ScRNA-seq analysis of human pituitary development

We obtained pituitaries from 21 human fetuses from 7 to 25 weeks postfertilization (11 females and 10 males) and performed a modified STRT-seq method on a total of 5181 cells, with 4113 high-quality cells being retained after filtration (Fig. 1a and Supplementary Fig. 1a). An average of 4506 genes and 86,497 transcripts (counted as unique molecular identifiers, UMIs) were detected in each cell (Supplementary Fig. 1c). The samples were detected with similar gene numbers and GAPDH expression across batches (Supplementary Fig. 1b, c). The morphology of the pituitary was verified (Supplementary Fig. 1g).

a Experimental flowchart for the scRNA-seq analysis of the human fetal pituitary. b UMAP plots showing the clusters of the cell cycle cells (CC), stem cells, the progenitor cells of PIT1 lineage (Pro.PIT1) or precursor cells of gonadotrope (Pre.Gonado) and the terminal endocrine cells (lower), and distribution of the fetal samples (upper). Dots: single cells. c Scatterplots showing expression of known markers, including TFs and hormone genes, projecting on the UMAP plot (b). Gray to red indicates no expression to high expression levels. d Bar plots showing the proportions of each cell type in each stage. Solid circles at the bottom indicate the earliest stages when a hormone producing cell type appears. e Heatmap showing the activated TFs predicated by SCENIC. For each cell type, the top 10 -log(P value) specific TFs being activated are shown, which ranked by number of cells. Columns are individual cells and rows are individual genes. White: not activated; red: activated. The expression levels of these TFs are shown in Supplementary Fig. 2c.

We used Seurat to identify cell clusters and Uniform Manifold Approximation and Projection (UMAP) for visualization (Fig. 1b and Supplementary Fig. 1d)¹². A total of 14 cell clusters identified with known marker genes (Fig. 1b, Supplementary Fig. 1d). We identified nine clusters of anterior pituitary endocrine cells including the stem cells (Stem), cycling cells (CC), corticotropes, progenitors of the PIT-1 lineage (Pro.PIT1), somatotropes, lactotropes, thyrotropes, precursors of gonadotropes (Pre.Gonado) and gonadotropes, comprising 2,388 cells (Fig. 1b and Supplementary Fig. 1a). PITX1 and PITX2 were expressed in all nine clusters, and SOX2 and PROP1 were specifically expressed in the stem cells (Fig. 1c and Supplementary Fig. 1f). Lineage-specific TFs, such as POU1F1, TBX19 and NR5A1, and hormone genes, including POMC, GH1, GH2, TSHB, PRL, FSHB and LHB, were expressed in special clusters of the endocrine cells. Mesenchymal cells, as marked by PDGFRA¹³, were the second most abundant cell type after endocrine cells, comprising 1,005 cells. A cluster of 90 cells was identified as posterior pituitary pituicyte cells (P) as marked by OTX2, LHX2, RAX and COL25A1^11,14 (Supplementary Fig. 1f). Other clusters included endothelial cells (PECAM1⁺), immune cells (IMM, PTPRC⁺) and red blood cells (RBC, HBQ1⁺). Each cell cluster was composed of multiple fetal samples, and the samples of similar stages, or different sexes, were largely mixedly distributed (Fig. 1b, Supplementary Fig. 1d, e and 2a).

We analyzed the timing of endocrine cell differentiation. The results showed that the corticotropes appeared first at 7 weeks, the earliest time point we analyzed; this was immediately followed by the somatotropes and the gonadotropes at 8 weeks, and the thyrotropes and lactotropes appeared later at approximately 10 weeks and 16 weeks, respectively (Fig. 1d). The timing of endocrine cell differentiation was validated by immunofluorescence staining for the hormone genes, which was consistent with previous studies (Supplementary Fig. 2b)¹⁵.

