A meta-analysis of two high-risk prospective cohort studies reveals autism-specific transcriptional changes to chromatin, autoimmune, and environmental response genes in umbilical cord blood.

Background: Autism spectrum disorder (ASD) is a neurodevelopmental disorder that affects more than 1% of children in the USA. ASD risk is thought to arise from both genetic and environmental factors, with the perinatal period as a critical window. Understanding early transcriptional changes in ASD would assist in clarifying disease pathogenesis and identifying biomarkers. However, little is known about umbilical cord blood gene expression profiles in babies later diagnosed with ASD compared to non-typically developing and non-ASD (Non-TD) or typically developing (TD) children.
Methods: Genome-wide transcript levels were measured by Affymetrix Human Gene 2.0 array in RNA from cord blood samples from both the Markers of Autism Risk in Babies-Learning Early Signs (MARBLES) and the Early Autism Risk Longitudinal Investigation (EARLI) high-risk pregnancy cohorts that enroll younger siblings of a child previously diagnosed with ASD. Younger siblings were diagnosed based on assessments at 36 months, and 59 ASD, 92 Non-TD, and 120 TD subjects were included. Using both differential expression analysis and weighted gene correlation network analysis, gene expression between ASD and TD, and between Non-TD and TD, was compared within each study and via meta-analysis.
Results: While cord blood gene expression differences comparing either ASD or Non-TD to TD did not reach genome-wide significance, 172 genes were nominally differentially expressed between ASD and TD cord blood (log2(fold change) > 0.1, p < 0.01). These genes were significantly enriched for functions in xenobiotic metabolism, chromatin regulation, and systemic lupus erythematosus (FDR q < 0.05). In contrast, 66 genes were nominally differentially expressed between Non-TD and TD, including 8 genes that were also differentially expressed in ASD. Gene coexpression modules were significantly correlated with demographic factors and cell type proportions. Limitations: ASD-associated gene expression differences identified in this study are subtle, as cord blood is not the main affected tissue, it is composed of many cell types, and ASD is a heterogeneous disorder. Conclusions: This is the first study to identify gene expression differences in cord blood specific to ASD through a meta-analysis across two prospective pregnancy cohorts. The enriched gene pathways support involvement of environmental, immune, and epigenetic mechanisms in ASD etiology.

Background

Autism spectrum disorder (ASD) is a neurodevelopmental disorder characterized by impaired social interaction and restricted and repetitive behaviors. Heritability of ASD risk has been well established with twin and family studies and is estimated at 52% [1-3]. While rare variants with large effects explain a relatively small proportion of all ASD cases, heritable common variants with individually minor effects contribute substantially to ASD risk [4]. Accumulating lines of evidence suggest that ASD arises from complex interactions between heterogeneous genetic and environmental risk factors. Gene expression levels are influenced by both genetic and environmental factors and determine the functional responses of cells and tissues. Postmortem brain gene expression studies have guided understanding of ASD pathophysiology and show evidence of changes in gene coexpression and enrichment in immune response and neuronal activity functions [5, 6]. Peripheral blood gene expression studies in children and adults using whole blood and in specific cell types (natural killer (NK) cell and lymphocytes) observed enrichment of immune and inflammatory processes in differential gene expression associated with ASD [7, 8]. Recent efforts have been focused on identifying how genetic risk factors converge into one or more unifying pathways and pathophysiological mechanisms [9, 10]. Yet, the majority of this work to date relies on postmortem or postsymptom timing of sample collection, rather than prospective assessment of gene expression.

Converging evidence suggests that most of the changes in the brain associated with ASD are initiated during prenatal brain development [11, 12], but the complete nature of these changes remains unknown. Umbilical cord blood captures fetal blood as well as the exchanges across the fetoplacental unit and provides a distinct insight into prenatal development. A unique cell mixture is represented in umbilical cord blood, including hematopoietic stem cells, B cells, NK cells, T cells, monocytes, granulocytes, and nucleated red blood cells [13]. Cord blood gene expression would reflect the immune response as well as endocrine and cellular communication essential for fetal development near the time of birth.

While several studies have previously examined child blood gene expression differences in ASD [8, 14–20], this is the first study to take advantage of cord blood samples collected from two prospective studies (Markers of Autism Risk in Babies-Learning Early Signs (MARBLES) and the Early Autism Risk Longitudinal Investigation (EARLI)) in order to assess the perinatal transcriptional changes that precede ASD diagnosis in high-risk children [21, 22]. The subjects in this study are all siblings of children with ASD and thus have a 13-fold increased risk for ASD compared to the general population [23]. They are also at a higher risk for non-typical neurodevelopment, including deficits in attention and behavior. We measured cord blood gene expression levels using the Affymetrix Human Gene 2.0 array and compared gene-level differential expression, gene set enrichment, and gene coexpression networks across ASD, non-typically developing (Non-TD) and neurotypical children (Additional file 1: Figure S1). Study-level results were then combined in a meta-analysis to investigate cord blood transcriptional dysregulation in ASD.

Methods

Sample population and biosample collection

MARBLES

The MARBLES study recruits Northern California mothers from lists of children receiving services through the California Department of Developmental Services who have a child with confirmed ASD and are planning a pregnancy or are pregnant with another child. Inclusion criteria for the study were as follows: (1) mother or father has one or more biological child(ren) with ASD; (2) mother is 18 years or older; (3) mother is pregnant; (4) mother speaks, reads, and understands English sufficiently to complete the protocol and the younger sibling will be taught to speak English; and (5) mother lives within 2.5 h of the Davis/Sacramento region at time of enrollment. As described in more detail elsewhere [21], demographic, diet, lifestyle, environmental, and medical information were prospectively collected through telephone-assisted interviews and mailed questionnaires throughout pregnancy and the postnatal period. Mothers were provided with sampling kits for cord blood collection prior to delivery. MARBLES research staff made arrangements with obstetricians/midwives and birth hospital labor and delivery staff to assure proper sample collection and temporary storage. Infants received standardized neurodevelopmental assessments beginning at 6 months, as described below, and concluding at 3 years of age. For this study, all children actively enrolled by March 1, 2017 (n = 347), with umbilical cord blood collected in a PAXgene Blood RNA tube (n = 262, 76%) were included.

EARLI

The EARLI study is a high-risk pregnancy cohort that recruited and followed pregnant mothers who had an older child diagnosed with ASD through pregnancy, birth, and the first 3 years of life. EARLI families were recruited at four EARLI Network sites (Drexel/Children’s Hospital of Philadelphia, Johns Hopkins/Kennedy Krieger Institute, University of California (UC) Davis, and Kaiser Permanente Northern California) in three distinct US regions (Southeast Pennsylvania, Northeast Maryland, and Northern California). In addition to having a biological child with ASD confirmed by EARLI study clinicians, to be eligible, mothers also had to communicate in English or Spanish and, at recruitment, meet the following criteria: be 18 years or older, live within 2 h of a study site, and be < 29 weeks pregnant. The design of the EARLI study is described in more detail in Newschaffer et al. [22]. EARLI research staff made arrangements with obstetricians/midwives and birth hospital labor and delivery staff to ensure proper cord blood sample collection and temporary storage. The development of children born into the cohort was closely followed through age 3 years. For this study, 212 infants born into EARLI as a singleton birth and followed to 1 year of age were considered for inclusion. Of the 212 infants, 97 were excluded because they were either missing umbilical cord blood samples or outcome measures at 36 months, leaving a final sample of 115.

Diagnostic outcomes

In both studies, development was assessed by trained, reliable examiners. Diagnostic assessments at 3 years included the gold standard Autism Diagnostic Observation Schedule (ADOS) [24, 25], the Autism Diagnostic Interview-Revised (ADI-R) [26] conducted with parents, and the Mullen Scales of Early Learning (MSEL) [27], a test of cognitive, language, and motor development. Participants were classified into one of three outcome groups, ASD, typically developing (TD), and Non-TD, based on a previously published algorithm that uses ADOS and MSEL scores [28, 29]. Children with ASD outcomes had scores over the ADOS cutoff and met DSM-5 criteria for ASD. The Non-TD group was defined as children with low MSEL scores (i.e., two or more MSEL subscales that are more than 1.5 standard deviations (SD) below average or at least one MSEL subscale that was more than 2 SD below average), elevated ADOS scores (i.e., within 3 points of the ASD cutoff), or both. Children with TD outcomes had all MSEL scores within 2 SD and no more than one MSEL subscale 1.5 SD below the normative mean and scores on the ADOS at least 3 or more points below the ASD cutoff.

