Genome-wide identification of transcriptional targets of RORA reveals direct regulation of multiple genes associated with autism spectrum disorder

Background

We have recently identified the nuclear hormone receptor RORA (retinoic acid-related orphan receptor-alpha) as a novel candidate gene for autism spectrum disorder (ASD). Our independent cohort studies have consistently demonstrated the reduction of RORA transcript and/or protein levels in blood-derived lymphoblasts as well as in the postmortem prefrontal cortex and cerebellum of individuals with ASD. Moreover, we have also shown that RORA has the potential to be under negative and positive regulation by androgen and estrogen, respectively, suggesting the possibility that RORA may contribute to the male bias of ASD. However, little is known about transcriptional targets of this nuclear receptor, particularly in humans.

Methods

Here we identify transcriptional targets of RORA in human neuronal cells on a genome-wide level using chromatin immunoprecipitation (ChIP) with an anti-RORA antibody followed by whole-genome promoter array (chip) analysis. Selected potential targets of RORA were then validated by an independent ChIP followed by quantitative PCR analysis. To further demonstrate that reduced RORA expression results in reduced transcription of RORA targets, we determined the expression levels of the selected transcriptional targets in RORA-deficient human neuronal cells, as well as in postmortem brain tissues from individuals with ASD who exhibit reduced RORA expression.

Results

The ChIP-on-chip analysis reveals that RORA1, a major isoform of RORA protein in human brain, can be recruited to as many as 2,764 genomic locations corresponding to promoter regions of 2,544 genes across the human genome. Gene ontology analysis of this dataset of genes that are potentially directly regulated by RORA1 reveals statistically significant enrichment in biological functions negatively impacted in individuals with ASD, including neuronal differentiation, adhesion and survival, synaptogenesis, synaptic transmission and plasticity, and axonogenesis, as well as higher level functions such as development of the cortex and cerebellum, cognition, memory, and spatial learning. Independent ChIP-quantitative PCR analyses confirm binding of RORA1 to promoter regions of selected ASD-associated genes, including A2BP1, CYP19A1, ITPR1, NLGN1, and NTRK2, whose expression levels (in addition to HSD17B10) are also decreased in RORA1-repressed human neuronal cells and in prefrontal cortex tissues from individuals with ASD.

Conclusions

Findings from this study indicate that RORA transcriptionally regulates A2BP1, CYP19A1, HSD17B10, ITPR1, NLGN1, and NTRK2, and strongly suggest that reduction of this sex hormone-sensitive nuclear receptor in the brain causes dysregulated expression of these ASD-relevant genes as well as their associated pathways and functions which, in turn, may contribute to the underlying pathobiology of ASD.

Autism spectrum disorder (ASD) is a neurodevelopmental disorder that is characterized by deficits in social understanding and interactions, aberrant communication, and repetitive, stereotyped behaviors, often with restricted interests[1–4]. With an overall prevalence of 1 in 88 individuals in the United States[5], autism is inexplicably biased towards males by a ratio of at least 4:1, although some recent studies[6, 7] have reported ratios closer to 2:1, depending on the population studied. The consistently observed male bias, however, has prompted theories that elevated fetal or neonatal testosterone levels may be a risk factor for ASD[8], and the recent association of various autism traits in individuals who presented prenatally with elevated testosterone levels in amniotic fluid has supported this hypothesis[9, 10]. However, the molecular and physiological mechanism(s) that lead to elevated testosterone levels in individuals with ASD, both prenatally and postnatally, remain essentially unknown.

We have recently identified a novel autism candidate gene, retinoic acid-related (RAR) orphan receptor-alpha (RORA)[11] which is regulated by male and female hormones in a manner that may provide an explanation for the higher testosterone levels and, possibly, sex bias in ASD[12]. RORA is a ligand-dependent orphan nuclear hormone receptor that, in combination with co-regulator proteins, serves as a transcriptional regulator. Although RORA has never before been associated with ASD, our recent studies have demonstrated: reduced expression of RORA in lymphoblastoid cell lines (LCL) derived from individuals with autism[13]; increased methylation leading to reduced expression of RORA in the LCL from cases vs. sibling controls[11]; and decreased expression of RORA protein in the prefrontal cortex and the cerebellum of individuals with autism[11]. Together, these results link these molecular changes in RORA in blood-derived peripheral cells to molecular pathology in the brain tissues of individuals with autism.

These findings are notable because studies on the Rora-deficient staggerer mouse model indicate that Rora is involved in several processes potentially relevant to autism, including Purkinje cell differentiation[14, 15], cerebellar development[16, 17], protection of neurons against oxidative stress[18], suppression of inflammation[19], and regulation of circadian rhythm[20]. Indeed, the involvement of Purkinje cells and cerebellar abnormalities as well as neuroinflammation and oxidative stress in the autistic brain has been comprehensively discussed in a consensus report on the pathological role of the cerebellum in autism[21]. Recently, the proposed circadian dysfunction in ASD[22, 23] has also been supported by both genetic studies that have identified polymorphisms in “clock” (circadian regulator) genes[24] as well as gene expression analyses that identified RORA as one of the 15 differentially expressed circadian genes in a phenotypic subgroup of individuals with ASD with severe language impairment[13]. The known functions of Rora in the mouse model thus appear to be relevant to the observed pathological findings in humans with ASD. Moreover, behavioral studies on the staggerer mouse, primarily used as a model to study ataxia and dystonia[16], show that RORA is also associated with restricted behaviors reminiscent of autism, such as perseverative tendencies[25], limited maze patrolling[26], anomalous object exploration[27], and deficits in spatial learning[28]. Although there are currently no reported studies on the social behaviors of staggerer mice, it is clear that RORA is associated with at least some of the symptomatology and pathology of autism.

