Expression Changes in Epigenetic Gene Pathways Associated With One-Carbon Nutritional Metabolites in Maternal Blood From Pregnancies Resulting in Autism and Non-Typical Neurodevelopment.

The prenatal period is a critical window for the development of autism spectrum disorder (ASD). The relationship between prenatal nutrients and gestational gene expression in mothers of children later diagnosed with ASD or non-typical development (Non-TD) is poorly understood. Maternal blood collected prospectively during pregnancy provides insights into the effects of nutrition, particularly one-carbon metabolites, on gene pathways and neurodevelopment. Genome-wide transcriptomes were measured with microarrays in 300 maternal blood samples in Markers of Autism Risk in Babies-Learning Early Signs. Sixteen different one-carbon metabolites, including folic acid, betaine, 5'-methyltretrahydrofolate (5-MeTHF), and dimethylglycine (DMG) were measured. Differential expression analysis and weighted gene correlation network analysis (WGCNA) were used to compare gene expression between children later diagnosed as typical development (TD), Non-TD and ASD, and to one-carbon metabolites. Using differential gene expression analysis, six transcripts (TGR-AS1, SQSTM1, HLA-C, and RFESD) were associated with child outcomes (ASD, Non-TD, and TD) with genome-wide significance. Genes nominally differentially expressed between ASD and TD significantly overlapped with seven high confidence ASD genes. WGCNA identified co-expressed gene modules significantly correlated with 5-MeTHF, folic acid, DMG, and betaine. A module enriched in DNA methylation functions showed a suggestive protective association with folic acid/5-MeTHF concentrations and ASD risk. Maternal plasma betaine and DMG concentrations were associated with a block of co-expressed genes enriched for adaptive immune, histone modification, and RNA processing functions. These results suggest that the prenatal maternal blood transcriptome is a sensitive indicator of gestational one-carbon metabolite status and changes relevant to children's later neurodevelopmental outcomes. LAY
SUMMARY: Pregnancy is a time when maternal nutrition could interact with genetic risk for autism spectrum disorder. Blood samples collected during pregnancy from mothers who had a prior child with autism were examined for gene expression and nutrient metabolites, then compared to the diagnosis of the child at age three. Expression differences in gene pathways related to the immune system and gene regulation were observed for pregnancies of children with autism and non-typical neurodevelopment and were associated with maternal nutrients.

Introduction

Autism spectrum disorder (ASD) is a group of neurodevelopmental disorders characterized by persistent impairment in social interactions, communication, restricted interests, or repetitive behaviors, and language deficits [Maenner et al., 2020]. Current data show that one in every 54 children in the United States has ASD [Maenner et al., 2020]. One major component of ASD risk is genetic heritability, based on studies of twins, siblings, and other family members [Ozonoff et al., 2011; Sandin et al., 2014; Wessels & Pompe van Meerdervoort, 1979]. Common genetic variants each having small effects dominate most ASD risk compared with rare gene variants with large effects [Weiner et al., 2017]. Large genome‐wide association studies (GWAS) support the role of common genetic variants in ASD with remaining challenges in ASD complexity and heterogeneity [Grove et al., 2019; Iossifov et al., 2014; Sanders et al., 2015]. Mutations in single genes can explain less than 1% of ASD cases [Bourgeron, 2015; Tsai & Bell, 2015].

Accumulating lines of evidence suggest that ASD risk arises from both genetic and environmental risk factors. In utero maternal exposures can contribute as ASD risk factors, including air pollution, fever, asthma, and nutrition, especially nutrients involved in the one‐carbon metabolic pathway [Raz et al., 2015; Schmidt et al., 2011, 2012; Zerbo et al., 2013]. Other studies suggest that one‐carbon metabolism is implicated in gene–environment interactions in ASD [L. Schaevitz, Berger‐Sweeney, & Ricceri, 2014; L. R. Schaevitz & Berger‐Sweeney, 2012]. Maternal prenatal nutritional supplements containing folic acid and additional B vitamins that play a role in one‐carbon metabolism are associated with ASD risk reduction [Schmidt et al., 2011; Schmidt, Iosif, Guerrero Angel, & Ozonoff, 2019; Suren et al., 2013]. A common genetic polymorphism affecting folic acid metabolism, MTHFR C677T, interacts with maternal nutrition, as association between folic acid and reduced ASD risk was strongest for mothers and children with MTHFR 677 C > T variant genotypes [Y. Li et al., 2020; Schmidt et al., 2012]. These findings suggest that additional gene–environment interactions relevant to ASD may be identified from investigations into maternal factors, since maternal and fetal metabolisms are shaped by both shared genetics and nutritional environment during pregnancy that can coordinately impact neurodevelopment.

Gene expression levels are also influenced by both genetic and environmental factors, especially by in utero maternal nutrition [Vucetic, Kimmel, Totoki, Hollenbeck, & Reyes, 2010; Yajnik & Deshmukh, 2012]. Maternal peripheral blood therefore offers a unique window to study transcriptome alterations during pregnancy that may reflect altered fetal development associated with nutrition [Costello et al., 2008; Croen et al., 2008]. Numerous environmental factors during pregnancy can alter gene expression levels [Haugen, Schug, Collman, & Heindel, 2015; Zerbo et al., 2013]. Other neurodevelopment disorders, such as schizophrenia, have also demonstrated a significant interaction of genetic risk with maternal perinatal environmental factors that affected the transcriptome [Ursini et al., 2018; Xu et al., 2014]. Postmortem brain gene expression studies revealed gene co‐expression modules enriched for immune response and neuronal development functions in ASD [Gupta et al., 2014; Voineagu et al., 2011]. Other studies using child peripheral blood and cord blood showed that differential gene expression in ASD was enriched for immune and inflammatory processes [Ansel, Rosenzweig, Zisman, Melamed, & Gesundheit, 2017; Mordaunt, Park, et al., 2019; Tylee et al., 2017].

