Transcriptome analysis reveals immune-related gene expression changes with age in giant panda (Ailuropoda melanoleuca) blood

The giant panda (Ailuropoda melanoleuca), an endangered species endemic to western China, has long been threatened with extinction that is exacerbated by highly contagious and fatal diseases. Aging is the most well-defined risk factor for diseases and is associated with a decline in immune function leading to increased susceptibility to infection and reduced response to vaccination. Therefore, this study aimed to determine which genes and pathways show differential expression with age in blood tissues. We obtained 210 differentially expressed genes by RNA-seq, including 146 up-regulated and 64 down-regulated genes in old pandas (18-21yrs) compared to young pandas (2-6yrs). We identified ISG15, STAT1, IRF7 and DDX58 as the hub genes in the protein-protein interaction network. All of these genes were up-regulated with age and played important roles in response to pathogen invasion. Functional enrichment analysis indicated that up-regulated genes were mainly involved in innate immune response, while the down-regulated genes were mainly related to B cell activation. These may suggest that the innate immunity is relatively well preserved to compensate for the decline in the adaptive immune function. In conclusion, our findings will provide a foundation for future studies on the molecular mechanisms underlying immune changes associated with ageing.


INTRODUCTION
The giant panda (Ailuropoda melanoleuca), an endangered species endemic to western China, has long been threatened with extinction due to human population expansion and habitat destruction [1,2].Recent research has shown that giant pandas are well adapted to a specialized bamboo diet via enigmatic gut microbiota and low energy expenditure [3][4][5].Contrary to previous studies, the latest genome sequencing and resequencing and genome-wide genetic studies have demonstrated that the giant panda has a relatively high genetic diversity [6][7][8].Over the years, extensive research has been conducted on genes of the major histocompatibility complex (MHC) in giant pandas and shown that DRB and DQA genes have low diversity but remain positive selection [9][10][11].However, most of these studies focused on the adaptations and genetic diversity of giant pandas and little is known about the gene expression of immune-related genes.AGING Giant pandas are vulnerable to a variety of carnivorespecific and mammalian infectious and parasitic diseases that seriously threaten the giant panda populations [12][13][14][15][16].For example, Baylisascaris schroederi is the most common parasite in wild and captive giant pandas and was responsible for half of panda deaths between 2001 and 2005 [15].Likewise, canine distemper virus (CDV) is highly contagious and fatal and affects a wide range of mammals, including the giant panda [17].In 2014-2015, an outbreak of CDV infection in giant pandas caused the death of five infected animals in two months [18].
The immune system is vulnerable to age-associated alterations, which accumulate to produce a progressive deterioration and lead to an increased incidence of infectious diseases [19].Although aging is an inevitable biological process and a powerful risk factor for many diseases, the underlying molecular mechanisms that lead to generalized disease susceptibility are largely unknown [20,21].Aging has proven difficult to dissect in part due to its interactions with environmental influences, other genetic factors, and a large number of age-related diseases [22].Immunosenescence, defined as immune changes with ageing, is an unavoidable life process and has been characterized in several species, such as humans [21,23], mice [24], zebra finches [25] and wolves [26].A comprehensive meta-analysis of age-related gene expression profiles indicated that signatures of aging most notably involve an overexpression of inflammation and immune response genes and an underexpression of genes associated with energy metabolism [27].However, the effects of aging on giant panda have not been characterized, especially the age-related changes of immune system.Moreover, understanding the principles of giant panda immune system is vital for the development vaccines that can elicit protective immunity [28].
The aims of this study were to determine which genes and pathways show differential expression with age in giant panda blood tissues and to understand age-related alterations of the immune system.We used RNA-seq technology to identify age-related differentially expressed genes (DEGs) in giant panda blood samples and performed functional enrichment analysis for these DEGs.The collective data generated in this study may represent a valuable resource to enable further advancements in immunological research in giant pandas.