As a transcriptional regulator, RORA is known to bind DNA as a monomer or as a homodimer to hormone response elements upstream of target genes to modulate expression of those genes. In staggerer mice that exhibit spontaneous disruption of the Rora gene, Rora deficiency has been reported to cause aberrant expression of several genes involved in Purkinje cell differentiation (for example, Shh), calcium-mediated signaling (for example, Itpr1, Calb1, Cals1, Pcp4), as well as excitatory neurotransmission at glutamatergic synapses of Purkinje cells in the cerebellum (for example, Slc1a6, Grm1, Pcp2)[29]. In mouse skeletal muscle, Rora influences genes associated with lipid and carbohydrate metabolism (for example, Acsl4, Cd36, Hif1a), LXR nuclear receptor signaling (for example, LXRa, Srebp-1c), and Akt and AMPK signaling (for example, Akt2)[30]. Liver tissues of Rora-deficient mice have been shown to exhibit altered expression of several genes involved in triglyceride synthesis and storage (for example, Cidec, Cidea, Mogat1)[31], thus demonstrating pleiotropic effects of Rora depending on tissue type. However, little is known about transcriptional targets of RORA in humans, particularly in the central nervous system. We therefore sought to identify, at the genome-wide level, putative transcriptional targets of RORA in human neuronal cells, and to validate a functionally relevant subset of targets that may play a role in ASD. Since we have previously demonstrated decreased RORA expression in the frontal cortex of individuals with autism relative to that of unaffected controls[11, 12], we also investigated mRNA expression of the confirmed RORA targets in postmortem brain tissues of individuals diagnosed with autism in comparison with the expression of those genes in the brain of sex-matched and age-matched unaffected individuals.

Cell culture

The human neuroblastoma cells SH-SY5Y (ATCC, Manassas,VA, USA) were cultured in 1:1 Modified Eagle’s Medium (MEM) and Ham’s F12 media (MediaTech, Manassas, VA, USA) supplemented with 15% (v/v) fetal bovine serum (Atlanta Biologicals, Lawrenceville, GA, USA) and 1% penicillin/streptomycin (MediaTech). Cells were maintained at 37°C with 5% CO₂, and split 1:2 every 3 or 4 days when the cells reached ~80% confluency. For harvesting, the cells were treated with trypsin-ethylenediamine tetraacetic acid (MediaTech) for 2 to 3 minutes to release them from the surface of the culture flask. Complete growth medium was then added to the flask containing suspension cells to inactivate trypsin. Cells were transferred to a sterile centrifuge tube and pelleted by spinning at 800 rpm for 5 minutes at 4°C and gently washed twice with ice-cold PBS.

Frozen human brain tissues

Frozen postmortem prefrontal cortex (BA9/10) tissues from male individuals with autism (n = 3) and from age-matched typically developing males (n = 3) were obtained through the Autism Tissue Program (San Diego, CA, USA). All frozen brain tissues were preserved by the Harvard Brain Tissue Resource Center, Harvard Medical School. The list of human brain tissues used in this study is shown in Additional file1, along with information on the age of the donor, the postmortem interval (PMI) between death and tissue collection, and the cause of death when known.

Chromatin immunoprecipitation (ChIP)-on-chip analysis

Chromatin immunoprecipitation (ChIP) for promoter array analysis was performed using the Millipore EZ-ChIP Chromatin Immunoprecipitation Kit (Millipore, Billerica, MA, USA) according to the manufacturer's protocol. Briefly, the human neuronal cell line SH-SY5Y was cultured in complete growth medium in a T-175 flask until ~80% confluency (approximately 1.5 × 10⁷ cells). The medium was then carefully removed without disturbing the cells. The cells were fixed with 37% formaldehyde for exactly 10 minutes to crosslink chromatin. The crosslinking reaction was terminated by addition of 10% glycine, after which ice-cold PBS supplemented with protease inhibitor cocktail was added to the culture flask to wash and chill the cells. The crosslinked cells were dislodged from the flask by scraping gently with a cell scraper and then transferred to a pre-chilled centrifuge tube. Crosslinked cells were pelleted and nuclear extraction was performed using the Active Motif Nuclear Extract Kit (Active Motif, Carlsbad, CA, USA) according to the manufacturer's protocol. Nuclear pellets were resuspended in SDS Lysis Buffer and sonicated on wet ice to shear crosslinked chromatin to 200–1,000 bp using a Heat Systems-Ultrasonics W-380 sonicator (Heat Systems-Ultrasonics/Misonix, Farmingdale, NY, USA) set to 30% of maximum output power for 20 × 10 seconds with 1-minute intervals. Sonicated chromatin was divided into several aliquots for immunoprecipitation reactions and stored at −80°C until use. Each immunoprecipitation reaction was conducted using 1 μg of goat anti-RORA1 or normal goat IgG antibody (Santa Cruz Biotechnology, Santa Cruz, CA, USA). Immunoprecipitated chromatin was then reverse-crosslinked by adding 5 M NaCl to the final concentration of 0.2 M and incubated at 65°C overnight. DNA from the immunoprecipitated chromatin was isolated and purified using Millipore DNA purification columns (Millipore) and then submitted to the Genomics Core Facility at The George Washington University for analysis on Affymetrix human promoter tiling arrays (Affymetrix, Santa Clara, CA, USA).

Analysis of RORA binding sites on Affymetrix promoter tiling arrays

Purified DNA from the chromatin immunoprecipitated with anti-RORA1 or control IgG antibodies was amplified using the Whole Genome Amplification 2 Kit (Sigma Aldrich, St. Louis, MO, USA) according to the manufacturer’s protocol. The amplification product of each immunoprecipitation reaction was purified and hybridized on Affymetrix GeneChip Human Promoter 1.0R microarrays. Probes significantly enriched in RORA1-immunoprecipitated DNA relative to IgG-immunoprecipitated DNA were identified using the Partek Genomics Suite software (Partek, St. Louis, MO, USA), as described in detail below.