We analyzed a number of TFs that are mutated in human pituitary genetic diseases¹⁶, showing that these TFs were expressed in a cell type-specific manner (Supplementary Fig. 1h). The SCENIC analysis identified activation of known TFs including SOX2, TBX19, POU1F1 and NR5A1 (Fig. 1e and Supplementary Fig. 2c)¹⁷. Taken together, these results indicated that our data provided comprehensive and precise information on human pituitary development.

Characterizing pituitary stem cells

Human pituitary stem cells have not been comprehensively characterized. We identified stem cell-specific genes including SOX2, PROP1, LHX3, HES1, ZFP36L1, ANXA1, NFIB, ZNF521 and NR2F2 (Fig. 2a). Gene set enrichment analysis (GSEA) showed that the TGF-β, Notch, Wnt and Hedgehog signaling pathways were enriched in the stem cells, and the “tight junction”, “cell cycle” and “ECM receptor interaction” pathways were also highly enriched (Fig. 2b). We identified potential ligand-receptor genes between stem cells and mesenchymal cells, which enriched Gene Ontology (GO) terms such as “ECM receptor interaction” and the Notch, Wnt, BMP and Eph signaling pathways, suggesting that the mesenchymal cells may involve in regulation of the stem cells (Supplementary Fig. 3a, b and Supplementary Data 1). Co-immunofluorescence staining for SOX2 and Collagen III (COL3) suggested that mesenchymal cells encompass stem cells in human fetal pituitaries (Supplementary Fig. 3c).

a DEGs between the stem cells and the differentiated endocrine cells shown by a z-scored heatmap. DEC: differentiated endocrine cells, which include the progenitor or precursor cells. b Bar plots showing differentially enriched KEGG pathways between the stem cells and the DEC detected by GSEA. Representative genes in each pathway are shown. Nominal P values are determined by two-sided Kolmogorov-Smirnov Test and adjusted by FDR. c Bar plots showing the stages of three subtypes of the stem cells. The Stem1 are mainly derived from the early stages (7 to 10 weeks), and the Stem3 cells are mainly derived from the late stages (19 to 25 weeks). Color: weeks postfertilization. d Violin plots of representative DEGs between the Stem1 and Stem3 cells. e Immunofluorescence staining for SOX2 and ASCL1 in the 8- and 23-week human fetal pituitaries. Triangles indicate representative cells coexpressing both genes. Scale bar, 50 μm. f The CC cluster are composed of different cell subtypes. The upper panel shows the average expression levels of the S and G2/M phase genes, and representative cell type markers. CC1 and CC2 are SOX2-positive stem cells with CC1 and CC2 being the S and G2/M phase, respectively. CC3 are non-stem cells expressing POU1F1, TBX19 or GATA2. The lower panel shows the proportions of the cycling to all stem cells as a function of weeks postfertilization. Data are presented as mean ± SEM. g Immunofluorescence staining of SOX2 and Ki67 in the 8-, 17- and 23-week human fetal pituitaries (left). Triangles indicate representative costaining cells. The right panel shows the ratios of double-positive and SOX2-positive cells in each sample (n ≥ 2). p-Value is determined by two-sided Wilcoxon test (*P = 0.023 which is less than 0.05, **P = 0.0024 which is less than 0.01). Data are presented as mean ± SEM. Scale bar, 50 μm.

Notably, reclustering the stem cells by Seurat revealed three subpopulations (Stem1, Stem2 and Stem3, Supplement Fig. 3d). Interestingly, Stem1 cells were found to be mainly derived from the early-stage pituitaries (7 to 10 weeks), while Stem3 cells were mainly derived from the late-stage pituitaries (19 to 25 weeks), indicating that these subpopulations represented time-dependent cell state shifts in the stem cells (Fig. 2c). The clustering result remained similar after regressing out the cell cycle genes (Supplement Fig. 3d).

