Altered expression of histamine signaling genes in autism spectrum disorder.

The histaminergic system (HS) has a critical role in cognition, sleep and other behaviors. Although not well studied in autism spectrum disorder (ASD), the HS is implicated in many neurological disorders, some of which share comorbidity with ASD, including Tourette syndrome (TS). Preliminary studies suggest that antagonism of histamine receptors 1-3 reduces symptoms and specific behaviors in ASD patients and relevant animal models. In addition, the HS mediates neuroinflammation, which may be heightened in ASD. Together, this suggests that the HS may also be altered in ASD. Using RNA sequencing (RNA-seq), we investigated genome-wide expression, as well as a focused gene set analysis of key HS genes (HDC, HNMT, HRH1, HRH2, HRH3 and HRH4) in postmortem dorsolateral prefrontal cortex (DLPFC) initially in 13 subjects with ASD and 39 matched controls. At the genome level, eight transcripts were differentially expressed (false discovery rate <0.05), six of which were small nucleolar RNAs (snoRNAs). There was no significant diagnosis effect on any of the individual HS genes but expression of the gene set of HNMT, HRH1, HRH2 and HRH3 was significantly altered. Curated HS gene sets were also significantly differentially expressed. Differential expression analysis of these gene sets in an independent RNA-seq ASD data set from DLPFC of 47 additional subjects confirmed these findings. Understanding the physiological relevance of an altered HS may suggest new therapeutic options for the treatment of ASD.

Introduction

Autism spectrum disorder (ASD) is a class of neurodevelopmental disorders (NDDs) characterized by social impairment, repetitive behavior and impaired communication.[1] Genetic studies indicate that the heritability of ASD is approximately 50–60%.[2, 3] The underlying genetic architecture is complex, with risk conferred by many genetic and environmental factors. Many of these factors appear to be shared by other NDDs,[4] consistent with a high level of comorbidity between ASD and other NDDs.[5] Tourette syndrome (TS) shares genetic risk factors with ASD[6, 7] and is among the most prevalently comorbid NDD with ASD.[5] This may be, in part, due to shared upregulation in neuroinflammation,[8, 9, 10] increased microglia activation[11] and/or abnormalities in language processing related to frontal/basal ganglia circuit dysfunction.[12, 13, 14] Indeed, boys with ASD show faster response times in procedural/grammar memory tasks similar to children with Tourette syndrome, which has been attributed to a reduction in the connectivity between the temporal lobe and the frontal/basal ganglia circuitry.[12] Interestingly, evidence suggests that overactivation of microglia may result in the reduced connectivity seen in ASD.[15]

Each of these shared abnormalities may be mediated in part by the histaminergic system (HS). Histamine has a prominent role in nueroinflammation[10, 16] and microglial regulation,[11, 17, 18, 19] and recent studies suggest a role for the HS in TS, possibly reflecting alterations within basal ganglia circuitry.[20, 21, 22] Critically, a single family study of TS demonstrated a causative role of a rare nonsense mutation in the histidine decarboxylase gene (HDC), which encodes the rate-limiting enzyme that converts histidine into histamine.[20] Interestingly one of the TS subjects in this study also had a diagnosis of Asperger’s syndrome. Subsequent reports showed that two SNPs in the HDC region also associated with TS in a family genetic association study.[23] Furthermore, recent studies of ASD also implicate potential involvement of HS in this disorder, as de novo deletions overlap histamine N-methyltransferase (HNMT), the gene encoding the enzyme that inactivates histamine.[24, 25]

The role of the HS in neuroinflammation, cognition, sleep, attention, sensory function and motor function has received considerable attention[26, 27, 28] and recent evidence suggests that the HS may have a role in regulating microglia activation,[11, 17] cytokine release[18] and migration.[19] Furthermore, the HS may influence microglia mediated phagocytosis and reactive oxygen species production through the activation of histamine receptors.[29] Within the central nervous system (CNS), this system consists of histaminergic neurons within the tuberomammillary nucleus (TMN) of the hypothalamus that have diffuse projections into the spinal cord, brain stem and many telencephalic brain regions. Histamine synthesis occurs predominantly in the TMN, where HDC is largely expressed. In addition, mast cells and microglia also express HDC and contribute to brain histamine, particularly in response to inflammation.[30, 31] Histamine inactivation by HNMT occurs throughout the brain. Four histamine G-protein-coupled receptors (H1R–H4R), encoded by the genes Histamine Receptor H1 (HRH1), Histamine Receptor H2 (HRH2), Histamine Receptor H3 (HRH3) and Histamine Receptor H4 (HRH4) are expressed throughout the CNS, although H4R expression is extremely low and its function within the CNS is not well characterized. Microglia express all four receptors and evidence suggests that histamine dose-dependently activates microglial activation and production of proinflammatory factors by binding H4R[11] and possibly H1R.[17]

