Regulatory Architecture of the LβT2 Gonadotrope Cell Underlying the Response to Gonadotropin-Releasing Hormone

The LβT2 mouse pituitary cell line has many characteristics of a mature gonadotrope and is a widely used model system for studying the developmental processes and the response to gonadotropin-releasing hormone (GnRH). The global epigenetic landscape, which contributes to cell-specific gene regulatory mechanisms, and the single-cell transcriptome response variation of LβT2 cells have not been previously investigated. Here, we integrate the transcriptome and genome-wide chromatin accessibility state of LβT2 cells during GnRH stimulation. In addition, we examine cell-to-cell variability in the transcriptional response to GnRH using Gel bead-in-Emulsion Drop-seq technology. Analysis of a bulk RNA-seq data set obtained 45 min after exposure to either GnRH or vehicle identified 112 transcripts that were regulated >4-fold by GnRH (FDR < 0.05). The top regulated transcripts constitute, as determined by Bayesian massive public data integration analysis, a human pituitary-relevant coordinated gene program. Chromatin accessibility [assay for transposase-accessible chromatin with high-throughput sequencing (ATAC-seq)] data sets generated from GnRH-treated LβT2 cells identified more than 58,000 open chromatin regions, some containing notches consistent with bound transcription factor footprints. The study of the most prominent open regions showed that 75% were in transcriptionally active promoters or introns, supporting their involvement in active transcription. Lhb, Cga, and Egr1 showed significantly open chromatin over their promoters. While Fshb was closed over its promoter, several discrete significantly open regions were found at −40 to −90 kb, which may represent novel upstream enhancers. Chromatin accessibility determined by ATAC-seq was associated with high levels of gene expression determined by RNA-seq. We obtained high-quality single-cell Gel bead-in-Emulsion Drop-seq transcriptome data, with an average of >4,000 expressed genes/cell, from 1,992 vehicle- and 1,889 GnRH-treated cells. While the individual cell expression patterns showed high cell-to-cell variation, representing both biological and measurement variation, the average expression patterns correlated well with bulk RNA-seq data. Computational assignment of each cell to its precise cell cycle phase showed that the response to GnRH was unaffected by cell cycle. To our knowledge, this study represents the first genome-wide epigenetic and single-cell transcriptomic characterization of this important gonadotrope model. The data have been deposited publicly and should provide a resource for hypothesis generation and further study.

