Blood Transcriptome Analysis reveals Age-associated changes in Expression Profile of Immune-Related Gene in Golden snub-nosed Monkey ( Rhinopithecus )

Background: Golden snub-nosed monkeys ( Rhinopithecus roxellana ) are endangered species of monkeys found in China. In this study, we provided the blood transcriptome sequences of golden snub-nosed monkeys obtained using RNA-Seq technology. The genomic annotation of these monkeys was useful to identify the polymorphisms and subpopulations, in order to understand age-related changes of immune system. This data may provide a valuable resource for further genetic and genomic research of golden snub-nosed monkeys. Results : 57.31 Gb high-quality sequencing data were obtained. The clean data of each sample were >5 Gb, and 86.17% to 94.48% of the reads of each sample could be compared to reference genome of snub-nosed monkey. After assembly, we obtained 24,992 genes, including 3,917 new genes. Many genes were up-regulated or down-regulated with age. Compared to the young group , there were 76 differential genes in adult group of R. r. roxellana , including 68 up-regulated and 8 down-regulated genes. While, compared with the adult group, there were 58 differential genes, including 25 up-regulated genes and 23 down-regulated genes in the old group of R. r. roxellana . In R. r. qinlingensis , compared with the young group, 117 differential genes were obtained, including 34 up-regulated and 83 down-regulated genes. Functional enrichment analysis indicated that the up-regulated genes were mainly related to innate immune response and T-cell activity, while the down-regulated genes were mainly involved in B-cell activity, suggesting that immune competence of adult group increased gradually compared to young group. However, the adaptive immune function declined gradually in the old group. Conclusions: Our findings can contribute to understanding of molecular mechanisms of age-related changes of immune system, which will provide a foundation for future studies of snub-nosed monkey.

snub-nosed monkeys were widely distributed in the central and southern China, including high altitude areas [3,4]. However, due to habitat changes and human activities, their population has been extinct in many areas. At present, the distribution of these monkeys is limited to mountain forests in three isolated temperate regions of China [5].
RNA sequencing (RNA-Seq) using NGS technology is becoming the predominant tool for large scale gene expression analyses. The transcriptome is the sum of expression of genes in the cell [16]. The spatiotemporal patterns of gene expression and transcription profiles are different at different stages of growth [17]. Transcriptions can reveal the dynamic changes of gene expression at the genomic level, including comprehensive details such as heredity information and influence of environmental factors [19][20][21]. Blood is the first line of defense of human body. Currently blood transcriptome has been performed in several mammalian species to detect their immune-expression profile, including in giant panda and pigs [22][23][24][25][26]. Studies have shown that the development and function of immune organs were different in different stages of growth. Age is associated with altered immune functions [27]. Age-related decline of immune function, referred to as "immunesenescence", is thought to contribute to the increased incidences of infectious disease [28,29]. However, the blood transcriptome data of these monkeys at different stages of growth are not available. Moreover, understanding the dynamics of age-related changes of immune system is important for the conversation of these monkeys.
Based on high-throughput sequencing technology, we sequenced the blood transcripts of golden snub-nosed monkey at different stages of growth. Subsequently, we filtered the off-line data to obtain high-quality sequences (Clean Data) and aligned it with the reference genome of snub-nosed monkeys. Simultaneously, we analyzed and explored variable splicing, optimized gene structure and found novel genes. The expression levels of these genes in each group were analyzed to understand the functions and pathway of differentially expressed genes (DEGs). It was useful to annotate the genome of golden snub-nosed monkey to identify polymorphisms and determine subpopulations. This helped us understand the age-related changes occurring in immune system, which provided a valuable insight for further genetic and genomic research of these monkeys.
Specifically, we aimed to (i) perform a functional annotation of the blood transcriptome of golden snub-nosed monkey and (ii) understand the changes occurring with age in immune-related gene expression profile. The data generated in this study may provide valuable resources for further genetic study of golden snub-nosed monkey.

