Using Blood Transcriptome Analysis to Determine the Changes in Immunity and Metabolism of Giant Pandas with Age

Simple Summary Giant pandas are considered a national treasure in China. Understanding the changes in immunity and energy metabolism of giant pandas with age might help develop more scientific guidelines for managing the health of giant pandas. Here, 48 giant pandas were collected, and their transcriptome was analyzed. The results showed that the immune system and energy metabolism of giant pandas changed significantly with age. Abstract A low reproductive rate coupled with human activities has endangered the giant panda, a species endemic to southwest China. Although giant pandas feed almost exclusively on bamboo, they retain carnivorous traits and suffer from carnivorous diseases. Additionally, their immune system is susceptible to aging, resulting in a reduced ability to respond to diseases. This study aimed to determine the genes and pathways expressed differentially with age in blood tissues. The differentially expressed genes in different age groups of giant pandas were identified by RNA-seq. The elderly giant pandas had many differentially expressed genes compared with the young group (3 years old), including 548 upregulated genes and 401 downregulated genes. Further, functional enrichment revealed that innate immune upregulation and adaptive immune downregulation were observed in the elderly giant pandas compared with the young giant pandas. Meanwhile, the immune genes in the elderly giant pandas changed considerably, including genes involved in innate immunity and adaptive immunity such as PLSCR1, CLEC7A, CCL5, CCR9, and EPAS1. Time series analysis found that giant pandas store glycogen by prioritizing fat metabolism at age 11, verifying changes in the immune system. The results reported in this study will provide a foundation for further research on disease prevention and the energy metabolism of giant pandas.


Introduction
The giant panda (Ailuropoda melanoleuca) is endemic to China. However, they are in danger of extinction due to habitat shrinkage as a result of the expansion of the human population [1,2]. The giant panda's diet consists almost entirely of bamboo, but it retains the digestive tract characteristics of a carnivore, making it difficult for the giant panda to obtain energy from food. Recent studies have shown that giant pandas can facilitate bamboo digestion through a unique gut microbiome [3]. A mutation in the DUOX2 gene allows pandas to maintain low thyroid hormone levels to reduce energy expenditure. The

1.
mRNA Isolation mRNA molecules were purified from total RNA using oligo(dT)-attached magnetic beads.
2. mRNA Fragment mRNA molecules were fragmented into small pieces using fragmentation reagent after reaction a certain period in proper temperature.

3.
First Strand cDNA Synthesis An appropriate amount of primers was added to the interrupted sample, mixed well, and reacted at a suitable temperature on a thermal cycler for a certain period of time to open the secondary structure and combine with the primers. Then, the first-strand synthesis reaction system was added, and the mix was prepared in advance. The first-strand cDNA was synthesized according to the corresponding procedure on the instrument.

Second Strand cDNA Synthesis
A second-strand synthesis reaction system was prepared (using dUTP instead of dTTP), then reacted on Thermomixer at an appropriate temperature for a certain amount of time to synthesize second-strand cDNA. The reaction product was purified by magnetic beads.

5.
End repair and Add 'A' Prepare the end repair & add "A" reaction system, react on a thermal cycler for a certain time, under the action of enzymes, repair the sticky ends of the cDNA doublestranded by reverse transcription, and add A base to the 3 end. 6.
Adaptor Ligation Prepare the linker connection reaction system. React on a thermal cycler at a suitable temperature for a certain time. Under the action of the enzyme, connect the linker to the A base. The reaction product was purified by magnetic beads.

PCR
Prepare the PCR reaction system, digest the U-labeled second-strand template with UDG enzyme and perform PCR amplification. PCR products were purified with XP Beads, and dissolved in EB solution. 8.
Library Quality Control The library was validated on the Agilent Technologies 2100 bioanalyzer (Agilent, Santa Clara, CA, USA).

Circularization
The double stranded PCR products were heat denatured and circularized by the splint oligo sequence. The single strand circle DNA (ssCir DNA) were formatted as the final library.

