Transcriptome Analysis of Testicular Aging in Mice

Male reproductive aging, or andropause, is associated with gradual age-related changes in testicular properties, sperm production, and erectile function. The testis, which is the primary male reproductive organ, produces sperm and androgens. To understand the transcriptional changes underlying male reproductive aging, we performed transcriptome analysis of aging testes in mice. A total of 31,386 mRNAs and 9387 long non-coding RNAs (lncRNAs) were identified in the mouse testes of diverse age groups (3, 6, 12, and 18 months old) by total RNA sequencing. Of them, 1571 mRNAs and 715 lncRNAs exhibited changes in their levels during testicular aging. Most of these aging-related transcripts exhibited slight and continuous expression changes during aging, whereas some (9.6%) showed larger expression changes. The aging-related transcripts could be classified into diverse expression patterns, in which the transcripts changed mainly at 3–6 months or at 12–18 months. Our subsequent in silico analysis provided insight into the potential features of testicular aging-related mRNAs and lncRNAs. We identified testis-specific aging-related transcripts (121 mRNAs and 25 lncRNAs) by comparison with a known testis-specific transcript profile, and then predicted the potential reproduction-related functions of the mRNAs. By selecting transcripts that are altered only between 3 and 18 months, we identified 46 mRNAs and 34 lncRNAs that are stringently related to the terminal stage of male reproductive aging. Some of these mRNAs were related to hormonal regulation. Finally, our in silico analysis of the 34 aging-related lncRNAs revealed that they co-localized with 19 testis-expressed protein-coding genes, 13 of which are considered to show testis-specific or -predominant expression. These nearby genes could be potential targets of cis-regulation by the aging-related lncRNAs. Collectively, our results identify a number of testicular aging-related mRNAs and lncRNAs in mice and provide a basis for the future investigation of these transcripts in the context of aging-associated testicular dysfunction.


Introduction
Male reproductive aging, or andropause, refers to age-related changes in male reproductive integrity, including sperm production and erectile function [1]. Unlike menopause resulting from female reproductive aging, which is relatively abrupt, male reproductive aging is accompanied by slow changes in function [1,2]. These gradual dysfunctional changes occur mostly in the testis, which is the main male reproductive organ and is responsible for producing sperm and hormones [1][2][3]. Slow and progressive decreases in the testicular mass and testosterone level, erectile dysfunction, and tubular sclerosis are observed during human andropause [2,4].
Spermatogenesis is the unique, complex, and tightly regulated process of male germ cell development [5]. It includes successive mitotic division, meiosis, and post-meiotic phases through which spermatogonial stem cells (SSCs) are differentiated into spermatocytes and spermatids respectively. A highly organized intrinsic genetic network is important for spermatogenesis [5]. To elucidate the molecular nature of the intrinsic program of spermatogenesis, various investigations have sought to identify and characterize relevant protein-coding genes [6] and noncoding RNAs, such as microRNAs (miRNAs) [7], piwi-interacting RNAs (piRNAs) [8], and long non-coding RNAs (lncRNAs) [9][10][11][12][13][14][15]. LncR-NAs could function as novel regulatory molecules which interact with DNA, RNA, and proteins. They are implicated in various biological process, including differentiation and development. Our group previously found that the testis harbors many tissue-specific lncRNAs [16]. We recently reported that the testicular germ cell-specific lncRNA, Teshl, functions to regulate the expression of sex chromosome genes and maintain the quality of Y chromosome-bearing sperm [17].
Testicular aging has been studied in humans [1][2][3][4][17][18][19], as well as in murine [20] and canine species [21]. In humans and canines, decreases in Leydig cell populations are commonly observed [19,21]. In the case of mouse testicular aging, plasma testosterone level was significantly decreased in aged groups (over age 450 days). It is unclear whether there are changes in semen parameters, such as sperm concentration, motility, and morphology, during aging [22]. An integrative study with genotype tissue expression (GTEx) data identified 22 and 7 aging-associated coding (mRNA) and non-coding (lncRNA) genes, respectively, in human testes [18]. However, although numerous aging-related studies have been performed on mouse tissues including testes and various tissues of other species, no previous study has undertaken aging-related transcriptomic profiling of whole mouse testes in a comprehensive manner.
Here, we identified and profiled mRNAs and lncRNAs associated with mouse testicular aging through total RNA sequencing analysis. To comprehensively investigate the expression changes of transcripts during all periods of aging, we analyzed the testicular transcripts obtained at postnatal months 3, 6, 12, and 18. We newly identified aging-related transcripts, observed that they typically exhibited slight and gradual expression changes, classified them into a number of different expression patterns, and assessed the potential features of the aging-related mRNAs and lncRNAs. To the best of our knowledge, this is the first study seeking to identify and characterize mRNAs and lncRNAs related to testicular aging in mice. Our study provides inclusive transcriptomic information that should facilitate future studies on the testicular dysfunction that occurs with age in mammalian species.