Sequencing and Assembly
In this study, 57.31 Gb clean data were obtained from 9 transcriptome sequencing libraries. The clean data of each sample exceeded 5Gb. Detailed results are listed in Table 1. The sequencing data of the samples were of high-quality. The number of clean data ranged from 19,923,730 to 26,157,724. The Q30 was higher than 85.79%. These results suggested that the quality of sequencing data was good for subsequent analysis.

Transcriptome Data and Reference Genome Sequence alignment
Clean reads of each sample were sequenced, assembled and quantified with the designated reference genome using HISAT2 and StringTie. The comparative efficiency of nine transcriptome libraries ranged from 86.17% to 94.48%. At least 36,987,254 mapped reads were obtained. Among them, 38,503,090 (80.68%) and 38,886,464 (89.10%) had unique alignment loci with genomic sequences (Table S1).
The ratio of reads of different regions to the reference genome was shown in Table S2. Reads from mature RNA were matched to exon regions, and nine transcript libraries were matched to exon regions from 47.10% to 76.81%. As variable splicing of the RNA precursor can cause intron retention, some reads could match the introns, and these libraries matched introns from 14.40% to 76.81%.

Transcriptome Annotation
Trinity software was used to predict the open reading frame (ORF) of the transcripts. To further understand the blood transcriptome of the golden snub-nosed monkey, we compared the predicted transcripts with 8 protein databases using BlastP and BlastX. We found that 24,102 (99.8%), 19,705 (81.6%), 9,098 (9.83%), 14,345 (59.4%) and 19,685 (81.5%) transcripts showed significant similarity matches with the Nr, GO, KEGG, COG, and PFAM databases, respectively ( Table 2).
In GO database, a total of 19705 transcripts were annotated, of which 12084 transcripts were larger than 1000 bp. GO analysis of the annotated information can be divided into three categories: biological processes, cellular components and molecular functions ( Figure 1). Among them, there were 9038 transcripts involved in metabolic process and 1233 transcripts involved in reproductive process. Whereas, there were 5139 transcripts involved in developmental process and 2154 transcripts involved in immune process.
In COG database, a total of 6668 (27.6%) transcripts were annotated, and 4627 transcripts were more than 1000 bp ( Figure 2). General function prediction represented the largest group, followed by translation, replication, recombination and repair, transcription, signal transduction mechanism, and post-translational modification.
In KEGG database, there were 14,345 transcripts (59.4%) with corresponding annotation information, and 9,239 transcripts of which were >1000 bp. Additionally, KEGG pathway analysis showed that transcripts were assigned to 330 pathways. Typical pathways included cancer, PI3K-AKT signaling, HTLV-I infection, MAPK signaling, actin cytoskeleton regulation, Epstein Barr virus infection, Rap1 signaling and Ras signaling pathway. In Nr database, 24102 transcripts (59.4%) were obtained with corresponding annotation information, and 14427 transcripts were >1000 bp.

Analysis of immune-related genes
The abundance of transcripts was expressed as FPKM. There were 2154 transcripts related to immune process in the blood of golden snub-nosed monkey.
Based on the quantitative results of expression value, 47 genes (Table S3) expressing >100 FPKM value were identified, and the expression value of these highly expressed immune-related genes varied from 107 to 6667. Compared with Nr database, the highest expression of these immune genes was for GOGO-C*0203a, ferritin heavy chain, major histocompatibility complex (MHC) class-I antigen, protein S100-A9, beta-actin, thioredoxin interacting protein, actin cytoplasmic 1.
Additionally, we also analyzed the MHC gene-family expressed in the blood transcriptome in golden snub-nosed monkey. Among them, Aime-128 of classical MHC class genes was the most highly expressed (Table S3). Comparing with Nr protein database, these MHC genes could be annotated in the database including MHC class-I and class-II proteins. Moreover, at least four subtypes of MHC class I, and ten subtypes of MHC class II were expressed in the blood transcription group of golden snub-nosed monkey.

