Identification of aberrant innate and adaptive immunity based on changes in global gene expression in the blood of adults with autism spectrum disorder

Autism spectrum disorder (ASD) is characterized as a neurodevelopmental disorder, and one of the main hypotheses regarding its cause is genetic factors. A previous meta-analysis of seven microarray studies and one RNA sequencing (RNA-seq) study using the blood of children with ASD identified dysregulation of gene expressions relevant to the immune system. In this study, we explored changes in global gene expression as the phenotype of ASD in the blood of adults with ASD. We recruited an RNA-seq cohort (ASD vs. control; n = 6 each) and a replication cohort (ASD vs. control; n = 19 each) and conducted RNA-seq to explore changes in global gene expression. We then subjected the significantly up- and downregulated genes to gene ontology (GO) and core analyses. Weighted gene correlation network analysis (WGCNA) was performed with all 11,617 genes detected in RNA-seq to identify the ASD-specific gene network. In total, 117 significantly up- and 83 significantly downregulated genes were detected in the ASD compared with the control group, respectively (p < 0.05 and q < 0.05). GO analysis revealed that the aberrant innate and adaptive immunity were more obvious in the 117 upregulated than in the 83 downregulated genes. WGCNA with core analysis revealed that one module including many immune-related genes was associated with the natural killer cell signaling pathway. In the results for the replication cohort, significant changes with same trend found in RNA-seq data were confirmed for MAFB (p = 0.046), RPSAP58 (p = 0.030), and G2MK (p = 0.004). The sample size was relatively small in both the RNA-seq and replication cohorts. This study examined the mRNA expression level, so the interaction between mRNA and protein remains unclear. The expression changes between children and adults with ASD were not compared because only adults with ASD were targeted. The dysregulated gene expressions confirmed in the blood of adults with ASD were relevant to the dysfunction of innate and adaptive immunity. These findings may aid in understanding the pathogenesis of ASD.


Background
Autism spectrum disorder (ASD) is characterized as a neurodevelopmental disorder that has three main traits: stereotyped behaviors, deficits in communication, and diminished social skills. As of 2012, the median prevalence of ASD was reported to be 0.62% [1], and a diagnosis of ASD is more frequent in males than in females [2]. However, the current prevalence has risen for several reasons, including increased recognition, awareness, and changes in clinical practice or service availability [1]. Although multiple factors, environmental toxins and stressors, mitochondrial dysfunction, and impaired immune responses are involved in the pathogenesis of ASD [3], the mainstream hypothesis regarding its cause is the individual's genetic background, such as singlenucleotide polymorphisms [4], rare genetic variants [5][6][7], and copy number variation [8][9][10]. As a nongenetic risk factor (in some cases, that interact with genetic factors), the most reliable and replicable association for ASD is advanced paternal age [11,12]. The de novo mutations in sperm increase with paternal age [13][14][15][16], and most of the de novo mutations found in individuals with ASD come from the paternal origin [17,18]. A growing body of evidence using postmortem brain studies has suggested the excessive production and increased density of microglia cells, especially in the dorsolateral prefrontal cortex [19], fronto-insular and visual cortex [20], cerebellum, midfrontal, and cingulate gyrus [21]. These findings are concordant with those from positron emission tomography studies [22]. At the morphology level, increased short-distance microglia-neuron interaction has been reported [23]. A transcriptome study revealed that these microglial changes were also found in gene expression levels, possibly induced by inflammatory cytokines [24]. Changes in inflammatory cytokines have also been found in the blood, including elevated proinflammatory cytokines such as TNF-α, interleukin (IL)-6, and IL-8, along with a decrease in anti-inflammatory cytokines such as TGF-β and IL-10 [25][26][27]. A metaanalysis of seven microarray studies using blood samples of children with ASD concluded that transcriptional cascades are typically elicited within circulating immune cells, contrary to the activated immune response in protein levels [28]. Recently, RNA sequencing (RNA-seq) has become a powerful method for investigating global gene expression as well as microarray data because it quantifies a large and dynamic range of expression levels with absolute rather than relative values. To our knowledge, only one transcriptome study has been conducted using RNA-seq of blood; the results revealed that immune dysregulation changes in gene expression levels were found based on genome-wide gene expression data using twin subjects with or without ASD [29]. Surprisingly, the one RNA-seq study and all seven microarray studies explored global gene expressions in the blood of only children with ASD. Considering that adults with ASD often experience more comorbidities as a result of facing social problems [30], gene expression among adults with ASD possibly differs from that among children with ASD because gene expressions generally change as a results of environmental factors such as smoking [31], alcohol consumption [32], and disease states. Given this background, the present study recruited adults with ASD and age-and sex-matched controls to investigate (1) global expression changes using RNA-seq of the blood, (2) the type of biological process to which differentially expressed genes (DEGs) belong based on the bio-computational method of gene ontology (GO) analysis, (3) ASD-specific gene networks using weighted gene correlation network analysis (WGCNA) with RNA-seq data, and (4) the expression profiles of selected genes in replication set using RNA-seq.

