Retinopathy of Prematurity and Bronchopulmonary Dysplasia are Independent Antecedents of Cortical Maturational Abnormalities in Very Preterm Infants

Very preterm (VPT) infants are at high-risk for neurodevelopmental impairments, however there are few validated biomarkers at term-equivalent age that accurately measure abnormal brain development and predict future impairments. Our objectives were to quantify and contrast cortical features between full-term and VPT infants at term and to associate two key antecedent risk factors, bronchopulmonary dysplasia (BPD) and retinopathy of prematurity (ROP), with cortical maturational changes in VPT infants. We prospectively enrolled a population-based cohort of 110 VPT infants (gestational age ≤31 weeks) and 51 healthy full-term infants (gestational age 38–42 weeks). Structural brain MRI was performed at term. 94 VPT infants and 46 full-term infants with high-quality T2-weighted MRI were analyzed. As compared to full-term infants, VPT infants exhibited significant global cortical maturational abnormalities, including reduced surface area (−5.9%) and gyrification (−6.7%) and increased curvature (5.9%). In multivariable regression controlled for important covariates, BPD was significantly negatively correlated with lobar and global cortical surface area and ROP was significantly negatively correlated with lobar and global sulcal depth in VPT infants. Our cohort of VPT infants exhibited widespread cortical maturation abnormalities by term-equivalent age that were in part anteceded by two of the most potent neonatal diseases, BPD and ROP.


Methods
Subjects. We prospectively enrolled 110 VPT infants from four level III NICUs from December 2014 to April 2016. These NICUs -Nationwide Children's Hospital, Ohio State University Medical Center, Riverside Hospital, and Mount Carmel St. Ann's Hospital -cared for approximately 80% of VPT births in the Columbus, Ohio region. We prospectively enrolled 51 healthy FT control infants from September 2014 to August 2015 from the well-baby nurseries of these hospitals. The inclusion criteria for full term infants included a GA of 38-42 weeks, appropriate weight for GA, and a normal well-baby nursery course. Full-term exclusions included congenital/ chromosomal anomalies of the brain, spine, heart, or lungs, significant maternal conditions (e.g. insulin-dependent diabetes or severe preeclampsia), intrauterine drug or alcohol exposure, and a history of perinatal distress or complications. Preterm inclusion criteria included a GA of ≤31 weeks, while exclusion criteria included congenital/chromosomal anomalies of the brain, spine, or heart (e.g. Dandy-Walker malformation, myelomeningocele, cyanotic heart disease). We excluded any infants hospitalized at 44 weeks postmenstrual age (unless they were at Nationwide Children's Hospital, the site of imaging). The Nationwide Children's Hospital Institutional Review Board approved the study (Protocol#: IRB13-00636). Reciprocity agreements with the Institutional Review Boards of Nationwide Children's Hospital, Ohio State University Medical Center, Riverside Hospital, and Mount Carmel St. Ann's Hospital allowed for approval at all sites. The study methods were carried out in accordance with the relevant guidelines and regulations. Written informed consent was obtained from a parent or guardian of all infants.
imaging methods. MRI data was acquired on a 3T Siemens Skyra scanner at Nationwide Children's Hospital at a mean (SD) postmenstrual age of 40.2 (0.6) weeks for the VPT cohort and 41.3 (0.8) weeks for the FT cohort. The majority of infants at Nationwide were imaged while inpatients, and all infants from the other NICUs were imaged after discharge. The scans were overseen by a skilled neonatal nurse and a neonatologist, who monitored heart rate and oxygen saturations throughout. We performed imaging without sedation using these procedures: the infant was fed 30 minutes prior to MRI, silicone earplugs were placed (Instaputty, E.A.R. Inc, Boulder, CO), and a blanket and vacuum immobilization device (MedVac, CFI Medical Solutions, Fenton, MI) helped promote natural sleep. We used these acquisition parameters: axial T2-weighted: echo time 147 ms, repetition time 9500 ms, flip angle 150°, resolution 0.93 × 0.93 × 1.0 mm 3 , scan time 4:09 min. Additional sequences such as T1and susceptibility-weighted images were acquired for qualitative assessment of abnormality.
MRi processing. All T2-weighted images were postprocessed using the developing Human Connectome Pipeline (dHCP) 32 developed for neonatal MRI, in order to automatically extract global and regional values for cortical features. The dHCP pipeline leverages FreeSurfer for tissue segmentation and volume estimation. It also calculates a value for the cortical features of interest, surface area, gyrification index, sulcal depth, and inner cortical curvature, for each of the 50 regions of the Gousias neonatal atlas 33 and for the whole brain. All tissue segmentations ( Fig. 1) and cortical surface results were visually inspected. Any scan with significant movement artifact, missing cortical boundaries, or moderate to severe ventriculomegaly was excluded at this point, as ventriculomegaly tended to interfere with proper tissue classification by the dHCP.
Statistical analysis. Our first major statistical analysis assessed significant differences in global and regional cortical features between the full-term and the very preterm groups. We began by transforming non-normal cortical feature values using either log transforms, power transforms, or the Yeo-Johnson transform 34 (for negative values), as appropriate. We then applied analysis of covariance (ANCOVA) to calculate significant group differences in the global and regional cortical features. We performed the ANCOVA analysis two ways: 1) with postmenstrual age (PMA) at MRI as the only covariate and 2) with PMA at MRI and total tissue volume (TTV) both as covariates. TTV was included to correct for the difference in brain size between the very preterm and the full term infants [35][36][37] . We did not use total cranial volume as a covariate in any analyses, because the very preterm group had significantly increased ventricular volume, despite the exclusion of subjects with ventriculomegaly. Following the ANCOVA, we used the 'adjust' command in Stata to center the covariates on their means and obtain the adjusted group means, which were then transformed back to their original units, if applicable. Because we analyzed so many cortical metrics (>200), we corrected for false discovery rate (FDR) using the Benjamini-Hochberg procedure with a chosen false discovery rate of 5%.
Our second major statistical analysis was a multivariable linear regression analysis examining the impact of BPD and ROP on all four global cortical metrics in the very preterm group, adjusted for covariates. We defined BPD (any severity) per Jobe and Bancalari's original definition 38 and ROP using the International Classification of ROP guidelines 39 . Because we hypothesized that BPD and ROP would predict cortical metric values independently of other known clinical or social demographic covariates, we took care to determine the best covariates for each regression. For each regression/global cortical metric, we started with 20 possible covariates (Table 1) from three distinct time periods -the antenatal, intrapartum, and postnatal periods -that had been associated with BPD, ROP, or NDI in previous research. We performed a sequential covariate selection method that has been published on several times previously 37,40,41 . Potential covariates were first tested with only other potential covariates from the same time period. For each individual period, covariates correlated with the cortical metric of interest at p < 0.25 were all entered into backward stepwise regression 42 , and were deemed significant for that period if p < 0.05 in multi-variable regression. Covariates deemed significant from all periods using the above criteria were included in the final model, even if they became non-significant with the addition of other covariates. Overall, this method 1) prevents later-occurring significant covariates from displacing earlier-occurring significant covariates and 2) helps avoid overfitting the model with too many covariates. Known highly-important covariates of brain development, PMA at MRI and injury score on structural MRI 43 , were forced into all final models. Next, we fit all resultant cortical surface metric models for each individual lobe, using multivariate regression in Stata, to determine whether the associations between BPD/ROP and cortical surface metrics were strongly regional.
Finally, we performed a third minor statistical analysis, using Pearson correlation analysis to assess the interrelation between all four global cortical metrics. All analyses were performed in STATA 15.1 (Stata Corp., College Station, TX), except the Benjamimi-Hochberg analysis, which was performed in python.