New gene analysis
Comparing with the reference genome, we found new genes using StringTie software. A total of 3814 new genes were found, of which 3004 genes could be compared to 8 protein databases ( Table 3). The new genes had 28% homology to human (823 transcripts), 17% to cynomolgus monkey (506 transcripts) and 12% to rhesus monkey (350 transcripts).
In these new genes, there were 2106 transcripts that can be annotated to GO database. Of which, 1263 transcripts were involved in metabolic process, 827 in developmental process and 87 in reproductive process. Additionally, according to GO results, 322 transcripts of these new genes were involved in the immune process. New genes belonging to the immune system were mainly involved in antigen processing and presentation, natural killer cell-mediated cytotoxicity, hematopoietic cell lineage, complement and coagulation cascade and NOD-like receptor signaling pathway.

The SNP detection
The reads of 9 golden Snub-nosed monkey were mapped to the reference genome of snub-nosed monkeys and the SNP was identified ( Figure 3). The number of SNP site ranged from 85,522 to 196,836, as shown in table 4. Among them, the distribution of SNPs in gene region varied from 75,053 to 149,387, while the total number of SNP sites in intergenic region ranged from 20,586 to 45,597.
The numbers of transitions and transversions at the SNP sites varied from 72.45% to 75.06% and 24.94% to 27.44% respectively. The SNP converts A to G and T to C most frequently, followed by G to A and C to T ( Figure 4).

Intercorrelations among the groups
According to the correlation thermograms of the nine samples ( Figure 5), they were divided into two subspecies, RRM1 RRM2, RRM3, RRM4 and RRM5 belonged to R. r. roxellana; while RRM6, RRM7, RRM8 and RRM9 belonged to R. r. qinlingensis. Further, these subspecies were divided into different groups according to their age ( Table 5). The differentially expressed genes of subspecies were compared to each other.

Identification of DEGs in golden snub-nosed monkey subspecies
DEGs between the subspecies of the golden snub-nosed monkey were analyzed. The results showed that there were 12 DEGs between the subspecies (Table S4). As per DEGs, there were 8 up-regulated genes and 4 down-regulated genes in the R. r. qinlingensis compared with the R.r. roxellana.
Compared to Nr protein database, the up-regulated genes included dynactin subunit 3, UPF0184 protein C9orf16, and down-regulated genes included histone H3.3, AT-rich protein-3 in the interaction domain, amongst other. The expression of DEGs was low, indicating that there was little difference in gene expression between subspecies of golden snub-nosed monkeys.

Identification of gender-related DEGs
In order to determine the gene expression patterns between the genders, we investigated DEGs between male and female of golden snub-nosed monkey. The results showed that there were 13 DEGs between different genders (Table S5). Among DEGs, there were 8 up-regulated genes and 5 down-regulated genes in the males (RRM1, RRM2, RRM3) compared with the females (RRM4, RRM5) in R. r. roxellana. However, a total of 176 DEGs were identified between males (RRM6, RRM9) and females (RRM7, RRM8) in R. r. qinlingensis (Table S6). Compared with the males, 89 genes were upregulated, and 87 genes were down-regulated in the females. Both GO and KEGG enrichment analysis showed that DEGs were mainly involved in reproductive processes and immune response.

Analysis of DEGs in young and old groups of R. r. roxellana
In young and old groups, 76 DEGs were identified (Table S7). Compared to the young group, 68 genes were up-regulated, and 8 genes were down-regulated in the old group. GO enrichment analysis showed that DEGs were mainly involved in immune responses, such as Actin, MHC class-I and MHC Class-II antigens were significantly up-regulated, while immunoglobulin, Cis-aconitic acid decarboxylase, and G0/G1 switch protein-2 were down-regulated in the old group. Some of these DEGs were detected only in the old group. Highly expressed genes included sperm-associated antigen 7, vasodilator-stimulated phosphoprotein isoform X5, and hCG2040310, which may be associated with the aging in golden snub-nosed monkey. While, some genes were only expressed in the young group, including major centromere autoantigen B, and proteinase-activated receptor 1, which may be associated with the development in young monkeys.

