Comprehensive analysis of circRNA expression profiles and circRNA-associated competing endogenous RNA networks in the development of mouse thymus.

The thymus plays an irreplaceable role as a primary lymphoid organ. However, the complicate processes of its development and involution are incompletely understood. Accumulating evidence indicates that non-coding RNAs play key roles in the regulation of biological development. At present, the studies of the circRNA profiles and of circRNA-associated competing endogenous RNAs (ceRNAs) in the thymus are still scarce. Here, deep-RNA sequencing was used to study the biological mechanisms underlying the development process (from 2-week-old to 6-week-old) and the recession process (from 6-week-old to 3-month-old) of the mouse thymus. It was found that 196 circRNAs, 233 miRNAs and 3807 mRNAs were significantly dysregulated. The circRNA-associated ceRNA networks were constructed in the mouse thymus, which were mainly involved in early embryonic development and the proliferation and division of T cells. Taken together, these results elucidated the regulatory roles of ceRNAs in the development and involution processes of the mouse thymus.

INTRODUCTION

The thymus is defined as a primary lymphoid organ with its inimitable role in T cell maturation, education and selection. It is the first formed lymphoid organ, but may also the most rapidly ageing tissue in the body. In response to the demand for large numbers of mature T cells, the thymus grows fast in the early stages of post‐natal life and reaches a plateau at 4‐6 weeks of age in the mouse. Afterwards, it starts the ageing process, of which the typical phenotype is thymic involution. The significantly decreased size and structural complexity of the thymus with age were manifested by the loss of thymocytes and thymic epithelial cells and the increase in the content of adipocytes. Actually, the processes about the rapid development and gradual involution of the thymus are still unclear.

MicroRNAs (miRNAs) are single‐stranded (average size 22 nt) non‐coding RNAs that regulate gene expression at the post‐transcriptional level. Through their broad effects on gene expression, miRNAs participate in the regulation of many cellular processes such as organismal development, , cellular differentiation and apoptosis , as well as mitochondrial metabolism. MiRNAs also take part in many aspects of T cell immunity such as T cell maturation, differentiation, activation and ageing. For example, some miRNAs including miR‐155, miR‐181c, miR‐9 and miR‐31 can regulate T cell activation by modulating the IL‐2 signalling pathway. Moreover, miR‐181a‐5p and miR‐195a‐5p expressed in thymic epithelial cells are involved in thymus involution, which directly target transforming growth factor β receptor 1 (Tgfbr1) and Smad family member 7 (Smad7) respectively. ,

Circular RNAs (circRNAs) are a class of single‐stranded closed RNA molecules without 3′‐poly (A) and 5′‐cap structures, which have been widely found in plants, animals and human beings. The cell‐type‐specific, tissue‐specific or developmental stage‐specific expression profiles of circRNAs suggest their regulatory functions in biological processes. , Numerous studies have reported that circRNAs are closely related to the tumorigenesis, , neurodegenerative diseases, cardiovascular diseases and immune diseases. , Some circRNAs might be involved in post‐transcriptional regulation by functioning as ‘sponges’ of miRNAs. Therefore, the circRNA‐miRNA‐mRNA networks may influence multiple biological pathways. At present, integrative analysis of circRNA‐associated competing endogenous RNAs (ceRNA) networks in the development and involution of the thymus are still scarce.

To gain further insight into the molecular events associated with thymic development and involution, RNA‐seq data were systematically analysed to identify aberrantly expressed circRNAs, miRNAs and mRNAs among mice thymuses at 2, 6 weeks and 3 months, respectively. In addition, the circRNA‐miRNA‐mRNA networks were constructed. This is the first comprehensive high‐throughput sequencing analysis of circRNA, miRNA and mRNA expression profiles in the thymus, which deepen our understanding of thymic development and involution.

MATERIALS AND METHODS

Animal tissues preparation

Male C57BL/6 mice were purchased from Beijing HFK Bioscience CO. LTD and maintained in a specific‐pathogen‐free environment. The mice were provided free access to standard diet until they met age requirements (2 weeks, 6 weeks and 3 months). Two biological replicates at each time point were used for sequencing. All the animal study protocols were approved by the Committee for the Ethics of Animal Experiments, Shenzhen Peking University‐The Hong Kong University of Science and Technology Medical Center (SPHMC; protocol number 2011‐004). The thymuses harvested from each group were placed in cryopreservation tubes and immediately immersed in liquid nitrogen (−196°C) for preservation.