Transcriptome sequencing and assembly
RNA prepared from blood tissues of four giant pandas were subjected to RNA sequencing using Illumina HiSeq™ 2000.Together with the three giant panda transcriptomes from our previous study, we totally acquired approximately 187 million raw paired-end reads and 162 million remained after removing adaptor sequences and discarding low quality reads.The total length of the reads was about 35.8 gigabases (Gb).We aligned each of the seven short-read libraries onto the giant panda Ensembl reference genome (ailMel1) separately and found that an average of 19.86 million high-quality reads (85.68%) could be successfully mapped to the genome per sample.This included a mean of 19.31 million uniquely mapped reads (97.23%) per sample, indicating that the majority of the paired reads aligned correctly (Table 1).These high-quality reads were assembled into 56,543 genes, giving rise to 88,071 transcripts containing 46,091 (52.33%) transcripts that had more than one exon (Supplementary Figure 1).The comparison of similarities between the assembled transcripts and the Ensembl transcripts showed that only 2943 (3.34%) matched exactly with annotated intron chain and a total of 33,600 (38.15%) were identified as potentially novel isoforms of a predicted Ensembl transcript with at least one splice junction shared.The majority of transcripts (34,660,39.35%) were annotated as intergenic transcripts and a small proportion of transcripts (7824, 8.88%) entirely fell into reference intron (Supplementary Table 1).

Identification of age-related differentially expressed genes
The samples were categorized by giant panda age into a young group (i.e. 2, 5, 6 years) and an old group (i.e. 18, 19, 21 years) according to an earlier study [29].The giant panda aged 12 years (i.e. a middle-aged adult) was not used to identify DEGs due to a lack of biological replicates.In total, 210 genes were found to show significant differential expression in giant panda blood samples between the young group and old group, after adjusting for confounding factors such as sex and library (Figure 1A, Supplementary Table 2).Among those DEGs, 146 genes were up-regulated while 64 genes were downregulated in the old group compared to young group, and the annotated up-and down-regulated genes were 113 and 49, respectively (Figure 1B).In addition, we conducted a hierarchical clustering of DEGs across six samples.Genes that were up-regulated clustered into one group and the genes that were down-regulated clustered into another group.Similarly, young samples were clustered into one group and old samples were clustered into another group (Figure 1C), which also highlights the expression level of DEGs in the old group compared to the young group.To determine the relationship of these DEGs with age, we compared all the detected DEGs to age-related genes in purified human immune cells, human peripheral blood and wolf leukocytes.We found that 23 DEGs were presented in human studies [21] and only one DEG was presented in wolf study [26].The 20 most significantly up-and down-regulated genes are shown in Tables 2 and 3. We also found that 2 MHC class I genes (LOC105855951, LOC102746195) and 3 MHC class II genes (LOC100464025, LOC101370503, LOC100464780) were both up-regulated in the old group.

Gene ontology enrichment analysis of differentially expressed genes
To gain insights into the biological roles of the DEGs, we performed a GO categories enrichment analysis.Only the annotated DEGs were selected and tested against the background set of all genes with GO annotation.We examined the GO categories separately.We found that for up-regulated genes with age, almost all the significantly enriched GO terms in biological process category were associated with the immune system, such as, immune response (GO:0006955), negative regulation of viral genome replication (GO:0045071), and defense response to virus (GO:0051607).MHC class II protein complex (GO:0042613) and GTP binding (GO:0005525) were the most significantly enriched GO term in the cellular component and molecular function categories, respectively (Figure 2, Supplementary Table 3).For down-regulated genes with age, the most significantly enriched biological process GO terms were complement activation, classical pathway (GO:0006958), B cell receptor signaling pathway (GO:0050853) and Fcepsilon receptor signaling pathway (GO:0038095) (Figure 3, Supplementary Table 4).Interestingly, several up-regulated genes were enriched in retina homeostasis (GO:0001895) and regulation of blood pressure (GO:0008217), which might be related to age.