Analysis of DEGs in adult and old groups of R. r. roxellana
A total of 58 DEGs were detected in the adult and old groups, and the difference between the group was statistically significant (Table S8). We found that 25 genes were up-regulated, while 23 genes were down-regulated in the old group compared to the adult group. GO enrichment analysis showed that DEGs mainly participated in immune responses, such as extension factor 1-alpha 1, MHC class-I antigen, transmembrane domain protein 1 isomer 1, high affinity immunoglobulin epsilon receptor alpha subunit, HLA I histocompatibility antigen Cw-3alpha, HLA II histocompatibility antigen DP beta isomer X1, HLA II histocompatibility antigen DP alpha, HLA II histocompatibility antigen DRB1 beta were significantly up-regulated, while immunoglobulin and protease-activated receptors were downregulated in the old group.

Analysis of DEGs in young and adult groups in R. r. qinlingensis
A total of 117 DEGs were identified in young and adult groups in R. r. qinlingensis (Table S9) Besides, immunoglobulin, and platelet factor 4 were down-regulated in the old group.

Discussion
In this study, blood samples were obtained from nine golden snub-nosed monkey and their transcriptomes were sequenced, to reflect their physiological status. We obtained 57.31GB highquality data and identified 584321 SNP loci. We also discovered 3,817 new genes, of which 3,004 were annotated by database. These data provided a valuable resource for improving genomic annotation and identification of biological pathways in the golden Snub-nosed monkey.
The major histocompatibility complex (MHC) plays a central role in adaptive immune system of the vertebrate [30]. The description of MHC gene sequences and expression in whole blood will be clearly useful in ongoing efforts to understand the role of this important multigene family of protein in ecology and evolution. The expression level of MHC genes was analyzed, and it was found that the MHC gene was highly diverse in golden snub-nosed monkey. The expression of Aime-128 MHC class-I was the highest in MHC genes, which was consistent with the studies in giant panda [17]. This suggested that the Aime-128 may play a crucial role in immune system of golden snub-nosed monkey.
These findings were also useful in annotation of the subpopulations. In wild, R. r. roxellana inhabits the western part of the Sichuan province and the southern part of the Gansu province, while R. r. qinlingensi inhabits only the southern Shaanxi province, which is a high altitude area of the southern Qinling Mountains [5]. As the monkeys used in this study all lived in the zoos, the environmental conditions were comparatively similar. Thus, the differences seen in gene expression of these monkeys were less aligned with characteristics of their native environment. There were few differentially expressed genes among subspecies, which could not be well-enriched in KEGG pathway.
However, among them, 40S ribosomal protein SA (up-regulated in R. r. qinlingensi) was involved in cell adhesion, migration and growth. A precursor of surface receptor of LAMR was widely distributed in tumor cells and some normal cells attached to the basement membrane [31], which means that R. r. qinlingensi subspecies were facing more severe physiological state in nature. Moreover, the expression of gene of dynactin subunit increased significantly in R. r. qinlingensi subspecies, which is thought to participate in the immune process, indicating that R. r. qinlingensi subspecies may adapt to worse natural environment in the wild.
At the same time, we found that gender had minimal effect on the gene expression patterns among R. r. roxellana, which was consistent with the recent studies on giant panda [17]. However, we identified DEGs that were significantly different between the males and females in R. r. qinlingensi.
Compared to females, many genes having essential role in male reproduction were highly expressed in males. For example, ATP-dependent RNA helicase DDX3 was encoded by the DDX3Y gene, which may play a role in spermatogenesis [32]. Lysine (K)-specific demethylase 5D (KDM5D) was located on the AZFb region of the Y chromosome and encoded a JmjC-domain-containing protein [33]. These male-specific genes may be crucial for male reproductive system.
Previous studies have shown that age has a broad impact on gene expression levels [17,34].
Therefore, we believe that it was a reasonable choice to use the blood transcriptome to explore agerelated differences in these monkeys. As protected endangered animals, the research on golden snub-nosed monkey had some limitations. The samples were obtained only from the dead snub-nosed monkeys in the most previous studies. However, studies done using dead ones could not best capture the complexities and intricacies of living. There are a large number of immune cells and factors in the blood. Transcriptions reveal the dynamic changes in the gene expression of these blood cells at the genome level, including comprehensive the effect of genetics and environmental factors. However, being an endangered animal, collection of the blood samples from golden snub-nosed monkey was often difficult. We could only collect the blood samples during routine physical examination of these monkeys in Chengdu Zoo. Thus, no duplicate samples were collected from the old group. The statistics showed the average life expectancy of snub-nosed monkeys was generally 16-18 years, while 20+ years old were considered to live a long life. Typically, golden snub-nosed monkeys are considered to be sexually mature at the age of five years, therefore in this study, the monkeys aged 0-5.0 years were considered "young", those between 5.1 and 13 years were considered "adult" and those beyond 14 years were considered "old". For R. r. roxellana, the difference among 4 years old, 5.72 years old, and 6.6 years old may be more subtle. We speculated that there might be weak power to meaningfully detect changes, especially after correcting sex and subspecies. So, in the study, We deleted the content of DEGs between the young and the adult groups R. r. roxellana. We only compared the DEGs between the young and the old groups, and the adult and the old groups in R. r.

