Single-cell transcriptome profiling reveals several LncRNAs differentially expressed in idiopathic germ cell aplasia

Abstract

Mechanisms underlying severe male infertility are still largely elusive. However, recently, a single-cell transcription study by our group identified several differentially expressed coding genes in all the somatic cell types in testes of patients with idiopathic germ cell aplasia (iGCA). Here, we leverage this work by extending the analysis also to the non-coding portion of the genome. As a result, we found that 43 LncRNAs were differentially expressed in the somatic cells of these patients. Interestingly, a significant portion of the overexpressed LncRNAs was found to be a target of TAF9B, a transcription factor known to be involved in germ cell survival. Moreover, several overexpressed LncRNAs were also found to be activated in a mouse model of Sertoli cells treated with bisphenol A, a widespread environmental contaminant, long suspected to impair male fertility. Finally, a literature search for MEG3, a maternally imprinted LncRNA overexpressed as well in our patients, found it to be involved, among other things, in obesity and inflammation, known comorbidities of iGCA, ultimately suggesting that our findings deepen the understanding of the molecular insights coupled not only to the pathogenesis, but also to the clinical course of this class of patients.

Introduction

The successful completion of spermatogenesis relies on the availability of spermatogonial stem cells, along with their capability to proliferate and to transform, first, into progenitor spermatogonia and, then, into spermatozoa (Oatley and Brinster, 2008). Sertoli-cell-only (SCO) syndrome, also known as germ cell aplasia (GCA), represents a condition of the testes where only Sertoli cells line the seminiferous tubules, ultimately leading to the most severe form of male infertility, the non-obstructive azoospermia (NOA) (Ramphul and Mejias, 2021).

While a genetic component has been determined for NOA (Nakamura et al., 2017), it is becoming clear that, in the vast majority of cases, there is not an evident underlying cause, leading to the broad definition of those cases as the wide group idiopathic NOA (iNOA) (European Association of Urology, 2022). Instead, there is overwhelming evidence that the male infertility status, and mainly in NOA cases, is linked with an augmented risk of diseases associated with aging, like type II diabetes, cardiovascular disease, autoimmune disease, obesity and cancers (Eisenberg et al., 2015). All this has brought to the development of a new conceptualization of the male fertility status as a proxy of the overall men’s health (Salonia et al., 2009).

Recently, single-cell RNA-sequencing (scRNA-seq) analysis from our group identified eight cell clusters in the testicular somatic cells populations of iGCA patients (Alfano et al., 2021). Thanks to the use of cell type marker genes on these clusters, the main somatic cell populations were recognized, Leydig (LEY), myoid (MYD), Sertoli (SRT) and endothelial (END) cells. Moreover, immune cells were also recognized, like macrophages (MCR) and T-cells (TCL). Finally, the stromal (STRO) cluster was assigned to pericytes or vascular smooth muscle cells, whereas one cluster, lacking clear marker genes, remained undetermined (UND). Identification of DE (DE) coding genes, lead to the identification of molecular pathways related to aging, inflammation and DNA damage, offering molecular insights into the pathogenesis of iGCA.

In the present study, we dissected the dysregulation of the non-coding transcripts in iGCA by focussing our attention to the analysis of the LncRNAs of our dataset. As a result, we found 43 DE LncRNAs in iGCA patients. Our analysis revealed a close connection between differential LncRNAs and several features of this disease, such as the survival of germ cells and the presence of comorbidities, like obesity and inflammation. Also, an environmental contribution to the disease insurgence was found to be highly consistent with our data.

Methods

Identification of DE LncRNAs

Tissue processing, ethical approval and the scRNA-seq procedure were already described (Alfano et al., 2021). DE LncRNAs in iGCA patients of this study were identified by mining DE transcript names in Supplementary Dataset 29 for the presence of the term “lncRNA” in the “gene_type” field of the “gencode.v39. annotation.gtf” file from the GENCODE project (www.gencodegenes.org). Normalized transcript expression data were extracted and used to build the heat map and the dot plot shown, respectively, in Figure 1A and Figure 1B. Normalized expression counts, along with pertinent sample metadata, are available as Supplementary Data respectively in the files “Normalized_LncRNA_counts.csv” and “Sample-phenotypes.csv”. Upstream and downstream transcripts of the DE LncRNAs located within 500 Kb of the beginning of each locus were obtained from the ‘gencode.v39. annotation.gtf’ file.