RNA extraction

Total RNA from each thymus sample was isolated with the TRIzol Reagent kit (Invitrogen) and RNA degradation and contamination were monitored on 1% agarose gels. NanoPhotometer spectrophotometer (IMPLEN) was used to check the RNA purity. RNA concentration was measured with the Qubit RNA Assay Kit in Qubit 2.0 Flurometer (Life Technologies), and RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies). All the RNA samples had an RNA Quality Index ≥8.

RNA‐sequencing and miRNA‐sequencing

Details of the RNA‐seq and miRNA‐seq methods are described in Table S15.

Differential expression analysis

The differentially expressed circRNAs (DEcircRNAs), mRNAs (DEmRNAs) and miRNAs (DEmiRNAs) were identified using DESeq Software Packages (http://bioconductor.org/packages/release/bioc/html/DESeq.html). P value and log2FC were used to screen differential transcripts with P < .05 and |log2FC|>1.

Integrated analysis of circRNAs‐miRNAs‐mRNAs

CircRNAs were blasted against the circBase for annotation. Those cannot be annotated were defined as novel circRNAs. For circRNAs that have been annotated in circBase, the target relationship with miRNAs can be predicted by StarBase (v2.0). For novel circRNAs, three softwares Mireap, Miranda (v3.3a) and TargetScan (v7.0) were used to predict targets. For the prediction of mRNAs interacting with circRNAs and miRNAs, miRTarBase (v6.1) was used to predict mRNAs targeted by miRNAs sponge. The triple network was finally built based on the ceRNA theory and the resulting correlation of circRNAs‐miRNAs‐mRNAs can be visualized by Cytoscape 3.01.

KEGG enrichment analysis of differentially expressed genes

To further understand the underlying biological mechanisms and pathways of ceRNA‐related genes, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was implemented using the clusterProfiler R package.

Data access

All raw and processed sequencing data have been submitted to the NCBI Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/) under accession number GSE139653.

RESULTS

Changes in thymus tissues

The age‐related changes were observed in the thymuses derived from the 2‐week‐old, 6‐week‐old and 3‐month‐old male mice (Figure S1). Consistent with the literature, the mouse thymus grows rapidly after birth, reaches the maximal size by sexual maturity and then gradually involutes. The volume of the thymus increases significantly in the 6‐week‐old mouse compared with that in the 2‐week‐old mouse. As expected, the thymus shows remarkable atrophy in 3‐month‐old mouse compared with that in 6‐week‐old mouse.

Identification of differentially expressed genes

A total of 94 differentially expressed circRNAs (DEcircRNAs, 45 up‐regulated and 49 down‐regulated), 154 differentially expressed miRNAs (DEmiRNAs, 38 up‐regulated and 116 down‐regulated) and 1836 differentially expressed mRNAs (DEmRNAs, 1018 up‐regulated and 818 down‐regulated) were identified in the 6‐week‐old mice relative to the 2‐week‐old mice (Figure 1A,D,G). In the 3‐month‐old mice, 102 DEcircRNAs (54 up‐regulated and 48 down‐regulated), 79 DEmiRNAs (47 up‐regulated and 32 down‐regulated) and 1971 DEmRNAs (1005 up‐regulated and 966 down‐regulated) were observed relative to the 6‐week‐old mice (Figure 1B,E,H). Cluster analyses and heat‐maps about the expression of these circRNAs, miRNAs and mRNAs were conducted respectively (Figure 1C,F,I). Detailed information is listed in Table [Link], [Link], [Link], [Link], [Link], [Link].

FIGURE 1

Expression profiles of distinct RNAs. A‐C, Expression profiles of circRNAs. A, B, Volcano plot of different circRNAs in the mouse thymus. A, 2W vs 6W thymus. B, 6W vs 3M thymus. Red, green and black points represent circRNAs that were up‐regulated, down‐regulated and not significantly different in the mouse thymus respectively. x‐axis: log2 ratio of circRNA expression levels, y‐axis: false‐discovery rate values (−log10 transformed) of circRNAs. C, Cluster analysis of expression of circRNAs. Red and blue: increased and decreased expression, respectively. Expression profiles are shown similarly for (D‐F) miRNAs and (G‐I) mRNAs

Construction of the ceRNA network