Sequencing
The library was amplified with phi29 to make DNA nanoball (DNB) which had more than 300 copies of one molecular. The DNBs were load into the patterned nanoarray and single end 50 (pair end 100/150) bases reads were generated in the way of combinatorial Probe-Anchor Synthesis(cPAS).

Quality Control
This project used the filtering software SOAPnuke (version: V1.4.0, The Beijing Genomics Institute, Guangdong, China), independently developed by UBC, for filtering. The specific steps were the following: (1) Remove reads containing joints (joint contamination); (2) Remove reads with unknown base N content greater than 5%; (3) Removal of low-quality reads (low-quality reads were defined as those with a mass value of fewer than 15 bases that accounted for more than 20% of the total number of bases in the reads).
The filtered "Clean Reads" were saved in FASTQ format.

Gene Enrichment Analysis
According to GO annotation results and official classification, differential genes were functionally classified, and phyper function in R software was used for enrichment analysis. According to KEGG annotation results and official classification, the differential genes were classified into biological pathways, and phyper function in R software was used for enrichment analysis. The p-value was then FDR corrected. Usually, a Q value ≤ 0.05 was considered a significant enrichment.

Transcriptome Sequencing and Assembly
In this project, 48 samples were tested using the DNBSEQ platform. The age distribution of the samples ranged from 3 to 30 years. Table 1 shows the sample ages and the biological repetitions. Each sample yielded an average of 6.44 Gb of data. The average rates of genome alignment and gene set alignment were 91.05% and 82.71%, respectively. A total of 22,013 genes were detected, including 21,104 known genes and 909 predicted new genes. Out of the 30,362 new transcripts detected, 29,415 belonged to new alternative splicing subtypes of known protein-coding genes, and 947 belonged to new protein-coding genes.

Identification of Age-Related Differentially Expressed Genes
Giant pandas have a natural life span of about 30 years. Although pandas typically live for about 20 years in the wild, they can live up to 30 years of age in captivity [17]. With the three-year-old group as the control, the statistics of differentially expressed genes in different age groups are shown in Figure 1. The age 3 group, the age 11 group, the age 18-19 group, and the age 27-28-30 group were selected, with at least three biological replicates for each group to facilitate the analysis. Differential genes (DEG) were identified and analyzed using DEseq2 software [18]. The results showed that 151 genes were significantly different in the age 11 group compared to the age 3 group, 83 genes were upregulated, and 68 genes were downregulated. The expression of 111 genes was significantly different in the age 18-19 group, 43 of which were upregulated and 68 of which were upregulated. Further, compared with the age 3 group, 949 genes were significantly different in expression level in the age 27-28-30 group, of which 548 genes were upregulated, and 401 genes were downregulated ( Figure 2A).

Gene Ontology Enrichment Analysis of Differentially Expressed Genes
GO enrichment analysis was performed on the DEGs to understand their biological role further. Analyzing the GO term, it was observed that compared to the age 3 group, the upregulated genes in the age 11 group were significantly enriched in response to another organism (GO:0051707), defense response (GO:0006952), and response to the virus (GO:0009615) and 78 GO terms (Figure 3, Supplementary Table S1). However, the downregulated genes were not significantly enriched in any of the GO terms. Compared with the age 3 group, the upregulated genes in the age 18-19 group were significantly enriched

Gene Ontology Enrichment Analysis of Differentially Expressed Genes
GO enrichment analysis was performed on the DEGs to understand their biological role further. Analyzing the GO term, it was observed that compared to the age 3 group, the upregulated genes in the age 11 group were significantly enriched in response to another organism (GO:0051707), defense response (GO:0006952), and response to the virus (GO:0009615) and 78 GO terms (Figure 3, Supplementary Table S1). However, the downregulated genes were not significantly enriched in any of the GO terms. Compared with the age 3 group, the upregulated genes in the age 18-19 group were significantly enriched in only two GO terms, defense response to a virus (GO:0051607) and response to a virus (GO:0009615). The downregulated genes were enriched in only two GO terms, immunoglobulin complex (GO:0019814) and negative regulation of Lipoprotein lipase activity (GO:0051005). GO enrichment analysis of differentially differentiated genes in the age 27-28-30 group found that Immune System Process (GO:0002376), Immune Response