Animals and RNA Preparation
We used C57BL/6 male mice representing four age groups: postnatal 3, 6, 12, and 18 months old. Three life phases (3-6 months of age for 'mature adult', 10-15 months of age for 'middle-aged', and 18-24 months of age for 'old') have been suggested for aging studies in mice (https://www.jax.org/research-and-faculty/research-labs/the-harrisonlab/gerontology/life-span-as-a-biomarker, accessed on 3 September 2019). All mice were provided by the Korea Basic Science Institute (KBSI, Gwangju, Korea). It is unclear whether mice of each age group were from the same litter. Each experiment was performed in triplicate. Testes were obtained from mice and immediately cryopreserved with liquid nitrogen and stored at −80 • C. Total RNA was isolated. Each experiment involving animals was conducted in accordance with the Korean Food and Drug Administration (KFDA) test guidelines. The protocols were reviewed and approved by the Institutional Animal Care and Use Committees (IACUC) of Gwangju Institute of Science and Technology (GIST) (permit number: GIST-2019-097).

Total RNA Sequencing
Sequencing libraries were prepared using a TruSeq RNA Sample Prep kit v2 (Illumina, San Diego, CA, USA) according to the provided protocol. Library preparation involved removing ribosomal RNA with a Ribo-Zero rRNA Removal Kit (Illumina), and then performing random fragmentation and cDNA synthesis. Each sample library was subjected to paired-end high-throughput RNA sequencing using a HiSeq 4000 system (Illumina). The raw reads obtained from high-throughput sequencing were subjected to quality control and preprocessing. Alignment was performed with HISAT2 against the UCSC mm10 reference genome. Expression profiles were obtained with the FPKM (fragments per kilobase of exon per million fragments) method. To determine the expression changes of transcripts, we used fold-change with pseudocount value (0.001). The sequencing data were uploaded to the Gene Expression Omnibus (GEO) database at the National Center for Biotechnology Information (NCBI) under GEO accession number GSE175633.

Tissue Expression Estimation
We previously performed microarray analysis of five different mouse tissues (brain, heart, kidney, liver, and testis) to identify tissue-specific expression of mRNAs and lncR-NAs [16]. Here, we identified age-related mRNAs and lncRNAs showing apparent testisspecific expression by comparing the previous tissue expression data with the newly identified age-related transcripts. We converted the transcript IDs of the age-related and testisspecific transcripts into ENSEMBL gene IDs using the DAVID gene ID conversion tool.

Investigating the Potential Features of Aging-Related Transcripts
We selected aging-related mRNAs and lncRNAs to investigate their potential features. A protein-protein network was constructed with the 46 mRNAs using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) [23]. To investigate cis-regulatory targets of the 34 aging-related lncRNAs, we used Genomic Regions Enrichment of Annotations Tool (GREAT) to identify their neighboring genes (Basal plus extension; proximal, 5 kb upstream and 1 kb downstream; distal, up to 1000 kb) [24].

