Statistics of clean reads alignment
After RNA sequencing and reads quality control, we obtained over 98 million 2×125 bp paired-end trimmed clean reads from each group, corresponding to an average of 12.3 Gb in size, of which over 79.4% were mapped to the Sus scrofa 10.2 assembly genome (Table 1). Approximately, 75.3%, 36.8%, 2.4%, 24.9%, 13.5%.24.7% of the clean reads were mapped to gene, coding, splicing, intron, no-coding and intergenic regions in Y-3 while 71.5%, 27.6%, 1.3%, 29.6%, 14.3%, 28.5% were mapped in Y-90 and 71.7%, 27.9%, 1.4%, 29.3%, 14.5%, 28.3% were mapped in Y-180 (Fig. S1). After aligning, we obtained 19,647 genes in nine Yorkshire spleens at 7, 90 and 180 days old in total.
Profile of distinct expressed genes
In order to understand how genes varied in spleen tissues of Yorkshire pigs at different development stages, we compared the expression levels of genes in the three age groups, Y-7, Y-90 and Y-180. Our ranking analysis (see Materials and Methods) yielded a list of 1,729 differential expressed genes, among which 1490 (617↑, 837↓) genes were in Y-90 vs Y-7, 1160 (556↑, 604↓) genes were in Y-180 vs Y-7 and 124 (58↑, 66↓) genes were in Y-180 vs Y-90. We found that the gene expression pattern in Y-90 and Y-180 showed a high degree of consistency (Fig. 1). Of the 1160 differently regulated genes in Y-180 vs Y-7, 957 were also detected in Y-90 vs Y-7. Only 27 differential expressed genes were at in common at trees stages, of which 11 gene were protein coding genes with characters, including IGKV-3, IGKV-7, IGKV-11, IGLV-2, IGLV-3, SLA-3, SLA-10, CXCL9, HDAC9, HSPA1A and MRC2. What was interesting was that some extremely differential expressed miRNAs were detected. For example, ssc-mir-199a-2 expressed in an extremely high level in Y-7 but was almost absent in Y-90 and Y-180 (↓4.7E-10, ↓4.7E-10) whereas ssc-mir-4332 expressed in an extremely high level in Y-7 and Y-180, but was absent in Y-90 (↓1.1E-8). In addition, ssc-mir-451 was down regulated both in Y-90 and Y-180 (↓0.2, ↓0.2).
GO and KEGG enrichment analysis of the distinct genes
After converting all the genes in pig into their corresponding H. sapiens Entrez gene IDs, 413, 27 and 497 differential expressed genes were remaining in Y-90 vs Y-7, Y-180 vs Y-90 and Y-180 vs Y-7. Go and KEGG enrichment analysis of these differential expressed genes were performed using Cytoscape. The overlapped up and down regulated genes and genes that shared the same enriched GO and KEGG ontology term(s) were shown in Circos plots (Fig. 2). Both up and down regulated gene expression profiles in Y-90 and Y-180 showed great similarity. Comparing with the up regulated genes, only a small number of down regulated genes that colored in light orange were unique to each of the three groups of samples.
In order to understand these enriched GO and KEGG ontology terms in detail, we extracted top 20 up and down regulated GO and KEGG clusters in the three groups and displayed them in the heatmaps (Fig. 3). No significant up or down regulated enriched term was found in Y-180 vs Y-90. As the age increased, various GO terms related to the immune system were accumulated in large numbers (Fig. 3a), such as GO: 0046649 and GO: 0002253 whereas the significant down regulated genes were almost enriched in cell replication and division (Fig. 3c), such as GO: 0006260 and GO: 0051301. For the up regulated genes enriched GO terms, 14 out of the top 20 were almost at equally significant level in both Y-90 vs Y-7 and Y-180 vs Y-7, among which the activation of lymphocytes and B cells and the regulation of various immune responses were at a relative higher significant level. The other six GO terms of the top 20 were all at a more significant level in Y-180 vs Y-7 than in Y-90 vs Y-7, which enriched in negative regulation of immune response and phagocytosis, such as GO:0050777 (negative regulation of immune response) and GO:0002683 (negative regulation of immune system process). For the down regulated GO terms, 18 of the top 20 were at almost the same level in Y-180 and Y-90.
With the age increased, disease infection and immune cell signaling pathways were enriched (Fig. 3b). Almost all terms had different expression levels in Y-90 and Y-180. In addition to Epstein-Barr virus infection and tuberculosis, the concentration of rheumatoid, staphylococcus and cytomegalovirus infection were all higher in Y-180. The concentration of antigen presentation and protein in endoplasmic reticulum and protein export were higher in Y-90, while the enrichment of chemokines, natural killer cells, T cells and other immune cells and intestinal immunity were higher in Y-180. The down-regulated pathway was basically the same between Y-90 and Y-180 (Fig. 3d). Cell cycle, alcoholism and DNA replication pathways were highly down-regulated in both groups. Signaling pathways including AMPK, mitophagy, adenine ribonucleotide biosynthesis, pyrimidine metabolism, insulin related PI3K-Akt and foxo, anemia, ferroptosis and fatty acid metabolism were less down-regulated. Down regulation of pyrimidine and fatty acid metabolism were higher in Y-180 vs Y-7 while PI3K-Akt signaling pathway was higher in Y-90 vs Y-7. Hematopoietic cell lineage and phagosome were only found in Y-90 vs Y-7.
Network analysis of known differentially expressed genes
To understand the differentially expressed gene patterns and the gene interactions involved in the significant expression modules at 7 days, 90 days and 180 days after birth in the Yorkshire spleens, we first used STEM to perform gene expression model analysis on the 1,729 genes and the results were shown in Fig. S2. A total of eight modules were obtained, of which three were significant and contained 1175 genes. Within the three significant modules, one included 666 genes that were up regulated from 7 days to 90 days and maintained the up regulation in 180 days, the other one involved 850 genes that were down regulated from 7 to 90 days and maintained the down regulation in 180 days, and the last one contained 52 genes that were up regulated from 7 days to 90 days and up regulated from 90 days to 180 days. Then, we used STRING to conduct PPI analysis on the 1175 module genes. A total of 756 interactions with combined scores no less than 0.97 were filtered to display using Cytoscape. Networks that combined no less than three edges, having average shortest pa We sequenced mRNA and lncRNA from spleen tissues of 7-day-old, 90-day-old and 180-day-old large white pigs, and obtained,,,, and lncRNA, among which,,,, and, were differentially expressed among the three age tissues.th length no less than three and with no less than three neighbors within distance two were finally displayed in Fig. 4. There was a total of 165 nodes and the center nodes genes were turned out to be CDK1, PCNA and PLK in Yorkshire spleen up and down regulation gene network.
To further capture how the relationships of the functional genes varied from one week to 180 days after birth, differentially expressed genes in Y-90 vs Y-7, Y-180 vs Y-90 and Y-180 vs Y-7 had been selected to analyze network enrichment using IPA (http://www.ingenuity.com/products/ipa) under default settings. The significant enriched networks of differentially expressed genes further supported the view of the similarity between Y-90 and Y-180 groups. Top one significant gene network of Y-90 vs Y-7, Y-180 vs Y-7 and Y-180 vs Y-90 were exhibited in Fig. 5. These networks were related to cellular function and maintenance, hematological system development and function, cell death and survival. Center genes of the networks of Y-90 vs Y-7 and Y-180 vs Y-7 were almost the same, such as CD40, REL and BNIP3L. Only SPP1 and HLA-A were involved in the gene network of Y-180 vs Y-90. In addition to that, top regulators in three groups were obtained, including BNIP3L, IL5, CD38, TGFβ1 and BCL3 in Y-90 vs Y-7, while BNIP3L, IL5, CD38, TGFβ1 and IL4 in Y-180 vs Y-7, and ERAP1, NLRC5 and IL2RG in Y-180 vs Y-90.
Expression profile of LncRNAs
After genome aligning, transcripts assembling and FPKM calculating (see Materials and Methods), we totally detected 219 known and 3219 putative lncRNAs, including 22 exonic_sense, 294 exonic_antisense, 232 intronic_sense, 232 intronic_antisense, 2435 intergenic, and 295 bidirectional lncRNAs (Fig. S3a). Using the same criteria for mRNA, our ranking analysis yielded a list of 64 differentials expressed LncRNAs that had RPKM value more than one, among which 51 (16↑, 35↓) were in Y-90 vs Y-7, 9 (4↑, 5↓) were in Y-180 vs Y-90, and 56 (20↑, 36↓) were found in Y-180 vs Y-7. We found that the lncRNA expression pattern in Y-90 and Y-180 also showed a high degree of consistency (Fig. S3b). Of the 51 different regulated lncRNAs in Y-90 vs Y-7, 43 were also found in Y-90 vs Y-7. This similarity for genes between Y-90 and Y-180 was also existed for lncRNAs. The 64-differential expressed lncRNAs were located on 13 genes, of which only three encode proteins with characters, those were SLA-3, SLA-8 and SLA-9. SLA-3 was up regulated in Y-90 vs Y-7 and Y-180 vs Y-7, then down regulated from the age of 90 days to 180 days. SLA-8 was up regulated in Y-180 vs Y-7. SLA-9 was up regulated in Y-180 vs Y-7 and Y-180 vs Y-90. The only one differential expressed lncRNA transcript in common for three groups was ENSSSCT00000001325 on SLA-3 gene.
In order to understand the characteristics of these 64 differentially expressed lncRNAs in Yorkshire spleen tissues at three time points and the biological processes involved, we predicted the target genes of the 64 lncRNAs and carried out KEGG pathways enrichment analysis of those genes at each time point. Enriched pathways were obtained from Y-90 vs Y-7 and Y-180 vs Y-7 which only involved in six genes including RRP7A, BMP6, NUP205, DDO, POLR2D, and COL5A2. Each of the pathways contained only one gene. In detail, from Y-90 vs Y-7, four genes enriched in six pathways, including RRP7A-hsa03008 (ribosome biogenesis), BMP6-hsa04340 and hsa04350 (hedgehog and TGF-β signaling pathway), NUP205-hsa03013 (RNA transport), DDO-hsa00250 and hsa04146 (alanine, aspartate and glutamate metabolism and peroxisome). From Y-180 vs Y-7, five genes were enriched in six pathways. Except for RRP7A, BMP6 and DDO genes and pathways in the above, POLR2D functions in five pathways including hsa00240 (pyrimidine metabolism), hsa00230 (purine metabolism), hsa03020 (RNA polymerase), etc, and COL5A2 functions in four pathways including hsa04512 (ECM receptor interaction), hsa04974 (protein digestion and absorption), hsa04510 (focal adhesion) and hsa05146 (amoebiasis) were detected.
Taking advantage of the mRNA expression profile in this study, we further analyzed the potential target differential expressed genes of the 64 lncRNAs. As a result, 32 of the 64 lncRNAs targeted 31 differentially expressed Ensembl genes, of which 23 genes had characters. In detail, 24, 24 and 9 lncRNAs potentially targeted 24, 21 and 6 Ensembl genes in Y-90 vs Y-7, Y-180 vs Y-7 and Y-180 vs Y-90, of which 18, 21 and 4 genes had characters (Table S2). Within the 23 characteristic genes, 17 were overlapped in Y-90 vs Y-7 and Y-180 vs Y-7 and with the same up & down regulated state. The four genes in Y-180 vs Y-90 were not unique. More than one lncRNA was found to target more than one gene, such as the target gene BMP6 and TRBC1, which were both targeted by two new lncRNAs in Y-90 vs Y-7 and Y-180 vs Y-7. Among the 23 target genes, SLA-3 and SLA-8 gene were acquired immunity related genes, SCARNA2 and SCARNA6 gene were members of small Cajal body-specific RNA family, NCAPD3 (Non-SMC Condensin II Complex Subunit D3) is a subunit of condensin II, a large protein complex involved in chromosome condensation, BMP6 (Bone morphogenetic protein 6) is a member of the TGFβ gene coded protein, ANP32B (Acidic Nuclear Phosphoprotein 32 Family Member B) functions in RNA polymerase binding, AHSP (Alpha-hemoglobin-stabilizing protein) is involved in hemoglobin assembly.
Real time PCR
Totally, six differential expressed genes and six LncRNA transcripts were randomly selected to perform RT-PCR in order to validate the sequencing data using GAPDH as the internal control gene for respectively. The sequencing and experiment result of these genes and lncRNAs were presented in Fig. 6. It was showed that there was a relatively good consistency between the sequencing and quantitative experiment results except for ALB gene and lncRNA transcript ENSSSCT00000022100.