Antagonists/inverse agonists of H3R have been investigated as a therapy for a number of CNS diseases, including Alzheimer’s disease, epilepsy, narcolepsy, attention deficit hyperactivity disorder and fetal alcohol exposure-induced learning deficits.[26, 32, 33, 34, 35] Unfortunately, clinical trial results have been largely disappointing; however, preclinical studies demonstrate promising results and early clinical trials of pitolisant showed encouraging results in the treatment of narcolepsy.[34] Pitolisant may also prove efficacious in epilepsy, a condition highly comorbid with ASD.[33] In addition, clinical studies demonstrate behavioral and sleep disturbance improvement in children and adolescents with ASD treated with H2R and H1R antagonists.[36, 37]

These therapeutic findings, combined with the recent genetic findings implicating HNMT in ASD and HDC in TS, the comorbidity of ASD and TS, and the overlap of a potential role for neuroinflammation, altered microglial activity and altered language processing between the two disorders raise the possibility that the HS is involved in ASD. We hypothesized, therefore, that the expression of histaminergic genes is altered in the brain in ASD. We focused on the expression of the following histaminergic genes: HDC, HNMT, HRH1, HRH2, HRH3 and HRH4 in the brains of subjects diagnosed with an ASD. Given the possibility of shared alterations in connectivity between the basal ganglia and frontal lobe in both TS and ASD, as well as the associated alterations in frontal lobe processing in general in ASD,[29] we selected tissue samples from the dorsolateral prefrontal cortex (DLPFC). Based on the positive treatment findings in patients and animal models of ASD using H1R,[37] H2R,[36] and H3R[35] antagonists, we predicted that the genes encoding these receptors would be more highly expressed in ASD subjects. We predicted that expression of HNMT would be altered, but we were unsure of the direction of change, given the overlap of HNMT with associated ASD deletions, and findings suggesting that polymorphisms of increased activity of HNMT are associated with modulating attention and hyperactivity behavior.[38, 39]

Materials and methods

A flow chart of the methods are outlined in Supplementary Figure 1.

Subjects and demographics

We evaluated the expression of histaminergic genes in the DLPFC of 13 individuals diagnosed with ASD and 39 non-psychiatric control subjects using deep RNA-seq coverage. Three control subjects were matched to each case subject by age (±6 years), gender and ethnicity when possible. Only five control subjects mismatched for ethnicity. A summary of subject demographic information is presented in Table 1 and more extensively in Supplementary Table 1.

Table 1

Subject demographics

Diagnosis	Gender (F/M)	Mean RIN (range)	Mean age (range)	Race (AA/CAUC)
ASD	3/10	8.57 (6.6 to 9.4)	22 (4 to 67)	6/7
Control	9/30	8.09 (6.0 to 10.0)	22 (2 to 69)	14/25

Abbreviations: AA, African American; ASD, autism spectrum disorder; CAUC, Caucasian; F, female; M, male; RIN, RNA integrity number.

This table indicates the demographic information for the subjects used in this analysis. Subjects with evidence of drug use, alcohol abuse or psychiatric illness were excluded from the control cohort.

The majority of samples were collected at the National Institute of Mental Health in the Section on Neuropathology, according to NIH protocol 90-M-0142, with informed consent of next-of-kin. Additional samples were consented through the University of Maryland Brain and Tissue Bank (formerly the NICHD Brain Bank), and one sample was consented through the Stanley Medical Research Institute. All the materials were transferred to the Lieber Institute for Brain Development under approved material transfer agreements. For details on National Institute of Mental Health brain collection, clinical characterization, neuropathological examinations and toxicological analyses, see Lipska et al.[40] For details on the University of Maryland Brain and Tissue Bank and the Stanley Medical Research Institute sample characterization, see their respective websites at: www.medschool.umaryland.edu/btbank/ and www.stanleyresearch.org. Subjects with evidence of drug use, alcohol abuse or psychiatric illness were excluded from the control cohort. All the analyzed samples had RNA integrity number (RIN) values >6.

RNA extraction and sequencing