Demographic characteristics

In both the MARBLES and EARLI studies, demographic information was prospectively collected through in-person and telephone-assisted interviews and mailed questionnaires throughout pregnancy and the postnatal period. Cotinine was measured in maternal urine during pregnancy, and maternal smoking was identified if the concentration of cotinine was > 50 ng/mL [30]. Within each study, demographic characteristics were stratified by diagnostic outcome and compared using Fisher’s exact test for categorical variables and one-way ANOVA for continuous variables.

RNA isolation and expression assessment

In both MARBLES and EARLI, umbilical cord blood was collected at the time of birth in PAXgene Blood RNA Tubes with the RNA stabilization reagent (BD Biosciences) and stored at − 80 °C. RNA isolation was performed with the PAXgene Blood RNA Kit (Qiagen) following the manufacturer’s protocol. RNA from 236 (90%) of the 262 MARBLES PAXgene blood samples and all of the EARLI PAXgene blood samples met quality control standards (RIN ≥ 7.0 and concentration ≥ 35 ng/uL) and volume requirements. Total RNA was converted to cDNA and in vitro transcribed to biotin-labeled cRNA, which was hybridized to Human Gene 2.0 Affymetrix microarray chips by the Johns Hopkins Sequencing and Microarray core. EARLI and MARBLES samples were measured separately and in multiple batches within each study. The manufacturer’s protocol was followed for all washing, staining, and scanning procedures. The raw fluorescence data (in Affymetrix CEL file format) with one perfect match and one mismatched probe in each set were analyzed using the oligo package in R.

Data preprocessing

Within each study, signal distribution was first assessed in perfect-match probe intensity and robust multi-chip average (RMA) normalized data [31]. During the quality control step, we identified outliers using the arrayQualityMetrics and oligo R packages [32, 33]. Outliers were excluded based on loading in principal component 1, the Kolmogorov-Smirnov test, median normalized unscaled standard error, and the sum of the distances to all other arrays. For the MARBLES study, 3 outlier samples were identified and excluded, and another 71 children had not yet received a diagnosis by April 12, 2018, so were excluded; 162 samples were normalized using RMA. For the EARLI study, 6 outliers were identified and excluded, then 109 samples were normalized using RMA. Probes were annotated at the transcript level using the pd.hugene.2.0.st R package [34], and those assigned to a gene (36,459 probes) were used in subsequent analyses.

Surrogate variable analysis

Surrogate variable analysis (SVA) was used to estimate and adjust for unmeasured environmental, demographic, cell type proportion, and technical factors that may have substantial effects on gene expression using the SVA R package [35, 36]. Twenty-one surrogate variables were detected in normalized expression data from MARBLES subjects for both the ASD versus TD and Non-TD versus TD comparisons. Specific factors associated with surrogate variables in MARBLES using linear regression included array batch, sex, maternal BMI, gestational age, delivery method, child ethnicity, and maternal education (false discovery rate (FDR) q < 0.1, Additional file 1: Figure S2a). In normalized expression data from EARLI subjects, 11 surrogate variables were detected for the ASD versus TD comparison, which were associated with sex, birth weight, gestational age, and paternal age (FDR q < 0.1, Additional file 1: Figure S3a). Twelve surrogate variables were detected for the Non-TD versus TD comparison, which were associated with sex and gestational age (FDR q < 0.1, Additional file 1: Figure S4a). Proportion of variance in expression of each gene explained by each surrogate variable was determined using the variancePartition R package [37]. Median variance explained by each surrogate variable ranged from 0.3 to 5.6% in MARBLES (Additional file 1: Figure S2b). In EARLI, median variance explained by each surrogate variable ranged from 0.8 to 7.1% for the ASD versus TD comparison and 0.5 to 7.2% for the Non-TD versus TD comparison (Additional file 1: Figures S3b and S4b).

Differential gene expression

Differential expression was determined using the limma package in R with diagnosis and all surrogate variables included in the linear model [38] (Additional file 1: Figures S5 and S6). ASD versus TD and Non-TD versus TD differential expression results were extracted from a single model with three levels for diagnosis for MARBLES, while two pairwise models were used for EARLI, although this did not affect the results (1 vs 2 model meta-analysis fold change ASD vs TD Pearson’s r = 0.97, Non-TD vs TD Pearson’s r = 0.99). Fold change and standard error from each study were input into the METAL command-line tool for meta-analysis using the standard error analysis scheme with genomic control correction [39]. In this approach, the fold changes from each study are weighted using the inverse of the standard error. Using the meta-analyzed data, differential probes were then identified as those with a nominal p value < 0.01 and an average absolute log2(fold change) > 0.1.

Gene overlap analysis

Gene overlap analysis by Fisher’s exact test was performed using the GeneOverlap R package [40]. Gene symbols annotated to differentially expressed probes were compared to autism-related or blood cell type-associated gene lists [41] for overlap relative to all genes annotated to probes on the array. Genes with variation previously associated with autism were obtained from the Simons Foundation Autism Research Initiative (SFARI) Gene database and a recent genome-wide association study meta-analysis [42, 43], while genes with expression previously associated with autism were obtained from multiple previous reports [6, 8, 44, 45]. Significant overlaps were those with an FDR q value < 0.05.

Overrepresentation enrichment analysis

Differential probes identified during meta-analysis were converted to Entrez gene IDs using the biomaRt R package [46]. Functional enrichment of only differential probes by hypergeometric test was relative to all probes on the array and was performed using the WebGestalt online tool with default parameters for the overrepresentation enrichment analysis method [47]. Enrichment databases included WebGestalt defaults and also a custom database of recently evolved genes obtained from [48]. WebGestalt default databases queried included Gene Ontology, KEGG, WikiPathways, Reactome, PANTHER, MSigDB, Human Phenotype Ontology, DisGeNET, OMIM, PharmGKB, and DrugBank. Significant enrichments were those with an FDR q value < 0.05.

Gene set enrichment analysis

All probes included in the analysis were ranked using meta-analysis log2(fold change) and input into the WebGestalt online tool using default parameters for the GSEA method [47]. GSEA assesses whether genes in biologically predefined sets occur toward the top or bottom of a ranked list of all examined genes more than expected by chance [49]. GSEA calculates an enrichment score normalized to the set size to estimate the extent of non-random distribution of the predefined gene set, and it then tests the significance of the enrichment with a permutation test. Enrichment databases included WebGestalt defaults (see above). Significant gene sets were called as those with an FDR q value < 0.05.

Weighted gene correlation network analysis

Weighted gene correlation network analysis (WGCNA) was performed using the WGCNA R package [50]. RMA-normalized expression data were adjusted for batch using ComBat due to large batch effects in MARBLES [36]. Samples were clustered with hierarchical clustering using the average method and excluded using static tree cutting with cut height set to 100, resulting in one outlier being removed for each study [50]. Expression data were renormalized with RMA and adjusted for batch after removing outliers. Final samples for WGCNA included 59 ASD (41 MARBLES/19 EARLI), 91 Non-TD (44 MARBLES/47 EARLI), and 119 TD (76 MARBLES/43 EARLI). Signed topological overlap matrices (TOMs) were obtained separately for each study in a single block using the biweight midcorrelation with the maximum percentile for outliers set to 0.1 and the soft-thresholding power set to 10. Study-specific TOMs were calibrated using full quantile normalization, and a consensus TOM was calculated as the parallel minimum of the study-specific TOMs. Modules were identified using dynamic hybrid tree cutting with deepSplit set to 4, and modules with a dissimilarity < 0.1 were merged. Module hub probes were determined as the probe in each module with the highest module membership. Study-specific module eigengenes were correlated with demographic factors or estimated cell type proportions using the biweight midcorrelation with the maximum percentile for outliers set to 0.1 and including only pairwise-complete observations. Study-specific correlation Z-scores were combined in a meta-analysis using Stouffer’s method with weights given by the square root of the sample n [51]. p values were adjusted for all comparisons using the FDR method. Significant correlations were called as those with an FDR q value < 0.05.