The Affymetrix Human Promoter 1.0R array contains over 4.6 × 10⁶ probes tiled over 25,500 promoter regions of annotated genes. Each gene promoter is thus interrogated on average by over 150 probes or “tiles”. For this study, DNA immunoprecipitated from SH-SY5Y cells with either a RORA1-specific antibody or a non-specific IgG as a control was hybridized on separate chips in triplicate, and the amount of probe hybridized to each element on the chips was determined. Partek Genomics Suite Software was used to analyze the intensities of the probe elements across the promoter array using the workflow recommended by Affymetrix for tiling arrays. In brief, the data normalization procedures include adjusting for probe sequence, RMA background correction, quantile normalization, and log (base 2) transformation of intensity data. Two-way analysis of variance was then used to determine differences between the RORA1-co-immunoprecipitated DNA and IgG-co-immunoprecipitated control DNA in hybridizing to specific probes. The MAT (model-based analysis of tiling arrays) peak-seeking algorithm[32] was subsequently used to detect enriched regions in the RORA1-co-immunoprecipated DNA versus IgG control samples. Nominal P value cutoff ≤0.05 (corresponding to false discovery rate (FDR) <7%), MAT score >0, minimum number of probes for each region >10, and average length for the ChIP fragment >600 bp were used as the peak detection parameters in this discovery phase of high-throughput screening for potential RORA targets. The microarray data from this study have been deposited into the Gene Expression Omnibus repository [GEO: GSE45756]. Positive MAT scores for a specific region indicate that RORA binding is enriched in that region relative to IgG controls. The genes mapping to the RORA-enriched promoter regions were identified using the gene annotation database (hg18) provided by Affymetrix. To further aid in the selection of gene promoters for confirmation analyses, an additional level of analysis was performed to determine the average intensity of probes across the regions hybridizing to the RORA-immunoprecipitated DNA relative to that hybridizing to the control DNA. These intensities were used to calculate the fold-change or enrichment across each gene promoter region for binding of the RORA-immunoprecipitated DNA over that of the IgG-precipitated DNA. Among the RORA-enriched promoter sites, the site corresponding to ITPR1 was used as a positive control since Itpr1 had been identified as a transcriptional target of Rora in mice[29]. Because the ITPR1 promoter region exhibited a fold-change of 1.29, we used this value as a minimum average fold-change in selecting target genes for further confirmation.

Rationale for selection of putative target genes for further confirmation

Of the 2,544 genes identified as putative targets of RORA by the ChIP-on-chip analysis (see Results), we selected ITPR1, CYP19A1, HSD17B10, A2BP1 (RBFOX1), NLGN1, and NTRK2 for further confirmation by ChIP-quanititative PCR analyses and by functional knockdown of RORA. Among these putative target genes, A2BP1, ITPR1, NLGN1, and NTRK2 have been previously identified as candidate genes for ASD[33–35]. The protein product of CYP19A1, also known as aromatase, is involved in the enzymatic conversion of androgen to estrogen via testosterone intermediates. Although CYP19A1 had previously been reported to be a target of RORA in human breast cancer cell lines[36], we recently demonstrated that its promoter region was a site of RORA binding in a human neuronal cell line[12]. Here, we wished to further confirm CYP19A1 as a transcriptional target of RORA by functional knockdown studies. Like CYP19A1, HSD17B10 also codes for an enzyme involved in the conversion of androgen to estrogen, but via 19-hydroxyandrostenedione, 19-oxyandrostenedione and estrone intermediates. Because of our earlier study demonstrating significant and correlated reductions in RORA and CYP19A1 protein levels in the autistic brain[12], we were interested in determining whether a deficiency in RORA could also impact this alternate biochemical pathway through down-regulation of HSD17B10, thus potentially reinforcing the buildup of androgen and reduction of estrogen in neuronal cells and tissues. The relevance of these genes to ASD is further elaborated in the Discussion.

Prediction of RORA binding elements

Putative binding elements of RORA in the promoter regions of A2BP1, CYP19A1, HSD17B10, ITPR1, NLGN1, and NTRK2 were predicted using PROMO 3.0[37, 38], JASPAR[39], and the EpiTect ChIP Search Portal of SABiosciences (Valencia, CA, USA). For each putative target gene, a total of three or four predicted transcription factor binding sites within 10 kb upstream of the transcription start site were selected for ChIP-quantitative PCR analyses.

ChIP-quantitative PCR analysis

Chromatin immunoprecipitation for ChIP-quantitative PCR analysis of RORA binding elements in the promoter region of putative RORA targets was conducted using the ChIP-IT Express Enzymatic Kit (Active Motif) following the manufacturer’s instructions. Briefly, confluent SH-SY5Y cells (approximately 1.5 × 10⁷ cells in a T-175 flask) were fixed with 10% formaldehyde for exactly 10 minutes and the fixation reaction was stopped by adding 10% glycine to the flask. The cells were washed with 10 ml ice-cold PBS for 5 seconds, then 6 ml ice-cold PBS supplemented with 0.5 mM (final concentration) phenylmethylsulfonylfluoride supplied in the kit was added to the culture flask to wash and chill the cells. The crosslinked cells were transferred from the flask to a pre-chilled 15 ml centrifuge tube by scraping gently with a cell scraper. Crosslinked cells were pelleted by centrifugation for 10 minutes at 2,500 rpm (720 relative centrifugal force) at 4°C, and then resuspended in 1 ml ice-cold Active Motif Lysis Buffer (Active Motif). The crosslinked cells were transferred to an ice-cold dounce homogenizer and the nuclei were released from the cells by douncing with a tight pestle. Optimal cell lysis was assessed under a phase contrast microscope using a hemacytometer. The nuclei were then transferred to an ice-cold 1.7 ml microcentrifuge tube and pelleted by centrifugation for 10 minutes at 5,000 rpm (2,400 relative centrifugal force) at 4°C. Chromatin was then isolated from the nuclear pellets and sheared into 150 to 1,000 bp fragments by incubating with 10 U/ml (final concentration) Enzymatic Shearing Cocktail (Active Motif) at 37°C for exactly 10 minutes. The enzymatic shearing reaction was stopped by adding ethylenediamine tetraacetic acid to a final concentration of 10 mM and chilling the reaction tube on ice for 10 minutes. Optimal shearing was assessed by agarose gel electrophoresis. For each ChIP reaction, enzymatically-sheared chromatin containing ~7 to 25 μg chromatin DNA was immunoprecipitated using 1 μg antibody and 25 μl Protein G Magnetic Beads (Active Motif). The list of antibodies for ChIP analyses is shown in Additional file2. Immunoprecipitated chromatin was reverse-crosslinked according to the ChIP-IT Express Enzymatic Kit protocol and DNA was isolated and purified from the chromatin using the ChIP DNA Purification Kit (Active Motif).

Real-time quantitative PCR analysis was conducted using the Applied Biosystems 7300 Real-Time PCR System (Applied Biosystems, Foster City, CA, USA) to determine the enrichment of each RORA binding element in immunoprecipitated DNA. Primers for quantitative PCR analysis were designed using Primer3 software[40] and were synthesized by Integrated DNA Technologies (Coralville, IA, USA). Input DNA was diluted into five 10-fold serial dilutions and included in quantitative PCR analyses. Relative enrichment values of RORA binding elements in each immunoprecipitated chromatin were calculated using standard curves obtained from the enrichment of RORA binding elements in the 10-fold serial dilutions of respective input DNA. The list of primers is shown in Additional file3.