Pathway Enrichment Analysis of Differentially Expressed Genes
KEGG pathway enrichment analysis was also performed to evaluate the biological significance of DEGs further. It was observed that upregulated genes were significantly enriched in cytosolic DNA sensing pathway (KO04623) in the age 11 group. The downregulated genes were not significantly enriched in any pathway. In the age 18-19 group, the upregulated genes were not significantly enriched in any pathway, but the downregulated genes were significantly enriched in the B Cell receptor signaling pathway (KO04662) and fat digestion and absorption pathway (KO04975). For the age 27-28-30 group, upregulated genes were significantly enriched in 15 pathways ( Figure 2D, Supplementary Table S4), while the downregulated genes were significantly enriched in 11 pathways ( Figure 2F, Supplementary Table S5).

Time Series Analysis
According to the dynamic gene expression changes over time, four giant panda models were identified, and the trend of change of partial DEG expression with age was summarized through soft clustering ( Figure 4A). The 630 genes in cluster 1 were highly expressed in 11-year-old giant pandas but were low in other age groups. The 1182 genes in cluster 2 were highly expressed in the young giant pandas, and the expression levels of these genes gradually decreased with an increase in age. The expression of the cluster 3

Pathway Enrichment Analysis of Differentially Expressed Genes
KEGG pathway enrichment analysis was also performed to evaluate the biological significance of DEGs further. It was observed that upregulated genes were significantly enriched in cytosolic DNA sensing pathway (KO04623) in the age 11 group. The downregulated genes were not significantly enriched in any pathway. In the age 18-19 group, the upregulated genes were not significantly enriched in any pathway, but the downregulated genes were significantly enriched in the B Cell receptor signaling pathway (KO04662) and fat digestion and absorption pathway (KO04975). For the age 27-28-30 group, upregulated genes were significantly enriched in 15 pathways ( Figure 2D, Supplementary Table S4), while the downregulated genes were significantly enriched in 11 pathways ( Figure 2F, Supplementary Table S5).

Time Series Analysis
According to the dynamic gene expression changes over time, four giant panda models were identified, and the trend of change of partial DEG expression with age was summarized through soft clustering ( Figure 4A). The 630 genes in cluster 1 were highly expressed in 11-year-old giant pandas but were low in other age groups. The 1182 genes in cluster 2 were highly expressed in the young giant pandas, and the expression levels of these genes gradually decreased with an increase in age. The expression of the cluster showed that these genes were also involved in antigen processing and presentation, thermogenesis, and intestinal IgA production. Cluster 3 expression pattern genes were involved in regulating cell cycle and gene expression, as well as B cell receptor and T cell receptor signaling pathways. Cluster 4 genes were involved in the immune (innate immune) response, inflammatory response regulation, and Toll-like receptor signaling. (A)

Discussion
In this study, many biological samples with a wide age distribution (3-30 years) and sufficient biological replicates were collected for scientific analysis. Differentially expressed genes based on age as a grouping basis were studied because previous studies have shown that age has a much more significant impact on gene expression patterns than gender [19]. Q < 0.05 was chosen as the threshold for all differences.  Table S7) enrichment analysis were performed to assess the function of DEG with different patterns during giant panda development. The analyses revealed that the genes with high expression in cluster 1 at the age of 11 were mostly related to energy metabolism and ATP generation. They were also involved in the differentiation of Th1 and Th2 cells and in heat production. Cluster 2 genes were involved in gene expression, peptide biosynthesis, translation, and synthesis of various organelles. KEGG analysis showed that these genes were also involved in antigen processing and presentation, thermogenesis, and intestinal IgA production. Cluster 3 expression pattern genes were involved in regulating cell cycle and gene expression, as well as B cell receptor and T cell receptor signaling pathways. Cluster 4 genes were involved in the immune (innate immune) response, inflammatory response regulation, and Toll-like receptor signaling.