FIGURE 1

Literature search

In order to evaluate the consistency of the MEG3 role in (IR) injury, we performed a literature search on 11 February 2022 using the keywords “MEG3 AND Reperfusion Injury”, “MEG3 AND infarction” and “MEG3 AND organ failure”. Both titles and abstracts were retrieved and, whenever it deemed to be necessary, also the full article was analyzed. A similar search strategy was used for evaluating the correlation between MEG3 expression and obesity/diabetes with Pubmed abstracts being retrieved using the keywords “MEG3 AND obesity” and “MEG3 AND diabetes”. Statistical significance was calculated by using a binomial test.

Identification of signatures associated to DE LncRNAs

All the datasets (C1—> C8 plus the hallmark gene set) from the MSigDB 7.5.1 database (Liberzon et al., 2015) were searched for overlaps with our DE LncRNAs using the “Investigate Gene Set” option with a FDR of 0.001. A custom R script was used to analyze a non-redundant, 147 libraries dataset available at the Enrichr resource (Xie et al., 2021) through the CRAN library “enrichr”. Differential expression of the LncRNAs was measured using the “FindMarkers” program from the R Seurat software package (Satija et al., 2015), using default parameters (Alfano et al., 2021). Differerential expression for the TAF9B transcript in the LEY and MYD cell types was measured using the ‘FindMarkers’ program as well, except that the “logfc.threshold” parameter, controlling the minimum differential expression fold change cutoff to be used between the two groups of cells, was set to zero instead of the default value of 0.25.

Results

Identification of DE LncRNAs in idiopathic germ cell aplasia patients

We searched our previously identified (Alfano et al., 2021) dataset of DE transcripts in iGCA patients for the presence of LncRNAs. This dataset profiled, by means of scRNA-seq, gene expression in each somatic cell type in testes of men with iGCA (Alfano et al., 2021). As a result, 43 LncRNAs were found to be DE in one or more cell types, as shown in the heatmap in Figure 1A and in Table 1. Further details on their genomic localization and on their differential expression are reported, respectively, in Supplementary Table S1 and in the dot plot comparison in Figure 1B, with the dot color intensity proportional to the average expression level for each group and the dot size related to the percentage of expressing samples. Among these transcripts, differential expression was very consistent, with 21 LncRNAs being downregulated in one or more cell types and 21 LncRNAs being instead upregulated in one or more cell types. Only the LncRNA NEAT1 presented a mixed expression, being downregulated in endothelial cells and upregulated in Sertoli cells (Figure 1B and Table 1).

Enrichment of DE LncRNAs in biological signatures

To evaluate the enrichment of the DE LncRNAs in specific lists or biological signatures, all the datasets from the Molecular Signature Database (Liberzon et al., 2015) were firstly searched. Since some of our cell types were not amenable to meaningful searches because they contained only a few DE LncRNAs, like for example STRO cells possessing only three downregulated LncRNAs (see Table 1), a pool of all cell types was used in our searches. To compensate for this added heterogeneity, a rather highly stringent FDR threshold, 0.001, was used. As a result, while no significant enrichment was found in the downregulated LncRNAs, the overexpressed ones were found to be significantly enriched in chromatin immunoprecipitation targets of the TATA-Box Binding Protein Associated Factor 9b (TAF9B), a basal component of the nuclear transcription machinery (Frontini et al., 2005), as shown in Table 2. Previous analysis (Alfano et al., 2021) was unable to detect any upregulation of TAF9B in cell types overexpressing a significant fraction of its target LncRNAs: LEY (6 target LncRNAs) and MYD (7 target LncRNAs). However, this analysis was performed using a relatively large log fold change (FC) cutoff, 0.25. Since for regulatory molecules, like TAF9B, even smaller FC can have a great significance in high throughput, transcriptional profiling experiments (Dalman et al., 2012), we repeated the analysis without employing this cutoff. As a result, it was found that TAF9B was significantly overexpressed, albeit with a lower log FC, both in LEY and MYD cell types, as shown in Supplementary Table S2, ultimately supporting its potential role as a transcriptional activator of DE LncRNAs in iGCA.We also looked for enrichments of our DE LncRNAs in the several datasets available at the Enrichr resource (Xie et al., 2021). Also in this case, no hit was found for the downregulated ones, but, as shown in Table 2, a highly significant enrichment was found with the genes overexpressed in a mouse model of Sertoli cells treated with Bisphenol A (BPA) (Tabuchi et al., 2006), a compound employed in the preparation of various plastics and known to be an endocrine disruptor, exhibiting weak estrogenic, anti-thyroid and anti-androgenic activities (Peretz et al., 2014). Similar results were obtained using single, i.e., no pooled, cell types harboring a high number of DE LncRNAs, like LEY (n = 16) and MYD (n = 18) cells (Supplementary Table S3).