The LβT2 mouse pituitary cell line has many characteristics of a mature gonadotrope and is a widely used model system for studying the developmental processes and the response to gonadotropin-releasing hormone (GnRH). The global epigenetic landscape, which contributes to cell-specific gene regulatory mechanisms, and the single-cell transcriptome response variation of LβT2 cells have not been previously investigated. Here, we integrate the transcriptome and genome-wide chromatin accessibility state of LβT2 cells during GnRH stimulation. In addition, we examine cell-to-cell variability in the transcriptional response to GnRH using Gel bead-in-Emulsion Drop-seq technology. Analysis of a bulk RNA-seq data set obtained 45 min after exposure to either GnRH or vehicle identified 112 transcripts that were regulated >4-fold by GnRH (FDR < 0.05). The top regulated transcripts constitute, as determined by Bayesian massive public data integration analysis, a human pituitary-relevant coordinated gene program. Chromatin accessibility [assay for transposase-accessible chromatin with high-throughput sequencing (ATAC-seq)] data sets generated from GnRH-treated LβT2 cells identified more than 58,000 open chromatin regions, some containing notches consistent with bound transcription factor footprints. The study of the most prominent open regions showed that 75% were in transcriptionally active promoters or introns, supporting their involvement in active transcription. Lhb, Cga, and Egr1 showed significantly open chromatin over their promoters. While Fshb was closed over its promoter, several discrete significantly open regions were found at −40 to −90 kb, which may represent novel upstream enhancers. Chromatin accessibility determined by ATAC-seq was associated with high levels of gene expression determined by RNA-seq. We obtained high-quality single-cell Gel bead-in-Emulsion Drop-seq transcriptome data, with an average of >4,000 expressed genes/cell, from 1,992 vehicle-and 1,889 GnRH-treated cells. While the individual cell expression patterns showed high cell-to-cell variation, representing both biological and measurement variation, the average expression patterns correlated well with bulk RNA-seq data. Computational assignment of each cell to its precise inTrODUcTiOn Gonadotropin-releasing hormone (GnRH) plays a key role in the control of reproduction in mammals. Secreted by the hypothalamus in a pulsatile fashion, GnRH acts via its receptor (GnRHR) to trigger the synthesis and release of the luteinizing hormone (LH) and follicle-stimulating hormone (FSH) by the pituitary gonadotropes. In turn, the gonadotropins regulate gametogenesis and steroidogenesis in the gonads. The gonadotropins are composed of a common glycoprotein hormone α subunit (CGA) and a specific β subunit (LHβ or FSHβ). The frequency of GnRH pulse release varies at different stages of reproductive life, e.g., during puberty and the female menstrual cycle. GnRH pulse frequency differentially regulates gonadotropin subunit gene expression and gonadotropin secretion (1). While Lhb gene expression is preferentially induced by high-frequency GnRH pulses, low-frequency pulses favor Fshb expression (2,3).
The immortalized LβT2 gonadotrope cells have been used extensively as an in vitro model for the study of gonadotropin gene regulation and GnRH signaling. The cell line was developed through targeted tumorigenesis in mice carrying the rat LHβ regulatory region linked to the SV40 T-antigen oncogene (4)(5)(6). LβT2 cells have some functional characteristics of mature gonadotropes, as they express Cga, Gnrhr, and Lhb. The cell line responds to pulsatile GnRH stimulation by upregulating Lhb and Gnrhr and secreting LH. In the presence of steroid hormones, LβT2 cells further increase the LH secretory response to GnRH pulses as well as the levels of Lhb and Gnrhr mRNAs (6). In addition, LβT2 cells induce Fshb under either activin A (7,8) or GnRH pulse stimulation (3), with the level of Fshb being influenced by both pulse frequency and average concentration of GnRH (9). While LβT2 cells exhibit an increase in intracellular calcium and exocytosis in response to GnRH stimulation (5, 6), they differ from mature anterior pituitary cells in that they lack a characteristic large-amplitude calcium oscillatory response to GnRH (10). In addition, continuous GnRH stimulation does not induce Gnrhr gene expression, which is in contrast with rat pituitary cells (11).
Previous studies in LβT2 cells showed that GnRH activates a complex cell signaling network that rapidly induces the expression of early genes such as Egr1, c-Fos, and c-Jun (12)(13)(14), whose products consecutively activate the transcription of gonadotropin subunit genes. Over the past two decades, a number of studies in the LβT2 cell line have implicated various pituitary factors in gonadotropin subunit gene regulation. These factors include secreted peptides such as bone morphogenetic proteins, pituitary adenylate cyclaseactivating polypeptide, growth differentiation factor 9, VGF nerve growth factor inducible (15)(16)(17)(18)(19) [for review, see Ref. (20)], as well as transcription factors (TFs) such as AP1 (Fos/Jun heterodimer), SF1, and Egr1 (14,(21)(22)(23). Nevertheless, the molecular mechanisms underlying the gonadotrope response to GnRH and the decoding of the GnRH pulse signal are not fully understood.
Recent advances in high-throughput sequencing technologies have enabled researchers to solve key questions about gene regulation both at the chromatin and at the transcriptome levels. Hence, mapping of "open" chromatin regions using the assay for transposase-accessible chromatin with high-throughput sequencing (ATAC-seq) allows the detection of putative DNA regulatory regions that are likely bound by TFs (24,25). Correlating the transcriptome (measured by RNA-seq) with a map of open chromatin may identify transcriptional regulatory elements that are involved in the GnRH response. Furthermore, as distinct cells within a cell population display significant variations in RNA expression, analysis of cell-to-cell variability of gene expression can deepen our understanding of cell population complexity and transcriptome dynamics by isolating transcriptomic heterogeneity (e.g., cell cycle status) that is concealed in cell population studies [for reference, see Ref. (26)] and providing insight into the cellular variation in gene expression levels and induction. Gel bead-in-emulsion (GEM) Drop-seq is a droplet-based single-cell (SC) RNA-seq method that can profile thousands of individual cells per sample with high sensitivity (27)(28)(29).
An unbiased epigenomic and SC transcriptomic study of LβT2 gonadotrope cells has not been reported. Here, we generated a global mRNA expression profile and a genome-wide atlas of accessible chromatin in GnRH-stimulated LβT2 cells by analyzing RNA-seq data (30) and ATAC-seq data (Figure 1A), respectively. An integrative analysis of the transcriptome and open chromatin data identified key GnRH-regulated genes along with putative cis-regulatory elements. We also analyzed cell-tocell variability of the transcriptome of GnRH-stimulated LβT2 cells and the effects of cell cycle state on this response using GEM Drop-seq ( Figure 1B).

MaTerials anD MeThODs
cell culture and Treatment  (9)] were harvested by trypsinization and resuspended in medium, as indicated. Cell samples were lysed, generating crude nuclei preparations. Next, the transposition reaction was performed in the presence of Tn5 transposase, which cuts and ligates adapters at DNA regions of increased accessibility. The quality of the transposed DNA fragment libraries was assessed using a Bioanalyzer. Libraries were sequenced, mapped to the genome, and processed bioinformatically, and accessible genomic regions were detected (peak calling). (B) Single-cell (SC) RNA-seq [gel bead-in-emulsion (GEM) Drop-seq] workflow. Following treatment with GnRH or vehicle, cultured LβT2 cells were harvested by trypsinization and resuspended in a PBS BSA-containing buffer. Cell samples were loaded on a microfluidic chip, combined with reagents and barcoded gel beads to form Gel beads in Emulsion (GEMs), and then mixed with oil. The resulting emulsions were collected and reverse transcribed in SC droplets. Oil was removed, and barcoded cDNA was amplified, purified, and subject to quality control (QC) assessment. Amplified cDNA was incorporated into libraries, which were subsequently subject to QC evaluation using Kapa, Bioanalyzer, and Qubit. Libraries were pooled for sequencing and demultiplexed for subsequent analyses.
Cells were frozen in freezing medium containing 70% DMEM, 20% FBS, and 10% DMSO (Sigma) and maintained in liquid nitrogen. Cell line authentication was achieved by comparing our cells with an early passage aliquot of LβT2 cells provided by Dr. Mellon and used as a standard reference (Idexx Bioresearch, Columbia, MO, USA). Our results confirmed that our LβT2 cells were Mycoplasma free, were of mouse origin, and had similar markers as the original cell line aliquot.
Protocols for GnRH stimulation were selected to provide wellcharacterized response patterns. Hence, early gene responses in LβT2 cells are sensitive to GnRH concentration, while Fshb transcript levels respond to variations in both GnRH concentration and pulse patterns (9,13). Because GnRH stimulation of LβT2 cells was previously shown to induce an early gene program, whose transcripts are upregulated within an hour of GnRH exposure (12,13), a 40-to 45-min GnRH treatment was used in this study to analyze the transcriptome response to GnRH. For the analysis of chromatin accessibility by ATAC-seq, a pattern of pulsatile GnRH exposure was employed to capture the chromatin state that accompanies maximal Fshb induction (9, 18).