Cell type proportion deconvolution

Estimation of cell type proportions was performed using CIBERSORT [41]. Final samples for cell type deconvolution were the same as for WGCNA: 59 ASD (41 MARBLES/19 EARLI), 91 Non-TD (44 MARBLES/47 EARLI), and 119 TD (76 MARBLES/43 EARLI). RMA-normalized expression data was adjusted for batch due to large batch effects in MARBLES. To correspond with identifiers used by CIBERSORT, array probes were matched to HUGO Gene Nomenclature Committee (HGNC) gene symbols using the biomaRt R package [46]. RMA-normalized expression data and the default LM22 adult blood signature gene files were input into the CIBERSORT web tool [41]. A similar signature gene file was not available for umbilical cord blood and so an adult reference panel was used. Relative and absolute modes were run together, with 100 permutations and without quantile normalization. Deconvolution goodness of fit p value was < 0.05 for all subjects. Estimated cell type proportions were correlated with demographic factors within each study using the biweight midcorrelation with the maximum percentile for outliers set to 0.1 and including only pairwise-complete observations. Study-specific correlation Z-scores were combined in a meta-analysis using Stouffer’s method with weights given by the square root of the sample n [51]. p values were adjusted for all comparisons using the FDR method. Significant correlations were called as those with an FDR q value < 0.05.

Results

Study sample characteristics

MARBLES subjects in the final analysis included 41 ASD (30 male, 11 female), 44 Non-TD (27 male, 17 female), and 77 TD subjects (40 male, 37 female). Paternal age and gestational age were nominally associated with diagnostic group in MARBLES, with slightly increased paternal age and gestational age for the ASD subjects (paternal age p = 0.02, gestational age p = 0.04, Table 1). Other demographic characteristics were not associated with diagnostic group among MARBLES subjects. EARLI subjects in the final analysis included 18 ASD (13 male, 5 female), 48 Non-TD (23 male, 25 female), and 43 TD subjects (19 male, 24 female). Child race and ethnicity and home ownership were nominally associated with diagnostic group in EARLI (race and ethnicity p = 0.02, home ownership p = 0.01, Table 2). Specifically, the ASD group included a lower proportion of white subjects and a lower rate of home ownership. Other demographic characteristics were not associated with diagnostic group among EARLI subjects. In the meta-analysis, which combined both the MARBLES and EARLI studies, gene expression was analyzed in 271 subjects, including 120 TD, 59 ASD, and 92 Non-TD subjects.

Table 1

Demographic characteristics of children and their parents in the MARBLES study, stratified by child diagnosis

	Child 36 month outcome
	ASD (n = 41)	Non-TD (n = 44)	TD (n = 77)	p valuea
Child male gender, n (%)	30 (73.2)	27 (61.4)	40 (51.9)	0.08
Child race and ethnicity, n (%)				0.90
White (European or Middle Eastern)	19 (46.3)	18 (40.9)	38 (49.4)
Black/African-American	2 (4.9)	2 (4.5)	1 (1.3)
Asian or Pacific Islander	4 (9.8)	3 (6.8)	7 (9.1)
Hispanic White	8 (19.5)	11 (25)	19 (24.7)
Hispanic Non-White	4 (9.8)	3 (6.8)	6 (7.8)
Multiracial	4 (9.8)	7 (15.9)	6 (7.8)
Mother age (years), mean (SD)b	35.4 (5.1)	34.1 (4.3)	33.3 (4.9)	0.09
Father age (years), mean (SD)c	38.7 (5.3)	36.1 (4.8)	35.7 (6.2)	0.02 k
Mother prepregnancy BMI, mean (SD)d	27.9 (7.4)	29.3 (8.3)	26.1 (6.2)	0.05
Mother type 1 diabetes, n (%)e	0 (0)	0 (0)	1 (1.3)	1.00
Mother type 2 diabetes, n (%)e	0 (0)	0 (0)	1 (1.3)	1.00
Mother gestational diabetes, n (%)f	12 (30)	6 (14.3)	11 (14.7)	0.12
Mother hypertension, n (%)e	1 (2.6)	1 (2.3)	1 (1.3)	1.00
Mother preeclampsia, n (%)g	2 (5.1)	4 (10)	2 (2.7)	0.29
Mother smoking by urine cotinine, n (%)h	4 (17.4)	1 (4.3)	1 (2.9)	0.12
Mother bachelor’s degree +, n (%)	18 (43.9)	19 (43.2)	40 (51.9)	0.57
Father bachelor’s degree +, n (%)	21 (51.2)	19 (43.2)	39 (50.6)	0.69
Own home, n (%)i	21 (55.3)	25 (58.1)	48 (63.2)	0.69
Cesarean delivery, n (%)j	11 (26.8)	17 (38.6)	31 (41.3)	0.28
Gestational age (weeks), mean (SD)	39.4 (1.1)	39 (1.3)	38.7 (1.4)	0.04 k
Birth weight (kg), mean (SD)	3.5 (0.4)	3.5 (0.4)	3.4 (0.5)	0.52

ASD autism spectrum disorder, Non-TD non-typically developing, TD typically developing

ap values from Fisher’s exact test for categorical variables and one-way ANOVA for continuous variables

bFrequency missing = 1 ASD

cFrequency missing = 2 ASD, 1 TD

dFrequency missing = 1 ASD, 1 TD

eFrequency missing = 2 ASD, 1 Non-TD, 1 TD

fFrequency missing = 1 ASD, 2 Non-TD, 2 TD

gFrequency missing = 2 ASD, 4 Non-TD, 2 TD

hFrequency missing = 18 ASD, 21 Non-TD, 42 TD

iFrequency missing = 3 ASD, 1 Non-TD, 1 TD

jFrequency missing = 2 TD

kNominally-significant association with child diagnosis

Table 2

Demographic characteristics of children and their parents in the EARLI study, stratified by child diagnosis

	Child 36 month outcome
	ASD (n = 18)	Non-TD (n = 48)	TD (n = 43)	p valuea
Child male gender, n (%)	13 (72.2)	23 (47.9)	19 (44.2)	0.13
Child race and ethnicity, n (%)b				0.02 h
White (European or Middle Eastern)	6 (35.3)	23 (48.9)	24 (55.8)
Black/African-American	3 (17.6)	8 (17)	1 (2.3)
Asian or Pacific Islander	2 (11.8)	5 (10.6)	0 (0)
Hispanic White	2 (11.8)	3 (6.4)	3 (7)
Hispanic Non-White	2 (11.8)	7 (14.9)	8 (18.6)
Multiracial	2 (11.8)	1 (2.1)	7 (16.3)
Mother age (years), mean (SD)	34.4 (4)	33.8 (4.3)	34.7 (5)	0.60
Father age (years), mean (SD)	35.4 (6.8)	35.7 (5)	36.2 (6.4)	0.88
Mother prepregnancy BMI, mean (SD)c	29.4 (7.1)	27.4 (7)	27.9 (7.2)	0.60
Mother smoking by urine cotinine, n (%)d	1 (6.3)	1 (2.3)	0 (0)	0.30
Mother bachelor’s degree +, n (%)	8 (44.4)	28 (58.3)	29 (67.4)	0.26
Father bachelor’s degree +, n (%)e	7 (38.9)	28 (62.2)	25 (58.1)	0.24
Own home, n (%)	8 (44.4)	31 (64.6)	36 (83.7)	0.01 h
Cesarean delivery, n (%)f	9 (56.3)	17 (39.5)	12 (35.3)	0.38
Gestational age (weeks), mean (SD)	39.2 (1.4)	39.4 (1.5)	39.1 (1.5)	0.70
Birth weight (kg), mean (SD)g	3.5 (0.6)	3.5 (0.6)	3.5 (0.6)	1.00