Discussion
In this study, many biological samples with a wide age distribution (3-30 years) and sufficient biological replicates were collected for scientific analysis. Differentially expressed genes based on age as a grouping basis were studied because previous studies have shown that age has a much more significant impact on gene expression patterns than gender [19]. Q < 0.05 was chosen as the threshold for all differences.
Genetic analysis of giant pandas of different age groups demonstrated no significant difference between the age 11 group and the age 3 group. GO enrichment analysis revealed that some upregulated genes such as CCL5, RSAD2, ACOD1, and STAT2 were involved in immune responses. CCL5 is a chemokine that is expressed initially in viral infection [20,21]. They are known to activate several cells directly involved in antiviral responses, including NK cells, T CD4+ lymphocytes, monocytes, mast cells, and dendritic cells [22]. RSAD2 is one of the few ISGs with direct antiviral activity. It also plays a vital role in regulating innate immunity [23]. STAT2, unlike others in the STAT family, is only involved in type I and III interferons (IFN-I/III) signaling pathways and also plays a role in innate immunity. Most of the genes involved in the immune response play essential roles in innate immune responses, indicating that the innate immunity of giant pandas is in continuous development during the growth. Additionally, CPT1A was highly expressed in the relatively young group. CPT1A is the first rate-limiting enzyme of fatty acid oxidation, and its high expression levels are conducive to fatty acid metabolism [24].
The number of differentially expressed genes between the age 27-28-30 group and the age 3 group was significantly different. Therefore, a hierarchical clustering of DEGs across seven samples was carried out. All the downregulated genes clustered in one group and all the upregulated genes clustered in the other. Similarly, the samples were also clustered into two groups by age ( Figure 2B). DEGs enrichment analysis showed that these genes were mostly related to immune function. In particular, pandas in the age 27-28-30 group were observed to have upregulated expressions of genes responsible for innate immunity, including CCL5, TLR2, OAS2, and ISG15s, compared with the age 3 group (Table 2). On the other hand, the genes related to adaptive immunity were significantly downregulated, including those related to B cell and T cell differentiation (Table 3). An aspect of aging is the decline in the integrity of the immune system, leading to an increase in the body's susceptibility to diseases [25]. Elderly giant pandas may compensate for the decline in adaptive immune function by improving their innate immunity. Further, the genes related to hematopoietic function were downregulated in older pandas. Interestingly, KEGG enrichment analysis of upregulated genes revealed a pathway associated with osteoclast differentiation, which has previously been reported to be associated with bone immunity and other bone diseases [26]. The downregulated genes were involved in the production of intestinal IgA, indicating that certain intestinal immune functions could be lacking in elderly giant pandas. Giant pandas 6-20 years of age were considered adults. The gene expression analysis in the two selected age groups, i.e., age 11 and age 18-19, revealed many DEGs. Compared with age 11, 121 genes were upregulated, and 31 genes were downregulated in the age 18-19 group. GO enrichment analysis showed that these upregulated genes were significantly enriched in membrane-enclosed lumen (GO:0031974), organelle lumen (GO:0043233), protein binding (GO:0005515), and 94 other GO terms (Supplementary Table  S8). Downregulated genes were not significantly enriched in any of the GO terms. The most significantly upregulated gene was HSP90AA1. This gene encodes a protein, Hsp90α, which has been reported to be associated with various diseases [27]. Other than that, there were no significant immune or metabolic changes in the adults.
In the aged giant panda group (21-30 years old), a comparative study of the change in expression levels of the genes with aging is crucial for preventing age-related diseases. The data of the 21-22 age group and the 27-28-30 age group were compared. It was observed that, in the age 27-28-30 group, 58 genes were upregulated, and 54 genes were downregulated. Enrichment analysis of DEGs revealed that the upregulated genes were significantly enriched in nine GO terms (Supplementary Table S9). However, the downregulated genes were not significantly enriched in any of the GO terms. GO term analysis of the upregulated genes showed that the upregulated genes, such as PLSCR1, CLEC7A, JUNB, and BATF3, were primarily involved in immune responses. In human studies, PLSCR1 has been found to play an essential role in IFN-dependent antiviral responses [28]. PLSCR1 directly interacts with and affects the function of several viral proteins, playing a pivotal role in PLSCR1-mediated antiviral activity [29][30][31]. The CLEC7A gene encodes dectin-1 and is widely expressed in the immune system. It is involved in biological processes such as allergies, cancer, autoimmune diseases, aseptic inflammation, and even aging [32]. JUNB mediates T cell development in the immune system [33]. BATF3 promotes the development of CD8α+ Conventional Dendritic cells (cDCs) by maintaining IRF8 and participating in the immune response [34]. MHC class II genes regulate Toll-like receptors, mediating innate immune responses in addition to their central role in adaptive immunity [35]. GO classification of downregulated genes was used to explain their functions. The analysis found that genes involved in immune response, such as CCL5, CCR9, and EPAS1, were downregulated. CCR9 is a G-protein-coupled receptor involved in biological processes such as immunity, inflammation, and tumors [36,37]. This result indicates that the elderly giant pandas enhanced the expression of some immune genes with age and compensated for the downregulated immune genes to resist the invasion of viruses, including innate immunity and adaptive immunity.
Time series analysis revealed that giant pandas had a high energy metabolism rate in the age 11 group. The genes strongly expressed at the age of 11 were enriched into the energy-related GO term. Protein misfolding is more likely to occur in high-salt environments, manifested in biological aging. ATP is an energy currency in animals, maintaining various metabolic networks of cells and also playing a pivotal role in protein stability through various mechanisms [38]. Almost all of the genes highly expressed in 11-year-old pandas are directly or indirectly related to ATP production. The analysis of cluster 2 genes revealed that giant pandas have a relatively vigorous metabolism when they are young. However, with age, the gene expression levels of adaptive immunity of the giant panda decrease. This decreased expression level of essential genes involved in antigen processing and presentation is bound to affect the integrity of the adaptive immune system. The downregulation of some of the genes engaged in heat production indicates that elderly pandas may be capable of generating less heat. The analysis of cluster 4 genes further proved that the expression of genes related to innate immunity is upregulated with aging in giant pandas. The genes related to regulating inflammatory response are also upregulated to a certain extent.

