HCAR3 May Mediate Autism Symptoms: Evidence From RNA-Sequencing of Peripheral Blood Mononuclear Cells
Article information
Abstract
Objective
Autism spectrum disorder (ASD) is a complex neurodevelopmental disorder for which there is currently no effective treatment, and no reliable diagnostic biomarkers are available for ASD.
Methods
The present study was designed to compare the gene expression profiles in children with ASD and typically developing children and validate differentially expressed genes (DEGs) to assist in the search for the pathophysiological markers of ASD.
Results
The results demonstrated the identification of 35 common DEGs with down-regulated expression and 21 with up-regulated expression. Gene ontology and Kyoto Encyclopedia of Genes and Genomes analyses demonstrated that these DEGs were predominantly involved in signaling and immunity pathways. Our findings revealed that hydroxycarboxylic acid receptor 3 (HCAR3) and tumor necrosis factor ligand superfamily member 15 (TNFSF15) were associated with ASD-related symbols. Following validation by external and internal cohorts, HCAR3 may be identified as a risk gene for ASD.
Conclusion
Collectively, our findings indicate that some signaling-related and immune-related genes are expressed abnormally in children with ASD, and suggest that HCAR3 plays a critical role in the ASD phenotype. These findings may offer promising avenues for developing effective diagnostic biomarkers for ASD.
INTRODUCTION
Autism spectrum disorder (ASD) is a complex neurodevelopmental disorder, characterized by social deficits and restricted and repetitive behaviors. The prevalence of ASD has increased over the past few decades, with a rate of 1 in 143 children in China, and a higher rate of 1 in 36 children in the United States [1,2]. At present, the clinical screening and diagnosis of ASD are predominantly based on scale scores and clinical interviews [3]. However, these methods were constrained by their subjective nature, and the resulting diagnoses were frequently inaccurate. A previous study indicated that the majority of ASD cases were not diagnosed in children until they reached the age of three or above [4]. Moreover, current therapeutic approaches focus on educational and behavioral interventions. The high prevalence and long-term treatment required for ASD place a heavy burden on individuals, families, and society. Biomarkers facilitated not only early objective diagnosis but also the identification of the pathophysiological mechanisms that influenced the development of the disease [5]. Consequently, the exploration of potential biomarkers of ASD may provide more significant insights for early diagnosis and treatment.
Transcriptomic studies serve as a valuable tool for characterizing complex human diseases and predicting associated molecular mechanisms, with RNA sequencing (RNA-Seq) being an efficient method for sequencing the transcriptome. This technique is capable of comprehensively and rapidly detecting a substantial number of samples, thereby obtaining transcript sequence data and expression levels [6]. RNA-Seq has become a widely utilized technique for elucidating the pathogenic mechanisms underlying various neurological diseases and identifying potential transcriptional biomarkers [7,8]. Due to their widespread use in transfusion and ease of acquisition, peripheral blood mononuclear cells (PBMC) have emerged as a significant focus within the field of ASD transcriptomic research [9]. The utilization of blood-based biomarkers holds the potential to facilitate large-scale screening and enhance the diagnostic process within clinical settings.
Therefore, PBMC samples were collected from ASD children and typical development (TD) children in this study. RNA-seq was conducted on the PBMC of all samples. A differential gene expression analysis, functional annotation, and pathway enrichment analysis were performed to explore potential targets in the development of ASD. Finally, differentially expressed transcripts detected by RNA sequencing were verified by external and internal validation cohorts. The objective of this study was to provide a basis for the clinical diagnosis and treatment of ASD.
METHODS
Study participants
In this study, we recruited a training cohort (children with ASD vs. sex- and age-matched TD children, n=10 each) for RNA sequencing and a validation cohort (children with ASD vs. sex- and age-matched TD children, n=10 each) for experimental verification. The children with ASD were recruited from the Child Development and Behavior Research Center of Harbin Medical University, Harbin, China, and designated rehabilitation institutions of the provincial Disabled Persons’ Federation, Heilongjiang Province, China. A diagnosis of ASD was made according to the Diagnostic and Statistical Manual of Mental Disorders, Fifth Edition [10] by at least two expert psychiatrists based on extensive clinical interviews and a review of medical records. Participants were excluded if they had other neuropsychiatric disorders or known genetic disorders, such as epilepsy, tuberous sclerosis, intellectual disability, or fragile X syndrome. Unrelated healthy children without a history of developmental delay or other neurological disorders were randomly selected from normal kindergartens in Harbin, China, as the typically developing children. Written informed consent was obtained from all parents before the study. The study protocol was approved by the Medical Ethics Committee of the Sixth Affiliated Hospital of Harbin Medical University (LC2023-002). All experiments on human subjects were conducted in accordance with the Declaration of Helsinki.
Clinical evaluation
The following assessment tools were used in this study: autism behavior checklist (ABC), Childhood Autism Rating Scale (CARS), and Developmental Diagnostic Scale of Children Aged 0–6 Years (DDSC). Please refer to Table 1 for detailed information on patient characteristics. All of these were administered by trained nurses. 1) ABC [11] is a parent- or caregiver-reported questionnaire. A cutoff total score of >67 indicates a high probability of autism, whereas a total score between 53 and 67 indicates questionable autism. 2) CARS [12] is a scale used to identify autism and quantitatively describe the severity of the disorder. 3) DDSC [13] can comprehensively evaluate the development level of children aged less than 6 years from five domains, including gross motor movements, fine motor movements, adaptability, language skills, and social skills. The assessment results are presented as a developmental quotient (DQ), with DQ ≤69 representing mental disability, DQ 70–84 representing low intelligence (borderline), and DQ ≥85 representing normal intelligence.
Demographic characteristics of the ASD and TD groups
Behavioral severity Z-score index
The integrated behavioral Z-score method allows the normalization of individual data observed in different behavioral parameters to decrease the intrinsic variability of single tests and enhance the sensitivity and reliability of the individual phenotyping. We calculated the integrated behavioral Z-score by the severity of behavior for the scores of ABC, CARS, and DDSC.
RNA extraction
The blood samples from children with ASD and control donors were drawn and collected in sterile ethylene diamine tetraacetic acid (EDTA) tubes. PBMCs were isolated by centrifugation over Histopaque 1077 density gradient (LTS1077, TBD). The final PBMCs were resuspended in TRIzol reagent (Invitrogen, Thermo Fisher Scientific). Total RNA was isolated using TRIzol reagent (Invitrogen) and RNAprep pure blood total RNA extraction kit (DP433, Tiangen) following the manufacturer’s instructions.
RNA sequencing
RNA-seq was carried out using Illumina based on high-throughput sequencing technology. The original disembarkation data were called raw data, which contained sequence information (reads) and corresponding sequencing quality information. The reads containing joint sequences, when the N content of a single-ended sequence read exceeded 5% of the read length, and the paired-end reads were filtered out. When the number of low-quality bases (≤10) in a single-ended read exceeded 20% of the read length, the reads and the paired-end reads were filtered out. SOAPnuke (v2.1.0; BGI Group) software was used for image recognition, deconfouling, and joint removal of sequencing results. FastQC (v0.11.9; Babraham Institute) and ht-stat (v0.21.0; German Cancer Research Center) software were used to calculate the number of sequencing reads, data yield, sequencing error rate, Q20 content, Q30 content, and GC content of clean data. Bowtie2 (v2.4.1; GitHub Release) software was used to compare clean data with the eukaryotic rRNA database to identify the rRNA contamination rate of each sample. The clean data were compared with reference genomes using hisat2 (v2.2.1; GitHub Release) software. Finally, the sample uniformity analysis was carried out. As expected, the sequencing data from all samples were of high quality.
Gene expression analysis
Stringtie (v2.1.5; Johns Hopkins University) software was used to predict transcripts of all samples, and then RSEM (v1.3.1; GitHub Release) software was used to call bowtie2 comparison results for statistics. The number of reads of each sample to each transcript was obtained. Then, they were converted using fragments per kilobase per million bases. The paired-end reads from the same fragment were counted as one fragment, and the expression levels of genes were obtained.
Differential expression analysis
After a comprehensive analysis of gene expression levels, the sample correlation and differential expression analyses were performed. The differential gene analysis (differentially expressed gene [DEG]) between sample groups was conducted using the R language package DEseq2 (v1.22.2; European Molecular Biology Laboratory), with screening thresholds of (false discovery rate) p<0.05, log2FC >1, or log2FC <-1. Next, the screened DEGs were analyzed by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses. GO and KEGG analyses were performed using hypergeometric distribution, with screening thresholds of p≤0.05.
GO analysis
A GO analysis of biological processes, molecular functions, and cellular components is the international standard classification system of gene function. A study on the distribution of DEGs in the GO can clarify the differences in gene function.
KEGG analysis
The most critical biochemical metabolic pathways and signal transduction pathways involved in DEGs can be identified through pathway enrichment analysis. KEGG is the main public database on the pathways. Pathway enrichment analysis was performed using KEGG as a unit. The hypergeometric test was used to find out the pathways that showed significant enrichment in DEGs compared with the whole-genome background. R software (R Foundation for Statistical Computing) and a self-written script were used to set parameters (BH correction was used) for pathway enrichment analysis.
The ASD validation set
In the validation cohort, RNA was reverse transcribed to cDNA using high-capacity cDNA reverse transcription kits (Applied Biosystem Inc.), using the following protocol: 10 minutes at 25°C, 2 hours at 37°C, and 5 minutes each at 85°C and 4°C. The quantitative polymerase chain reaction (qPCR) was performed with SYBR Green PCR Master Mix (Applied Biosytems Inc.) on a Light Cycler 96 system (Roche Applied Science). The thermal cycling conditions were as follows: 95°C for 10 minutes and 40 cycles at 95°C for 15 seconds and 60°C for 1 minute. The corresponding primers are shown in Table 2. All transcript levels were normalized to GAPDH expression. Fold-change expression was calculated by the 2-ΔΔCT method.
Primer pairs for quantitative polymerase chain reaction
In biometric verification, the Series Matrix File data file of GSE18123 was downloaded from the NCBI GEO public database. The GSE18123 dataset (platform GPL6244) included peripheral blood samples of 82 control subjects and 104 ASD children. The “limma” package of the R software was applied for differential gene analysis.
Statistical analysis
SPSS 25 software (IBM Corp.) was used for the statistical analysis. All qPCR results were presented as mean±standard error of the mean. The group means were compared using Student’s t-test as indicated. A p-value <0.05 (two-tailed) indicated a statistically significant difference for all tests.
RESULTS
DEG analysis
The DEG analysis results in the two groups were presented in the follow. The volcano plots showed the differential DEGs in Figure 1A, which were 21 up-regulated and 35 down-regulated. Subsequently, a Z-normalization was applied to all behavioral parameters, with the objective of combining data from all tests into a single index (Figure 1B). A correlation was conducted between all of the differential genes and the Z-scores. The results only demonstrated a negative correlation between the up-regulated genes, hydroxycarboxylic acid receptor 3 (HCAR3), and the down-regulated genes, tumor necrosis factor ligand superfamily member 15 (TNFSF15), and the Z-scores (p<0.05) (Figure 1C).
The analysis between DEGs and behavioral severity Z-score index. A: Volcano plot of DEGs. B: Heat map for behavioral severity Z-score index. C: Spearman analysis between DEGs and behavioral severity Z-score index. DEG, differentially expressed gene; ABC, autism behavior checklist; CARS, Childhood Autism Rating Scale; DDSC, Developmental Diagnostic Scale of Children Aged 0–6 Years; HCAR3, hydroxycarboxylic acid receptor 3; TNFSF15, tumor necrosis factor ligand superfamily member 15.
GO analysis
The results of the GO and KEGG enrichment analysis based on DEGs were shown in Figure 2. As illustrated in Figure 2A and B, the high-enrichment GO terms included the following: cellular process, biological regulation, metabolic process, signaling, response to stimulus, cellular anatomical entity and binding. The principal enrichment pathways identified in the KEGG enrichment analysis were the signal transduction, signal molecules and interaction, transport and catabolism, immune disease, and development and regeneration pathways (Figure 2C and D).
GO and KEGG pathway analysis of common DEGs. Bar charts of downregulated (A) and upregulated (B) DEG GO classification in the cohort. Bar charts of downregulated (C) and upregulated (D) KEGG pathway classification analysis in the cohort. The X-axis refers to the number of genes, and the Y-axis refers to GO pathway name or the KEGG metabolic pathway. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; DEG, differentially expressed gene.
The ASD validation set
We used qPCR to verify the two DEGs derived from RNAseq in the validation group (children with ASD vs. sex- and age-matched TD children, n=10 each). As illustrated in Figure 3, the validation cohort exhibited a notable alteration with a similar trajectory as observed in the RNA-seq data for HCAR3 (p=0.014), whereas TNFSF15 demonstrated no change (p=0.822). We employed bioinformatics techniques to examine transcriptional data from human specimens with a view to elucidating its potential involvement further. The analysis was conducted using the GSE18123 dataset, a high-throughput sequencing dataset generated from the peripheral blood of patients with ASD. After preprocessing, the gene expression level of HCAR3 in ASD patients was observed to be higher than that in controls, with a similar result obtained by qPCR (p<0.01) (Figure 3B).
The expression of HCAR3 and TNFSF15. A: The mRNA expression in the validation group. B: Gene expression levels in GSE18123 dataset. *p<0.05; **p<0.01. HCAR3, hydroxycarboxylic acid receptor 3; TNFSF15, tumor necrosis factor ligand superfamily member 15; ASD, autism spectrum disorder; TD, typical development; ns, no significance.
DISCUSSION
In this study, we detected mRNA expression profiles in PBMC and identified DEGs by RNA-seq. A total of 56 significantly differentially expressed mRNAs were found, including 21 which were upregulated and 35 which were downregulated mRNAs. GO and KEGG functional enrichment analyses show that these 56 DEGs were mainly localized in cellular processes, signaling, immune and so on. Moreover, we found that HCAR3 and TNFSF15 were related to ASD symptoms. Finally, our findings suggested that HCAR3 may be a risk gene for ASD.
The findings of our study indicated that HCAR3, a protein involved in signaling pathways and one of the main receptors for nicotinic acid, was elevated in ASD children. Previous research pointed out that many signaling proteins were also demonstrated to be involved in the pathogenesis of ASD [14,15]. Meanwhile, we found that signaling pathways were strikingly associated with ASD by KEGG and GO analyses. The signaling of molecular events involved in a series of pathological processes, such as synaptic signaling transmission, has been identified as a potential target for interventions in ASD [16]. Through RNA sequencing of hippocampal tissue from ASD patients, DEGs were found to be enriched in “synaptic regulation,” and molecular biological experiments also revealed reduced expression of synaptic markers [17]. Moreover, some studies demonstrated that there was a strong correlation between the dysregulation of synaptic signaling transmission and the manifestation of social deficits and memory loss in ASD rats [18]. What’s more, Shank3 (a synaptic-related gene) mutant mice have been widely utilized as an ASD model in exploring the etiological mechanisms of ASD [19,20]. In addition, HCAR3 was located in the 12q24.31 chromosomal region, where patients with 12q24.31 mutations were diagnosed with ASD [21,22]. Elevated expression of the HCAR3 gene influenced the niacin response in patients with psychiatric disorders, and HCAR3 may serve as a risk gene for these conditions [23,24]. Furthermore, mice deficient in HCAR1, which is the homologous protein of HCAR3, exhibited a reduction in social behavior accompanied by an increase in repetitive and anxiety-like behaviors [25]. Notably, the expression of HCAR3 was regulated by nicotinic acid and a deficiency in the range of 50%–60% for nicotinic acid was detected in ASD children by a cross-sectional study comprising 285 ASD children and 224 relatives of these children [26,27]. The mRNA for HCAR3, also known as HM74, was found to be increased in the anterior cingulate cortex of subjects with other neurological diseases when compared to controls [28,29]. All of the above evidence indicates that HCAR3 plays an important role in maintaining the integrity of the biological functions associated with ASD.
In addition, enrichment results showed that the immune pathway was abnormally enriched in this study. There is evidence to suggest that dysfunction of the immune system may be associated with ASD. Several pro-inflammatory molecules were elevated in the serum and cerebrospinal fluid of patients diagnosed with ASD [30,31]. Similarly, some studies have indicated that there is evidence of immunological dysregulation in animals with ASD, with an elevation in pro-inflammatory factors [32,33]. Notably, the maternal immune activation mouse has been employed extensively as a model for ASD in mice [34,35]. Furthermore, immunomodulatory drugs could markedly enhance the frequency and duration of social sniffing in ASD mice [36]. In the present study, we also found that immune-related protein (TNFSF15) was associated with behavioral symptoms in ASD children. TNFSF15, also called TL1A, is a member of the tumor necrosis factor family, expressed in various immune cells [37]. TNFSF15-induced expression of TNF-α served to regulate the intestinal microenvironment, thereby exacerbating intestinal inflammation [38]. Some studies reported that the level of TNF-α was elevated in ASD patients and the regulation of TNF-α enabled the amelioration of ASD-like phenotypes [39,40]. In other neurological disorders, suppressing the expression of TNFSF15 could reduce inflammatory responses, decrease neuronal apoptosis, and thereby ameliorate ASD-related behavioral phenotypes [41,42]. Nevertheless, no significant difference was observed in the validation group with respect to TNFSF15. Further experiments are needed to determine whether TNFSF affects the development of ASD.
Although our study has some important results, it still has some limitations. First, it would be beneficial to validate our findings with a larger sample size, to control for intraindividual variations. Second, the cohort in this study is exclusively male, which may limit the generalizability of the findings to females. Third, the methodology presented in this study is solely based on RNA-seq. Future studies should incorporate additional omics analyses, such as those of the proteome and metabolome, to gain a more comprehensive understanding [43,44].
In conclusion, we obtained 56 DEGs using RNA-seq data from PBMC of ASD children and TD children. The GO and KEGG pathway analyses indicated that the target genes identified in this study may be involved in signaling and immune pathways. Following both extrinsic and intrinsic validation, we confirmed that HCAR3 may be a risk gene associated with ASD. The current findings may offer preliminary evidence for the potential biomarkers and provide novel insights into the genetic etiology of ASD.
Notes
Availability of Data and Material
The data analyzed during the current study are available from the corresponding author on reasonable request.
Conflicts of Interest
The authors have no potential conflicts of interest to disclose.
Author Contributions
Conceptualization: Mingyang Zou, Caihong Sun. Formal analysis: Sen Yang, Jingting Xu. Funding acquisition: Mingyang Zou, Caihong Sun. Methodology: Feng Wang. Supervision: Caihong Sun. Validation: Feng Wang, Jiao Zuo. Writing—original draft: Feng Wang. Writing—review & editing: Mingyang Zou.
Funding Statement
This work was supported by the University Nursing Program for Young Scholars with Creative Talents in Heilongjiang Province (Grant No.UNPYSCT-2020168).
Acknowledgments
None