All RNA extraction and sequencing was performed at the Lieber Institute for Brain Development. RNA was extracted from postmortem homogenate brain tissue of gray matter dissected from Brodmann area 46 and 9. Total RNA was extracted using the RNeasy kit (Qiagen, Germantown, MD, USA) from approximately 100 mg of homogenate tissue. RNA-seq was performed using the TruSeq Stranded Total RNA Library Preparation kit with Ribo-Zero Gold ribosomal RNA depletion from Illumina (San Diego, CA, USA). One hundred base pair paired-end sequencing was run on the HiSeq 2000. The Illumina Real Time Analysis (RTA) module was used to perform image analysis and base calling and the BCL Converter (CASAVA v1.8.2) was used to generate the sequence reads. Sequencing depth was 40–60 million paired-end reads, yielding 80–120 million reads per sample. TopHat[41] (v2.0.4) was used to enforce strand specificity and to align the sequencing reads to known transcripts of the Ensembl Build GRCh37.67, (a total of 57 659 genes with 196 495 known transcripts). Visual inspection of reads at each gene of interest was evaluated for each subject and no structural variants were found.

Differential expression analysis

We used featureCounts[42] (v1.4.4) to count the total number of reads overlapping each gene using the default settings, with paired-end and reverse-stranded counting specified using the Ensembl Build GRCh37.67 gtf file. These counts were merged from both reads of the paired-end sequencing and normalized by library size and coding gene length to form Reads Per Kilobase of Gene per Million mapped reads (RPKM) values. These RPKM values were transformed using log2 after applying an offset of 1 to each count to stabilize the variance among lowly expressed genes and to prevent negative normalized counts, resulting in normalized gene expression levels: log2(RPKM+1). These normalized expression values were then filtered for genes with mean normalized expression values less than 0.5, reducing the total number of genes to 13 011. The influence of diagnosis status on all normalized and filtered gene abundance estimates was evaluated using a linear regression analysis covarying for principal components determined by the transcriptome-wide RNA-seq expression values for the 13 011 expressed Ensembl genes to control for potential technical artifacts and latent biological effects.

The appropriate number of principal components (PCs) to include in the model to correct for known and unknown confounds, yet still maintain variance related to diagnosis, was determined by the num.sv() function[43] of the Bioconductor package ‘sva’ in R.[44] This method uses parallel analysis to calculate the number of components that significantly contribute to overall variance, based on how often the eigenvalue for each PC is larger than eigenvalues calculated for random permutations of the data. This process is performed after removing variation due to specified variables of interest, in our case: diagnosis, RIN, sex, age, race and exonic mapping rate. This function determined that nine PCs significantly contributed to the overall variance. The percent variance contributed by each PC was as follows: 58.5, 9.2, 7, 3.8, 2.7, 2.3, 1.8, 1.5 and 1.2%, accounting for 88% of the overall variance. See Supplementary Figure 2 for a plot of the first two PCs and Supplementary Figure 3 for a heatmap indicating the association of each component with various known confounders.

After the model optimization, linear regression analysis was performed on the normalized gene expression levels across the transcriptome, modeling the effect of diagnosis while adjusting for the first nine PCs, as well as for age, RIN, sex, race and exonic mapping rate. Resulting t-statistics and corresponding P-values were calculated using the empirical Bayes method in the Bioconductor package ‘limma’ in R.[45, 46] Multiple testing correction was performed using the Benjamini–Hochberg false discovery rate method. The histamine genes (HRH1, HRH2, HRH3, HRH4, HDC and HNMT) were then evaluated in the context of this genome-wide differential expression analysis. We also considered only male subjects (10 ASD and 30 controls) as sensitivity analyses to determine how much heterogeneity across the genders of the cases and controls mattered in the differential expression analysis, as recent studies indicate the potential existence of gender-specific protective and risk genetic factors.[47] In the male subjects, we found that adjusting for the first seven PCs was the optimal number to include in the differential expression model to control for latent heterogeneity.

Gene set enrichment analysis