Literature survey and meta-analysis of the LncRNA MEG3

A bibliographical search, aimed at evaluating possible relationships of the DE LncRNAs to disease features, was also undertaken. We started our analysis with the MEG3 LncRNA because it had a relatively high number of abstracts available in Pubmed (#931, as of 11 February 2022). MEG3, overexpressed in LEY, MYD and STRO cell types of our patients, was an interesting molecule also because of its role as a maternally imprinted transcript (Schmidt et al., 2000). Moreover, given the relevant role played by inflammation in the testes of our patients (Alfano et al., 2021), it was i intriguing that, recently, MEG3 was found to promote pyroptosis, i.e., a form of cell-death associated with inflammatory signals, in testicular ischemia-reperfusion (IR) injury (Ning et al., 2021). In order to evaluate the consistency of the MEG3 role in IR injury, we performed a literature search using pertinent keywords (see Methods). A total of 29 experimental, non-redundant articles were retrieved and the direction of MEG expression (up or down) was evaluated. Results, shown in Table 3, indicated that MEG3 expression was upregulated in all cases, strongly validating it as a potential therapeutic target in our patients. Given the relevant role of the body-mass index (BMI), obesity and diabetes as comorbidities in iGCA (Boeri et al., 2022), we also performed a literature search in order to determine MEG3 direction of expression in these circumstances. Also in this case, albeit not exclusively, MEG3 upregulation was found to be correlated to these diseases (p = 0.02), as shown in Supplementary Table S4.

Discussion

Here, we present a follow-up study to our previous characterization of the transcriptional atlas of the testicular somatic cell populations linked to human iGCA (Alfano et al., 2021). Previous work had indeed determined the contribution of several coding genes to the pathological status of these patients, underlying their relevance for frequently observed features in iGCA, like senescence of the testicular somatic cells, immaturity of LEY cells, persistent inflammation, DNA damage and defects in gene imprinting. The current work extends the previous analysis by looking also at the non-coding portion of the genome. While the inevitably low depth of the scRNA seq approaches somehow limits the identification of weakly expressed molecules like ncRNAs (Pokhilko et al., 2021), we still were able to identify 43 LncRNAs DE in iGCA patients. Also, in light of the extensive spectrum of comorbidities observed in unfertile males (Salonia et al., 2009; Shiraishi and Matsuyama, 2018; Boeri et al., 2022), the study of regulatory transcripts, like LncRNAs, in this class of patients appears to be of paramount importance.

Overexpressed LncRNAs in our patients were found to overlap a small, but biologically very relevant number of lists. For instance, the overlap with the target genes of TAF9B, a component of TFIID upregulated during the germline development (Gura et al., 2020), is rather intriguing. Indeed, it is in keeping with recent findings showing that cells of the germline use a highly specialized form of TFIID, the transcription initiation complex of nuclear genes that once was thought to be, instead, rather general and ubiquitous in all tissues (Hoey et al., 1990). Interestingly, TAF9B is also within the network of proteins regulated by the gene Deleted in Azoospermia Like (DAZL), whose function is essential for spermatogenesis and meiosis completion (Zagore et al., 2018; Li et al., 2019). Summing up, while the formal proof of the transcriptional activation of the identified LncRNAs by TAF9B must await further studies, present data are already suggestive of an involvement of this transcription factor in iGCA.

Of interest also appears the overlap with the transcripts stimulated by BPA exposure in mouse Sertoli cells, since the impact of the environment on reproductive health has been reported as an explanation for idiopathic infertility (Shi et al., 2017). It will be interesting to see if BPA can induce also the expression of the identified LncRNAs in human testicular somatic cells in ex-vivo experiments. BPA, along with its analogs, is present in several daily use products, like manufactured plastics, cans and paper (Wong and Durrani, 2017). It is known to possess a variety of deleterious effects on human health, including male reproduction and fertility (Siracusa et al., 2018). A replacement of BPA with more chemically stable analogs has been attempted in the manufacturing process. However, these analogs have proved to be both worse in terms of biodegradability and better in terms of dermal penetration (Liao et al., 2012). BPA, along with its analogs, contaminate the environment, including air, water, food and house dust (Wu et al., 2018). Therefore, the primary intake of bisphenols in humans is likely to occur through the diet, mainly through the consumption of canned foods or drinking bottled water (Eladak et al., 2015). All this ultimately suggests a limitation of their use, also in consideration of the wide spectrum of comorbidities being associated to BPA usage (Ma et al., 2019; Barra et al., 2022), even through the maternal route (Blaauwendraad et al., 2022).

A literature survey also found MEG3 consistently linked to IR injury and obesity/diabetes, validating it, or molecules in its pathway, as a therapeutic target in IGCA. MEG3 contribution to the inflammatory phenomena of this disease could be, indeed, rather important. Recently, inflammation induced by lipopolysaccharides (LPS) in Leydig cells was found to decrease both testosterone levels and cell viability, along with a MEG3 transcriptional activation (Zhou et al., 2021). Ablation of MEG3 expression in this setting attenuated the inflammatory response and increased testosterone levels through the sponging effect on miR-93-5p (Zhou et al., 2021). The sponging capabilities of MEG3, mainly related to the RI phenomenon, have been shown also in several other works (reviewed in Zhao et al., 2022). Therefore, miRNAs in its pathway are likely to be therapeutic targets as well in inflammatory phenomena and should be considered for therapeutic intervention.

References

• 1

  AlfanoM.TasciniA. S.PederzoliF.LocatelliI.NebuloniM.GianneseF.et al (2021). Aging, inflammation and DNA damage in the somatic testicular niche with idiopathic germ cell aplasia. Nat. Commun.12, 5205–5221. 10.1038/s41467-021-25544-0

  • CrossRef
  • Google Scholar

• 2

  BarraN. G.KwonY. H.MorrisonK. M.SteinbergG. R.WadeM. G.KhanW. I.et al (2022). Increased gut serotonin production in response to bisphenol A structural analogs may contribute to their obesogenic effects. Am. J. Physiol. Endocrinol. Metab.323, E80–E091. 10.1152/ajpendo.00049.2022

  • CrossRef
  • Google Scholar

• 3

  BlaauwendraadS. M.GaillardR.SantosS.SolC. M.KannanK.TrasandeL.et al (2022). Maternal phthalate and bisphenol urine concentrations during pregnancy and early markers of arterial health in children. Environ. Health Perspect.130, 47007. 10.1289/EHP10293

  • CrossRef
  • Google Scholar

• 4

  BoeriL.VentimigliaE.CazzanigaW.PederzoliF.FallaraG.PozziE.et al (2022). Risk of health status worsening in primary infertile men: A prospective 10-year follow-up study. Andrology10, 128–136. 10.1111/andr.13090

  • CrossRef
  • Google Scholar

• 5

  DalmanM. R.DeeterA.NimishakaviG.DuanZ. H. (2012). Fold change and p-value cutoffs significantly alter microarray interpretations.BMC Bioinforma.13, S11. 10.1186/1471-2105-13-S2-S11

  • CrossRef
  • Google Scholar

• 6

  EisenbergM. L.LiS.BehrB.PeraR. R.CullenM. R. (2015). Relationship between semen production and medical comorbidity. Fertil. Steril.103, 66–71. 10.1016/j.fertnstert.2014.10.017

  • CrossRef
  • Google Scholar

• 7

  EladakS.GrisinT.MoisonD.GuerquinM. J.N'Tumba-BynT.Pozzi-GaudinS.et al (2015). A new chapter in the bisphenol A story: bisphenol S and bisphenol F are not safe alternatives to this compound. Fertil. Steril.103, 11–21. 10.1016/j.fertnstert.2014.11.005

  • CrossRef
  • Google Scholar

• 8

  European Association of Urology (2022). http://uroweb.org/guideline/male-infertility/ (Accessed February 8, 2022).

  • Google Scholar

• 9

  FrontiniM.SoutoglouE.ArgentiniM.Bole-FeysotC.JostB.ScheerE.et al (2005). TAF9b (formerly TAF9L) is a bona fide TAF that has unique and overlapping roles with TAF9. Mol. Cell. Biol.25, 4638–4649. 10.1128/MCB.25.11.4638-4649.2005

  • CrossRef
  • Google Scholar

• 10

  GuraM. A.MikedisM. M.SeymourK. A.de RooijD. G.PageD. C.FreimanR. N. (2020). Dynamic and regulated TAF gene expression during mouse embryonic germ cell development. PLoS Genet.16, e1008515. 10.1371/journal.pgen.1008515

  • CrossRef
  • Google Scholar

• 11

  HoeyT.DynlachtB. D.PetersonM. G.PughB. F.TjianR. (1990). Isolation and characterization of the Drosophila gene encoding the TATA box binding protein, TFIID. Cell.61, 1179–1186. 10.1016/0092-8674(90)90682-5

  • CrossRef
  • Google Scholar

• 12

  LiH.LiangZ.YangJ.WangD.WangH.ZhuM.et al (2019). DAZL is a master translational regulator of murine spermatogenesis. Natl. Sci. Rev.6, 455–468. 10.1093/nsr/nwy163

  • CrossRef
  • Google Scholar

• 13

  LiaoC.LiuF.MoonH. B.YamashitaN.YunS.KannanK. (2012). Bisphenol analogues in sediments from industrialized areas in the United States, Japan, and korea: Spatial and temporal distributions. Environ. Sci. Technol.46, 11558–11565. 10.1021/es303191g

  • CrossRef
  • Google Scholar

• 14

  LiberzonA.BirgerC.ThorvaldsdóttirH.GhandiM.MesirovJ. P.TamayoP. (2015). The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell. Syst.1, 417–425. 10.1016/j.cels.2015.12.004

  • CrossRef
  • Google Scholar

• 15

  MaY.LiuH.WuJ.YuanL.WangY.DuX.et al (2019). The adverse health effects of bisphenol A and related toxicity mechanisms. Environ. Res.176, 108575. 10.1016/j.envres.2019.108575

  • CrossRef
  • Google Scholar

• 16

  NakamuraS.MiyadoM.SaitoK.KatsumiM.NakamuraA.KoboriY.et al (2017). Next-generation sequencing for patients with non-obstructive azoospermia: Implications for significant roles of monogenic/oligogenic mutations. Andrology5, 824–831. 10.1111/andr.12378

  • CrossRef
  • Google Scholar

• 17

  NingJ. Z.HeK. X.ChengF.LiW.YuW. M.LiH. Y.et al (2021). Long non-coding RNA MEG3 promotes pyroptosis in testicular ischemia-reperfusion injury by targeting MiR-29a to modulate PTEN expression. Front. Cell. Dev. Biol.9, 671613. 10.3389/fcell.2021.671613

  • CrossRef
  • Google Scholar

• 18

  OatleyJ. M.BrinsterR. L. (2008). Regulation of spermatogonial stem cell self-renewal in mammals. Annu. Rev. Cell. Dev. Biol.24, 263–286. 10.1146/annurev.cellbio.24.110707.175355

  • CrossRef
  • Google Scholar

• 19

  PeretzJ.VroomanL.RickeW. A.HuntP. A.EhrlichS.HauserR.et al (2014). Bisphenol a and reproductive health: Update of experimental and human evidence, 2007-2013. Environ. Health Perspect.122, 775–786. 10.1289/ehp.1307728

  • CrossRef
  • Google Scholar

• 20

  PokhilkoA.HandelA. E.CurionF.VolpatoV.WhiteleyE. S.BøstrandS.et al (2021). Targeted single-cell RNA sequencing of transcription factors enhances the identification of cell types and trajectories. Genome Res.31, 1069–1081. 10.1101/gr.273961.120

  • CrossRef
  • Google Scholar

• 21

  RamphulK.MejiasS. G. (2021). Sertoli-cell-only syndrome. Treasure Island (FL): In: StatPearls.

  • Google Scholar

• 22

  SaloniaA.MatloobR.GallinaA.AbdollahF.SaccàA.BrigantiA.et al (2009). Are infertile men less healthy than fertile men? Results of a prospective case-control survey. Eur. Urol.56, 1025–1031. 10.1016/j.eururo.2009.03.001

  • CrossRef
  • Google Scholar

• 23

  SatijaR.FarrellJ. A.GennertD.SchierA. F.RegevA. (2015). Spatial reconstruction of single-cell gene expression data. Nat. Biotechnol.33, 495–502. 10.1038/nbt.3192

  • CrossRef
  • Google Scholar

• 24

  SchmidtJ. V.MattesonP. G.JonesB. K.GuanX. J.TilghmanS. M. (2000). The Dlk1 and Gtl2 genes are linked and reciprocally imprinted. Genes. Dev.14, 1997–2002. 10.1101/gad.14.16.1997

  • CrossRef
  • Google Scholar

• 25

  ShiM.SekulovskiN.MacLeanJ. A.2ndHayashiK. (2017). Effects of bisphenol A analogues on reproductive functions in mice. Reprod. Toxicol.73, 280–291. 10.1016/j.reprotox.2017.06.134

  • CrossRef
  • Google Scholar

• 26

  ShiraishiK.MatsuyamaH. (2018). Effects of medical comorbidity on male infertility and comorbidity treatment on spermatogenesis. Fertil. Steril.110, 1006–1011. 10.1016/j.fertnstert.2018.07.002

  • CrossRef
  • Google Scholar

• 27

  SiracusaJ. S.YinL.MeaselE.LiangS.YuX. (2018). Effects of bisphenol A and its analogs on reproductive health: A mini review. Reprod. Toxicol.79, 96–123. 10.1016/j.reprotox.2018.06.005

  • CrossRef
  • Google Scholar

• 28

  TabuchiY.TakasakiI.KondoT. (2006). Identification of genetic networks involved in the cell injury accompanying endoplasmic reticulum stress induced by bisphenol A in testicular Sertoli cells. Biochem. Biophys. Res. Commun.345, 1044–1050. 10.1016/j.bbrc.2006.04.177

  • CrossRef
  • Google Scholar

• 29

  WongK. H.DurraniT. S. (2017). Exposures to endocrine disrupting chemicals in consumer products-A guide for pediatricians. Curr. Probl. Pediatr. Adolesc. Health Care47, 107–118. 10.1016/j.cppeds.2017.04.002

  • CrossRef
  • Google Scholar

• 30

  WuL. H.ZhangX. M.WangF.GaoC. J.ChenD.PalumboJ. R.et al (2018). Occurrence of bisphenol S in the environment and implications for human exposure: A short review. Sci. Total Environ.615, 87–98. 10.1016/j.scitotenv.2017.09.194

  • CrossRef
  • Google Scholar

• 31

  XieZ.BaileyA.KuleshovM. V.ClarkeD.EvangelistaJ. E.JenkinsS. L.et al (2021). Gene set knowledge discovery with Enrichr. Hoboken, New Jersey: Current protocols, e90.

  • Google Scholar

• 32

  ZagoreL. L.SweetT. J.HanniganM. M.Weyn-VanhentenryckS. M.JobavaR.HatzoglouM.et al (2018). DAZL regulates germ cell survival through a network of PolyA-proximal mRNA interactions. Cell. Rep.25, 1225–1240. e6. 10.1016/j.celrep.2018.10.012

  • CrossRef
  • Google Scholar

• 33

  ZhaoY.LiuY.ZhangQ.LiuH.XuJ. (2022). The mechanism underlying the regulation of long non-coding RNA MEG3 in cerebral ischemic stroke. Cell. Mol. Neurobiol.10.1007/s10571-021-01176-2

  • CrossRef
  • Google Scholar

• 34

  ZhouX.HeJ.ChenJ.CuiY.OuZ.ZuX.et al (2021). Silencing of MEG3 attenuated the role of lipopolysaccharides by modulating the miR-93-5p/PTEN pathway in Leydig cells. Reprod. Biol. Endocrinol.19, 33. 10.1186/s12958-021-00712-5

  • CrossRef
  • Google Scholar