Examination of differentially expressed genes (DEGs) among the three subpopulations identified 114 and 165 genes that were highly expressed in Stem1 and Stem3, respectively (Supplementary Fig. 3e). Stem1 cells were enriched in genes for “stem cell proliferation” (e.g. HMGA2), while the Stem3 cells were strongly enriched in genes for “negative regulation of cell proliferation” (e.g., CDKN1A, Fig. 2d and Supplementary Fig. 3e). Reclustering of the CC cluster showed that most cycling cells were SOX2-positive stem cells, and a small portion were POU1F1 or TBX19-positive cells, which were verfied by immunofluorescence staining (Fig. 2f and Supplementary Fig. 3f). The proportion of cycling stem cells in the total stem cell population gradually decreased from the early to late stages (20% in the early stages, 10% in the middle stages, and less than 5% in the late stages; Fig. 2f), which was validated by coimmunostaining for SOX2 and Ki67 in the pituitaries of 8-, 17- and 23-week fetuses (Fig. 2g). These results together suggest that Stem3 cells enter a quiescent or lowly proliferative state.

Interestingly, ASCL1 was specifically expressed in Stem1 (Fig. 2d). This expression was consistent with the expression of Ascl1 in the early stage of pituitary development in mouse and zebrafish^18,19. We performed immunofluorescence staining to show that a large portion (24%) of SOX2-positive stem cells were also positive for ASCL1 at 10 weeks, while there were nearly no double-positive cells at 23 weeks (Fig. 2e). ASCL1 is required for differentiation of all anterior pituitary endocrine cell types in zebrafish, while in mice, mutation of Ascl1 affects corticotropes and gonadotropes, and maybe also thyrotropes, at different regulatory steps^19,20,21. The downregulation of ASCL1 expression in the late pituitary stem cells may play a role in the state shift. The Wnt signaling modulator SFRP2 and the non-canonical WNT gene WNT5A was prominently upregulated in Stem2 and Stem3 (Fig. 2d), which was consistent with the enrichment of the GO term “negative regulation of canonical Wnt signaling pathway” in Stem3 (Supplementary Fig. 3c).

The marker of the folliculostellate cells (FSCs) S100B was expressed in a small fraction (3.4%) of the stem cells of both the early and late stages (Fig. 1c). Comparison between the S100B-postive and S100B-negative stem cells revealed no DEGs relating to FSCs. The results suggested that FSCs were not present in the human fetal pituitary during 25 weeks postfertilization, which was consistent with previous immunostaining studies showing that S100-positive FSCs are basically not detected in the prenatal rat pituitary²².

In summary, we identified an embryonic state shift of the pituitary stem cells accompanied by negative regulation of cell proliferation and alternations in TFs and signaling pathways.

Pituitary stem cells in hybrid epithelial/mesenchymal state

Previous studies have suggested that mouse pituitary stem cells undergo an epithelial-mesenchymal transition (EMT)-like process during differentiation^23,24. We explored EMT dynamics during differentiation of the human pituitary stem cells. Among the Kyoto Encyclopedia of Genes and Genomes (KEGG) terms enriched in the stem cells, “tight junction” is related to the epithelial state, and “ECM receptor signaling pathways” is related to the mesenchymal state (Fig. 2b). By principal component analysis (PCA), all endocrine cells were ordered in a gradient transition from the stem cells to the differentiated cells on the principal component 1 (PC1) axis, so we used the PC1 axis as a trajectory for analyzing the stem cell differentiation (Fig. 3a). Examination of the expression of epithelial and mesenchymal markers showed that both the epithelial markers, such as CDH1, and the mesenchymal markers, such as VIM and CDH2, were highly expressed in the stem cells but lowly expressed in the differentiated cells; EPCAM was highly expressed in both groups (Fig. 3b and Supplementary Fig. 4b). This pattern suggested that the stem cells existed in a state expressing both the epithelial and mesenchymal markers, fitting the term of a hybrid E/M state as previously identified during mouse organogenesis and tumor transition^25,26.

