Effects of Castration on miRNA, lncRNA, and mRNA Profiles in Mice Thymus

Thymic degeneration and regeneration are regulated by estrogen and androgen. Recent studies have found that long non-coding RNAs (lncRNAs) and microRNAs (miRNAs) are involved in organ development. In this study, RNA sequencing (RNA-seq) results showed that ovariectomy significantly affected 333 lncRNAs, 51 miRNAs, and 144 mRNAs levels (p < 0.05 and |log2fold change| > 1), and orchiectomy significantly affected 165 lncRNAs, 165 miRNAs, and 208 mRNA levels in the thymus. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis showed that differentially expressed genes (DEGs) were closely related to cell development and immunity. Next, we constructed two lncRNA–miRNA–mRNA networks using Cytoscape based on the targeting relationship between differentially expressed miRNAs (DEMs) and DEGs and differentially expressed lncRNAs (DELs) analyzed by TargetScan and miRanda. Besides, we screened DEGs that were significantly enriched in GO and in ceRNA networks to verify their expression in thymocytes and thymic epithelial cells (TECs). In addition, we analyzed the promoter sequences of DEGs, and identified 25 causal transcription factors. Finally, we constructed transcription factor-miRNA-joint target gene networks. In conclusion, this study reveals the effects of estrogen and androgen on the expression of miRNAs, lncRNAs, and mRNAs in mice thymus, providing new insights into the regulation of thymic development by gonadal hormones and non-coding RNAs.


Introduction
As a primary immune organ, the thymus regulates adaptive immunity by producing naïve T lymphocytes [1]. Thymus degeneration is closely related to age, and the main age-related manifestations are thymus size becomes smaller, reduced thymic epithelial cells, and gradually decreased naive T lymphocyte output [2,3]. When naive T lymphocytes in the peripheral circulation system decrease with aging, effector-memory T lymphocytes will proliferate to maintain a constant number of total lymphocytes [4][5][6]. However, the reactivity of old lymphocytes to new infections is weak. Therefore, the reduction in the number of naïve T lymphocytes is considered to be one of the causes of the decline in resistance to infections and cancer in older adults [7,8]. Moreover, Gui found that the number of thymocytes in female and male mice decreased sharply at 3 months of age [9].
Gonadal hormone levels are also closely related to aging [10]. Estrogen and androgen can regulate thymus development, which has been widely reported [11,12]. Studies have found that pregnancy or estrogen injection can induce thymus atrophy and inhibit the expression of chemokines, such as

Animals and Sample Collection
For the experiment, 48 one-month-old CD1 mice (half male and half female) were purchased from the Laboratory Animal Center of Guangzhou University of Chinese Medicine (license key: SCXK (Yue) 2013-0034). The mice were fed in a specific pathogen-free environment (12/12 h light/dark cycle, 22-24 • C, 40-60% humidity) and randomly subdivided into ovariectomy group (F3x), female control group (F3), orchiectomy group (M3x), and male control group (M3) (12 mice per group). Mice were allowed to adapt for 7 days before the experiment. On the eighth day, the mice in the treatment groups were anesthetized for ovary or testicle removal, and the control groups received the same treatment but were not neutered. At 3 months old [9], the thymuses of each group were collected aseptically, frozen immediately with liquid nitrogen, and stored at -80 • C until use.

Total RNA Extraction and Sequencing
Referring to known research and economic factors, we mixed the tissues of each group equally for sequencing by a service provider (LC-BIO Bio-tech ltd, Hangzhou, China). According to the manufacturer's instructions, total RNA was extracted from mice thymus using Trizol reagent (Invitrogen, Carlsbad, CA, USA). A Bioanalyzer 2100 (Agilent, Santa Clara, CA, USA) and agarose gel electrophoresis were used to assess the quantity, purity, and integrity of total RNA. Total RNA with RNA integrity number greater than 7 was selected for further experiments. Then a TruSeq Small RNA Sample Prep Kit (Illumina, San Diego, CA, USA) was used to construct miRNA libraries according to the manufacturer's instructions. Finally, we performed single-ended sequencing on an Illumina HiSeq 2500 (Illumina, San Diego, CA, USA).
Ribosomal RNA (rRNA) was depleted from DNA-free RNA by an Epicentre Ribo-Zero Gold Kit (Illumina, San Diego, CA, USA). Then rRNA-free RNA was reverse-transcribed to create the final complementary DNA (cDNA) libraries by mRNA-Seq sample preparation kit (Illumina, San Diego, CA, USA). Finally, the libraries were paired-end sequenced on an Illumina HiSeq 4000 (Illumina, San Diego, CA, USA).