Identification of mRNAs and lncRNAs in Mouse Testes during Aging
To comprehensively identify transcripts expressed in mouse testes at different ages, we prepared testes from C57BL/6 mice representing four age groups: postnatal 3 months old (3 M), 6 months old (6 M), 12 months old (12 M), and 18 months old (18 M). We performed total RNA sequencing analysis of testicular tissues using an Illumina HiSeq 4000. The raw data comprised 65 to 90 million total reads, over 95% of which were clean reads (Table 1). These data were deposited to the Gene Expression Omnibus (GEO) under GEO accession number GSE175633. The testicular transcripts obtained from the RNA sequencing data were processed into individual transcripts using StringTie [25]. A total of 71,967 transcripts were assembled from the entire set of raw RNA sequences. From the total set of individual transcripts, 31,386 transcripts were identified as mR-NAs based on our analysis using the MGI (Mouse Genome Informatics) database and GffCompare [26]. Total RNA sequencing can also provide information on the expression profiles of non-coding RNAs. To identify lncRNAs from the transcript assembly, we designed an in silico pipeline ( Figure 1A). First, based on the definition of an lncRNA, we selected transcripts longer than 200 nucleotides (nt). We then filtered out single-exon transcripts, which yielded 64,957 multi-exon transcripts, to eliminate experimental artifacts and background noise. As a consequence, single-exon lncRNAs were excluded from our lists. Finally, we assessed the coding potential of these transcripts using CPC (Coding Potential Calculator) [27], CPAT (Coding Potential Assessment Tool) [28], txCDSPredict (provided by kentUtils), and HMMSearch (provided by HMMER) against the pfam database [29] ( Figure 1B). Of the transcripts, 9387 were commonly evaluated as non-coding sequences by these tools and were thus considered to be testicular lncRNAs. Further classification of these lncRNAs revealed that 2152 were known non-coding transcripts, 1734 were novel isoforms of known transcripts, and the remaining 5274 were novel unannotated transcripts (NUTs) ( Figure 1C).

Global Expression and Transcriptomic Features of mRNAs and lncRNAs Expressed in Mouse Testes during Aging
We characterized the expression and transcriptomic features of the mRNA and lncRNA transcripts identified by our total RNA sequencing. The average expression levels of mRNAs were modestly higher than those of lncRNAs (average 1.32-fold) in both young (3M) and old (18M) age groups ( Figure 2A). Most of the lncRNAs (69%) varied in length from 200 to 2000 nt, and the majority of the mRNAs (74.2%) ranged from 200 to 4000 nt ( Figure 2B). The expression levels of total transcripts, mRNAs, and lncRNAs were found to be similar among the age groups, from 3M to 18M, for all mouse chromosomes except for the Y chromosome (Supplementary Figure S1): the average expression levels of total transcripts and lncRNAs, but not mRNAs, increased slightly from 3M to 18M for the Y chromosome ( Figure 2C).