ASD autism spectrum disorder, Non-TD non-typically developing, TD typically developing

ap values from Fisher’s exact test for categorical variables and one-way ANOVA for continuous variables

bFrequency missing = 1 ASD, 1 Non-TD

cFrequency missing = 1 Non-TD

dFrequency missing = 2 ASD, 4 Non-TD, 3 TD

eFrequency missing = 3 Non-TD

fFrequency missing = 2 ASD, 5 Non-TD, 9 TD

gFrequency missing = 1 TD

hNominally-significant association with child diagnosis

ASD-associated differential gene expression in cord blood

We examined differential expression of single genes in cord blood samples in association with ASD diagnosis status at 36 months. In meta-analysis, no transcripts were differentially expressed at a conservative FDR q value < 0.05. Under the thresholds of log2(fold change) > 0.1 and nominal p value < 0.01, 172 transcripts were differentially expressed between ASD and TD cord blood (ASD n = 59, TD n = 120, Fig. 1a, Additional file 2: Table S1). Among these differential transcripts, 87 were upregulated and 85 were downregulated, and the median absolute log2(fold change) was 0.12. The differential transcript with the greatest absolute fold change was TUBB2A (log2(fold change) = 0.35, standard error = 0.12, p = 4.8E−3, Fig. 1b, Table 3). Additionally, the estimated fold changes for differentially expressed genes were strongly correlated between the two studies (Pearson’s r = 0.80, p < 2.2E−16), although the fold changes of all transcripts were weakly correlated (Pearson’s r = 0.02, p = 4.6E−4, Additional file 1: Figure S7a). Many of the differentially expressed genes were non-coding or uncharacterized transcripts; however, the median expression of differentially expressed genes was not lower than non-differentially expressed genes on the array (MARBLES: differential = 4.70, non-differential = 4.64, p = 0.74; EARLI: differential = 4.34, non-differential = 4.19, p = 0.52; Additional file 1: Figure S8).

Fig. 1

Identification and function of ASD-associated differentially expressed genes in cord blood from two high-risk prospective studies. Gene expression in umbilical cord blood samples from subjects with typical development (n = 120, 59 male/61 female) or those diagnosed with ASD at age 36 months (n = 59, 43 male/16 female) was assessed by expression microarray. SVA was performed to control for technical and biological variables including sex and array batch. Differential expression analysis was carried out separately by study and combined in a meta-analysis. a Identification of 172 differentially expressed genes in meta-analysis (178 probes, log2(fold change) > 0.1, p < 0.01). b Fold change in gene expression for top 20 differentially expressed genes sorted by meta-analysis log2(fold change) and plotted for individual studies and meta-analysis. c ASD-associated differentially expressed genes were analyzed for functional enrichment with WebGestalt using the hypergeometric test and compared to all genes annotated to the array. Significantly enriched ontology terms are shown (FDR q < 0.05). d ASD-associated differentially expressed genes belonging to significantly enriched ontology terms are listed

Table 3

Top 20 ASD-associated differentially expressed genes by log2(fold change) in meta-analysis

Probe ID	Gene symbol	MARBLES		EARLI		Meta-analysis
		log2(FC)	p value	log2(FC)	p value	log2(FC)	SE	p value	q value
17004333	TUBB2A	0.36	9.6E−03	0.33	2.1E−01	0.35	0.12	4.8E−03	0.75
17094496	MIR1299	0.28	6.1E−03	0.15	2.3E−01	0.23	0.08	4.3E−03	0.75
17065047	LOC389831	0.26	6.1E−03	0.16	2.4E−01	0.23	0.08	4.2E−03	0.74
17063807	TRBV11-2	0.22	1.7E−02	0.19	1.7E−01	0.21	0.08	7.0E−03	0.79
16812322	RNU6-667P	0.13	1.4E−01	0.34	4.1E−03	0.20	0.07	3.7E−03	0.73
16662296	RNA5SP42	0.22	3.1E−03	0.12	4.4E−01	0.20	0.07	3.0E−03	0.73
17079513	RNU6-748P	0.15	1.3E−01	0.26	1.5E−02	0.20	0.07	6.4E−03	0.78
16837063	SNORA76C	0.20	2.5E−02	0.19	1.5E−01	0.20	0.08	9.2E−03	0.83
16755265	MIR331	0.15	2.1E−02	0.26	7.0E−04	0.20	0.05	6.4E−05	0.47
16698226	LINC01136	0.17	3.0E−02	0.30	4.9E−02	0.20	0.07	5.4E−03	0.76
16738477	RP11-872D17.8	0.22	2.6E−03	0.12	3.2E−01	0.19	0.06	2.3E−03	0.72
16797869	HERC2P7	− 0.17	4.7E−02	− 0.22	5.4E−02	− 0.19	0.07	7.0E−03	0.79
16778413	TPTE2P5	− 0.25	4.9E−04	− 0.01	9.2E−01	− 0.19	0.06	2.3E−03	0.72
16865104	MIR519A2	− 0.28	3.5E−05	− 0.07	3.3E−01	− 0.19	0.05	1.8E−04	0.70
17113511	RNU6-154P	− 0.22	1.9E−03	− 0.13	2.6E−01	− 0.20	0.06	1.4E−03	0.72
16810171	RP11-507B12.1	− 0.23	3.1E−03	− 0.16	1.8E−01	− 0.21	0.07	1.6E−03	0.72
17062251	RNA5SP241	− 0.15	5.8E−02	− 0.45	3.4E−03	− 0.21	0.07	2.6E−03	0.72
16788697	SNORD114-28	− 0.29	9.1E−03	− 0.17	2.8E−01	− 0.25	0.09	6.7E−03	0.78
16928422	CTA-390C10.10	− 0.17	9.8E−02	− 0.55	4.3E−04	− 0.29	0.09	7.2E−04	0.72
16733787	RP11-304M2.6	− 0.28	1.5E−02	− 0.40	5.3E−02	− 0.31	0.10	2.6E−03	0.72

ASD autism spectrum disorder, FC fold change, SE standard error

Several differentially expressed genes in cord blood overlapped with genes previously associated with ASD in genetic or gene expression studies, although the overlap was not statistically significant (FDR q > 0.05, Additional file 1: Figure S9, Additional file 2: Table S2). Specifically, SLC7A3, VSIG4, and MIR1226 have been associated with ASD in genetic studies [42, 43], while SNORD104, OR2AG2, and DHX30 were differentially expressed in the same direction in ASD in gene expression studies [8, 44]. Further, GFI1, GPR171, KIR2DL4, PTGIR, and TRPM6 were differentially expressed in ASD cord blood and are also differentially expressed in specific blood cell types, including natural killer cells and T cells, although a significant enrichment was not observed (q > 0.05, Additional file 1: Figure S10, Additional file 2: Table S3) [41].

Overrepresentation enrichment analysis, which looks for overlap of only differentially expressed genes with biologically predefined gene lists, revealed that ASD differential transcripts were significantly enriched for functions in the response to toxic substances (TRPM6, CYP1A1, FOS, GCH1, AC012476.1, RAD51, and AQP10, fold enrichment = 9.5, q = 0.027) and ultraviolet radiation (CDO1, CYP1A1, FOS, and GCH1, fold enrichment = 7.6, q = 0.037, Fig. 1c, Additional file 2: Table S4). Both of these functional enrichments included the genes CYP1A1, FOS, and GCH1. Additionally, downregulated transcripts were enriched for functioning in blood coagulation (GNG12, MAFF, PF4, and PLG, fold enrichment = 12.5, q = 0.009) and xenobiotic metabolism (CDO1, CYP1A1, GCH1, and PLG, fold enrichment = 8.6, q = 0.019), but no significant enrichments were observed for upregulated transcripts alone.