Given our interest in the impact of the HS as a pathway in ASD, we also evaluated these genes using gene set enrichment tests (GSETs). In addition, we considered five other curated histamine-related gene sets not based on our a priori hypotheses. HDC and HRH4 were more lowly expressed than the expression filter (<0.5) and were therefore not included in this analysis. We analyzed the expression of our target gene set including HRH1, HRH2, HRH3 and HNMT using the Bioconductor limma ROAST test, a rotation gene set test that is a self-contained gene set test and therefore tests whether the genes within the gene set are differentially expressed.[48] The advantage of this test is that it does not rely on the assumption of gene independence[48] and it has been found to be robust to small sample sizes and sample heterogeneity.[49] Furthermore, this test is robust to small gene set sizes, as it does not rely on the hypergometric test, which tests for enrichment of the testing pre-specified gene set against all genes identified to have statistics below a significance threshold. Many gene set tests implement the hypergometric test, which can lead to inflated significance of small gene sets when one gene is highly differentially expressed. However, the ROAST test instead combines all single gene tests statistics into one set-level statistic, which is compared with calculated rotated gene set statistics. This test was performed using the normalized expression values and design matrix as used in the previous analyses. This test estimates significance based on random rotations of the orthogonalized residuals and provides output for three directional hypotheses: up, down and mixed. A mixed hypothesis indicates that the gene set includes extreme t-statistics and thus that genes within the gene set represent a large increase or decrease in expression.[48] The other gene sets evaluated included three from gene ontology according to http://amigo.geneontology.org/: histamine receptor activity (GO:0004969), histamine secretion (GO:0001821) and histamine production involved in inflammatory response (GO:0002349), as well as two from GeneRIF identified using Harminizome:[50] histamine and histaminergic. Genes within these gene sets that were expressed above the minimum expression threshold were evaluated in the GSETs. For a list of these genes, see Supplementary Table 2. We hypothesized that all gene set results would either be in the up or mixed direction, based on the positive treatment findings in patients and animal models of ASD using H1R,[37] H2R[36] and H3R[35] antagonists and the inclusion of these genes or related genes in each gene set. In this analysis, 10 000 random rotations were run to obtain more precise statistical results. Significance values for the mixed and up hypothesis tests were corrected for multiple testing correction using the Benjamini–Hochberg false discovery rate method. The effect sizes of the q-values were obtained using the pes() function of the compute.es package[51] in R.

Replication analysis

We performed a replication analysis using data recently published with independent samples in another recent ASD RNA-seq study.[52] We received raw data from the Geschwind lab at UCLA for Ribo-Zero RNA-seq samples from 57 subjects including controls and cases with ASD and RIN >6 from the DLPFC, (Brodmann area 9). We excluded six ASD samples that were identified to have a 15q11.2-13.1 duplication, as the manifestation of ASD in these subjects may be due to a different pathogenesis. In addition, we excluded four samples because they overlapped with our original data set. Therefore, we used the raw data for a total of 47 subjects (28 controls and 19 patients). See Supplementary Table 3 for more extensive demographic data for these subjects.

Using these subjects allowed for a strict replication analysis, as the differential expression and GSETs were performed identically as in our data. The only difference was that these data were not strand-specific, and thus the setting used for counting was not strand-specific. Again, the num.sv() function[43] of the Bioconductor package ‘sva’ in R was used to determine the optimal number of PCs to include in the model. This process was performed after removing variation due to specified variables of interest, in our case: diagnosis, RIN, sex, age and the percentage of assigned reads to Ensembl genes. This function determined that six PCs significantly contributed to the overall variance. The percent variance contributed by each PC was as follows: 32.9, 14.8, 12.8, 6.3, 5.2 and 3.8%, accounting for 75.8% of the overall variance. See Supplementary Figure 4 for a plot of the first two PCs and Supplementary Figure 5 for a heatmap indicating the association of each component with various known confounders. Only genes or gene sets with significant findings in the original data set were evaluated in the replication data. TRPC3 was not expressed in the replication data set and thus not tested when testing the GeneRIF histamine gene set.

Meta-analysis

Meta-analysis was performed using the Fisher’s method and the P-value statistics from both the original and the replication data sets using the sumlog() function of the ‘matap’ package[53] in R. The histamine genes of interest, the genes that showed differential expression in the discovery analysis and the gene sets of interest that showed positive results in the discovery analysis were tested. Multiple testing correction was done using the Benjamini–Hochberg method. Importantly, this evaluation provides a measure for the statistical significance of the absolute value of gene expression changes across both the data sets but does not account for the directionality of the changes in expression.

Code availability

The code used to generate the results of this study is accessible at https://github.com/LieberInstitute/Wright_2017_ASD_Histamine. The information about data access can also be found at this site.

Results

Differential expression analyses reveal no difference in the expression of individual histamine genes

