Prenatal Exposure to Metabolism-Disrupting Chemicals, Cord Blood Transcriptome Perturbations, and Birth Weight in a Belgian Birth Cohort

Prenatal exposure to metabolism-disrupting chemicals (MDCs) has been linked to birth weight, but the molecular mechanisms remain largely unknown. In this study, we investigated gene expressions and biological pathways underlying the associations between MDCs and birth weight, using microarray transcriptomics, in a Belgian birth cohort. Whole cord blood measurements of dichlorodiphenyldichloroethylene (p,p’-DDE), polychlorinated biphenyls 153 (PCB-153), perfluorooctanoic acid (PFOA), perfluorooctane sulfonic acid (PFOS), and transcriptome profiling were conducted in 192 mother–child pairs. A workflow including a transcriptome-wide association study, pathway enrichment analysis with a meet-in-the-middle approach, and mediation analysis was performed to characterize the biological pathways and intermediate gene expressions of the MDC–birth weight relationship. Among 26,170 transcriptomic features, we successfully annotated five overlapping metabolism-related gene expressions associated with both an MDC and birth weight, comprising BCAT2, IVD, SLC25a16, HAS3, and MBOAT2. We found 11 overlapping pathways, and they are mostly related to genetic information processing. We found no evidence of any significant mediating effect. In conclusion, this exploratory study provides insights into transcriptome perturbations that may be involved in MDC-induced altered birth weight.