Demographic and clinical data
We recruited an RNA-seq cohort (ASD vs. control; n = 6 each) and a replication cohort (ASD vs. control; n = 19 each) and conducted RNA-seq to explore changes in global gene expression. The demographic and clinical data of the RNA-seq and replication cohorts are shown in Tables 1 and 2, respectively. A diagnosis of ASD was made according to the Diagnostic and Statistical Manual of Mental Disorders, Fifth Edition, by at least two expert psychiatrists on the basis of extensive clinical interviews and a review of medical records. The control subjects without psychiatric signs or a past history of mental disorders, and who were diagnosed as neurotypical normal based on clinical interviews by at least two expert psychiatrists, were recruited. The recruited individuals with ASD and controls were of unrelated Japanese origin and provided written informed consent before the study began in accordance with the ethics committee of Ehime University Graduate School of Medicine (No. 31-K8). The severity of ASD symptoms in the RNA-seq cohort were assessed by using the Autism Diagnostic Observation Schedule, Second Edition [33] or Autism Diagnostic Interview-Revised [34].

RNA isolation and synthesis of complementary DNA
Blood was collected into PaxGene Blood RNA Systems tubes (BD, Tokyo, Japan), and RNA was isolated according to the manufacturer's protocol. RNA concentration and quality were calculated by using the NanoDrop 1000 system (Thermo Fisher Scientific, Yokohama, Japan). RNA samples indicating a 260/280 ratio between 1.8 and 2.0 were assumed to be pure, and 1.0 μg of RNA was used to synthesize 40-μL reaction mixtures of complementary DNA (cDNA) using the High-Capacity cDNA Reverse Transcription Kit (Applied Biosystems, Foster City, CA, USA).

Quantitative PCR
Reverse transcription quantitative PCR (RT-qPCR) was conducted to measure mRNA expression levels using the StepOnePlus Real-Time PCR System (Applied Biosystems). for GAPDH. RT-qPCR was conducted by using the Prime-Time Gene Expression Master Mix (Integrated DNA Technologies, Inc., Coralville, IA, USA) in a final volume of 10 μL. mRNA expression levels were measured in duplicate, and the same sample was used in each plate to remove errors between plates. The relative expression value was calculated by using Livak's ΔΔCt method [35].

RNA-sequencing
The RNA integrity number was measured with an Agilent 2100 Bioanalyzer (Agilent Technologies Inc., Santa Clara, CA, USA) and an Agilent RNA 6000 Nano Kit (Agilent Technologies Inc.), as shown in Table 1

Bioinformatic analysis of RNA sequencing data
Raw data files in FASTQ format were generated from the MiSeq system (Illumina). The reads were aligned to the reference genome (hg38) with TopHat software [36]. The expression levels (metric fragments per kilobase of transcript per million mapped reads [FPKM value]) of the known genes were estimated using Cufflinks [37]. The number of identified genes per group was calculated based on the average FPKM values ≥ 1.0 in each group. Genes including an FPKM value of 0 in either the ASD or control group were excluded. A total of 11,617 genes were analyzed according to the criteria. DEG analysis was performed using the edgeR package [38]. Statistical significance for the DEGs was set at p < 0.05 and q < 0.05. A volcano plot of DEGs was created using DEseq2 [39]. A heatmap based on significantly changed DEGs was generated by using heatmap.2.