We identified 1463 genes differentially expressed between ASD subjects and healthy controls across all moderately expressed Ensembl genes (13 011) at marginal statistical (P<0.05) significance. Three control subjects were matched to each ASD subject for age and sex and nearly all subjects were also matched for ethnicity when possible. See Supplementary Table 4 for a full list of the test statistics for all evaluated genes. Following multiple testing correction, eight remained significantly differentially expressed at false discovery rate <5%: SNORA74A, SNORA53, SNORD17, TUBE1, SNORA54, SNORA74B, a lncRNA HSA-LNCG003387, also called RP6-206I17.3 or BX284650.2, and SNORD114-23, (Supplementary Table 5). None of the histamine genes of interest were significantly differentially expressed between cases and controls (Table 2); however, HNMT was nominally significantly differentially expressed, showing increased expression in ASD (P=0.003). HDC and HRH4 were excluded from the analysis as their mean estimated expression value was below the lower limit threshold (<0.5), as expected from prior knowledge of their expression. The analysis of only male subjects yielded 1027 differentially expressed genes before multiple testing correction, and three survived multiple testing correction including: SNORA54, SNORA74A and SNORA53 (false discovery rate <5%, Supplementary Table 6). The histaminergic gene results are presented in Supplementary Table 7.

Table 2

Statistical summary for the histaminergic genes of interest in the transcriptome-wide differential expression analysis

Gene symbol	Log2 fold change	P-value	q-value	Control mean expression	ASD mean expression
HNMT	0.080	0.003	0.245	0.74	0.82
HRH1	0.095	0.072	0.484	0.12	0.22
HRH2	0.029	0.445	0.819	1.01	1.04
HRH3	0.061	0.166	0.626	0.94	1.00

Statistical summary of the histamineric genes of interest that were sufficiently expressed for differential expression analysis between autism spectrum disorder (ASD) and control subjects in the dorsolateral prefrontal cortex. Log fold change is indicated with the controls as the reference group; therefore, each of these genes is more highly expressed in ASD subjects than in control subjects. The P-value indicates the nominal significance for each gene in the differential expression analysis. The q-value indicates the multiple testing corrected significance value for each gene. Mean expression is expressed as log2(RPKM+1) values and reflect correction of confounding effects by principal components.

We performed differential expression analysis using the same methods with the replication sample of 47 subjects. This yielded 1090 significantly differentially expressed genes before multiple testing correction (out of 13 473 that passed the minimum expression filter), and zero after multiple testing correction. Again, none of the histamine genes of interest were significantly differentially expressed individually following multiple testing correction, (Supplementary Table 8); however, HRH3 (P=0.013) was nominally significantly differentially expressed. HRH1 and HRH2 were found to be more lowly expressed in the patient group of the replication sample, and therefore showed the opposite direction of change as that of the original sample. None of the eight genes identified as significantly differentially expressed in our data were significantly differentially expressed in the replication data set at even nominal significance and roughly only half of these genes showed the same direction of fold change, as would be expected by chance, (Supplementary Table 9). The lack of consistency between the two data sets may be due to heterogeneity among samples and the relatively small size of our samples. Supplementary Table 10 shows the meta-analysis results for each of the histamine genes of interest and the genes identified to be significantly differentially expressed in the original data set. All differentially expressed genes identified in the original data set were also significant in the meta-analysis. HNMT (q=0.012), HRH1(q=0.033) and HRH3 (q=0.019) were also significant. However, HRH1 showed the opposite change of direction between the two data sets. HRH2 (q=0.296) was not significant in the meta-analysis.

Gene set enrichment discovery and replication analysis reveals that collectively the histamine genes show altered expression in the ASD subjects

In contrast to the individual gene differential expression analyses, our target GSET of HNMT, HRH1, HRH2 and HRH3 using the ROAST test revealed that this gene set is significantly overexpressed (q-value =0.01), and significantly differentially expressed regardless of directionality (‘mixed’ test, q=0.02) in ASD patients, (Table 3). Evaluation of the fold changes within just the male subjects and our full data set in a sensitivity analysis indicates that the full data set was not biased by potential gender differences between diagnosis groups, (Supplementary Figure 6). Given the strength of the effect sizes, our results provide compelling evidence that the expression of the evaluated HS-related genes may be collectively altered in the DLPFC of ASD subjects. Furthermore, the replication analysis using the sample of 47 subjects confirmed that this gene set is similarly altered in expression in ASD prefrontal cortex. The gene set was found to be significantly different regardless of directionality also in this independent sample of ASD subjects, with the ROAST ‘mixed’ test (q-value =0.03, Table 4).

Table 3

Statistical summary for gene set enrichment analysis of the histaminergic genes of interest

