All experiments were performed with permissions from the Animal Welfare Committee of the Institute of Experimental Medicine Hungarian Academy of Sciences (Permission Number: A5769-01) and in accordance with legal requirements of the European Community (Decree86/609/EEC). All animal experimentation described was conducted in accordance with accepted standards of humane animal care and all efforts were made to minimize suffering.

Adult, gonadally intact female mice were used from local colonies bred at the Medical Gene Technology Unit of the Institute of Experimental Medicine (IEM). They were housed in light (12:12 light-dark cycle, lights on at 06:00 h)—and temperature (22 ± 2°C) controlled environment, with free access to standard food and tap water. GnRH-green-fluorescent protein (GnRH-GFP) transgenic mice (Suter et al., 2000) bred on a C57BL/6J genetic background were used. In this animal model, a GnRH promoter segment drives selective GFP expression in the majority of GnRH neurons. The estrous cycle was monitored daily between 9 and 10 a.m. by microscopic evaluation of vaginal cytology (Nelson et al., 1982; Byers et al., 2012; Cora et al., 2015). Proestrous (n = 6) and metestrous (n = 6) female mice with at least two consecutive, regular estrous cycles were used. In order to avoid the possible circadian effect, animals were sacrificed at the same period of the day, between 16:00 and 18:00 h. Those animals were considered to be in the proestrus stage that fulfilled the following criteria: (1) vaginal smear staining with predominance of nucleated epithelial cells (Byers et al., 2012); (2) LH serum concentrations >5 mg/L (15.11 ± 3.4 mg/L); (3) uterus wet weights >0.15 g (0.19 ± 0.01 g). Accordingly, the following criteria were applied for the metestrous cycle phase: (1) vaginal smears consisting of the three cell types: leukocytes, cornified and nucleated epithelial cells (Byers et al., 2012); (2) serum LH levels < 0.5 mg/L (0.35 ± 0.02 mg/L); (3) uterus wet weights < 0.1 g (0.08 ± 0.01 g). Serum LH concentrations and uterus weight data are presented in Supplementary Figure 1.

Blood samples were collected from the heart of deeply anesthetized mice immediately before the brain fixation step. The samples were chilled on ice, centrifuged at 1300 g for 3 min at 4°C. Plasma was aspirated then frozen and stored at −80°C until further use. Serum LH concentrations were measured with a rodent LH ELISA kit #ERK R7010 from Endocrine Technologies Inc. (Newark, CA, USA) according to manufacturers' instructions.

Brain fixation, preparation of sections for the subsequent laser capture microdissection (LCM) and microarray profiling were performed as reported elsewhere (Khodosevich et al., 2007; Vastagh et al., 2015). Briefly, metestrous (n = 6) and proestrous female (n = 6) mice were deeply anesthetized and perfused transcardially with 80 ml 0.5% paraformaldehyde followed by 20% sucrose. For microdissection, 7 μm thick coronal brain sections were cut. Sections were mounted on PEN-membrane slides (Zeiss, Jena, Germany), processed further for laser microdissection. Uniform and representative sampling of the entire GnRH neuronal population was performed using LCM performed on a PALM Microbeam system (Carl Zeiss Microimaging Gmbh, Jena, Germany) which was equipped with an epifluorescent setup. About 250 GFP-positive neurons were dissected per animal from 80 to 100 consecutive sections to generate GnRH cell samples from each brain.

GnRH cell samples (metestrous: n = 6, proestrous: n = 6) collected with LCM were incubated in 200 ml lysis buffer at 56°C for 3 h. RNA was isolated from the lysate by proteinase K/acid phenol method (Khodosevich et al., 2007), then precipitated by adding isopropanol (Sigma-Aldrich) and 20 μg glycogen (Thermo Fischer Scientific, Waltham, MA, USA) at −20°C for 30 min followed by centrifugation at 14,000 g at 4°C. The pellet was washed in 70% ethanol, air-dried and resuspended in water. The remaining genomic DNA was eliminated by treatment the mixture with 1U of RNase-free DNase-I (Thermo Fischer Scientific). RNA was purified using RNeasy MinElute Cleanup kit (Qiagen, Hilden, Germany). Total RNA was eluted with 14 μl of ribonuclease-free water. For the analysis of RNA integrity, RNA was isolated from the tissue of the medial preoptic area, and measured using RNA Pico Chip on the 2100 Bioanalyzer (Agilent, Santa Clara, CA, US; Supplementary Figure 2A).