While numerous studies have investigated specific genes or pathways in children with ASD, none have focused on the maternal transcriptome during pregnancy. Further, most previous ASD transcriptome studies used data from specimens collected postmortem or after childbirth, as opposed to prospective studies to help understand potential etiologic changes that occur before behavioral symptoms. Other large epidemiology studies examined environmental effects in ASD, but how the environment influences alterations at the molecular level remains to be understood. The goal of this study was to examine maternal prenatal gene expression profiles associated with both maternal serum one‐carbon metabolites and the child outcome (ASD, Non‐TD, TD) to shed light on molecular changes during pregnancy.

Methods

Study Design

The MARBLES study recruited mothers in Northern California with at least one child with ASD who were pregnant or planning another pregnancy. Due to a shared genetic background, 24% of the next children met the criteria for ASD. A previous publication detailed the study design of MARBLES [Hertz‐Picciotto et al., 2018; Mordaunt, Park, et al., 2019]. In order to enroll into MARBLES, all five of the following criteria needed to be met: (a) the prospective child has one or more first or second degree relatives with ASD; (b) mother is 18 years or older; (c) mother is pregnant or able to become pregnant; (d) mother is able to speak, read, and understand English and plans to raise the child with English spoken at home; (e) mother lives within a 2.5‐hr drive from Davis/Sacramento, California. Demographic, diet, environmental, and medical information were collected by telephone interviews or questionnaires throughout the pregnancy. Infants received standardized neurodevelopmental assessments from 6 months until 3 years old [Hertz‐Picciotto et al., 2018]. At 3 years old, the child was assessed clinically using the gold standard Autism Diagnostic Observation Schedule (ADOS) [Lord et al., 2000], the Autism Diagnostic Interview—Revised (ADI‐R) [Rutter, LeCouteur, & Lord, 2015], and the Mullen Scales of Early Learning (MSEL) [Mullen, 1995]. Based on a previously published algorithm using ADOS and MSEL scores [Mordaunt, Park, et al., 2019; Schmidt et al., 2019], participants were classified into three outcome groups including ASD, TD, and Non‐TD [Chawarska et al., 2014; Ozonoff et al., 2014]. Children with ASD had scores over the ADOS cutoff and fit ASD DSM‐5 criteria. Children with Non‐TD outcomes were defined as children with low MSEL scores (two or more MSEL subscales with more than 1.5 standard deviations (SD) below averages or at least one MSEL subscale more than 2 SD below average) and elevated ADOS scores. Children with TD outcome had all MSEL scores within 2.0 SD and no more than one MSEL subscale that is 1.5 SD below the normative mean and scores on the ADOS at least three points lower than the ASD cutoff.

Isolation and Expression Microarray

Maternal peripheral blood was collected at study visits during all three trimesters of pregnancy in PAXgene Blood RNA tubes with the RNA stabilization reagent (BD Biosciences) and stored frozen at −80°C. The first timepoint sample was used for mothers who had multiple blood draws (n = 12) during pregnancy. RNA was isolated using the PAXgene Blood RNA Kit (Qiagen) according to the default protocol. Total RNA was converted to cDNA and biotin labeled. Expression was measured using Human Gene 2.0 Affymetrix microarray chips with three batches by the John Hopkins Sequencing and Microarray core following washing, staining, and scanning procedures based on manufacturer's protocol.

Data Preprocessing and Normalization

Robust Multi‐Chip Average (RMA) [Bolstad, Irizarry, Astrand, & Speed, 2003; R. A. Irizarry, Hobbs, et al., 2003; Rafael A. Irizarry, Bolstad, et al., 2003] from the oligo R package was used for normalization of Affymetrix CEL files. For quality control, we used the oligo and ArrayQualityMetrics R packages [Carvalho & Irizarry, 2010; Kauffmann, Gentleman, & Huber, 2009]. No samples were identified as outliers by principal component analysis, the Kolmogorov–Smirnov test, or Euclidean distance to other arrays. Probes were mapped at the transcript level using the pd.hugene.2.0.st R package, and those annotated to genes (36,459) were used in subsequent analyses.

One‐carbon Nutrient Metabolite Measurements

Serum and plasma samples from the same blood draw as specimens used for RNA expression analysis were used to measure one‐carbon and nutrient metabolites. Metabolites were selected for their role in one‐carbon metabolism and may be relevant to gene expression. S‐adenosylmethionine (SAM) and S‐adenosylhomocysteine (SAH), adenosine, homocysteine, cystine, cysteine, glutathione (GSH), and glutathione disulfide (GSSG) were measured in the James' laboratory at the Arkansas Children's Research Institute using HPLC with electrochemical detection as previously described [Hollowood et al., 2018; Melnyk et al., 2012]. Serum pyridoxal phosphate (PLP), the biologically active form of vitamin B6 (Vit B6), was measured by HPLC using fluorescence detection in the Green‐Miller laboratory at the UC Davis Medical Center (inter‐assay coefficient of variation (CV) = 4.8%) [Talwar et al., 2003]. Total serum vitamin B12 (Vit B12) was measured using automated chemiluminescence in the CLIA‐approved Medicine Clinical Laboratories at UC Davis Medical Center (inter‐assay CV = 6.2%). Plasma choline, betaine, and dimethylglycine (DMG) were measured using LC–MS/MS stable isotope dilution methods in the Caudill laboratory [Holm, Ueland, Kvalheim, & Lien, 2003; Yan et al., 2011] with modifications to include measurements of trimethylamine N‐oxide (TMAO) and methionine [Wang et al., 2008; Yan et al., 2012]. Intra‐assay and inter‐assay CVs of the in‐house controls were 3.0% and 3.6% for choline; 1.5% and 1.7% for betaine, 2.5% and 2.4% for DMG; 2.6% and 2.6% for methionine; and 3.1% and 3.4% for TMAO. Serum 5‐methyltetrahydrolate (5‐MeTHF) and folic acid were quantified in the Caudill laboratory using LC–MS/MS stable‐isotope dilution methods [Pfeiffer, Fazili, McCoy, Zhang, & Gunter, 2004] with modifications based on the instrumentation [West et al., 2012]. Intra‐assay and inter‐assay CVs of in ‐house controls were 1.8% and 1.9% for 5‐methyltetrahydrofolate; and 4.9% and 8.5% for folic acid.

