Participant

Our participant (E.R.C.) was a healthy 38-year-old primiparous woman who underwent in-vitro fertilization (IVF) to achieve pregnancy. Previous studies reported no observable differences in neural changes from prepregnancy to postpregnancy between women who conceived naturally versus women who conceived via IVF13, and doing so provides a controlled way of monitoring pregnancy status. The participant experienced no pregnancy complications (for example, gestational diabetes and hypertension), delivered at full term via vaginal birth, nursed through 16 months postpartum, and had no history of neuropsychiatric diagnosis, endocrine disorders, prior head trauma or history of smoking. The participant gave written informed consent and the study was approved by the University of California, Irvine Human Subjects Committee.

Study design

The participant underwent 26 MRI scanning sessions from 3 weeks before conception through 2 years postpartum (162 weeks), during which high-resolution anatomical and diffusion spectrum imaging scans of the brain were acquired. Scans were distributed throughout this period, including prepregnancy (four scans), first trimester (four scans), second trimester (six scans), third trimester (five scans) and postpartum (seven scans; Fig. 1c). The first 6 sessions took place at the UCSB Brain Imaging Center (BIC), the final 20 sessions took place at the UCI Facility for Imaging and Brain Research (FIBRE). The majority of scans took place between 9 AM and 2 PM, limiting significant AM–PM fluctuations49. The MRI protocol, scanner (Siemens 3T Prisma) and software (version MR E11) were identical across sites. Each scanner was checked weekly for the duration of the study and passed all QC reports indicating no significant alterations in the geometry. To ensure the robustness of the findings, after the final study session, the participant completed back-to-back validation scans at UCI and UCSB within a 12-h window to assess reliability between scanners. Intraclass correlation coefficients (two-way, random effects, absolute agreement, single rater) reveal ‘excellent’ test–retest reliability between scanners, including ROI-level GMV (ICC = 0.97, 95% CI: 0.80–0.99), ROI-level CT (ICC = 0.96, 95% CI: 0.90–0.98), MTL subfield volume (ICC = 0.99, 95% CI: 0.97–0.99) and ROI-level QA (ICC = 0.94, 95% CI: 0.91–0.97). Furthermore, when examining the relationship between gestation week and GMV among UCI-only gestational sessions, findings were consistent (Supplementary Fig. 12), indicating that site differences are highly unlikely to have contributed meaningfully to the observed effects. Although not applicable here, we note that having a control participant scanned over a similar duration within the same scanner is critical for estimating how much variation in the brain can be attributed to within-scanner variability.

To monitor state-dependent mood and lifestyle measures, the following scales were administered on each experiment day: Perceived Stress Scale50, Pittsburgh Sleep Quality Index51, State-Trait Anxiety Inventory for Adults52 and Profile of Mood States53. Correlation analyses between state-dependent measures, summary brain metrics and gestation week revealed little to no relationships. The only exception to this was a moderate negative association between global QA and state anxiety (Spearman’s correlation (ρ) = −0.65, q = 0.04; baseline—36 weeks, n = 16). By making this data openly accessible, we encourage a more nuanced approach toward exploring mood and lifestyle measures in relation to brain changes over pregnancy.

Endocrine procedures

The participant underwent a blood draw (n = 19; Fig. 1c) before MRI scanning. Sex steroid concentrations were determined via ultra-sensitive liquid chromatography–mass spectrometry at the Brigham and Women’s Hospital Research Assay Core (BRAC). Assay sensitivities, dynamic range and intra-assay coefficients of variation were as follows: estradiol—1.0 pg ml−1, 1–500 pg ml−1, <5% relative s.d. (RSD); progesterone—0.05 ng ml−1, 0.05–10 ng ml−1, 9.33% RSD. Serological samples were not acquired in five sessions due to scheduling conflicts with UC Irvine’s Center for Clinical Research.

MRI acquisition