According to the ceRNA hypothesis, competing endogenous RNAs (ceRNAs) members can regulate each other through competing for the same miRNA response elements (MREs). RNA transcripts communicate with the ceRNA language. The ceRNA networks in the mouse thymus were constructed based on our RNA‐seq data. As the Venn diagram shown in the Figure 2, we divided the differentially expressed transcripts (circRNAs, miRNAs and mRNAs) into two groups: (a) 2W/6W (+) 6W/3M (−): differential expression in 2W vs 6W group but not in 6W vs 3M group. We consider these transcripts in the group as the factors that play a vital role in the rapid development of the thymus. (b) 2W/6W (−) 6W/3M (+): no differential expression in 2W vs 6W group, but differential expression in 6W vs 3M group. The differential transcripts in this group are thought to be factors that play an important role in the early decline of the thymus.

FIGURE 2

Venn diagrams of differentially expressed transcripts. A, circRNA; B, miRNA; C, mRNA

A total of 40 circRNAs, 136 miRNAs and 523 mRNAs were significantly dysregulated in the 2W/6W (+) 6W/3M (−) group, which were selected to construct ceRNA networks (Tables S7 and S8). The ceRNA networks included both positive and negative links (Figure 3A,B). Figure 3A shows the decreased circRNAs, increased miRNAs and decreased mRNAs and Figure 3B shows the increased circRNAs, decreased miRNAs and increased mRNAs. In the 2W/6W (−) 6W/3M (+) group, 45 DEcircRNAs, 65 DEmiRNAs and 407 DEmRNAs were used to construct the ceRNA networks (Table S9 and S10). As mentioned above, the ceRNA networks also include two types of links in this group (Figure 4A,B). Figure 4A shows the decreased circRNAs, increased miRNAs and decreased mRNAs and Figure 4B shows the increased circRNAs, decreased miRNAs and increased mRNAs. Information about the differentially expressed transcripts used to construct the networks can be found in Tables [Link], [Link], [Link], [Link]. And the top15 DEcircRNAs between groups was shown in Tables 1 and 2.

FIGURE 3

CircRNA‐associated ceRNA networks. CeRNA networks were constructed based on identified circRNA–miRNA and miRNA–mRNA interactions. The networks include increased circRNAs‐decreased miRNAs‐increased mRNAs and decreased circRNAs ‐ increased miRNAs‐decreased mRNAs in 2W/6W (+) and 6W/3M (−) group. (A) down‐up‐down (B) up‐down‐up

FIGURE 4

CircRNA‐associated ceRNA networks. CeRNA networks were constructed based on identified circRNA–miRNA and miRNA–mRNA interactions. The networks include increased circRNAs‐decreased miRNAs‐increased mRNAs and decreased circRNAs ‐ increased miRNAs‐decreased mRNAs in 2W/6W (−) and 6W/3M (+) group. (A) down‐up‐down (B) up‐down‐up

TABLE 1

Top15 DEcircRNAs in the thymus between 2‐wk‐old and 6‐wk‐old mouse

CircRNA ID	T‐2W_RPM	T‐6W_RPM	log2(FC)	P value
novel_circ_002600	6282.377	2524.3654	−1.31539	5.80E‐08
mmu_circ_0000130	1739.733	411.06649	−2.08142	9.51E‐05
novel_circ_002601	3203.865	1161.2558	−1.46413	.000193
novel_circ_002657	637.3261	22.816464	−4.80388	.000193
novel_circ_005640	959.8191	167.44411	−2.51908	.000606
novel_circ_006249	0.001	414.93092	18.66251	.000805
mmu_circ_0000909	0.001	399.84332	18.60908	.001563
novel_circ_002221	596.1799	98.994722	−2.59032	.002377
novel_circ_004311	455.5065	0.001	−18.7971	.002433
novel_circ_002628	0.001	365.43356	18.47925	.002507
novel_circ_002392	0.001	342.61709	18.38624	.003053
novel_circ_004939	3891.912	1720.4442	−1.1777	.003094
novel_circ_001377	412.4453	26.680896	−3.95032	.003707
mmu_circ_0001023	0.001	357.70469	18.44841	.004018
novel_circ_001409	46.89112	517.79007	3.464981	.004374

Abbreviation: FC, fold change; RPM, reads per million mapped reads.

TABLE 2

Top15 DEcircRNAs in the thymus between 6‐wk‐old and 3‐mo‐old mouse