Pathway enrichment analysis of differentially expressed genes
To further evaluate the biological significance of the DEGs, we also performed KEGG pathway enrichment analysis.Hypergeometric tests with a P value cutoff of 0.05 was used as the criteria for pathway detection.We found that up-regulated genes were significantly enriched in 35 pathways, among which 22 pathways were disease related and 7 pathways were immune related.However, down-regulated genes were not significantly enriched to any pathways and this may be because many downregulated genes were not annotated in KEGG pathways.Figure 4 illustrates the result of pathway enrichment of up-regulated genes, in which antigen processing and presentation (ko04612), cytosolic DNA-sensing pathway (ko04623) and RIG-I-like receptor signaling pathway (ko04622) were most significantly enriched.

Protein-protein interaction network of differentially expressed genes
We have performed a protein-protein interaction analysis of all DEGs.A total of 122 interactions between 36 DEGs were extracted from the STRING database.Interferonstimulated gene 15 (ISG15), signal transducer and active tor of transcription 1 (STAT1), interferon regulatory factor 7 (IRF7) and DExD/H-Box Helicase 58 (DDX58) that play important roles in response to pathogen invasion were at key position of the interaction network (Figure 5).

DISCUSSION
In a previous study, we have characterized giant panda blood transcriptome and identified 15 immune pathways where more than 70% of the total known genes were mapped by assembled transcripts [30].In the present study, we still applied the widely used RNA-seq approach to highlight significant DEGs in blood between young and old giant pandas.The previous work has proved that age had a broad impact on gene expression levels, whereas sex had very minimal effects on gene expression patterns [26].Therefore, it is appropriate to investigate age-related changes in giant panda by sequencing the blood transcriptome.To the best of our knowledge, this is the first study to assess the impact of age on immune-related gene expression.
The assembly results unfortunately underlined the incompleteness of the giant panda genome annotation of current version, and the need for further studies to decipher its complexity.We have provided a valuable resource for improvement of genome assembly and annotation.For example, a recent study has improved the genome assembly and identified novel transcripts by using transcriptomes from 12 giant panda tissues, excluding blood [31].Another major outcome of transcrip-  tome sequencing is the ability to identify DEGs between different conditions [32].We identified 210 DEGs, which is much lower than the 625 and 1497 DEGs found in wolves [26] and humans [21], respectively.This may suggest that age-related alterations of gene expression are species-specific [33].The most significantly up-regulated gene SERPING1, as a signature of mammalian aging, encodes the complement component 1 (C1) inhibitor [27].The C1 inhibitor's main function is the inhibition of several different proteases in the complement, contact, coagulation and fibrinolytic systems and has been reported to be related to an increased risk of several age-related human diseases [34][35][36].The most significantly down-regulated gene C4BPA encodes the alpha-subunit of complement 4 binding protein.C4BPA is a key regulator of the complement system and has been considered a major factor in the agerelated decline of the immune response due to decreased expression [37].
We identified age-related DEGs that may be critical in the progression of the giant panda aging process, such as ISG15, STAT1, IRF7.ISG15 is an innate immune and interferon-induced protein that plays a central role in response to viral infection [38].ISG15 is also a stress-response gene that has been implicated as a tumor suppressor and contributor to inflammatory responses [39].ISG15 expression was regulated by telomere length and increased with shortening telomeres caused by aging [40].STAT1 carries tumor suppressor functions and regulates various cellular activities, such as apoptosis, angiogenesis, invasion and evasion of the immune system [41,42].Previous studies indicated that loss of STAT1 expression has been implicated in the pathobiology of various types of human cancer [41].IRF7 has been discovered as the crucial regulator of type I interferons against pathogenic infections and plays important roles in immune modulation and regulation of oncogenesis and apoptosis [43].Our results showed ISG15 can interact Furthermore, the invasive viruses can lead to activation of IRFs and induction of interferons (IFNs) resulting in antiviral IFNs signal through the JAK/STAT pathway to induce ISG production [44].
Aging is associated with functional integrity of the immune system declines leading to increased susceptibility to aging diseases [45].Both GO and KEGG enrichment analysis showed that up-regulated genes were mainly involved in the immune response, particularly defense against viruses.The involved genes, such as IFIT5, BST2 and OAS1, were closely related to the innate immunity response [46][47][48].In contrast, the down-regulated genes were implicated in the B cell receptor signaling pathway and positive regulation of B cell activation.In accordance with observations in wolf blood transcriptome with age, down-regulated genes were mainly associated with B cells and up-regulated genes were significantly associated with innate immunity [26].
Although, both innate and adaptive immune responses are affected by the aging process, the adaptive response seems to be more affected by age-related changes [49].However, contrary to the previous results observed in the wolf [26], MHC class I and II genes were up-regulated with age in our study.In addition to their central role in adaptive immunity, both MHC class I and II genes can regulate Toll-Like Receptor to mediate the innate immune response [50,51].Taken together, these suggest the possibility that the aging immune system maintains a persistently activated innate immune response to compensate for the decline in the adaptive immune function [52].Figure 3. GO enrichment analysis of down-regulated DEGs.
Vaccines appear as one of the most efficient medical interventions against infectious diseases [53] and have been also used to assist giant pandas to defend against CDV infections [54,55].However, the vaccines for giant pandas do not elicit consistent antibody titers [56].One explanation may be that the attenuated vaccine is inadequate to stimulate giant pandas to produce a high level of antibodies against CDV [55].Another possible explanation might be aging-associated declines in immune system function.It is clear that the decline in adaptive immunity leads to dramatically reduced vaccine responses and vaccine longevity in older individuals [57].We have noted that the immune response genes, such as IGHV3-23 and IGKC, that participate in the antigen recognition and B cell activation were downregulated with ageing.This may lead to a decline of adaptive immunity resulting in decreased response to vaccines in giant pandas.In addition, the decreased ability of the adaptive immune system is closely related with the decline in production of naive lymphocytes and the accumulation of functionally impaired memory lymphocytes during aging [58].
In conclusion, we identified 210 genes which differentially expressed between old and young giant pandas.The upregulated genes with age were mainly involved in the innate immune response, particularly in defense against virus.Instead, the down-regulated genes with age were mainly related to B cell activation.This suggests that innate immunity is relatively well preserved in aged giant pandas, although adaptive immunity deteriorates due to a decline of naive lymphocytes and the expansion of incompetent memory lymphocytes.Moreover, ISG15, STAT1, IRF7 and DDX58 were identified as the hub genes that play a pivotal role in response to viral infection and their expression were up-regulated with ageing.This work provides insight into gene expression changes associated with aging and will provide a foundation for future studies on the molecular mechanisms underlying immune changes associated with ageing.