Figure 1.
An example T2-weighted MRI scan from one VPT subject is shown in sagittal, coronal, and horizontal views (top row). A segmentation from the Developing Human Connectome Pipeline (dHCP) is overlaid on the original image (bottom row). The dHCP performs cortical and sub-cortical volume segmentation, cortical surface extraction, and cortical surface inflation and was specifically designed for neonatal T2-weighted MRI brain scans.

Results
Of the 110 eligible VPT infants, one was missing a T2 image, 10 had significant ventriculomegaly that adversely affected brain segmentation. Four T2 images had significant artifact, and one had missing cortical boundaries. The remaining 94 infants were included in the analysis. Of these, only one infant had significant macrostructural injury on MRI. Of the 51 eligible FT subjects, four T2 images had missing boundaries, and one had poor alignment with the dHCP atlas. The remaining 46 were included in the analysis. Table 1 shows the baseline characteristics of our cohorts.
The results of our two ANCOVA group analyses were highly similar, therefore, we will present results controlled for both PMA at MRI and TTV. (Results controlled for PMA only are in the supplement.) Globally, compared to the full-term infants, the very preterm infants showed a 6 to 7% reduction in both surface area and gyrification index and a 6% increase in curvature. Whole-brain sulcal depth was not significantly different between the groups (Table 2).
Regionally, the very preterm cohort had reduced surface area bilaterally in the frontal, parietal, and temporal lobes and in the insula and the right occipital lobe compared to the full-term group. (Range: 3.1% to 9.0%). VPT surface area increased (13.5%) only in the left fusiform gyrus (Fig. 2A). The VPT infants had reduced gyrification in bilateral regions of all lobes (Range: 6.1% to 24.0%.) compared to the full-term group, with the largest decreases in the temporal lobes. (Fig. 2B). Very preterm curvature increased in bilateral regions of the frontal, parietal, and temporal lobes (7.1% to 85.7%.) compared to the full-term group, with extreme values in the bilateral fusiform gyrus (586.0% right; 238.6% left) (Fig. 2C), however given the small size of this structure, its magnitude may be spurious. Curvature decreased in only one region in the VPT group, the bilateral parahippocampal gyrus (16.5% right; 13.8% left). Although there was no significant difference between groups for global sulcal depth, our preterm group had many significant regional sulcal depth differences compared to the full-term group (Fig. 2D) in the frontal, parietal, and temporal lobes. (Range: −100.4% to 41.5%).
Our preterm only regression analysis uncovered correlations between BPD and surface area and ROP and sulcal depth, in bivariable analysis and in multivariable regression adjusted for covariates. For surface area, the significant clinical covariates from the sequential variable selection procedure 37,40,41 were injury score on MRI, PMA at MRI, sex, intrauterine growth restriction, head circumference at birth, and maternal breastmilk duration within the first 28 postnatal days. In a multivariable model with these covariates, BPD was negatively correlated with surface area (p = 0.009) ( Table 3). The significant clinical covariates chosen for sulcal depth were injury score on MRI, PMA at MRI, GA, and birth weight z-score. In a multivariable model with these covariates, ROP was www.nature.com/scientificreports www.nature.com/scientificreports/ negatively correlated with sulcal depth (p = 0.032) ( Table 3). ROP was not significantly correlated with surface area and BPD was not significantly correlated with sulcal depth in multi-variable regression with covariates. Likewise, BPD and ROP were both not significantly correlated with either gyrification index or curvature in multi-variable regression with covariates. Regionally, BPD was significantly negatively correlated with surface area in all lobes of the brain (Table 4). ROP was significantly negatively correlated with sulcal depth in the parietal lobe (Table 5).
Our four global cortical metrics were strongly and significantly interrelated. Sulcal depth had a significant negative relationship with curvature (r = −0.37) and a positive relationship with surface area (r = 0.34). Curvature was significantly negatively correlated with both gyrification index (r = −0.20) and surface area (r = −0.40). Gyrification index was positively correlated with surface area (r = 0.28). Sulcal depth and gyrification index were not significantly correlated.