To clarify the correlation between the hybrid E/M state and stemness, we defined an epithelial score (E.score), a mesenchymal score (M.score) and a stemness score (S.score) (Supplementary Data 2). We found that all three scores decreased along the differentiation axis (Fig. 3c and Supplementary Fig. 4c). This pattern was exemplified by the expression levels of known stemness markers (SOX2, SOX4, SOX9, NOTCH2 and HES1), epithelial markers (EPCAM, CDH1, KRT8, CLDN4 and GRHL2) and mesenchymal markers (CDH2, VIM, COL1A1, COL1A2 and SNAI1) which decreased along the timeline (Fig. 3c). We did not observe significant EMT changes in the stem cell substates (Supplementary Fig. 4d). The stem cells also specifically expressed CDH3, which plays roles in both maintaining stemness and promoting the hybrid E/M state in development and cancers (Fig. 2a)²⁷.

The SCENIC analysis suggested that CDH1, CDH2, CDH3 and VIM were potential targets of SOX2 and SOX4; VIM and CDH3 were potential targets of PROP1 whose targeted genes enriched in the “epithelial to mesenchymal transition” term, being consistent with the previous chromatin immunoprecipitation sequencing (ChIP-seq) study²⁴ (Supplementary Fig. 4e, f).

In sum, the data suggested that the pituitary stem cells existed in a hybrid E/M state which is associated with their stemness characters.

Reconstructing developmental trajectories

Next, we reconstructed the developmental progression of five endocrine cell lineages and identified transient precursors by applying the RNA velocity²⁸ and Slingshot²⁹ for the pseudotime trajectory analysis (Fig. 4a). The Slingshot results revealed lineage-shared and specific TFs being downregulated or upregulated during the pseudotime developmental progression, with some TFs showing peak expression in intermediated stages (Fig. 4b, Supplementary Data 3). A group of 29 TFs were identified as downregulated genes being shared among all five lineages, including PROP1, SOX2, LHX3, HES1, TCF7L1 and TGIF1 (Fig. 4c, d). Mutations in TCF7L1 and TGIF1 have been reported in patients with hypopituitarism recently^30,31. POU1F1, TBX19 and NR5A1 were strongly upregulated in the PIT-1, corticotrope and gonadotrope lineages, respectively. GATA2, FOXL2 and ISL1 were shared by the thyrotrope and gonadotrope lineages. NEUROD1 showed peak expression at the intermediated states of both the gonadotrope and corticotrope lineages (Fig. 4c).

a Pseudotime analyses of the endocrine cells shown in the UMAP plot. Left: the RNA velocity result with the arrows predicting directions of the pseudotime. Right: the Slingshot result with the lines indicating the trajectories of lineages and the arrows indicating manually added directions of the pseudotime. Dots: single cells; colors: cell types. Yellow circle in the right panel represents the start point of the trajectories which was set as the Stem1 subcluster. Since the corticotrope was already a separate cluster in the earliest sample (7 weeks), the corticotrope trajectory may start from stem cells earlier than Stem1. The CC cluster was omitted in the Slingshot pseudotime analysis as it contained several cell types and caused wrong trajectories. b Heatmap showing the relative expression of TFs displaying significant changes (P ≤ 1E−5) along the pseudotime axis of each lineage. The progenitor or precursor cells of each lineage were enclosed in frames of dashed lines. Colors: loess-smoothed expression (red, high; blue, low). The columns represent the cells being ordered along the pseudotime axis with the cell type and fetal week informations being shown above. Rows represent genes being ordered by their peak expression along the pseudotime axis. P values are determined by one way ANOVA test. c Scatterplots showing the expression levels of representative TFs along the pseudotime axis. Dots: single cells; colors: cell types. d Venn diagram of downregulated and upregulated TFs along the pseudotime axis of each lineage. Representative TFs were listed in the right boxes.

Corticotropes experience two subtates