Sample collection
Blood samples were collected during a routine examination from four captive giant pandas at the Chengdu Research Base of Giant Panda Breeding in Chengdu, China.The other three giant panda blood transcriptomes were obtained from our previous study in NCBI SRA (project Accession no.SRP041998).Samples collection and utility protocols were approved by the Chengdu Institute of Biology Animal Use Ethics, which is responsible for the Chengdu Research Base of Giant Panda Breeding.The fresh blood samples were immediately stored in RNAlater (Ambion Inc., Austin, TX, USA).Total RNA was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer's protocol and treated with RNase-free DNase I.

Library preparation and sequencing
A total amount of 3 μg RNA per sample was used as input material for the RNA sample preparations.Sequencing libraries were generated using NEBNext ® Ultra™ RNA Library Prep Kit for Illumina ® (NEB, USA) following manufacturer's recommendations and index codes were added to attribute sequences to each sample.Briefly, mRNA was purified from total RNA using poly-T oligo- First strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase (RNase H − ).Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities.After adenylation of 3' ends of DNA NEBNext Adaptor with hairpin loop structure were ligated to prepare for hybridization.
In order to select cDNA fragments of 150~200 bp in length, the library fragments were purified with AMPure XP system (Beckman Coulter, Beverly, USA).Then 3 μl USER Enzyme (NEB, USA) was used with size-selected, adaptor-ligated cDNA at 37° C for 15 min followed by 5 min at 95° C before PCR.Then PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers and Index (X) Primer.Lastly, PCR products were purified (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system.The library preparations were sequenced on an Illumina Hiseq 2000 platform and 100 bp paired-end reads were generated.

Quality control and transcriptome assembly
The giant panda genome assembly ailMel1 and reference annotation were downloaded from Ensembl v83 (www.ensembl.org).In order to obtain high quality reads, we employed NGS QC Toolkit [59] with stringent criteria (high-quality paired reads with more than 90% of bases with Q-value ≥ 20 were retained) to remove the lowquality paired-end reads or reads containing adaptors.High-quality paired-end reads from seven libraries were mapped on giant panda reference genome using HISAT [60], separately.Each of the alignment output files were assembled into separate transcriptomes using StringTie [61], which produces a transcript GTF file.The transcript GTF files produced after assembling were merged together to generate a complete transcript annotation file using StringTie with-merge option.

Functional annotation
The merged assemblies were compared using the Cufflinks [62] inclusive utility Cuffcompare with the Ensembl annotations to find known and novel genes and isoforms, as well as transcripts expressed from intergenic regions.Only the genes with "=" and "j" class code were considered to be known genes.The GO annotations of these genes were extracted from Ensembl BioMart (http:// www.ensembl.org/biomart/martview).The unknown gene sequences were searched against Swiss-Prot and NCBI non-redundant (NR) database using BLASTX [63] with a typical cutoff E-value of 1E-5.Gene ontology (GO) was applied with the Blast2GO [64] to obtain annotation of unknown genes.

Gene expression estimation and differentially expressed genes analysis
We used the merged GTF file as the reference annotation file to guide the assembly process.StringTie was performed to estimate the genes expression level with -G, -b and -e options.The expression value of TPM (transcripts per million) and raw read counts for each gene were extracted from above result.Differential expression analysis was performed using DESeq2 [65] that takes the raw read count as the input.The results of all statistical tests were corrected for multiple testing with the Benjamini-Hochberg false discovery rate (FDR ≤ 0.05) and an absolute value of log 2 fold change ≥ 1 was used to determine significant differences in gene expression.

Gene enrichment analysis
To functional classify the DEGs, we performed GO enrichment analysis using clusterProfiler [66] and searched for over-represented GO terms in three categories, namely biological process, molecular function and cellular component.We also performed KEGG pathway enrichment analysis using KOBAS [67].An FDR cut-off of 0.05 was considered to be significantly enriched.

Protein-protein interaction network analysis
Protein-protein interaction networks of DEGs were constructed on the basis of the STRING [68] protein interaction database.Cytoscape [69] was used to visualize the protein-protein interaction network for DEGs, and find the hub genes.

Figure 1 .
Figure 1.Differentially expressed genes in old group compared to young group.(A) MA plot showing the distribution of gene expression plotted against log2 fold change for each gene.Red dots indicate differentially expressed genes (FDR ≤ 0.05), black dots indicate non-differentially expressed genes.(B) Number of up-and down-regulated DEGs and annotated DEGs.(C) Heat map plot of DEGs using TPM expression value of genes by adopting hierarchical clustering method.The expression values of six individual are presented after being centered and scaled in the row direction.Each column represents a specimen and each row represents a gene.Red color indicates genes that were up-regulated and blue color indicates genes that were down-regulated.

Figure 5 .
Figure 5. Protein-protein interaction network of differentially expressed genes.Size of the node is proportional to the number of DEGs interacted with it, and color of node represents Log2FoldChange in expression levels of DEGs between old and young giant pandas.