roxellana.
Functional enrichment analysis indicated that up-regulated genes with age were mainly related to innate immune response and T-cell activity. For example, elongation factor 1-alpha 1, a protein that promotes protein synthesis during cell cycle and elongation [35]. Whereas, apoptotic cells can secrete elongation factor 1-alpha 1 and other elongation factors to induce autoimmune response during cancer and help kill cancer cells [36]. MHC class-I gene and MHC class-II gene play a pivotal role by presenting peptide antigens and leading to CD4 + T cell-dependent cellular and humoral immune responses [37,38]. The up-regulated expression of these genes implied that the enhancement of immunity in adult group compared to the young group.
We also found that the down-regulated immune-related genes with age were mainly involved in antimicrobial response of B-cell activity. The expression of various immunoglobulin genes and B-cell antigen receptor decreased with age. Immunoglobulin helps the immune system recognize "self" and "non-self" substances, and binds to foreign substances such as bacteria and viruses, found in the blood of mammals and on the surface of B-cells The B-cell antigen receptor is a heterogenous oligomer complex composed of membrane-bound immunoglobulin and Ig-alpha/Ig-beta [40].
Membrane immunoglobulin recognizes the antigen and Ig-alpha/Ig-beta transduces antigen stimulus, resulting in cell proliferation and differentiation [41]. B-cell antigen receptor and immunoglobulin were closely related to acquired immunity. The down-regulation of B-cell antigen receptor and immunoglobulin may indicate the decline of acquired immune function in adult R. r. roxellana.
We also found that the expression of other genes decreased with age, e.g. platelet factor 4, a small cytokine of the CXC chemokine family. Platelet factor 4 can promote blood coagulation by regulating the role of heparin-like molecules, presumably contributing to wound repair and inflammation [42].
The G0/G1 switch protein 2 (G0S2) kills the cancer cell by promoting apoptosis [43]. However, the expression of G0/G1 switch protein 2 decreased with age, indicating that the immune function declines with age organisms.

Conclusions
This study investigated the blood transcriptome profile of the golden snub-nosed monkey using Illumina RNA-seq technology. Functional enrichment analysis indicated that up-regulated genes were mainly related to innate immune response and T-cell activity, while the down-regulated genes were mainly involved in B-cell activity. Our findings may contribute to understanding the molecular mechanisms of immune changes associated with ageing, which will provide a foundation for future study in golden snub-nosed monkey.