Short hairpin RNA (shRNA) transfection

SH-SY5Y cells were cultured in a six-well culture plate containing complete growth medium without antibiotics (approximately 2.5 × 10⁵ cells per well) until cells were ~70% confluent. For each well, the cells were transfected with 2.50 μg RORA1 shRNA (Santa Cruz Biotechnology) or 2.50 μg negative control shRNA (Santa Cruz Biotechnology) using Lipofectamine LTX and PLUS reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's protocol. Briefly, shRNA (2.50 μg) was diluted in 500 μl Opti-MEM I Reduced Serum Medium (Invitrogen) and 1.25 μl PLUS reagent was added to the diluted shRNA solution. Lipofectamine LTX (25 μl) was then added to the shRNA-PLUS solution, incubated for 25 minutes at room temperature to form shRNA-Lipofectamine-LTX-PLUS complexes, and added to the cells. The cells were then incubated with the shRNA-Lipofectamine-LTX-PLUS complexes at 37°C and 5% CO₂ for 24 hours before harvesting. The list of shRNAs is shown in Additional file2.

RNA isolation and quantitative RT-PCR analysis

Quantitative RT-PCR analyses were performed as previously described[41]. Total RNA from the shRNA-transfected cells was isolated using TRIzol (Invitrogen) and purified using the RNeasy Mini Kit (Qiagen, Gaithersburg, MD, USA) following the manufacturers' instructions. Human brain tissues were homogenized in the Bullet Blender Homogenizer (Next Advance, Averill Park, NY, USA) using nuclease-free glass beads, after which total RNA from homogenized brain tissues was isolated using the RNeasy Mini Kit (Qiagen). A total of 1 μg purified total RNA was used for cDNA synthesis using the iScript cDNA Synthesis Kit (BioRad, Hercules, CA, USA) according to the manufacturer’s protocols. The reaction (20 μl) was incubated at 25°C for 5 minutes, followed by 42°C for 30 minutes, and ending with 85°C for 5 minutes. After reverse transcription, the cDNA reaction mixture was diluted to a volume of 50 μl with nuclease-free water and used as a template for quantitative PCR analyses. Real-time PCR analyses were conducted using the Applied Biosystems 7300 Real-Time PCR System (Applied Biosystems). Primers for quantitative RT-PCR analyses designed by Primer3 software are listed in Additional file3. The relative quantity of transcripts in each sample was calculated using standard curves based on the relative quantity of 18S RNA transcript in 10-fold serial dilutions of the respective sample.

Statistical analysis

A paired, two-sided, Student’s t test was performed to determine significance of the differences in numerical data obtained by quantitative PCR analyses. P <0.05 was considered statistically significant.

Hypergeometric distribution analyses[42] were used to determine the statistical significance of enrichment in autism candidate genes among the transcriptional targets identified by ChIP-on-chip analysis, relative to the genes present in AutDB[34], AutismKB[35], or a combination of the two databases of autism-associated genes.

Pathway and gene ontology analyses

Network prediction, functional, and gene ontology analyses were accomplished using licensed Ingenuity Pathway Analysis (IPA) and Pathway Studio7.0 software, as well as the open-access DAVID Bioinformatics Resources 6.7[43, 44]. For functional and network analyses, Fisher exact P values were calculated using the entire set of genes in the Ingenuity Knowledge Base as the reference set.

Identification of RORA transcriptional targets using whole-genome promoter array analysis and gene ontology analyses

ChIP-on-chip analysis revealed that a total of 2,764 probes (corresponding to 2,544 unique genes) were significantly enriched in RORA-immunoprecipitated chromatin relative to IgG-immunoprecipitated chromatin (P <0.05; FDR <7%). The complete list of gene-associated regions from the ChIP-on-chip analyses is shown in Additional file4. Gene ontology analysis of this complete list of genes using DAVID Bioinformatics Resources 6.7 shows that neuron differentiation and development as well as axonogenesis are significantly over-represented among biological processes shown in Annotation cluster 1 which has a highly significant Enrichment Score (ES) of 5.107 (Table 1). The dataset is also significantly enriched for genes involved in synaptic transmission and plasticity (Annotation cluster 4 with enrichment score 3.456) and post-synaptic density (Annotation cluster 5 with enrichment score 3.176). Furthermore, it is remarkable that all of genes in annotation clusters 1 to 5 are contained in AutDB and/or AutismKB databases. (See Additional file5 for a full list of genes and functions associated with annotation clusters 2 to 5.) This is noteworthy because, although this nuclear hormone receptor is known to have pleiotropic functions in different tissues, neurological functions are clearly enriched among the putative transcriptional targets of RORA within the context of a human neuronal cell line.

To further mine the dataset of potential transcriptional targets for higher level biological functions, disorders, and canonical pathways, we conducted network prediction and functional analysis on the RORA-bound genes using Ingenuity Pathway Analysis (IPA) software. Table 2 summarizes the functional and pathway analysis which reveals that nervous system development and function, neurological disease, and axonal guidance signaling are among the top five most significantly over-represented biological functions, disorders, and canonical pathways, respectively, associated with the RORA-enriched gene dataset. More detailed investigation into the specific neurological functions associated with the genes showed significant enrichment of genes involved in development of the brain and nervous system, axonogenesis, cell-cell adhesion of neurons, long-term potentiation of granule cells, neuritogenesis, and development of the cerebellum (Table 3). The genes associated with neurological disorders/diseases were enriched for schizophrenia, Huntington’s disease, movement disorders, dyskinesia, and seizure disorder, the latter three of which often present with the most severe subtype of ASD, which we have found to be associated with RORA deficiency[13]. Interestingly, a number of behaviors that are often disrupted or impaired in ASD, such as cognition, learning, circling behavior, emotional behavior, memory and spatial learning, are also significantly over-represented in the dataset of putative transcriptional targets of RORA. The complete list of genes associated with each of the biological functions, diseases/disorders, and behaviors described in Table 3 are provided in Additional file6.

To determine whether the ChIP-on-chip-identified putative targets of RORA were enriched in autism candidate genes, hypergeometric distribution analyses were performed to calculate P values for over-representation of genes from two autism gene databases (AutDB and AutismKB) within the target gene list. The results of the hypergeometric distribution analyses (Table 4) indicate that autism candidate genes, identified through various other studies, are indeed enriched among our transcriptional targets.

