Global analysis of gene expression mediated by OX1 orexin receptor signaling in a hypothalamic cell line

The orexins and their cognate G-protein coupled receptors have been widely studied due to their associations with various behaviors and cellular processes. However, the detailed downstream signaling cascades that mediate these effects are not completely understood. We report the generation of a neuronal model cell line that stably expresses the OX1 orexin receptor (OX1) and an RNA-Seq analysis of changes in gene expression seen upon receptor activation. Upon treatment with orexin, several families of related transcription factors are transcriptionally regulated, including the early growth response genes (Egr), the Kruppel-like factors (Klf), and the Nr4a subgroup of nuclear hormone receptors. Furthermore, some of the transcriptional effects observed have also been seen in data from in vivo sleep deprivation microarray studies, supporting the physiological relevance of the data set. Additionally, inhibition of one of the most highly regulated genes, serum and glucocorticoid-regulated kinase 1 (Sgk1), resulted in the diminished orexin-dependent induction of a subset of genes. These results provide new insight into the molecular signaling events that occur during OX1 signaling and support a role for orexin signaling in the stimulation of wakefulness during sleep deprivation studies.

In addition to efforts demonstrating the behavioral effects of the orexin system, a number of studies have addressed the intracellular molecular signaling events that occur in response to orexin receptor activation [26][27][28][29][30]. Upon [44][45][46]. While these signaling studies have provided much insight, they are largely focused on the upstream, canonical GPCR signaling pathways, leaving the detailed downstream signaling cascades unclear. Furthermore, most of these signaling studies have been done in heterologous expression systems that are non-neuronal, resulting in varied responses and raising questions of applicability to orexin receptor signaling in a neural context [37,[47][48][49][50][51][52][53][54][55].
In this study, we report the generation a of a neuronal cell line that stably expresses OX 1 and the transcriptional profile seen upon receptor activation, as determined by RNA-Seq. OX 1 activation resulted in the differential regulation of a large set of genes, several of which have previously been shown to be similarly regulated by sleep deprivation (SD), in vivo. Additionally, the downstream role of one of the more highly regulated genes, Sgk1, was further characterized.

Analysis of mRNA transcripts by qPCR
Cells were grown to near confluence in 60mm culture dishes. For orexin treatments, culture media was replaced with fresh, warm media containing 100nM orexin A (OxA, Tocris) and incubated for 3h at 37˚C. For OX 1 inhibition, cells were pretreated with media containing 3μM SB-334867 (Tocris) for 10 min prior to adding OxA. RNA was purified from cells with the RNeasy Plus Mini Kit (Qiagen). First strand cDNAs were synthesized with the iScript cDNA Synthesis Kit (BioRad) in a 20μl reaction using 1μg RNA. The qPCR reactions were done in triplicate on the StepOnePlus real time PCR system (Applied Biosystems) with Taq-Man gene expression assays (Applied Biosystems) using the following conditions: 1μl cDNA, 1μl TaqMan probe (S1 Table), 10ul TaqMan Gene Expression Master Mix (Applied Biosystems), and 8ul nuclease-free water. Thermal cycling conditions were 95˚C for 10m, then 40 cycles of 95˚C for 15s and 60˚C for 60s.
For the OX 1 inhibition assay, 4-fold serial dilutions of SB-334867 were done in DMSO, starting from 2.5mM. 7μl of each serial dilution was then added to 293μl 1X Stimulation Buffer. These dilutions were then added in a 1:1 ratio (v/v) to each well of a low-volume 384-well plate (white) already containing 20,000 GT1-7-OX 1 cells per well in 6μl 1X Stimulation Buffer. After a 30 minute incubation at 37˚C, 2μl OxA (varying concentrations in 1X Stimulation Buffer) was added to each well (final [DMSO] = 1%). After 45 minutes at 37˚C, the HTRF reagents were added (3μl each) and the plate was incubated for 1 hour at room temperature. Data were acquired on a Tecan Infinite M1000 Pro plate reader. All experiments were done in triplicate.

Generation of GT1-7-OX 1 stable cells
The gene for human OX 1 (Genecopoeia, EX-U0062-M02) was subcloned into pCDH-CMV-MCS-EF1-copGFP (Systems Biosciences) at the XbaI and BamHI restriction sites. Lentiviral particles were generated in HEK293T cells (ATCC) by co-transfecting the pCDH-OX 1 plasmid with the pPACKH1 HIV Lentivector Packaging Kit (Systems Biosciences) using FUGENE HD Transfection Reagent (Promega). The lentiviral particles were then concentrated with PEG-it Virus Precipitation Solution (Systems Biosciences), resuspended in 400μl DPBS + 25mM HEPES, and stored at -80˚C until ready for use. For the viral transduction, GT1-7 cells were cultured in 6-well plates to 60% confluence. Various volumes of lentivirus (50, 100, 200μl) were added directly to the culture media and mixed. After 72h, reporter gene expression and cell viability were analyzed via microscopy and the cells that demonstrated the highest levels of GFP expression with minimal cell toxicity were expanded for further analysis.
RNA-Seq library construction and sequencing GT1-7-OX 1 cells were plated in 75cm 2 culture flasks and grown to near confluence. For cell treatments, culture media was replaced with fresh, warm media containing vehicle (H 2 O) for 8 hours, or 200nM OxA for 3 or 8 hours. RNA was isolated with TRIzol Reagent (Life Technologies) according to manufacturer's protocol, including the addition of 10μg RNase-free glycogen (ThermoFisher). The RNA samples were then treated with DNase (New England Biolabs) to remove any genomic DNA contamination and then cleaned up with the Purelink RNA Micro Kit (Invitrogen). This process was repeated twice for n = 3 per condition. The DNasetreated Total RNA (250ng) was depleted of ribosomal RNA using the TruSeq Stranded Total RNA kit (Illumina) and quality assessed on an Agilent 2100 Bioanalyzer to confirm that 18S and 28S rRNA peaks were depleted. The rRNA-depleted RNA was converted to dsDNA libraries by following the TruSeq Stranded Total RNA sample prep kit user guide. Briefly, the RNA was fragmented, converted to cDNA, and ligated with adaptors. The adaptor-ligated DNA was then PCR amplified using 11 cycles to generate the final libraries. The final libraries were size selected and purified using 1.0 x Ampure XP beads (Beckman Coulter) then validated by the Bioanalyzer and qPCR quantified using primers that recognize the Illumina adaptors. The libraries were then pooled at equimolar ratios, quantified using qPCR (quantification of only the adaptor-ligated libraries) and loaded onto the NextSeq 500 flow cell (Illumina) at 1.8pM final concentration. Demultiplexed and quality filtered raw reads (fastq) generated from the NextSeq 500 were trimmed (adaptor sequences) using Trimmomatic, version 0.35 [56] and aligned to the reference genome (UCSC-mm10) using STAR, version 2.5.2a [57]. HTSeq-count (version 0.6.0) was used to generate gene counts and differential gene expression analysis was performed using DESeq2 (version 1.10.1, R version: 3.2.3) [58], comparing the OxAtreated samples to those treated with vehicle. The principle component analysis was performed via the plotPCA function in DESeq2, using the regularized log-transformed values of the 500 genes that were the most variable across all samples. In order to more closely identify the relationship between each sample and every other sample, the Euclidean distance between each pair of samples was calculated using the log-transformed values of the complete data set. Complete linkage clustering was then used to generate a sample-to-sample distance heatmap, via the pheatmap package in R source. For statistical analyses, raw counts for the two conditions of interest were imported into DESeq2 and transformed using the negative binomial Wald test. Adjusted p-values were generated via the Benjamini-Hochberg procedure.

Promoter analysis
Gene symbols of the differentially expressed genes identified in the RNA-Seq data were entered into the DAVID Gene ID Conversion tool (https://david.ncifcrf.gov/) [59,60] under the settings "OFFICIAL_GENE_SYMBOL" (input) and "REFSEQ_MRNA" (output). From the complete set of 332 genes that were differentially regulated 2-fold or greater by OX 1 signaling at 3 or 8h (adj. p-values <0.05, log 2 Fold Change (log 2 FC) >1.0 or < -1.0), there were 31 official gene symbols that either were not recognized by DAVID, or could not be converted into a RefSeq_mRNA ID recognized by PSCAN, and could not be included in the PSCAN analysis. The complete list of gene ID conversions, including the unrecognized gene ID's, can be found in S2 Table. The remaining gene ID's were entered into the PSCAN user interface (http://159.149.160.88/pscan/), and run with the following settings: Mus musculus (organism), -450 +50 (region), Jaspar 2016 (Descriptors).
Comparison between GT1-7-OX 1 RNA-Seq and SD microarray metaanalysis The RNA-Seq data from this work were cross-compared to results from an SD microarray meta-analysis [61]. Of the 91 SD-related mouse genes presented in the Wang, et al. study, five did not have Gene ID's that correlated to our data set (2310076G05Rik, 3110003A17Rik, 4932442K08Rik, C330006P03Rik, and D930028F11Rik) and were not included in the comparison. Notably, our data set does not distinguish between the long and short isoforms of Rbm3, which were oppositely regulated by SD in the meta-analysis. Therefore, we also excluded Rbm3 for this comparison.
Sgk1 inhibition assay GT1-7-OX 1 cells were plated in 75cm 2 culture flasks and grown to 70-90% confluence. Growth media was replaced with fresh, warm media containing 1.0μM GSK-650394 (Apexbio Technology), or DMSO vehicle. After 30 minutes at 37˚C, either H 2 O vehicle or OxA was added to the media at 200nM (final). After another 3 hours of incubation at 37˚C, RNA was purified from cells with the RNeasy Plus Mini Kit (Qiagen). First strand cDNAs were synthesized with the iScript cDNA Synthesis Kit (BioRad) in 400μl reactions using 20μg RNA. The qPCR reactions were done with PowerUp SYBR Green Master Mix (Applied Biosystems) and PrimeTime qPCR primer pairs (IDT, S3 Table) in 20μl reactions (10μl SYBR, 2μl primer pair (500nM, final), 0.5μl cDNA), in triplicate, on the StepOnePlus real time PCR system (Applied Biosystems). Cycling conditions were 50˚C for 2 minutes, 95˚C for 2 minutes, then 40 cycles of 95˚C for 15 seconds, 50˚C for 1 minute. In addition to 89 genes of interest, primer pairs targeting seven housekeeping genes were included (Actb, B2m, Gusb, Polr2a, Ppia, Rplp0, and Tbp). As it demonstrated the strongest stability amongst treatments, with a geNorm M value <0.2 (as determined via qbase+ software, version 3.1), B2m was used as the endogenous control for data analysis. Data were analyzed by the 2 -ΔΔC T method and represented as fold-change over control samples.

Characterization of orexin receptor-expressing cell lines
The initial goal of this study was to identify a cell line that would be a reasonable model in which to analyze OX 1 signaling. To that end, several cell lines, originating from different species, that have been reported to express one or both of the orexin receptors, endogenously, were acquired [51, [62][63][64]. Each cell line was screened for the presence of OX 1 and OX 2 mRNA by qPCR with a set of probes designed to span various exons (Table 1). While some amplification was observed sporadically, the high Ct values and inconsistency between probes did not clearly demonstrate the presence of orexin receptor transcripts in any of the cell lines tested. In order to look for the presence of functional orexin receptors, the IP-One HTRF assay was employed. This assay is a FRET-based immunoassay that measures accumulation of inositol monophosphate (IP1) upon activation of the phospholipase C pathway and is often used as a measure of G q -coupled GPCR activation. At baseline levels, the kit components (FRET donor and acceptor) are bound together, resulting in a high HTRF ratio. Upon ligand-induced GPCR activation, production of cellular IP1 is stimulated. This native, unlabeled IP1 then displaces acceptor-labeled IP1, disrupting the proximity of the donor and acceptor molecules and resulting in decreased HTRF ratios. While this assay is commonly used to quantify IP1 production (via a standard curve), we utilized it as a simple measure of receptor functionality and did not perform this quantitation. As shown in Fig 1, none of the cell lines tested demonstrated detectable GPCR activation in response to increasing concentrations of the OxA peptide, whereas a CHO-based cell line stably expressing OX 1 (CHO-OX 1 ) provided a robust response.

Generation of GT1-7 cells stably expressing OX 1
In the absence of an existing cell line that expresses OX 1 endogenously, we turned to the construction of a neuronal cell line that would express recombinant OX 1 . As the orexin receptors are known to be highly expressed in the hypothalamus [65,66] and GT1-7 is a mouse cell line derived from hypothalamic neurons, it seemed reasonable to use this as the parental cell line. A lentiviral transduction system was used to stably incorporate human OX 1 into GT1-7 cells. The presence of the transcript was verified by qPCR (Fig 2A) and the presence of functional receptor was demonstrated via IP-One HTRF assay ( Fig 2B). The functional response to OxA in these cells was inhibited with SB-334867 [67], a selective inhibitor of OX 1 (Fig 2C).

RNA Seq analysis
In order to assess changes in gene expression brought about by OX 1 signaling, GT1-7-OX 1 cells were treated with OxA for 3 or 8 hours. RNA-Seq analysis was then used to identify transcripts regulated differentially compared to cells treated with vehicle. To evaluate the overall relationships between samples and test for batch effects, a principal component analysis was conducted (Fig A in S1 Fig). In addition, for a more complete analysis of how each sample compares to every other sample, a clustering analysis was performed (Fig B in S1 Fig). In each case, strong clustering among replicates and treatment groups was evident. Additionally, more than 3 x 10 7 mapped reads were obtained per sample, representing approximately 90% of total reads, while only 0.4% of reads mapped to a ribosomal RNA reference, indicating minimal rRNA contamination. The data discussed in this publication have been deposited in NCBI's Gene Expression Omnibus [68] and are accessible through GEO Series accession number GSE99690 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE99690).
As an initial characterization of the recombinant model, the basal expression of known neural marker genes was analyzed. Using the average transcripts per million (TPM) values of the vehicle-treated control samples as a measure of expression level (Avg. TPM >2.0), GT1-7-OX 1 cells expressed a number of neuronal marker genes [69], but not glial marker genes, indicating  (Table 2). Of note, as GT1-7 cells were isolated from mouse hypothalamic tumor cells designed to express the SV-40 T-antigen under the control of the gonadotropin releasing hormone (GnRH) promoter sequence [70], we expected to see elevated levels of GnRH in these cells. As expected, these cells express GnRH at very high levels (Avg. TPM = 3924).
In order to validate the changes in gene expression demonstrated by the RNA-Seq data, the experiment was repeated and a subset of genes was analyzed by qPCR. The levels of gene expression measured were consistent with the RNA-Seq data. (Fig 4A). In a second, smaller validation experiment, SB-334867 was used to inhibit the orexin-dependent differential expression of another subset of genes, demonstrating that OX 1 activation was required for these changes in transcription (Fig 4B).
The results were also compared to a microarray study of OX 1 -expressing HEK293 cells published previously [71]. Of the genes that were up-regulated 2-fold or greater by OX 1 activation in HEK293 cells, 346 had gene symbols that were present in the GT1-7-OX 1 RNA-Seq data. Of these, only 24 (6.9%) were similarly up-regulated in GTI-7-OX 1 cells (2-fold or greater at 3 or 8h, adj. p <0.05). Of the genes that were down-regulated 2-fold or greater by OX 1 activation in HEK293 cells, 370 had gene symbols that were present in the GT1-7-OX 1 RNA-Seq data, with only 1 (0.3%) that was similarly down-regulated (2-fold or greater at 3 or 8h, adj. p <0.05).
Conversely, 5 of the 370 (1.4%) down-regulated genes were actually up-regulated 2-fold or greater at either 3 or 8h in GT1-7-OX 1 cells. So, while some similarities were observed, the  OX 1 orexin receptor signaling in a hypothalamic cell line data were poorly replicated in the two systems, highlighting the importance of using an appropriate model for studying signal transduction, in vitro. Nevertheless, we did identify a set of genes comparably regulated by OX 1 in both HEK293 and GT1-7 cells, strongly supporting a role for these genes in OX 1 signal transduction (S4 Table).
Transcription factors regulated by OX 1 signaling GPCR-mediated signaling usually results in the activation or suppression of transcription factors, which then drive downstream changes in global gene expression. The data were thus searched for OX 1 -regulated genes known to encode transcription factors. Several families of related transcription factors, mostly immediate early genes, were highly regulated by orexin. These include the early growth response genes (Egr1, Egr2, Egr3), AP-1 (Fos, Fosb, Junb, Fosl2), high mobility group superfamily A (Hmga1, Hmga2), the Nr4a family of nuclear hormone receptors (Nr4a1, Nr4a3), the Id family of transcriptional repressors (Id1, Id2, Id3, Id4)  In order to identify the transcription factors likely responsible for driving the more downstream changes in gene expression, we performed an in silico analysis with PSCAN, an automated program that examines the promoter regions of regulated genes in order to identify common transcription factor binding sites [72]. The PSCAN analysis identified several putative transcription factors but, when combined with the expression data, strongly suggested a role for two distinct families of transcription factors in the regulation of downstream genes, the early growth response (EGR) proteins and Kruppel-like factors (KLF), in addition to the Myc transcription factor (Table 4). Of note, the EGR genes are heavily associated with neural plasticity and memory [73] as well as enhanced long-term potentiation that impacts relapse in drug-related reward behaviors [74,75], processes in which orexin signaling is known to be involved [12,[76][77][78].

Sleep deprivation-related genes
To assess the physiological relevance of the data presented above, we searched the literature to compare our results with those from animal studies focused on orexin-related behaviors. Interestingly, there are striking similarities between our data set and those from an SD computational meta-analysis aimed at identifying highly conserved SD-related genes [61]. In this study, the authors combined and analyzed data from all available SD microarray studies spanning four species (mouse, rat, sparrow, and fruit fly). They found that SD resulted in a highly conserved (across at least 3 species) induction of Egr1, Nr4a1, and Arc, all of which were induced strongly by OxA in GT1-7-OX 1 cells. Furthermore, the authors reported a set of 91 mouse genes that were differentially regulated during short-term SD (zeitgeber time 0-6h). Of the 90 SD-related genes that could be cross-referenced to our data (see Materials and Methods), 45.6% (41/90) were differentially regulated in response to 3h OX 1 signaling at statistically significant levels (adj. p-value <0.05, Table 5). Additionally, 33.3% of these genes (30/90) were regulated in the same fashion (up or down) in both studies. These genes notably included Egr1, Egr2, Egr3, Nr4a1, Nr4a3, Arc, and Sgk1. Since it is known that SD results in elevated levels of OxA in the locus coeruleus and hypothalamus [79], increased hypothalamic OxA immunoreactivity [80,81], and increased expression of c-fos in orexinergic neurons [82], these parallels between the data presented here and gene regulation in SD suggest strongly that the events occurring in OX 1 -expressing GT1-7 cells are of significant physiological relevance.

Role of Sgk1 in orexin signaling
One of the most highly orexin-responsive genes in GT1-7-OX 1 cells, serum/glucocorticoidregulated kinase 1 (Sgk1), is also highly up-regulated by SD. This transcript was of particular interest due to its roles in neuronal excitation and synaptic plasticity [83][84][85][86]. Sgk1 is expressed in all tissues of the body, including the brain, and regulates numerous ion channels, molecular transporters, and signaling proteins [87][88][89][90][91]. Transcription of Sgk1 mRNA has been shown to be induced by several stimuli and, as it relates to orexin, by exposure to drugs of abuse [92,93] or fasting conditions [94][95][96]. The Sgk1 transcript was highly induced in response to OxA at both 3 and 8 hours (9-fold and 4-fold, respectively), with adj. p-values approaching zero. In order to examine the potential role of Sgk1 in orexin-regulated gene expression, an Sgk1 inhibitor, GSK-650394 [97], was added to GT1-7-OX 1 cells prior to the addition of OxA. The effects of the Sgk1 inhibition on transcription of 89 of the most highly OX 1 -regulated genes were determined by qPCR (S2 Fig). From this set of genes, eleven showed reduced levels of induction when pretreated with GSK-650394 ( Fig 5). These data argue that Sgk1 regulates these genes, possibly through phosphorylation of a transcription factor(s) that targets them. In an attempt to identify this putative factor(s), the PSCAN analysis was repeated with this set of OX 1 orexin receptor signaling in a hypothalamic cell line genes (S5 Table). As seen in the earlier results, the Egr family of transcription factors was again identified. Interestingly, one of the other transcription factors identified, Sp1, is a known substrate of Sgk1 [98] and has been shown to regulate the transcription of Ldlr [99,100], Dok7 [101], Rara [102,103], and Cldn4 [104,105].

Discussion
Because the orexin system has been shown to regulate behavior primarily via its actions in the central nervous system, studying orexin receptor signaling in a neuronal context is of particular interest. To that end, a recombinant model stably expressing OX 1 was generated in GT1-7 cells, a mouse cell line derived from GnRH-expressing neurons of the hypothalamus. With this recombinant model established, RNA-Seq was used to identify a large set of genes regulated by OX 1 signaling including several immediate early genes, transcription factors, kinases, and phosphatases. Results from qPCR validation experiments correlated well with the RNA-Seq data, affirming many of the transcriptional changes observed. In order to identify the primary transcriptional regulators or OX 1 signaling, an in silico transcription factor binding site analysis was performed. While the occurrence of false positives is a legitimate concern when using these types of analyses, correlating these results with the gene expression data further supports roles for several of the transcription factors identified in the promoter analysis. Notably, these included related transcription factors such as Egr1, Egr2, and Egr 3 as well as a number of Kruppel-like factors and the Myc transcription factor. Of particular importance, though, is whether or not these data are physiologically relevant. To address this, the data were compared to those from similar studies focused on orexinrelated behaviors, such as sleep and wakefulness. These comparisons demonstrated that OX 1 signaling shows similarities with transcriptional profiles seen in in vivo SD microarray studies, with several genes being similarly regulated. This suggests that orexin signaling may be responsible for the SD-induced changes in expression of these genes. Indeed, SD has been shown to result in the activation of orexigenic neurons and increased OxA peptide levels in the brain, further supporting this hypothesis.
One of the genes highly regulated by both orexin signaling in GT1-7-OX 1 cells and SD in the mouse brain is Sgk1, a kinase that is known to regulate a number of cellular proteins and is strongly associated with neuronal excitability, synaptic plasticity, and memory formation. As orexin also has strong associations with all of these processes, the role of Sgk1 in OX 1 signaling seemed to be particularly interesting. These analogous roles imply that Sgk1 may be a crucial mediating factor in transducing the biological effects of orexin. A likely scenario is that, when induced by OX 1 activation, Sgk1 phosphorylates cellular ion channels, altering the excitability of OX 1 -expressing neurons, thereby influencing associated behaviors.
In order to evaluate the downstream effects of Sgk1 on OX 1 signaling, we analyzed the impact of Sgk1 inhibition on the orexin-dependent differential regulation of a subset of genes by qPCR. Inhibition of Sgk1 resulted in diminished induction of a small set of genes that included Gja3, Ldlr, Dok7, Rara, Kcnk3, and Cldn4, amongst others ( Fig 5). Promoter analysis identified several putative transcription factors that could regulate these genes. Of these, Sp1 seems to be the most interesting, as it is known to be phosphorylated by Sgk1 and regulates transcription of many of the genes whose OX 1 -dependent transcription was impaired by Sgk1 inhibition, although the role of Sp1 was not empirically addressed in this study.
In summary, the orexin system sits at the crossroads of a diverse set of related behaviors including sleep, memory, synaptic plasticity, and addiction. To give insight into the molecular mechanisms influencing the role of orexin in regulating these behaviors, a neuronal, recombinant cell line was utilized to identify a set of candidate genes involved in OX 1 signaling. Corroborating evidence in the literature strongly supports the physiological relevance of this data set, as several of the genes regulated by OX 1 in the model cell line are similarly regulated by SD, in vivo. While this study brought focus to the role of Sgk1 in OX 1 signaling, the data set is rich with other candidate genes whose capacity in OX 1 signaling commands further study. Some examples include Homer1, Arc, Nr4a1, and Ldlr, each of which were heavily induced by OX 1 signaling and have been associated with numerous orexin-related behaviors [73,93,[106][107][108][109][110].  Table. Gene ID conversions. The complete list of genes that were regulated 2-fold or greater, at 3h or 8h, by OxA in GT1-7-OX 1 cells were entered into the DAVID Gene ID Conversion Tool to generate RefSeq mRNA IDs (See Methods). The RefSeq mRNA ID's were then entered into PSCAN for promoter analysis. (XLSX) S3 Table. List of primer pairs used for the Sgk1 inhibition qPCRs. The complete set of 96 primer pairs includes 89 of the most highly OX 1 -regulated genes plus 7 housekeeping genes (Actb, B2m, Gusb, Polr2a, Ppia, Rplp0, and Tbp).