Gene set	Up P-value	Up q-value	Effect size g	Mixed P-value	Mixed q-value	Effect size g	Number of genes	Up active prop	Mixed active prop
Hypothesis-driven gene set	0.001	0.010	0.84 (0.19, 1.5)	0.006	0.024	0.73 (0.08, 1.38)	4	0.50	0.50
GO: histamine receptor activity	0.015	0.036	0.68 (0.03, 1.33)	0.096	0.193	0.42 (−0.22, 1.06)	3	0.33	0.33
GeneRIF: histaminergic	0.461	0.564	0.18 (−0.45, 0.82)	0.011	0.032	0.70 (0.05, 1.35)	5	0.20	0.60
GeneRIF: histamine	0.373	0.564	0.18 (−0.45, 0.82)	0.003	0.016	0.79 (0.13, 1.44)	41	0.15	0.34
GO: histamine secretion	0.452	0.564	0.18 (−0.45, 0.82)	0.976	0.976	0.01 (−0.62, 0.64)	6	0	0
GO: histamine production in inflammation	0.470	0.564	0.18 (−0.45, 0.82)	0.911	0.976	0.01 (−0.62, 0.64)	4	0	0

Abbreviation: GO, Gene Ontology.

Summary of results for the gene set enrichment tests of the histaminergic genes. The results are shown with the controls as the reference group. Therefore, the hypothesis of ‘Up’ refers to the gene set being increased in autism spectrum disorder (ASD) patients as compared with controls, whereas the hypothesis of ‘Mixed’ refers to the gene set including genes with extreme t-statistics in both directions. Thus, as the results indicate that several of the gene sets show altered expression in ASD, the active prop is the proportion of genes in the set contributing meaningfully to significance, defined as those with squared z-values greater than 2.

Table 4

Statistical summary for gene set enrichment analysis of the histaminergic genes of interest in the replication data set

Gene set	Up P-value	Up q-value	Effect size g	Mixed P-value	Mixed q-value	Effect size g	Number of genes	Up active prop	Mixed active prop
Hypothesis-driven gene set	0.441	0.588	0.16 (−0.43, 0.75)	0.009	0.026	0.67 (0.07, 1.28)	4	0.25	0.50
GO: histamine receptor activity	0.656	0.749	0.09 (−0.49, 0.68)	0.009	0.026	0.67 (0.07, 1.28)	3	0.33	0.67
GeneRIF: histaminergic	0.904	0.904	0.04 (−0.55, 0.62)	0.008	0.026	0.67 (0.07, 1.28)	5	0.20	0.60
GeneRIF: histamine	0.080	0.150	0.43 (−0.17, 1.02)	0.170	0.269	0.33 (−0.27, 0.92)	40	0.20	0.25

Abbreviation: GO, Gene Ontology.

Summary of results for the gene set enrichment tests of the histaminergic genes in the replication data set. The results are shown with the controls as the reference group. Therefore, the hypothesis of ‘Up’ refers to the gene set being increased in autism spectrum disorder (ASD) patients as compared with controls, whereas the hypothesis of ‘Mixed’ refers to the gene set including genes with extreme t-statistics in both directions. Thus, as the results indicate that several of the gene sets show altered expression in ASD, the active prop is the proportion of genes in the set contributing meaningfully to significance, defined as those with squared z-values greater than 2.

The Gene Ontology (GO) receptor gene set (GO:0004969) also was significantly upregulated (q=0.04), and both the GeneRIF gene sets showed significantly mixed expression (q=0.03 and q=0.02; Table 3) in the Lieber Institute for Brain Development sample. However, the GO histamine secretion (GO:0001821) and histamine production in inflammation (GO:0002349) gene sets did not show differential expression between the two subject groups. The GO receptor gene set and the GeneRIF histaminergic gene set also showed significantly mixed expression, (q=0.03 and q=0.03, Table 4) in the independent replication sample. See Figure 1 for a graphical representation of our target GSET in both the data sets. See Supplementary Figure 8 for a plot of the difference of the absolute log2 fold changes of these genes between the two data sets. Supplementary Table 10 shows the results of meta-analysis for the gene sets statistics across the two data sets. Our hypothesis-driven gene set using the ROAST up (q=0.004) and mixed (q=0.001) statistics were significant in the meta-analysis. In addition, the curated gene sets (evaluated in both data sets) using the ROAST mixed statistics were also significant.

Figure 1