Bulk rna-seq assay
The generation of the RNA-seq data set that we analyze in this study was previously described (30). Briefly, LβT2 cells were serum starved overnight and stimulated with either 5 nm GnRH or vehicle for 45 min. Each group consisted of four independent replicates. Total RNA (2.5 µg) from each replicate was sequenced at the Mount Sinai Genomics Core Facility using an Illumina platform (Illumina, Inc., San Diego, CA, USA) and a HiSeq 2000 sequencing system (100-nucleotide length, single read type, multiplexing three samples per lane). The RNA-Seq data are deposited in GEO (GSE42120). assay for Transposase-accessible chromatin with high-Throughput sequencing Assay for transposase-accessible chromatin with high-through put sequencing was performed as previously described (25) on two replicate samples of LβT2 cells treated with 2 nM GnRH pulses every 2 h for a duration of 6 h and 45 min (4 pulses in total; cells harvested 45 min after last pulse) in a high-throughput GnRH pulse system (9). Briefly, 4,500 cells were washed with cold PBS at 4°C and lysed for 10 min at 4°C. Cell pellets were resuspended in the transposase reaction mix [5 µl 2× TD buffer, 0.5 µl transposase (Illumina) and 4.5 µl nuclease-free water] and incubated at 37°C for 30 min. DNA from the transposase reaction was purified with a DNA Clean & Concentrator-5 Kit (Zymo Research). PCR amplification was performed using Nextera PCR primers. The optimal number of cycles was determined via quantitative real-time PCR (qPCR) to stop the amplification before saturation. Libraries were purified with AMPure beads and then quantified using a KAPA Library Quantification Kit (Kapa Biosystems) and High-Sensitivity DNA Bionalyzer kit (Agilent) and sequenced on an Illumina HiSeq 2500 to >164M reads with 50 bp read length, paired-end. The ATAC-seq data are deposited in GEO (GSE102480).

