Transcriptome-scale spatial gene expression in rat arcuate nucleus during puberty

A variety of neurons in hypothalamus undergo a complicated regulation on transcription activity of multiple genes for hypothalamic–pituitary–gonadal axis activation during pubertal development. Identification of puberty-associated cell composition and characterization of the unique transcriptional signatures across different cells are beneficial to isolation of specific neurons and advanced understanding of their functions. The hypothalamus of female Sprague–Dawley rats in postnatal day-25, 35 and 45 were used to define the dynamic spatial atlas of gene expression in the arcuate nucleus (ARC) by 10× Genomics Visium platform. A surface protein expressed selectively by kisspeptin neurons was used to sort neurons by flow cytometric assay in vitro. The transcriptome of the isolated cells was examined using Smart sequencing. Four subclusters of neurons with similar gene expression signatures in ARC were identified. Only one subcluster showed the robust expression of Kiss1, which could be isolated by a unique membrane surface biomarker Solute carrier family 18 member A3 (SLC18A3). Moreover, genes in different subclusters presenting three expression modules distinctly functioned in each pubertal stage. Different types of cells representing distinct functions on glial or neuron differentiation, hormone secretion as well as estradiol response precisely affect and coordinate with each other, resulting in a complicated regulatory network for hypothalamic–pituitary–gonadal axis initiation and modulation. Our data revealed a comprehensive transcriptomic overview of ARC within different pubertal stages, which could serve as a valuable resource for the study of puberty and sexual development disorders.

Kiss1 neurons have been mapped within the arcuate nucleus (ARC) and anteroventral periventricular nucleus (AVPV) with dense continuum of cells as well as the dorsomedial nucleus and posterior hypothalamus with low-density group of scattered cells [4]. Besides that, the projection patterns of various Kiss1 neuron populations are also different. Kiss1 neurons at preoptic area project to the cell body and proximal dendrites of GnRH neurons and are responsible for the luteinizing hormone (LH) surge required for ovulation, while Kiss1 neurons in the ARC project to the distal dendron of GnRH neurons and generate the release of pulsatile GnRH [5]. Recently, researches have reported the extensive colocalization of neurokinin B (NKB) (encoded by TAC3) and dynorphin (encoded by Pdyn) with Kiss1 expression. Therefore, Kiss1 neurons were also identified as kisspeptin/neurokinin B/dynorphin (KNDy) neurons [6]. Additionally, the coordination of neurotransmitters secreted from other glutamatergic and GABAergic neurons also seems to be strongly interconnected and synchronised with Kiss1 neurons around puberty [7,8].
A fundamental step in understanding the neuroendocrine regulation of the reproductive axis requires accurate mapping of synaptic connectivity to KNDy neurons. Given the current research progress of the transgenic Kiss1 animal models, the transgenic GFP or β-galactosidase were too weak to be observed [9]. Moreover, the CRE-activated reporter usually drives the expression of transgenic target at the ectopic sites where kisspeptin protein is not regularly expressed [10]. The accurate characterization and isolation of KNDy neurons or other important neurons in vivo for molecular biological experiments (especially the next generation sequencing) was still missing.
Our previous studies have determined that the rats with postnatal day (PND)-25, 35 and 45 represent childhood, adolescent and adult stages via the fertile phenotypes as well as the dynamic changes of Kiss1 and GnRH expression [11,12]. Barcoding-based spatial transcriptomic sequencing [13] that provides the visualization of high-quality RNA-sequencing data in tissue with twodimensional positional information gives us a clue for advanced understanding of the distribution and dynamic changes of cell-specific neurons in wild type rodents during pubertal process without any transgenic symbols. In this study, we used 10× Genomics Visium platform to generate transcriptome-scale spatial gene expression in ARC regions of three individual female Sprague-Dawley (SD) rats with the different ages of PND-25, 35 and 45. By comparison among the different hypothalamic cell clusters, our data uncovered important novel hallmarks for specific Kiss1 neurons isolation. We provide these data as a significant resource for the neuroscience community to augment current molecular profiling and spatial transcriptomics efforts in the neuroendocrinology for pubertal development.

Experimental animals
Female Sprague-Dawley rats in PND-25, 35 and 45 purchased from Shanghai SLAC Laboratory Animal Co., Ltd. (Shanghai, China) were allowed free access to food and water. Rats were sacrificed via cervical dislocation. The whole brains and ovaries were isolated immediately. Ovaries were sectioned by 5 μm and conducted for hematoxylin-eosin (HE) staining (C0105S, Beyotime, Shanghai, China). The area of active corpus luteum was assessed for the sex maturation. Ten rats in each pubertal stage were collected for flow cytometric assay. All the procedures were followed by the Institutional Animal Care and Use Committee of Shanghai Children's Hospital.