Magnitude of altered expression of histamine-related genes in autism spectrum disorder (ASD). (a) The absolute log2 fold changes of gene expression between ASD subjects and control subjects within the original data set are shown on the x axis for all Ensembl genes evaluated in this differential expression analysis, while the y axis shows the frequency, or number of genes with each log2 fold change value. The x axis shown is limited to changes between 0 to 0.2 to allow easier evaluation of the genes of interest, however, the full range of absolute log2 fold changes was 4.9 × 10−6 to 1.2. This indicates the extent to which genes showed differential expression, regardless of the direction of change. The red lines indicate the absolute log2 fold changes of the histamine-related genes of interest. (b) The absolute log2 fold changes of gene expression between the ASD subjects and the control subjects within the replication data set are shown on the x axis for all Ensembl genes evaluated in this differential expression analysis, while the y axis shows the frequency, or number of genes with each log2 fold change value. The x axis shown is limited to changes between 0 to 0.2 to allow easier evaluation of the genes of interest, however, the full range of absolute log2 fold changes was 8.7 × 10−7 to 6.02. Again, the red lines indicate the absolute log2 fold changes identified for the genes of interest.

Discussion

Given recent HS findings in TS, the shared comorbidity between TS and ASD, the role of the HS in neuroinflammation and preliminary clinical trial data of histamine receptor inverse agonists/antagonists in ASD, we hypothesized that gene expression of histaminergic genes may be altered in postmortem brains of ASD patients. We evaluated the expression of individual HS genes as well as gene sets comprised of the genes encoding the more highly expressed histamine receptors in the CNS and the gene encoding the enzyme that inactivates and controls the activity of histamine.

None of our genes of interest were individually identified to be differentially expressed in the DLPFC of our sample or in our replication data set. Given the sample sizes involved, this is not inconsistent with what might have been expected. However, collectively as a gene set, our results suggest that the expression of these genes as representing a histaminergic signaling system, do show differential expression in the DLPFC of subjects with ASD compared with controls. Although this gene set had significant up and mixed ROAST results in the discovery sample, only the mixed result (of overall differential expression regardless of directionality) was replicated in the replication sample, an independent cohort of 47 subjects with RNA-seq data derived from Brodmann area 9. The directionality of a subset of specific genes varied across the discovery and replication data sets, yet the significance of overall differential expression suggests that the HS may be perturbed in ASD at a system level. Given that these genes regulate the response to histamine and therefore the influence of this neurotransmitter system on other neurotransmitters, such changes could be consequential for brain function and behavior. In addition, we identified several other histamine-related curated gene sets that showed significant differential expression regardless of directionality. Of note, we did not find significant differential expression of gene sets relating to histamine secretion or histamine production in response to inflammation. This suggests that perhaps the response to histamine is influenced in ASD, rather than its secretion or production in the brain due to inflammation. However, the GO production of histamine in response to inflammation gene set does not include HDC. Furthermore, our analyses did not evaluate the expression of this gene, as it was too lowly expressed. Therefore, production of histamine may still be altered, but this may not be due to a change of the tested inflammatory genes within this gene set. The HS system may still influence neuroinflammation rather than the other way around. Further research with additional and larger studies are necessary to characterize HS gene expression in ASD and to determine the potential consequence of altered gene expression of this signaling system. In addition, more research is required to characterize the HS and its interplay with neuroinflammation in ASD. Although the role of the HS is well described in systemic inflammation, the interaction of the HS with neuroinflammation in particular, is still being characterized. Evidence suggests that HRH1 may particularly be involved, however, the interplay may be highly contextually dependent.[29]

At the transcriptome level, our analysis yielded eight significantly differentially expressed individual genes, a number lower than one previous report[54] yet comparable to another.[55] These variable results are likely due to differences in the stringency of confounder correction. Our conservative method may obscure some true findings, yet as demonstrated in Supplementary Figure 7, the potential role of population confounders (for example, race) and of implicit RNA quality differences is not fully addressed by adjusting for typical demographic or RNA-seq quality control measures. Given our goal of testing differential expression of specific genes of interest, we opted for a conservative approach to minimize Type I error. Less stringent confounder correction may be more appropriate when testing more global differential expression hypotheses. Further work is needed to resolve the trade-offs in approaches to control for confounders in RNA-seq.

The top differentially expressed genes included small nucleolar RNAs (snoRNAs), which are interesting given recent findings of differentially methylated regions in paternal sperm and postmortem ASD human brain mapping to genes in the SNORD family.[56] This warrants further investigation, but is beyond the scope of the current work. Interestingly, a deletion syndrome that included TUBE1 was associated with developmental delay and intellectual disability, suggesting the importance of this gene in brain development.[57] To our knowledge, associations of SNORA74A, SNORA53, SNORD17, SNORA54, SNORA74B, RP6-206I17.3 and SNORD114-23 have not been previously reported in ASD.