Next, we investigated development of each cell lineage. The corticotrope is the first hormone-producing cell type that appears and TBX19 is crucial for its development³. The pseudotime and reclustering analysis identified two subclusters: Corticotrope1, which comprised cells mainly derived from the early stage (7 to 9 weeks), and Corticotrope2, which comprised cells mainly derived from the middle or late stages (10 to 25 weeks) (Fig. 5a). TBX19 and POMC were highly expressed in both clusters, which was consistent with the previous immunofluorescence staining study showing that all TBX19-positive cells also express POMC (Fig. 5b, c)³². Androgen receptor (AR) was prominently expressed in both corticotrope clusters, which was predicted to be a potential TF regulating POMC by SCENIC, and the expression was verified by immunofluorescence (Figs. 1e, 5b and c). AR mutant mice have been shown to exhibit an increased expression of POMC³³.

We identified DEGs between Corticotrope1 and Corticotrope2. DEGs in Corticotrope1 enriched GO terms such as “regulation of mitotic cell cycle” (e.g., HMGA2), suggesting that they could potentially proliferate (Fig. 5d, e, and Supplementary Data 4). Indeed, a few TBX19-positive cycling cells were captured (Fig. 2g). Interestingly, follistatin (FST), which is an antagonist of ACTVIN, was specifically expressed in Corticotrope1, suggesting that this intermediate cell type may also play some paracrine roles for regulating the gonadotropes. DEGs in Corticotrope2 enriched GO terms such as “response to corticosterone” (e.g., NR4A1 and NR4A2), suggesting that these cells were more matured (Fig. 5d, e)³⁴.

Another POMC-expressing lineage is the melanotrope in the intermediate lobe, which appears after 14 weeks of gestation but is scarce in the human pituitary¹⁵. The melanotrope coexpresses TBX19 and PAX7, both of which are required for its differentiation³⁵. Only 4 cells, which were 7- or 8-week stem cells, coexpressed TBX19 and PAX7 in our data (Fig. 1c). Comparison between these cells and other stem cells revealed only PAX7 and three other DEGs that were not expressed in the mouse and rat melanotropes. These results suggested that we have not captured the melanotrope.

In summary, we defined corticotrope differentiation from an early intermediate state to a maturing state with upregulation of the expression of genes to establish the cortisol feedback.

PIT-1 lineages segregate from a common progenitor

The PIT-1 lineage is comprised of three endocrine cell types: the somatotrope, the lactotrope and the thyrotrope, all of which are governed by POU1F1^2,36. In addition to these three hormone producing cell types, the pseudotime trajectories analysis identified three intermediated progenitor or precursor cell populations: the Pro.PIT1_all cells as a common progenitor for all three hormone producing cell types (146 cells), the Pre.Thy as a precursor for the thyrotrope (74 cells, GATA2-positive) and the Pre.Som as a potential precursor for the somatotrope (19 cells, Fig. 6a, b).

a Venn diagram of progenitor cells distributed in each PIT1 sublineage as shown in the pseudotime analysis. Pro.PIT1_all: the common progenitor cells of all three sublineages; Pre.Thy, the precursor cells of thyrotropes; Pre,Som: the potential precursor cells of somatotropes. b Distribution of the subtypes of the Pro.PIT1 on the UMAP plot. Dots: single cells; colors: weeks postfertilization (left) or cell types (right). c Violin and inside box plots showing the gene expression dynamics of known TFs for each PIT1 lineage. ^#P > 0.1, **P < 0.01 ***P < 0.001. P-values are determined by two-sided Wilcoxon test and adjusted by Bonferroni correction. d Venn diagram showing the intersections among the differentially expressed TFs between any two of the three PIT1 lineages, and the 13 TFs that were significantly upregulated in the Pro.PIT1_all cells in comparison with the stem cells, for which the names of genes are shown. e Violin with inside box plots showing representative identified TFs for each PIT1 lineage. f Heatmap showing the relative expression levels of the significantly upregulated TFs of thyrotropes along the pseudotime axis. Colors: loess-smoothed expression (red, high; blue, low). The columns represent cells being ordered along the pseudotime axis, and cell type information is shown above the heatmap. Rows represent genes being ordered by their peak expression along the pseudotime axis. g The scatterplot (left) expression and immunofluorescence staining of NKX2-2 as a somatotrope-specific TF. Scale bar, 50μm. Arrows indicate representative cells co-expressing GH and NKX2-2.