Spatial transcriptomic sequencing
OCT-embedded whole brains were cryosectioned coronally to a thickness of 10 μm (bregma: − 2.52 to − 2.92, interaural: 6.08 to 6.48 mm) using a CryoStar NX70 cryostat (Thermo Fisher Scientific, Waltham, MA, USA). Tissue sections were layered onto the visium spatial tissue optimization slide containing oligonucleotides for mRNA capture (10× Genomics, Pleasanton, CA, USA). Each capture area per slide has 5000 spatially barcoded with a diameter of 55 μm and a center-to-center distance of 100 μm, over an area of 6.5 mm by 6.5 mm. A Master Mix containing reverse transcription (RT) reagents and fluorescently labeled nucleotides is added on top of the tissue sections, resulting in fluorescently labeled cDNA synthesis. Tissue is enzymatically removed, leaving behind fluorescent cDNA covalently linked to oligonucleotides on the slide. Fluorescent cDNA is visualized under fluorescence imaging conditions verified using the Visium Imaging Test Slide. Fluorescent cDNA is visualized under fluorescence imaging conditions verified using the Visium Imaging Test Slide. H&E and fluorescence images are compared. The permeabilization time that results in maximum fluorescence signal with the lowest signal diffusion is optimal. The libraries were sequenced with paired-end 150 bp sequencing (PE150) by NovaSeq 6000 platform.
Space ranger showed the capture area of the tissue in the slide and differentiates reads for each spot based on spatial barcode information. STAR was used to assess the sample quality by total number of spots, the number of pairs of reads in each Spot, the number of detected genes, and the number of unique molecular identifier (UMIs). Sctransform was used to normalize data, construct regularized negative binomial model of gene expression, and detect high variance characteristics. The method of t-distributed stochastic neighbor embedding (t-SNE) and uniform manifold approximation and projection for dimension reduction (UMAP) was used for dimension reduction. K-means was used to cluster the same type spots together. Seurat was used to identify the genes with spatial expression tendency through unsupervised clustering or prior knowledge. SPOTlight was used to identify the cell types based on deconvolution algorithm for nonnegative matrix factorization. MAST difference test method was used to evaluate differential expressed genes (DEGs) under adjusted p-values < 0.001. Gene ontology (GO) was used to analyze the functions and associated signaling pathways of DEGs.
Cells were analyzed on FACSCanto II (BD Bioscience, Franklin Lakes, NJ, USA) and screened on AriaIII (BD Bioscience). The following gating strategy was adopted to obtain SLC18A3 positive cells: first, FSC A/H gating was used to remove adherent cells, and then 7-AAD positive were removed to exclude dead cells. After that, FSA/ SSA gating was used to remove fragments, and two cell population were obtained, namely FSA large or small cell. Finally, four groups of cells were sorted to obtain SLC18A3-positive large cells, SLC18A3-positive small cells, SLC18A3-negative large cells, SLC18A3-negative small cells.