geM Drop-seq assay
Gel bead-in-emulsion Drop-seq was performed as described [10× Genomics, Pleasanton, CA, USA (29)]. Briefly, LβT2 cells were seeded at 350,000 cells per well in 12-well plates for 48 h and treated on day 3 with either 2 nM GnRH or vehicle for 40 min. Cells were then trypsinized and resuspended in medium before being washed and resuspended in 1× PBS/0.04% BSA. Following filtration of the cell suspension, cells were counted on a Countess instrument, and viability was assessed to be above 90% using Trypan Blue. Final concentration was set at 1,000 cells/μl in 1× PBS/0.04% BSA. As a starting point, ~8,000 cells from each sample were loaded into the fluidics chip. Reverse transcription was performed in the emulsion, and cDNA was amplified for 12 cycles before library construction. Quality control (QC) and quantification of the amplified cDNA were assessed using the High-Sensitivity DNA Bioanalyzer kit. Library quality control and quantification were evaluated. The SC data set is deposited in GEO (GSE102480).

Quantification and Qc of rna and libraries
RNA concentrations were determined with Quant-iT RiboGreen RNA reagent (Invitrogen, Carlsbad, CA, USA) using a fluorescence microplate reader (SpectraMax M3, Molecular Devices, Sunnyvale, CA, USA). RNA quality was assessed by determining the RNA Integrity Number using Bioanalyzer.

Quantitative real-time Pcr
Following total RNA isolation, 1 µg of RNA was reverse transcribed with the Affinity Script reverse-transcriptase (Agilent, Santa Clara, CA, USA). Next, samples were diluted 1:20 in molecular biology grade H2O (Cellgro, Manassas, VA, USA). SYBR Green qPCR assays were performed (40 cycles) in an ABI Prism 7900HT thermal cycler (Applied Biosystems, Foster City, CA, USA) using 5 µl of cDNA template and 5 µl of master mix containing the specific primers for the targeted gene, Platinum ® Taq DNA polymerase, and the required qPCR buffer, following the manufacturer's recommendations. Three technical qPCR replicates were run for each biological replicate. Results were exported as cycle threshold (Ct) values, and Ct values of target genes were normalized to that of Rps11 in subsequent analysis. Data were expressed as arbitrary units by using the formula, E = 2,500 × 1.93( rps11 CT value − gene of interest CT value) , where E is the expression level in arbitrary units. Primer sequences were as previously described (9,12).

Bulk rna-seq Data analysis
The RNA-seq data generated about 36-45 million reads per sample. The RNA-seq reads were aligned using STAR (31) v2.5.1b with the mouse genome (GRCm38 assembly) and gene annotations (release M8, Ensembl version 83) downloaded from the https://www.gencodegenes.org/web site. 91-93% of the reads were uniquely mapped to the mouse transcriptome. The matrix counts of gene expression for all eight samples were computed by featureCounts v1.5.0-p1 (32). Differentially expressed genes (5% FDR and at least 2 log2 fold change) were identified using the voom method (33) in the Bioconductor (34) package Limma (35). When comparing the bulk RNA-seq analysis with SC RNA-seq or ATAC-seq data, the transcripts per million (TPM) computed by RSEM (36) was used for the comparison.