Network prediction of selected transcriptional targets to assess relevance to autism spectrum disorder (ASD)

We selected six potential gene targets, ITPR1, CYP19A1, A2BP1, HSD17B10, NLGN1, and NTRK2, for confirmation by independent ChIP-quantitative PCR and functional analyses based on the reasons given in Methods. The probe enrichment data and genomic locations for these genes are shown in Table 5. Functional analysis of these six genes using the Pathway Studio 7.0 network prediction program revealed an association of these genes with neurological disorders, including autism and ataxia, as well as autism-related neurological functions, including synaptogenesis, synaptic transmission, long-term potentiation, learning, and memory (Figure 1). Interestingly, other autism candidate genes, including NLGN3, NRXN1, RELN, and GABAA, were also included in the predicted gene network, thus supporting the investigation of our selected genes as ASD-relevant transcriptional targets of RORA.

Pathway analysis of selected potential transcriptional targets of RORA. A biological network was created using the Pathway Studio 7.0 program to identify biological functions and disorders associated with potential RORA transcriptional targets selected for confirmation (that is, ITPR1, CYP19A1, A2BP1, HSD17B10, NLGN1, and NTRK2 which are highlighted with a blue halo).

Full size image

Confirmation of RORA binding to selected transcriptional targets in human neuronal cells

To confirm that RORA protein binds to the promoter region of these six potential transcriptional targets, we conducted ChIP using anti-RORA1 or nonspecific IgG antibody, followed by quantitative PCR analysis to determine the enrichment of RORA binding sites in the respective promoter regions of the selected genes. Figure 2 shows that RORA binding sites in the promoter region(s) of each of these potential targets were significantly enriched in chromatin samples immunoprecipitated by anti-RORA1 antibody relative to IgG-immunoprecipitated chromatin, indicating that RORA protein indeed binds to promoter regions of A2BP1, CYP19A1, HSD17B10, ITPR1, NLGN1, and NTRK2.

Chromatin immunoprecipitation-quantitative PCR of RORA transcriptional targets. Chromatin immunoprecipitation followed by quantitative PCR analysis was conducted to determine whether RORA protein binds to A2BP1, CYP19A1, HSD17B10, ITPR1, NLGN1, and NTRK2 promoters. Chromatin was isolated from SH-SY5Y cells and immunoprecipitated with anti-RORA1 or IgG antibody. DNA was purified from the immunoprecipitated chromatin and quantitative real-time PCR analysis (n = 3) was performed to determine the enrichment of each RORA binding element in the promoter region of the selected potential RORA targets. Error bars indicate standard error of the mean. * P <0.05; ** P <0.01; ^#undetectable.

Full size image

RORA suppression reduces expression of selected transcriptional targets in human neuronal cells

To further validate that RORA regulates expression of A2BP1, CYP19A1, HSD17B10, ITPR1, NLGN1, and NTRK2, and to examine whether reduction of RORA expression leads to reduced expression of these selected RORA transcriptional targets, shRNA-mediated knockdown of RORA was performed in SH-SY5Y cells using RORA1 shRNA, and quantitative RT-PCR analyses were conducted to measure expression of RORA1 and the potential targets. Expression of RORA1 in SH-SY5Y cells was significantly down-regulated by approximately 50% in comparison with cells transfected with negative control shRNA, indicating that shRNA-mediated knockdown of RORA1 was successful. Notably, expression of all selected RORA transcriptional targets was significantly reduced (Figure 3), indicating that expression of these genes is regulated by RORA and that these genes indeed are RORA transcriptional targets in human neuronal cells. As a negative control, the expression level of GAPDH, which is not a transcriptional target of RORA, is unaffected by transfection of the SH-SY5Y cells with RORA1 shRNA.

Quantitative RT-PCR analysis of RORA short hairpin RNA-transfected SH-SY5Y cells. SH-SY5Y cells were transfected with RORA1 short hairpin RNA (shRNA) or negative control shRNA for 24 hours and quantitative RT-PCR analysis (n = 3) was performed to determine expression of RORA1 and potential transcriptional targets in the transfected cells. GAPDH was included as a negative control. Error bars indicate standard error of the mean. *P <0.05, **P <0.01.

Full size image

Expression of transcriptional targets of RORA is relatively reduced in the frontal cortex of individuals with autism in comparison with age-matched controls

We have previously reported reduction of RORA transcript and/or protein in four independent cohorts using LCL as well as tissues from the prefrontal cortex and the cerebellum of individuals diagnosed with autism[11–13]. To examine whether RORA reduction in the brain of individuals with autism may be associated with aberrant expression of the transcriptional targets identified in this study, a pilot study involving quantitative RT-PCR analysis of frozen postmortem prefrontal cortex tissues from individuals with autism (n = 3) and age-matched controls (n = 3) was performed. As shown in Table 6, the average expression of each of the six gene targets in the combined cases is reduced relative to the average expression levels of the respective genes in the combined controls, although the P values (all P >0.05) indicate that the differences are not statistically significant.

Based on reduced expression of RORA in LCL and in postmortem brain tissues of individuals with ASD versus unaffected controls, coupled with its known functions in cerebellar development and neuroprotection against inflammation and oxidative stress in mice, we postulated that this nuclear hormone receptor may be responsible for at least some of the pathobiology associated with ASD. In particular, our recent finding that RORA specifically binds the promoter region of aromatase, whose protein expression in the human brain is highly correlated with that of RORA, suggests a molecular explanation for the increased levels of testosterone observed in some individuals with ASD[12]. Specifically, we anticipate that downregulation of aromatase would lead to a buildup of its substrate, testosterone, with a corresponding reduction in estrogen synthesis. Moreover, the downregulation of RORA expression in response to androgen and upregulation in response to estrogen suggest a mechanism for the introduction of sex bias in ASD.