Smart sequencing
Cells harvested from FACS were centrifuged at 10,000g for sheath fluid removal, and suspended by 10  The raw data was trimmed adaptors and filter out low quality reads using Trimmomatic (non-default parameters: SLIDINGWINDOW:4:15 LEADING:10 TRAIL-ING:10 MINLEN:35), and checked the quality of clean reads using Fastqc. Next, clean reads were aligned to the latest mouse genome assembly mm10 using Hisat2 v2.0.5 (non-default parameters: -rna-strandness RF -dta). The transcripts were assembled and the expression levels were estimated with FPKM values using the String-Tie algorithm (non-default parameters: -rf ). Differential mRNA and lncRNA expression among the groups were evaluated using an R package Ballgown, and the significance of differences (log 2 |FC|> 1, q value < 0.05) by the Benjamini & Hochberg (BH) adjustment method were computed. Gene annotation was described by Ensembl genome browser database (http:// www. ensem bl. org/ index. html). The R package ClusterProfiler was used to annotate the differential genes with GO analysis.

Data deposits
The raw sequencing data was deposited to ArrayExpress assigned with the accession number E-MTAB-10975 and E-MTAB-10976.

Overview of hypothalamic cell clusters identified by spatial transcriptomic sequencing
Female SD rats aged PND-25, 35 and 45 were examined for the ovarian maturation. The area of active corpus luteum (≥ 30 mm) occupied in ovarian structure indicated the different reproductive ability during these pubertal stages (Additional file 1: Fig. S1). The whole brains were performed 10 μm serial tissue sections transversely to expose ARC regions (bregma: − 2.52 to − 2.92 mm, interaural: 6.08 to 6.48 mm) according to Allen Brain Atlas [15] and processed by spatial transcriptomic sequencing (Fig. 1A). The total 12,436 spots containing 10 cells average were sequenced to a median depth of 497.9 M reads with 5179 genes per spot and 29,277 UMIs (Additional file 2: Table S1). Quality control for the ratio of mitochondrial gene counts, expressed gene counts and unique molecular identifiers in each spot ensure the qualified sequencing data of tissues without cellular apoptosis or cytolysis (Additional file 1: Fig.  S2A-C).
Fourteen cell clusters in the detected regions of integrating three samples have been characterized by t-SNE (Fig. 1B) and UMAP (Fig. 1C) dimension reduction clustering analysis. We observed a series of clear-cut neurons in hypothalamus and adjacent regions distinguished by the global differentially expressed signatures in spots ( Fig. 1D-F), and the areas, structures and positions of all these neurons were generally similar in three samples. We observed that the clusters of ARC (Cluster 9), ventromedial hypothalamus (VMH) (Cluster 1) and paraventricular nucleus (PVN) (Cluster 13) gathered closely to each other compared to other clusters, indicating that the nucleus distributed in the adjacent regions displayed more similar gene expression profiles. Given the hallmarks of each cell type (Additional file 3: Table S2), we noticed that many specific marker genes such as Ghrh, Npvf and Ces1d in ARC as well as Trh, Dlk1 and Nkx2-1 in VMH were also co-expressed in PVN (Fig. 1G, Additional file 1: Fig. S3), which suggested that genes exhibited a gradient expression profiles across the transition zones among adjacent neuronal clusters. Moreover, these fourteen neuronal clusters could be further divided into GABAergic and glutamatergic neuronal cells, and nine non-neuronal cell types ( Fig. 1H-J, Additional file 1: Fig.  S4) by the combined analysis of the resources of single cell sequencing data of rat brain from "Allen Brain Atlases and Data" (https:// portal. brain-map. org/) [16]. From the perspective of the interplay of different cell types made up in each spot in spatial transcriptomic sequencing, the overall change of proportions across different pubertal stages was obvious. For example, the occurrence of glutamatergic and GABAergic neurons as well as GABAergic neurons and oligodendrocytes were gradually decreased from PND25 to PND45. In turn, the occurrence of glutamatergic neurons and astrocytes, glutamatergic neurons and oligodendrocytes, as well as astrocytes and oligodendrocytes were increased (Fig. 1K-M). Our initial glance indicated that ARC, PVN, VMH and other hypothalamic regions closely associated with pubertal development were all shown in our spatial transcriptomics data.

Identification and characterization of KNDy neurons in rat hypothalamus
Due to the crucial roles of KNDy neurons in puberty initiation, we majorly focused on the expressions of KNDy (Kiss1, Tac3 and Pdyn) in whole hypothalamic regions. The Kiss1 positive spots were concentrated on ARC, PVN and scattered in other regions, whereas Tac3 and Pdyn were robustly enriched not only in the neurons above, but also in VMH (Cluster 1), dorsomedial hypothalamic nucleus, dorsal part (DMD) (Cluster 4), posterior hypothalamic area (PH) (Cluster 5), medial tuberal nucleus (MTu) (Cluster 8), VMH ventrolateral part (VMHVL) (Cluster 10), ventral reuniens thalamic nucleus (VRe) (Cluster 11) and reuniens thalamic nucleus (Re) (Cluster 12) ( Fig. 2A, Additional file 1: Fig. S5). The widely acknowledged functions of NKB and dynorphin on mental controlling, pain repression and addiction indicated that these two peptides extensively existed in multiple types of neurons and were not the appropriate biomarkers for KNDy neurons tracing and isolation. Top ten biomarkers (Fezf1, Kiss1, Ghrh, Sox14, Npvf, S100g, Ces1d, Gck, Crhr2, Six6) of Cluster 9 distinguishing ARC from other thirteen clusters were provided to address whether there were any other hallmarks being used for KNDy neurons identification besides Kiss1 (Fig. 2B), and we observed that Sox14, S100g, Gck and Six6 were uniquely expressed in Cluster 9 more specific than Kiss1 (Additional file 1: Fig. S6).
Next, ARC containing abundant KNDy positive neurons was further divided into more subclusters according to the comprehensive gene expression. Four subtypes of neurons in ARC were shown (Fig. 3A). Herein, Kiss1 (FC = 1.417, p = 3.58 × 10 -24 ) and Tac3 (FC = 2.563, p = 5.02 × 10 -72 ) were highly expressed in the Subcluster 1 of Cluster 9 compared to other three subclusters, while Pdyn was highly expressed in Subcluster 2 (FC = 1.065, p = 5.12 × 10 -59 ) (Additional file 4: Table S3), indicating that population in Subcluster 1 was likely to contain a great deal of Kiss1-expressed neurons. KNDy neurons The thicker lines represent that these two cell types share more spots might be not equal with Kiss1-expressed neurons because of the differential localization of Pdyn expression. In Cluster 9, top ten genes of Tgfbi, Nr5a2, Spp1, Galp, Fndc3c1, Slc18a3, Dlx5, Glp1r, Tekt3, and Kiss1 could be considered as the appropriated candidate hallmarks to distinguish the largest number of Kiss1 positive neurons in Subcluster 1 from other three subclusters (Fig. 3B). Tgfbi, Nr5a2, Spp1 and Slc18a3 were picked up to examine their expression compared to Kiss1 in ARC by IF (Fig. 3C-G), the specificity of Nr5a2 and Slc18a3 in ARC beyond other regions was observed. In our data, PVN had too few spots to perform subclustering analysis, and although VMH could be also divided into four subclusters (Additional file 1: Fig. S7), the expression of Kiss1 (no statistical significance), Tac3 (FC = 1.269, p = 4.41 × 10 -57 in Subcluster 3) and Pdyn (FC = 1.3, p = 1.03 × 10 -68 in Subcluster 4) (Additional file 5: Table S4) implied that KNDy neurons were seemingly scattered but not centrally distributed in VMH. Taken together, the spatial transcriptomics sequencing data provided applicable biomarkers for KNDy neurons characterization from wild type rat hypothalamus.

Kiss1-expressed neuron isolation by Slc18a3
Kisspeptin, as a secreted protein, is not suitable for KNDy neuron isolation despite the best expression specificity. In top ten DEGs in Cluster 9, the cytomembrane surface expressed marker Slc18a3 was highly expressed in ARC, and central medial thalamic nucleus (Cluster 12 in Fig. 1C-E) (Fig. 4A) was used to examine the practicability of alive Kiss1-expressed neuron isolation from the hypothalamus by flow cytometric assay. Here, we obtained two populations according to the cell size by FSC/SSC signals, and further harvested a small proportion of Slc18a3 + cells from hypothalamus of PND-25, 35 and 45 (Fig. 4B). Slc18a3 + large cell counts seemed to gradually increase from PND-25 to PND-45 (Fig. 4B). To verify the effectiveness of this methods for isolation of living neurons with high expression of Kiss1, smartseq was conducted to detect the transcriptome of four population namely large Slc18a3 ± cells as well as small Slc18a3 ± cells. The transcriptome of large Slc18a3 + cells displayed highest expression of Kiss1, Tac3 and Pdyn compared to other three population (Fig. 4C, Additional file 6: Table S5). Concurrently, candidate genes including Nr5a2, Tgfbi and Slc18a3 were highly expressed in both large and small Slc18a3 + cells. Interestingly, Gfap, Aldh1l1, Sox10, and Olig2 were robustly expressed in small Slc18a3 + cells compared to large ones. Combined with the cell sizes, we speculated that the large ones were probably neurons while small ones were glial cells. As expected, GO analysis on DEGs compared between large and small Slc18a3 + cells showed that large population presented function on axon and dendrites development as well as hormone transport, whereas small population presented glial differentiation (Fig. 4D). Taken together, we provided a reliable method to isolate living Kiss1expressed neurons from hypothalamus of wild type rats.

Dynamic changes of gene expression in ARC in different stages of puberty
Although kisspeptin was tightly connected with puberty, the precise and clear dynamics of gene expression profile of Kiss1-expressed neurons in hypothalamus had never been described yet. The spatial transcriptomic sequencing of ARC at different pubertal stages provided an opportunity for understanding the transcriptional program of Kiss1-expressed neurons in vivo. To this end, an unsupervised pseudotime analysis was conducted in four subclusters of ARC. We noticed that the spots of ARC in PND-45 (Fig. 5A) were largely in accordance with the spots of Subcluster 1 (Fig. 5B), whereas the spots of ARC in PND-25 majorly overlapped with the spots of Subcluster 3, and the spots of ARC in PND-35 were involved in all subclusters, suggesting that the pseudotime axis of cell programming well mimicked the pubertal development (Fig. 5C, Additional file 1: Fig. S8). Nevertheless, it was noteworthy that the overall transcriptome of a proportion of ARC neurons in PND-35 were different from in PND-45 (Fig. 5C), indicating two parallel developmental routes of ARC neurons from PND-25 to PND-35 as well as from PND-25 to PND-45, but not a one-way route from PND-25, PND-35 to PND-45. In other words, the distinct subclusters between PND-35 and PND-45 might reflect the differences of genotype between puberty and complete sex maturation. Moreover, three modules with the distinct expression landscape of 2,948 hypervariable genes in each module were characterized along the puberty process based on the gene expression tendency (Fig. 5D, Additional file 7: Table S6). Given the related biological processes among different modules of gene expression by GO analysis, genes in Module 1 contributed to enhancing glial cells but repressing neuron cells proliferation induced by estradiol (Fig. 5E), and genes in Module 3 were responsible for neurons differentiation and signals transmission (Fig. 5G). Interestingly, only genes in Module 2 substantially played a regulatory role in hormone secretion (Fig. 5F). Correspondingly, we exemplified the expression of Pdyn and Fto (Module 3) as well as Kiss1, Tac3, Dlk1 and Pou2f2 (Module 2) in different modules by pseudotime analysis (Fig. 5H), indicating that multiple genes in ARC participated in pubertal onset and development actually via distinct functions.

Neuron type alteration across puberty
In ARC, our clustering analysis identified a number of GABAergic neuron subtypes, glutamatergic neuron subtypes, and non-neuronal subtypes, which were intercrossed with each other (Fig. 6A-C). In the annotated neuron types, GABAergic neurons with high expression of Meis2, glutamatergic neurons as well as astrocytes were largely coordinated with a variety of other cells throughout the entire route of sexual development. The components of these neuron mixture dynamically changed across the puberty (Fig. 6D-F, Additional file 1: Fig. S9). We found that GABAergic and glutamatergic neurons were both gradually replaced by non-neuron cells appearing in these four clusters from PND-25 to PND-45. Collectively, these results demonstrated that our unbiased data analysis were able to reveal the particular cell types in ARC that contributed to puberty initiation and development.

Discussion
The current study provides the spatial transcription atlas of multiple neuronal nuclei within ARC across different stages of puberty. Owing to the limitation of spots in our spatial transcriptome data, many specific neurons such as Pomc-, Ghrh-and Agrp-highly-expressed neurons are all classified as Subcluster 1, which displays a relatively rough description on the coordination of neurons in ARC compared to single-cell RNA-seq. Additionally, the data of AVPV as the other main group of Kiss1-expressed neuron bodies, median eminence as well as ARC in sexual development related disease model such as central precocious puberty is not included in this study. Nevertheless, compared to previous single-cell sequence studies on hypothalamus [17][18][19], our data uniquely extend the dynamic spatial identity of neurons nucleus subsets within ARC specifically for the event of female mammalian puberty. Our findings will be an abundant source for evaluation of transcriptional profiles in central precious puberty and sexual dysplasia, and provide comparative neuroendocrine data of for future studies. Previous studies on single-cell RNA-seq have determined that ARC can be divided into eleven clusters in embryonic E15 brain [20], and twenty clusters in adult brain [21] by gene expression signatures. The signaling of kisspeptin-GPR54 at GnRH neurons are critical for maintaining GnRH release needed for fertility especially in puberty. In general, kisspeptin-expressed neurons are majorly located in the rostral periventricular area of the third ventricle of rodents or the preoptic area of other mammals and the ARC of the hypothalamusare sites of the highest density of kisspeptin expression within the mammalian brain [6]. Due to the high degree of colocalization of the other two peptides (Neurokinin B and Dynorphin) with kisspeptin, and for simplicity, this population of kisspeptin-expressed neurons is also abbreviated as KNDy neurons. Although more detailed and comprehensive information of neurons in ARC have been characterized, the spatial position of KNDy neurons are not fully elucidated in ARC in pubertal event. Our data fills this gap to some extent, and provides the distribution and the alteration of these four clusters in ARC during pubertal development (Fig. 3A). From our results, we notice that the cell clusters with the highest expression of Pdyn are not identical with Kiss1 and Tac3 in ARC (Additional file 4: Table S3), suggesting that Pdyn may play a more extensive role in neuroendocrine beyond the reproductive function, and seem to be not a specific biomarker for KNDy neurons.
In turn, the top ten highly expressed genes in Subcluster 1 indicate a series of potential hallmarks for Kiss1 neurons. For example, transforming growth factor beta induced (Tgfbi) is never reported in puberty, and first discovered in Subcluster 1 of ARC. TGFBI exerted as a tumor suppressor or enhancer in cancer progression is an extracellular matrix (ECM) protein that is associated with a number of ECM proteins such as fibronectin, biglycan, decorin, and several types of collagen and functions as a ligand to modulate cell adhesion and migration via various types of integrins [22]. The high-expressed Tgfbi is supposed to be of benefit to the microenvironment for the signals communication among different neurons in ARC. Another important example is nuclear receptor subfamily 5 group A member 2 (Nr5a2), a critical regulator for adult neurogenesis [23]. NR5A2 is sufficient to reduce cell proliferation, govern fate specification and differentiation of adult neural stem cells, and promote axon outgrowth. Moreover, Nr5a2 can bind to the promoter of Kiss1 for stimulating Kiss1 transcription particularly in Kiss1 neurons of ARC but not AVPV [24] and regulates gonadotrope function [25], which suggests a molecular basis for the distinct regulation of basal kisspeptin expression between ARC and AVPV neurons.
In previous studies of the epigenetic mechanism on central puberty, specimens used for next generation sequencing or microarray are almost harvested from the entire ARC [26] and even hypothalamus tissue [27], which suggests an obscure landscape of gene expression profiles. Isolation of living Kiss1 neurons is necessary and benefit to culturing primary cell line, and carrying on more subsequent genetic experiments in vitro for puberty study. In our study, Slc18a3 is first discovered and examined the specificity of Kiss1 neurons by flow cytometric assay. Although Slc18a3 + cells only account for a very small proportion of the whole ARC tissues, this population can be further divided according to the cell size and density of intracellular particles. The cell number is gradually increased from PND-25 to 45 in the large group, while is relatively stable in the small group (Fig. 4B). The gene profiles by smart-seq indicate that the Slc18a3 + large population shows the higher expression of Kiss1, Nr5a2, Tgfbi, and Glp1r compared to negative population, which is closer to Subcluster 1 than the Slc18a3 + small population. Although the small population has been found the robust expression of Gfap, Aldh1l1, Sox10, and Olig2, which are the classic biomarkers of astrocyte, microglia and oligodendrocyte, the phenotype and function of Slc18a3 + cells still need to be further explored before we can truly define them. Nevertheless, we are able to isolate the Kiss1-expressed cells with high purity from wild type rat by our methods. Slc18a3 encodes for the vesicular acetylcholine transporter for loading newly synthesized acetylcholine from the neuronal cytoplasm into synaptic vesicles. It is a significant discovery that Kiss1 neurons have partial property of acetylcholinergic neurons, in which the role of Slc18a3 in pubertal regulation or other function needs to be figured out in future research.
Finally, it is significant to investigate the alteration of gene signatures from a single-cell perspective during puberty. Three different modules of gene expression in four subclusters in ARC summarize the dynamic changes at the specific pubertal stages (Fig. 5C, D), indicating a novel rule of transcriptional regulation in ARC. For example, Rbm38, Ptbp1, Mapkapk2, Ybx1, Samd4b, Carhsp1, Vim and Pabpc1 contributing to RNA stabilization and 5′ or 3′ UTR binding are all classified in Module 1, in which genes expression are increased by pseudotime analysis. The post-transcriptional regulation seems to be a new field to study puberty. These genes expressed in different types of cells representing distinct functions on glial or neuron differentiation, hormone secretion as well as estradiol response precisely coordinate and affect with each other, and result in a complicated regulatory network for hypothalamic-pituitary-gonadal axis initiation and modulation.

Conclusion
In summary, our data revealed a comprehensive transcriptomic atlas of ARC with different pubertal stages, which could serve as a valuable resource for puberty and sexual development disorders study.