MRI scanning sessions at the University of California, Santa Barbara and Irvine were conducted on 3T Prisma scanners equipped with 64-channel phased-array head/neck coil (of which 50 coils are used for axial brain imaging). High-resolution anatomical scans were acquired using a T1-weighted (T1w) magnetization prepared rapid gradient echo (MPRAGE) sequence (repetition time (TR) = 2,500 ms, time to echo (TE) = 2.31 ms, inversion time (TI) = 934 ms, flip angle = 7°, 0.8 mm thickness) followed by a gradient echo field map (TR = 758 ms, TE1 = 4.92 ms, TE2 = 7.38 ms, flip angle = 60°). A T2-weighted (T2w) turbo spin echo scan was also acquired with an oblique coronal orientation positioned orthogonally to the main axis of the hippocampus (TR/TE = 9,860/50 ms, flip angle = 122°, 0.4 × 0.4 mm2 in-plane resolution, 2-mm slice thickness, 38 interleaved slices with no gap, total acquisition time = 5 min and 42 sec). The Diffusion Spectrum Imaging (DSI) protocol sampled the entire brain with the following parameters: single phase, TR = 4,300 ms, echo time = 100.2 ms, 139 directions, b-max = 4,990, FoV = 259 × 259 mm, 78 slices, 1.7986 × 1.7986 × 1.8 mm voxel resolution. These images were linearly registered to the whole-brain T1w MPRAGE image. A custom foam headcase was used to provide extra padding around the head and neck, as well as to minimize head motion. Additionally, a custom-built sound-absorbing foam girdle was placed around the participant’s waist to attenuate sound near the fetus during second-trimester and third-trimester scanning.

Image processing Cortical volume and thickness

CT and GMV were measured with Advanced Normalization Tools54 version 2.1.0 (ANTs). We first built a subject-specific template (SST) (antsMultivariateTemplateConstruction2) and tissue priors (antsCookTemplatePriors) based on our participant’s two preconception whole-brain T1-weighted scans to examine neuroanatomical changes relative to the participant’s prepregnancy baseline. We used labels from the OASIS population template, provided by ANTs, as priors for this step. For each session, the structural image was processed and registered to the SST using the ANTs CT pipeline (antsCorticalThickness). This begins with an N4 bias field correction for field inhomogeneity, then brain extraction using a hybrid registration/segmentation method55. Tissue segmentation was performed using Atropos54 to create tissue masks of CSF, gray matter, white matter and deep gray matter. Atropos allows prior knowledge to guide the segmentation algorithm, and we used labels from our SST as priors to minimize warping and remain in native participant space. CT measurements were then estimated using the DiReCT algorithm56, which estimates the gray–white matter interface and the gray matter–CSF interface and computes a diffeomorphic mapping between the two interactions, from which thickness is derived. Each gray matter tissue mask was normalized to the template and multiplied to a Jacobian image that was computed via affine and nonlinear transforms. Using MATLAB (version 2022a), summary, regional-level estimates of CT, GMV and CSF for each scan were obtained by taking the first eigenvariate (akin to a ‘weighted mean’57) across all voxels within each parcel of the Schaefer 400-region atlas58. We then averaged ROIs across networks, which were defined by the 17-network Schaefer scheme58,59. Global measures of CT, GMV and CSF were computed for each session by summing across all voxels within the respective output image; total brain volume was computed by summing across all voxels within each session’s brain extraction mask. Our findings held when using an SST derived from all 26 MRIs (prepregnancy through postpartum), as well as when estimating the mean (versus weighted mean) of all voxels within each parcel. The ANTs CT pipeline is highly validated with good test–retest reproducibility and improved ability to predict variables such as age and gender from region-wise CT measurements compared to surface-based FreeSurfer55. However, to reproduce our findings across software packages, we also ran the T1w data through the longitudinal FreeSurfer (v.7) CT pipeline60,61, which corroborated our findings using both the Schaefer-400 (Supplementary Fig. 2 and Supplementary Tables 1 and 4) and popular Desikan–Killiany62 (Supplementary Table 3) cortical parcellations. Whole-brain T1w-based subcortical volume estimates (including cerebellum and lateral ventricles) were also derived using this FreeSurfer pipeline, wherein we derived 28 region-of-interest estimates via the commonly used ‘aseg’ parcellation scheme63 (Supplementary Fig. 6a). A complete reporting of findings can be found in Supplementary Data 1.