Discussion
Our results extend published findings of significant cortical maturation abnormalities in very preterm infants at TEA. Overall, our very preterm group showed decreased global and regional gyrification, decreased global and regional surface area, and increased global and regional curvature compared to the full-term group. Our preterm group also showed regional alterations in sulcal depth compared to the full-term group, with decreases tending to dominate. Unlike prior studies, we also identified two important risk factors -BPD and ROP -for abnormal  Table 2. Whole-brain mean values and percent differences for four global cortical metrics in very preterm and full-term infants (corrected for PMA at MRI scan and TTV).

Figure 2. Mean Percent Difference in Cortical
Metrics (Full-term to Very Preterm). Regional percent differences in adjusted group means for surface area (panel A), gyrification index (panel B), curvature of the white matter surface (panel C), and sulcal depth (panel D). For regions with significant differences between groups (after false discovery rate correction), percent difference values (VPT − FT)/FT * 100 are projected onto a representative subject brain from the very preterm group. These values have been corrected for postmenstrual age at MRI scan and total brain tissue volume. For each panel,  www.nature.com/scientificreports www.nature.com/scientificreports/ changes in surface area and sulcal depth, respectively. In prior neuroimaging studies, preterm infants with BPD and ROP have not exhibited consistent structural abnormalities to suggest that ROP and BPD-associated NDI are mediated through qualitative abnormalities visible on structural MRI. Although we cannot determine a definitive causal link, our results suggest that aberrant cortical maturation could be in the causal pathway between severe neonatal illnesses like ROP and BPD and future NDI.
Reduced surface area in preterm subjects is well-documented at TEA 10,13 and in childhood 11,12 , and reduced gyrification index in preterm subjects is also well-documented at TEA 13,14,44 . Engelhardt et al. reported approximately 13% reduced gyrification and 19% reduced surface area in preterm infants compared to full-term infants at TEA, which is larger than our finding (7% reduced gyrification and 6% reduced surface area, globally). However, we excluded subjects with ventriculomegaly or injuries such as intraventricular hemorrhage (27% of the preterm subjects in Engelhardt's study had intraventricular hemorrhage), and we studied twice as many infants. Zhang reported 4% reduced gyrification index and 10% reduced cortical surface area in preterm children at seven years of age, which is closer to our result, albeit in a different age group. Zhang normalized their gyrification and surface area measures by convex hull area and total cortical surface area, respectively, while Engelhardt reported unnormalized group differences. Both normalization methods are fundamentally different from our analysis, in which we corrected for PMA at MRI and TTV. The differences in magnitude between our study and these previous studies can likely be explained by the injury status of the preterm cohort and the type of normalization performed.
Other works examining sulcal depth in very preterm subjects have reported just decreases 14 or both increases and decreases at TEA 13 and in childhood 15 , with decreases tending to dominate. Although there was no global difference between our two groups, there were more significant regional sulcal depth decreases (12) than increases (five) in the very preterm group, which is consistent with previous work.
Curvature was increased at TEA in our VPT group, which corroborates Makropoulos et al. 's findings 16 . Many different curvature formulations have been utilized, making this a challenging metric to compare to previous work. Curvature has been reported on both the outer cortical surface and on the inner cortical/outer white matter surface, as in the current work. It has been formulated as the mean of the principal curvatures as in the current work, Gaussian curvature, etc. (See Pienaar 45 for an in-depth review). Given the larger number of regions with significant sulcal depth decrease in the very preterm group, immature sulci may still be driving increased curvature. Shallower sulci are missing sulcal wall area (low curvature), but have the same amount of sulcal trough area (high curvature), which results in higher average curvature in that local region 14 . The negative relationship between 1) sulcal depth and curvature and 2) gyrification index and curvature reinforces the idea that decreased folding may increase curvature in the very preterm group.
Our preterm only regression analyses confirmed BPD and ROP as clinical risk factors of abnormal cortical maturation in very preterm infants. A significant relationship persisted after controlling the final multivariable  www.nature.com/scientificreports www.nature.com/scientificreports/ models for carefully selected clinical covariates. The negative relationship between BPD and surface area was widespread, involving all lobes of the brain (Table 4). Our lobar analysis also revealed that parietal lobe sulcal depth was the main driver of the negative relationship between this cortical metric and ROP (Table 5). Although BPD is a well-established, independent risk factor for NDI, no structural abnormalities except ventriculomegaly consistently explain the poor outcomes associated with BPD 27,46 . BPD and ROP both predicted NDI independently of structural injury in a large cohort of extremely preterm infants 21 , but our study is the first to link both to cortical maturational abnormalities. A few studies linked these antecedents to decreased white matter maturation, but none using the objective MRI measures of our study. Both BPD and ROP were associated in univariable analyses with lower overall maturation score (evaluated visually/qualitatively) on TEA MRI in preterm infants 47 , but only BPD remained significant in multivariable analyses. Both BPD and ROP share systemic inflammation and hyperoxia as common mechanisms that could disrupt brain development. It is challenging to determine whether the link between BPD/ROP and cortical dysmaturation is causal or a result of the most developmentally delayed infants requiring more intervention and support. However, in a rat pup model, hyperoxia caused BPD-like lung changes, ROP-like hypervascularization, and reduced brain surface area 48 , suggesting that BPD and ROP are interrelated and cause brain changes beyond those induced by prematurity. More work is needed to confirm similar causal effects in preterm infants.
Our study's strengths include a population-based cohort, assessment of quantitative cortical maturation measurements, and many covariates collected and temporally examined. However, several limitations should be acknowledged. Despite the large number of possible covariates collected, residual confounding could have introduced bias. Nevertheless, the association we identified with BPD and cortical surface area was especially strong. Furthermore, because we assessed numerous metrics, some false discovery is inevitable, even after applying FDR correction. Given that we excluded severely brain injured subjects, including those with ventriculomegaly, the brain changes reported here may be less severe than in other brain injured preterm cohorts. Finally, any automated segmentation is inherently imperfect. The dHCP pipeline represents a substantial improvement over labor-intensive manual segmentation, however its creators acknowledge a 2% rate of significant error 32 .
We identified widespread cortical maturation abnormalities at TEA in very preterm infants compared to full-term infants. These cortical abnormalities, specifically surface area and sulcal depth, appear to be exacerbated by clinical risk factors, BPD and ROP, respectively. Larger studies are needed to validate the relationship of these antecedents with cortical dysmaturation. Further work should examine the predictive value of cortical dysmaturation for NDI in VPT infants. This will allow clinicians to assess risk of impairment early and provide targeted early interventions.

Data availability
The datasets generated during the current study are available from the corresponding author on reasonable request.