To explore how three different lineages segregate from the common progenitor cell state, we first investigated expression dynamics of known lineage-enriched genes (Fig. 6c). Previous studies have shown that mutation of Neurod4 in mice leads to almost completely lack of GHRHR expression and markedly decrease in GH expression^21,37. Mutation of Foxo1 in mice also results in delayed somatotrope differentiation³⁸. Notably, we found that NEUROD4 was clearly activated in Pro.PIT1_all cells in comparison with the stem cells (logFC = 1.6, P = 1.3E–18, Fig. 6c). NEUROD4 was further significantly upregulated in the Pre.Som cells and the differentiated somatotropes. FOXO1 was not expressed in the Pro.PIT1_all cells before upregulation in the Pre.Som cells and a portion of the somatotropes. Thus, the results demonstrated that NEUROD4 played roles earlier than FOXO1 in the somatotropes.

ZBTB20 has been recently demonstrated to be crucial for lactotrope specification in mice^39,40. This gene was expressed in all pituitary cell types including the stem cells. Notably, its expression level was significantly upregulated in the Pro.PIT1_all cells in comparison with the stem cells (logFC = 0.6, P = 9.8E–08), and the level was similar between the Pro.PIT1_all cells and the lactotropes (P > 0.1, Fig. 6c). Due to downregulation in other two lineages, the expression level of ZBTB20 was slightly but significantly higher in the lactotropes comparing with the somatotropes and the thyrotropes (Lactotrope versus Somatotrope: logFC = 0.7, P = 9E-25; Lactotrope versus Thyrotrope: logFC = 0.4, P = 1.1E-08). This expression pattern was consistent with the essential role of ZBTB20 in lactotrope specification.

GATA2, and possible ASCL1, have been implicated in thyrotrope development in mouse^18,21,41. In zebrafish, sox4b has been shown to be required for thyrotrope development by activating gata2a expression⁴². We found that ASCL1 and SOX4 were prominently expressed in the Pro.PIT1_all cells. GATA2 was not significantly upregulated in the Pro.PIT1_all cells (P > 0.1), but were prominently activated in the Pre.Thy and further upregulated in the terminal differentiated thyrotropes. Interestingly, SOX11, a member of SoxC family genes with SOX4, was significantly upregulated in the Pro.PIT1_all cells in comparison with the stem cells (logFC = 1.0, P = 4.2E-12, Fig. 6d). The SCENIC analysis also suggested that both SOX4 and SOX11 bind to the regulatory region of GATA2. ASCL1 displayed significantly higher expression levels in the thyrotropes than two other lineages (Thyrotrope versus Somatotrope: logFC = 1.3, P = 6.6E-15; Thyrotrope versus Lactotrope: logFC = 0.7, P = 3.7E-09).

Then, we investigated a comprehensive view of lineage-enriched genes and identified a total of 1,277 DEGs, including 107 TFs, between any two of three lineages (logFC ≥ 0.5 and adjusted P ≤ 0.01, Fig. 6d, Supplementary Fig. 5a, b and Supplementary Data 5). Thirteen TFs were significantly upregulated in the Pro.PIT1_all cells comparing with the stem cells. Among these TF, remarkably, NEUROD4 and ZBTB20 were identified as the primary lineage-enriched TFs for the somatotropes and the lactotropes, respectively (Fig. 6c, d). RXRG and DACH1 were identified as two thyrotrope-enriched genes, both of which were prominently upregulated in the Pro.PIT1_all cells and kept the high expression levels in the thyrotropes while downregulated in the somatotropes and the lactotropes (Fig. 6c, d).