CircRNA ID	T‐6W_RPM	T‐3M_RPM	log2(FC)	P value
novel_circ_002601	1161.256	4700.284	2.017062	3.25E‐14
novel_circ_002602	1047.173	4015.974	1.939249	2.34E‐09
mmu_circ_0001100	22.81646	944.118	5.37082	4.86E‐07
novel_circ_002600	2524.365	5384.594	1.092917	5.23E‐07
novel_circ_004338	936.5854	97.4279	−3.265	.000121
novel_circ_003454	841.085	64.95194	−3.69481	.000436
novel_circ_001994	247.4868	1066.482	2.107436	.000486
novel_circ_002939	1138.069	219.7915	−2.37238	.000543
novel_circ_004395	521.6545	32.47597	−4.00565	.000775
novel_circ_001902	388.25	0.001	−18.5666	.000947
novel_circ_003367	978.3539	219.7915	−2.15422	.001128
novel_circ_004936	0.001	479.5993	18.87147	.001142
novel_circ_001285	2189.107	981.5215	−1.15725	.001571
mmu_circ_0000861	2337.969	879.166	−1.41105	.001759
novel_circ_004267	22.81646	387.0989	4.084555	.002271

Abbreviations: FC, fold change; RPM, reads per million mapped reads.

Functional enrichment analysis

Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis was used to determine the underlying biological mechanisms and pathways in the mouse thymus. Bubble charts were used to represent the top 20 pathways in 2W/6W (+) 6W/3M (−) group (Figure 5A,B) and 2W/6W (−) 6W/3M (+) group (Figure 5C,D), respectively. Each group included both positive and negative regulation. In the 2W/6W (+) 6W/3M (−) group, ceRNA gene‐related KEGG pathways were involved in pyruvate metabolism, lysine degradation, PI3K‐Akt signalling pathway, Ras signalling pathway, cytokine‐cytokine receptor interaction, insulin signalling pathway and so on (Tables S11 and 12). The ceRNA genes in the 2W/6W (−) 6W/3M (+) group were associated with biological processes including fatty acid biosynthesis and metabolism, regulation of lipolysis in adipocyte, pyruvate metabolism, insulin secretion, T cell receptor signalling pathway, regulation of autophagy, etc (Table S13 and S14).

FIGURE 5

KEGG pathways analysis in the thymus. KEGG pathway enrichment analysis of up and down‐regulated mRNAs with top twenty Enrichment score. (A) 2W/6W (+) 6W/3M (−) group: up‐regulated mRNAs. (B) 2W/6W (+) 6W/3M (−) group: down‐regulated mRNAs. (C) 2W/6W (−) 6W/3M (+) group: up‐regulated mRNAs. (D) 2W/6W (−) 6W/3M (+) group: down‐regulated mRNAs. The Y‐axis label represents the pathway and the X‐axis label represents the rich factor. The colour and size of the bubble represent enrichment significance and the amount of differentially expressed genes enriched in the pathway, respectively

DISCUSSION

The thymus is a vital immune organ in humans and other mammals which is the primary site of T cell development. The thymus provides a specialized environment to facilitate T cell maturation and generate an extremely diverse T cell repertoire. It also establishes its self‐tolerance to prevent autoimmune diseases. The thymus begins the involution process and apparently exhibits symptoms of age‐related changes at an accelerated rate relative to other organs. Previous studies have found that some miRNAs play an important role in thymic development and involution. , , However, the role of circRNA‐associated ceRNA networks in these processes is unknown. We designed this experiment to further investigate the molecular mechanisms underlining these complex processes of the post‐natal rapid development and subsequent gradual decline in the thymus.

Consistent with the literature, the mouse thymus undergoes rapid growth after birth. After reaching its maximum size at 6 weeks, the mouse thymus begins to decline. To identify the key molecules involved in these processes and understand their possible mechanisms, RNA‐seq was performed to compare the expression of circRNA, miRNA and mRNA in the thymus of 2‐week‐old, 6‐week‐old and 3‐month‐old mice.