Library preparation and amplification were performed according to the manufacturer's instructions for the WTA2 kit (Sigma-Aldrich). In the first step of the WTA2 protocol, the RNA was reverse transcribed implementing non self-complementary and quasi-random 3′ and universal 5′ primers. The newly synthetized single strands were the templates for annealing and extension steps. When the SYBR Green signal reached a plateau, the reaction was stopped. The resultant library contained cDNA fragments with size distribution between 100 and 1000 base pairs (data obtained from Agilent Bioanalyzer using DNA 1000 microfluidic chip; Supplementary Figure 2B). The amplified double-stranded cDNA was purified and quantified on a Nanodrop ND-1000 spectrophotometer (Thermo Fisher Scientific).

Eight microgram of cDNA was fragmented by DNase I and biotinylated by terminal transferase obtained from the GeneChip Mapping 250K Nsp Assay Kit (Affymetrix Inc., Santa Clara, CA, USA). Hybridization, washing, staining, and scanning of Affymetrix Mouse Genome 430 PM Strip arrays were performed following the manufacturer's recommendations. The Mouse Genome 430 PM Strip array allows the analysis of 34,325 well-annotated genes using 45,123 distinct probe sets. Scanned images (DAT files) were transformed into intensities (CEL files) using the AGCC software (Affymetrix). RMA analysis was performed by means of the statistical analysis software Partek Genomics Suite (Partek Inc., St. Louis, MO, USA) to obtain probe set level expression estimates.