Using fold changes to rank all transcripts for gene set enrichment analysis (GSEA), we observed significant enrichment for upregulation of gene sets involved in chromatin regulation (q < 0.05, Fig. 2, Additional file 2: Table S5). In other words, genes associated with chromatin regulation tended to be ranked toward the top of the distribution of fold change in ASD cord blood. Chromatin gene sets upregulated in ASD included DNA methylation (23 leading edge (LE) transcripts, normalized enrichment score (NES) = 2.16, q = 0.009), condensation of prophase chromosomes (24 LE transcripts, NES = 2.11, q = 0.007), nucleosome assembly (24 LE transcripts, NES = 1.96, q = 0.026), histone deacetylase (HDAC)-mediated deacetylation (30 LE transcripts, NES = 1.90, q = 0.040), and polycomb-repressive complex 2 (PRC2)-mediated methylation (22 LE transcripts, NES = 1.89, q = 0.002). Additionally, the gene set for the autoimmune disease systemic lupus erythematosus was significantly upregulated (45 LE transcripts, NES = 2.13, q = 0.003). Most of the genes associated with these sets compose a cluster of histone genes located at the 6p22.2 locus, which was also enriched (27 LE transcripts, NES = 2.15, q = 0.007). The above findings of differential expression across two prospective cohorts suggest transcriptional dysregulation in environmentally responsive genes is present at birth in cord blood of high-risk subjects later diagnosed with ASD.

Fig. 2

Chromatin and autoimmune gene sets are upregulated in cord blood from ASD subjects. Meta-analysis log2(fold change) for ASD versus TD gene expression was used to rank probes for gene set enrichment analysis (GSEA) with WebGestalt. GSEA assesses whether genes in biologically defined sets occur toward the top or bottom of a ranked list more than expected by chance. Significantly enriched gene sets are shown (FDR q < 0.05)

Non-TD-associated differential gene expression in cord blood

To assess the specificity of ASD-associated transcriptional differences in cord blood, we also examined differential expression between cord blood samples from infants later classified as Non-TD compared to TD at 36 months. Meta-analysis results showed no transcripts differentially expressed at a conservative FDR q value < 0.05. Under the thresholds of log2(fold change) > 0.1 and nominal p value < 0.01, 66 transcripts were differential, with 38 upregulated and 28 downregulated (Non-TD n = 92, TD n = 120, Fig. 3a, Additional file 2: Table S6). The median absolute log2(fold change) was 0.12. The gene with the greatest fold change between Non-TD and TD subjects was TAS2R46 (log2(fold change) = 0.37, standard error = 0.12, Fig. 3b, Table 4). Further, the estimated fold changes of Non-TD-associated differentially expressed genes were highly correlated between the individual studies (Pearson’s r = 0.80, p = 3.9E−16); however, fold changes of all transcripts were weakly correlated (Pearson’s r = 0.01, p = 0.10, Additional file 1: Figure S7b). Additionally, median expression of differentially expressed genes was not different from other genes on the array (MARBLES: differential = 4.48, non-differential = 4.64, p = 0.65; EARLI: differential = 4.15, non-differential = 4.20, p = 0.90; Additional file 1: Figure S11).

Fig. 3

Identification and functional enrichment of Non-TD-associated differentially expressed genes in cord blood. Gene expression in umbilical cord blood samples from subjects with typical development (TD, n = 120, 59 male/61 female) or those diagnosed with Non-TD at age 3 (n = 92, 50 male/42 female) was assessed by expression microarray. SVA was performed to control for technical and biological variables including sex and array batch. Differential expression analysis was carried out separately by study and combined in a meta-analysis. a Identification of 66 differentially expressed genes in meta-analysis (66 probes, log2(fold change) > 0.1, p < 0.01). b Fold change in gene expression for top 20 differentially expressed genes by meta-analysis log2(fold change) for individual studies and meta-analysis. c Non-TD-associated differentially expressed genes were analyzed for functional enrichment in recently evolved genes with WebGestalt using the hypergeometric test. Significantly enriched clades are shown (FDR q < 0.05). d Non-TD-associated differentially expressed genes specific to primates are listed

Table 4

Top 20 Non-TD-associated differentially expressed genes by log2(fold change) in meta-analysis

Probe ID	Gene symbol	MARBLES		EARLI		Meta-analysis
		log2(FC)	p value	log2(FC)	p value	log2(FC)	SE	p value	q value
16761521	TAS2R46	0.43	6.7E−03	0.29	9.3E−02	0.37	0.12	1.7E−03	0.78
16846019	FAM215A	0.08	4.5E−01	0.32	9.5E−04	0.22	0.07	2.5E−03	0.85
17063807	TRBV11-2	0.30	9.4E−04	0.14	9.3E−02	0.21	0.06	4.5E−04	0.74
17094496	MIR1299	0.17	7.7E−02	0.22	2.1E−02	0.20	0.07	3.8E−03	0.91
16842217	SNORD3D	0.13	7.0E−02	0.25	2.9E−03	0.18	0.05	8.3E−04	0.74
16818580	RNA5SP415	0.16	6.8E−02	0.20	5.1E−02	0.18	0.07	7.9E−03	0.91
16846016	LRRC37A2	0.09	3.5E−01	0.24	6.4E−03	0.17	0.06	7.6E−03	0.91
16902656	RNA5SP103	0.22	3.0E−03	0.06	5.2E−01	0.16	0.06	5.5E−03	0.91
16975969	RNU6-310P	0.20	1.6E−03	0.08	2.5E−01	0.15	0.05	1.6E−03	0.78
16845957	LRRC37A4P	0.21	9.3E−03	0.11	1.1E−01	0.15	0.05	3.8E−03	0.91
16788576	SNORD113-3	0.14	1.1E−02	0.15	7.2E−03	0.14	0.04	2.0E−04	0.74
16840210	SLC52A1	0.19	3.4E−05	0.08	8.6E−02	0.14	0.03	1.7E−05	0.34
16774442	SMIM2-AS1	0.14	1.6E−03	0.14	4.8E−03	0.14	0.03	1.9E−05	0.34
16985844	PMCHL2	0.11	2.9E−02	0.17	1.4E−02	0.14	0.04	1.3E−03	0.74
16805216	ASB9P1	− 0.15	1.6E−02	− 0.13	3.6E−02	− 0.14	0.04	1.4E−03	0.74
17078751	REXO1L1P	− 0.14	2.3E−03	− 0.12	1.0E−01	− 0.14	0.04	5.2E−04	0.74
16728518	FAM86C1	− 0.21	6.6E−03	− 0.08	3.1E−01	− 0.15	0.06	7.8E−03	0.91
17061129	RASA4	− 0.20	6.3E−03	− 0.09	2.7E−01	− 0.16	0.06	5.0E−03	0.91
16836858	RNU6-446P	− 0.26	3.2E−03	− 0.06	5.7E−01	− 0.18	0.07	7.1E−03	0.91
16733848	RP11-326C3.7	− 0.07	4.4E−01	− 0.33	9.6E−04	− 0.20	0.07	3.8E−03	0.91

Non-TD non-typically developing, FC fold change, SE standard error

Several of the 66 nominally differentially expressed genes between Non-TD and TD cord blood samples have been previously associated with genetic variation or gene expression in ASD, although the overlap was not statistically significant (q > 0.05, Additional file 1: Figure S9, Additional file 2: Table S2). Genetic deficiencies in MIR4269 have been previously associated with a reduced risk for ASD [43], while DHCR24, GNAO1, and TYMS were differentially expressed in ASD in other studies [8, 44]. Additionally, none of the Non-TD differentially expressed genes were known cell type-specific genes (Additional file 1: Figure S10, Additional file 2: Table S3) [41]. Differentially expressed genes in Non-TD that overlap genes previously associated with ASD likely function in general neurodevelopment.