A total of 19 somatotrope-enriched TFs, including FOXO1, CEBPD and NKX2-2, were upregulated during the differentiating process of the somatotrope. NKX2-2, which has been shown to be essential for development of neuroendocrine, gastrointestinal tract and pancreas⁴³, was specifically upregulated in a portion of terminal somatotropes (Fig. 6g), and it is highly expressed in the human adult pituitary in the GTEx database. We validated the coexpression of NKX2-2 with GH by immunofluorescence (Fig. 6g). A total of 8 lactotrope-enriched TFs (e.g. SIX6) and 38 thyrotrope-enriched (e.g., GATA2, ISL1 and FOXL2) were upregulated during the differentiation of the lactotropes and the thyrotropes, respectively (Fig. 6d, e). The pseudotime analysis of thyrotrope-lineage differentiation showed that expression of SOX4 and SOX11 peaked before activation of GATA2 in the Pre.Thy cells (Fig. 6f). SOX4 and SOX11 were significantly downregulated in the thyrotrope comparing with the Pre.Thy (SOX4: logFC = 1.0, P = 2.7E-7; SOX11: logFC = 1.3, P = 4.9E-13).

A surprising finding was that the lactotropes and the thyrotropes were close to each other relative to the somatotropes. This was displayed in UMAP, and fewer DEGs were identified between these two lineages comparing with the somatotropes, and the transcriptomes of the lactotropes and the thyrotropes were closely correlated (Figs. 4a, 6d and Supplementary Fig. 5a, b, c). In contrast, the somatotropes and lactotropes were more close to each other comparing with the thyrotropes in the mouse and rate adult pituitaries (Supplementary Fig. 5c). We speculated that this may be due to that both the lactotropes and the thyrotropes were in less matured states comparing with the somatotropes before 25 weeks. In the human fetus, GH begins to secrete before 10 weeks and peaks at approximate 22 weeks, while PRL begins to secrete at 25 weeks and peaks at birth⁴⁴. Supporting this, we found that TRHR and ESR1 were lowly expressed in the fetal thyrotropes and lactotropes, respectively.

Furthermore, despite of many evidences suggesting the existence of mammosomatotropes coexpressing GH and PRL in the human and rodent pituitary^45,46, our analysis did not identify a distinct cell cluster that corresponds to a common precursor of the somatotropes and the lactotropes, even after removing potential batch effects of the somatotropes and the lactotropes; the result was similar to that of the rat scRNA-seq study⁹ (Supplementary Fig. 5d).

Taken together, our results characterized the transcriptome dynamics during specification of the PIT-1 lineages, in which a common progenitor coexpresses lineage-enriched TFs prior to activation of divergent TF networks.

Gonadotropes exhibit two developmental trajectories

Gonadotropes mainly secrete two types of hormones, LH and FSH, which are essential for reproduction in both sexes. The pseudotime and reclustering analyses interestingly identified five subclusters (Fig. 7a). There was a clear intermediate precursor cell state, Pre.Gonado, which expressed GATA2 and FOXL2, but not NR5A1. The other four clusters (Gonadotrope1, 2, 3 and 4) comprised cells that expressed NR5A1 and GNRHR with different expression patterns of LHB, FSHB and the primate-specific hormone chorionic gonadotropin (CGBs). NR5A1, GNRHR and LHB were expressed in the four cell clusters at a similar level, while CGBs were expressed solely in Gonadotrope2 (LHB^highCGB^highFSHB^low), and FSHB was more highly expressed in Gonadotrope4 (LHB^highCGB^lowFSHB^high).

a UMAP plots (left) showing that Gonadotropes are comprised of two sublineages with scatterplots (right) showing the expression patterns of representative marker genes. Dots: single cells; colors: cell types (upper) or weeks postfertilization (lower). Gray to red indicate no to high expression levels. b RNA velocity analysis of gonadotropes projected onto the UMAP plot showing two sublineages. Dots: single cells; colors: cell types; arrows: predicted directions of the pseudotime. c Volcano plots (left) of DEGs between Gonadotrope2 and Gonadotrope4, with bar plots (right) showing GO terms enriched in each subtype. Colored dots indicate significant DEGs in each subtype. P-values are determined by two-sided Wilcoxon test and adjusted by Bonferroni correction. d Violin with inside box plots showing representative DEGs of two sublineages. e Diagram of the development of human fetal anterior pituitary. For a given stage, representative enriched TFs are highlighted. The Pre.Som is labeled with a question mark as this intermediated substate is supported by sufficient number of cells.