Mean framewise displacement (FWD) estimates from gestation sessions with a 10-min resting-state scan (n = 18) were used to indirectly assess whether motion increased throughout pregnancy. Average FWD (mm) was extremely minimal across the entire experiment (M = 0.13, s.d. = 0.02, range = 0.09–0.17) and varied only slightly by pregnancy stage (pre, M = 0.11 and s.d. = 0.004; first, M = 0.11 and s.d. = 0.01; second, M = 0.13 and s.d. = 0.02; third, M = 0.16 and s.d. = 0.007; post, M = 0.13 and s.d. = 0.01). While mean FWD did correspond with gestation week (r = 0.88, P < 0.001), controlling for this did not alter our main findings (for example, total GMV remained negatively associated with gestation after partial correlation with FWD (r = −0.87 and P < 0.001) because motion differences between stages were minuscule (Supplementary Fig. 4a).

As a further test of the robustness of the dataset, we ran QC assessments on all T1w images using the IQMs pipeline64 from MRIQC (version 23.1). Assessments of interest included (1) coefficient of joint variation (CJV), (2) signal-to-noise ratio for gray matter (SNR) and (3) contrast-to-noise ratios (CNR). All QC metrics fell within expected standard ranges65 (Supplementary Fig. 4b–d). Although relationships existed between gestation week and QC measures (CJV, r = 0.70 and P < 0.001; SNR and CNR, r = −0.83 and P < 0.001), including these variables in the regression models did not detract from our finding suggesting cortical GMV reductions occur over gestation, especially within regions belonging to attention and somatosensory networks (Supplementary Fig. 5). When looking across all MRIQC outputs, discrepancies were noted in session seven (gestation week nine, first trimester). Removing this day from the analyses only strengthened observed relationships between cortical volume and gestation; however, for completeness, data from this day is included in the main findings. These QC outputs for each session of the experiment can be found in Supplementary Data 1. Finally, we used FreeSurfer’s Eueler number to evaluate a field-standard quantitative assessment of each T1w structural image66. We observed no significant relationships between the Euler number and gestation week or summary brain metrics. A discrepancy (for example, two s.d. below average) was noted in session eight; however, again, removing this session did not detract from our main findings showing reductions in GMV over gestation.

Hippocampal segmentation

T1- and T2-weighted images (n = 25) were submitted to the automatic segmentation of hippocampal subfields package (ASHS67, version July 2018) for parcellation of seven MTL subregions: CA1, CA2/CA3, dentate gyrus, subiculum, perirhinal cortex, entorhinal cortex and PHC (Supplementary Fig. 6b). The ASHS segmentation pipeline automatically segmented the hippocampus in the T2w MRI scans using a segmented population atlas, the Princeton Young Adult 3T ASHS Atlas template68 (n = 24, mean age = 22.5 years). A rigid-body transformation aligned each T2w image to the respective T1w scan for each day. Using ANTs deformable registration, the T1w was registered to the population atlas. The resulting deformation fields were used to resample the data into the space of the left and right template MTL ROI. Within each template ROI, each of the T2w scans of the atlas package was registered to that day’s T2w scan. The manual atlas segmentations were then mapped into the space of the T2w scan, with segmentation of the T2w scan computed using joint label fusion69. Finally, the corrective learning classifiers contained in ASHS were applied to the consensus segmentation produced by joint label fusion. The output of this step is a corrected segmentation of the T2w scan. Further description of the ASHS protocol can be found here67. T2w scans and segmentations were first visually examined using ITK-SNAP70 for quality assurance and then subjected to manual editing in native space using ITK-SNAP (v.3.8.0-b; C.M.T.). One session (scan 15, third trimester) was discarded due to erroneous scan orientation. The anterior extent of the segmented labels was anchored 4 mm (two slices) anterior to the appearance of the limen insulae, and the posterior extent was anchored to the disappearance of hippocampal gray matter from the trigone of the lateral ventricle. Boundaries between perirhinal, entorhinal and parahippocampal cortices were established in keeping with the Olsen–Amaral–Palombo (OAP) segmentation protocol71. In instances where automatic segmentation did not clearly correspond to the underlying neuroanatomy, such as when a certain label was missing several gray matter voxels, manual retouching allowed for individual voxels to be added or removed. All results are reported using the manually retouched subregion volumes to ensure the most faithful representation of the underlying neuroanatomy. Scans were randomized and segmentation was performed in a random order, blind to pregnancy stage. To assess intrarater reliability for the present analyses, two days underwent manual editing a second time. The generalized Dice similarity coefficient72 across subregions was 0.87 and the intraclass correlation coefficient was 0.97, suggesting robust reliability in segmentation.

White matter microstructure

Diffusion scans were preprocessed using the automation software QSIprep (version 0.16.1) compiled using a singularity container73 and run primarily with the default parameters, with the exceptions ‘–output resolution 1.8’, ‘–dwi denoise window 5′,–force-spatial-normalization’, ‘–hmc model 3dSHORE’, ‘–hmc-transform Rigid’ and ‘–shoreline iters 2’. Twenty-three sessions were preprocessed and analyzed, with the remaining three scans excluded due to missing DSI scans (sessions 9 and 15) or corresponding field map for distortion correction (session 7). Despite passing QC assessments during preprocessing, visual inspection of the field maps in session 10 revealed a slight artifact. However, removal of this session had minimal impact on the overall results and remained in the final analyses. T1w images were corrected for intensity nonuniformity (N4BiasFieldCorrection) and skull-stripped (antsBrainExtraction). The images underwent spatial normalization and registration to the ICBM 152 Nonlinear Asymmetric template. Finally, brain tissue segmentation of CSF, GM and WM was performed on each brain-extracted T1w using FMRIB’s Automated Segmentation Tool (FAST). Preprocessing of diffusion images began by implementing MP-PCA denoising with a 5-voxel window using MRtrix3’s dwidenoise function. B1 field inhomogeneity was corrected using dwibiascorrect from MRtrix3 with the N4 algorithm. Motion was corrected using the SHORELine method. Susceptibility distortion correction was based on GRE field maps. Preprocessed Nifti scans were prepared for tractography using DSI Studio via singularity container version Chen-2022-07-31 (ref. 74). Diffusion images were converted to source code files using the DSI Studio command line ‘--action=src’ and a custom script to convert all images. The diffusion data were reconstructed in MNI space using q-space diffeomorphic reconstruction75 with a diffusion sampling of 1.25 and output resolution of 1.8 mm isotropic. The following output metrics were specified to be included in the output FIB file: QA and mean diffusivity (MD). The quality and integrity of reconstructed images were assessed using ‘QC1: SRC Files Quality Control’. First, the consistency of image dimension, resolution, DWI count and shell count was checked for each image. Second, each image was assessed for the ‘neighboring DWI correlation’ which calculates the correlation coefficient of low b DWI volumes that have similar gradient direction. Lower correlation values may indicate issues with the diffusion signal due to artifacts or head motion. Finally, DSI Studio performed an outlier check, labeling images as a ‘low-quality outlier’ if the correlation coefficient was >3 s.d. from the absolute mean. None of our scans were flagged as outliers. The reconstructed participant files were aggregated into one connectometry database per metric.

Day2Day control dataset

To compare our findings against a control group of nonpregnant densely-sampled individuals, we used the Day2Day dataset23 which offered comparable whole-brain T1 and T2 MTL scans for eight participants (two male) scanned 12–50 times over 2–7 months. Each participant was run through the ANTs CT and ASHS processing pipelines as outlined above (‘Cortical volume and thickness’ and ‘Hippocampal segmentation’). To note, for each participant, we created an SST based on their first two sessions for consistency with the primary dataset; subfield volumes for the T2 MTL scans did not undergo manual retouching. Due to missing header information on the publicly available diffusion scans, we were unable to benchmark our white matter changes with the Day2Day dataset.

Statistical analysis

Statistical analyses were conducted using R (sMRI; version 3.4.4) and DSI Studio (dMRI; Chen-2022-07-31).

Summary brain metrics

To reflect the existing literature, we first explored brain metrics across the entire study duration (prepregnancy through postpartum, n = 26 scans). When including all sessions, total brain volume, GMV, CT, global QA, ventricle volume and CSF displayed nonlinear trends over time; therefore, we used generalized additive models (GAM; cubic spline basis, k = 10, smoothing = GCV), a method of nonparametric regression analysis (R package, mgcv76), to explore the relationship between summary brain metrics (outcome variables) and gestation week (smooth term). Each model underwent examination (gam.check function) to ensure it was correctly specified with regards to (1) the choice of basis dimension (k) and (2) the distribution of model residuals (see mgcv documentation in ref. 76). The general pattern of results held after toggling model parameters; however, we note the risk of overinterpreting complex models with small sample sizes77. To address overfitting and cross-validate our basis type selection, we also fit the data using nonpenalized general linear models (GLM) with both linear and polynomial terms for gestation week. We compared the performance of each GLM (that is, models using only a linear term versus models with polynomial terms) via the Akaike information criterion (AIC), which revealed that cubic models consistently outperformed both linear and quadratic models (AICdiff > 3), providing additional evidence for nonlinear changes in structural brain variables over time. Determining whether these patterns replicate in larger cohorts and whether complex models are better suited to capture data patterns across individuals will be a necessary next step.

Cortical GMV and CT

We then narrowed our analyses to the first 19 sessions (baseline—36 weeks gestation) to assess novel brain changes occurring over the gestational window. We first computed Pearson’s product-moment correlation matrices between the following variables: gestation week, estradiol, progesterone and the 17 network-level average GMV values. We then ran a multivariate regression analysis predicting ROI-level GMV changes by gestation week. To identify which regions were changing at a rate different from the global decrease, we then ran the analyses again to include total GMV in the regression model (Supplementary Table 2). This was extended to the network level, where we ran partial correlations accounting for total GMV. These same analyses were then run with CT measures. Globally-corrected results provided in Supplementary Tables 1–5. Percent change at the network level was computed by subtracting the final pregnancy value (36 weeks pregnant) from the first prepregnancy baseline value, then dividing that difference by said first prepregnancy baseline value. All analyses underwent multiple comparisons testing (false discovery rate (FDR)-corrected at q < 0.05).

Subcortical GMV

A similar statistical approach was taken for subcortical volume estimates. We ran a multivariate regression analysis predicting GMV changes over gestation in 28 ROIs (Supplementary Fig. 6a) by gestation week (FDR-corrected at q < 0.05).

To evaluate the relationship between gestation week and MTL subregion volume over pregnancy (n = 7 bilateral subregions and n = 18 MTL scans), we used a combination of linear and nonlinear models based on individual subregion data patterns. Models were compared for best fit with each subregion via AIC from the GLM output (as described in ‘Summary brain metrics’). A linear regression model was most appropriate for PHC (AICdiff < 3), whereas a quadratic model performed best for CA1 and CA2/CA3. As a control, we repeated the analyses with MTL subregion volumes after proportional volume correction of total GMV calculated by ASHS. Finally, we evaluated the relationship between endogenous sex hormones (estrogen and progesterone) and subregion volumes using linear regression. Relationships were considered significant only if they met FDR correction at q < 0.05.

White matter microstructure

DSI Studio’s correlational tractography74 was used to analyze the relationship between white matter structure and gestational week (n = 16). A truncated model was run to examine the relationship between white matter and sex steroid hormones (n = 14) for the subset of diffusion scans with paired endocrine data during gestation. A nonparametric Spearman’s correlation was used to derive the correlation between gestational week and endocrine factors and our metrics of interest (QA and MD; see Supplementary Table 9 and Supplementary Fig. 10 for MD results) because the data were not normally distributed. Statistical inference was reached using connectometry, a permutation-based approach that tests the strength of coherent associations found between the local connectome and our variables of interest. It provides higher reliability and replicability by correcting for multiple comparisons. This technique provides a high-resolution characterization of local axonal orientation. The correlational tractography was run with the following parameters: t score threshold of 2.5, four pruning iterations and a length threshold of 25 voxel distance. To estimate the FDR, a total of 4,000 randomized permutations were applied to obtain the null distribution of the track length. Reported regions were selected based on FDR cutoff (FDR < 0.2, suggested by DSI Studio), and contained at least ten tracts. For visualization of global and tract QA at each gestational stage, mean QA values were extracted using DSI Studio’s whole-brain fiber tracking algorithm and ROI-based tracking using the default HCP842 atlas78.

Day2Day dataset: measurement variability

To establish a marker of normative variability over half a year, we computed metrics of measurement variability using the Day2Day dataset23, which provided both whole-brain T1 and high-resolution T2 MTL scans. For each region, j, of the Schaefer parcellation, we assessed across-session variability, ε, as

$$_=100\times }\left(\frac_-\hat\right|}}\right)$$