All statistical and data mining works were performed in R-environment (R Core Team, 2016) with Bioconductor packages (Huber et al., 2015). Quality assessment of microarrays (n = 12) was performed using affyQCReport (Kauffmann et al., 2009). For quality control data regarding RNA integrity, see Supplementary Figure 3. Raw microarray data were pre-processed for analysis using RMA (Robust Multi-Array Average (Irizarry et al., 2003). Fold change (FC) estimation and difference analysis of gene expression were based on linear models combined with Bayesian methods using limma package (Ritchie et al., 2015). FC was calculated from normalized and log₂ transformed gene expression microarray data for each probe sets. The obtained p-values were adjusted by the FDR-based method. The following cut-off criteria were applied on the differentially expressed genes (DEG): |FC| >1.6; and adjusted p-value (p_adj) < 0.05. The top 10 differentially expressed genes (up- and down-regulated in proestrus, respectively) are shown in Supplementary Table 1.

For over-representation analysis, the list of DEGs was analyzed using the web-based functional annotation tool DAVID Bioinformatics Resources 6.7 (https://david.ncifcrf.gov) at default settings (Huang da et al., 2009a,b).

For visualization of biological pathways, the data matrix of log₂ transformed and RMA-normalized microarray data were analyzed using a public web server Graphite Web (Sales et al., 2013). Of the available analytical methods, Signaling Pathway Impact Analysis (SPIA) was used. Over the classical probabilistic components, the SPIA also applies statistical systems biology approaches: it takes into consideration (1) the extent of the expression changes of each gene (2) the position of the DEG within a pathway, (3) the topology of the pathway, (4) the type of interaction between the genes to identify significantly impacted pathways where the total net accumulated perturbation in the pathway (tA) can be calculated (Tarca et al., 2009) The level of significance was 0.05 using FDR correction. Interacting genes of the significant neurotransmitter signaling pathways (pGFDR < 0.05) were visualized using Cytoscape (ver. 3.2.1) open source software platform.

For quantitative real-time PCR (qPCR) investigations of LCM-derived GnRH samples (proestrous females n = 6, metestrous females n = 5) RNA isolation and WTA were performed as described in the previous section. Amplified and column-purified cDNA was used as template for qPCR. Whole transcriptome-amplified cDNA from LCM samples were diluted in 0.1x TE buffer for qPCR investigation. Inventoried TaqMan assays were used to confirm microarray results by qPCR. Each assay consisted of a FAM dye-labeled TaqMan MGB probe and two PCR primers. Thermal cycling conditions of the qPCR were as follows: 2 min at 50°C and 10 min at 95°C, followed by 40 cycles of 15 s at 95°C and 1 min at 60° C using ViiA 7 real-time PCR platform (Thermo Fisher Scientific). Differential expression of genes was calculated by the 2^−ΔΔCt method (Livak and Schmittgen, 2001) using GAPDH gene as a reference. Student's t-test were used as a statistical method in comparison of gene expression of the two groups (proestrus: n = 6; metestrus: n = 5).

The functional annotation tool DAVID Bioinformatics Resources was used to discover the most relevant gene ontology (GO) annotations. GO categories were differentially represented (FDR adjusted p < 0.05) as follows: GO:0005874~microtubule, GO:0000166~nucleotide binding; GO:0016071~mRNA metabolic process; GO:0031988~membrane-bounded vesicle; GO:0007268~synaptic transmission; GO:0044456~synapse part.

Analysis of microarray data revealed differentially expressed genes associated with various neurotransmitter signaling mechanisms. Receptor subunits of GABAergic, glutamatergic, and cholinergic communication, as well as receptors participating in adrenergic, serotonergic, dopaminergic, glycinergic, adenosinergic, and purinergic neurotransmission exhibited significant changes (FC > 1.6; adjusted p < 0.05) in gene expression levels. In addition, genes encoding G-proteins and downstream effectors as well as neurotransmitter/amino acid transporters of the solute carrier family also showed differential expression (Figure 1 and Table 1).

Significant downregulation was found in gene expression of adenosinergic (Adora2a, Adora2b), adrenergic (Adra1b, Adra2a, Adra2c) serotonergic (Htr1b), purinergic (P2rx7), and cholinergic (Chrm4, Chrne) receptors or their subunits, respectively, by microarray data analysis. Similarly, the dopaminergic receptor Drd3 and Drd4 exhibited lower mRNA level in proestrus. The Gabbr1 coding for GABA-B receptor 1 and the GABA-A receptor subunit Gabra3, Gabrb2, Gabrb3, Gabrd, Gabrg2 were upregulated. The gene of the GABA-A receptor associated protein, Gabarap were also upregulated. Other GABA-A receptor subunit, the Gabrd showed downregulated gene expression. Receptor subunits of the glutamatergic signaling were upregulated (Grin3a, Gria1, Gria2, and Grin1) whereas Grin2d, the metabotropic glutamate receptor Grm1 and glutaminase (Gls) were downregulated. The expression level of the nicotinic cholinergic receptor beta polypeptide 2 (Chrnb2), as well as several heterotrimeric G-protein alpha subunits (Gnai1, Gnai2, and Gnas) were increased. Genes of downstream primary effector proteins were expressed at either higher (Adcy5) or lower level (Adcy3, Plcb2) in proestrus. The secondary effector Prkaca, Prkacb, and Prkca were upregulated. As far as the solute carrier families are concerned, the glutamate Slc1a4, the vesicular glutamate transporter Slc17a6 and sodium-dependent vesicular transporter Slc6a17 exhibited higher expression level in GnRH neurons in proestrus. Expression of Clcn6 and Glra1 genes, both involved in cellular chloride ion transport, were up- or down-regulated, respectively.

TaqMan real-time PCR was used for validation of microarray data (Figure 2A). Out of the 9 selected target genes, representing GABA, glutamate and acetylcholine neurotransmission, the differential expression of 8 genes was confirmed (Slc17a6, Grin1, Gria1, Chrnb2, Gabrb2, Gabbr1, Gabrb3, Gabra3) and 6 of them exhibited FC >1.5 (Slc17a6, Grin1, Gria1, Chrnb2, Gabrb2, and Gabbr1). We found significant differences in gene expression of GABAB1 and Slc17a6 genes by qPCR. Differences between other targets were not significant due to the high standard deviation of individual samples. The correlation coefficient (Pearson's r = −0.8445) indicates the quantitative interrelation between microarray and real-time PCR data (Figure 2B).

High throughput analysis of microarray data from 6 proestrous and 6 metestrous females by Signaling Pathway Impact Analysis (False Discovery Rate, FDR of DEGs < 0.05) identified 32 KEGG and 34 REACTOME pathways (adjusted pG using FDR, pGFdr < 0.05). Three of the KEGG pathways (GABAergic, cholinergic, and dopaminergic signaling) and two of the REACTOME pathways (Neurotransmitter receptor binding and downstream transmission in the postsynaptic cell; Activation of NMDA receptor upon glutamate binding and postsynaptic events) were strongly associated with neurotransmitter signaling and/or its postsynaptic events (Table 2).

In the previous decades, it was shown that hypothalamic GnRH neurons receive functional GABAergic primary afferent inputs (Leranth et al., 1985; Adler and Crowley, 1986; Herbison and Dyer, 1991; Kimura and Jinnai, 1994; Jung et al., 1998; Sim et al., 2000; DeFazio et al., 2002; Temple and Wray, 2005; Zhang et al., 2009). The estrous cycle is characterized by profound changes in GABA transmission to GnRH neurons (Herbison, 1997): prior the LH surge the mRNA level of the GABA synthetizing enzyme GAD67 is significantly reduced (Herbison et al., 1992) which exhibits E2-dependent temporal pattern in the territory of the AVPV (Curran-Rauhut and Petersen, 2002). The release of GABA also falls significantly prior to and during the time of estrogen-induced LH surges (Jarry et al., 1992; Tin-Tin-Win-Shwe et al., 2004). This drop of GABA level before the surge seems essential for the positive feedback (Herbison and Dyer, 1991; Kimura and Jinnai, 1994).

Expression of the α1-3, α5, β1-3, and γ2 subunits of the GABA_A receptor were detected in the GnRH neurons of adult female mice earlier (Sim et al., 2000). We found an increased mRNA expression level of the α1, β1-3, and γ2 subunits and in the proestrous female. Although a GnRH neuron-specific γ2 knock-out mice exhibited normal fertility (Lee et al., 2010), higher level of this subunit can increase the conductance and sensitivity of the GABA_A receptors (Brickley et al., 1999). The increased mRNA level of the GABA_B receptor 1 and that of several GABA_A subunits might also indicate the upregulation of these receptors in response of the dropping GABA level in proestrus. However, the observed negative value of the perturbation score (tA = −4.55) indicates a repressed status of the complete GABAergic synapse pathway in GnRH neurons, probably due to downregulation of genes (effector proteins i.e., Adcy4, Adcy3, Plcb2), downstream to the GABA receptors in the signaling pathway. Recent functional data support the view that net GABA effects might be determined by the balance of the excitatory GABA_A (Herbison and Moenter, 2011) and inhibitory GABA_B (Zhang et al., 2009; Liu et al., 2011) receptor-mediated tones on GnRH neurons, even in vivo (Constantin et al., 2013). We showed higher expression of both type of the GABA receptors and a repressed status of the GABAergic pathway suggesting an increased sensitivity for the ligand and inhibited signaling cascade that might play a role in the facilitation of the GnRH surge.

GnRH neurons express ionotropic AMPA (α-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid) and NMDA (N-methyl-D-aspartate) receptors (Gore et al., 1996; Spergel et al., 1999; Ottem et al., 2002) and their antagonists block the LH surge (López et al., 1990, 1992; Brann and Mahesh, 1991; Ping et al., 1997) indicating a prominent role of glutamate in the positive feedback regulation. At the time of the GnRH/LH surge, the number of contacts that express both glutamate and GABA transporters (Ottem et al., 2004), as well as the level of glutamate increases in the neighborhood of the GnRH cell bodies (Ping et al., 1994; Jarry et al., 1995). The expression of AMPA receptors changes at the level of the somata in steroid-induced LH surge model of young female rats, although, very few (1%) GnRH neurons co-synthetize AMPA receptors and Fos in the morning of proestrus. In contrast, the increase in the number of Gria1/Fos positive (25%) and Gria2/Fos positive (71%) GnRH neurons coincides with the increase and peak of LH (Bailey et al., 2006). In concert with these data, we found that expression of the AMPA1 and AMPA2 receptor subunits (Gria1 and Gria2, respectively) increased by the time of LH surge. Though the percentage of NMDAR1-immunopositive GnRH neurons does not change in the mouse (Adjan et al., 2008), we found that the expression level of NMDA subunits Grin1 and Grin3a increased significantly in proestrus. The higher expression level of these receptor subunits and enrichment of genes in the ionotropic glutamate receptor complex suggest an enhanced glutamate receptor signaling in GnRH neurons during the preovulatory GnRH/LH surge.

In the past decades, acetylcholine (ACh) has been shown to exert a stimulatory effect on the release of the GnRH peptide in vitro (Fiorindo and Martini, 1975; Richardson et al., 1982). Cholinergic afferents to hypothalamic GnRH neurons have been detected at both light and electron microscopic levels in the rat (Turi et al., 2008). Studies using perifused hypothalamic and immortalized GnRH (GT1-7) neurons indicated that selective activation of the nicotinic ACh receptor stimulated, in contrast, muscarinic receptor-specific activation inhibited GnRH release (Krsmanovic et al., 1998). In vivo pharmacological experiments in rats suggested an estrous cycle-dependent stimulation of GnRH release by selective muscarinic antagonists (Koren et al., 1992) induced via the M4 receptor subtype encoded by the Chrm4 gene. We found a significantly downregulated Chrm4 expression in proestrus and an upregulated level of Chrnb2. The altered regulation of the two subunits and the diverse action of the nicotinic and muscarinic-type ACh receptor activation can promote GnRH release and surge in proestrus.

GnRH neurons receive tyrosine-hydroxylase (TH) positive terminals demonstrated both at light microscopic (Jennes et al., 1983) and ultrastructural levels (Leranth et al., 1988; Chen et al., 1989; Horvath et al., 1993). Dual phenotype kisspeptin-TH positive fibers of AVPV origin apposed to the somata of GnRH neurons can be one of the sources of these dopaminergic inputs (Clarkson and Herbison, 2011). Dopamine inhibited the firing and AVPV-evoked GABA/glutamate postsynaptic currents in ~50% of the GnRH neurons in vitro mediated by D1 and D2-like receptors in male and female mice (Liu and Herbison, 2013) showing no evidence of an estrous cycle-dependent modulation when diestrous, proestrous, and estrous female mice were compared in the morning period. In the current study, D2-like dopamine receptors encoded by Drd3 and Drd4 genes were found to be down-regulated in proestrous GnRH neurons, suggesting that expression of these genes are affected by the estrous cycle. However, further functional studies are needed to clarify if down-regulation in the expression of the D3 and D4 receptor genes can decrease the inhibitory effect of dopamine on the excitability of the GnRH neurons in the late afternoon period of proestrus.

Numerous studies have supported the view for a long time that norepinephrine (NE) exerts a stimulatory effect on GnRH release in the presence of E2 in OVX+E replacement models (Gallo and Drouva, 1979; Ferris et al., 1984; Wise, 1984). Later, NE was found to suppress the excitability of the GnRH neurons acting through both alpha 1 and beta adrenergic receptors with reduced rate of responsive GnRH neurons during proestrus, otherwise exerting a similar effect across the estrous cycle (Han and Herbison, 2008). We found that adrenergic receptors Adra1b, Adra2a, and Adra2c are significantly downregulated in proestrus compared to metestrus. The decreased mRNA expression level of these receptors might indicate an attenuated catecholamine sensitivity of GnRH neurons specifically in the late afternoon period; prior to and during LH surge.

Earlier investigations provided evidence that 5-HT neurons project directly to GnRH neurons in rodents (Jennes et al., 1982; Kiss and Halász, 1985; Wada et al., 2006; Campbell and Herbison, 2007) In the mouse, serotonin (5-HT) exerts a biphasic effect on the excitability of GnRH neurons: both inhibitory and excitatory responses have been demonstrated, mediated by activation of 5-HT1A and 5-HT2A receptors, respectively (Bhattarai et al., 2014). The serotonergic input to GnRH neurons showed an estrous cycle-dependent regulation in the adult female mouse: during proestrus, the serotonergic inhibition was reduced, accompanied by the absence of the biphasic responses. The 5-HT1B receptor also mediates inhibitory neurotransmission (Mathur et al., 2011; Huang et al., 2013) in the CNS. Although the function of the 5-HT1B receptor was not investigated in GnRH neurons to date, we found that it was down-regulated in proestrous GnRH neurons. Compared to other stages of the estrous cycle, 5-HT1 receptor expression is decreased up to 40% in the basal forebrain during proestrus (Biegon et al., 1980) when peaking level of E2 significantly down-regulates 5-HT receptors (Biegon and McEwen, 1982). The current study further supports that a suppression of the inhibitory 5-HT tone takes place in GnRH neurons in proestrus.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.