By systematically analysing the sequencing data, we identified DEcircRNAs, DEmiRNAs and DEmRNAs in the 2W vs 6W group and the 6W vs 3M group (fold change ≥2.0 and P < .05), and mainly constructed the ceRNA networks of the 2W/6W (+) 6W/3M (–) group and the 2W/6W (−) 6W/3M (+) group. These circRNA‐associated ceRNA networks are potentially involved in the development and early decline of thymus. For instance, mmu_circ_0001023, novel_circ_003800 were identified as ceRNAs of mmu‐miR‐434‐5p, which targets Ids. The expression of Ids was higher in 6‐week‐old mice than that in 2‐week‐old mice. Ids are necessary for skeletal myogenesis, cardiac development and also play key roles in the differentiation and proliferation of T cells. , Moreover, mmu_circ_0000909 was found to be a ceRNA of mmu‐miR‐92a‐1‐5p, which targets KMT2D KMT2D is a major enhancer of H3K4me1/2 methyltransferase which is essential for early embryonic development in mice. There is evidence that KMT2D is required for cell‐type‐specific gene expression and cell differentiation in model systems for adipogenesis and myogenesis. Its mutation is closely related to the occurrence of T and B cell‐associated lymphoma. , In the 2W/6W (−) 6W/3M (+) group, mmu_circ_0000459, novel_circ_000862, novel_circ_001902 and novel_circ_002116 were identified as ceRNAs of mmu‐miR‐125a‐5p, which targets Erap1. Erap1 is an M1 zinc metalloprotease family member which has a significant impact on peptide processing function and the repertoire of peptides presented. In addition, the deregulation of ERAP1 can affect CD8+ T cell response. , Several researches have shown that the absence of ERAP1 expression in mice results in increased production of proinflammatory cytokines. The occurrence of inflammation is closely related to thymic ageing. , The down‐regulation of Erap1 in 3‐month‐old mice may be associated with the early thymic decline. Besides, mmu_circ_0001431, novel_circ_004960, novel_circ_005000 and novel_circ_005893 were identified as ceRNAs of mmu‐miR‐135a‐5p, which targets Celf1. Celf1, a member of the Celf family, belongs to RNA‐binding proteins. Previous research has indicated that Celf1 may affect the embryonic development process , and involve in the proliferation and division of T cells. Other links in these ceRNA networks are listed in Tables [Link], [Link], [Link], [Link].

In order to clarify ceRNA functions in mouse thymus, KEGG pathway enrichment analysis was used to determine the underlying biological pathways. The metabolism pathways such as pyruvate metabolism, fatty acid metabolism were significantly enriched pathways in both groups. Metabolism interconnects with signalling events to influence cell cycle, differentiation, cell death and immunological function. Yang et al reported that mTORC1 coordinated multiple metabolic programmes in T cells including glycolysis, lipid synthesis and oxidative phosphorylation to affect lymphocyte activation and fate decisions in adaptive immunity. In addition, Ras signalling pathway was found to be enriched in the 2W/6W (+) 6W/3M (−) group. Some evidence indicates that the activation of the Ras signalling pathway is critical for T cell development in the thymus. , , In the 2W/6W (−) 6W/3M (+) group, the expression of DEmRNAs related to the Insulin secretion pathway was significantly down‐regulated. Insulin is a key hormone which controls the metabolism of carbohydrates, proteins and lipids. Previous studies indicated that insulin receptor signalling controlled T cell proliferation, cytokine production and regulated adaptive immunity. , It suggests that insulin signalling pathway may modulate the thymus ageing through controlling the development of the thymocytes.

In summary, we demonstrated the circRNA‐associated ceRNA profiles of the mouse thymus at different time points by using deep‐RNA‐seq analysis for the first time. These results shed light on the possible contribution of the ceRNAs to the development and involution of the thymus.

CONFLICT OF INTEREST

The authors declare that there are no conflicts of interest.

AUTHORS CONTRIBUTIONS

WZ, JW and YZ designed the study. WL, NM, BY, YZ, YY, QL and T.YW. participated in the animal experiments including tissue collection and RNA/protein extraction. WL, NM, XF.C., JW and WZ analysed the data and wrote the paper. The final manuscript has been seen and approved by all authors.

Fig S1

Click here for additional data file.

Table S1

Click here for additional data file.

Table S2

Click here for additional data file.

Table S3

Click here for additional data file.

Table S4

Click here for additional data file.

Table S5

Click here for additional data file.

Table S6

Click here for additional data file.

Table S7

Click here for additional data file.

Table S8

Click here for additional data file.

Table S9

Click here for additional data file.

Table S10

Click here for additional data file.

Table S11

Click here for additional data file.

Table S12

Click here for additional data file.

Table S13

Click here for additional data file.

Table S14

Click here for additional data file.

Table S15

Click here for additional data file.