Because genes recently evolved in primates have been hypothesized to play a role in human neurodevelopment, differentially expressed genes in Non-TD cord blood were assessed for enrichment in recently evolved genes by vertebrate lineages, ranging from tetrapods to homininae using overrepresentation enrichment analysis [48]. Non-TD-associated genes were significantly enriched for genes recently evolved in mammalia, theria, eutheria, boreoeutheria, euarchontoglires, primates, catarrhini, and hominidae, with the greatest enrichment in primate-specific genes (fold enrichment = 7.5, q = 2.1E−5, Fig. 3c, Additional file 2: Table S7). Of genes that evolved in primates, SLC52A1, SPANXN5, and TRIM49B were upregulated in Non-TD cord blood, while FAM86C1, RASA4, RASA4B, and TRIM73 were downregulated (Fig. 3d). In contrast, ASD differentially expressed genes were not significantly enriched in recently evolved genes from any of the vertebrate lineages (q > 0.05).

After GSEA with all probes ranked by fold change in Non-TD compared to TD subjects, we observed significant enrichment for upregulation of sensory perception gene sets (q < 0.05, Fig. 4a, Additional file 2: Table S8). Taste receptor activity (14 LE transcripts, NES = 2.30, q < 1E−4), metabotropic glutamate receptors (17 LE transcripts, NES = 2.21, q = 4.9E−3), and olfactory receptor activity (105 LE transcripts, NES = 1.96, q = 0.018) gene sets were all upregulated in cord blood from Non-TD subjects. Additionally, gene sets that interact with the compounds quinine (19 LE transcripts, NES = 2.30, q = 2.0E−3) and citric acid (22 LE transcripts, NES = 2.17, q = 2.5E−3) were significantly upregulated, while those interacting with indomethacin (18 LE transcripts, NES = − 2.02, q = 0.037) and H2-receptor antagonists (6 LE transcripts, NES = − 2.03, q = 0.047) were downregulated. Taste receptor genes included in these gene sets and the top Non-TD-upregulated gene, TAS2R46, are located at the 12p13.2 locus, which was also enriched (16 LE transcripts, NES = 2.11, q = 8.3E−3).

Fig. 4

Sensory perception gene sets are dysregulated in Non-TD subject cord blood. Meta-analysis log2(fold change) for Non-TD versus TD gene expression was used to rank probes for GSEA with WebGestalt. GSEA assesses whether genes in biologically defined sets occur toward the top or bottom of a ranked list more than expected by chance. Normalized enrichment score was plotted for significantly enriched gene sets (FDR q < 0.05) for default WebGestalt databases

Comparison of ASD and Non-TD differentially expressed genes

Eight genes were differentially expressed in both ASD and Non-TD compared to TD subjects, which was more than expected by chance (odds ratio = 28.3, p = 1.67E−9, Fig. 5a). Specifically, IGLV1-40, LRRC37A4P, MIR1299, PMCHL2, and TRBV11-2 were upregulated, while RNU4ATAC11P, TPTE2P5, and TRIM74 were downregulated in both ASD and Non-TD subjects (Fig. 5b). LRRC37AP, MIR1299, PMCHL2, and TRBV11-2 were among the top upregulated genes in Non-TD subjects (Fig. 3b). Additionally, the fold changes across all transcripts were moderately correlated between the ASD versus TD and Non-TD versus TD comparisons both within study and in the meta-analysis (meta-analysis log2(fold change) Pearson’s r = 0.38, p < 2.2E−16, Additional file 1: Figure S7c-e). These findings suggest that some ASD-associated transcriptional alterations in cord blood are also present in Non-TD subjects.

Fig. 5

A subset of ASD-associated differentially expressed genes are also differentially expressed in Non-TD subjects. a Overlap of ASD- and Non-TD-associated differentially expressed genes from meta-analysis by direction. Significance was assessed by Fisher’s exact test. b Meta-analysis log2(fold change) in gene expression for ASD- and Non-TD-associated differentially expressed genes sorted by average log2(fold change)

Coexpression network analysis and cell type deconvolution in cord blood from high-risk children

As a complementary approach to differential gene expression analysis, we performed WGCNA to identify consensus gene coexpression modules and their correlations with ASD or Non-TD diagnosis and other variables. Seventy-nine consensus coexpression modules were identified, which ranged in size from 20 to 4137 transcripts (Additional file 1: Figure S12, Additional file 2: Table S9, S10). Overall the module eigengene networks were highly preserved between MARBLES and EARLI subjects, indicating the identified consensus modules are robust (overall preservation = 0.93, Additional file 1: Figure S13). Module eigengenes were correlated with diagnosis and demographic factors within each study, and these results were combined in meta-analysis (Fig. 6, Additional file 1: Figures S14-S16, Additional file 2: Table S11). Across MARBLES and EARLI subjects, modules were identified that significantly correlated with diagnostic group, sex, gestational age, birth weight, ethnicity, paternal age, delivery method, and maternal smoking (FDR q < 0.05, Fig. 6, Additional file 1: Figure S16). In particular, gestational age, birth weight, and paternal age were associated with more than 20 modules each, suggesting that these factors have major effects on gene expression in cord blood from high-risk children. The skyblue1 module, which mostly includes genes on chromosome Y, was significantly upregulated in ASD compared to TD subjects (Z-score = 3.4, FDR q = 0.046, Additional file 1: Figure S17a-b, Additional file 2: Table S11). Two of these genes, TTTY10 and ZFY, were nominally upregulated in ASD compared to TD subjects in this study, which is a significant enrichment (Additional file 1: Figure S17d, odds ratio = 33.6, p = 0.002). However, the association of skyblue1 with ASD is likely driven by increased expression of skyblue1 genes in males and a greater frequency of males with ASD. Skyblue1 was significantly correlated with sex (Z-score = − 28.1, FDR q = 2.7E−172), paternal age (Z-score = 2.2, FDR q = 0.047), and maternal smoking (Z-score = − 3.4, FDR q = 0.015). Probes in skyblue1 map to 21 genes, 18 of which are located only on chromosome Y, explaining the strong upregulation in males (Additional file 1: Figure S17c).

Fig. 6

Consensus modules are correlated with diagnosis and demographic factors in meta-analysis. Heatmap of meta-analysis biweight midcorrelation Z-scores of module eigengenes with sample covariates (MARBLES: ASD n = 41, Non-TD n = 44, TD n = 76; EARLI: ASD n = 19, Non-TD n = 47, TD n = 43). Z-scores from the individual studies were combined using Stouffer’s method with weights given by the square root of sample n. p values were adjusted for all comparisons shown in Additional file 1: Figure S16 using the FDR method (*q < 0.05). Modules and covariates with significant correlations are shown

Because cord blood is a heterogeneous tissue composed of many cell types with distinct transcriptional profiles, we examined the proportions of 22 cell types estimated with the CIBERSORT web tool and their correlation with diagnostic group and demographic variables. The most prevalent cell types detected in cord blood overall were neutrophils and naïve CD4+ T cells, which made up 26% and 24% of whole cord blood, respectively (Additional file 2: Table S12). When cell type proportions were compared to diagnostic group and demographic factors within each study and combined in meta-analysis, no cell types differed between diagnostic groups (Fig. 7a). However, cell fractions were significantly correlated with demographic factors including sex, gestational age, birth weight, delivery method, and maternal smoking (FDR q < 0.05, Fig. 7b, Additional file 2: Table S13). To determine the correspondence of consensus modules with the proportions of specific cell types, the module eigengenes were compared to the cell type fractions followed by meta-analysis. The eigengenes for the majority of consensus modules were strongly and significantly correlated with the proportion of multiple cell types (FDR q < 0.05, Additional file 1: Figure S18, Additional file 2: Table S14). The cell types significantly correlated with the most modules were neutrophils, naïve CD4+ T cells, CD8+ T cells, M0 macrophages, and plasma cells, which were associated with more than 30 modules each. In contrast, skyblue1, which was correlated with diagnosis, along with sex, paternal age, and maternal smoking, was not correlated with the proportions of any examined cell type (Additional file 1: Figures S16 and S18, Additional file 2: Table S14). As an example, the grey60 module was positively correlated with the proportions of naïve B cells, plasma cells, and activated dendritic cells, but negatively correlated with resting NK cells and neutrophils (FDR q < 0.05, Additional file 1: Figure S18). The hub gene for grey60 was CD79A, which encodes Igα, the B cell antigen receptor component necessary for signal transduction [52]. Interestingly, the grey60 module eigengene, CD79A expression, and the proportion of naïve B cells were all significantly upregulated during cesarean delivery (FDR q < 0.05) and nominally downregulated in ASD compared to Non-TD subjects (p < 0.05, Additional file 1: Figure S19, Additional file 2: Tables S11 and S13). These results from coexpression and cell type deconvolution analyses suggest that biological factors including cell type, gestational age, birth weight, and paternal age are major drivers of interindividual variation in cord blood gene expression.