aTac-seq Data analysis
Primary data analysis involved read mapping via bowtie2 (37) to GRCm38 (mm10), followed by duplicate read removal and peak calling via MACS2 (38) with parameters "-g dm -nomodelshiftsize 100 extsize 200. " For examining correspondence between chromatin accessibility and gene expression (Figure 4A), all peaks were first annotated to the nearest gene, and the most significant peak score for each gene was selected.

sc rna-seq Data analysis
Single-cell RNA-seq data were processed using the Cell Ranger pipeline v1.3, which provides a data matrix of expression for all genes and all cells. Differentially expressed genes were analyzed using the sSeq method (41), as implemented in the R package cellrangerRkit v1.1. The cell phase computation for the single cells follows the ideas described in the Supplementary Material of the study by Macosko et al. (28) with our own customized R script implementation. A schematic of the cell phase score computation is described in Figure 6A. The t-SNE analysis (42,43) was performed using the implementation from the Cell Ranger pipeline.

Transcriptome Profiling of lβT2 cells
To characterize the transcriptome response to GnRH in LβT2 cells for comparison with the chromatin accessibility and SC studies described below, we analyzed an RNA-seq experiment in which LβT2 cells were exposed to either GnRH (5 nM) or vehicle for 45 min [(30); n = 4 for each group]. We identified 112 differentially expressed genes relative to the control (>4fold at FDR < 0.05, see Figure 2A; Figure S1 in Supplementary Material and Table S1 in Supplementary Material), including a large number of early gene transcripts that are known to be regulated by GnRH in gonadotropes, such as Egr1, Fos (c-Fos), and Jun (c-Jun) (12,45,46). As is evident from the asymmetrical expression change versus significance volcano plot (Figure 2B), the majority of regulated transcripts were upregulated at this time point. The plot highlights the most highly upregulated genes, which include Jun, Fos, Fosb, Egr1, Egr2, Egr3, Egr4, Nr4a3, and Nr4a1.
We investigated the functional relationships of the 21 highest fold-change annotated GnRH-regulated genes in the context of a human pituitary gene network using Genome-scale Integrated Analysis of gene Networks in Tissues [GIANT; Figure 2C; (39)]. GIANT uses a massive public data compendium to infer tissuespecific data-driven functional gene-gene relationships. The output of a GIANT analysis is a graph indicating the strength of the functional relationship of each pair of the input genes in the pituitary, as well as the inference of additional highly related genes. The edges connecting the nodes indicate the statistical strength of the evidence for that relationship in the human pituitary. Notably, the highest regulated genes form highly interrelated subgroups (e.g., Jun/Fos/Ier2), and many highly related inferred genes are regulated by GnRH (e.g., Dusp1, Nab2) or known to be central to key gonadotropin developmental or regulatory processes [e.g., Pitx2; Table S2 in Supplementary Material; (12,14,(47)(48)(49)(50)]. Included among the statistically significant gene enrichment set were MAPK signaling and SMAD protein signaling (Table S3 in Supplementary Material). Overall, this analysis reveals a coordinated gene program activated by short-term GnRH exposure.

genome-Wide Mapping of chromatin accessibility in lβT2 cells
To map open chromatin regions, we carried out ATAC-seq in LβT2 cells treated with GnRH ( Figure 1A). ATAC-seq uses the hyperactive Tn5 transposase, loaded with adapters for high-throughput DNA sequencing, to integrate into regions of accessible chromatin. The resulting DNA fragments, generated from locations of open chromatin, are amplified, sequenced, and computationally mapped to the genome to obtain a genome-wide accessibility landscape (24,25). The ATAC-seq libraries and the sequence data showed a characteristic ~200 bp size distribution periodicity (Figure 3A), reflecting individual nucleosome occupancy patterns and confirming specific transposase activity and assay accuracy. Primary data analysis (see Materials and Methods) identified more than 58,000 statistically significant regions of open chromatin (peaks). The open chromatin map was found to be reproducible across independent samples and libraries. We note that, while determining whether chromatin accessibility changes with GnRH exposure is an important question, the present study was intended only to provide a baseline analysis of open chromatin structure. Performing a reliable comparative analysis of changes in specific regions with GnRH stimulation would require assaying a large number of samples, which was not feasible for this investigation.
We examined the 2,000 most prominent open chromatin regions (showing the highest peak scores; see Materials and Methods) with respect to their location relative to annotated genomic features. Approximately 75% of these peaks were located in immediate gene promoters or introns, which is consistent with open chromatin in the proximal regions of transcriptionally active genes (Figure 3B). Focusing on the sequence fragments in genomic areas flanking the transcription start sites (TSS, −3 to 3 kb), we observed a distribution highly preferential to the regions in close vicinity to the TSS (Figure 3B). We next examined chromatin state in the proximal promoter regions of several key genes ( Figure 3C). While open chromatin peaks were detected at genes that are either constitutively expressed in LβT2 cells (Gnrhr) and/or regulated by GnRH (Lhb, Cga, and Egr1), chromatin was closed at the Tshb gene that is not expressed in gonadotropes. To gain insight into the overall transcriptional regulatory state of these cells, we used the HOMER tool to determine TF binding motifs showing enrichment among all open chromatin regions (40). The enrichment results show extremely high statistical significance (see Table S4 in Supplementary Material), reflecting the high power of this global analysis. Finding highly enriched global binding motifs for the TFs Smad2 and Six6, which have  been implicated in Fshb gene expression (51,52), suggests that this analysis is likely to have generated useful leads for further study. In these high-resolution ATAC-seq data, we were also able to observe TF footprints, which create a notch in a region of otherwise open chromatin due to the presence of a bound TF. As an example, Figure 3D shows a notch detected in the otherwise open Jun promoter that precisely matches an SP1 consensus site.
This finding is consistent with the known role of SP1 in mediating induction of this gene (53).
To examine the relationship of the epigenetic landscape and the global pattern of mRNA expression, we compared the ATAC-seq and RNA-seq data. While the GnRH treatment conditions used for ATAC-seq versus RNA-seq experiments differed, they were both compatible with early gene induction. Indeed,  (40) calculates a data-driven motif model from our ATAC-seq peak data, shown in forward and reverse configurations, and then associates this motif with known motifs. In this example, the known SP1 consensus motif is nearly identical to the data-driven motif.  previous analysis of temporal responses of the early genes Fos and Egr1 in response to GnRH pulse stimulation demonstrated that they are highly expressed within 40 min after the fourth pulse (9). In this analysis, we identified the open chromatin peak showing the highest ATAC-seq peak score associated with each annotated gene. This estimate of chromatin accessibility for each gene was plotted against its level of expression following GnRH treatment using the RNA-seq data ( Figure 4A). The results show a general pattern of increased chromatin accessibility for expressed genes. Genes found to be highly induced by GnRH (e.g., Egr1, Egr2, Egr4, Jun, Fosb, Atf3) tend to show high chromatin accessibility. The gonadotropin subunit genes Lhb and Fshb showed very high levels of chromatin accessibility at specific locations. The Lhb area of high accessibility, shown in Figure 3C, is centered around the proximal promoter. In contrast, Fshb does not show statistically significant open chromatin in the region of the proximal promoter. Notably, several significant open chromatin peaks were found between −40 and −90 kb upstream of the Fshb TSS. The segment of highest chromatin accessibility, at −90 kb, contains a binding motif for AP1, a known Fshb regulator [ Figure 4B; (14)]. This pattern raises the possibility that this region may represent a distal Fshb enhancer.
sc Transcriptome analysis of lβT2 cells and characterization of early gene response to gnrh To investigate cell-to-cell variability in gene expression in the gonadotrope response to GnRH, we performed a SC transcriptome analysis of LβT2 cells exposed to either GnRH or vehicle for 40 min using GEM Drop-seq (see Figure 1B). This assay measures the entire transcriptome in thousands of individual cells from each sample. Specifically, we sought to assess cell-tocell heterogeneity with respect to the individual cell response to GnRH.
We started with ~8,000 LβT2 cells in each sample. The resulting SC RNA-seq libraries all exhibited the expected electropherogram traces on the Agilent Bioanalyzer High Sensitivity Chip, thus passing the QC assessment. In addition, qPCR assays of several early response genes in the library provided further evidence that our SC RNA-seq libraries were suitable for the detection of individual gene expression and regulation (data not shown). We sequenced to a depth of ~300 million reads in the SC libraries from each sample (vehicle-and GnRH-treated LβT2 cells). We obtained good sequence data for 1,992 vehicle-treated and 1,889 GnRH-treated cells, with ~98,000 mean reads/cell and >4,000 median genes/cell above detection threshold. Analysis of the GEM Drop-seq data demonstrated even coverage, with the slow drop-off of expression/cell in the top nearly 2,000 cells being similar in untreated versus GnRH-treated cells (Figure 5A). These features are indicative of a high-quality SC transcriptome data set.
A heat map of selected genes expressed in all 3,881 analyzed cells revealed a global pattern of differential expression between GnRH-treated and vehicle-treated cells ( Figure 5B). We identified 95 differentially expressed genes that included known GnRH-regulated immediate-early genes such as Egr1, Fos, Fosb, Jun, Btg2, Junb, and Nr4a1 (Table S5 in Supplementary  Material). Notably, we observed high cell-to-cell heterogeneity in response to GnRH, and several GnRH-treated cells even exhibited a gene expression pattern similar to that of untreated cells (Figure 5B; Figure S2 in Supplementary Material). Comparison of SC gene expression measurements averaged across two randomly selected subgroups of 996 cells highlighted the high reproducibility and consistency of the data, as the correlation coefficient was >0.99 in either vehicle-or GnRH-treated cells (Figure 5C). To further evaluate the performance of the digital transcriptome, we compared Drop-seq gene expression measurements averaged across all SC with the bulk RNA-seq measurements analyzed above, in either vehicle-or GnRHtreated cells. Despite the SC and bulk sequencing data sets coming from different experiments, the aggregated Drop-seq data showed high correlation with the bulk RNA-seq data (correlation coefficient >0.90; Figure 5D). The high cell-to-cell variation in the response to GnRH and in gene expression levels results from a combination of technical measurement variation and true cell-to-cell biological variation. To gain a sense of the degree of cell expression measurement resulting from technical noise, we modeled the technical measurement variation by a Poisson distribution [(54); Figure 5E]. Notably, at high levels of expression, the contribution of technical noise is relatively small and supports the presence of high levels of true biological cell-to-cell expression variation in these cells.
analysis of cell cycle Dependence of the sc Transcriptome of gnrh-stimulated lβT2 cells One limitation in many studies of response to GnRH performed by our group and others in LβT2 cells is that the effects of cell division on GnRH-induced gene regulation are not controlled. LβT2 cultures comprise asynchronously dividing cells. A SC transcriptome experiment provides the possibility of identifying the cell cycle stage of each individual cell and determining whether cell cycle stage influences the response to GnRH. Previous studies in yeast indicate that cell cycle has global effects on protein and RNA synthesis, thus affecting the transcriptional activity (55)(56)(57). RNA levels are controlled by transcriptional bursting that can vary during different cell cycle stages (58).
We examined cell cycle state of LβT2 cells as well as the influence of cell cycle phase on gene expression. Comparison of each cell's average expression with gene sets known to be enriched in one of the five cycle phases [G1/S, S, G2/M, M, and M/G1; Table  S6 in Supplementary Material; for reference, see Ref. (44)] enabled alignment of each cell by the cell cycle stage (Figures 6A,B). Consistent with this alignment, a gene associated with the M cycle phase (Cenpf) exhibited the highest expression at the M phase and high expression at the G2/M phase ( Figure 6C). In contrast, Cenpf was poorly expressed at G1/S and moderately expressed at S and M/G1. Likewise, the expression levels of other genes known to be associated with a specific phase of the cell cycle (e.g., Pcna with G1/S and S, Top2a with both S and G2/M) were highest at that cell cycle phase ( Figure S3 in Supplementary Material). On the other hand, immediate-early genes induced by GnRH (e.g., Egr1, Fos, Junb) and other GnRH-regulated genes (e.g., Gdf9) showed no significant expression change with cell cycle phase ( Figure 6C). The overall cell distribution with respect to cell cycle phase was comparable in GnRH-treated cells versus control (vehicle-treated cells ) cells (data not shown). Additional analysis showed that the five cycle phases were represented throughout all analyzed cells, with individual cells partly forming clusters based on their cell cycle phase (Figure 6D).

DiscUssiOn
In this study, we analyzed the global transcriptional, epigenetic, and single-cell transcriptional landscapes of the LβT2 gonadotrope cell line using RNA-seq data, ATAC-seq data and GEM Drop-seq data, respectively. To our knowledge, this represents the first genome-wide epigenetic characterization and the first SC transcriptome study performed in any gonadotrope experimental system. Our results provide insight into the global transcriptional regulatory processes of these cells and provide data sets and hypotheses to guide further work in this field.
A recent characterization of epigenetics in gonadotrope models reported CpG methylation, DNAse hypersensitivity, and histone modification at regions of several specific genes in several cell lines, including LβT2 cells (59). Their finding of open chromatin correlating with expression at the specific genes investigated corresponds with our genome-wide characterization of this relationship. The global ATAC-seq epigenetic approach that we have pursued opens the avenue to entirely unexpected discovery, such as the putative new Fshb enhancer identified as highly open chromatin (see Figure 4B). Further study is needed to evaluate the functional role of this novel putative enhancer.
Analysis of SC transcriptome in vehicle-and GnRH-treated cells demonstrates, for the first time, that there is no influence of cell cycle stage on the gene response to GnRH. The SC variability observed at the level of basal gene expression and gene induction by GnRH is high and largely not explained by the technical   variation to which SC transcriptome analysis is prone (60). High levels of SC expression and gene induction variation in LβT2 cells could be anticipated based on previous studies of response variation and noise scale of selected transcripts (61,62). The cause of this high level of variation in expression and response to GnRH in LβT2 cells and whether it accurately models the expression and  Table S2 in Supplementary Material mapped to mouse gene using Homologene. The  gene sets for the five cell cycle phases (G1/2 response patterns of the intact mouse gonadotrope are unknown. While recapitulating many properties of a mature gonadotrope, LβT2 cells are a transformed cell line generated by tumorigenesis. It is conceivable that the heterogeneity of these cells has been augmented by the process of transformation (63,64). Cell-to-cell variation in expression and response to stimulation can also result from normal biological variation and stochastic mechanisms (65). Thus, alternatively, this line may faithfully reflect SC expression and response variation and represent an accurate model of the behavior of the primary gonadotrope. Individual gonadotropes largely function as a SC processor in controlling the reproductive axis and understanding the role of cell-to-cell variation in the engineering of this system is a relevant question to approach through SC biology. Evaluation of chromatin accessibility variation at the SC level [SC ATAC-seq (66,67)] is also an interesting area for future investigation. The generation of high-quality global ATAC-seq chromatin accessibility and SC GEM Drop-seq transcriptome data should provide a useful resource for the research community.
aUThOr cOnTriBUTiOns FR-Z and MF designed and performed research, analyzed, and interpreted data; YG, EZ, GN, ST, and HW contributed analytic tools and analyzed data; VN performed research and analyzed data; HP interpreted data and drafted the work; JT analyzed and interpreted data; SC conceived research, analyzed data, and drafted the work. All authors drafted or revised the work critically and approved the final version to be submitted.

acKnOWleDgMenTs
Sequencing for the ATAC-seq and GEM Drop-seq assays were performed at the Genomics and Epigenomics Core Facility of Weill Cornell Medical College. qPCR assays were performed at the Quantitative PCR CoRE at the ISMMS.