Conclusions
Giant pandas are known to undergo significant changes in immune-related genes with age. This study showed that congenital immunity developed faster in giant pandas as they grew from juvenility to adulthood. This was consistent with the results of a previous study showing that the immune system of giant pandas gradually improved from infancy to adulthood [39]. Our study found that innate immunity might develop faster than adaptive immunity. There was also an upregulation in innate immunity and downregulation in adaptive immunity in old giant pandas relative to the young. Additionally, previous studies have not proven the humoral immunity and change in T cell immunity [40], these factors supplementing the prevention of aging-related diseases. The immune system of the elderly giant panda is still changing. In the process of aging, the body undergoes the adjustment of some immune genes to compensate for the loss of other immune functions. It was also shown that giant pandas at the age of 11 promoted fat consumption to meet the storage of liver glycogen and high ATP production. Further studies are required to understand the biological significance of these processes.

Supplementary Materials:
The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/vetsci9120667/s1, Table S1: GO enrichment analysis of upregulated DEGs in age groups 11, Table S2: GO enrichment analysis of upregulated DEGs, Table S3: GO enrichment analysis of down-regulated DEGs, Table S4: KEGG pathway enrichment analysis of upregulated DEGs, Table S5: KEGG pathway enrichment analysis of down-regulated DEGs, Table  S6: Gene Ontology, Table S7: Kyoto Encyclopedia of Genes and Genomes, Table S8: GO enrichment analysis of upregulated DEGs in age groups 18-19 compared to the age groups 11, Table S9: GO enrichment analysis of upregulated DEGs in age groups 27-28-30 compared to the age groups 21-22.