Where ts is the morphometric measurement of a parcel for session s and \(\hat\) is the mean of t across sessions55,79. Thus, we defined variability as the mean absolute percent difference between each individual and the mean across sessions. Across-session variability estimates for all 400 regions were then averaged across eight participants, and a global measure of cortical GMV variability was computed by averaging across the 400 regions. This approach was repeated independently for the T2 hippocampal scans, wherein we computed across-session variability for each parcel of the ASHS parcellation scheme (n = 7 bilateral subfields). However, it is important to note that raw subfield values (that is, no manual retouching) were used for Day2Day variability assessments and should be interpreted with caution. Finally, to better compare against our own data, we repeated this approach using our participant’s first two baseline scans (that is, preconception) to derive within-participant variability estimates.

Benchmarking our data in this way allows us to capture the degree of change expected due to factors such as image processing and instrumentation variability or other day-to-day changes that could potentially modulate brain size and shape (see ref. 80 for review). The percent change observed over pregnancy (baseline versus 36 weeks gestation) far exceeds the expected variability estimated using both the Day2Day dataset (Supplementary Fig. 11) and our within-participant control data. This was quantified by dividing the observed percent change in GMV metrics (baseline versus 36 weeks) by the global measure of GMV percent variability of each control group (that is, Day2Day, within-participant control), independently for cortex and subcortex.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.