Transcript Assembly
Cutadapt [32] and FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) were used to remove reads containing adaptors and low-quality bases, and to verify sequence quality ( Figure 1). Then clean data were mapped to the mice genome using Bowtie 2 [33] and Tophat 2 [34]. Finally, StringTie [35] and Ballgown [36] were used to merge all transcriptomes from this experiment and estimate their expression levels by calculating fragments per kilobase per million reads.

Total RNA Extraction and Sequencing
Referring to known research and economic factors, we mixed the tissues of each group equally for sequencing by a service provider (LC-BIO Bio-tech ltd, Hangzhou, China). According to the manufacturer's instructions, total RNA was extracted from mice thymus using Trizol reagent (Invitrogen, Carlsbad, CA, USA). A Bioanalyzer 2100 (Agilent, Santa Clara, CA, USA) and agarose gel electrophoresis were used to assess the quantity, purity, and integrity of total RNA. Total RNA with RNA integrity number greater than 7 was selected for further experiments.

Transcript Assembly
Cutadapt [32] and FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) were used to remove reads containing adaptors and low-quality bases, and to verify sequence quality ( Figure 1). Then clean data were mapped to the mice genome using Bowtie 2 [33] and Tophat 2 [34]. Finally, StringTie [35] and Ballgown [36] were used to merge all transcriptomes from this experiment and estimate their expression levels by calculating fragments per kilobase per million reads.

lncRNA Identification
First, we removed all transcripts that overlapped with known mRNAs and were shorter than 200 bp. Then CPC [37], CNCI [38], and Pfam [39] were used to predict transcripts with coding potential ( Figure 6). All transcripts with CPC score <-1 and CNCI score <0 were deleted. Finally, the

lncRNA Identification
First, we removed all transcripts that overlapped with known mRNAs and were shorter than 200 bp. Then CPC [37], CNCI [38], and Pfam [39] were used to predict transcripts with coding potential.
All transcripts with CPC score <-1 and CNCI score <0 were deleted. Finally, the transcripts annotated as i, j, o, u, and x were retained. These letters represent transfrags falling entirely within a reference intron; potentially novel isoform; generic exonic overlap with a reference transcript; and unknown, intergenic, and exonic overlap with reference on the opposite strand, respectively.