As a nuclear hormone receptor, RORA, in combination with various co-regulator proteins, can potentially regulate the transcription of a large number of gene targets. This study was therefore conducted in order to identify additional transcriptional targets of RORA, specifically within the context of a human neuronal cell model. To our knowledge, this is the first ChIP-on-chip study directed towards the comprehensive identification of transcriptional targets of RORA at the genome-wide level in any species. Consequently, in this pilot study aimed at discovery of genes that may be regulated by RORA, we chose to use less stringent P values (P ≤0.05; FDR <7%) for identification of RORA-enriched regions in order to capture as many potential gene targets of RORA as possible. Not surprisingly, our ChIP-on-chip analysis identified 2,764 promoter regions enriched for RORA binding sites, which corresponded to 2,544 unique genes. Interestingly, gene ontology analysis of this complete set of putative target genes revealed a strikingly high enrichment in genes associated with neuron differentiation and development, neuron projection morphogenesis, axonogenesis, and axon guidance in the top functional annotation cluster (enrichment score 5.107). Intriguingly, all of the genes associated with these processes are listed in the AutDB (SFARI Gene) and/or AutismKB databases of autism-associated genes, indicating relevance of our target gene set to ASD. Two additional highly significant functional annotation categories of genes are related to synaptic transmission and plasticity (enrichment score 3.456) and postsynaptic density (enrichment score 3.176). Again, all of the genes in both annotation categories are contained within one of the two aforementioned databases. The statistically significant enrichment of autism candidate genes within our dataset of genes identified by ChIP-on-chip analysis was also confirmed by hypergeometric distribution analyses which used the number of genes in either or both of the autism databases as the total number of interesting markers and the overlap between the genes in our dataset (that is, the selected markers) and the genes in either database as the number of selected interesting markers. The total number of annotated genes on the microarray represents the total number of general markers (that is, the population).

Biological network and pathway analysis using Ingenuity Pathway Analysis software further revealed that these potential transcriptional targets of RORA are significantly associated with nervous system development and function including development of cortex and cerebellum, axonogenesis, adhesion of neuronal cells, neuronal migration, neuritogenesis, neurotransmission, and synaptic density, all of which have been associated with autism[45–50]. It is also noteworthy that neurological disorders known to be co-morbid with autism, including schizophrenia[51, 52], movement disorder[53–57], and seizure disorder[58–61], were also over-represented among the genes in our dataset. Of particular relevance to behavioral deficits in ASD, genes related to cognition, learning, repetitive behaviors (for example, circling), memory, and spatial learning are also associated with these potential targets of RORA.

Among the selected targets of RORA investigated in this study, we confirmed ITPR1 and CYP19A1 as transcriptional targets in human neuronal cells by both ChIP-quantitative PCR and functional knockdown of RORA with shRNA. As mentioned earlier, both genes had been implicated as targets of RORA by earlier studies with Rora-deficient mice[29] as well as in human breast cancer[36]. Our current study, however, also reveals four novel transcriptional targets of RORA, thus expanding the repertoire of genes and pathways that may be impacted by RORA deficiency in humans. These novel target genes are A2BP1 (RBFOX1), HSD17B10, NLGN1, and NTRK2. Furthermore, our pilot study examining the expression of these six transcriptional targets of RORA in postmortem brain tissues from individuals with ASD and that of age-matched controls shows an average reduction of all targets in the autism samples, although the differences between combined cases and combined controls were not statistically significant, all exhibiting P >0.05. Surprisingly, the least significant differences in average expression levels were observed for RORA and CYP19A1. Inasmuch as we had previously detected significant differences (P <0.05) between the protein expression levels of both RORA and CYP19A1 in the frontal cortex of 22 cases versus 12 controls by confocal immunofluorescence analyses[12], we suggest that the lack of statistical significance for reduced expression of these two transcripts as well as the other targets in this study is probably due to the limitation in sample size (n = 3 for each group) coupled with the natural variability of gene expression among the individuals. Nevertheless, the data from this limited sampling of cases and controls suggest a trend towards reduced expression of transcriptional targets in autism brain tissues exhibiting reduced RORA. As discussed below, an independent review of the literature revealed that these six RORA targets were reliably associated with autism.

Relevance of A2BP1, CYP19A1, HSD17B10, ITPR1, NLGN1, and NTRK2 to the pathobiology of autism

A2BP1 (ataxin 2-binding protein 1), also known as RBFOX1, is an RNA-binding protein that regulates neuron-specific alternative splicing[62]. Several studies, including genetic and gene expression analyses, have established the link between this gene and autism. Using fluorescent in situ hybridization and quantitative PCR analyses, Martin and colleagues found a cryptic deletion of the A2BP1 gene in a female with autism, epilepsy, and global development delay[33]. Reduction of the A2BP1 transcript level was also observed in the lymphocytes of this individual, suggesting that the deletion causes aberrant expression of this gene in this autism case. By genotyping 27 SNPs across this gene in 206 parent–child trios, they identified two regions exhibiting a nominal association with autism. Moreover, a recent noise-reduction genome-wide association study (GWAS) of two autism GWAS datasets (with 597 and 696 families) from the Autism Genetic Resource Exchange revealed 1,535 significant linkage disequilibrium blocks overlapping 431 genes[63]. Interestingly, regions in the A2BP1 gene were identified among the most significant linkage disequilibrium blocks (P = 3.6×10^-5). Furthermore, Voineagu and colleagues also conducted transcriptomic analysis of postmortem frontal and temporal cortex tissues from 19 individuals with autism and 17 controls using microarrays[64]. They found that A2BP1 expression was significantly reduced in both frontal and temporal cortex tissues from individuals with autism relative to controls. Using high-throughput RNA sequencing and semi-quantitative RT-PCR analyses, they also demonstrated that splicing of A2BP1-dependent alternative exons in the brain of individuals with autism was dysregulated[64], suggesting that aberrant expression of A2BP1 results in dysregulation of alternative splicing in autism. Findings from several copy number variation (CNV) studies[65–70] and linkage studies[71–75] have also reported that A2BP1 is associated with autism. Aside from its association with ASD, recent studies on the physiological function of A2BP1/RBFOX1 demonstrate its involvement in synaptic transmission and neuronal excitation[76] as well as its role in the regulation of transcriptional networks involved in neuronal development[77].

CYP19A1 encodes aromatase, a key enzyme that converts androstenedione to estrone and testosterone to estradiol. Although its neuronal function is most often associated with the regulation of reproductive behaviors through the sex hormones, more recent studies have revealed unexpected functions of aromatase in the brain, including neurogenesis, neuronal differentiation, synaptic activity and plasticity, and protection against oxidative stress[50, 78–82]. Another interesting observation is the colocalization of aromatase with oxytocin in several regions of the rat brain, including the periventricular nucleus of the hypothalamus and the zona incerta, which suggests a role for this gene in oxytocinergic neurons in the limbic system[83], which may have implications for social cognition. CYP19A1 has also been associated with autism by genetic studies. Allen-Brady and colleagues conducted genome-wide screening of 70 families (192 individuals with autism and 461 unaffected relatives) and identified three regions shown to be highly significant in the linkage analysis[84]. Interestingly, one of the three regions is 15q21.1-15q22.2 (heterogeneity logarithm (base 10) of the odds = 5.31) where the CYP19A1 gene is located. Moreover, Chakrabarti and colleagues studied SNPs in 68 candidate genes for autism to identify common genetic variations associated with autistic traits and Asperger syndrome, a high-functioning subgroup of autism. A case–control association analysis of individuals with Asperger syndrome (n = 174) revealed regions in 14 genes, including CYP19A1, that showed a nominally significant association with Asperger syndrome[85].