Fig. 7

Cell type proportions are correlated with demographic factors in meta-analysis. a Barplot of mean estimated proportion of each cell type by diagnosis and study from CIBERSORT using default adult peripheral blood reference panel. b Heatmap of meta-analysis biweight midcorrelation Z-scores of cell type proportions with sample covariates (MARBLES: ASD n = 41, Non-TD n = 44, TD n = 76; EARLI: ASD n = 19, Non-TD n = 47, TD n = 43). Z-scores from the individual studies were combined using Stouffer’s method with weights given by the square root of sample n. p values were adjusted using the FDR method (*q < 0.05)

Discussion

Perinatal transcriptional alterations in ASD

Based on meta-analysis across two high-risk pregnancy cohorts of similar design, we examined gene expression differences in cord blood between high-risk children who went on to develop ASD, were Non-TD, or were TD at 36 months. Significant differential gene expression in cord blood at individual genes was not observed in individuals who developed ASD compared to TD after adjusting for multiple comparisons. However, significant gene set enrichment was seen in toxic substance response, blood coagulation, chromatin regulation, and autoimmune response functions among genes differentially expressed at a nominal level in ASD.

Environmental factors are thought to contribute to ASD risk, especially during the perinatal period, a sensitive window for neurodevelopment [53]. Genes nominally differentially expressed in cord blood from ASD subjects, including CYP1A1 and GCH1, were significantly enriched for functions in xenobiotic metabolism and response to both toxic substances and ultraviolet radiation. Notably, CYP1A1 was downregulated in ASD cord blood and has been previously found to be transcriptionally regulated in blood by toxicants that affect neurodevelopment, including polychlorinated biphenyls [54-56]. GCH1—which is the rate-limiting enzyme in the synthesis of tetrahydrobiopterin, a precursor to folate, dopamine, and serotonin [57]—was also downregulated in cord blood from ASD subjects. GCH1 expression increases in response to valproic acid, an anti-epileptic drug associated with behavioral deficits in mice and increased risk of autism in humans after in utero exposure [58, 59]. Interestingly, GCH1 is genetically associated with ASD subphenotypes and is downregulated in peripheral blood from ASD children, and its product tetrahydrobiopterin is decreased in cerebrospinal fluid from ASD subjects [17, 60, 61]. These environmentally responsive genes may have altered expression due to increased genetic susceptibility and/or early life exposure to toxicants in patients with ASD [62, 63]. It is important to note that the enrichment analysis identified gene sets with significant overrepresentation, and so the differential expression of individual genes should be interpreted cautiously.

Epigenetic modifications, such as those to DNA and histone protein subunits, are affected by both genetic and environmental factors and are thought to play a role in mediating ASD risk [64, 65]. Immune dysregulation has also been found in children with ASD, and immune cells rely on epigenetic regulation for lineage commitment and cell activation in response to infection [66, 67]. A cluster of histone genes at 6p22.2 was significantly enriched for upregulated genes in ASD cord blood. Genes associated with the autoimmune disease systemic lupus erythematosus (SLE) were also upregulated, including histone-encoding, complement pathway, and antigen presentation genes. Epigenetic dysregulation is a feature of SLE, including global DNA hypomethylation, histone H3 and H4 hypoacetylation, and H3K9 hypomethylation in CD4+ T cells [68-70]. Notably, maternal SLE increases risk for ASD in offspring, suggesting an association between these two disorders [71]. Together, this implicates both epigenetic and immune mechanisms in ASD pathobiology.

Perinatal transcriptional alterations in Non-TD

To assess the specificity of cord blood gene expression changes in ASD compared to other neurodevelopmental disorders, this analysis examined transcriptional differences in Non-TD compared to TD subjects across the two studies. While no single genes were significantly differentially expressed in Non-TD subjects after adjusting for multiple comparisons, sets of functionally related genes were enriched among nominally differentially expressed genes. Significantly enriched gene sets included sensory perception and primate-specific genes. The top upregulated gene in Non-TD cord blood was TAS2R46, encoding a member of the taste 2 receptor (TAS2R) family, which was included in the top upregulated gene set, taste receptor activity. Upregulated genes in this gene set were primarily other TAS2Rs located at the 12p13.2 locus. TAS2Rs are G protein-coupled receptors (GPCRs) highly expressed in taste receptor cells and are associated with the perception of bitter taste [72]. Interestingly, individuals with attention-deficit/hyperactivity disorder and epilepsy have previously been found to perceive bitter tastes more intensely than healthy controls [73, 74]. TAS2Rs are also expressed in blood leukocytes, where they function in chemosensation of foodborne flavor compounds and response to food uptake [75]. Further, TAS2Rs are upregulated in leukocytes from asthma patients and levels of lipopolysaccharide-induced pro-inflammatory cytokines are decreased by TAS2R agonists [76]. Taste receptor upregulation may reflect altered chemosensation in the immune and nervous systems in Non-TD subjects.

Differentially expressed genes in cord blood from Non-TD subjects included genes recently evolved in primates and genes that function in neurodevelopment. Primate-specific genes originated at a similar evolutionary time that the neocortex expanded and have biased expression for the fetal neocortex in humans [48]. RASA4 is a fetal-biased primate-specific gene that was also one of the top downregulated genes in Non-TD subjects. RASA4 is a GTPase-activating protein in the Ras signaling pathway that functions in the activation of T cells, mast cells, and macrophages [77-79]. Children with Non-TD were also observed to have downregulation of GNAO1, encoding a G protein alpha subunit important for neurodevelopment and synaptic signaling. Mutations in GNAO1 are associated with epileptic encephalopathy, involuntary movements, and intellectual disability [80, 81]. Additionally, missense mutations and downregulation of GNAO1 in lymphocytes occur in individuals with schizophrenia [82, 83]. In individuals with ASD, GNAO1 is upregulated in postmortem cortex [6]. Further, GNAO1 is required in mast cells for toll-like receptor 2-mediated pro-inflammatory cytokine release, suggesting GNAO1 functions in cell signaling in both the nervous and immune systems [84].

Cord blood as a window into transcriptional alterations specific to ASD

Umbilical cord blood gene expression offers a unique snapshot of molecular differences in perinatal life, a critical window for neurodevelopment [85]. Hematopoietic cells in the blood are actively differentiating and responding to environmental cues, such as pathogens and xenobiotics [86, 87]. Epigenetic marks written during this period, which reflect short-term transcriptional activity, have the potential to have long-term effects on gene regulation and cell function [88, 89]. Signals from the immune system cross the blood-brain barrier during gestation and influence the development of the brain [90]. Toxicant exposure during gestation can also impact brain development [91, 92]. In this study, genes involved in toxic substance response, xenobiotic metabolism, and chromatin regulation were altered in cord blood from subjects diagnosed with ASD at 36 months. Transcriptional differences in cord blood from ASD and Non-TD subjects compared to TD subjects were largely independent, with only 8 genes in common. Enriched gene sets associated with Non-TD expression included sensory perception and primate-specific genes and did not overlap with ASD expression gene sets. Further, genes associated with ASD in previous studies of genetic variation and gene expression had few overlaps with ASD-associated genes in cord blood [6, 8, 42–44]. Instead, cord blood ASD genes likely represent tissue-specific environmentally responsive genes that may reflect in utero exposures and long-term altered neurodevelopmental outcomes.

Major factors contributing to transcriptional variability in whole cord blood