Aging-Related Expression Patterns of mRNAs and lncRNAs
Although the overall expression levels of transcripts were similar among the age groups, we expected that sets of transcripts would show changes during aging. To characterize aging-related mRNAs and lncRNAs in depth, we analyzed their expression patterns during testicular aging. For data preprocessing of the identified mRNAs and lncRNAs, we set the expression level criteria to select transcripts with FPKM ≥ 1 in at least one age group. This analysis yielded 13,797 mRNAs and 6230 lncRNAs. In addition, biological aging is accompanied with long and gradual accumulation of genetic damage [30]. To detect the result of accumulated damage and/or physiological changes, we investigate the continuous expression changes in testicular aging. We found that 1571 mRNAs and 715 lncRNAs showed a continuous increase or decrease, respectively, from 3M to 18M. Evaluation of the expression changes was based on the average values of three testicular samples in each age and statistical significance was not considered in this analysis ( Figure 3A and Supplementary Table S1).
To further analyze the expression patterns of transcripts exhibiting continuous gradual increases or decreases during testicular aging (3M, 6M, 12M, and 18M), we classified the mRNAs and lncRNAs by their degree of expression change across three consecutive agecomparison groups: 3M to 6M, 6M to 12M, and 12M to 18M. We considered transcripts with log 2 (|Fold change|) ≥ 1.2 to exhibit a "substantial expression change", as opposed to transcripts exhibiting a "slight change" over that period. From this analysis, we classified the mRNAs and lncRNAs into each eight groups representing increasing or decreasing expression patterns ( Figure 3B,C and Table 2). Total 1571 715 2286 1 The "+" mark represents a substantial expression change between age groups. 2 The numbers of aging-related mRNAs, lncRNAs, and total transcripts are indicated.  Table 2 and Supplementary Tables S2-S5. The classified expression patterns revealed several features. The majority (over 80%) of the transcripts [93.1% of mRNAs (1463 of 1571) and 85.8% of lncRNAs (614 of 715)] were classified as type 1, the members of which commonly showed very slight changes in expression level between the age groups ( Figure 3B,C and Table 2). The remaining transcripts (108 mRNAs and 111 lncRNAs) distributed to types 2-8, each of which included mRNAs and lncRNAs with increasing or decreasing patterns. The types can be divided into two groups, one comprising types 2, 3, and 5 and one comprising types 4, 7, and 8. The members of types 2, 3, and 5 showed a substantial expression change in a single age-comparison group (3M-6M, 6M-12M, or 12M-18M) ( Figure 3B,C and Table 2). For example, members of type 2 exhibited a substantial expression change only between 3M and 6M (28 mRNAs and 41 lncRNAs) ( Table 2). In contrast, members of types 4, 7, and 8 exhibited substantial expression changes in more than one age-comparison group ( Figure 3B,C and Table 2). For instance, members of type 4 showed substantial expression changes at both 3M-6M and 6M-12M groups (5 mRNAs and 7 lncRNAs) ( Table 2). Our analysis revealed that more transcripts showed a substantial expression change in only a single age-comparison group (174 transcripts) than in multiple age-comparison groups (36 transcripts). Interestingly, we observed that major expression changes during aging tend to be limited to the 6 month period (maximum 12M to 18M). We obtained information on gene ontology (GO) enriched in each classified expression pattern (Supplementary Tables S6 and S7). Spermatogenesis (GO:0007283) and single fertilization (GO:0007338) were found to be enriched in increasing expression patterns with substantial expression changes (Types 2-8) (Supplementary Table S6). GO terms with extracellular exosome (GO:0070062) and hydrolase activity (GO:0016787) were enriched in decreasing expression patterns having substantial expression changes (Types 2-8) (Supplementary Table S7).