HSD17B10 encodes 3-hydroxyacyl-CoA dehydrogenase type-2, a mitochondrial enzyme involved in mitochondrial integrity and oxidation of fatty acids and steroids. As mentioned earlier, like CYP19A1, HSD17B10 is involved in the conversion of androgens to estradiol via androstenedione and estrone intermediates. A reduction in HSD17B10 would thus be expected to lead to increased androgen and reduced estrogen, either of which may have serious consequences for brain development[86]. Indeed, reduced expression of HSD17B10 causes X-linked mental retardation, choreoathetosis, language impairment, and abnormal behavior[87–89]. These associated disorders are particularly interesting since we have identified RORA deficiency (specifically, reduced expression) in the phenotypic subtype of ASD associated with severe language impairment[13]. Using chromosome microarray analysis, Edens and colleagues reported duplications involving Xp11.22-p11.23, a region where HSD17B10 gene is located, in two females with autism and epilepsy[90].

ITPR1 encodes inositol 1,4,5-triphosphate receptor type 1, a ligand-gated ion channel that is activated by cytosolic calcium and inositol triphosphate, and that is involved in synaptogenesis and formation of dendritic contacts[91]. Computational modeling also suggests that the biochemical and biophysical properties of ITPR1 coupled with its high density and lower sensitivity to IP₃ in cerebellar Purkinje cells are critical to its postulated role in long-term depression[92]. With respect to ASD, transcriptomic analysis of frontal and temporal cortex tissues from 16 individuals with autism and 16 controls revealed down-regulation of ITPR1[64]. Bremer and colleagues used high-resolution whole genome array-based comparative genomic hybridization to screen 223 individuals with autism to identify gene dose alterations associated with autism susceptibility. They found regions in chromosome 3p25.3, where ITPR1 is located, significantly enriched in individuals with autism[93]. Another CNV study of 1,124 families, each of which included a single proband, his/her unaffected parents and, in most cases, an unaffected sibling, revealed as many as 234 CNV regions in several genes, including ITPR1, to be associated with autism[69]. Moreover, Levy and colleagues studied genomic copy-number variation in a large cohort of families with a single affected child and at least one unaffected sibling. Among significant regions identified to be de novo CNVs were regions containing the ITPR1 gene[94].

NLGN1 encodes neuroligin-1, a neuronal cell surface adhesion molecule that interacts with neurexin molecules. NLGN1 is involved in formation and remodeling of synapses in the central nervous system[95]. In addition to synaptogenesis, a recent study has also demonstrated a role for NLGN1, acting together with NRXN1β, in neuritogenesis and the formation of neuronal circuits upon activation of the fibroblast growth factor receptor-1[96]. In the context of ASD, Glessner and colleagues conducted a whole-genome CNV analysis of 859 individuals with autism and 1,409 typically developing children of European ancestry who were genotyped with ~550,000 SNP markers to comprehensively identify CNVs conferring susceptibility to autism[97]. Positive findings were evaluated in an independent cohort of 1,336 individuals with autism and 1,110 controls of European ancestry. Interestingly, a CNV region in the NLGN1 gene was significantly enriched in individuals with autism. Moreover, a family-based association analysis of 100 families with autism found a modest significant association at a region in NLGN1 gene[98]. A high-density SNP, genome-wide linkage study also showed that the chromosomal region 3q26.31-q27.3 was significantly associated with autism[99]. Deletion of this gene was also reported to cause impaired spatial memory and increased repetitive behavior[100].

NTRK2 encodes neurotrophic tyrosine kinase receptor type 2 (also known as TRKB) that is a catalytic receptor for several neurotrophins, including BDNF, NT4, and NT3. NTRK2 is involved in axon guidance signaling as well as synapse formation and plasticity[101–104]. Its role in ASD is implicated by genetic association analysis of 146 children with autism and 50 typically developing controls which identified six SNPs and multiple haplotypes in the NTRK2 gene[105]. This study also showed that BDNF levels were significantly increased in children with autism relative to that in controls, supporting a role of NTRK2/BDNF signaling as a susceptibility factor for the disorder.

Study limitations and future directions

Because little is known with regard to transcriptional targets of RORA, particularly in human neuronal cells, this ChIP-on-chip analysis was undertaken as a discovery tool to screen for potential transcriptional targets at the genome-wide level. As a result, we used a less stringent P-value of 0.05 as a filter to select probes enriched for RORA binding without applying correction for multiple testing (which is not part of the standard Partek analytical workflow for tiling arrays). However, all of the regions identified using P ≤0.05 are retained in the dataset even after imposing FDR ≤7%. Thus, although the list of 2,544 genes that are identified as putative targets of RORA should be interpreted with caution pending independent validation, we suggest that this list may be used as a basis for hypothesis generation with regard to the genes, pathways and functions that may be impacted by RORA deficiency. Another caveat is that the list of RORA targets does not replicate all of the several targets identified in the mouse cerebellum by ChIP-quantitative PCR analyses[29] other than ITPR1. We do not know the precise reason for this difference, but it is possible that the target selection may be a function of the specific tissues and cell types as well as species, since the mouse targets were identified in primary cerebellar cells, and the human targets were identified in the SH-SY5Y human neuroblastoma cell line. Indeed, co-regulator availability is also tissue-specific, and transcriptional regulation by RORA requires association with gene-specific and tissue-specific co-regulator proteins, thus adding another layer of complexity to target selection.