Ethics Statement
Blood samples were collected during routine physical examination by an experienced veterinarian of Chengdu Zoo, who has been taking care of golden snub-nosed monkeys for years. All samples collection and utility protocols were approved by the Chengdu Institute of Biology Animal Use Ethics.
Our experimental procedures complied with the current laws on animal welfare and research in China.

Sample Collection
Blood samples were collected during routine examinations from nine captive golden snub-nosed monkeys (five male and four females) of different age groups in Chengdu zoo, China (Table 5).
Samples were preserved on ice and promptly treated in the laboratory.

RNA isolation
Total RNA was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer's protocol. Residual genomic DNA was extracted from total RNA using DNase I (RNasefree). The RNA samples were sent to Biomarker Technologies Corporation (Beijing, China) for cDNA library construction.

Preparation of cDNA library and transcriptome sequencing
For library construction, 1 μg RNA of each sample was used. Total RNA was purified using oligo (dT) magnetic beads. The first strand of DNA was synthesized by random primers and reverse transcriptase using RNA as template. Next, dNTPs, RNase H and DNA polymerase I were used to synthesize the second strand of DNA.
The double-stranded DNA was purified, and the DNA with single-stranded ends was repaired by 3 to 5' terminal nucleic acid exonuclease and polymerase. A base tail was added to the 3' end of the repaired DNA. The library fragment was purified by kit in order to select the suitable length of the DNA fragment. The cDNA fragments were enriched by PCR to construct the final cDNA library [17]. The cDNA library was sequenced on the Illumina 2000 sequencing platform (Novogene Bioinformatics Technology Co., Ltd.) using the paired-end technology [44], by Beijing Genomics Institute (BGI), China.

Transcriptome assembly
The original sequence contained some low-quality sequences with sequencing joints, which could interfere with subsequent data analysis. To ensure the accuracy of the results; it was necessary to filter the original data and get clean reads [17].
High-quality clean reads were obtained by deleting low-quality data from the original data. These clean reads were mapped to the reference genome sequence of snub-nosed monkeys. According to the reference genome, only the complete matches were further analyzed and annotated. The reference genome was mapped using HISAT2 software 41].

Functional annotation
The assembled transcripts were searched against the Nr (NCBI non-redundant protein sequences database) KEGG (Kyoto Encyclopedia of Genes and Genomes), GO (Gene Ontology), COG (Cluster of Orthologous Groups) and Pfam (Protein family) by BLASTX to obtain protein functional annotation based upon sequence similarity [45,46].

Single nucleotide polymorphism (SNP) Analysis
Each sample was classified using the software picard-tools V1.41 and samtools v0.1.18, then duplicate data were deleted. The SNP sites were analyzed using the software Genome Analysis Toolkit (GATK). The original variant cell format (VCF) was filtered by GATK standard filtering method and other parameters, and only the SNP with distance >5 were retained.

Gene expression quantification
For each transcript, fragments per kilobase of transcript per million mapped reads (FPKM) were calculated to obtain the true level of transcript or gene expression in different samples.

Differentially expressed genes (DEG)
For biologically repetitive samples, DEseq2 [47] was used to analyze the differentially expressed genes between the two groups. Owing to the problem of false-positive results, Benjamini-Hochberg method was used to correct the significant p value of the original hypothesis test. The error detection rate was used as one of the criteria for screening differentially expressed genes. The criterion for screening differentially expressed genes by DEseq2 was that the false discovery rate (FDR) value must be < 0.05 and log2 (fold change) > 1. For the samples without biological replication, EBseq was used to analyze the differentially expressed genes between the two samples.
Declarations request sent to the corresponding author.

Ethics approval and consent to participate
Ethics approval: The collection of blood samples was performed in strict accordance with the guidelines of Chengdu Institute of Biology Animal Use Ethics.

Consent for publication
Not applicable. Tables   Table 1 statistics