Differential Gene Expression

After normalization, surrogate variable analysis (SVA) was used to estimate and adjust for hidden confounding variables on gene expression [Leek, Johnson, Parker, Jaffe, & Storey, 2012]. Differential gene expression was identified using the limma R package by a linear model that included the children's diagnosis outcome (ASD, Non‐TD, and TD) and all surrogate variables [Ritchie et al., 2015]. Differential gene expression analysis with children diagnosed as ASD, Non‐TD, and TD were included in the same model with three levels of diagnosis using the F‐test [Ritchie et al., 2015]. Pairwise fold change, standard error, and P‐value between ASD versus TD, and Non‐TD versus TD were extracted from the same model using the limma R package. Differentially expressed transcripts were identified as those with an unadjusted P‐value <0.05. Genome‐wide significant differentially expressed transcripts were classified as those with a false discovery rate (FDR) adjusted P‐value (q‐value) <0.05. To complement adjustment by SVA, we used an alternative approach that adjusted for known confounders, including batches, trimesters, child gender, and cell types using limFit function in limma R package. Differentially expressed transcripts between ASD versus TD, and non‐TD versus TD were extracted using the same approach as SVA adjustment and compared.

Gene Overlap Analysis

Gene overlap analysis was performed by Fisher's exact test using the GeneOverlap R package [Li, 2019]. In each comparison, the null distribution was generated from 1000 random samples of all genes annotated to transcripts on the array. Gene symbols annotated to differentially expressed transcripts were compared to 943 genes in the Simons Foundation Autism Research Initiative (SFARI) Gene database [Abrahams et al., 2013]. As an alternative approach to using the whole probe set on the Affymetrix arrays (36,459) as background, an analysis of filtered probe values was performed based on removing probe intensity on the lowest 5% intensity probes (34,636) to remove genes with the lowest expression levels.