While the purpose of this study was to identify changes in cord blood gene expression associated with ASD and Non-TD diagnosis, we also identified biological factors associated with transcription variation in cord blood, including gestational age, birth weight, paternal age, and cell type composition. Gestational age was correlated with many consensus gene coexpression modules in cord blood, and it likely reflects ongoing differentiation and expansion of hematopoietic cells captured at different developmental timepoints. Specifically, gestational age was positively correlated with the estimated proportion of neutrophils and negatively correlated with the proportions of CD8+ T cells, naïve CD4+ T cells, and resting mast cells, suggesting alterations in these cell types are contributing to variability in gene expression related to gestational age. Notably, birth weight was correlated with many of the same gene coexpression modules as gestational age and is itself correlated with gestational age, so modules associated with birth weight likely reflect gestational age as well.

A mostly independent set of coexpression modules were correlated with paternal age. Increased paternal age has previously been associated with global changes in transcription in peripheral blood lymphocytes, including a downregulation of genes involved in transcriptional regulation and a decrease in the distribution of gene expression variance [93]. Interestingly, both increased paternal age and decreased gene expression variance have also been associated with ASD [93, 94]. The mechanism for the effect of paternal age on transcription is still unclear, but previous studies have observed increases in telomere length and de novo mutations and alterations in DNA methylation at specific genes in offspring from older fathers [95-97].

The factor with the strongest association with cord blood gene expression was cell type composition, as observed previously [98]. The cell types whose proportions correlated with the most modules were neutrophils, naïve CD4+ T cells, CD8+ T cells, M0 macrophages, and plasma cells. Neutrophils and naïve CD4+ T cells were also the most prevalent cell types in cord blood. The strong correlation of estimated cell type proportions with coexpression modules suggests that variability between samples in cell type composition is contributing a large portion of the transcriptional signal measured in whole cord blood. Overall, the large impact of biological factors on gene expression in a bulk tissue such as cord blood makes detecting differentially expressed genes associated with a heterogeneous disorder such as ASD practically difficult. Future studies should take care to reduce these sources of noise, by isolating specific cell populations, selecting subjects with similar gestational and paternal age, or focusing on more narrow ASD endophenotypes [99].

Limitations

Cord blood samples from two enriched autism risk pregnancy cohort studies were analyzed together to improve the power to detect differential gene expression between ASD and TD individuals. Meta-analysis substantially increased the total number of ASD and Non-TD subjects; however, the study was not adequately powered to detect moderate differences at a single gene level after correcting for multiple comparisons. We have included a supplemental power calculation for future studies based on our findings (Additional file 1: Figure S20). Differential genes used in the enrichment analysis also did not meet multiple comparison-corrected statistical significance. Lack of statistical evidence of individual differential genes does not eliminate the potential of identifying a biologically significant gene set enrichment across diagnostic groups; however, the results should be interpreted cautiously. ASD is a heterogeneous disorder, and this may mask differential expression. The two cohort studies coordinated classifying participants into ASD, Non-TD, and TD diagnostic groups using the same metrics to improve consistency between the studies. Nonetheless, it is possible that heterogeneity within ASD remains in this study without further breaking down the group into subtypes.

The two cohorts used the same platform to measure the RNA expression levels and were subject to the limitations in transcript coverage and measurement precision of the microarray platform. Transcripts not covered by a probe on the array were not analyzed in this study, and targeted quantitative analysis in the general population would be needed to validate specific transcriptional changes as ASD risk biomarkers. A cord blood cell type gene expression reference panel was not available at the time of analysis, and so samples were deconvoluted using an adult panel. Future studies should incorporate a cord blood-specific cell type deconvolution method to improve accuracy. Additionally, genetic, epigenetic, and environmental variation is known to impact both gene expression and ASD risk, but this was not investigated in this study. Future studies that integrate cord blood gene expression with genetic, epigenetic, and environmental factors will be important to improve understanding of ASD etiology.

Conclusions

In the first study to investigate gene expression in cord blood from high-risk newborns later diagnosed with ASD, we identified nominally statistically significant transcriptional alterations specific to ASD, which were enriched for toxic substance response and epigenetic regulation functions. Differentially expressed genes in ASD had few overlaps with those identified in cord blood from newborns with other non-typical neurodevelopmental outcomes in this high-risk population. Instead, Non-TD-associated genes were enriched for sensory perception functions and primate-specific genes. Further, coexpression and cell type analyses revealed that gestational age, birth weight, paternal age, and cell type composition have large impacts on variability in gene expression measured in whole cord blood.

A strength of this high-risk prospective pregnancy cohort design was the observation of gene expression at birth, prior to the onset of symptoms, diagnosis, and treatment. Perinatal life is a critical period for neurodevelopment, where environmental stressors could have long-term impact. Additionally, ASD-associated differential expression was meta-analyzed across two independent studies with harmonized methods and compared with expression changes in other subjects with non-typical neurodevelopment. Finally, cord blood is an accessible tissue that reflects the perinatal environment, and ASD-associated gene expression changes in cord blood may have potential as a predictive biomarker.

Additional file 1: Figure S1. Overview of study design. Figure S2. Surrogate variable analysis in MARBLES subjects. Figure S3. Surrogate variable analysis in EARLI subjects for the ASD versus TD comparison. Figure S4. Surrogate variable analysis in EARLI subjects for the Non-TD versus TD comparison. Figure S5. Identification of ASD-associated differentially-expressed genes in cord blood within each study. Figure S6. Identification of Non-TD-associated differentially-expressed genes in cord blood within each study. Figure S7. Correlations between ASD and Non-TD expression differences in MARBLES and EARLI subjects. Figure S8. Expression level distribution of meta-analysis ASD versus TD differential probes is similar to non-differential probes. Figure S9. Cord blood differentially-expressed genes are not enriched for ASD-associated gene sets. Figure S10. Cord blood differentially-expressed genes are depleted for blood cell-specific genes. Figure S11. Expression level distribution of meta-analysis Non-TD versus TD differential probes is similar to non-differential probes. Figure S12. Consensus coexpression modules identified in MARBLES and EARLI. Figure S13. Consensus module eigengene networks are preserved between MARBLES and EARLI subjects. Figure S14. Consensus modules are correlated with diagnosis and demographic factors in MARBLES subjects. Figure S15. Consensus modules are correlated with demographic factors in EARLI subjects. Figure S16. Consensus modules are correlated with diagnosis and demographic factors in meta-analysis. Figure S17. Skyblue1 module is specifically expressed in males and is enriched for genes upregulated in ASD. Figure S18. Consensus modules are strongly correlated with cell type proportions in meta-analysis. Figure S19. The B cell-associated grey60 module is upregulated during cesarean delivery and is downregulated in ASD. Figure S20. Anticipated power curve calculated by ssize.fdr R package with projected proportions of non-differentially expressed genes as 0.8, 0.85, 0.9, 0.95 and 0.98.

Additional file 2: Table S1. Differential expression analysis of ASD compared to TD subjects for individual studies and meta-analysis. Table S2. Gene overlap analysis of ASD- and Non-TD-associated differentially-expressed genes with genes previously associated with ASD. Table S3. Gene overlap analysis of ASD- and Non-TD-associated differentially-expressed genes with blood cell-type specific genes. Table S4. Significant terms from overrepresentation enrichment analysis of ASD-associated differentially-expressed genes. Table S5. Significant terms from GSEA of probes ranked by ASD versus TD differential expression. Table S6. Differential expression analysis of Non-TD compared to TD subjects for individual studies and meta-analysis. Table S7. Significant terms from overrepresentation enrichment analysis of Non-TD-associated differentially-expressed genes. Table S8. Significant Terms from GSEA of probes ranked by Non-TD versus TD differential expression. Table S9. Consensus coexpression module membership in MARBLES subjects. Table S10. Consensus coexpression module membership in EARLI subjects. Table S11. Correlation meta-analysis of consensus module eigengenes with demographic factors. Table S12. Cell type proportions in MARBLES and EARLI subjects estimated with CIBERSORT. Table S13. Correlation meta-analysis of cell type proportions with demographic factors. Table S14. Correlation meta-analysis of consensus module eigengenes with cell type proportions.