The pseudotime analysis revealed two developmental trajectories: one trajectory was from the Pre.Gonado to the Gonadotrope1 and terminated at the Gonadotrope2 (Type I trajectory), while the other trajectory was from the Pre.Gonado to the Gonadotrope3 and Gonadotrope4 (Type II trajectory, Fig. 7a, b). For the Type I trajectory, the intermediate Gonadotrope1 was solely in the early stages (8 to 14 weeks) while Gonadotrope2 was comprised of both early and late stages; all Type II trajectory cells were in the late stages (15 to 25 weeks). These results indicated that the Type I and II trajectories represented an early and a late gonadotrope lineage, respectively. Among all five clusters, Gonadotrope2 had the most DEGs, and comparing between Gonadotrope2 and Gonadotrope4 identified 265 and 30 DEGs, respectively (Fig. 7c, Supplementary Fig. 6 and Supplementary Data 6). GO analysis showed that Gonadotrope2 DEGs enriched “regulated exocytosis” and “C21 − steroid hormone biosynthetic process” (e.g. the steroidogenic enzyme CYP11A1), suggesting that Gonadotrope2 actively secreted hormones (Fig. 7c). WNT4 and GATA2 were more highly expressed in Gonadotrope2 and Gonadotrope4, respectively (Fig. 7d). Other DEGs included HIF3A and MC2R in the Type I trajectory and folate receptor FOLR1 and secretoglobin SCGB2A1 in the Type II trajectory (Fig. 7d and Supplementary Fig. 6).

Together, these data determined two gonadotrope developmental trajectories including a previously unappreciated LHB^highCGB^highFSHB^low subcluster.

Species comparison between human and rodent pituitaries

Next, we compared our scRNA-seq data with two recently published scRNA-seq datasets of mouse and rat adult pituitaries^9,10. The human scRNA-seq data, which used a plate-based method, recovered higher number of genes per cell comparing with the rodent data using the droplet-based 10X genomic method (Supplementary Fig. 7a). Integrating three datasets showed that most pituitary cell types were conserved among human and rodent (Fig. 8a). Anterior pituitary known markers and new markers, including SOX2, POMC, GH1, PRL, TSHB, GNRHR, ALDH1A2, NR3C1, DLK1, OLFM1, DIO2 and KCNK3, showed similar cell-type-specific expression patterns (Fig. 8b, Supplementary Fig. 7b).

Notably, no progenitor or precursor cell types were found in the rodent datasets, indicating that differentiation rarely occurs in these adult pituitaries. Also, in the human fetal pituitaries, the proliferating cells were mainly the stem cells, while in the adult rodent pituitaries, the proliferating cells were mainly the somatotropes and the lactotropes. The rat FSCs were clustered closely to human and mouse stem cells; all expressed SOX2, but only the FSCs expressed S100B.

We then attempted to find species-specific genes, despite that any differences between the fetal human dataset and the adult rodent datasets could reflect the species or stage differences. We recognized cell-type-specific genes for each species and then identified human (or fetal)-specific genes and rodent (or adult)-specific genes (Supplementary Data 7). Among the identified genes, CGB and GH2 are primate-specific and thus do not exist in the rodent data (Fig. 8c). Trhr and Esr1 were highly expressed in the rodent thyrotropes and lactotropes, respectively, but lowly expressed in the human fetal pituitary, thus reflecting stage differences. Interestingly, AR was highly expressed in the human fetal corticotropes while lowly expressed in the rodent adult corticotropes instead with high expression in the gonatropes (Fig. 8c). Also, PRLHR was highly expressed in the human fetal gonatropes while nearly not expressed in rodent⁴⁷. IL17RB and GAD2 also showed specific expression in the human fetal pituitary. Together, these results revealed a small number of potential species-specific genes.