Gene Ontology (

Transcripts with significant expression levels or selected gene lists were exported to DAVID bioinformatics software with default settings for GO analysis [Ashburner et al., 2000; The Gene Ontology Consortium, 2017]. The analysis was done using the GO ontology database and Fisher's exact test with multiple test correction by the FDR method [The Gene Ontology Consortium, 2017]. GO term enrichments were presented with hierarchical terms. GO terms with an FDR adjusted P‐value <0.05 were considered statistically significant. GO terms enrichment analyses was performed separately for both all genes and filtered probes representing expressed genes.

Weighted Gene Network Analysis (

A weighted gene co‐expression network was built using the WGCNA R package [Langfelder & Horvath, 2008; Zhang & Horvath, 2005] with normalized expression levels after adjustment for batch effects (3 batches) using the ComBat function from the sva R package [Johnson, Li, & Rabinovic, 2007]. The correlation matrix included all probes and all samples. To construct a signed adjacency matrix, estimated soft thresholding power (power = 6) was used to achieve approximately scale‐free topology (R^2 fit >0.8). Adjacency values were transformed into a signed topological overlap matrix (TOM). Co‐expression modules were identified from the dissimilarity matrix (1‐TOM) with a minimum module size of 30 probes using Pearson's coefficient. Module eigengenes were clustered based on correlation. Similar modules were merged based on a cut height of 0.25 to generate co‐expression modules. Each module's expression profile was summarized into a module eigengene (ME) using the matched module's first principal component. The correlation between each gene in the module with the ME was represented as intramodule connectivity (kME). Module hub probes were defined as the probe in each module with the highest module membership. Hierarchical clustering was done using the standard R function hclust with the default setting using ward's agglomeration method [Gentle, Kaufman, & Rousseuw, 1991]. Pearson's correlation coefficient was used to measure the correlation between traits and modules. Highly correlated modules were defined as those with an FDR adjusted P‐value <0.05.

Cell Type Proportion Deconvolution

CIBERSORT was used to estimate the proportions of each cell type using the default settings and the LM22 adult peripheral blood signature gene expression profiles [Newman et al., 2015]. Normalized expression levels adjusted for batch effects were used to estimate cell type proportions. Both relative and absolute modes were performed together with 100 permutation tests. P‐values were calculated using FDR multiple test adjustment. Significant associations were defined based on FDR adjusted P‐value <0.05.

Results

Study Sample Characteristics and Nutrient Measurements

High quality RNA was isolated from 300 maternal peripheral blood samples collected during pregnancy from the MARBLES high risk ASD cohort (Table S1). Children from MARBLES pregnancies were diagnosed at 3 years old as ASD (67, including 47 male and 20 female), Non‐TD (79, including 46 male and 33 female), and TD (154, including 79 male and 75 female) (Table S2).

Nutrients in the one‐carbon metabolism pathway, including methionine, SAM, SAH, adenosine, homocysteine, 5‐MeTHF, folic acid, Vit B6, Vit B12, choline, DMG, betaine, cystine, cysteine, GSH, and GSSG were directly measured from maternal blood in 14%–62% of all samples (Table S3). None of these metabolites in maternal blood were significantly associated with clinical outcomes of children (Table S3). Measurements for one‐carbon metabolites and transcriptomes were conducted on samples collected throughout pregnancy (Fig. S1).

Differential Gene Expression Analyses by Child Outcome

Expression was measured using the Human Gene 2.0 Affymetrix microarray and adjusted for all surrogate variables, followed by differential gene expression analysis for child diagnosis (ASD, Non‐TD, TD) on 36,459 transcripts. There were 28 surrogate variables (SVs) identified, including 5 SVs significantly associated with batch effect and 2 SVs significantly associated with gestational age of maternal blood draw (Fig. S2). Six transcripts located at four genes (TGR‐AS1, SQSTM1, HLA‐C, and RFESD) were associated with child outcomes (ASD, Non‐TD, TD) with genome‐wide significance by F‐test (FDR adjusted P‐value <0.05) (Table S4). Three out of these six transcripts mapped to HLA‐C (Major Histocompatibility Complex, Class I, C) (FDR adjusted P‐value <0.05).

Comparing the maternal blood transcriptome between ASD and TD outcomes revealed 2,012 differentially expressed transcripts at a nominal confidence level (unadjusted P‐value <0.05) that mapped to 1,912 genes, including 980 up‐regulated and 1,032 down‐regulated transcripts, with none significant after FDR adjustment (Fig. 1A, Table S5). There was a significant overlap between these 1,912 differentially expressed genes and a list of strong ASD candidate genes from the Simons Foundation Autism Research Initiative (SFARI Gene, including TRIO, GRIA1, SMARCC2, SPAST, DIP2C, FOXP1, and CNTN4, Fisher's exact test, P‐value <0.05) [Abrahams et al., 2013].

Figure 1

Identification and function of ASD associated and Non‐TD associated differentially expressed genes in maternal peripheral blood. Differential expression analysis was performed in maternal peripheral blood transcriptomes (n = 300) after adjustment for surrogate variables. (A) Identification of 1,912 differentially expressed genes (2,012 transcripts, P‐value <0.05) compared between children diagnosed as ASD (n = 67) and TD (n = 154). (B) Identification of 1,919 differential expressed genes (1,987 transcripts, P‐value <0.05) compared between children diagnosed as Non‐TD (n = 79) and TD (n = 154). Two transcripts located at RFESD and TRG‐AS1 were genome‐wide significant in the Non‐TD to TD comparison (Table S2). (C) Venn diagram represents the overlap in differentially expressed transcripts (unadjusted P <0.05) identified in ASD to TD versus Non‐TD to TD comparisons, which was greater than expected by random using a Fisher's exact test (P‐value <0.001***). (D) Gene ontology (GO) and pathway analysis was performed on the 218 transcripts differentially expressed in both ASD‐TD and Non‐TD‐TD comparisons, with significant enrichments (Fisher's exact test, FDR P‐value <0.05). In contrast, the differentially expressed transcripts uniquely associating with either ASD or Non‐TD were not significantly enriched for any GO terms.

Comparing the maternal blood transcriptome between Non‐TD and TD outcomes revealed 1,987 differentially expressed transcripts at a nominal confidence level (unadjusted P‐value <0.05) that mapped to 1,919 genes, including 1,044 up‐regulated and 943 down‐regulated transcripts (Fig. 1B, Table S6). Two of these transcripts, RFESD and TRG‐AS1, also passed genome‐wide significance (FDR adjusted P‐value <0.05). Unlike the ASD versus TD comparison, however, no significant overlap was observed between Non‐TD versus TD differentially expressed genes and SFARI gene lists. An alternative approach of adjustment for known confounders revealed a significant overlap with SVA adjusted differentially expressed transcripts for both ASD versus TD and Non‐TD versus TD comparisons (Fisher's exact test, P‐value <2.2E‐16) (Fig. S3, Table S4).

Differential gene expression analysis using SVA resulted in a significant overlap of 218 transcripts between ASD versus TD differentially expressed transcripts and Non‐TD versus TD differentially expressed transcripts (Fisher's exact test, P‐value <2.2E‐16) (Fig. 1C). Gene ontology (GO) analysis of these 218 transcripts revealed significant enrichment for the interferon‐gamma mediated signaling pathway, apoptosis in muscle, response to interferon gamma, and metal ion transport (Fig. 1D; Fig. S4). CaMK (calmodulin‐dependent protein kinase) families (CAMK2A, CAMK2B, CAMK2D, and CAMK2G) and HLA (human leukocyte antigen) systems (HLA‐B, HLA‐C, and HLA‐E) were included in those significant signaling pathways (Fig. 1D). In contrast, neither list of ASD‐ or Non‐TD‐specific differentially expressed transcripts were significantly enriched for any GO terms. To control for potential bias towards expressed genes in the GO terms, we also used an alternative approach of filtering the probes based on intensity that removed the lowest 5% of expressed transcripts (Fig. S5), which resulted in an identical list of significant GO terms (Table S7).

Weighted Gene Co‐expression Network Analysis (

WGCNA was performed as a complementary bioinformatic approach that incorporates the independent and inter‐related associations of transcript levels with measured concentrations of maternal nutrients. First, expression values were adjusted for batch effects, then correlation patterns among all transcripts were analyzed across all 300 samples. WGCNA identified 27 co‐expressed gene modules in our dataset, representing 17,049 transcripts, distinguished from 19,410 transcripts without evidence of co‐expression were grouped into the “gray” module (Fig. 2A,B; Fig. S6, Table S8). For each module, the number of transcripts, as well as the hub gene, defined as the gene with the highest correlation with the module eigengene, were determined (Fig. 2B; Table S9). Out of those 27 co‐expression modules, 23 modules showed associations between eigengene expression level and at least one variable related to demographics, diagnosis, or maternal nutrients, after FDR correction (FDR adjusted P‐value <0.05) (Fig. 2A; Fig. S6). All 27 modules were significantly associated with one or more traits, including child clinical outcome, demographic factors, and maternal blood metabolite concentrations at unadjusted P‐value <0.05 (Fig. 2A; Fig. S7).

Figure 2

Co‐expression network modules with demographic factors and maternal peripheral blood one‐carbon metabolites. (A) Heatmap of Z‐scores of modules eigengenes with sample covariates with 27 co‐expression network modules on all 300 maternal blood samples. Each row represents a different module eigengene and each column is the associated trait, which include child clinical outcome, demographic factors, and maternal blood metabolite concentrations. The matrix was calculated by Pearson correlation and P‐values adjusted for the total number of comparisons. Color represents the direction (red, positive correlation; blue, negative correlation) and intensity reflects the significance. (^ unadjusted P‐value <0.05 and FDR adjusted P‐value >0.05; * FDR adjusted P‐value <0.05). (B) Number of transcripts and hub genes from all 27 co‐expressed modules are listed.

Multiple co‐expression modules were significantly correlated (FDR adjusted P‐value <0.05) with gestational age at blood draw and four maternal metabolites, including 5‐MeTHF, folic acid, DMG, and betaine (Fig. 2A). None of the additional measured variables was significantly associated with any co‐expression gene modules, including clinical outcome. However, the module “greenyellow” showed a nominally significant positive correlation with outcome (unadjusted P‐value = 0.02, FDR adjusted P‐value = 0.14) and a negative significant correlation with both 5‐MeTHF (FDR adjusted P‐value = 0.02) and folic acid levels (FDR adjusted P‐value = 0.02) (Figs. 2A and 3A,B; Figs. S6 and S7; Table S10). Interestingly, the greenyellow module eigengene was correlated in opposite directions with ASD and 5‐MeTHF, consistent with a putative 5‐MeTHF protective effect in ASD (Fig. 3A,B).

Figure 3

“Greenyellow” module was positively associated with diagnosis and negatively associated with folic acid and 5‐MeTHF. (A) “Greenyellow” module eigengene was significantly associated with child diagnosis (one‐way ANOVA, unadjusted P‐value <0.05). Greenyellow eigengene values were higher in maternal blood from ASD pregnancies than TD or Non‐TD pregnancies. (B) “Greenyellow” module eigengene level was significantly negatively associated with 5‐MeTHF concentrations in maternal blood (ANOVA, P‐value <0.05). (C) Bar graph shows gene ontology (GO) and pathway significant enrichments from the 224 transcripts in “greenyellow” module (Table S8). (D) Transcripts (16857547, 16867905, and 16867910) from MBD3L3‐5 genes encoding proteins involved in methylation‐CpG binding functions were significantly negatively associated with “greenyellow” module eigengene.

This “greenyellow” module contained 224 transcripts with TRNAI2 as the hub gene (Table S10). These 224 transcripts showed a significant enrichment for gene ontology functions in methylation‐CpG binding, methyl‐dependent chromatin silencing, and keratinocyte differentiation (Fisher's exact test, FDR adjusted P‐value <0.05) (Fig. 3C, Table S11). The three known genes with methyl‐binding functions included MBD3L3, MBD3L4, and MBD3L5, represented by 16857547, 16867905, and 16867910 transcripts (Fig. 3C). Normalized expression of those three transcripts was also significantly associated with the “greenyellow” module eigengene, supporting their membership in the module (Fig. 3D).

There were 12 modules in the maternal transcriptome significantly associated with gestational age, including 8 positively correlated modules (darkred, lightgreen, cyan, darkgrey, brown, magenta, orange, and green) and 4 negatively correlated modules (lightyellow, lightcyan, red, tan) (Fig. 2A, Table S12). Genes within modules with a positive correlation with gestational age were significant enriched for functions in RNA binding, chromatin binding, and ATP binding, among others (Table S12). In contrast, genes within modules negatively correlated with gestational age were significantly enriched in functions related to granulocyte activation, neutrophil mediated immunity, coagulation, and other blood functions (Table S12).

Eight Co‐expression Modules Strongly Clustered with Betaine and

Among the 27 identified co‐expression modules, eight modules (darkred, lightgreen, cyan, darkgrey, brown, magenta, orange, and white) were highly correlated with each other and clustered based on unsupervised hierarchical clustering, representing a total of 2,582 transcripts (Fig. 4; Fig. S8, Table S13). Betaine and DMG were significantly associated and clustered together with this distinct block of co‐expression modules.

Figure 4

Eight weighted gene co‐expression modules associated with maternal betaine and DMG concentrations were strongly clustered. (A) Unsupervised hierarchical clustering dendrogram was performed with module eigengenes, betaine and DMG. The height of each node represents the intergroup dissimilarity. Similar nodes clustered together under one branch. (B) Unsupervised hierarchical clustering adjacency heatmap, with color and intensity representing the degree of correlation (dark, high; light, low correlation). Black box indicates the block of eight weighted gene co‐expression modules associated with betaine and DMG concentrations.

Transcripts inside these eight clustered co‐expression modules associated with betaine and DMG showed significant enrichment for 18 gene pathways involved in adaptive immune response, RNA processing, histone modification, inflammatory response, and Rett syndrome (Fisher's exact test, FDR adjusted P‐value <0.05) (Fig. S9). Network analysis using GeneMANIA [Warde‐Farley et al., 2010] identified a network with EVL in the center, linked with other hub genes (Fig. S10).

Cell Type Composition in Maternal Peripheral Blood was Associated With Maternal Metabolites But Not Child Clinical Outcomes

In order to determine the effects of cell composition differences on the findings associated with maternal transcriptomes, cell type specific information from 22 immune cell types was deconvoluted using the CIBERSORT web tool. Maternal peripheral blood samples reflected a mixture of cell types, with neutrophils as the largest and most variable population ranging from 17% to 48% (Fig. 5A; Table S14). The eigengenes for 21 out of 27 modules were significantly correlated with at least one cell type (FDR adjusted P‐value <0.05) (Fig. S11). No significant difference was observed in cell type composition between child diagnosis outcomes or gender (Fig. 5A,B; Fig. S12). Furthermore, neither the “greenyellow” module, nor the betaine and DMG variables were significantly associated with cell type proportions, suggesting that the associations identified with these modules were largely cell type independent (Fig. 5B; Figs. S11 and S12).

Figure 5

Imputed cell type proportions in material peripheral blood associated with demographic factors and maternal nutrients. (A) Barplot of each cell type mean estimated proportion separated by children diagnosis outcomes using peripheral blood reference panel in CIBERSORT. (B) Heatmap of correlation between sample demographic factors and maternal nutrients with cell type proportions. Each row represents a cell type proportion and columns represent traits, including child diagnostic outcome, demographic factors, and maternal blood nutrient concentrations. P‐values adjusted for the total number of comparisons. Color represents the direction (red, positive correlation; blue, negative correlation) and intensity reflects the significance, *P‐value <0.05 after FDR correction.

In contrast, some cell type proportions were significantly correlated with some maternal metabolites. Vit B6, 5‐MeTHF, choline, cysteine, the ratio of DMG/betaine, and the ratio of cystine/cysteine were separately associated with six cell types (FDR adjusted P‐value <0.05) (Fig. 5B; Fig. S12). Vit B12, folic acid, the ratio of DMG/betaine, and the ratio of SAM/SAH were associated with more than one cell type (FDR adjusted P‐value <0.05) (Fig. 5B; Fig. S12). The most significant association was between vit B12 and memory B cells (FDR adjusted P‐value = 0.0001) (Fig. 5B; Fig. S12).

Discussion

Maternal blood collected during pregnancy can provide molecular insights into the in utero environment relevant to the etiology of ASD. This was the first study to our knowledge to examine gene expression differences together with one‐carbon metabolites in peripheral blood during pregnancy from mothers of children that went on to develop ASD, Non‐TD, or TD at 36 months. Using complementary bioinformatics approaches, we identified several genes and gene pathways consistent with proinflammatory and oxidative stress responses in mothers of children with adverse neurodevelopmental outcomes. We also identified eight novel co‐regulated gene modules associated with maternal blood betaine and DMG concentration.

Genes and Gene Patterns Common to Mothers of Children with and Non‐typical Neurodevelopment

Using differential gene expression analysis of individual genes, we describe four genes (SQSTM1, HLA‐C, TRG‐AS1, and RFESD) that were differentially expressed in mothers by child diagnosis outcome of either TD, ASD, or Non‐TD. SQSTM1 encodes the p62 sequestosome that acts as a receptor for ubiquitinated cargo in the selective autophagy response induced by oxidative stress [Sánchez‐Martín & Komatsu, 2020], and also links the mTOR and GABA signaling pathways in brain [Hui & Tanaka, 2019]. RFESD, encoding an iron–sulfur cluster binding protein with oxidoreductase activity, is located on 5q15, a hotspot for copy number variants in intellectual and developmental disabilities [Kaminsky et al., 2011; Sajan et al., 2013]. TRG‐AS1, T‐cell receptor gamma locus antisense RNA 1, is located on 7p14.1, another locus previously associated with developmental delay, intellectual disability, and ASD [Kaminsky et al., 2011; Klamt et al., 2016; Wenger et al., 2016]. HLA‐C belongs to the HLA (human leukocyte antigen) polymorphic loci encoding major histocompatibility class I (MHC I) proteins involved in antigen presentation to CD8+ T cells and NK cells. HLA‐C is important for both tolerance to fetal allo‐antigens and viral immunity during pregnancy [Papúchová, Meissner, Li, Strominger, & Tilburgs, 2019]. Proinflammatory cytokines such as interferon gamma (IFNγ) induce HLA‐C expression in both lymphocytes and placental trophoblasts. In the WGCNA modules, HLA‐C, SQSTM1, and TGR‐AS1 were all members of the same module, as evidence of co‐expression. A number of previous studies have shown that the HLA locus is associated with ASD [Saresella et al., 2009; Torres, Maciulis, Stubbs, Cutler, & Odell, 2002; Torres et al., 2006] or HLA locus activation in ASD children and their mothers [Guerini et al., 2014; Guerini et al., 2009; Torres, Westover, Gibbons, Johnson, & Ward, 2012], which is consistent with our findings at HLA‐C. Furthermore, two additional class I loci, HLA‐B and HLA‐E, were also differentially expressed in mothers of children classified as ASD and Non‐TD compared to TD children in this study, providing further evidence of an MHC I response in pregnancies of non‐typical neurodevelopment. These findings are consistent with the known maternal effect of viral infection during pregnancy of increasing susceptibility to ASD and schizophrenia during pregnancy [al‐Haddad et al., 2019; Zerbo et al., 2013]. They also suggest how common genetic polymorphisms at these HLA loci may interact with common environmental factors such as viral infection during pregnancy to impact diverse neurodevelopmental outcomes.

Furthermore, gene pathway analysis of differentially expressed genes common between ASD and Non‐TD revealed enrichment for the interferon‐gamma mediated signaling pathway, which has been previously found to be elevated in mothers of children with ASD and other neurodevelopmental disorders [Goines et al., 2011; Krakowiak et al., 2017]. In one such study, elevated interferon‐gamma levels in maternal midgestation peripheral blood was associated with a 50% increased risk of offspring ASD risk [Goines et al., 2011]. A second enriched pathway included CaMK family members which play an important role in neuronal connectivity and synaptic plasticity [Bourgeron, 2015; Stamou, Streifel, Goines, & Lein, 2013; Zimmerman, Pessah, & Lein, 2008] as well as immune response and inflammation [Racioppi & Means, 2008]. Prior ASD studies have implicated the CaMK pathway in dendritic growth and local connectivity alterations related to gene–environment interactions in ASD [Bourgeron, 2015; Stamou et al., 2013; Zimmerman et al., 2008].

Although genome‐wide significance of individual differentially expressed genes was not observed between samples from mothers whose children developed ASD compared to TD after adjusting for multiple comparisons, seven nominally‐significant genes were also on the SFARI list of strong ASD candidate genes. TRIO, Trio Rho guanine nucleotide exchange factor, promotes exchange of GDP for GTP and provides necessary support for cell migration and cell growth related to Alzheimer disease and other neurological conditions [De Rubeis et al., 2014; Katrancha et al., 2017]. GRIA1, encoding a receptor for glutamate, the predominant excitatory neurotransmitter in brain, is activated by normal neurophysiologic processes [De Rubeis et al., 2014; Geisheker et al., 2017]. SMARCC2 encodes a chromatin remodeling protein with helicase and ATPase activities that has been implicated in altering chromatin structure in ASD [Devlin et al., 2012]. CNTN4 functions in neuronal network formation and plasticity, and is associated with nervous system development at the transcriptome level [Fernandez et al., 2004; Yoshihara et al., 1995]. Mutations in FOXP1, a developmental transcription factor, are observed in rare cases of intellectual disability with ASD [Ferland, Cherry, Preware, Morrisey, & Walsh, 2003; Teramitsu, Kudo, London, Geschwind, & White, 2004]. Dysregulated expression of these ASD risk genes in maternal blood could reflect an underlying shared genetic risk for ASD in these high‐risk families.

Methylation and Methyl‐binding Functions in a Gene Module Oppositely Associated with Folic Acid and Risk

The complementary co‐expression network analysis further revealed a module of 224 co‐expressed genes in maternal blood showing an association with folic acid and 5‐MeTHF levels in the opposite direction from ASD risk that could not be explained by cell type differences. Interestingly, these ASD and nutrient associated genes were functionally enriched for DNA methylation binding and methylation‐dependent chromatin silencing, consistent with prior DNA methylation changes observed in ASD [Coulson et al., 2018; Mordaunt et al., 2020; Vogel Ciernia et al., 2019; Zhu et al., 2019] as well as ASD‐like syndromes associated with methyl binding proteins [Cukier et al., 2012; Cukier et al., 2010]. MBD3L, which has methyl‐binding function, is predicted to assist with demethylation reactions and functions as a transcriptional repressor [Fouse, Nagarajan, & Costello, 2010; Mungall, 2002; Zhou et al., 2019]. Folic acid, the synthetic form of folate that contributes the substrate for one‐carbon metabolism, and 5‐MeTHF, one of the active biological forms of folate that plays a critical role in one‐carbon metabolism, have also been shown to be inversely associated with developmental delay.

One‐carbon metabolites associated with changes in gene expression in this study have also been associated with the prevention of neurodevelopmental conditions [Afman & Müller, 2006; Mordaunt, Kieffer, et al., 2019; Waterland & Jirtle, 2004; Zhu et al., 2019]. The co‐regulated block of betaine and DMG co‐expression modules contained genes enriched in the adaptive immune system and chromatin modification functions, as well as Rett syndrome, a known syndromic form of ASD [Amir et al., 1999; Craig, 2004; Ducker & Rabinowitz, 2017; Paparo et al., 2014]. Choline is metabolized to betaine, which converts homocysteine to form methionine, generating DMG in the one‐carbon pathway [Ueland, Holm, & Hustad, 2005; Zeisel & Blusztajn, 1994]. A previous study of maternal peripheral blood collected at term showed that changes in betaine and DMG were in the opposite direction from choline when compared with nonpregnant women [Friesen, Novak, Hasman, & Innis, 2007]. EVL (Enah/Vasp‐like) is involved in actin cytoskeleton remodeling and is crucial for central nervous system processes and immune system functions [Gardiner et al., 2013; Krause, Dent, Bear, Loureiro, & Gertler, 2003; Tsunoda et al., 2015]. One study also showed EVL as a differentially expressed gene in schizophrenia in peripheral blood [Gardiner et al., 2013].

Previous studies in ASD have been focused on post mortem brain tissue [Ginsberg, Rubin, Falcone, Ting, & Natowicz, 2012; Voineagu et al., 2011], as a tissue relevant to the condition, but collected after diagnoses were made, raising concerns about reverse causation in determining etiologically‐relevant expression changes. Few studies have focused on prospective transcriptomic profiles collected prior to the presentation of the condition [Glatt et al., 2012; Mordaunt, Park, et al., 2019], and none have examined maternal gene expression profiles during pregnancies at high risk for developing ASD. In addition, few studies have integrated maternal transcriptome and one‐carbon metabolite data within biospecimens. Furthermore, most studies of ASD expression biomarkers have not considered the roles of nutritional factors during pregnancy that could be relevant to fetal development. Although maternal one‐carbon metabolites were not significantly associated with child outcome in this study, they were associated with transcriptomic differences that were relevant to child neurodevelopmental outcomes. This could be explained by differential maternal responses to similar concentrations of maternal metabolites, or if there were ongoing transcriptomic changes in response to nutrient metabolite concentrations, but that there were timing‐specific metabolite concentrations (e.g., earlier in gestation, as indicated in some studies [Schmidt et al., 2011; Schmidt et al., 2019; Schmidt et al., 2012; Suren et al., 2013]) that are relevant to neurodevelopmental outcomes and we did not capture all potentially critical time points in the metabolite measures.

A limitation of using maternal peripheral blood to examine expression is that it contains multiple cell types, and proportions can differ across samples. However, estimated cell type composition of maternal blood was not significantly associated with the child's clinical outcomes, the “greenyellow” module, betaine, or DMG, which suggests that our main findings were not driven by differential cell type proportions. After correcting for multiple comparisons, this study did not identify any individual differentially maternally expressed genes specifically associated with ASD, although 6 transcripts from 4 genes reached genome‐wide significance with diagnosis of either ASD or Non‐TD. Another potential limitation is the relatively low number of samples with certain metabolite measurements available. Future studies on maternal nutritional factors and child outcomes would be beneficial in confirming our data. Furthermore, lack of genome‐wide evidence of individual differentially expressed genes specific to a pairwise comparison of ASD versus TD is likely due to the relatively small sample size that is inherent to a prospective ASD study, that is underpowered to detect small differences in transcript levels. However, this does not eliminate the importance of identifying and understanding the biologically significant gene set enrichments and co‐expression network modules using differential gene expression and WGCNA analysis. Additionally, other factors, including genetics, epigenetics, gestational age, and other environmental factors can influence the transcriptome and ASD risk. Approaches incorporating those factors will be important in future studies.

In summary, genome‐wide gene expression analysis of maternal peripheral blood samples revealed transcriptome changes associated with maternal one‐carbon metabolites and child neurodevelopmental outcomes implicating maternal immune, apoptotic, and epigenetic mechanisms in the development of ASD in offspring. In addition, folic acid and 5‐MeTHF were associated with expression of genes involved in methylated‐CpG binding in an opposite direction to that of ASD, consistent with prior evidence of protection. Finally, maternal betaine and DMG levels clustered with co‐expressed genes related to immune, chromatin modification, and development functions. These results therefore provide important biological insights into maternal gene pathways associated with adverse neurodevelopment in the child, as well as the suggested protective association with one carbon metabolites in the complex etiology of ASD.

Conflict of Interest

The authors declare that there is no conflict of interest.

Declarations

Ethics Approval and Consent to Participate

The UC Davis Institutional Review Board and the State of California Committee for the Protection of Human Subjects approved this study and the MARBLES Study protocols. All data and specimens were collected after parent given written informed consent form.

Availability of data and material

Data are shared in the Gene Expression Omnibus (GEO) accession number (GSE148450) based on participates consent. Code and scripts for this study are available on GitHub (https://github.com/Yihui‐Zhu/AutismMaternalBloodExpression). Other related data and information are included in supplementary tables.

Figure S1. Density plot of metabolites and transcriptome on gestational age of maternal blood draw.

Figure S2. Surrogate variable analysis in MARBLES subjects.

Figure S3. Differential expressed transcripts on using known confounders and SVA adjustment.

Figure S4. Gene ontology and pathway analysis using directed acyclic graph (DAG) on the 218 transcripts common to ASD versus TD and Non‐TD versus TD differentially expressed gene lists.

Figure S5. Probe intensity density plot based all 300 maternal blood Affymetrix data.

Figure S6. Co‐expression network modules with diagnosis, demographic factors, and maternal blood nutrient concentrations. The values in the cells represent Pearson r (adjusted P‐value). P‐value were adjusted for all comparisons.

Figure S7. Co‐expression network modules with diagnosis, demographic factors, and maternal blood nutrient concentrations. The values in the cells represent Pearson r (P‐value). P‐value shown here were unadjusted P‐value without adjustment.

Figure S8. Unsupervised hierarchical clustering adjacency heatmap correlation and P‐value.

Figure S9. Gene ontology and pathway analysis for the block of eight weighted gene co‐expression modules associated with betaine and DMG.

Figure S10. Association network on hub genes from the block of eight weighted gene co‐expression modules associated with betaine and DMG.

Figure S11. Heatmap of correlation between module eigengenes and cell type proportions with FDR adjusted P‐value.

Figure S12. Heatmap of correlation between sample demographic factors and nutrients and cell type proportions with FDR adjusted P‐value.

Click here for additional data file.

Table S1. Sample variables and RNA quality in MARBLES subjects.

Table S2. Demographic characteristics of mother participants and their children in MARBLES, stratified by child diagnosis outcomes.

Table S3. Descriptive statistics of maternal peripheral blood nutrients level in MARBLES, stratified by children diagnosis.

Table S4. Differential expression analysis on maternal gene expression and diagnosis on adjustment using either SVA or known confounders.

Table S5. ASD related significant genes from differential expression analysis.

Table S6. Non‐TD related significant genes from differential expression analysis.

Table S7. Gene Ontology terms on all transcripts and filtered transcripts as background.

Table S8. Weighted gene co‐expression network module memberships for each transcript on the array (MM: Pearson correlation coefficient for module membership; p.MM: P value for the preceding relationship).

Table S9. Weighted gene co‐expression network module features, including number of transcripts and hub genes characters.

Table S10. “Greenyellow” gene co‐expression network module memberships on 224 transcripts (MM: Pearson correlation coefficient for module membership; p.MM: P value for the preceding relationship).

Table S11. “Greenyellow” gene co‐expression network module 224 transcripts gene ontology terms and gene lists.

Table S12. Modules significant associated with trimester during pregnancy and their enriched gene ontology terms.

Table S13. Eight weighted gene co‐expression modules block memberships on 2,582 transcripts.

Table S14. Cell type proportions in all 300 maternal peripheral blood samples estimated with CIBERSORT.

Click here for additional data file.