Testis-Specific Expression of Aging-Related mRNAs and lncRNAs
To investigate the potential effect of aging on the transcriptome related to testis-specific features, such as spermatogenesis, we identified aging-related mRNAs and lncRNAs that are thought to be exclusively expressed in testes. In a previous study, we identified testisspecific mRNAs and lncRNAs from microarray analysis [16]. Using Ensembl Gene IDs, we herein compared aging-related transcripts found in the present study (1571 mRNAs and 715 lncRNAs of types 1-8) with the previously identified testis-specific transcripts. We identified at least 121 mRNAs and 25 lncRNAs as putative testis-specific aging-related transcripts (Supplementary Table S8). For the mRNAs, we used DAVID to perform functionenrichment analysis and predict their functional roles. We obtained information on gene ontology (GO) terms with p-value ≤ 0.05, and found that the 121 mRNAs were related to male reproductive components and processes, such as acrosomal vesicle (GO:0001669), male meiosis (GO:0007140), spermatogenesis (GO:0007283), and sperm-egg recognition (GO:0035036) (Supplementary Table S9). To provide more precise information on mR-NAs with testis specificity, we cross-checked testis-specific aging-related mRNAs with mouse ENCODE RNA sequencing data (http://genome.ucsc.edu/ENCODE, accessed on 13 October 2021) containing global tissue expression data. With the TissueEnrich tool [31], we narrowed down 121 putative testis-specific mRNAs to 77 mRNAs showing further testis-enriched expression (Supplementary Figure S3). Thus, 4.9% (77 of 1571) of the agingrelated mRNAs were predicted to be expressed in a testis-specific manner, suggesting that their expression changes could affect male reproductive functions during aging.

Potential Features of Aging-Related Transcripts Showing Changes between 3M and 18M
To further investigate the potential functional features of aging-associated transcripts in testes, we selected transcripts that were more stringently related to the terminal stage of aging. Unlike our analysis of transcripts whose expression levels showed continuous increases or decreases across all age groups ( Figure 3B,C and Table 2), we focused on transcripts that showed changes only between 3M (youngest age) and 18M (oldest age). These transcripts were chosen using DESeq2 with the criteria of a significant expression change (p-value ≤ 0.05 for three biological replicates) and a particular fold-change (log 2 (|Fold change|) ≥ 1) between the 3M and 18M age groups. This strategy selected 46 mRNAs and 34 lncRNAs (Figure 4), 55% of which (44/80) were members of types 2-8 and showed substantial expression changes ( Figure 3B,C and Table 2). The remaining transcripts (36/80) did not show substantial expression changes between the studied age-group pairs (type 1), but showed significant changes between 3M and 18M. We investigated the potential functional features of these transcripts (Table 3). For the mRNAs, 71% (33/46) were members of types 2-8 and showed substantial expression changes; in contrast, only 32.4% (11/34) of the lncRNAs were members of these types. We constructed a protein-protein interaction (PPI) network of the mRNAs predicted to be strongly associated with testicular aging [23]. Through a STRING-based PPI network analysis, we identified potential hub genes that were predicted to have more than one connection with other aging-related genes (Supplementary Figure S4). Notably, some of these hub genes are known to be involved in male reproduction. For example, progesterone receptor membrane component 1 (Pgrmc1) is implicated in the proliferation of male germ cells, as found in a conditional deletion study [32]. Cytochrome P450, family 7, subfamily a, polypeptide 1 (Cyp17a1) and cytochrome P450, family 11, subfamily a, polypeptide 1 (Cyp11a1) are cytochrome P450 family genes [33]. They encode proteins that are located in the inner mitochondrial membrane and endoplasmic reticulum, respectively, and participate in steroidogenesis. Finally, the expression level of kallikrein 1-related peptidase b27 (Klk1b27) is known to be regulated by testosterone in testis [34].

Potential Cis-Regulatory Targets of Aging-Related lncRNAs
Many functional lncRNAs are known to act in a cis-regulatory manner. To investigate the potential cis-regulatory targets of the 34 aging-related lncRNAs described above (Figure 4), we analyzed nearby genomic regions using Genomic Regions Enrichment of Annotations Tool (GREAT) [24]. This strategy identified 56 protein-coding genes that were located in genomic regions near the identified aging-related lncRNAs, and could thus be potential cis-regulatory targets. Indeed, five of these neighboring genes were found to show changes during aging: The mRNA expression levels of four decreased (Agtpbp1, Suds3, 2410089E03Rik, and Tex14) with age, whereas one increased (4921509C19Rik) (Figure 4). LncRNAs for three mRNAs (Suds3, 2410089E03Rik and Tex14) also showed decreases in level during aging. In contrast, the expression pattern (increase or decrease in level) of lncRNAs for the other two mRNAs (Agtpbp1 and 4921509C19Rik) was opposite to that of the mRNAs. We also examined the testis expression of the neighboring genes using the ENCODE database, and found that 21 of them are known to be expressed in testes (Table 4). Of these testicular genes, 13 exhibit testis-predominant (6 genes) or -specific (7 genes) expression. Finally, we performed functional enrichment tests with the potential testis-specific and -predominant cis-regulatory target genes. We found that three of the genes were related to the GO terms of cellular response to hormone stimulus (GO:0032870: 4921509C19Rik, and Gm14717) and protein phosphorylation (GO:0006468: 4921509C19Rik, Gm14717, and Tex14). Table 4. Testis-expressed potential cis-regulatory targets of aging-related lncRNAs identified by GREAT analysis.

Discussion
In the present study, we comprehensively investigated transcriptional changes associated with biological aging in adult mouse testis. In previous reports, sperm parameters and the morphology of testicular somatic cells were analyzed in consecutive age groups of various mammalian species, including humans [1][2][3][4]19], canines [21], and murine species [17]. Further studies investigated aging-related testicular changes in the context of epigenetic modification and gene expression [18,20]. Here, we used total RNA sequencing to reveal the profiles of transcripts, including both mRNAs and lncRNAs, in mouse testes of multiple age groups (3M, 6M, 12M, and 18M). We identified a total of 31,386 mRNAs and 9387 lncRNAs in the testes. Notably, the number of testicular lncRNAs identified by our total RNA sequencing analysis is similar to those found previously using microarray analysis (8265~14,256 lncRNAs) and single-cell RNA sequencing analysis (9431 lncRNAs) in mouse testes [16,[35][36][37]. Nonetheless, library preparation with TruSeq RNA Library Prep Kit v2, which is non-stranded, limited our analysis in identifying antisense RNAs exhaustibly. Further study using stranded RNA library may extend our profiling of antisense lncRNAs. Importantly, we newly identified 1571 mRNAs and 715 lncRNAs associated with testicular aging in mice (Figures 1 and 2).
Our genome-level analysis revealed that the total transcripts and lncRNAs expressed from the Y chromosome increased slightly during testicular aging (Figure 2). It was previously reported that histone modifications were altered on the peri-chromocenter, predicted to be putative sex chromosomes, of aged human spermatogenic cells [20]. It is possible that the observed increase in lncRNAs reflects transcriptional noise derived from cellular senescence [30,38]. Further studies are warranted to examine aging-related transcriptional changes at the chromosome level.
We investigated the expression patterns of the aging-associated transcripts (1571 mRNAs and 715 lncRNAs) in depth, seeking to elucidate the age(s) at which major transcriptional changes occur in the testis. The degree of differential gene expression during aging in mice has been found to vary by tissue type [39]. Here, we observed that the majority of aging-associated genes in testes showed only slight changes between the age groups. It should be noted that these changes were not statistically significant. Thus, the expression changes could be artefacts due to differences among animals and/or deviation between RNA sequencing analyses. Alternatively, this may represent the nature of aging that shows gradual and slight expression changes hard to detect experimentally. To confirm the gene expression changes observed in this study requires further investigation. In this regard, gene expression analysis of specific testicular cell types, instead of whole testes, is necessary, since cell-type specific expression differences may be hidden. Previously, gene expression analysis of spermatocytes in rats during aging revealed alteration of genes related to cell adhesion [40].
The analyses showing larger (herein termed "substantial") changes tended to be found in the 3M-6M and 12M-18M groups ( Table 2). This is in line with previous reports that the largest numbers of expression-altered genes were observed during periods considered middle-to-old age (12M-18M) in gonadal adipose tissue (GAT) and subcutaneous adipose tissue (SCAT) in mice [30]. Perhaps gene expression alterations that emerge from middle age could affect testicular function in late life. Transcriptional analysis of an age group older than 18M might be necessary to explore this possibility.
Testis-specific genes play important roles in male reproduction [5,6,16]. Interestingly, the testis contains the largest number of tissue-specific mRNAs and lncRNAs. We recently reported that the testis-specific lncRNA, Teshl, promotes the expression of genes on the Y chromosome and thereby regulates the offspring sex ratio in mice [17]. In the present study, we identified 121 mRNAs and 25 lncRNAs as being testis-specific aging-related transcripts. Function-enrichment analysis of these mRNAs revealed that some are related to male reproductive functions. Regarding the potential cis-regulatory targets of agingrelated lncRNAs, one of the candidates is Tex14. This gene is required for intercellular bridge formation in spermatogenic cells and essential for male mouse fertility [41]. The observed decrease in level of Tex14 in our analysis suggests a possible link of the gene to abnormal spermatogenesis during testicular aging. Another interesting novel candidate is 4921509C19Rik, which showed an increase in expression level during testicular aging. This gene has been predicted to play a role in sperm motility (https://www.uniprot.org/ uniprot/Q8C0X8, accessed on 13 October 2021). It should be noted that search of potential cis-regulatory targets was performed with mild criteria (1 Mb distance), requiring further investigation for cis-regulatory relationship. These testis-specific aging-related transcripts could be causative genes for aging and/or reflect male reproductive dysfunction due to testicular aging.
The decline of male reproductive hormone levels is a main cause of andropause, and aging in mammals is associated with a decrease in the testosterone level [2,4,19]. Here, we observed aging-related expression changes of genes related to steroid hormone synthesis and reactions, such as Cyp17a1, Cyp11a1, and Klk1b27 [33,34]. In the context of male reproductive aging, altered hormonal regulation is likely to affect the intrinsic regulation of gene expression in male germ cells.
Notably, most of the gene expression changes observed herein during testicular aging were slight and gradual. The majority of identified transcripts were members of type 1, showing small changes in expression level between the age groups. This is consistent with the gradual, physiological dysfunction seen in andropause, which contrasts with the more abrupt alterations seen during female menopause. The nature of aging appears to involve transcriptional changes of small proportions of genes, as shown in human kidney [42], human brain [43,44], monkey skeletal muscle [45] (~4%), fruit fly heart (~3%), human Achilles tendon [46], human retina [47,48], rodent brain [49], rodent skeletal muscle [50], rodent liver [51], and rodent heart (≤2%) [39,50,52,53]. In the present study, the agingrelated testicular transcripts accounted for 5.6% (2292) of the total identified transcripts (40,773). These observations suggest that the transcriptional landscape of aging should be elucidated in a sophisticated and precise way. In this regard, it could be important to investigate a possibility that some gene expression changes observed during aging are attributable to changes in cellular composition. In fact, a recent study showed that aging could affect differentiation of mouse SSCs, resulting in transcriptomic alterations [54].
In conclusion, we herein analyzed the transcriptional changes of mRNAs and lncRNAs during mouse testicular aging, and provide a comprehensive profile of the aging-associated transcripts. Most of the identified transcripts showed modest and gradual transcriptional changes. We classified the aging-related testicular transcripts into eight types based on their expression patterns. Further analyses provided additional insights, including potential features of aging-related mRNAs and potential cis-regulatory targets of aging-related lncRNAs. The present findings improve our understanding of the molecular mechanisms underlying testicular dysfunction in aging and should facilitate future investigations into the transcriptional signature of male reproductive aging.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/cells10112895/s1: Figure S1. Global expression level distribution of whole transcripts, mR-NAs and lncRNA in whole chromosomes. Expression levels are presented as log2(FPKM + 0.001); Figure S2. PCA plot between samples by DESeq2. (A) Aging-related mRNAs. (B) Aging-related lncR-NAs; Figure S3. Tissue expression of 77 testis-enriched genes identified among testis-specific agingrelated mRNAs by TissueEnrich, with mouse ENCODE expression data (BioProject, PRJNA66167); Figure S4. Protein-protein (PPI) network constructed by STRING database. Width of edge indicates the number of interactions detected in STRING and the color of each node presents each protein; Table S1. Aging-related mRNAs and lncRNAs showed a continuous increase or decrease from 3M to 18M. Class codes were identified with Cuffcompare (http://cole-trapnell-lab.github.io/cufflinks/ cuffcompare/, accessed on 1 June 2020): known (class code '='), novel (class code 'j'), intergenic (class code 'u'), intronic (class code 'i'), and antisense (class code 'x'); Table S2. Increasing mRNAs in classified expression patterns; Table S3. Decreasing mRNAs in classified expression pattern; Table S4. Increasing lncRNAs in classified expression patterns; Table S5. Decreasing lncRNAs in classified expression patterns; Table S6. Functional enrichment of mRNAs classified as increasing expression patterns. Gene ontology (GO) IDs and terms with p-value ≤ 0.05 and annotated ratio ≥ 5% are described. Category indicates GO aspects: molecular function (MF), cellular compartment (CC), and biological process (BP). GO enrichment data that do not satisfy criteria were omitted; Table S7. Functional enrichment of mRNAs as classified decreasing expression patterns. GO IDs and terms with p-value ≤ 0.05 and annotated ratio ≥ 5% are shown. Category indicates GO aspects; molecular function (MF), cellular compartment (CC), biological process (BP). GO enrichment data that do not satisfy criteria were omitted; Table S8. Putative testis-specific aging-related mRNAs and lncRNAs; Table S9. Functional enrichment of putative testis-specific aging-related mRNAs.