Skip to main content
Download PDF

Postnatal developmental trajectory of sex-biased gene expression in the mouse pituitary gland

Biology of Sex Differences volume 13, Article number: 57 (2022) Cite this article

Abstract

Background

The pituitary gland regulates essential physiological processes such as growth, pubertal onset, stress response, metabolism, reproduction, and lactation. While sex biases in these functions and hormone production have been described, the underlying identity, temporal deployment, and cell-type specificity of sex-biased pituitary gene regulatory networks are not fully understood.

Methods

To capture sex differences in pituitary gene regulation dynamics during postnatal development, we performed 3’ untranslated region sequencing and small RNA sequencing to ascertain gene and microRNA expression, respectively, across five postnatal ages (postnatal days 12, 22, 27, 32, 37) that span the pubertal transition in female and male C57BL/6J mouse pituitaries (n = 5–6 biological replicates for each sex at each age).

Results

We observed over 900 instances of sex-biased gene expression and 17 sex-biased microRNAs, with the majority of sex differences occurring with puberty. Using miRNA–gene target interaction databases, we identified 18 sex-biased genes that were putative targets of 5 sex-biased microRNAs. In addition, by combining our bulk RNA-seq with publicly available male and female mouse pituitary single-nuclei RNA-seq data, we obtained evidence that cell-type proportion sex differences exist prior to puberty and persist post-puberty for three major hormone-producing cell types: somatotropes, lactotropes, and gonadotropes. Finally, we identified sex-biased genes in these three pituitary cell types after accounting for cell-type proportion differences between sexes.

Conclusion

Our study reveals the identity and postnatal developmental trajectory of sex-biased gene expression in the mouse pituitary. This work also highlights the importance of considering sex biases in cell-type composition when understanding sex differences in the processes regulated by the pituitary gland.

Highlights

  • Male and female mouse pituitary gland gene and miRNA expression was profiled across five postnatal ages spanning pubertal development.

  • Sex differences in pituitary gene expression exist prior to puberty and become more prominent upon puberty.

  • Combining expression data from genes and miRNAs revealed 18 putative sex-biased gene targets of 5 sex-biased miRNAs.

  • Sex differences in the proportions of somatotropes, lactotropes, and gonadotropes are predicted to occur prior to puberty.

Background

The pituitary gland plays a central role in regulating growth, lactation, reproduction, metabolism, stress responses, and puberty. These physiological processes are mediated by hormones released from five main anterior pituitary cell-types: growth hormone (GH) from somatotropes, prolactin (PRL) from lactotropes, follicle-stimulating hormone (FSH) and luteinizing hormone (LH) from gonadotropes, thyroid-stimulating hormone (TSH) from thyrotropes, and adrenocorticotrophic hormone (ACTH) from corticotropes [1]. Unlike the more glandular anterior pituitary, the more neural posterior pituitary consists primarily of astroglial-like pituicytes as well as axonal projections from the hypothalamus which store and release oxytocin and vasopressin into the systemic circulation [2, 3]. Non-hormone producing pituitary cells, including stem cells, and folliculostellate cells, are also present in the pituitary gland and support hormone-producing cells by functioning as progenitor cells and facilitating intercellular signaling within the pituitary [4, 5].

Both pituitary hormone production and many physiological processes regulated by the pituitary gland are sex-biased. For example, GH is secreted in a sexually different pattern in rodents—more pulsatile in males compared to females, and regulates sex-biased gene expression in liver [6]. Moreover, the clinical presentation and prevalence of pituitary-related disorders can also differ between sexes. For example, the prevalence of prolactinoma is significantly higher in women, but men are more likely to present with macroadenomas [7, 8]. While the sex differences in pituitary function and disease are well known, the gene regulatory networks underlying such differences remain elusive.

Several studies in rodents and humans have clearly highlighted sex differences in pituitary gland gene regulation. These studies include: targeted qPCR profiling of genes encoding for the main pituitary hormones in rat anterior pituitaries [9], serial analysis of gene expression (SAGE) in whole adult mouse pituitaries [10], and RNA-sequencing of adult human pituitaries as part of the Genotype-Tissue Expression (GTEx) project [11,12,13]. Most recently single-cell RNA-seq (scRNA-seq) has been performed in male and female adult mouse and rat pituitary glands revealing genes with sex-biased expression within specific cell types [14,15,16]. While most gene expression studies have focused on the adult pituitary, sex differences in pre-pubertal gene expression have been revealed using qPCR in mouse and rat pituitary glands [9, 17] and by RNA-seq in juvenile mouse gonadotropes [18]. While these studies suggest that some sex differences in pituitary gene regulation are established prior to puberty, we still lack a comprehensive view of pituitary gene regulation during postnatal development.

Another essential aspect of gene regulation is post-transcriptional regulation by microRNAs (miRNAs). While miRNA expression has been explored in pituitary glands of several mammalian species, these studies were not focused on postnatal development in males and females [19,20,21,22,23,24]. Clear examples of miRNA regulation of sex-biased gene expression have been reported in the neonatal hypothalamus and pubertal liver [25, 26], but evidence is currently limited in the pituitary gland. In one example, gonadotrope-specific Dicer knock-out (KO) in the mouse results in male-specific loss of follicle-stimulating hormone β-subunit (Fshb) gene expression, suggesting a sex-specific DICER-dependent post-transcriptional regulation of Fshb. This example points to the need for more comprehensive study of miRNA regulation of sex-biased gene expression in the whole pituitary.

The objective of this study was to characterize male and female pituitary gene expression during postnatal development and investigate the role of miRNAs in regulating pituitary sex differences. To achieve this, we profiled gene (3’ untranslated region sequencing, 3’UTR-seq) and miRNA (small RNA sequencing, sRNA-seq) expression in male and female mice at multiple postnatal days spanning pubertal transition to identify genes and miRNAs exhibiting known or novel sex differences. The resulting temporal gene and miRNA expression data from this study can be queried and visualized at https://wilsonlab-sickkids-uoft.shinyapps.io/pituitary_gene_mirna_shiny/. By combining this data with published single-cell RNA-seq (scRNA-seq) datasets, we provide evidence that sex differences in cell-type proportions emerge prior to the onset of puberty and likely contribute to sex biases in bulk gene expression.

Methods

Animal and tissue collection

All studies and procedures were approved by the Toronto Centre for Phenogenomics (TCP) Animal Care Committee (AUP 09-08-0097) (see Ethics approval for details). Conditions in which C57BL/6J mice were maintained, killed, and dissected are described previously in [17] as pituitary samples from the same animals were used in the current study.

Physical markers of puberty, vaginal opening (VO) and preputial separation (PS), were assessed every day after weaning to determine the pubertal stage of our female and male mice, respectively [27,28,29]. For assessing VO, the mouse was held by her base and a sterile pipette tip was used to brush the fur and assess the opening of the vagina. For assessing PS, the mouse was held on the hopper by his base and the degree of separation was assessed by gently pushing the prepuce using a sterile pipette tip [28].

Upon dissection, the pituitary gland was directly moved to RNAlater (containing 10% w/v sodium citrate tribasic dihydrate and 60% w/v ammonium sulphate) and stored at − 20 °C until RNA extraction.

RNA extraction

Prior to RNA extraction, pituitary tissue samples were placed into bead mill tubes containing six 1.4-mm ceramic beads (MoBio Laboratories) and homogenized for 30 s at 6.5 m/s at 4 °C using an Omni Bead Ruptor 24 bead mill. RNA was extracted with the NucleoSpin® miRNA kit (Macherey Nagel) in combination with TRIzol lysis (Invitrogen) following the manufacturers’ protocols, allowing for collection of small RNA (< 200 nt) and large RNA (> 200 nt) simultaneously into separate tubes from total RNA. RNA quantity was determined using Nanodrop and the quality was assessed by an Agilent 2100 Bioanalyzer.

Library preparation and sequencing

To construct RNA-seq libraries, we established an automated 3’UTR-seq (QuantSeq 3’mRNA-seq; Lexogen GmbH, Vienna) using the Agilent NGS Workstation (Agilent Technologies, Santa Clara) at The Centre for Applied Genomics (TCAG) (Toronto, Canada) as per the manufacturer’s protocol [30]. Briefly, 250 ng of RNA, from the large RNA fraction, was used to generate cDNA. cDNA was amplified with 17 PCR cycles as determined by qPCR analysis using the PCR Add-on kit (Lexogen). ERCC RNA spike-in Mix 1 was added following the manufacturer’s instructions. The resulting libraries were quantified with Qubit DNA HS (Thermo Fisher, Waltham). Libraries with insufficient library quantity determined by Qubit (below detection range) were excluded from further downstream processing. The remaining libraries had fragment sizes analyzed on the Agilent Bioanalyzer using the High Sensitivity DNA assay prior to sequencing. Sequencing was performed at TCAG on the HiSeq 2500 v4 flow cell (Illumina, San Diego) with SR50 bp with cycles extended to 68 bp.

Small RNA-seq libraries were generated from the small RNA fraction of the same samples (except for PD37 which was excluded from this experiment) by TCAG using NEBNext Small RNA Library Kit (New England BioLabs) according to the manufacturer’s protocol in two batches: 20 ng of RNA was used for batch 1 (replicates 1–3) and 10 ng of RNA was used for batch 2 (replicates 4–6). Ages and sexes were equally represented in the two batches. Batch effects were corrected for using RUVSeq (see Methods—mRNA and miRNA normalization and differential analysis). Sequencing was performed at TCAG on a HiSeq 2500 v4 flow cell with SR50 bp.

mRNA sequencing reads processing

FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) was used to examine the quality of sequenced reads. Next, a customized script was used to trim both the polyAs and adapter sequences at the end of the reads. A subset of sequencing reads obtained from 3’UTR-seq show a mixture of polyAs and sequencing adapters towards the end of the reads, which are not effectively trimmed by available read-trimming tools. We developed a trimming strategy which can identify and trim off polyA sequences embedded in the adapter sequences. Only reads longer than 36 bp after trimming were used. In addition, the first 12 nucleotides were trimmed based on the manufacturer’s recommendations. After trimming, FastQC was performed again to examine read quality. At this step, we found that ribosomal reads were overrepresented through priming by oligo(dT) binding to A-rich regions in the ribosomal RNA loci. In addition, the sequence of a brain-enriched small RNA, BC1, was also overrepresented. Thus, reads that map to these overrepresented transcripts were removed. Trimmed and filtered reads were aligned to the mouse genome using a splice-aware aligner, STAR (version 2.7.0f) [31], with parameters “–outWigType wiggle –twopassMode Basic –outFilterMismatchNoverLmax 0.05”. To build the genome index used by STAR, mouse genome was obtained from the UCSC database (mm10) and combined with ERCC sequences while mouse gene annotation was obtained from GENCODE (VM21) and combined with ERCC annotations. ERCC sequences and annotations were downloaded from the manufacturer’s website (https://assets.thermofisher.com/TFSAssets/LSG/manuals/ERCC92.zip). Quality control of mapped RNA-seq reads was performed using Qualimap (version 2.2.1) [32]. To visualize QuantSeq signal on the UCSC genome browser, wiggle files generated by STAR were converted to bigWig file format using “wigToBigWig” obtained from UCSC after removing ERCC genome coordinates.

PolyA site identification and gene annotation modification

GENCODE version M21 was the primary annotation used. To achieve a more comprehensive annotation of the 3’UTRs, we also incorporated the 3’UTRs annotated in RefSeq, which is obtained from the UCSC database (mm10). In addition, we identified potential polyA (pA) sites from the data. To do this, only the 3’ most nucleotide of each read is used to build a signal track for each sample. R package “derfinder” [33] was used to identify expressed regions (ERs) from these signal tracks. Specifically, an average read pile-up cutoff of 1 RPM (reads per million mapped reads) was used. ERs are annotated to gene annotations (3’UTR, 5’ UTR, exons, introns, and intergenic regions) based on GENCODE version M21, allowing for overlapping categories. ERs mapped to introns and intergenic regions are further analyzed to identify novel polyA sites. To filter for potential internal polyA priming events, sequence composition around ERs is examined. ERs with (a) matches 18-mer polyAs (with up to 6 mismatches) within 150 bp downstream from the ends; (b) matches 7-mer polyAs (with up to 1 mismatches) within 20 bp downstream from the ends; or (c) more than 50% of the polyAs within 20 bp downstream from the ends, are removed. In addition, ERs that overlap more than 20 bp with an annotated repeat region are also excluded. Filtered ERs are first mapped to RefSeq 3’UTR annotations (obtained from UCSC) and are associated with the corresponding genes. The rest of the unmapped ERs are then annotated to (a) the corresponding gene if it is intronic, or (b) the nearest gene upstream if it is within 5 kb from the gene ends. The intronic ERs are extended for 5 bp in each direction and the ERs downstream of genes are used to extend the gene’s 3’UTR annotation. The novel intronic polyA sites and extended gene annotations are then added to the gtf file used for gene counting. In total, additional internal polyA sites were added to 228 genes, and extended 3’UTRs were added to 476 genes, and 28 genes have both. Given the complexity of transcripts, the assignment to intronic polyA sites or 3’UTR extension may be impossible to distinguish in some cases. In addition, 2 of the 676 novel ERs did not make a difference to gene quantification as they overlapped with annotations from other genes.

Gene quantification

Trimmed and filtered reads were assigned to genes using featureCounts (v1.6.2) [34] with parameters “ -s 1 -Q 255” for 3’UTR-seq.

Processing of small RNA-sequencing reads

FastQC was used to examine the quality of sequenced reads. BBDuk (BBMap suite v37.90) was used to trim adapter sequences from reads with reference adapter sequences provided by BBMap suite and settings “hdist = 1 mink = 11” for small RNA-seq reads [35]. For miRNA size specificity, only reads less than 23 nucleotides in length were retained. Following trimming, FastQC was used to examine the quality of trimmed sequenced reads. miRDeep2 mapper.pl was used with default parameters to map reads of at least 18 nucleotides in length to the mouse genome (mm10) [36]. Known and novel miRNAs were identified using miRDeep2 main algorithm (miRDeep2.pl) with default parameters. For known miRNAs, the mature miRNA sequences in mouse were obtained from miRBase (v21) [37]. For novel miRNAs, only those with miRDeep score ≥ 2 and a sequence not matching previously reported small RNAs (rfam alert = FALSE) were retained for downstream analysis.

mRNA and miRNA normalization and differential analysis

Low-count mRNAs and miRNAs were filtered out prior to analysis. Only mRNAs and miRNAs with normalized read count (counts per million mapped reads, CPM) > 2 in at least 10 samples were retained for downstream analysis. CPM was used because of the consistent read mapping with UTR-seq.

For mRNAs, quantification of mitochondrial genes was not considered in this study and ERCC transcripts were removed prior to differential analysis. Remove Unwanted Variation from RNA-Seq Data (RUVSeq, v1.18.0) was used with RUVg() function with empirically detected negative genes to estimate unwanted variations in mRNA data based on the previously shown superior performance of this method compared to methods using ERCC for library normalization [38]. Empirical negative-control genes were identified with an ANOVA-like test comparing all conditions (FDR < 0.1). For miRNAs, RUVSeq was used with replicates (RUVs) to normalize and remove variation between the two batches from miRNA counts [38].

Quasi-likelihood F-test method was used to test for differential expression of mRNAs and miRNAs with a significance cutoff of absolute fold-change (FC) > 1.5 and false discovery rate corrected (FDR) < 0.05 using edgeR (v3.26.5) [39, 40]. Both the mRNA and miRNA analyses were performed in R version 3.6.

Novel miRNA identification

The mature sequences of novel miRNAs which were included in the miRNA differential expression analysis were used to identify homologous miRNAs in other species reported in miRBase v21 and precursor genome coordinates of the novel miRNAs were used to gain insight into their mechanism of biogenesis [41,42,43]. Specifically, miRNA identification with “Single sequence search” function on miRBase [44] (https://www.mirbase.org/) was used with the following parameters: “Search sequences: Mature miRNAs”, “Search method: “BLASTN”, “E-value cutoff: 10”, “Maximum no. of hits: 100”, “Show results only from specific organisms: Mouse”, “Word size: 4”, “Match score: + 5”, “Mismatch penalty: -4”. For miRNAs with more than one result, the miRNA with the best alignment (lowest E-value) is reported. If no miRNAs were identified in mouse, BLASTN was rerun with no species filter and the best alignment was reported. miRNAs with no results reported from any species are denoted with “NA”.

Predicting transcriptional regulation using Lisa

Transcriptional regulator (TR) binding sites were predicted using Lisa [45] (http://lisa.cistrome.org/) with genes which were female- or male-biased in 2 or more ages between PD27, PD32, PD37. The full Lisa model was applied (“TR ChIP-seq Peak-RP (regulatory potential)” and “ISD-RP (in silico deletion-regulatory potential) for both motif and ChIP-seq” methods) using the DNase-seq and H3K27Ac ChIP-seq data and 3000 genes which were randomly selected as the background gene set. Results combined from H3K27ac-ChIP-seq and DNase-seq ISD models, and TR ChIP-seq peak-only models using the Cauchy combination test are shown for the ChIP-seq model.

Pathway enrichment analysis

All pathway enrichment analyses were performed using gProfileR (v0.6.7) in R (v3.6). For differentially expressed gene pathway enrichment, all detected genes in this dataset were used as background. For miRNA–gene target pathway enrichment, all gene targets detected in this dataset were used as background with parameters “min_set_size = 3, min_isect_size = 2”.

Co-expression module identification

Gene co-expression modules were identified using CEMitool [46] (version 1.8.2) using log2-transformed, normalized read counts with default settings. Briefly, CEMitool first uses an unsupervised method to filter for genes with sufficient variation. By default, CEMitool models the variants of the genes as an inverse gamma distribution and chooses genes with a p value < 0.1. Next, it automatically determines the similarity criteria before it separates genes into modules using the dynamic tree cut method. Hub genes are identified by ranking the summed similarities between a certain gene and all other genes in the same module. This analysis was performed in R version 3.6.

miRNA–gene target correlation

Computationally predicted gene targets were curated from TargetScanMouse (v7.2) [47] and experimentally validated gene targets were curated from miRTarBase (v8.0) [48]. Only miRNA–gene target pairs from TargetScan with “Cumulative weight context score” < − 0.1 were used. “Context score” for a specific target site is defined by [47] as the summed contribution from 14 features which likely influence miRNA targeting a given gene, including “site type”, “local AU”, “3’ UTR length”, and “Probability of conserved targeting” (full feature list http://www.targetscan.org/vert_70/docs/context_score_totals.html). Spearman’s correlation coefficient (rho) was calculated for each pair using log2-transformed normalized counts and p-values were adjusted with FDR to account for multiple testing. Pairs were considered negatively correlated if rho < 0 and FDR-adjusted p-value < 0.1. This analysis was performed in R version 3.6.

Integration of single-nuclei RNA-seq dataset from [16]

Single-nuclei RNA-seq (snRNA-seq) of 10–12 week-old snap-frozen male and female C57BL/6 mouse pituitary gland was obtained from GSE151961 and processed using Seurat 4.1.0 [49]. Data were merged from replicates for each sex separately and cells with mitochondrial gene content > 15% and ribosomal gene content > 3% were removed. The percent of mitochondrial gene content and percent of ribosomal gene content were variables regressed out of the merged data for each sex during data normalization using SCTransform [50]. Data integration between sexes was performed using Seurat with the normalized merged data [51]. Cell types were labeled by comparing gene markers identified in cell clusters calculated by Seurat using the “FindAllMarkers” functions to gene markers identified in the original paper [16]. For all downstream analyses, the “Debris” cluster, which was also identified in the original paper, was removed. This analysis was performed using R version 4.1.2. For detailed methods, our code is available at https://github.com/wilsonlabgroup/pituitary_transcriptome_analyses.

Cell-type enrichment of co-expression module genes

A one-sided Kolmogorov–Smirnov (KS) test was performed to compare the distribution of expression of a given co-expression module gene within each cell type versus its distribution of expression values in all other cell types based on the expression data from the [16] snRNA-seq dataset (see Methods—Integration of single-nuclei RNA-seq dataset from [16] for data processing details). Resulting p-values were FDR-adjusted for multiple testing. Only co-expression module genes with FDR-adjusted p-value ≤ 0.05 in at least one cell type comparison were plotted. A one-sided hypergeometric test was then used to determine if there was enrichment for a group of module genes with statistically greater expression based on the KS test (FDR-adjusted p-value ≤ 0.05) in each cell type. Resulting p-values were FDR-adjusted for multiple testing and a group of module genes was considered to be significantly enriched if FDR-adjusted p-value < 0.05. This analysis was performed using R version 4.1.2.

Proportions in admixture RNA-seq deconvolution

Proportions in admixture RNA-seq deconvolution was run as part of a wrapper script in scMappR (v1.0.7) [52] using the “compare_deconvolution_methods” function which calls the ADAPTS package (v1.0.21) [53] with RUVSeq normalized bulk counts and a custom signature matrix. In scMappR, the Proportions in Admixture method is called “WGCNA”. The custom signature matrix was generated using the “generes_to_heatmap” function from scMappR with gene markers identified from sex-integrated Ruf-Zamojski snRNA-seq data (see Integration of single-nuclei RNA-seq dataset from [16] for data processing) and selecting only the top 3000 genes with the greatest variance. This analysis was performed using R version 4.1.2.

Identifying cell-type-specific sex-biased genes using scMappR

Cell-weighted fold-change (cwFC) of PD37 sex-biased genes was calculated using scMappR (v1.0.7) (https://cran.r-project.org/package=scMappR) with sex-integrated adult C57BL/6 mouse pituitary transcriptome [16] using the “WGCNA” (Proportions in Admixture) deconvolution method. Genes which were not detected in the single-cell reference dataset and genes which were flagged as a “false positive” by scMappR (cwFoldchange_gene_flagged_FP) were filtered out. Genes were further filtered to include outliers as determined by scMappR based on their cwFC (cwFoldchange_gene_assigned). Finally, only genes with an absolute gene-normalized cwFC > 0.5 for a given cell-type were considered cell-type-specific sex-biased genes and shown in the heatmaps. This analysis was performed using R version 4.1.2.

Results

Profiling postnatal mouse pituitary gland development with 3’UTR-seq and small RNA-seq

To assess changes in the mouse pituitary transcriptome across postnatal development, we profiled pituitary RNA expression at 5 postnatal days spanning the pubertal transition (PD: 12, 22, 27, 32 and 37). We observed physical markers of pubertal onset, preputial separation (PS) and vaginal opening (VO) occurring on average at PD27 and PD29 in male and female mice, respectively, in our C57BL/6J colony ([54], Fig. 1A). In all analyses performed, pubertal onset refers to ages at which PS and VO were recorded (see Fig. 1B for a summary of analysis workflow).

Fig. 1
figure 1

Overview of the pituitary transcriptome during postnatal development in male and female samples. A Schematic of experimental design. Marks (purple) on the timeline denote the age in which the pituitary gland was collected. Vertical arrows denote average ages for onset of puberty of the specified sex in our colony (determined by preputial separation for males or vaginal opening for females). Fraction of pubertal mice out of total mice for males (blue) and females (red) at each age is shown. B Schematic of analysis workflow. Summary of miRNA expression analyses (yellow), gene expression analyses (blue), miRNA–gene target identification (green), processing of single-nuclei RNA-seq (snRNA-seq) data from [16] (red), and combining snRNA-seq data with bulk gene expression data (purple). C Genome browser screenshots showing QuantSeq signal at Fshb and Prop1. X-axis: genomic coordinates; y-axis: reads per million mapped reads (RPM); PD postnatal day. Gene name and gene model are shown on the bottom of each panel. Each track represents overlapping signal from 5–6 biological replicates. D Scatter plot showing the correlation between gene quantification measured by qPCR and by 3’UTR-seq in one pituitary sample (PD37M4). X-axis: ΔCt values obtained by qPCR; y-axis: log2-transformed normalized counts (log2(normCounts)) values obtained by 3’UTR-seq. Sample name, Spearman correlation coefficient, and number of genes included are labeled on the plot. PCA plot for pituitary gland samples based on E gene expression and F miRNA expression. Principal component analysis (PCA) was performed using log2(normCounts) after filtering for low-count genes/miRNAs and normalization using RUVSeq. Only scores of the first 2 PCs are shown. Age is indicated by shape while sex is indicated by colors

Full size image

To measure mRNA expression in a genome-wide, cost-effective, and relatively high-throughput manner, we first automated the QuantSeq 3’ mRNA-Seq protocol which profiles 3’UTR of mRNA transcripts (Methods—Library preparation and sequencing). Of the 60 libraries generated, 55 libraries were of sufficient library quantity for sequencing which resulted in 4 to 6 biological replicates for each sex at each of the five postnatal days (see Additional file 2: Table S1 for quality control metrics).

It was shown previously that single-cell-based 3’UTR profiling in the pituitary missed the expression of Prop1 due to a lack of tissue-specific 3’UTR gene annotation [55]. To avoid similar issues with our 3’UTR profiling, we refined gene 3’-end annotations by identifying clusters of sequencing reads from our data and re-annotating 3’UTRs (Methods—PolyA site identification and gene annotation modification and Additional file 1: Fig. S1A). Improved pituitary-specific 3’UTR annotations were generated for 676 genes, allowing for assignment of significantly more reads to them, including important pituitary genes such as Pou1f1, Ghrhr, Fshb, and Prop1 (Fig. 1C, Additional file 1: Fig. S1B).

To profile miRNA expression, we performed sRNA-seq in the same male and female samples used for 3’UTR-seq besides PD37 which was excluded from this experiment (PD: 12, 22, 27, and 32; n = 6 biological replicates; 48 libraries total). Using the miRDeep2 workflow, we identified 273 known mouse miRNAs (miRBase v21) and 19 novel miRNAs (Additional file 3: Table S2).

For both 3’UTR-seq and sRNA-seq experiments, the biological replicates were well correlated (Pearson’s correlation coefficient: 3’UTR-seq 0.95–0.97; sRNA-seq 0.86–0.90) (Additional file 1: Fig. S2A, B). Furthermore, gene expression level quantified by 3’UTR-seq correlated well with microfluidic qPCR data previously generated from the same 55 RNA samples (178 puberty-related genes plus 5 control genes) [17] (median Spearman correlation coefficient: 0.74) (Fig. 1D, Additional file 1: Fig. S1C). Using principal component analyses (PCA), we observed a separation of PD12 samples from PD22 and older for gene expression profiles along PC1 (Fig. 1E). Although more subtle, we observed samples distributed by age along PC1 for miRNA expression profiles (Fig. 1F). We also observed separation between male and female samples along PC2 at all ages based on gene expression profiles, which became more pronounced across postnatal ages (Fig. 1E). In contrast, no obvious sex differences were observed in our miRNA expression data (Fig. 1F).

Sex-biases in the transcriptome occur prior to puberty

We quantified sex differences in the pituitary transcriptome by comparing the expression of genes and miRNAs between male and female samples at each age. Across all the profiled ages, we observed an increase in the numbers of sex-biased genes and miRNAs, with the most dramatic increase occurring at PD27, when all the males and half of the females had gone through puberty as measured by PS and VO (Figs. 1A, 2A, B see Tables 1, 2 for list of significant sex-biased genes and miRNAs; see Additional file 4: Table S3, Additional file 5: Table S4 for full list of differential analysis results).

Fig. 2
figure 2

Pituitary transcriptome is increasingly sex-biased across postnatal development. Barplot showing the number of A intersecting sex-biased genes or B intersecting sex-biased miRNAs between each age. Horizontal bars on the bottom left side of each plot show the numbers of male- (blue) or female-biased (red) genes/miRNAs at each age (absolute FC > 1.5; FDR < 0.05). Different intersection combinations between sex-biased genes/miRNAs identified at each age are represented by the dotplot. The number of genes/miRNAs which intersect in the indicated combination of sex comparisons is shown by the vertical barplots (# overlapping sex-biased genes/miRNAs). Expression plots of example C pre-pubertal (PD12-22) sex-biased genes and genes with sex-by-age effect between PD12 and PD22, D peri-/post-pubertal sex-biased genes, and E sex-biased miRNAs. Log2-transformed normalized counts (log2(normCounts)) are plotted for each gene/miRNA. Expression changes are shown across ages (x-axis). Large, filled points represent median expression at each age and unfilled points represent each biological replicate. Blue: male samples; red: female samples. Black asterisks highlight ages at which the corresponding genes/miRNAs are detected as sex-biased and red asterisks highlight genes with significant sex-by-age effect between PD12 and PD22. Network representation of pathways enriched for F male-biased and G female-biased differentially expressed genes. Each node represents a pathway and nodes are connected based on similarity in genes found enriching for the connected pathways (Jaccard distance). Node size represents the number of differentially expressed genes enriching for the given pathway and nodes are colored based on the differential expression comparison in which the genes were identified. Pathways that share similar genes are circled (dashed lines) and labeled manually based on pathway functions. Specific pathways and their manual labels can be found in Additional file 6: Table S5

Full size image
Table 1 Sex-biased genes at each profiled age
Full size table
Table 2 Sex-biased miRNAs at each profiled age
Full size table

Although we observed most sex-biased gene expression at peri- and post-pubertal ages, sex differences in the pituitary begin to manifest earlier in postnatal development at PD12 (Table 1). At PD12, 12 genes, including seven sex chromosome-linked genes: Ddx3y, Uty, Kdm5d, Gm29650, Eif2s3y (Y-linked), Xist and Kdm6a (X-linked); and five autosomal genes, Chrna4, Kcna4, Lhb, Th, and Drd4, are identified as significantly sex-biased (FDR < 0.05, absolute fold-change > 1.5) (Fig. 2C, Additional file 1: Fig. S3A). The seven sex chromosome-linked genes were male- or female-biased across all five profiled ages (PD12 to PD37). Other than Kdm5d, Ddx3y, Eif2s3y, Uty, Xist, and Lhb [9, 18, 56], sex differences in the expression of other sex-biased genes we detected at PD12 have not been reported previously in pituitary.

We identified 25 male-biased and 16 female-biased genes at PD22 (Table 1), an age preceding puberty, including puberty-related genes Dgkk, Fshb, Osbp2, and Pcsk1, which we previously identified as sex-biased using microfluidic qPCR [17] (Fig. 2C, Additional file 1: Fig. S3A). Of these 41 sex-biased genes at PD22, we found 18 genes which start to exhibit sex-biased expression at PD22 and maintain the same sex-biased trend at all older ages we profiled, (e.g., Asb4 and Serpine2) (Fig. 2A, C, Additional file 1: Fig. S3A), demonstrating that there is establishment of sex-biased expression in the pituitary prior to puberty.

We noticed some genes exhibited sex-biased changes between PD12 and PD22 (e.g., Fshb). We then specifically tested for sex-by-age interaction effect, and identified 13 genes (Fshb, Pcsk1, Serpine2, Lhb, Chrna4, Nupr1, Dgkk, Steap3, Timp1, Kcna4, Gpr101, 2010007H06Rik, and Th) with significant sex-by-age interaction effect between PD12 and PD22 (FDR < 0.05). All of these genes showed similar or increased expression in females compared to males at PD12 but displayed decreased expression in females compared to males at PD22 (Fig. 2C, Additional file 1: Fig. S3A).

For miRNAs, we identified 3 miRNAs with sex-biased expression prior to puberty (PD12 or PD22). The one sex-biased miRNA at PD12 displayed female-biased expression and was one of 20 miRNAs (after normalization and filtering out low-count miRNAs; Additional file 3: Table S2) not previously identified in mice based on miRBase v21 (we tentatively named it novel46; Fig. 2E, Table 2, Additional file 1: Fig. S3B). We did not observe any alignments of novel46 to miRNAs in other species found in miRBase v21. Based on its miRNA precursor coordinates, novel46 is a mirtron (a miRNA which is spliced from a host gene intron) expressed from the last intron of calcium voltage-gated channel subunit alpha1 G (Cacna1g), located on chromosome 11 (Additional file 3: Table S2, Additional file 1: Fig. S3C).

At PD22, we found that miR-499-5p displayed male-biased expression, and based on its miRNA precursor coordinates, miR-449-5p is a mirtron expressed from Myh7b (Table 2, Additional file 1: Fig. S3B). We also found that miR-342-5p displayed female-biased expression at PD22, and based on its miRNA precursor coordinates, miR-342-5p is a mirtron expressed from Evl (Fig. 2E, Table 2, Additional file 1: Fig. S3B). However, neither Myh7b nor Evl were detected as significantly sex-biased at any profiled age.

Peri- and post-pubertal sex differences in gene expression reflect sex differences in pituitary endocrine functions

At peri- and post-pubertal stages, we detected similar numbers of male- and female-biased genes. The number of sex-biased genes roughly doubled across the pubertal transition (140, 253, and 351 male-biased genes and 156, 232, and 294 female-biased genes at PD27, PD32, and PD37, respectively; Fig. 2A, Table 1, Additional file 4: Table S3). Many of these genes (43 male and 73 female-biased) showed sex-biased gene expression throughout the pubertal transition (across PD27, PD32, and PD37) (Fig. 2A, D, Additional file 1: Fig. S3A). While we recovered genes with previously known sex-biased expression in pituitary, such as Dlk1 and Prl [9, 17, 57], many other genes whose sex-biased expression has not been previously identified in the pituitary were found (Fig. 2D, Table 1), including the male-biased expression of an abundant transcript in the pituitary, neuronatin (Nnat) [58], at PD27, PD32 and PD37.

Male-biased genes are enriched for pathways related to hormone synthesis and secretion (Fig. 2F, Additional file 6: Table S5), including “peptide hormone biosynthesis” (PD27, p = 7.22e−06), “regulation of secretion” (PD37, p = 7.87e−03), and “secretory vesicle” (PD37, p = 4.9e−04); and pathways related to reproduction, including “male sex differentiation” (PD37, p = 4.76e−02) and “reproductive process” (PD27, p = 4.14e−02). In addition, male-biased genes are enriched for pathways associated with signaling receptors, ion channels, extracellular region, and plasma membrane, like “G-protein coupled receptor signaling pathway”, which is enriched at all three ages (PD27: p = 4.7e−02; PD32: p = 1.22e−03; PD37: p = 8.67e−03). Finally, pathways related to endopeptidase inhibitor activity are also enriched, including several serpin family genes (Serpine2, Serpina3c, and Serpinb1a).

Female-biased genes, including Stat5a, Cckbr, Slit2, Robo2, Nrip1, Nhlh2, Prl, and Pgr, are enriched for “ovulation cycle” (PD37: p = 1.03e−02), linking female-biased pituitary genes to female-specific physiological processes. Notably, other female-biased pathways are predominantly neuron-related (Fig. 2G, Additional file 6: Table S5); this could be attributed to genes expressed in the posterior pituitary, which contains axons extended from the hypothalamus, or genes expressed in neuroendocrine cells, which are known to be activated in a neuron-like manner [59,60,61]. Particularly, female-biased genes at all three ages, including Crhbp, Prl, Calb1, Pgr, Pak7, Reln, Dmrta1, Prkce, Ar, Nhlh2, Grm5, and Cacna1c, are enriched for “behavior” (PD27: p = 5.32e−03; PD32: p = 4.25e−02; PD37: p = 5.32e−03). Several of these genes, Crhbp [62], Prl [63], Calb1 [64], Reln [65], Prkce [66], Ar [67], Nhlh2 [68], Grm5 [69], and Cacna1c [70], are related to the regulation of stress, which is sex-biased in its activity [71].

For miRNAs, we detected 3 and 4 male-biased miRNAs and 3 and 8 female-biased miRNAs at PD27 and PD32, respectively (Fig. 2B, Table 2, Additional file 1: Fig. S3B, Additional file 5: Table S4). miR-224-5p, a mirtron of Gabre, and miR-383-5p, a mirtron of Sgcz, respectively, displayed male- and female-biased expression levels at both PD27 and PD32 (Fig. 2E, Table 2).

Connecting sex-biased miRNAs to target genes with sex-biased expression

To evaluate post-transcriptional regulation of sex-biased changes in the pituitary gland by miRNAs, we identified computationally predicted or experimentally validated miRNA–target gene pairs that both exhibit sex-biased expression. Since miRNAs are usually predicted to repress target gene expression, we additionally required miRNA and target gene pairs to be significantly negatively correlated in expression across matched samples (Spearman’s rho < 0, FDR < 0.1) (see Methods—miRNA–gene target correlation for details) (Fig. 3A, Additional file 7: Table S6).

Fig. 3
figure 3

Regulation of sex-biased pituitary gene expression by miRNAs and transcriptional regulators (TRs). A Interaction network of sex-biased miRNAs and their sex-biased target genes for all ages. Each node represents a gene (gray) or a miRNA (orange) with the node size representing its number of connections with other nodes. Outlines on nodes indicate gene/miRNA is female-biased (red) or male-biased (blue). Only miRNAs with 2 or more connections are shown. Edges between each node show a predicted interaction between the miRNA and gene. Edge thickness indicates Spearman’s correlation coefficient (rho) calculated for the given pair. B Heatmap of sex-biased pituitary TR gene expression. TRs which are sex-biased in the same direction at a minimum of two ages between PD27, PD32, and PD37 are plotted. Colors of the heatmap represent row-scaled and centered expression levels of each gene. Column annotation bars indicate sample age and sex. Row annotation bars indicate age at which the gene was found to be sex-biased. C Scatterplot comparing Lisa TR rankings with combined P-values predicted to regulate female- and male-biased genes. Each TR is represented by a point. The combined − log10(P-value) is plotted for each TR based on female-biased and male-biased gene sets which are sex-biased in the same direction at a minimum of two ages between PD27, PD32 and PD37. Colored points show TRs which have change in − log10(p-value) between sexes which is two standard deviations greater than the mean change in − log10(p-value). Red points indicate TRs which are enriched for regulating female-biased genes; blue points indicate TRs which are enriched for regulating male-biased genes

Full size image

At pre-pubertal ages, PD12 and PD22, we did not identify any negatively correlated miRNA–gene pairs where the miRNA and its target gene were sex-biased. This is likely due to the relatively lower number of sex-biased genes and miRNAs identified at these earlier profiled ages.

At peri- and post-pubertal ages, we found 18 putative sex-biased gene targets of 5 sex-biased miRNAs. For example, the male-biased miRNAs, miR-181a-5p (at PD27) and miR-224-5p (at PD27 and PD32), showed significant negative correlation with female-biased target genes, Ammecr1 (target of miR-181a-5p and miR-224-5p), Ank1 (target of miR-181a-5p), and Prkce (target of miR-181a-5p) (Fig. 3A, Additional file 7: Table S6). In comparison, female-biased miR-205-5p (at PD32) was negatively correlated with Inhba (Fig. 3A, Additional file 7: Table S6).

Predicting transcriptional regulators associated with sex-biased pituitary gene expression

Sex-biased expression of transcriptional regulators (TRs) is a principal mechanism by which sex-biased regulatory networks can be generated. Of the sex-biased genes we detected, 5 male-biased and 16 female-biased genes are annotated as TRs in mouse by AnimalTFDB3 [72] (Fig. 3B). Known molecular functions in the pituitary and pituitary-related knockout phenotypes of these 21 sex-biased TRs are summarized in Table 3.

Table 3 Sex-biased TRs and their pituitary-related phenotypes and functions
Full size table

To gain insights into how TR proteins might regulate sex-biased pituitary gene expression across pubertal transition, we used the bioinformatic tool Lisa (epigenetic Landscape In Silico deletion Analysis; [45]). Lisa takes a list of genes and builds transcriptional regulatory models based on both publicly available DNase and ChIP-seq datasets and returns predictions of candidate TRs that regulate them. To focus on the major increase in sex-biased gene expression occurring at peri-pubertal and post-pubertal ages (PD27, PD32, and PD37), we provided Lisa with pituitary genes classified as having sex-biased gene expression at two or more postnatal days between PD27, PD32, and PD37 in males (n = 183) or females (n = 174). We identified 6 and 22 TRs (Δ-log10(p-value) > mean ± 2SD) including gonadal nuclear hormone receptors, ESR1 and AR, that were predicted by Lisa to regulate female-biased and male-biased genes, respectively. In addition, CBX7, a component of Polycomb repressive complex 1 (PRC1), as well as SUZ12 and EZH2, components of Polycomb repressive complex 2 (PRC2), were predicted by Lisa to regulate female-biased genes. In contrast, a heterodimeric TR, CLOCK:BMAL1 (ARNTL), was predicted by Lisa to regulate male-biased gene expression (Fig. 3C).

Gene co-expression analysis reveals dynamic modules enriching for pituitary cell types

To further characterize developmental changes in the pituitary transcriptome, we utilized our temporal and sex-dependent gene expression profiling of the postnatal pituitary to build gene co-expression networks because co-expressing genes tend to share related functions or regulatory pathways. We used CEMiTool [46] to identify transcriptome-wide gene expression correlation networks in the postnatal pituitary transcriptome (Methods—Co-expression module identification). In total, nine co-expression gene modules were identified based on the expression of 1205 genes (Fig. 4, Additional file 1: Figs. S4, S5, Additional file 8: Table S7).

Fig. 4
figure 4

Co-expression network analysis identifies gene modules underlying pituitary transcriptome changes. Heatmap shows genes selected for co-expression analysis, separated into 9 modules. Column annotation bars indicate sample age and sex. Colors of the heatmap represent row-scaled and centered expression levels of each gene. A summary of the expression profiles of each module is shown in the left column of the right panel. The solid line represents the median expression (scaled and centered as shown in the heatmap) at each age and sex (red: female samples, blue: male samples) for all the genes in the corresponding module. Dash lines represent the scaled expression profiles of each gene in the module. Module names and number of genes included in the module are labeled. Top 10 hub genes (calculated based on genes’ connectivity within the modules) are listed for each module. All co-expressing module genes are listed in Additional file 8: Table S7. Cell types in which module genes are enriched in based on single-nuclei RNA-seq expression from [16] are shown on the right (see Additional file 1: Fig. S6C and Methods—Cell-type enrichment of co-expression module genes for details). One-sided hypergeometric test: * FDR < 0.05, ** FDR < 0.01, *** FDR < 0.001

Full size image

The temporal expression profiles of the 9 co-expression gene modules showed three general patterns: (1) decreasing expression, particularly between PD12 and PD22 (Modules 1 (M1), M4, M7, and M8); (2) increasing expression (M2 and M3); and (3) sex-biased expression (M5, M6, and M9) (Fig. 4, Additional file 1: Fig. S5A-B).

To further understand the nature of the co-expression modules, we used three independent approaches: (1) pathway enrichment analysis; (2) hub gene identification; and (3) an enrichment test where we ask if a set of genes in a module show enhanced expression in specific adult pituitary cell types (previously determined by single-nuclei RNA-seq (snRNA-seq) [16]).

Using this approach, we found that our co-expression modules contain clear cell-type-enriched gene expression signatures. For example: (1) the M7 module, shows decreased expression during postnatal development in both males and females and is enriched for cell cycle-related pathways (Additional file 1: Fig. S5C); (2) the top M7 module hub genes include Top2a, Mki67, Cdca3, Cenpf, and Prc1, all of which are canonical markers of proliferating cells (Additional file 1: Fig. S4); and (3) using published snRNA-seq datasets of adult female and male mouse pituitary gland [16] (Additional file 1: Fig. S6A-B), we found that there is a significant enrichment of M7 genes in the proliferating cell population (Hypergeometric test, FDR = 2.60e-41) (Fig. 4, Additional file 1: Fig. S6C, Methods—Cell-type enrichment of co-expression module genes).

Similarly, the top hub genes for M2 and M5 include Gh and Dlk1 (Fig. 4), markers of somatotropes [57], and there was a significant enrichment of both M2 (FDR = 6.84e−11) and M5 (FDR = 2.82e−10) module genes in somatotropes (Additional file 1: Fig. S6C). The top hub gene of M3 is Prl, a marker of lactotropes, and there is a significant enrichment of M3 module genes (FDR = 3.18e−24) in the lactotropes. Finally, the top hub genes in M8 consist of pituicyte markers, including Col13a1, Scn7, and Col25a1 [73], there is also a significant enrichment of M8 module genes (FDR = 1.61e−17) within the pituicytes (Additional file 1: Fig. S6C). Overall, we observe a clear association between cell types and specific gene co-expression modules identified in the pituitary gland.

Sex differences in cell proportions emerge prior to and across puberty

Although sex differences in the number of somatotropes, lactotropes and gonadotropes have been documented in the adult pituitary gland [15, 16, 74], much less is known about when these sex differences arise during postnatal development. To assess whether sex differences in cell proportions of somatotropes, lactotropes and gonadotropes are dynamic across pubertal development, we estimated cell-type proportions in our temporal bulk RNA-seq using the Proportions in Admixture RNA-seq deconvolution algorithm [75] from Automated Deconvolution Augmentation of Profiles for Tissue Specific cells (ADAPTS) [53] (Fig. 5A, B, Additional file 1: Fig. S6D). All cell types detected in the sex-integrated single-nuclei adult female and male pituitary transcriptome were used as our reference [16] (Additional file 1: Fig. S6A, B).

Fig. 5
figure 5

Estimating sex differences in pituitary cell types by leveraging single-nuclei RNA-seq. A Schematic of analysis workflow for estimating sex differences in pituitary cell types using bulk gene expression with single-nuclei RNA-seq data from [16]. B Estimated cell-type proportions by RNA-seq deconvolution using Proportions in Admixture (WGCNA) changes across profiled ages of cell types with previously established sex-biased proportions in the adult pituitary. Estimated cell-type proportions are plotted across postnatal ages along the x-axis. Large circles and triangles represent the mean cell-type proportion at each age and small circles and triangles represent each biological replicate. Lighter color, solid line, circle points: female samples; dark color, dotted line, triangle points: male samples. Wilcoxon test was performed to compare cell proportions between both sexes at each age (*p < 0.05, **p < 0.01). See Additional file 1: Fig. S6D for all other pituitary cell type proportions. C Heatmap of cell-weighted fold-changes (cwFC) for sex-biased genes at PD37 in somatotropes, lactotropes, and gonadotropes. Color gradient indicates gene-normalized cwFC value calculated by scMappR; red: more female-biased; blue: more male-biased. Genes with |gene-normalized cwFC|> 0.5 in at least one cell type are plotted. Genes shown to be sex-biased in the same direction in the same cell-type by [14] are indicated by the purple text and an asterisk within the cell of the corresponding cell-type

Full size image

We found that we were able to recapitulate known sex biases in adult pituitary cell-type proportions [15, 16, 74] at our oldest profiled age, PD37, where the estimated cell-type proportions for somatotropes were significantly male-biased, and lactotropes were significantly female-biased (Fig. 5B). In addition, we identified sex differences in estimated cell-type proportions which were dynamic across our earlier profiled ages. For example, we observed significant sex differences in estimated gonadotrope cell-type proportions between PD12 and PD22, ages prior to puberty, and these proportions show a sex-by-age trend (from PD12 to PD22, decreasing proportions in females and increasing proportions in males), which could contribute to some of the sex-by-age bulk gene expression we observed (Fig. 5B). At PD27, when puberty has occurred in a subset of our mice (Fig. 1A), we observed the emergence of male-bias in somatotrope cell proportions and female-bias in lactotrope cell proportions (Fig. 5B).

Inferring sex-biased gene expression in pituitary cell types

We next examined if any sex-biased genes identified remained sex-biased after adjusting for cell-type proportion sex differences in somatotropes, lactotropes, and gonadotropes. To do this, we used a bioinformatic pipeline, Single-cell mapper (scMappR) [52] to calculate a cell-weighted fold-change (cwFC) for each sex-biased gene identified in our bulk RNA-seq data (Fig. 5A). Briefly, cwFCs were calculated for each sex-biased gene in somatotropes, lactotropes, and gonadotropes, by readjusting the bulk gene expression fold-change with the cell-type specificity of the gene and the ratio of deconvolution-estimated cell-type proportions (determined by snRNA-seq). To minimize potential confounding effects from differences in ages, we focused on sex-biased genes identified at PD37, which is the closest age to the snRNA-seq reference dataset from adult male and female mice [16]. We identified 75 sex-biased genes which remain sex-biased after adjusting for sex differences in cell-type proportions (|cwFC|> 0.5) across somatotropes, lactotropes, and gonadotropes (Fig. 5C, Additional file 9: Table S8). Of these 75 sex-biased genes, 4 genes (Gata2, Rab3b, Vegfa, and Lama1) were found to be concordantly sex-biased in the corresponding cell type in the rat anterior pituitary (genes highlighted in Fig. 5C, [14]). Thus, by combining our temporal postnatal bulk pituitary 3’UTR-seq with a snRNA-seq dataset we can infer postnatal pituitary cell-type specificity of sex-biased genes and gain insights into sex-biased temporal gene regulation in a cell-type-specific manner.

Discussion

The pituitary gland displays sex differences in its regulated physiological functions, including stress response, somatic growth, reproduction, and pubertal timing. In our previous work, we showed by qPCR profiling that selected genes associated with the onset of puberty were increasingly sex-biased with pubertal development in the mouse pituitary gland [17]. In this study, we aimed to further understand if sex bias exists beyond puberty-related genes and identify potential regulatory mechanisms. By comparing transcriptome-wide pituitary gland gene and miRNA expression profiles between male and female mice during postnatal development, we have identified sex-biased genes and miRNAs that contribute to sex differences in postnatal pituitary development and function.

While the majority of sex differences in pituitary gene expression were observed at puberty and beyond, we observed robust sex-biased expression of genes and miRNAs prior to onset of physical markers of puberty (VO/PS) at PD12 and PD22. These sex differences may be attributed to irreversible organizational effects [76] that are established by gonadal hormone exposures during fetal development and perhaps in the first week after birth in mice [64]. Unlike activational effects, sex differences established by organizational effects persist even when gonadal hormones are present at low levels in the system, such as at pre-pubertal ages. We found that many sex-biased genes identified at PD12 or PD22 or genes showing sex-by-age effect between PD12 and PD22 are expressed in the gonadotropes, including Fshb and Lhb (which encode for β-subunits of gonadotrope-secreted hormones), Chrna4 (which is a marker for gonadotropes in rats [14]), Nupr1 (which is involved in embryonic gonadotrope development [77]), as well as Gpr101 and Serpine2 (which display sex-biased gene expression in gonadotropes isolated from juvenile mice but not adult mice [18]). Our RNA-seq deconvolution analyses also suggest that gonadotrope cell proportions are male-biased prior to puberty. As gonadotropes produce LH and FSH to regulate gonadal maturation during puberty, these pre-pubertal sex differences in gonadotrope cell proportions and sex-biased expression suggest that these sex differences may be necessary to prime the pituitary gland for its regulation of gonadal maturation during puberty.

By identifying a list of sex-biased genes, we were able to predict trans regulators that are likely to underlie their expression. Expectedly, we predicted by Lisa that ESR1 and AR are involved in regulating female- and male-biased gene expression. However, since Lisa does not differentiate between activating and repressing effects of TRs on gene expression, Lisa-predicted regulation of male-biased genes by AR suggests that AR is either activating their expression in males or is repressing their expression in females. Our observation of female-biased Ar gene expression suggests the latter. Pituitary-specific loss of AR in male mice result in increased Prl gene expression and serum PRL levels [78], and gonadotrope-specific loss of AR in female mice result in decreased Fshb gene expression and serum FSH and LH levels [79]. In addition, global AR ablation affects its repression of sex-biased genes in the liver by altering methylation levels at the promoter regions [80], whether this has a similar effect on pituitary sex-biased gene expression remains to be examined. We also predicted by Lisa that female-biased gene expression involved regulation by polycomb complexes PRC1 (CBX7) and PRC2 (SUZ12 and EZH2). This observation is consistent with results obtained from the GTEx consortium who compared adult human male and female pituitary gene expression [13]. Mechanistically, sex-biased deposition of H3K27me3 marks by PRC2 through its catalytic subunits Ezh1/Ezh2 has been shown to repress female-biased expression in the male adult mouse liver [81] and is also involved in the regulation of pubertal onset in female rats by inhibiting Kiss1 expression in the arcuate nucleus prior to the initiation of puberty [82]. Our study highlights that future studies of these mechanisms in the mouse pituitary would be ideal to study prior to puberty.

miRNAs have potentially important roles in responding to sex-biased hormonal release across multiple hypothalamic–pituitary axes. For example, multiple miRNAs have been shown to be estrogen-responsive in the neonatal mouse hypothalamus by the administration of an aromatase inhibitor [26]. miR-1948 and miR-802, identified as sex-biased in adult mouse liver, are known to be regulated by sex-biased pattern of growth hormone release from the pituitary gland [25]. Here we identified two male-biased miRNAs, miR-181-5p and miR-224-5p with estrogen-responsive predicted targets: protein kinase C​​ε (encoded by Prkce) whose signaling in response to epinephrine-induced inflammatory pain was found to be suppressed by estrogen in rats [83] and Ammecr and Ank1 whose gene expression is upregulated by estradiol in mouse pituitary [84]. We also identified a female-biased novel miRNA at PD12, novel46, whose host gene, Cacna1g, has been shown in anterior pituitary primary cells to promote LH secretion by estradiol signaling through estrogen receptor 1 (ESR1) upon GnRH-induction [84]. While additional functional studies of these miRNAs are needed, these putative miRNA–mRNA connections represent potential gene regulatory networks underlying sex differences in postnatal pituitary gland development.

Sex differences in the number of somatotropes and lactotropes in the adult pituitary have been known for more than two decades [74]. These findings have recently been expanded on by single-cell genomic analyses of the adult pituitary [14,15,16]. In this study, we combined our temporal bulk RNA-seq data with snRNA-seq data from adult pituitaries to estimate sex differences in cell-type proportions during postnatal development. We suggest that these sex differences in cell proportions occur prior to puberty. It was previously demonstrated that the gonadotrope population displays plasticity in response to reproductive processes, for example, the density of gonadotropes increases in post-pubertal female mice (8–15 weeks) relative to pre-pubertal ages (3 weeks) and the localization of gonadotropes change in lactating mice [85]. In addition, neonatal exposure to testosterone is known to influence the number of somatotropes and lactotropes in the adult male and female rat pituitary [86, 87]. There are several proposed mechanisms by which cell composition changes arise in the pituitary gland: differentiation from the adult pituitary stem cell niche [88], transdifferentiation from differentiated cell-types [15, 89], and self-proliferation from existing cell-types [90, 91]. Future work such as single-cell multiomic profiling of female and male pituitaries during early postnatal development will be needed to pinpoint the gene regulatory networks that govern pituitary cell type composition.

By using bulk mRNA and miRNA profiling of multiple replicates of male and female mice at multiple postnatal timepoints, we identified sex biases in gene expression which are not explained by differences in cell type proportion. Several of these genes have plausible links to regulating sex biases in specific pituitary cell types. For example, the male-biased expression of Rab3b in lactotropes is consistent with its inhibitory role of PRL secretion in males [92]. In addition, the male-biased expression of Gata2 in gonadotropes is consistent with male-specific growth deficiency observed with pituitary-specific Gata2 KO mice [93] and the male-specific reduction in serum FSH in gonadotrope-specific Gata2 KO mice [94]. Our analysis also reveals potential novel sex-biased roles for genes. For example, Nr4a1 may be female-biased in gonadotropes, and this gene has previously been shown to be upregulated in response to GnRH [95].

Early evidence of the power of single-cell genomics for studying gene regulatory networks active in postnatal pituitary comes from Ruf-Zamojski et al. who used snRNA-seq and snATAC-seq to characterize the sex-biased specific regulatory landscape of the male and female adult mouse pituitary (10–12 weeks) [16]. Particularly, Ruf-Zamojski et al. highlighted a latent variable (LV) showing increased expression/accessibility in females. We found that 10 of the top 30 genes associated with this LV also show evidence of female-biased expression in our study (Ankra2, Crhbp, Ddx21, Ern1, Gadd45g, Greb1, Npr2, Nrg4, Rps6ka2, Stat5a). Technological advances in single-cell small RNA sequencing (reviewed in [96]), and in particular single-cell miRNA–mRNA co-sequencing techniques [97], will likely permit a higher-throughput way to assign miRNAs to specific cell types and study relationships with their target genes.

Perspectives and significance

We performed comprehensive profiling of pituitary gene and miRNA expression across the pubertal transition, a major postnatal developmental milestone, in both male and female mice. We have made this resource freely available so that others can use this data to design experiments to further understand the biology of sex differences in the pituitary gland. Our initial analysis of these data identified novel pituitary-expressed sex-biased genes and miRNAs that cannot be explained by differences in cell-type proportions alone. We also revealed dramatic developmental changes between PD12 and PD22, the mechanisms of which remain to be identified. By providing evidence for pre-pubertal sex differences in pituitary cell-type proportions, another specific challenge is dissecting the mechanisms which generate these differences. The unique role of the pituitary gland in reproductive and stress-related processes that differ between sexes highlights the importance of such future work.

Availability of data and materials

The datasets generated and analyzed during the current study are available in the ArrayExpress repository under accession numbers: E-MTAB-9459, E-MTAB-9460. 3’UTR-seq bigWig files are available in E-MTAB-9459 for genome browser visualization. All datasets and scripts used to perform analyses in the paper and reproduce figures can be found here: https://github.com/wilsonlabgroup/pituitary_transcriptome_analyses.

References

  1. Ooi GT, Tawadros N, Escalona RM. Pituitary cell lines and their endocrine applications. Mol Cell Endocrinol. 2004;228:1–21. https://doi.org/10.1016/j.mce.2004.07.018.

    Article  CAS  PubMed  Google Scholar 

  2. Baylis PH, Ball S. The neurohypophysis: endocrinology of vasopressin and oxytocin. In: De Groot LJ, Beck-Peccoz P, Chrousos G, Dungan K, Grossman A, Hershman JM, Koch C, McLachlan R, New M, Rebar R, Singer F, Vinik A, Weickert MO, Editors. Endotext. MDText.com, Inc., South Dartmouth (MA); 2000.

  3. Bucy PC. The pars nervosa of the bovine hypophysis. J Comp Neurol. 1930;50:505–19. https://doi.org/10.1002/cne.900500209.

    Article  Google Scholar 

  4. Fauquier T, Lacampagne A, Travo P, Bauer K, Mollard P. Hidden face of the anterior pituitary. Trends Endocrinol Metab. 2002;13:304–9.

    Article  CAS  Google Scholar 

  5. Yoshida S, Kato T, Yako H, Susa T, Cai LY, Osuna M, Inoue K, Kato Y. Significant quantitative and qualitative transition in pituitary stem / progenitor cells occurs during the postnatal development of the rat anterior pituitary. J Neuroendocrinol. 2011;23:933–43. https://doi.org/10.1111/j.1365-2826.2011.02198.x.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Waxman DJ, O’Connor C. Growth hormone regulation of sex-dependent liver gene expression. Mol Endocrinol. 2006;20:2613–29. https://doi.org/10.1210/me.2006-0007.

    Article  CAS  PubMed  Google Scholar 

  7. Agustsson TT, Baldvinsdottir T, Jonasson JG, Olafsdottir E, Steinthorsdottir V, Sigurdsson G, Thorsson AV, Carroll PV, Korbonits M, Benediktsson R. The epidemiology of pituitary adenomas in Iceland, 1955–2012: a nationwide population-based study. Eur J Endocrinol. 2015;173(5):655–64. https://doi.org/10.1530/EJE-15-0189.

    Article  CAS  PubMed  Google Scholar 

  8. Mindermann T, Wilson CB. Age-related and gender-related occurrence of pituitary adenomas. Clin Endocrinol (Oxf). 1994;41(3):359–64. https://doi.org/10.1111/j.1365-2265.1994.tb02557.x.

    Article  CAS  PubMed  Google Scholar 

  9. Bjelobaba I, Janjic MM, Kucka M, Stojilkovic SS. Cell type-specific sexual dimorphism in rat pituitary gene expression during maturation. Biol Reprod. 2015;93:21. https://doi.org/10.1095/biolreprod.115.129320.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Nishida Y, Yoshioka M, St-Amand J. Sexually dimorphic gene expression in the hypothalamus, pituitary gland, and cortex. Genomics. 2005;85:679–87. https://doi.org/10.1016/j.ygeno.2005.02.013.

    Article  CAS  PubMed  Google Scholar 

  11. Gershoni M, Pietrokovski S. The landscape of sex-differential transcriptome and its consequent selection in human adults. BMC Biol. 2017;15:7. https://doi.org/10.1186/s12915-017-0352-z.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Lopes-Ramos CM, Chen C-Y, Kuijjer ML, Paulson JN, Sonawane AR, Fagny M, Platig J, Glass K, Quackenbush J, DeMeo DL. Sex differences in gene expression and regulatory networks across 29 human tissues. Cell Rep. 2020;31: 107795. https://doi.org/10.1016/j.celrep.2020.107795.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Oliva M, Muñoz-Aguirre M, Kim-Hellmuth S, Wucher V, Gewirtz ADH, Cotter DJ, Parsana P, Kasela S, Balliu B, Viñuela A, Castel SE, Mohammadi P, Aguet F, Zou Y, Khramtsova EA, Skol AD, Garrido-Martín D, Reverter F, Brown A, Evans P, Gamazon ER, Payne A, Bonazzola R, Barbeira AN, Hamel AR, Martinez-Perez A, Soria JM, Pierce BL, Stephens M, Eskin E, Dermitzakis ET, Segrè AV, Im HK, Engelhardt BE, Ardlie KG, Montgomery SB, Battle AJ, Lappalainen T, Guigó R, Stranger BE, GTEx Consortium. The impact of sex on gene expression across human tissues. Science. 2020. https://doi.org/10.1126/science.aba3066.

    Article  PubMed  PubMed Central  Google Scholar 

  14. Fletcher PA, Smiljanic K, Maso Prévide R, Iben JR, Li T, Rokic MB, Sherman A, Coon SL, Stojilkovic SS. Cell type- and sex-dependent transcriptome profiles of rat anterior pituitary cells. Front Endocrinol (Lausanne). 2019;10:623. https://doi.org/10.3389/fendo.2019.00623.

    Article  Google Scholar 

  15. Ho Y, Hu P, Peel MT, Chen S, Camara PG, Epstein DJ, Wu H, Liebhaber SA. Single-cell transcriptomic analysis of adult mouse pituitary reveals sexual dimorphism and physiologic demand-induced cellular plasticity. Protein Cell. 2020;11:565–83. https://doi.org/10.1007/s13238-020-00705-x.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Ruf-Zamojski F, Zhang Z, Zamojski M, Smith GR, Mendelev N, Liu H, Nudelman G, Moriwaki M, Pincas H, Castanon RG, Nair VD, Seenarine N, Amper MAS, Zhou X, Ongaro L, Toufaily C, Schang G, Nery JR, Bartlett A, Aldridge A, Jain N, Childs GV, Troyanskaya OG, Ecker JR, Turgeon JL, Welt CK, Bernard DJ, Sealfon SC. Single nucleus multi-omics regulatory landscape of the murine pituitary. Nat Commun. 2021;12:2677. https://doi.org/10.1038/s41467-021-22859-w.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Hou H, Uusküla-Reimand L, Makarem M, Corre C, Saleh S, Metcalf A, Goldenberg A, Palmert MR, Wilson MD. Gene expression profiling of puberty-associated genes reveals abundant tissue and sex-specific changes across postnatal development. Hum Mol Genet. 2017;26:3585–99. https://doi.org/10.1093/hmg/ddx246.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Qiao S, Nordström K, Muijs L, Gasparoni G, Tierling S, Krause E, Walter J, Boehm U. Molecular plasticity of male and female murine gonadotropes revealed by mRNA sequencing. Endocrinology. 2016;157:1082–93. https://doi.org/10.1210/en.2015-1836.

    Article  CAS  PubMed  Google Scholar 

  19. Bak M, Silahtaroglu A, Møller M, Christensen M, Rath MF, Skryabin B, Tommerup N, Kauppinen S. MicroRNA expression in the adult mouse central nervous system. RNA. 2008;14:432–44. https://doi.org/10.1261/rna.783108.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Bottoni A, Zatelli MC, Ferracin M, Tagliati F, Piccin D, Vignali C, Calin GA, Negrini M, Croce CM, Degli Uberti EC. Identification of differentially expressed microRNAs by microarray: a possible role for microRNA genes in pituitary adenomas. J Cell Physiol. 2007;210:370–7. https://doi.org/10.1002/jcp.20832.

    Article  CAS  PubMed  Google Scholar 

  21. Ye J, Yao Z, Si W, Gao X, Yang C, Liu Y, Ding J, Huang W, Fang F, Zhou J. Identification and characterization of microRNAs in the pituitary of pubescent goats. Reprod Biol Endocrinol. 2018;16:51. https://doi.org/10.1186/s12958-018-0370-x.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Ye R-S, Li M, Qi Q-E, Cheng X, Chen T, Li C-Y, Wang S-B, Shu G, Wang L-N, Zhu X-T, Jiang Q-Y, Xi Q-Y, Zhang Y-L. Comparative anterior pituitary miRNA and mRNA expression profiles of Bama Minipigs and Landrace Pigs reveal potential molecular network involved in animal postnatal growth. PLoS ONE. 2015;10: e0131987. https://doi.org/10.1371/journal.pone.0131987.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Yuan B, Han D-X, Dai L-S, Gao Y, Ding Y, Yu X-F, Chen J, Jiang H, Chen C-Z, Zhang J-B. A comprehensive expression profile of micrornas in rat’s pituitary. Int J Clin Exp Med. 2015;8:13289–95.

    CAS  PubMed  PubMed Central  Google Scholar 

  24. Zhang H, Qi Q, Chen T, Luo J, Xi Q, Jiang Q, Sun J, Zhang Y. Age-related changes in microRNA in the rat pituitary and potential role in GH regulation. Int J Mol Sci. 2018. https://doi.org/10.3390/ijms19072058.

    Article  PubMed  PubMed Central  Google Scholar 

  25. Hao P, Waxman DJ. Functional roles of sex-biased, growth hormone-regulated microRNAs miR-1948 and miR-802 in young adult mouse liver. Endocrinology. 2018;159:1377–92. https://doi.org/10.1210/en.2017-03109.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Morgan CP, Bale TL. Sex differences in microRNA-mRNA networks: examination of novel epigenetic programming mechanisms in the sexually dimorphic neonatal hypothalamus. Biol Sex Differ. 2017;8:27. https://doi.org/10.1186/s13293-017-0149-3.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Danilovich N, Wernsing D, Coschigano KT, Kopchick JJ, Bartke A. Deficits in female reproductive function in GH-R-KO mice; role of IGF-I. Endocrinology. 1999;140:2637–40. https://doi.org/10.1210/endo.140.6.6992.

    Article  CAS  PubMed  Google Scholar 

  28. Korenbrot CC, Huhtaniemi IT, Weiner RI. Preputial separation as an external sign of pubertal development in the male rat. Biol Reprod. 1977;17:298–303. https://doi.org/10.1095/biolreprod17.2.298.

    Article  CAS  PubMed  Google Scholar 

  29. Sánchez-Garrido MA, Castellano JM, Ruiz-Pino F, Garcia-Galiano D, Manfredi-Lozano M, Leon S, Romero-Ruiz A, Diéguez C, Pinilla L, Tena-Sempere M. Metabolic programming of puberty: sexually dimorphic responses to early nutritional challenges. Endocrinology. 2013;154:3387–400. https://doi.org/10.1210/en.2012-2157.

    Article  CAS  PubMed  Google Scholar 

  30. Yuki KE, Eyck TT, Bannister S, Kyriakopoulou L, Shlien A, Wilson MD. Automation of the Lexogen QuantSeq3’ mRNA Kit on the Agilent NGS workstation produces high-qualitysequencing libraries. Agilent Technologies. 2018.

  31. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, Batut P, Chaisson M, Gingeras TR. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29:15–21. https://doi.org/10.1093/bioinformatics/bts635.

    Article  CAS  PubMed  Google Scholar 

  32. Okonechnikov K, Conesa A, García-Alcalde F. Qualimap 2: advanced multi-sample quality control for high-throughput sequencing data. Bioinformatics. 2016;32:292–4. https://doi.org/10.1093/bioinformatics/btv566.

    Article  CAS  PubMed  Google Scholar 

  33. Collado-Torres L, Nellore A, Frazee AC, Wilks C, Love MI, Langmead B, Irizarry RA, Leek JT, Jaffe AE. Flexible expressed region analysis for RNA-seq with derfinder. Nucleic Acids Res. 2017;45: e9. https://doi.org/10.1093/nar/gkw852.

    Article  CAS  PubMed  Google Scholar 

  34. Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014;30:923–30. https://doi.org/10.1093/bioinformatics/btt656.

    Article  CAS  PubMed  Google Scholar 

  35. Bushnell B. BBMap: a fast, accurate, splice-aware aligner. SourceForge. 2014.

  36. Friedländer MR, Mackowiak SD, Li N, Chen W, Rajewsky N. miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 2012;40:37–52. https://doi.org/10.1093/nar/gkr688.

    Article  CAS  PubMed  Google Scholar 

  37. Kozomara A, Griffiths-Jones S. miRBase: annotating high confidence microRNAs using deep sequencing data. Nucleic Acids Res. 2014;42:D68-73. https://doi.org/10.1093/nar/gkt1181.

    Article  CAS  PubMed  Google Scholar 

  38. Risso D, Ngai J, Speed TP, Dudoit S. Normalization of RNA-seq data using factor analysis of control genes or samples. Nat Biotechnol. 2014;32:896–902. https://doi.org/10.1038/nbt.2931.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. McCarthy DJ, Chen Y, Smyth GK. Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation. Nucleic Acids Res. 2012;40:4288–97. https://doi.org/10.1093/nar/gks042.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26:139–40. https://doi.org/10.1093/bioinformatics/btp616.

    Article  CAS  PubMed  Google Scholar 

  41. Kim VN. MicroRNA biogenesis: coordinated cropping and dicing. Nat Rev Mol Cell Biol. 2005;6:376–85. https://doi.org/10.1038/nrm1644.

    Article  CAS  PubMed  Google Scholar 

  42. O’Brien J, Hayder H, Zayed Y, Peng C. Overview of microrna biogenesis, mechanisms of actions, and circulation. Front Endocrinol (Lausanne). 2018;9:402. https://doi.org/10.3389/fendo.2018.00402.

    Article  Google Scholar 

  43. Ruby JG, Jan CH, Bartel DP. Intronic microRNA precursors that bypass Drosha processing. Nature. 2007;448:83–6. https://doi.org/10.1038/nature05983.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Kozomara A, Birgaoanu M, Griffiths-Jones S. miRBase: from microRNA sequences to function. Nucleic Acids Res. 2019;47:D155–62. https://doi.org/10.1093/nar/gky1141.

    Article  CAS  PubMed  Google Scholar 

  45. Qin Q, Fan J, Zheng R, Wan C, Mei S, Wu Q, Sun H, Brown M, Zhang J, Meyer CA, Liu XS. Lisa: inferring transcriptional regulators through integrative modeling of public chromatin accessibility and ChIP-seq data. Genome Biol. 2020;21:32. https://doi.org/10.1186/s13059-020-1934-6.

    Article  PubMed  PubMed Central  Google Scholar 

  46. Russo PST, Ferreira GR, Cardozo LE, Bürger MC, Arias-Carrasco R, Maruyama SR, Hirata TDC, Lima DS, Passos FM, Fukutani KF, Lever M, Silva JS, Maracaja-Coutinho V, Nakaya HI. CEMiTool: a Bioconductor package for performing comprehensive modular co-expression analyses. BMC Bioinform. 2018;19:56. https://doi.org/10.1186/s12859-018-2053-1.

    Article  CAS  Google Scholar 

  47. Agarwal V, Bell GW, Nam J-W, Bartel DP. Predicting effective microRNA target sites in mammalian mRNAs. Elife. 2015. https://doi.org/10.7554/eLife.05005.

    Article  PubMed  PubMed Central  Google Scholar 

  48. Chou C-H, Shrestha S, Yang C-D, Chang N-W, Lin Y-L, Liao K-W, Huang W-C, Sun T-H, Tu S-J, Lee W-H, Chiew M-Y, Tai C-S, Wei T-Y, Tsai T-R, Huang H-T, Wang C-Y, Wu H-Y, Ho S-Y, Chen P-R, Chuang C-H, Hsieh P-J, Wu Y-S, Chen W-L, Li M-J, Wu Y-C, Huang X-Y, Ng FL, Buddhakosai W, Huang P-C, Lan K-C, Huang C-Y, Weng S-L, Cheng Y-N, Liang C, Hsu W-L, Huang H-D. miRTarBase update 2018: a resource for experimentally validated microRNA-target interactions. Nucleic Acids Res. 2018;46:D296–302. https://doi.org/10.1093/nar/gkx1067.

    Article  CAS  PubMed  Google Scholar 

  49. Hao Y, Hao S, Andersen-Nissen E, Mauck WM, Zheng S, Butler A, Lee MJ, Wilk AJ, Darby C, Zager M, Hoffman P, Stoeckius M, Papalexi E, Mimitou EP, Jain J, Srivastava A, Stuart T, Fleming LM, Yeung B, Rogers AJ, McElrath JM, Blish CA, Gottardo R, Smibert P, Satija R. Integrated analysis of multimodal single-cell data. Cell. 2021;184:3573-3587.e29. https://doi.org/10.1016/j.cell.2021.04.048.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Hafemeister C, Satija R. Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression. Genome Biol. 2019;20:296. https://doi.org/10.1186/s13059-019-1874-1.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  51. Stuart T, Butler A, Hoffman P, Hafemeister C, Papalexi E, Mauck WM, Hao Y, Stoeckius M, Smibert P, Satija R. comprehensive integration of single-cell data. Cell. 2019;177:1888-1902.e21. https://doi.org/10.1016/j.cell.2019.05.031.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Sokolowski DJ, Faykoo-Martinez M, Erdman L, Hou H, Chan C, Zhu H, Holmes MM, Goldenberg A, Wilson MD. Single-cell mapper (scMappR): using scRNA-seq to infer the cell-type specificities of differentially expressed genes. NAR Genom Bioinform. 2021;3:lqab011. https://doi.org/10.1093/nargab/lqab011.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Danziger SA, Gibbs DL, Shmulevich I, McConnell M, Trotter MWB, Schmitz F, Reiss DJ, Ratushny AV. ADAPTS: Automated deconvolution augmentation of profiles for tissue specific cells. PLoS ONE. 2019;14: e0224693. https://doi.org/10.1371/journal.pone.0224693.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Corre C, Shinoda G, Zhu H, Cousminer DL, Crossman C, Bellissimo C, Goldenberg A, Daley GQ, Palmert MR. Sex-specific regulation of weight and puberty by the Lin28/let-7 axis. J Endocrinol. 2016;228:179–91. https://doi.org/10.1530/JOE-15-0360.

    Article  CAS  PubMed  Google Scholar 

  55. Cheung LYM, George AS, McGee SR, Daly AZ, Brinkmeier ML, Ellsworth BS, Camper SA. Single-cell RNA sequencing reveals novel markers of male pituitary stem cells and hormone-producing cell types. Endocrinology. 2018;159:3910–24. https://doi.org/10.1210/en.2018-00750.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. Eckstrum KS, Weis KE, Baur NG, Yoshihara Y, Raetzman LT. Icam5 expression exhibits sex differences in the neonatal pituitary and is regulated by estradiol and bisphenol A. Endocrinology. 2016;157:1408–20. https://doi.org/10.1210/en.2015-1521.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  57. Cheung LYM, Rizzoti K, Lovell-Badge R, Le Tissier PR. Pituitary phenotypes of mice lacking the notch signalling ligand delta-like 1 homologue. J Neuroendocrinol. 2013;25:391–401. https://doi.org/10.1111/jne.12010.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  58. Nishida Y, Yoshioka M, St-Amand J. The top 10 most abundant transcripts are sufficient to characterize the organs functional specificity: evidences from the cortex, hypothalamus and pituitary gland. Gene. 2005;344:133–41. https://doi.org/10.1016/j.gene.2004.09.007.

    Article  CAS  PubMed  Google Scholar 

  59. Robinson AG, Verbalis JG. Posterior pituitary. In: Williams textbook of endocrinology. Elsevier; 2011. pp. 291–323. doi:https://doi.org/10.1016/B978-1-4377-0324-5.00010-9.

  60. Stojilkovic SS, Bjelobaba I, Zemkova H. Ion channels of pituitary gonadotrophs and their roles in signaling and secretion. Front Endocrinol (Lausanne). 2017;8:126. https://doi.org/10.3389/fendo.2017.00126.

    Article  Google Scholar 

  61. Stojilkovic SS, Tabak J, Bertram R. Ion channels and signaling in the pituitary gland. Endocr Rev. 2010;31:845–915. https://doi.org/10.1210/er.2010-0005.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. Stinnett GS, Westphal NJ, Seasholtz AF. Pituitary CRH-binding protein and stress in female mice. Physiol Behav. 2015;150:16–23. https://doi.org/10.1016/j.physbeh.2015.02.050.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  63. Torner L. Actions of prolactin in the brain: from physiological adaptations to stress and neurogenesis to psychopathology. Front Endocrinol (Lausanne). 2016;7:25. https://doi.org/10.3389/fendo.2016.00025.

    Article  Google Scholar 

  64. Li R, Vannitamby A, Yue SSK, Handelsman D, Hutson J. Mouse minipuberty coincides with gonocyte transformation into spermatogonial stem cells: a model for human minipuberty. Reprod Fertil Dev. 2017;29:2430–6. https://doi.org/10.1071/RD17100.

    Article  CAS  PubMed  Google Scholar 

  65. Schroeder A, Buret L, Hill RA, van den Buuse M. Gene-environment interaction of reelin and stress in cognitive behaviours in mice: implications for schizophrenia. Behav Brain Res. 2015;287:304–14. https://doi.org/10.1016/j.bbr.2015.03.063.

    Article  CAS  PubMed  Google Scholar 

  66. Hodge CW, Raber J, McMahon T, Walter H, Sanchez-Perez AM, Olive MF, Mehmert K, Morrow AL, Messing RO. Decreased anxiety-like behavior, reduced stress hormones, and neurosteroid supersensitivity in mice lacking protein kinase Cepsilon. J Clin Invest. 2002;110:1003–10. https://doi.org/10.1172/JCI15903.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  67. Heck AL, Handa RJ. Sex differences in the hypothalamic-pituitary-adrenal axis’ response to stress: an important role for gonadal hormones. Neuropsychopharmacology. 2019;44:45–58. https://doi.org/10.1038/s41386-018-0167-9.

    Article  CAS  PubMed  Google Scholar 

  68. Terenina EE, Cavigelli S, Mormede P, Zhao W, Parks C, Lu L, Jones BC, Mulligan MK. Genetic factors mediate the impact of chronic stress and subsequent response to novel acute stress. Front Neurosci. 2019;13:438. https://doi.org/10.3389/fnins.2019.00438.

    Article  PubMed  PubMed Central  Google Scholar 

  69. Shin S, Kwon O, Kang JI, Kwon S, Oh S, Choi J, Kim CH, Kim DG. mGluR5 in the nucleus accumbens is critical for promoting resilience to chronic stress. Nat Neurosci. 2015;18:1017–24. https://doi.org/10.1038/nn.4028.

    Article  CAS  PubMed  Google Scholar 

  70. Moon AL, Haan N, Wilkinson LS, Thomas KL, Hall J. CACNA1C: association with psychiatric disorders, behavior, and neurogenesis. Schizophr Bull. 2018;44:958–65. https://doi.org/10.1093/schbul/sby096.

    Article  PubMed  PubMed Central  Google Scholar 

  71. Oyola MG, Handa RJ. Hypothalamic-pituitary-adrenal and hypothalamic-pituitary-gonadal axes: sex differences in regulation of stress responsivity. Stress. 2017;20:476–94. https://doi.org/10.1080/10253890.2017.1369523.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  72. Hu H, Miao Y-R, Jia L-H, Yu Q-Y, Zhang Q, Guo A-Y. AnimalTFDB 3.0: a comprehensive resource for annotation and prediction of animal transcription factors. Nucleic Acids Res. 2019;47:D33–8. https://doi.org/10.1093/nar/gky822.

    Article  CAS  PubMed  Google Scholar 

  73. Chen Q, Leshkowitz D, Blechman J, Levkowitz G. Single-cell molecular and cellular architecture of the mouse neurohypophysis. eNeuro. 2020. https://doi.org/10.1523/ENEURO.0345-19.2019.

    Article  PubMed  PubMed Central  Google Scholar 

  74. Sasaki F, Iwama Y. Sex difference in prolactin and growth hormone cells in mouse adenohypophysis: stereological, morphometric, and immunohistochemical studies by light and electron microscopy. Endocrinology. 1988;123:905–12. https://doi.org/10.1210/endo-123-2-905.

    Article  CAS  PubMed  Google Scholar 

  75. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559. https://doi.org/10.1186/1471-2105-9-559.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  76. McCarthy MM, Arnold AP. Reframing sexual differentiation of the brain. Nat Neurosci. 2011;14:677–83. https://doi.org/10.1038/nn.2834.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  77. Million Passe CM, White CR, King MW, Quirk PL, Iovanna JL, Quirk CC. Loss of the protein NUPR1 (p8) leads to delayed LHB expression, delayed ovarian maturation, and testicular development of a sertoli-cell-only syndrome-like phenotype in mice. Biol Reprod. 2008;79:598–607. https://doi.org/10.1095/biolreprod.108.068304.

    Article  CAS  PubMed  Google Scholar 

  78. O’Hara L, Curley M, Tedim Ferreira M, Cruickshanks L, Milne L, Smith LB. Pituitary androgen receptor signalling regulates prolactin but not gonadotrophins in the male mouse. PLoS ONE. 2015;10: e0121657. https://doi.org/10.1371/journal.pone.0121657.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  79. Wu S, Chen Y, Fajobi T, DiVall SA, Chang C, Yeh S, Wolfe A. Conditional knockout of the androgen receptor in gonadotropes reveals crucial roles for androgen in gonadotropin synthesis and surge in female mice. Mol Endocrinol. 2014;28:1670–81. https://doi.org/10.1210/me.2014-1154.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  80. AlOgayil N, Bauermeister K, Galvez JH, Venkatesh VS, Zhuang QK-W, Chang ML, Davey RA, Zajac JD, Ida K, Kamiya A, Taketo T, Bourque G, Naumova AK. Distinct roles of androgen receptor, estrogen receptor alpha, and BCL6 in the establishment of sex-biased DNA methylation in mouse liver. Sci Rep. 2021;11:13766. https://doi.org/10.1038/s41598-021-93216-6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  81. Lau-Corona D, Bae WK, Hennighausen L, Waxman DJ. Sex-biased genetic programs in liver metabolism and liver fibrosis are controlled by EZH1 and EZH2. PLoS Genet. 2020;16: e1008796. https://doi.org/10.1371/journal.pgen.1008796.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  82. Lomniczi A, Loche A, Castellano JM, Ronnekleiv OK, Bosch M, Kaidar G, Knoll JG, Wright H, Pfeifer GP, Ojeda SR. Epigenetic control of female puberty. Nat Neurosci. 2013;16:281–9. https://doi.org/10.1038/nn.3319.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  83. Dina OA, Aley KO, Isenberg W, Messing RO, Levine JD. Sex hormones regulate the contribution of PKCε and PKA signalling in inflammatory pain in the rat. Eur J Neurosci. 2001;13:2227–33. https://doi.org/10.1046/j.0953-816x.2001.01614.x.

    Article  CAS  PubMed  Google Scholar 

  84. Kim HJ, Gieske MC, Trudgen KL, Hudgins-Spivey S, Kim BG, Krust A, Chambon P, Jeong J-W, Blalock E, Ko C. Identification of estradiol/ERα-regulated genes in the mouse pituitary. J Endocrinol. 2011;210:309–21. https://doi.org/10.1530/JOE-11-0098.

    Article  CAS  PubMed  Google Scholar 

  85. Alim Z, Hartshorn C, Mai O, Stitt I, Clay C, Tobet S, Boehm U. Gonadotrope plasticity at cellular and population levels. Endocrinology. 2012;153:4729–39. https://doi.org/10.1210/en.2012-1360.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  86. González-Parra S, Argente J, García-Segura LM, Chowen JA. Cellular composition of the adult rat anterior pituitary is influenced by the neonatal sex steroid environment. Neuroendocrinology. 1998;68:152–62. https://doi.org/10.1159/000054361.

    Article  PubMed  Google Scholar 

  87. Gonzázalez-Parra S, Argente J, García-Segura LM, Chowen JA. Effect of neonatal and adult testosterone treatment on the cellular composition of the adult female rat anterior pituitary. J Endocrinol. 2000;164:265–76. https://doi.org/10.1677/joe.0.1640265.

    Article  PubMed  Google Scholar 

  88. Andoniadou CL, Matsushima D, Mousavy Gharavy SN, Signore M, Mackintosh AI, Schaeffer M, Gaston-Massuet C, Mollard P, Jacques TS, Le Tissier P, Dattani MT, Pevny LH, Martinez-Barbera JP. Sox2(+) stem/progenitor cells in the adult mouse pituitary support organ homeostasis and have tumor-inducing potential. Cell Stem Cell. 2013;13:433–45. https://doi.org/10.1016/j.stem.2013.07.004.

    Article  CAS  PubMed  Google Scholar 

  89. Scheithauer BW, Sano T, Kovacs KT, Young WF, Ryan N, Randall RV. The pituitary gland in pregnancy: a clinicopathologic and immunohistochemical study of 69 cases. Mayo Clin Proc. 1990;65:461–74. https://doi.org/10.1016/s0025-6196(12)60946-x.

    Article  CAS  PubMed  Google Scholar 

  90. Oishi Y, Okuda M, Takahashi H, Fujii T, Morii S. Cellular proliferation in the anterior pituitary gland of normal adult rats: influences of sex, estrous cycle, and circadian change. Anat Rec. 1993;235:111–20. https://doi.org/10.1002/ar.1092350111.

    Article  CAS  PubMed  Google Scholar 

  91. Taniguchi Y, Yasutaka S, Kominami R, Shinohara H. Proliferation and differentiation of rat anterior pituitary cells. Anat Embryol (Berl). 2002;206:1–11. https://doi.org/10.1007/s00429-002-0271-8.

    Article  CAS  Google Scholar 

  92. Perez F, Lledo PM, Karagogeos D, Vincent JD, Prochiantz A, Ayala J. Rab3A and Rab3B carboxy-terminal peptides are both potent and specific inhibitors of prolactin release by rat cultured anterior pituitary cells. Mol Endocrinol. 1994;8:1278–87. https://doi.org/10.1210/mend.8.9.7838160.

    Article  CAS  PubMed  Google Scholar 

  93. Charles MA, Saunders TL, Wood WM, Owens K, Parlow AF, Camper SA, Ridgway EC, Gordon DF. Pituitary-specific Gata2 knockout: effects on gonadotrope and thyrotrope function. Mol Endocrinol. 2006;20:1366–77. https://doi.org/10.1210/me.2005-0378.

    Article  CAS  PubMed  Google Scholar 

  94. Schang G, Ongaro L, Brûlé E, Zhou X, Wang Y, Boehm U, Ruf-Zamojski F, Zamojski M, Mendelev N, Seenarine N, Amper MA, Nair V, Ge Y, Sealfon SC, Bernard DJ. Transcription factor GATA2 may potentiate follicle-stimulating hormone production in mice via induction of the BMP antagonist gremlin in gonadotrope cells. J Biol Chem. 2022;298: 102072. https://doi.org/10.1016/j.jbc.2022.102072.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  95. Ruf-Zamojski F, Fribourg M, Ge Y, Nair V, Pincas H, Zaslavsky E, Nudelman G, Tuminello SJ, Watanabe H, Turgeon JL, Sealfon SC. Regulatory architecture of the LβT2 gonadotrope cell underlying the response to gonadotropin-releasing hormone. Front Endocrinol (Lausanne). 2018;9:34. https://doi.org/10.3389/fendo.2018.00034.

    Article  Google Scholar 

  96. Hücker SM, Fehlmann T, Werno C, Weidele K, Lüke F, Schlenska-Lange A, Klein CA, Keller A, Kirsch S. Single-cell microRNA sequencing method comparison and application to cell lines and circulating lung tumor cells. Nat Commun. 2021;12:4316. https://doi.org/10.1038/s41467-021-24611-w.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  97. Wang N, Zheng J, Chen Z, Liu Y, Dura B, Kwak M, Xavier-Ferrucio J, Lu Y-C, Zhang M, Roden C, Cheng J, Krause DS, Ding Y, Fan R, Lu J. Single-cell microRNA-mRNA co-sequencing reveals non-genetic heterogeneity and mechanisms of microRNA regulation. Nat Commun. 2019;10:95. https://doi.org/10.1038/s41467-018-07981-6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  98. Tong Y, Zhou J, Mizutani J, Fukuoka H, Ren S-G, Gutierrez-Hartmann A, Koeffler HP, Melmed S. CEBPD suppresses prolactin expression and prolactinoma cell proliferation. Mol Endocrinol. 2011;25:1880–91. https://doi.org/10.1210/me.2011-1075.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  99. Pyczek J, Buslei R, Schult D, Hölsken A, Buchfelder M, Heß I, Hahn H, Uhmann A. Hedgehog signaling activation induces stem cell proliferation and hormone release in the adult pituitary gland. Sci Rep. 2016;6:24928. https://doi.org/10.1038/srep24928.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  100. Vila G, Papazoglou M, Stalla J, Theodoropoulou M, Stalla GK, Holsboer F, Paez-Pereda M. Sonic hedgehog regulates CRH signal transduction in the adult pituitary. FASEB J. 2005;19:281–3. https://doi.org/10.1096/fj.04-2138fje.

    Article  CAS  PubMed  Google Scholar 

  101. Kober P, Boresowicz J, Rusetska N, Maksymowicz M, Paziewska A, Dąbrowska M, Kunicki J, Bonicki W, Ostrowski J, Siedlecki JA, Bujko M. The role of aberrant DNA methylation in misregulation of gene expression in gonadotroph nonfunctioning pituitary tumors. Cancers (Basel). 2019. https://doi.org/10.3390/cancers11111650.

    Article  Google Scholar 

  102. Miyamoto J, Matsumoto T, Shiina H, Inoue K, Takada I, Ito S, Itoh J, Minematsu T, Sato T, Yanase T, Nawata H, Osamura YR, Kato S. The pituitary function of androgen receptor constitutes a glucocorticoid production circuit. Mol Cell Biol. 2007;27:4807–14. https://doi.org/10.1128/MCB.02039-06.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  103. García IA, Torres Demichelis V, Viale DL, Di Giusto P, Ezhova Y, Polishchuk RS, Sampieri L, Martinez H, Sztul E, Alvarez C. CREB3L1-mediated functional and structural adaptation of the secretory pathway in hormone-stimulated thyroid cells. J Cell Sci. 2017;130:4155–67. https://doi.org/10.1242/jcs.211102.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  104. Greenwood M, Paterson A, Rahman PA, Gillard BT, Langley S, Iwasaki Y, Murphy D, Greenwood MP. Transcription factor Creb3l1 regulates the synthesis of prohormone convertase enzyme PC1/3 in endocrine cells. J Neuroendocrinol. 2020;32: e12851. https://doi.org/10.1111/jne.12851.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  105. Konishi H, Ogawa T, Nakagomi S, Inoue K, Tohyama M, Kiyama H. Id1, Id2 and Id3 are induced in rat melanotrophs of the pituitary gland by dopamine suppression under continuous stress. Neuroscience. 2010;169:1527–34. https://doi.org/10.1016/j.neuroscience.2010.06.030.

    Article  CAS  PubMed  Google Scholar 

  106. Zhu X, Zhang J, Tollkuhn J, Ohsawa R, Bresnick EH, Guillemot F, Kageyama R, Rosenfeld MG. Sustained Notch signaling in progenitors is required for sequential emergence of distinct cell lineages during organogenesis. Genes Dev. 2006;20:2739–53. https://doi.org/10.1101/gad.1444706.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  107. Cogliati T, Delgado-Romero P, Norwitz ER, Guduric-Fuchs J, Kaiser UB, Wray S, Kirsch IR. Pubertal impairment in Nhlh2 null mice is associated with hypothalamic and pituitary deficiencies. Mol Endocrinol. 2007;21:3013–27. https://doi.org/10.1210/me.2005-0337.

    Article  CAS  PubMed  Google Scholar 

  108. Gordon A, Garrido-Gracia JC, Aguilar R, Sánchez-Criado JE. Understanding the regulation of pituitary progesterone receptor expression and phosphorylation. Reproduction. 2015;149:615–23. https://doi.org/10.1530/REP-14-0592.

    Article  CAS  PubMed  Google Scholar 

  109. Turgeon JL, Waring DW. Progesterone regulation of the progesterone receptor in rat gonadotropes. Endocrinology. 2000;141:3422–9. https://doi.org/10.1210/endo.141.9.7688.

    Article  CAS  PubMed  Google Scholar 

  110. Li J-T, Xie X-M, Yu J-Y, Sun Y-X, Liao X-M, Wang X-X, Su Y-A, Liu Y-J, Schmidt MV, Wang X-D, Si T-M. Suppressed calbindin levels in hippocampal excitatory neurons mediate stress-induced memory loss. Cell Rep. 2017;21:891–900. https://doi.org/10.1016/j.celrep.2017.10.006.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

We would like to thank: Theresa Ten Eyck (Agilent), Jekaterina Aleksejeva and Stephanie Bannister (Lexogen), Karen Ho and Sergio Pereira (The Centre for Applied Genomics (TCAG)) for their efforts and guidance setting up the automation of the QuantSeq protocol; and Sunyun Lee for assistance with setting up the Shiny app.

Funding

This work was supported by CIHR grants: 312557 (MRP/MDW/AG) and 437197 (Melissa Holmes/MDW/MRP). MDW is supported by the Canada Research Chairs Program. RQ, C.Chan and DS were supported in part by NSERC grant RGPIN-2019-07014 to MDW. C.Chan and MH were supported by a SickKids RESTRACOMP scholarship. DS is supported by NSERC CGS M, PGS D and Ontario Graduate Scholarships. HH is supported by the Genome Canada Genomics Technology Platform, The Centre for Applied Genomics. MFM is supported by NSERC PGS D and the association computing machinery special interest group on high performance computing (ACM/SIGHPC) Intel Computational and Data Science Fellowship. LU was supported by the CRS Scholarships for the Next Generation of Scientists.

Author information

Author notes

  1. Huayun Hou and Cadia Chan contributed equally

Authors and Affiliations

  1. Genetics and Genome Biology, SickKids Research Institute, Toronto, ON, Canada

    Huayun Hou, Cadia Chan, Kyoko E. Yuki, Dustin Sokolowski, Anna Roy, Liis Uusküla-Reimand, Mariela Faykoo-Martinez, Matt Hudson, Christina Corre, Anna Goldenberg, Mark R. Palmert & Michael D. Wilson

  2. Department of Molecular Genetics, University of Toronto, Toronto, ON, Canada

    Cadia Chan, Dustin Sokolowski, Matt Hudson, Zhaolei Zhang & Michael D. Wilson

  3. Donnelly Centre for Cellular and Biomolecular Research, Toronto, ON, Canada

    Cadia Chan

  4. Interdepartmental Program of Computational Biology and Bioinformatics, Yale University, New Haven, CT, USA

    Rihao Qu

  5. Department of Pathology, Yale School of Medicine, New Haven, CT, USA

    Rihao Qu

  6. Department of Cell and Systems Biology, University of Toronto, Toronto, ON, Canada

    Mariela Faykoo-Martinez

  7. Department of Computer Science, University of Toronto, Toronto, ON, Canada

    Anna Goldenberg & Zhaolei Zhang

  8. Division of Endocrinology, The Hospital for Sick Children, Toronto, ON, Canada

    Christina Corre & Mark R. Palmert

  9. Institute of Medical Science, University of Toronto, Toronto, ON, Canada

    Mark R. Palmert

  10. Departments of Pediatrics and Physiology, University of Toronto, Toronto, ON, Canada

    Christina Corre & Mark R. Palmert

Authors
  1. Huayun Hou
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Cadia Chan
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Kyoko E. Yuki
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Dustin Sokolowski
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Anna Roy
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Rihao Qu
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Liis Uusküla-Reimand
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Mariela Faykoo-Martinez
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Matt Hudson
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. Christina Corre
    View author publications

    You can also search for this author in PubMed Google Scholar

  11. Anna Goldenberg
    View author publications

    You can also search for this author in PubMed Google Scholar

  12. Zhaolei Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

  13. Mark R. Palmert
    View author publications

    You can also search for this author in PubMed Google Scholar

  14. Michael D. Wilson
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

Experimental and analysis study design was conceived by HH, C.Chan, DS, MFM, AG, ZZ, MRP, and MDW. Bioinformatic analysis was performed by HH, C.Chan and DS. Experiments were performed by KEY, AR, LUR, MH, C.Corre. Pipeline for preprocessing 3’UTR-seq data was developed by RQ and HH. The manuscript was written by HH, C.Chan, and MDW with support from all authors. AG, ZZ, MRP and MDW supervised and obtained funding for this work. All authors have read and approved the manuscript for publication.

Corresponding authors

Correspondence to Mark R. Palmert or Michael D. Wilson.

Ethics declarations

Ethics approval and consent to participate

All studies and procedures were approved by the Toronto Centre for Phenogenomics (TCP) Animal Care Committee (AUP 09-08-0097) in accordance with recommendations of the Canadian Council on Animal Care, the requirements under Animals for Research Act, RSO 1980, and the TCP Committee Policies and Guidelines.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Figure S1.

Illustration of 3’UTR-seq method and quality control. A. Overview of the QuantSeq data analysis pipeline. B. Genome browser screenshot of extended 3’UTRs for genes Pou1f1 (left) and Ghrhr (right). Gene name and gene model are shown on the bottom of each panel. Each track represents overlapping signal from 5-6 biological replicates. C. Summary of qPCR vs. 3’UTR-seq comparisons across samples in all samples. Bottom panel: bar plots showing the numbers of genes detected using qPCR (white bars) and 3’UTR seq (red bars). The Spearman correlation coefficients between two experiments are shown for each sample in the dot plot (middle panel) and as well as in a density plot (top panel). Figure S2. Correlation heatmaps for pituitary gland samples. Pearson correlation between samples is calculated based on (A) gene expression and (B) miRNA expression. Hierarchical clustering is performed based on pairwise Pearson correlation coefficients (PCCs). Heatmap color intensity represents PCCs (orange: lower correlation; purple: higher correlation). Sample conditions are shown in top bars. Green shades: different ages, blue: male samples; red: female samples. Figure S3. Characterization of sex-biased mRNAs and miRNAs. A. Gene expression heatmap of sex-biased genes at PD12 and PD22, and genes with significant sex-by-age effect between PD12 and PD22. Each row represents a gene and each column represents a sample. Column annotation bars indicate sample age and sex. Colors represent row-scaled log2 (normalized counts). Whether a gene is female-biased (red) or male-biased (blue) at each corresponding age is summarized by the row annotation bars on the left. Sex chromosome-linked genes are labeled with asterisks. Genes with significant sex-by-age interaction effect between PD12 and PD22 are bolded. B. Expression plots of all sex-biased miRNAs. Log2(normCounts) are plotted for each miRNA across ages. Large filled points represent median expression at each age and unfilled points represent each biological replicate. Red: female samples; blue: male samples. C. novel46 is a mirtron of Cacna1g. The precursor of novel46 is expressed from the last intron of Cacna1g (highlighted in blue). The canonical seed region (nucleotide positions 2-8) in the mature sequence of novel46 is bolded. Figure S4. Summary of gene co-expression modules. Left panel. Module expression profile (log2-transformed normalized counts (log2(normCounts)) scaled per gene across all samples) is plotted for genes within each co-expression module across profiled postnatal ages. Dotted lines: individual gene profiles, solid line: median expression profile of module genes. Number of genes within each module is labeled at the top of the plots. Red: female samples; blue: male samples. Right panel. Expression profile of top five hub genes for each module. log2(normCounts) is plotted across profiled postnatal ages. Large, filled points represent median expression at each age and unfilled points represent each biological replicate. Blue: male samples; red: female samples. Figure S5. Characterization of co-expression gene modules. A. Correlation heatmap based on module eigengene (first PC, obtained using “mod_summary()” function from “CEMitool”). Color scale represents pairwise Pearson correlation coefficients. B. Barplots showing the eigengene for each sample in each module. Samples are grouped by age and colored by sex (blue: male; red: female). C. Barplots showing pathway enrichment results for genes in each module. Each bar represents a pathway. Color represents pathway categories (BP: Biological Process; CC: Cellular Component; MF: Molecular Function). Number of module genes in the pathway is labeled. X axis: -log10(P-value) of each pathway. Figure S6. Co-expression module gene enrichment in single-nuclei RNA-seq data. Uniform Manifold Approximation and Projection (UMAP) dimension reduction representation of adult female and male mouse pituitary gland single-nuclei transcriptome integrated by sex from Ruf-Zamojski et al. 2021 (A) colored by cell type and (B) colored by sample sex. Samples derived from snap-frozen pituitaries were first merged between replicates for each sex (n=3/sex). Merged samples were then integrated between sexes. C. Enrichment heatmap of co-expression module genes within cell types. One-sided Kolmogorov–Smirnov (KS) test was performed to test for enrichment of each co-expression module gene within a given cell type compared to all other cell types. Color gradient represents KS test FDR-adjusted P-value for each gene (dark purple: high enrichment; gray: low enrichment). Only genes with FDR ≤ 0.05 in at least one cell type are plotted and only FDR < 0.05 are shown. Each column represents a cell type as labeled at the bottom of the heatmap. Each row represents a co-expression module gene and the genes are grouped by the co-expression module in which the gene was identified (labeled on the right). Breaks were added in the heatmap between co-expression modules. Colored boxes represent the level of significance based on a one-sided hypergeometric test to determine if a group of module genes with KS test FDR ≤ 0.05 was significantly enriched for a given cell type. D. Estimated cell-type proportions by RNA-seq deconvolution using Proportions in Admixture (WGCNA) changes across profiled ages of pituitary cell types without known sex differences in their proportions. Estimated cell-type proportions are plotted across postnatal ages along the x-axis. Large circles and triangles represent the mean cell-type proportion at each age and small circles and triangles represent each biological replicate. Lighter color, solid line, circle points: female samples; dark color, dotted line, triangle points: male samples. Wilcoxon test was performed to compare cell proportions between both sexes at each age (*P<0.05, **P<0.01). See Figure 5B for estimated proportions of cell types with known sex biases.

Additional file 2: Table S1.

Quality control metrics of 3’UTR-seq and small RNA-seq libraries.

Additional file 3: Table S2.

Novel miRNAs and sequences.

Additional file 4: Table S3.

Differential analysis results for gene expression.

Additional file 5: Table S4.

Differential analysis results for miRNA expression.

Additional file 6: Table S5.

gProfiler pathway enrichment of sex-biased genes.

Additional file 7: Table S6.

Correlation of miRNA–gene targets.

Additional file 8: Table S7.

Genes in co-expression modules.

Additional file 9: Table S8.

Sex-biased genes identified by scMappR in cell types.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Hou, H., Chan, C., Yuki, K.E. et al. Postnatal developmental trajectory of sex-biased gene expression in the mouse pituitary gland. Biol Sex Differ 13, 57 (2022). https://doi.org/10.1186/s13293-022-00467-7

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13293-022-00467-7

Biology of Sex Differences

ISSN: 2042-6410

Contact us