Weighted correlation network analysis
Weighted gene correlation networks based on the RNAseq data were generated using the "WGCNA" R package [40]. Weighted gene correlation network analysis (WGCNA) was performed to identify gene modules and summarize clusters by the module eigengene or an intramodular hub gene in consideration of external traits according to the WGCNA tutorial (https://horvath. g ene tic s. uc la .e du / ht m l / C o ex p re ss ion Net w o rk / Rpackages/WGCNA/). A "signed network adjacency" was applied to construct a network from the RNA-seq data (12 subjects: n = 6 each in the ASD and control groups) by setting a soft thresholding power of "6," which was determined according to the result of scalefree topology and mean connectivity ( Figure S1). The adjacency matrix was converted into a "topological overlap matrix (TOM)" to minimize the effects of noise and spurious associations. The "dynamicTreeCut" R package was used to detect groups of highly correlated genes.
The corresponding dissimilarity was calculated as "1 -TOM." Then, the modules were merged by setting a cut height value of 0.25′, a deep split of 2′, and a minimum module size of 30′ ( Figure S2). For the created modules, the module eigengene was calculated by using the "Mod-uleEigengenes" function. Correlation analysis between a module and each parameter (trait [ASD or control subjects], age, and sex) was performed using the "corPvalueStudent" function, with statistical significance set at p < 0.05. The genes of the interested module were then subjected to functional annotation and core analysis.

Functional annotation
The gene ontology (GO) of biological process (BP), molecular function (MF), and cellular component (CC) were conducted by using the ClueGO plugin [41] in the Cytoscape program (ver. 3.8.0). Statistical significance was set at p < 0.05 and corrected by suing the Benjamini-Hochberg method. Graphical networks of significant GO terms were created by using GO term fusion based on the following criteria: visual style = Groups; GO Term/Pathway Network Connectivity = Medium (kappa score = 0.50).

Core analysis (canonical pathway, diseases, and functions)
Ingenuity Pathway Analysis software (Qiagen, Valencia, CA, USA) was used to consider a module for the functional enrichment of target genes and decipher their role in canonical and disease pathways by using Fisher's exact test, with the p value threshold set at < 0.05.
Statistical analysis SPSS 22.0 software (IBM Japan, Tokyo, Japan) was used for the statistical analysis. The assessment of a normal distribution was considered by using the Shapiro-Wilk test. Average differences in age and mRNA levels were assessed by using Student's t test or the Mann-Whitney U test. Differences in sex were analyzed with Fisher's exact test. Correlations between gene expressions and covariates were conducted by suing the Pearson correlation coefficient or Spearman's rank correlation coefficient. The average differences of mRNA expression in qPCR by sex were tested by using Student's t test or the Mann-Whitney U test. Statistical significance was defined at the 95% level (p = 0.05).

Demographic data
No significant differences were found in the demographic data between the ASD and control subjects in the RNA-seq cohort (age: p = 1.0, sex: p = 1.0) or the replication cohort (age: p = 0.76, sex: p = 1.0).

RNA-seq
A total of 11,617 genes were used for subsequent analysis under stringent criteria, as mentioned in the "Methods" section; these genes are depicted in the volcano plot in Fig. 1a. Of these genes, 117 significantly upregulated and 83 significantly downregulated genes were detected in the ASD compared with the control group (p < 0.05 and q < 0.05). A heatmap of significantly changed DEGs is shown in Fig. 1b. The high fold-change genes are shown in Table 3, and all significantly changed genes are shown in Table S1.
Functional annotation and core analysis in the RNA-seq data All 117 upregulated genes were subjected to analysis because the fold change of 117 genes was more than 1.5. As a result, 90 BP terms reached significant levels (Table  S2). Of those BP terms, immune-related BP terms were abundant and had functionally interacted with each other (Fig. 2a). In terms of CC and MF terms, 18 CC and four MF terms reached statistical significance. The terms of ficolin-1-rich granule lumen (p = 0.00045) and ether hydrolase activity (p = 0.00012) were most significant in CC and MF, respectively. For the same reason, all 83 downregulated genes were used to conduct functional annotation. Eight BP terms, including immune response, and six CC terms were found to be significant terms, but no significant result was found in MF (Table  S3). When considering canonical pathways, several inflammatory pathways and immune responses were relevant to both up-and downregulated genes. Those immune-related pathways and responses were more abundant in upregulated than in downregulated genes (Fig. 2b). This trend was similar to the results of diseases and functions. In addition, neurological disease, nervous system development and function, and development disorder overlapped between up-and downregulated genes, but psychological disorders were found only in upregulated genes (Fig. 2c). were associated with traits (ASD or control subjects), as shown in Fig. 3a. Subsequently, we focused on the MEbrown4 module, which was the most significant (479 genes, as shown in Table S4).

Functional annotation and core analysis in the MEbrown4 module
All 479 genes were subjected to analysis. As a result, 78 BP terms reached significant levels (Table S5). Numerous immune-related BP terms functionally interacted with each other (Fig. 3b). Moreover, 11 MF terms reached significant levels, with core promoter sequencespecific DNA binding being the most significant (p = 0.022). No significant CC terms were found. Core analysis (canonical pathways, diseases, and functions) also revealed a relationship between the MEbrown4 module and immune-related terms (Fig. 3c).
qPCR validation for RNA-seq data and replication set The same samples used in RNA-seq (ASD vs. control subjects = 6 vs. 6) were tested in a validation qPCR experiment. We picked up three most significant (p = 5.00E-05) and two highest fold change genes about up and downregulated genes, respectively. Additionally, three genes associated with the immune system (TLR1, IL1B, and TNFAIP6) which were significantly upregulated in RNA-seq result were measured. The trends in expression and fold changes of the qPCR results were similar to the RNA-seq data for all 13 genes (Fig. 4a). Significant changes were found for MAFB Next, the mRNA expression levels measured by using qPCR in the replication cohort were evaluated in terms of their association with covariates such as age and sex. No changes were associated with covariates except for ALDH2 (age in ASD, p = 0.019, r = 0.534), RPSAP58 (sex in ASD, p = 0.031: male vs. female = 1.14 ± 0.44 vs. 1.63 ± 0.43), and TNFAIP6 (sex in ASD, p = 0.008: male vs female = 1.14 ± 0.56 vs. 2.11 ± 0.0.55), as shown in Table S6.

Discussion
To our knowledge, this is the first study to conduct RNA-seq using the blood of adults with ASD and unrelated control subjects. Both up-and downregulated DEGs were associated with the immune system based on the results of bio-computational analyses. Additionally, WGCNA revealed that one module (gene network) was significantly correlated with ASD pathogenesis through numerous immune-related genes.
Most of the selected genes were successfully validated in the qPCR experiment within the same samples used in RNA-seq. However, only a few had significant and identical direction changes in the replication set; the other genes showed no significant changes, although the directions were almost the same. In addition to ASD pathogenesis, individual backgrounds (e.g., lifetime event. parenting, comorbidities) could affect gene expression changes [42]. In this regard, many samples are needed to verify this point.
In addition, a larger number of immune-related BP terms was confirmed when subjecting the 117 upregulated DEGs to GO analysis compared with the 83 downregulated DEGs, and these genes also showed higher interaction with each other. Especially, many BP terms were related to innate immunity (e.g., myeloid leukocyte mediated immunity, neutrophil degranulation, mast cell degranulation), which is one of the two main types of immunity in vertebrates (the other is the adaptive immune system). Children with ASD have specific reactions toward benign factors such as common illnesses and environmental challenges [43]. These children often suffer from inflammation of the gastrointestinal mucosa and nodular hyperplasia, and in turn, present with innate immune abnormalities [14,44]. Numerous studies have reported immune dysregulation, which is reported as increased levels of pro-inflammatory cytokines secreted by peripheral blood mononuclear cells [25,45,46]. Indeed, the dysregulation of cytokine production was revealed in upregulated DEGs from both BP terms and the canonical pathway. On the other hand, the canonical pathway revealed that the upregulated DEGs were involved in the Fig. 2 Immune-related BP results and core analysis based on RNA-sequencing data. a The network was constructed by immune-related BP terms based on 117 upregulated genes (p < 0.05 and q < 0.05), which were revealed by using ClueGo in the Cytoscape program. Canonical pathways, disease, and function were analyzed using 117 up-and 83 downregulated genes (p < 0.05 and q < 0.05). Significant results for the canonical pathway (b) and disease and function (c) are shown in red (up) and blue (down) circles. Spot/circle size is the function of -log(base=10) of Fisher's exact test enrichment p value. BP biological process adaptive immune system (e.g., Th1 pathway, Th2 pathway, Th1 and Th2 activation pathway), in addition to BP terms such as CD4-positive, alpha-beta T cell activation, positive regulation of T-helper 1 type immune response, and T cell differentiation involved in the immune response. Replicative findings have shown atypical adaptive T cell responses in ASD subjects [47]. One possible reason for this is that the dysregulation of Th1 and Th2 immune responses induces the elevated pro-inflammatory cytokine TNF-α in ASD subjects. Interestingly, Ashwood et al. reported that increased TNF-α was associated with stereotypical behaviors in ASD subjects [47]. In fact, the gene expression changes found in the replication cohort were affected by age and sex among ASD subjects. This suggests that the expression change of immune-related genes could be affected by ASD symptoms, age, and sex as well as ASD pathogenesis. Compared with elevated pro-inflammatory cytokines such as TNF-α (Th1), IL-6 (Th2), GMCSF, and IL-8 (multiple immune cells), anti-inflammatory cytokines (e.g., TGF-β, IL-10) are commonly decreased in the blood [25][26][27]48]. Collectively, our results suggest that both innate and adaptive immunity are dysregulated in adults with ASD and possibly induce increased pro-inflammatory and decreased antiinflammatory cytokines.
Moreover, WGCNA identified one module (MEb-rown4, 479 genes) that is inversely correlated with ASD. These genes are functionally highly interacted with each other and relevant to the immune system, as verified by BP terms and core analyses. This result suggests that dysregulated immune-related genes are possibly affected by ASD pathogenesis together, not individually, which is revealed by WGCNA. Central nervous system (CNS) activity may generally impact immunological functions via neurotransmitters and glucocorticoids. Conversely, proinflammatory cytokines, monocytes, macrophages, and T or B lymphocytes from the peripheral immune system act on the CNS through microglial cell activation [49, Fig. 3 The results of WGCNA, gene ontology, and core analysis based on the enriched immune-related genes (MEbrown4) module. a The figure represents only the significant module eigengene correlated with phenotypic traits (ASD or control). Each row represents the module eigengene or ME (the correlation matrix of the module and sample, labeled by color). b The network was constructed by immune-related BP terms based on 479 genes involved in MEbrown4, which were revealed by using ClueGo in the Cytoscape program. c Significant results for the canonical pathway and disease and function are based on 479 genes involved in MEbrown4. ASD autism spectrum disorder 50]. Interestingly, the most significant canonical pathway was the "natural killer (NK) cell signaling" pathway, which is responsible for innate immunity. The association between ASD and NK cells has been frequently discussed [51]. A recent study suggested that immune dysregulation was implicated in the higher expression of NK cytotoxicity genes in the blood of children with ASD [52]. Of the seven genes that were highly annotated for NK activity, four (SPON2, IL2RB, PRF1, and CX3CR1) were found in this module. The authors also investigated CD56 + NK cells under both stimulated and unstimulated situations and found higher resting but reduced cytolytic activity compared with the control subjects. A more recent study involving 104 children with ASD indicated significantly lower levels of CD57 + NK cells in children with ASD despite having a normal level of CD56 + NK cells [53]. In essence, the specific gene network revealed by WGCNA is implicated in the lower activity of NK cells, even in adults with ASD.

Limitations
This study had several limitations. First, four of six ASD subjects in the RNA-seq cohort have intellectual disability, and could not take Autism Diagnostic Observation Schedule. We did not investigate how the effect of ASD severity impacts global gene expressions. Second, the sample size was relatively small in both the RNA-seq cohort (ASD vs. control; n = 6 each) and the replication cohort (ASD vs. control; n = 19 each). In this regard, a type I error may be present; however, most of the significant genes found in RNA-seq were successfully validated in the qPCR experiment. Third, most preceding studies have verified the dysregulation of the immune system from the perspective of protein levels. The present study Fig. 4 Validation qPCR in the RNA-seq and replication cohorts. a The y-axis represents the ratio of the relative expression value (ASD/Ct) in qPCR and RNA-seq, and the x-axis represents the gene names selected by the results of RNA-seq. This experiment was done in the RNA-seq cohort. b The y-axis represents the ratio of relative expression values of ASD and control subjects, and the x-axis represents the gene names selected by the results of RNA-seq. This experiment was done in the replication cohort. *p < 0.05; †p < 0.01. ASD autism spectrum disorder, Ct control explored the mRNA expression level, so the interaction between mRNA and protein remains unclear. Lastly, we targeted only adults with ASD. Basically, the three main traits of ASD are found at birth and last a lifetime [54]. From this viewpoint, the pathogenesis of ASD lasts from birth to adulthood, but several factors, such as comorbidities, lifestyle, and drugs, could affect gene expression. In the future, expression changes in both mRNA and protein levels need to be investigated in larger cohorts of children and adults with ASD.

Conclusion
The dysregulated gene expressions confirmed in the blood of adults with ASD were relevant to the dysfunction of innate and adaptive immunity, and those possibly induce the increased pro-inflammatory and decreased anti-inflammatory cytokines. These findings may aid in understanding the pathogenesis of ASD.