Theory and implementation of the gauge-including atomic orbital (CIAO) ansatz for the gauge-invariant calculation of nuclear magnetic resonance chemical shifts are described for the coupled-cluster singles and doubles (CCSD) approach. Results for the shielding constants of the hydrides HF, H2O, NH3, and CH4 as well as for a few multiply bonded systems such as CO, N-2, and HCN demonstrate the importance of higher-order correlation corrections, as good agreement with experiment is only obtained at the CCSD level and to some extent at partial fourth-order many-body perturbation theory [SDQ-MBPT(4)] with the latter slightly overestimating correlation effects due to single and double excitations. For relative chemical shifts, GIAO-CCSD calculations provide in difficult cases (e.g., CO and CF4) more accurate results than previous GIAO-MBPT(2) calculations. But, it seems that it is often more important to include rovibrational effects (as well as possible molecule-solvent interactions) than higher-order correlation corrections. Despite that, GIAO-CCSD proves to be a powerful tool for the accurate calculation of NMR chemical shifts. Its capabilities as well as its limitations are demonstrated in shielding calculations for formaldehyde, diazomethane, and ozone. At least for the latter, the description provided by the CCSD ansatz is not sufficient and even higher excitations need to be considered.