There are important limitations of our study beyond the sample size, which though small in contrast to genetic population studies, is comparable to prior brain tissue studies of autism. First, many factors can influence the quantification of gene expression in postmortem tissue, including cellular composition, RNA quality, sequencing depth and so on. Evaluating cell-type-specific methylation signatures, using the methods as in Jaffe et al.[58] indicated that there was no significant difference in the neuronal and microglial composition between our sample groups (data not shown), which is in agreement with results of another ASD DNA methylation study.[59] In addition, the inclusion of transcriptome-wide principal components in our analyses reduces the influence of these known confounders as well as other unknown confounders. Furthermore, ASD is a very heterogeneous disorder, with potentially varied causative factors among different subjects. To consider this, we evaluated the copy number variations in a subset of our samples and identified that only 1 out of 10 evaluated ASD samples had a large copy number variation, which was located in chromosome 1. Importantly, none of the evaluated histamine genes of interest are located on this chromosome. We also excluded samples with the 15q11.2-13.1 duplication syndrome in our replication analysis. Our use of a gene set test that is robust to small sample size and sample heterogeniety[49] should help mitigate some of the bias that may be introduced by each of these confounding effects. It would be especially important for future studies to investigate the expression of these histamine genes in neurons and microglia separately, to determine whether these genes are differentially expressed in primarily neurons, or microglia, or both. Deeper evaluation of the expression of each of the transcripts for these genes would also be helpful, as the shorter transcript of the HRH3, acts as an autoreceptor and therefore functions very differently than the longer isoforms. Deeper evaluation of the exon expression of HRH3 in our data indicate that only the longest transcripts were expressed. In the replication sample, the short autoreceptor transcript may have been expressed, although this is difficult to determine because it is unclear which exons are uniquely specific to the short form of the receptor. However, importantly all exons were more abundantly expressed in the ASD samples. Furthermore, expression studies cannot delineate between expression alterations due to causation or epiphenomena associated with the illness state of the condition studied.

HDC and HRH4 were not found to be highly expressed in the cortical regions evaluated in this study. Although this was expected, it would be interesting to evaluate the expression of these genes elsewhere. This is especially relevant for HDC, to determine whether it is under-expressed in the TMN, as it appears to be in particular cases of TS.[7, 20, 21, 60] This possibility of altered receptor expression due to altered HDC expression is consistent with unpublished findings indicating that expression of HRH3 is elevated in the basal ganglia in Hdc knockout mouse models of TS.[60] Another study of Hdc knockout mice found an increase in expression of HRH1 and HRH2 that was not significant, a significant decrease in expression of HRH3 in the hippocampus, and a significant increase in HRH3 expression in the TMN.[61] Also, a human transcriptome study of the striatum found an increase in HRH3 that was not significant in TS subjects.[62] Therefore, reduced HDC may induce regionally specific compensatory responses, although this requires further study. To our knowledge, the expression of HDC in the TMN has not yet been explored in ASD and it remains to be determined whether the expression of this gene is altered. Interestingly, higher plasma levels of histidine, the precursor of histamine, was identified in several ASD studies,[63, 64] suggestive of peripheral alterations of histidine catabolism and histamine production. Histidine readily passes the blood–brain barrier; therefore, an overabundance of histidine peripherally could lead to increased histamine production in the CNS by HDC. A high level of histamine in the brain is associated with neuroinflammation, which is also observed in ASD and is predicted to be involved in the pathogenesis of many neurological diseases, including ASD.[10, 65] Indeed TS may also have an inflammatory basis, as a recent study suggests enterovirus infection is associated with an increased incidence of Tic disorders;[66] and similar to ASD,[67] the incidence of allergies and asthma are strongly correlated with increased risk for TS.[68] Further research is required to determine what role the HS may or may not have in the enhanced neuroinflammation observed in each of these disorders. Indeed changes in either neuroinflammation or the HS may simply be reflective of changes in the other, and potentially due to compensatory effects, rather than indicative of a causative association.

In conclusion, the HS is involved in modulating cognition and behavior and has a role in neuroinflammation. Our study provides, to our knowledge, the first specific evaluation of the expression of a set of histaminergic genes representing histamine signaling in the brain of ASD subjects. Many questions remain about the involvement of this system in ASD and the potential overlap in pathogenesis of ASD with TS. However, our results, which are replicated in an independent ASD postmortem data set, provide a starting point for the investigation of this system in ASD.