RNA Expression and Functional Analysis
Differentially expressed lncRNAs (DELs), miRNAs (DEMs), and DEGs were screened by the R package Ballgown with log2 (fold change (FC)) > 1 or log 2 (FC) < -1 and statistical significance p < 0.05. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) in the DAVID online tool (https://david.ncifcrf.gov/summary.jsp) were used to analyze the functions of DEGs; p < 0.05 represents significant difference.

lncRNA-miRNA-mRNA Network Analysis
Recent studies have shown that lncRNAs and mRNAs containing the same miRNA binding site can regulate each other's expression levels by competitively binding miRNAs. To determine the interactions between DELs, DEMs, and DEGs after ovariectomy or orchiectomy, miRWalk, miRanda, RNAhybrid, and Targetscan were used for screening. Then, the intersections of co-expressed lncRNAs, miRNAs, and mRNAs in the 4 software programs were used to construct a lncRNA-miRNA-mRNA network. Finally, the lncRNA-miRNA-mRNA network was visualized using Cytoscape v3.7 software (https://cytoscape.org/). As competing endogenous RNA (ceRNA) analysis is based on miRNAs, the relationship between lncRNAs and mRNAs needs further investigation.

Cell Isolation and Culture
Thymocytes were isolated from 1-month-old mice. According to the manufacturer's instructions, thymocytes were obtained by centrifugation using the lymphocyte separation solution (Dakewe, Shenzhen, China). Red blood cell lysate (Solarbio, Beijing, China) was used to remove red blood cells. Then, cells were diluted with RPMI-1640 (Gibco, Waltham, MA, USA) containing 10% fetal bovine serum (Gibco, Waltham, MA, USA) and 0.2% penicillin / streptomycin (Invitrogen, Carlsbad, CA, USA) and culture for 24 h in 37 • C at 5% CO 2 to remove adherent cells. TECS were isolated from mice thymus in our previous experiments [40]. TECS were cultured in DMEM (Gibco, Waltham, MA, USA) containing 10% fetal bovine serum (Gibco, Waltham, MA, USA) and 0.2% penicillin / streptomycin (Invitrogen, Carlsbad, CA, USA) at 37 • C in 5% CO 2 . Finally, thymocytes and TECS were collected and immediately placed in liquid nitrogen and stored at −80 • C until use.

Quantitative Real-Time PCR Analysis
Total RNA in mice thymus, thymocytes and TECS were extracted with Trizol reagent (Takara, Kusatsu, Japan). First-strand cDNA was synthesized with a ReverTra Ace quantitative real time-PCR (qRT-PCR) RT Kit (Toyobo, Osaka, Japan), with random hexamers (for mRNAs and lncRNAs) and stem-loop RT primers (for miRNAs, purchased from RiboBio) according to the manufacturer's instructions. Primers of mRNAs and lncRNAs (Table S1) were designed by Primer Premier 5.0 software (Premier Biosoft International, USA), and miRNA primers were purchased from RiboBio (Guangzhou, China). β-actin (for mRNAs and lncRNAs) and U6 (for miRNAs) were used to normalize the data. qRT-PCR was performed by a Bio-Rad CFX96 Real-Time PCR system (Bio-Rad, Hercules, CA, USA) using SYBR Green Real-Time PCR Master Mix (Toyobo) following the manufacturer's instructions. All qRT-PCR assays were conducted in triplicate, and the relative levels were measured in terms of threshold cycle (Ct) and calculated using the formula 2 −∆∆Ct [41].

Analysis of Transcription Factor Binding Sites of DEGs
The putative promoter sequence was obtained from Browser1 of the UCSC genome. TESS software v6.0 was used to searched Position-weigh matrices in the TRANSFAC database [42]. The relative fraction line was set to 0.9. Next, the internal PERL script was used to conduct hypergeometric testing. p < 0.05 was defined as enriched transcription factor.

Identification of Transcription Factor-Related miRNAs and Their Joint Target Genes
According to the analysis of DEGs-related transcription factors, we used CircuitsDB (http: //biocluster.di.unito.it/circuits/) to obtain information about miRNAs and joint target genes related to these transcription factors. Cytoscape was used to build transcription factor-miRNAs-joint target genes networks

Expression Profiles of miRNAs, lncRNAs, and mRNAs in Thymus
To identify putative transcripts in the thymus, five mice each from the F3x group, F3 group, M3x group, and M3 group were mixed together for sequencing. In this experiment, we obtained a total of 380 million uniquely mapped lncRNA reads (Table 1). In addition, we got 56,673 unique miRNA readings in the Rfam database (Table 2, Supplementary Figure S1). In total, we identified 17,498 candidate lncRNAs, 1528 miRNAs, and 46,595 mRNAs in all chromosomes (Supplementary Tables  S2-S4). Venny analysis showed that 13,343 lncRNAs, 878 miRNAs, and 37,643 mRNAs were expressed in all groups, accounting for 76.3%, 57.5%, and 80.8% of the total readings, respectively (Figure 2a-c). Chromosome analysis showed that lncRNAs were unevenly distributed on the chromosome and were basically consistent with the trends of chromosome length (Figure 2d,e). Interestingly, the minimum length, mean length, median length, and N50 of lncRNAs were very close to those of mRNAs in this experiment ( Figure 2f, Table 3). In addition, the distribution of miRNAs showed that miRNA length was mainly concentrated in 20-24 nt (Figure 2g).

Identification of DELs, DEMs, and DEGs
Data analysis showed that ovariectomy significantly affected the levels of 333 lncRNAs, 51 miRNAs, and 144 mRNAs in mice thymus (p < 0.05, FC > 2 or FC < 0.5) (Figure 3a Tables S8-S10). To further explore the effects of gender on thymic degeneration, we analyzed differential RNAs in the F3x and M3x groups. The results showed that 416 lncRNAs, 68 miRNAs, and 225 mRNAs were remarkably differentially expressed after ovariectomy and orchiectomy (Supplementary Figure S2a-c, Supplementary Tables S11-S13). In addition, the levels of 136 lncRNAs, 129 miRNAs, and 124 mRNAs were remarkably different between the F3 and M3 groups (Supplementary Figure S2d-f, Supplementary Tables S14-S16). The analysis results show that most DEGs were upregulated after ovariectomy when compared with the F3 group, whereas most miRNAs and lncRNAs were downregulated (Figure 3a-c). It is worth noting that most DELs were upregulated after orchiectomy when compared with the M3 group (Figure 3d

Identification of DELs, DEMs, and DEGs
Data analysis showed that ovariectomy significantly affected the levels of 333 lncRNAs, 51 miRNAs, and 144 mRNAs in mice thymus (p < 0.05, FC > 2 or FC < 0.5) (Figure 3a (d-f) Volcano plots of DELs, DEMs, and DEGs, respectively, in orchiectomy and male control groups.
(g-i) Heat maps of specific DELs, DEMs, and DEGs, respectively, in ovariectomy and female control groups. (j-l) Heat maps of specific DELs, DEMs, and DEGs, respectively, in orchiectomy and male control groups. Each row represents one RNA. Green and red represent relatively lower and higher gene expression in castration group, respectively (p < 0.05 and |log2FC| > 1).
In order to verify our RNA-seq data, we randomly selected 3 DEGs, 3 DELs, and 3 DEMs for qRT-PCR analysis (Figure 4). The results showed that the differentially expressed RNAs had the same expression trends in qRT-PCR and RNA-seq (R = 0.747, p = 0.003). (d-f) Volcano plots of DELs, DEMs, and DEGs, respectively, in orchiectomy and male control groups. (g-i) Heat maps of specific DELs, DEMs, and DEGs, respectively, in ovariectomy and female control groups. (j-l) Heat maps of specific DELs, DEMs, and DEGs, respectively, in orchiectomy and male control groups. Each row represents one RNA. Green and red represent relatively lower and higher gene expression in castration group, respectively (p < 0.05 and |log2FC| > 1).
In order to verify our RNA-seq data, we randomly selected 3 DEGs, 3 DELs, and 3 DEMs for qRT-PCR analysis (Figure 4). The results showed that the differentially expressed RNAs had the same expression trends in qRT-PCR and RNA-seq (R = 0.747, p = 0.003).

Functional and Pathway Annotation
GO terms and KEGG pathway annotations of DEGs were enriched via DAVID. The enriched GO terms were divided into 3 categories: biological process (BP), cellular component (CC), and molecular function (MF). In this study, the BP enrichment of DEGs after ovariectomy includes extracellular matrix organization (p = 0.015), peptide transport (p = 0.024), extracellular structure organization (p = 0.040), and cellular cation homeostasis (p = 0.042) (Figure 5a). KEGG enrichment of DEG after ovariectomy showed no significant enrichment pathway. Meanwhile，the BP enrichment of DEGs after orchiectomy includes immune response (p = 0.00000054), response to virus (p = 0.000025), defense response (p = 0.000036), inflammatory response (p = 0.0017), regulation of nitric oxide biosynthetic process (p = 0.0091), positive regulation of defense response (p = 0.0094), cellular alkene metabolic process (p = 0.011), positive regulation of response to stimulus (p = 0.014), production of nitric oxide during acute inflammatory response (p = 0.015). KEGG enrichment of DEG after orchiectomy includes arachidonic acid metabolism (p = 0.014), viral myocarditis (p = 0.020), and intestinal immune network for IgA production (p = 0.045).

lncRNA-miRNA-mRNA Network Construction and Visualization
As a member of ceRNAs, lncRNA can affect the expression of miRNA target genes by adsorbing miRNA. In this study, we screened co-expressed genes from differentially expressed RNAs to construct ceRNA networks. As shown in Figure 6a, the female lncRNA-miRNA-mRNA network included 37 lncRNAs, 10 miRNAs, and 18 mRNAs. The male lncRNA-miRNA-mRNA network included 27 lncRNAs, 14 miRNAs, and 36 mRNAs (Figure 6b). The results of the node degree distribution analysis in the ceRNAs network showed that the slope of the power law distribution in female was −1.115 and R 2 = 0.723, and the slope in the male was -1.306 and R 2 = 0.795, which indicates that these two ceRNA networks have typical bio-network scale-free feature (Figure 6c,d).

lncRNA-miRNA-mRNA Network Construction and Visualization
As a member of ceRNAs, lncRNA can affect the expression of miRNA target genes by adsorbing miRNA. In this study, we screened co-expressed genes from differentially expressed RNAs to construct ceRNA networks. As shown in Figure 6a, the female lncRNA-miRNA-mRNA network included 37 lncRNAs, 10 miRNAs, and 18 mRNAs. The male lncRNA-miRNA-mRNA network included 27 lncRNAs, 14 miRNAs, and 36 mRNAs (Figure 6b). The results of the node degree distribution analysis in the ceRNAs network showed that the slope of the power law distribution in female was −1.115 and R 2 = 0.723, and the slope in the male was -1.306 and R 2 = 0.795, which indicates that these two ceRNA networks have typical bio-network scale-free feature (Figure 6c,d).

Expression of DEGs in Different Cells of Thymus
In order to further explore the function of DEGs in thymus development, we screened DEGs that were significantly enriched in GO and in the ceRNA network for verification. Then, qRT-PCR was used to detect the expression of these DEGs in thymocytes and TECS. As shown in Figure 7, H2-M2 was highly expressed in thymocytes. SLC7A2, CAV1, CACNB4, and PTGS2 were highly expressed in TECS. The expression trend of POLR3K in thymocytes and TECS was consistent.

Expression of DEGs in Different Cells of Thymus
In order to further explore the function of DEGs in thymus development, we screened DEGs that were significantly enriched in GO and in the ceRNA network for verification. Then, qRT-PCR was used to detect the expression of these DEGs in thymocytes and TECS. As shown in Figure 7, H2-M2 was highly expressed in thymocytes. SLC7A2, CAV1, CACNB4, and PTGS2 were highly expressed in TECS. The expression trend of POLR3K in thymocytes and TECS was consistent.

Expression of DEGs in Different Cells of Thymus
In order to further explore the function of DEGs in thymus development, we screened DEGs that were significantly enriched in GO and in the ceRNA network for verification. Then, qRT-PCR was used to detect the expression of these DEGs in thymocytes and TECS. As shown in Figure 7, H2-M2 was highly expressed in thymocytes. SLC7A2, CAV1, CACNB4, and PTGS2 were highly expressed in TECS. The expression trend of POLR3K in thymocytes and TECS was consistent.

DEGs Transcription Factor-related DEMs and Joint Target Genes prediction
We use CircuitsDB to predict the transcription factors of miRNAs associated with DEGs, and predict the joint target genes. We found that one transcription factor (PEA3), which was related to DEGs after ovariectomy, was associated with mmu-miR-135a and 15 joint target genes (Figure 9a). In addition, we found that 5 DEGs-related transcription factors are associated with DEMs after orchiectomy. PBX-1 was associated with 3 DEMs, including mmu-miR-342, mmu-miR-135a, and mmu-miR-135b, and 5 joint target genes (Figure 9b). HOXA4 was associated with mmu-miR-135b and 2 joint target genes (Figure 9c). IRF was associated with mmu-miR-15a and 2 joint target genes (Figure 9d). SREBP-1 was associated with mmu-miR-200b and 3 joint target genes (Figure 9e). NF-γ

Discussion
The thymus, as a central immune organ, is the main site for generation of naive T lymphocytes [43]. Thymus atrophy is considered to be a biomarker of immune system aging and is influenced by many factors, such as age, hormones, and gender [44]. Studies have extensively demonstrated that ovariectomy or orchiectomy can result in transient thymic regeneration. However, the mechanisms by which gonadal hormones regulate thymic development remains unclear [17,18]. In addition, the role of non-coding RNA in regulating thymic development is revealed. In this study, 333 DELs, 51 DEMs, and 144 DEGs, and 165 DELs, 165 DEMs, and 208 DEGs were identified after ovariectomy and orchiectomy, respectively. In order to identify the accuracy of the RNA-seq data in this experiment, we randomly selected some differentially expressed RNAs for qRT-PCR identification. The results showed that DELs, DEMs, and DEGs had the same expression trends in qRT-PCR and RNA-seq (R = 0.747, p = 0.003), which parts supported the high reliability of our RNA seq data.
Next, GO and KEGG in DAVID were used to enrich the function of DEGs after ovariectomy and orchiectomy. DEGs after ovariectomy were mainly concentrated in extracellular structure organization, peptide transport, cellular cation homeostasis, and extracellular matrix organization in BP. Upregulation of peptide transport may help regulate lymphocyte proliferation [45]. Besides, cellular cation homeostasis is necessary for normal life activities of cells. For instance, K + is used to compensate negative charges and activate protein translation processes [46]; Mg 2+ is an important component of cell proliferation-related cofactors [47]. It is worth noting that the thymic extracellular matrix can increase thymocyte output in vivo and promote TEC differentiation in vitro [48]. Meanwhile, DEGs after orchiectomy were mainly concentrated in immune response and inflammatory response in BP. Castration or androgen deficiency may result in thymus enlargement, which may be related to androgen inhibiting immature thymocyte proliferation and accelerating its apoptosis [49]. In addition, androgens can affect the production of immune-related cytokines in the thymus. Upregulation of TGF-β and activation of TGF-β-induced Smad signaling pathway induces accumulation of connective fibers in the thymus [50]. As a proinflammatory cytokine, IL-6 can also stimulate the proliferation of thymocytes in vitro [51]. Besides, Oxidative stress affects the homeostasis of the immune system, causing the irreversibility of thymus degradation [52].
lncRNAs can participate in proliferation and aging physiological processes in different ways [53,54]. With the development of non-coding RNA research, the theory of ceRNA was proposed. In this theory, acting as key nodes, miRNAs can bind to lncRNAs to regulate the expressions of other non-coding RNAs or target genes with the same binding sites. In this experiment, we hypothesized that lncRNA is involved in regulating thymus development after mice castration through ceRNA. Next, we performed ceRNA analysis on the differentially expressed RNAs between castration and Transcription factor-miRNAs-joint target genes networks.
(a) Transcription factor-miRNAs-joint target genes network after ovariectomy. (b-f) Transcription factor-miRNAs-joint target genes network after orchiectomy. The purple triangles, brown circles, and green diamonds represents transcription factors, miRNAs, and joint target genes, respectively.

Discussion
The thymus, as a central immune organ, is the main site for generation of naive T lymphocytes [43]. Thymus atrophy is considered to be a biomarker of immune system aging and is influenced by many factors, such as age, hormones, and gender [44]. Studies have extensively demonstrated that ovariectomy or orchiectomy can result in transient thymic regeneration. However, the mechanisms by which gonadal hormones regulate thymic development remains unclear [17,18]. In addition, the role of non-coding RNA in regulating thymic development is revealed. In this study, 333 DELs, 51 DEMs, and 144 DEGs, and 165 DELs, 165 DEMs, and 208 DEGs were identified after ovariectomy and orchiectomy, respectively. In order to identify the accuracy of the RNA-seq data in this experiment, we randomly selected some differentially expressed RNAs for qRT-PCR identification. The results showed that DELs, DEMs, and DEGs had the same expression trends in qRT-PCR and RNA-seq (R = 0.747, p = 0.003), which parts supported the high reliability of our RNA seq data.
Next, GO and KEGG in DAVID were used to enrich the function of DEGs after ovariectomy and orchiectomy. DEGs after ovariectomy were mainly concentrated in extracellular structure organization, peptide transport, cellular cation homeostasis, and extracellular matrix organization in BP. Upregulation of peptide transport may help regulate lymphocyte proliferation [45]. Besides, cellular cation homeostasis is necessary for normal life activities of cells. For instance, K + is used to compensate negative charges and activate protein translation processes [46]; Mg 2+ is an important component of cell proliferation-related cofactors [47]. It is worth noting that the thymic extracellular matrix can increase thymocyte output in vivo and promote TEC differentiation in vitro [48]. Meanwhile, DEGs after orchiectomy were mainly concentrated in immune response and inflammatory response in BP. Castration or androgen deficiency may result in thymus enlargement, which may be related to androgen inhibiting immature thymocyte proliferation and accelerating its apoptosis [49]. In addition, androgens can affect the production of immune-related cytokines in the thymus. Upregulation of TGF-β and activation of TGF-β-induced Smad signaling pathway induces accumulation of connective fibers in the thymus [50]. As a proinflammatory cytokine, IL-6 can also stimulate the proliferation of thymocytes in vitro [51]. Besides, Oxidative stress affects the homeostasis of the immune system, causing the irreversibility of thymus degradation [52].
lncRNAs can participate in proliferation and aging physiological processes in different ways [53,54]. With the development of non-coding RNA research, the theory of ceRNA was proposed. In this theory, acting as key nodes, miRNAs can bind to lncRNAs to regulate the expressions of other non-coding RNAs or target genes with the same binding sites. In this experiment, we hypothesized that lncRNA is involved in regulating thymus development after mice castration through ceRNA. Next, we performed ceRNA analysis on the differentially expressed RNAs between castration and control groups and constructed two scale-free biological networks of ceRNAs. qRT-PCR preliminarily confirmed the expression of these genes in the ceRNA network were consistent with RNA-seq. In addition, the scale-free distribution of these two ceRNA networks indicates that there are important nodes in the regulation of gonadal hormones on thymus development.
In order to further clarify the function and expression site of DEGs in the ceRNA network. we screened DEGs (H2-M2, POLR3K, SLC7A2, CAV1, CACNB4, and PTGS2 that were significantly enriched in GO and in the ceRNA network to verify their expression in thymocytes and TECS. Cav1 is a channel for Ca2 + influx required for T cell development [55]. Cyclooxygenase 2 (COX-2) is encoded by the PTGS2 gene, which inhibits T cell proliferation by regulating NF-kappaB [56]. This experiment found that Cav1 and PTGS2 were significantly overexpressed in TECS, which may be related to the TEC providing a microenvironment for thymocyte development [23]. The qRT-PCR showed that H2-M2 and POLR3K were highly expressed in thymocytes, and SLC7A2 and CACNB4 were highly expressed in TECS. However, the roles of these genes in the thymus have not been reported. These genes may be worthy of further study.
In order to explore the regulatory mechanism of DEGs, we predicted the causal transcription factors of DEGs through enrichment experiments. We found that the binding sites of KROX, Lmo2, PEA3, STAT1, and STAT4 were significantly enriched in downregulated DEGs after ovariectomy, while AP-1, AP-2rep, c-Myc, E2F-1, DP-2, ETF, MOVO-B, Pax, and ZF5 binding sites were significantly enriched in upregulated DEGs. Meanwhile, the binding sites of NF-kappaB, Pbx-1, SMAD, and SREBP-1 were significantly enriched in downregulated DEGs after orchiectomy, while the binding sites of Hoxa4, ICSBP, IRF, ISGF-3, MOVO-B, NF-γ, Oct-4, and TEF-1 were significantly enriched in upregulated DEGs. The main function of Krox-20 is to inhibit the c-Jun NH2-terminal protein kinase-c-Jun pathway, whose activation is necessary for both proliferation and death [57]. Overexpression of LMO2 leads to delayed development of thymus progenitor cells and increases the number of T lymphocytes in lymphoid organs [58]. As members of the STAT family, TAT1 and STAT4 mainly promote the differentiation of Th1 cells [59]. AP-1 is one of the most important target members in the MAPK pathway, which can regulate cell division and apoptosis [60]. AP-2 rep is a transcriptional repressor of AP-2 and plays an important role in embryo development [61]. c-Myc makes cells sensitive to TNF stimulation and promotes apoptosis of thymocytes [62]. E2F-1 regulates apoptosis and proliferation by stabilizing p53 tumor suppressor [63]. The NF-kappaB pathway is considered to be a typical transduction pathway of proinflammatory cytokines represented by the IL-1 and TNF receptor families, which was required for T cell proliferation and emigration [64]. Smad is at the core of the TGF-β pathway [65]. HOXA4 expression is inversely related to cell cycle, metastasis, and Wnt signaling pathway, and its overexpression can inhibit tumor cell proliferation, migration and invasion [66]. Oct-4 belongs to the POU domain transcription factor family that are involved in regulating the proliferation and differentiation of various cells in tissues [67]. The function of some transcription factors in thymus development has been partially verified, but the role of PEA3, DP-2, ETF, MOVO-B, Pax, ZF5, Pbx-1, SREBP-1, ICSBP, IRF, ISGF-3, MOVO-B, NF-γ, and TEF-1 in regulating thymus development needs further investigation.

Conclusions
In summary, we identified the expression profiles of lncRNAs, miRNAs, and mRNAs in mice thymus from control groups and ovariectomy or orchiectomy group, and constructed two scale-free ceRNA networks. Bioinformatics analysis showed that DEGs were mainly enriched in extracellular matrix organization and immune response-related physiological functions. This study provides new insights into the effects of estrogen and androgen on thymic development, and the specific mechanisms need to be further explored.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4425/11/2/147/s1, Figure S1: Reading miRNAs mapped to the reference genome; Figure S2: Volcano plots of estrogen-and androgen-related co-expression RNAs in mice thymus; Table S1: Specific primer sequences for qRT-PCR; Table  S2: List of expression of lncRNAs in all groups; Table S3: List of expression of miRNAs in all groups; Table S4: List of expression of mRNAs in all groups; Table S5: List of DELs in mice thymus between ovariectomy and female control groups; Table S6: List of DEMs in mice thymus between ovariectomy and female control groups; Table S7: List of DEGs in mice thymus between ovariectomy and female control groups; Table S8: List of DELs in mice thymus between orchiectomy and male control groups; Table S9: List of DEMs in mice thymus between orchiectomy and male control groups; Table S10: List of DEGs in mice thymus between orchiectomy and male control groups; Table S11: List of DELs in mice thymus between ovariectomy and orchiectomy groups; Table S12: List of DEMs in mice thymus between ovariectomy and orchiectomy groups; Table S13: List of DEGs in mice thymus between ovariectomy and orchiectomy groups; Table S14: List of DELs in mice thymus between female control and male control groups; Table S15: List of DEMs in mice thymus between female control and male control groups; Table S16: List of DEGs in mice thymus between female control and male control group.

Conflicts of Interest:
The authors declare no conflict of interest.