Introduction
Metabolism-disrupting chemicals (MDCs) have been defined as natural or anthropogenic endocrine-disrupting chemicals (EDCs) that can promote metabolic changes and ultimately lead to obesity, type 2 diabetes and/or non-alcoholic fatty liver disease (NAFLD) [1]. In line with the Developmental Origins of Health and Disease (DOHaD) hypothesis [1], the prenatal period is a highly sensitive and vulnerable phase during which stressors, such as MDCs, can alter cell numbers and fate, gene expression, and protein levels that may lead to changes in tissue and organ function and contribute to increased susceptibility to a variety of non-communicable diseases later in life [2]. This may be the result of differences in toxicokinetics between children and adults and from time-dependent programming during early development [3].
Both high and low birth weight (HBW and LBW) are considered important predictors of later perturbed metabolic outcomes in children and adults [4][5][6]. While some observational studies have demonstrated associations between exposure to MDCs [including dichlorodiphenyldichloroethylene (p,p'-DDE), polychlorinated biphenyl-153 (PCB-153), perfluorooctanoic acid (PFOA), and perfluorooctane sulfonic acid (PFOS)] and birth weight [7][8][9][10], the molecular mechanisms of action remain poorly understood. The field of omics, based on high-throughput biochemical data, provides promising opportunities to advance and enhance our understanding of the impact of MDCs on child health, including by revealing changes in the gene expression using transcriptome profiling [11,12].
Assessing the effects of various chemical exposures on gene expression may help to uncover cellular mechanisms through which exposures influence the development of metabolic disorders in human populations. Several recent epidemiological studies using transcriptomics data have increased our understanding of how exposure to MDCs may perturb gene expression, and have identified regulatory pathways that may be affected by these exposures [13][14][15], as well as links between gene expression and birth weight [16][17][18][19]. However, to our knowledge, a study assessing the transcriptome in relation to both MDCs and birth weight in the same study population has not been performed.
Based on results from our previous birth cohort study [15], several MDCs (p,p'-DDE, PCB-153, PFOA, and PFOS) were suggested to play a role in transcriptional changes which are related to metabolic health outcomes. This led us to hypothesize that prenatal exposure to MDCs induces transcriptional modifications that, in turn, affect birth weight and have adverse effects on human health. Here, we aim to identify transcriptomic alterations in the cord blood of Belgian mother-child pairs that are associated with both prenatal MDC levels and birth weight in order to better understand the molecular effects and the underlying mechanisms.

Population Characteristics
Demographic and exposure information for participants are shown in Tables 1 and S1. The median gestational age was 40 weeks. Most children (98%) had a birth weight at or more than 2500 g, with a median of 3540 g. The median concentrations were 75.9 ng/g lipid, 28.7 ng/g lipid, 1600 ng/L, and 2700 ng/L for p,p'-DDE, PCB-153, PFOA, and PFOS, respectively (Table S1). The majority of the mothers had completed a high level of education (59%), had a normal pre-pregnancy body mass index (BMI) between 18.5 and 25 kg/m 2 (71%), and did not smoke during pregnancy (85%). In addition, 38% of mothers were nulliparous and 57% were above 30 years of age at delivery.

Gene Expression Associated with MDCs and Birth Weight
Using a transcriptome-wide association study (TWAS) approach, we failed to select any features from models (1) or (2) with significance levels of false discovery rate (FDR) <0.05 or 0.20, and selected only a few features with a stringent p-value < 0.01 (Table 2). In order to avoid excluding weak but possibly relevant features, we used a relatively lenient p-value < 0.05 to select features for further analyses as an exploratory study. With p-value < 0.05, we found that 2110 out of 26,170 features were associated with one or more MDCs (777, 623, 333, and 624 for p,p'-DDE, PCB-153, PFOA, and PFOS, respectively; Table 2), and 775 features were associated with birth weight. A similar number of associated features were found in the sensitivity analyses of gestational age-unadjusted MDC-transcriptome associations (Table S2). In addition, as shown in the volcano plots (Figure S1a-e), the significance and directionality of gene expression obtained with and without adjustment for gestational age were consistent in the TWAS models for MDCs and features.

Pathways Associated with MDCs and Birth Weight
The MDC-or birth weight-associated pathways at FDR < 0.05 with at least five genes involved are represented in Table S4. There were 17, 3, 27, 20, and 33 pathways associated with p,p'-DDE, PCB-153, PFOA, PFOS, and birth weight, respectively; most of them were related to genetic information processing and organismal systems. Notably, one metabolic pathway [glycosaminoglycan biosynthesis] was linked to p,p'-DDE; three [metabolism of xenobiotics by cytochrome P450, drug metabolism, and type 1 diabetes mellitus (T1D)] were linked to PFOA, two [amino sugar and nucleotide sugar metabolism, type I diabetes mellitus] were linked to PFOS, and five [oxidative phosphorylation (OXPHOS), non-alcoholic fatty liver disease (NAFLD), cysteine and methionine metabolism, sulfur metabolism, and valine, leucine, and isoleucine degradation] were linked to birth weight.
At FDR < 0.05, we found that four, three, six, and three pathways associated with birth weight ovrlapped with p,p'-DDE, PCB-153, PFOA, and PFOS, respectively (Figure 1). They mostly belong to the "genetic information processing" category in the Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway Database [22], and none of them were metabolismrelated pathways ( Table 4). The PC1 scores used to represent pathways in the mediation analysis explained 37-57% of the variance in the involved genes, and given the insignificant average causal mediation effects (ACMEs) with large FDR values, we did not observe any pathway that mediated both MDC and birth weight (Table 4).

Discussion
Transcriptome changes in early life may act in response to environmental exposures and subsequently lead to adverse health outcomes later in life; however, epidemiological studies are scarce. This is the first paper that evaluated the cord blood transcriptome with MDC exposures and birth weight. We examined differences in transcriptomics at the gene and pathway levels.
The five gene expressions that are metabolism-related and were found to be associated with both an MDC and birth weight are BCAT2, IVD, SLC25a16, HAS3, and MBOAT2. Birth weight may be altered by an MDC through one of these gene expressions, although we did not find a mediating effect to be significant. Branched-chain amino acids (BCAAs) are associated with the progression of obesity-related metabolic disorders [23]; additionally, BCAA catabolism is suggested to play a role in the pathogenesis of metabolic disturbances, and BCAT2 is an important enzyme that catalyzes the initial step of BCAA catabolism [24,25]. In a recent human study [26], BCAT2 variants were detected in Spanish infants suspected of having maple syrup urine disease-a rare metabolic disorder that some babies are born with. In line with our finding on higher BCAT2 expression with high birth weight, LBW pigs were found to express less BCAT2 mRNA in the longissimus dorsi muscle compared to normal birth weight pigs [27]. Similarly, higher BCAT2 mRNA was revealed in the blastocysts of diabetic rabbits compared to control blastocysts [28]. We observed an inverse association of IVD expression with birth weight. It has been demonstrated that the deficiency of the mitochondrial enzyme IVD may lead to isovaleric acidemia (IVA), an inherited metabolic disorder that may cause problems with the breakdown of the amino acid leucine [29]. Children with this condition may fail to gain weight and often experience developmental delays [30]. SLC25a16 has been considered as a carrier of Grave's disease, which causes hypothyroidism [31]. On the other hand, hypothyroidism is thought to cause HBW [32,33], which may explain the association we found between SLC25a16 and higher birth weight, but this needs to be further explored. Our results on gene expression suggest new insights into birth weight changes indirectly caused by MDCs, and also provide some support, albeit weak signals, for the existing evidence from transcriptomics-birth weight research. However, it is also important to note that none of the features selected for further analyses from the TWAS models passed the FDR correction threshold; our gene expression results should therefore be viewed as exploratory and hypothesis-generating.
Metabolism-related pathways linked to both an MDC and birth weight were not observed in this study. However, some results on the metabolism-related pathways associated with an MDC or birth weight are noteworthy. For PFOA, we have observed positive associations with the metabolism of xenobiotics by cytochrome P450 and drug metabolism and inverse association with T1D. In a mouse study, PFOA was found to induce the cytochrome P450 enzyme by activating constitutive androstane receptor (CAR) nuclear receptors [34]. Another mouse study has shown that PFOA may induce drug metabolism, and then lead to liver injury [35]. For PFOS, we observed inverse associations with amino sugar and nucleotide sugar metabolism, and T1D. Consistently, PFOS-induced altered amino sugar and nucleotide sugar metabolism were found in a recent zebrafish study, as well as in Hispanic children [36,37]. In a large U.S. study, PFOA and PFOS were associated with a reduced risk of T1D in adults [38], but in a recent Finnish study, they both were associated with an increased risk of T1D in newborns [39]. For birth weight, lower birth weight was found to be associated with six metabolism-related pathways, comprising OXPHOS, NAFLD, cysteine and methionine metabolism, sulfur metabolism, valine, leucine and isoleucine degradation, and fatty acid biosynthesis. Consistent with our findings, LBW was shown to be associated with OXPHOS in the skeletal muscle and myotubes of Danish individuals [40,41]. A study investigating the relationship between birth weight and NAFLD, in 538 children, also showed an overrepresentation of LBW in those with NAFLD compared with the general U.S. population [42]. This inverse relationship between birth weight and NAFLD occurrence was also confirmed in a large French prospective cohort study of 55,034 adults [43]. Likewise, a recent systematic review and meta-analysis demonstrated that excess methionine and cysteine led to lower birth weight [44]. The effect of branched chain amino acids (valine, leucine and isoleucine) on birth weight was not yet clear, and most of the existing studies were animal studies [44]. In addition, there is growing evidence that there may be an association between high fatty acid levels and LBW [45][46][47][48].
The strengths of our study include the well-defined sampling frame and the use of omics techniques, which allow for the investigation of multiple genes and pathways simultaneously, in order to explore the impact of MDCs on the transcriptome perturbations and the subsequent impact on the birth weight. We also acknowledge several limitations of this study. First, the relatively small sample size (n = 193 mother-child pairs) of our study population was prone to modest statistical power in detecting associations. Also for this reason, we did not perform sex-specific analysis despite that EDCs have been shown to exert different adverse effects in males and females, both in laboratory animals and in humans [49]. Second, it should also be noted that the concentrations of p,p'-DDE, PCB-153, and PFOS in our study population were relatively low compared with the median exposure levels observed in other studies that found associations with birth weight [7,50,51], and they may not have been high enough to have a measurable effect, or the limited contrast in exposures may have limited statistical power to detect associations; PFOA levels were more comparable with levels in other studies. Third, the cross-sectional design of the study precluded establishing a temporal or causal relationship between MDC concentrations, transcriptome, and birth weight. Last, as with any other observational epidemiological study, there may be residual confounding bias due to uncontrolled unmeasured confounders, but we expected these to be minimal, as we carefully adjusted for a set of covariates that have been shown to be important with the help of directed acyclic graphs (DAGs).
In addition, the mechanisms are complex and sensitive windows, for exposure to MDCs may vary depending on the specific chemical. Alterations at the molecular level caused by MDCs may also differ according to the specific outcome being studied. Therefore, different exposure windows and outcomes should be assessed in further studies investigating the metabolism-disrupting effects of chemicals.

Study Population
We used data from the second cycle of the Flemish Environment and Health Study (FLEHS II, 2008-2009), whose design and recruitment have been previously described in detail [52]. In short, 255 mother-child pairs were recruited from Flanders, Belgium, using a two-stage sampling procedure, with provinces as the primary sampling unit and maternity units as the secondary sampling unit. Mothers who had lived for at least 10 years in Flanders and were able to fill in Dutch questionnaires were invited to participate. The number of participants in each province was proportional to the number of inhabitants. Among the mother-child pairs, 195 were randomly selected for transcriptome profiling. We restricted our analyses to the 193 term births (gestational age ≥37 weeks) in this study because preterm birth is a potential mediator of the effects of chemical exposures on birth weight [53].

Exposure Assessment
Several classes of environmental chemicals were measured in cord blood samples. Here, we have focused our analyses on MDC exposures that could be detected in at least 60% of the cord blood samples [54]: p,p'-DDE, PCB-153, PFOA, and PFOS (see Supplementary Material, Table S1 for detection rates, which ranged from 97 to 100% for these four chemicals). Samples were collected immediately after birth and stored at −80 • C until the measurements. MDC concentrations were measured using gas chromatographyelectron capture negative ionization mass spectrometry (for p,p'-DDE and PCB-153) and high-performance liquid chromatography with tandem mass spectrometry detection (for PFOA and PFOS), as previously described [55,56]. All of the samples had quantifiable concentrations of p,p'-DDE, PFOA, and PFOS, while for PCB-153, 3% of the samples had values below the limit of quantification (LOQ, 300 ng/L). These values were then imputed using maximum likelihood estimation, assuming a censored log-normal distribution for values above the LOQ and conditional on the observed values for other biomarkers [54,57]. Lipidstandardized p,p'-DDE and PCB-153 concentrations were calculated based on estimated total lipids [total lipids = 50.49 + 1.32 × (cholesterol + triglycerides) (mg/dL)] and expressed as ng/g lipids for subsequent analyses [15]. All MDC concentrations were log 2 -transformed in order to reduce the potential influence of extreme values.

Transcriptome Profiling and Processing
As previously described [15,58], total RNA was extracted from the cord blood samples and stored at −80 • C. Amplified and labeled cRNA were then hybridized to 4 × 44 K Agilent Whole Human Genome Microarray (design 014850, one-color experimental setup with Cy3-labeling; Santa Clara, CA, USA), according to the manufacturer's protocol. Preprocessing, quality assessment, and normalization of the microarray data were performed as described previously [15]. Briefly, the arrays were scanned with an Agilent scanner (G2565BA) and were subjected to primary quality control using the Agilent Feature Extraction Software (Version 10.7; Santa Clara, CA, USA). Furthermore, for each feature on the array, the quantile-normalized and log2-transformed signal intensity derived from Cy3 fluorescent dye was used for subsequent analyses. For replicated features on the array, the mean of signal intensities was calculated. After control and noise filtering by removing features with signal intensity below 3, 33,543 features retained. Thereafter, we used the R package Combat to eliminate possible batch effects related to different hybridization dates (28 dates from 14 September 2011 to 11 January 2012) [59,60]. Lastly, 26,170 (78.02%) features were annotated to a total of 17,880 unique gene symbols according to the Molecular Signatures Database (MSigDB) and were subjected to further statistical and functional analyses [61].

Outcome Assessment and Covariates
We considered birth weight (g) as our outcome of interest. DAGs were used to guide the selection of covariates (Figures S1b and S2a). The set of minimally sufficient covariates included sex of the child (girl, boy), smoking during pregnancy (smoking, non-smoking), parity (0, 1, ≥2), maternal education (low, medium, high), maternal age at delivery (<27, 27 < 30, 30 < 33, ≥33 years), pre-pregnancy BMI (<18.5, 18.5 < 25, 25 < 30, ≥30 kg/m 2 ), and gestational age (weeks). Birth weight and child sex were collected from maternity medical records. Other covariate data was obtained from questionnaires. Missing data in covariates and exposures that were completely missing (1-3% and 1% of participants had one of more missing values, respectively) were singly imputed using the R package mice [62].

TWAS
TWASs were conducted in order to investigate the association of global transcriptomics with (1) MDCs and (2) birth weight. We used the following multivariable linear models to evaluate the effects of MDC exposures and potential predictors of birth weight, for each feature and MDC separately: log 2 (feature intensity i ) = β 0 + β 1 log 2 (MDC i ) + β 2 sex i + β 3 smoking during pregnancy i + β 4 parity i + β 5 education i + β 6 age at delivery i + β 7 pre-pregnancy BMI i + β 8 gestational age i + ε 1i (1) birth weight i = γ 0 + γ 1 log 2 (feature intensity i ) + γ 2 sex i + γ 3 smoking during pregnancy i + γ 4 parity i + γ 5 education i + γ 6 age at delivery i + γ 7 pre-pregnancy BMI i + γ 8 gestational age i + ε 2i (2) where i indexes the study subjects and Model (1) describes the association between a single transcriptomic feature and a single MDC, while Model (2) describes the association between birth weight and a single transcriptomic feature. Parameters β 0 and γ 0 are the model intercepts, while β 1 and γ 1 refer to the effect estimates (slopes) for a single MDC on a single transcriptomic feature, and for a single transcriptomic feature on birth weight, respectively. Parameters β 2-8 and γ [2][3][4][5][6][7][8] are coefficients corresponding to other covariates in the model, and ε 1i and ε 2i represent the residual errors, which are assumed to follow a normal distribution. According to observed p-values for β 1 and γ 1 , we estimated FDR using the method of Benjamini and Hochberg to correct for multiple testing and to select significant features [63].

Enrichment Pathway Analysis
In order to find pathways associated with MDC exposures and birth weight, we carried out Gene Set Enrichment Analyses (GSEA) using the WEB-based GEne SeT AnaLysis Toolkit (WebGestalt; Los Angeles, CA, USA) tool with pathway gene sets from the KEGG database [22,64,65]. First, we generated the respective ranked lists of all 26,170 features, sorted by their degree of differential expression (log 2 -fold change) in cord blood in relation to MDCs and birth weight, i.e., β 1 and γ 1 obtained from Models (1) and (2) [66,67]. Subsequently, the normalized enrichment scores were calculated, reflecting the degree to which pathways were enriched by ranked genes, where positive and negative values represent positive and inverse associations of pathways with MDCs or birth weight, respectively [68]. We restricted to pathways with at least five genes involved, and estimated the statistical significance using 1000 gene set permutations with FDR correction for multiple testing. Pathways with FDR < 0.05 were considered significant. Figure 2 outlines the workflow of the meet-in-the-middle approach used in this study [69]. The overlapping selected features and pathways observed in association with any of the four MDCs and birth weight were further explored by mediation analysis using the R package mediation [70] to explore potential biological mechanisms and mediating effects linking exposure and outcome. When assessing an overlapping feature as a mediator, we included it in the mediation model, and computed ACMEs (also known as indirect effects) using 1000 bootstrapped samples with FDR correction, and it was considered as a potential mediating feature if the FDR < 0.05. When examining an overlapping pathway as a mediator, we first performed a principal component analysis (PCA) on the genes belonging to that pathway, and then used the first principal component score (PC1) to represent that pathway in the mediation model [71,72]. ACMEs with FDR < 0.05 were generated to identify potential mediating pathways.

Sensitivity Analysis
Recognizing that gestational age could be associated with MDC exposure and impact transcriptome levels [15], combined with several other studies also showing that the transcriptome was substantially influenced by gestational age [73,74], gestational age was included as a control variable in our primary regression model of MDCs and transcriptome (TWAS Model (1)). On the other hand, the causal direction of the association between gestational age and MDC is not entirely clear, and it is possible that gestational age mediates the outcome [75]. Therefore, in a sensitivity analysis, we assessed MDC and transcriptome associations without adjusting for gestational age in order to avoid adjustment for a potential mediator [53].
All statistical analyses were performed in R version 4.1.0 [76].

Conclusions
In summary, we integrated cord blood TWASs in order to identify gene expressions and pathways associated with MDCs and birth weight. Taken together, our study suggested five gene expressions associated with at least one MDC and birth weight. This provides insight into the etiology of higher and lower birth weight and possible later metabolic disorders, but again, this is an exploratory study with weak signals. In order to validate our results and further understand the potential link between MDC exposures and birth weight, and to elucidate the underlying mechanisms, studies with larger sample sizes and prospective study designs combined with advanced omics techniques are warranted.