Another potential complication is the possibility that target selection may be directed by circadian rhythm, as RORA is a regulator of the circadian cycle. Thus, it would be interesting to conduct ChIP-on-chip analysis on cells synchronized according to stage in the circadian cycle. Finally, as noted earlier, although the average expression levels of the RORA targets were reduced in the frontal cortex of individuals with ASD relative to that of controls as expected, the differences were not statistically significant. Notably, the protein expression levels of the genes exhibiting the least significant expression differences (that is, highest P values) in this study, RORA and CYP19A1, were found to be significantly different in our previous study in which 12 autism tissue samples were compared with 22 controls[12], thus emphasizing the need for an expanded investigation of RORA targets involving additional brain tissues from individuals with ASD and controls, including samples from different brain regions.

In summary, using a high-throughput genome-wide method for identifying potential transcriptional targets of RORA followed by independent validation of six targets with functional relevance to ASD, we have demonstrated that RORA deficiency can lead to a cascade of transcriptional deregulation that impacts a number of genes known to be associated with ASD. Figure 4 summarizes the principal findings from this study.

Possible downstream consequences of deregulation of the six confirmed transcriptional targets of RORA. This figure summarizes findings from this study. ChIP, chromatin immunoprecipitation; qPCR, quantitative PCR.

Full size image

Using ChIP-on-chip analyses, we identify at the genome-wide level over 2,500 potential transcriptional targets of RORA, a nuclear hormone receptor, and confirm six of these targets (four of them novel with respect to regulation by RORA) in human neuronal cells by ChIP-quantitative PCR and functional analyses. Moreover, we demonstrate that the transcript levels of these target genes are reduced in the postmortem frontal cortex of autistic males who exhibit reduced RORA expression relative to that of control males, further implicating RORA as an upstream regulator of these genes in human brain. The significance of these findings is the likelihood that any mechanism that contributes to reduced RORA expression, including elevated levels of male hormones[12] and epigenetic modifications that may be the result of as yet unknown environmental factors[11], may lead to increased risk for autism based on the aberrant transcription of downstream targets of RORA, which we show are significantly enriched for autism candidate genes based on hypergeometric distribution analyses of our dataset in comparison with established databases of autism genes.

A2BP1:

Ataxin 2 binding protein 1

Acsl4:

Acyl-CoA synthetase long-chain family member 4

Akt2:

V-akt murine thymoma viral oncogene homolog 2

ASD:

Autism spectrum disorder

BA:

Brodmann area

BDNF:

brain-derived neurotrophic factor

bp:

base pairs

Calb1:

Calbindin 1, 28kDa

Cals:

Callose synthase

Cd36:

Trombospondin receptor

ChIP:

chromatin immunoprecipitation

ChIP-on-chip:

chromatin immunoprecipitation followed by microarray analysis

Cidea:

Cell death-inducing DFFA-like effector a

Cidec:

Cell death-inducing DFFA-like effector c

CNV:

Copy number variant(s)

CYP19A1:

Cytochrome P450, family 19, subfamily A, polypeptide 1

FDR:

False discovery rate

GABA:

Gamma aminobutyric acid

Grm1:

Glutamate receptor, metabotropic 1

Hif1a:

Hypoxia inducible factor 1, alpha subunit

HSD17B10:

Hydroxysteroid (17-β) dehydrogenase 10

IPA:

Ingenuity Pathway Analysis

ITPR1:

Inositol 1,4,5-trisphosphate receptor, type 1

LCL:

Lymphoblastoid cell line

LXRa:

Liver X receptor alpha

MEM:

modified Eagle’s medium

Mogat1:

Monoacylglycerol O-acyltransferase 1

NLGN1:

Neuroligin 1

NLGN3:

Neuroligin 3

NRXN1:

Neurexin 1

NT3:

Neurotransmitter transporter 3

NT4:

Neurotransmitter transporter 4

NTRK2:

Neurotrophic tyrosine kinase, receptor, type 2

PBS:

phosphate-buffered saline

Pcp2:

Purkinje cell protein 2

Pcp4:

Purkinje cell protein 4

PCR:

polymerase chain reaction

RBFOX1:

RNA-binding protein, fox-1 homolog

RELN:

Reelin

RORA:

Retinoic acid receptor-related orphan receptor alpha

RORA1:

Retinoic acid receptor-related orphan receptor alpha, isoform 1

RT:

Reverse transcription

Shh:

Sonic hedgehog

shRNA:

Short hairpin RNA

SH-SY5Y:

Human neuroblastoma cell line

Slc1a6:

Solute carrier protein 1, member 6

SNP:

Single-nucleotide polymorphism

Srebp-1c:

Sterol regulatory element-binding protein 1c

The authors wish to thank the families of the brain donors for their invaluable gift to autism research and the Autism Tissue Program (San Diego, CA, USA) for making frozen human brain tissues available for this work through the Harvard Brain Tissue Resource Center (Harvard University, Boston, MA, USA). They also thank Ms. Yi Hong for her assistance with the Affymetrix promoter array hybridization and analysis, and Ms. Jing Jing for help with RNA isolation. The authors are also grateful to Dr. Yinglei Lai (Department of Statistics at The George Washington University) for conducting the hypergeometric distribution analyses.

This study was supported in part by a generous gift from the LIFE Foundation (Aspen, CO, USA) and by an intramural grant from The George Washington University (Medical Faculty Associates award). Neither of the funding sources played any role in the study design, collection, analysis, and interpretation of data, writing of the manuscript, or decision to submit this study for publication. TS was a predoctoral student in the Institute for Biomedical Sciences at the George Washington University, who was supported by the Higher Educational Strategic Scholarship for Frontier Research Network (SFR scholarship) from the Office of the Commission on Higher Education of the Royal Thai Government, Thailand, through the Faculty of Allied Health Sciences, Chulalongkorn University.

Author notes

1. • Tewarit Sarachana

   Present address: Center for Biologics Evaluation and Research, Food and Drug Administration, 8800 Rockville Pike, Bethesda, MD, 20892, USA

Affiliations

1. Department of Biochemistry and Molecular Medicine, The George Washington University School of Medicine and Health Sciences, 2300 I Street NW, Washington, DC, 20037, USA

   • Tewarit Sarachana
   •  & Valerie W Hu

2. Center for Biologics Evaluation and Research, Food and Drug Administration, 8800 Rockville Pike, Bethesda, MD, 20892, USA

   • Tewarit Sarachana

Corresponding author

Correspondence to Valerie W Hu.

Competing interests

The authors declare that they have no competing interests.

Authors' contributions

TS and VWH conceived of the study, contributed to the study design, performed the data analyses, and prepared the manuscript. TS conducted the experiments. All authors read and approved the final manuscript.

Reprints and Permissions