Research Paper Volume 12, Issue 21 pp 21777—21797

The immune and metabolic changes with age in giant panda blood by combined transcriptome and DNA methylation analysis

Xiaoyu Huang1,2,3, *, , Qingyuan Ouyang1, *, , Mingxia Ran1, *, , Bo Zeng1, , Linhua Deng2,3, , Shenqiang Hu1, , Mingyao Yang1, , Guo Li2,3, , Tao Deng2,3, , Ming He2,3, , Ti Li2,3, , Haidi Yang2,3, , Guiquan Zhang2,3, , Heming Zhang2,3, , Changjun Zeng1, , Jiwen Wang1, ,

  • 1 Farm Animal Genetic Resources Exploration and Innovation Key Laboratory of Sichuan Province, Sichuan Agricultural University, Chengdu 611130, Sichuan, China
  • 2 China Conservation and Research Center for the Giant Panda, Dujiangyan 611830, Sichuan, China
  • 3 Key Laboratory of State Forestry and Grassland Administration on Conservation Biology of Rare Animals in The Giant Panda National Park, Dujiangyan 611830, Sichuan, China
* Equal contribution

Received: June 9, 2020       Accepted: August 14, 2020       Published: November 7, 2020
How to Cite

Copyright: © 2020 Huang et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.


Giant panda (Ailuropoda melanoleuca) is an endangered mammalian species. Exploring immune and metabolic changes that occur in giant pandas with age is important for their protection. In this study, we systematically investigated the physiological and biochemical indicators in blood, as well as the transcriptome, and methylation profiles of young, adult, and old giant pandas. The white blood cell (WBC), neutrophil (NEU) counts and hemoglobin (HGB) concentrations increased significantly with age (young to adult), and some indicators related to blood glucose and lipids also changed significantly with age. In the transcriptome analysis, differentially expressed genes (DEGs) were found in comparisons of the young and adult (257), adult and old (20), young and old (744) groups. Separation of the DEGs into eight profiles according to the expression trend using short time-series expression miner (STEM) software revealed that most DEGs were downregulated with age. Functional analysis showed that most DEGs were associated with disease and that these DEGs were also associated with the immune system and metabolism. Furthermore, gene methylation in giant pandas decreased globally with age, and the expression of CCNE1, CD79A, IL1R1, and TCF7 showed a highly negative correlation with their degree of methylation. These results indicate that the giant panda’s immune function improves gradually with age (young to adult), and that changes in the methylation profile are involved in the effects of age on immune and metabolic functions. These results have important implications for the understanding and conservation of giant pandas.


The giant panda (Ailuropoda melanoleuca) is a first-class national protected species in China. Over the past 20 years, extensive research on the protection of giant pandas has been conducted all over the world. Previous studies showed very little difference in the physiological and biochemical indexes of giant panda blood between males and females, while several indicators, including the number of immune cells, were found to be linearly related to age [1]. The immune system undergoes profound transformations with age, and these changes are globally known as ‘immunosenescence’ in mammals [2]. Immune cells and circulating factors, including chemokines, cytokines, and other soluble molecules, are generally considered to be critical factors that affect the immune system during the aging process [3]. Exploring the changes in the giant panda’s immune cells with age and the underlying mechanisms is essential to developing strategies that delay or prevent the age-related decline in immunity. Aging is arguably the universal contributor to the etiologies of metabolic decay and related diseases, including type 2 diabetes mellitus, cardiovascular disease, and insulin resistance [4, 5]. Substantial evidence suggests that aging organisms exhibit homeostasis disorders of carbohydrate, amino acid, and fatty acid metabolism [6, 7]. Enhancing mitochondrial function or inhibiting glycolysis can extend lifespan and promote healthy aging in many species [8]. However, research on age-related changes in immune and metabolic functions of giant pandas is still lacking.

Two recent transcriptomic studies [9, 10] of the changes in gene expression with age in young and old giant pandas showed that the differentially expressed genes (DEGs) identified in the blood were associated with immunity. Nevertheless, aging involves two processes (young to adult and adult to old), and separate studies of these two processes are still lacking. Furthermore, changes in gene expression may be caused by DNA methylation [11], a process that has been widely reported to be involved in the aging process in humans and mouse [12, 13]. Approximately 0.5% of the genome in pandas is methylated [14] and DNA methylation has been shown to be involved in the regulation of cataract disease in old giant pandas [15]. However, the molecular mechanism of age-related immune and metabolic changes in giant pandas has not yet been systematically investigated using multi-omics technology.

In this study, we systematically analyzed the physiological and biochemical changes in the blood of giant pandas with age and investigated the underlying molecular mechanisms. We identified age-related DEGs and used the short time-series expression miner (STEM) to classify the age-related DEGs (upregulated and downregulated). Each type of DEG was used for functional enrichment analysis, and the expression and degree of methylation of the functional DEGs were used for correlation analysis. Furthermore, we constructed the first age-related dynamic methylation profile of panda blood. Identification and functional analysis of differentially methylated genes (DMGs) was also conducted. This research is of great significance as the basis of further investigations aimed at developing strategies for the protection of the giant panda.


Physiological and biochemical indicators in blood change with age

As shown in Figure 1, the physiological and biochemical indexes in the blood of the giant pandas of different ages were systematically determined and analyzed (Supplementary Table 1). The white blood cell (WBC) and neutrophil (NEU) counts and hemoglobin (HGB) concentrations increased significantly (P < 0.05) young to adult, while decrease significantly from adult to old (Table 1). Furthermore, the levels of chloride ions (Cl-) in the blood decreased significantly with age (P < 0.05), while the levels of glutamate transaminase (AST) and glutamyl transpeptidase (GGT) increased significantly (P < 0.05). Some physiological and biochemical indicators related to metabolism in blood were also found to change significantly with age. Blood glucose (GLU) levels decreased with age, while triglycerides (TG) and apolipoprotein B (APOB) levels increased (Table 1). Furthermore, cholesterol (CHOL) levels were significantly higher in old giant pandas compared with those in adult giant pandas (P < 0.05).

Schematic diagram of experimental process design.

Figure 1. Schematic diagram of experimental process design.

Table 1. Blood physiological and biochemical indicators closely related to age.

ProjectYoungAdultOldPYoung vs. AdultYoungvs. OldAdultvs. Old
WBC (109/L)6.14±0.217.55±0.396.04±0.320.002**0.002**0.7970.001**
NEU (109/L)4.34±0.215.71±0.344.84±0.230.002**0.000**0.1310.024*
NEUT (%)66.55±2.8475.59±1.9378.06±0.830.000**0.007**0.000**0.456
HGB (U/L)118.29±1.52125.44±3.77115.04±2.000.013*0.037*0.2740.003**
Na+ (mmol/L)126.41±0.44125.85±0.68124.80±0.510.0680.4800.022*0.195
Cl- (mmol/L)96.09±0.4494.96±0.7093.25±0.480.000**0.1510.000**0.035*
AST (mmol/L)56.43±1.9763.31±2.6370.65±2.540.000**0.0590.000**0.047*
GGT (mmol/L)6.57±0.275.75±0.5110.65±1.120.000**0.4870.000**0.000**
GLU (mmol/L)4.671±0.124.581±0.274.01±0.140.006**0.7070.002**0.022*
TG (mmol/L)1.544±0.0921.911±0.0972.098±0.0970.000**0.015*0.000**0.215
CHOL (mmol/L)5.594±0.1834.841±0.285.728±0.3580.1150.0860.7230.047*
HDL (mmol/L)3.769±0.0762.945±0.1113.019±0.0960.000**0.000**0.000**0.602
LDL (mmol/L)3.312±0.1372.846±0.2233.463±0.2580.1610.1490.5860.060
APOA1 (g/L)0.696±0.0210.632±0.030.66±0.0220.1870.0770.2450.447
APOB (g/L)0.02±0.0020.032±0.0050.033±0.0030.008**0.019*0.004**0.867
Note: Most commonly used hematologic and biochemical parameters were tested. The marked * was a significant difference (P<0. 05), and the marked ** was an extremely significant difference (P<0.01). The unmarked letter indicated no statistical difference.

Identify the DEGs

To identify of the DEGs between different ages, we constructed cDNA libraries derived from blood samples to obtain the normalized expression of each gene, using |log2FoldChange| > 1 and false discovery rate (FDR) <0.05 as the standard for screening DEGs. Eight cDNA libraries were constructed using total RNA from blood samples obtained from giant pandas, and a total of 657,404,608 reads were obtained. Reads from the samples covered an average of 89.80% (88.45%–90.86%) of the giant panda reference genome [16], and an average of 86.33% (84.0%–87.8%) of the reads were compared only once (Supplementary Table 2). We used the expression levels of all genes for principal component analysis (PCA) and clustering heat map analysis. As shown in Figure 2A, all samples from the young group and two samples from the adult group were clustered together, and one sample from the adult group and all samples from the old group were clustered into a separate cluster. PCA results showed that all samples in the young group were separated from the samples in the other groups, and the samples in the adult and old groups could not be distinguished effectively (Figure 2B). In total, we detected 257, 20, and 744 DEGs in the Young vs. Adult, Adult vs. Old, and Young vs. Old groups (Figure 2C). The Young vs. Adult and Young vs. Old groups had the highest proportion of overlapping DEGs, while the Young vs. Old group has the most unique DEGs (Figure 2D).

Profile of the giant panda blood transcriptome. (A) Heat map plot of all genes using TPM expression value of genes by adopting hierarchical clustering method. (B) PCA of all genes using TPM expression value of genes by adopting hierarchical clustering method. (C) The histogram of the number of DEGs in each group. The purple box represents the number of upregulated DEGs, and the blue box represents the number of downregulated DEGs. (D) Venn diagram of the number of annotated DEGs in each group.

Figure 2. Profile of the giant panda blood transcriptome. (A) Heat map plot of all genes using TPM expression value of genes by adopting hierarchical clustering method. (B) PCA of all genes using TPM expression value of genes by adopting hierarchical clustering method. (C) The histogram of the number of DEGs in each group. The purple box represents the number of upregulated DEGs, and the blue box represents the number of downregulated DEGs. (D) Venn diagram of the number of annotated DEGs in each group.

Functional analysis of DEGs

Clarification of the function of DEGs is of considerable significance in understanding the physiological changes that occur with age in giant pandas. The protein interaction network (PPI) constructed based on the DEGs showed that BUB1 (nodes = 26) and CCNB1 (nodes = 26) had the most nodes (Supplementary Figure 1). As shown in Figure 3A, all the DEGs were divided into eight profiles by STEM according to the transcripts per million (TPM) in the three stages. However, only Profile7, Profile0, Profile1, and Profile6 showed significantly enrichment (P < 0.05). Profile7 contained the largest number of DEGs; however, due to the scattered functions of the DEGs, no significant KEGG pathway was found to be enriched in Profile7 (Supplementary Table 3). As shown in Figure 3B, KEGG pathways associated with metabolism, genetic information processing, environmental information processing, cellular processes, organismal systems, and human disease categories were enriched by the DEGs in the different profiles. To be more specific, the porphyrin and chlorophyll metabolism pathways belonging to the metabolism category were significantly upregulated in old giant pandas. All the DEGs enriched in the environmental information category pathways showed the lowest expression in the young giant pandas, and the DEGs enriched in the PI3K-Akt signaling and ECM-receptor interaction pathways showed the highest expression in the old giant pandas. DEGs enriched in most pathways of the cellular processes category showed the highest expression in the old giant pandas, while the expression levels of DEGs in the focal adhesion pathway increased with age. All pathways in the category of organismal systems belonged to the immune system subcategories. The DEGs in the B cell receptor signaling pathway showed the highest expression levels in the young pandas, and the other immune system-related pathways were found to be activated in the old giant pandas. The largest number of KEGG pathways belonged to the human disease category, although the expression trends varied among the different disease pathways.

The functional analysis of DEGs. (A) Identify different profile through STEM. Trend blocks with color: Trends with significant enrichment, trend blocks with similar trends have the same color. Trend blocks without color: Trends with no significant enrichment. (B) KEGG pathways significantly enriched by DEGs in each profile.

Figure 3. The functional analysis of DEGs. (A) Identify different profile through STEM. Trend blocks with color: Trends with significant enrichment, trend blocks with similar trends have the same color. Trend blocks without color: Trends with no significant enrichment. (B) KEGG pathways significantly enriched by DEGs in each profile.

Analysis of DNA methylation patterns

We constructed the first age-related dynamic methylation map of the giant panda with age and obtained a total of 830.13 Gb of original data (Supplementary Table 4). Comparison of the methylation sites in each sample revealed that the cytosine (C) methylation rate of the eight genomic DNA samples was approximately 5%, with the methylation rate at the CG site, while methylation at the C sites of CHG and CHH accounted for a smallest proportion among the C sites in the genome (Supplementary Figure 2A). In addition, the ratio of methylated CG:CG in all methylated cytosine sites and the ratio of CG-methylated cytosine in all methylated cytosine sites ranged from 78.21%–81.01%, while less than 1% of methylated cytosines were mCHG and mCHH, indicating that methylation of the giant panda genome mainly occurs at CG sites (Supplementary Table 5). The levels of cytosine methylation observed in each functional element in the eight samples were similar in three sequence environments (Supplementary Figure 2B). DNA methylation levels in introns and the 3'-untranslated region (3'-UTR) were higher in the CG sequence environment, followed by the CGI shore, repeat, and exon regions, while DNA methylation levels were lowest in promoter regions. Furthermore, in the promoter region, the proximal region had the lowest DNA methylation level, while the distal region had the highest DNA methylation level. In the CHG and CHH sequences, the level of DNA methylation in the CGI region was the highest, followed by the CGI region and the 5'-UTR region, and the methylation levels of the remaining functional elements were almost unchanged. We also analyzed the density of methylation in the gene coding region and the regions 2 kilobases (2 Kb) upstream and downstream of the gene (Supplementary Figure 2C). In the CG sequence environment, all samples showed the highest methylation levels in the region 2 K downstream of the gene, and the lowest in the region 2 K upstream. In the mCHG/mCHG sequence environment, we found that most individuals had a methylation peak in the region 2 K downstream of the gene. In contrast, the methylation levels in the other regions remained unchanged and were maintained at a low level. In the mCHH/mCHH environment, the methylation level was also high in the region 2 K downstream of the gene and remained low in the other regions.

Identification and functional analysis of DMGs

The differently methylated regions (DMRs) in each group were identified using the ‘methylKit’ package. The number of DMRs in the Young vs. Adult, Adult vs. Old, and Young vs. Old groups were 4,273, 6,544, and 8,834, respectively. In addition, the number of DMGs in the Young vs. Adult, Adult vs. Old, and Young vs. Old groups were 466, 823, and 991, respectively. In all groups, there were greater numbers of hypomethylated than hypermethylated DMRs and DMGs. The proportions of demethylated DMGs in each group were 82.5%, 93.0% and 94.5%, respectively, suggesting that methylation patterns decrease globally with age. The KEGG pathways significantly enriched by DMGs in each group are shown in Supplementary Table 6. The KEGG pathways enriched in the three groups included pathways belonging to the categories of environmental information processing, organismal systems, and human diseases. There were also pathways belonging to cellular process in the Adult vs. Old group.

Combined analysis of transcriptome and Whole Genome Bisulfite Sequencing (WGBS) data

To determine whether DNA methylation is involved in the regulation of DEG expression, we first identified the intersection of DEGs and DMGs in each age group. As shown in Figure 4A, there were three, one, and 36 genes that were both differentially methylated and differentially expressed in the Young vs. Adult, Adult vs. Old, and Young vs. Old groups, respectively (Supplementary Table 7). Furthermore, we found that the KEGG pathways enriched by DMGs in each group overlapped with the KEGG pathways enriched by DEGs in all profiles (Figure 4B). In addition, 49 DEGs that were significantly enriched in the KEGG pathways were selected from each profile based on previous research. We obtained the expression levels, CpG methylation values in the promoter region and CpG methylation values for each of these DEGs. Then, we analyzed potential correlations between the expression and CpG methylation levels of the 49 DEGs. The expression levels of CCNE1, CD79A1, IL1R1, and TCF7 were highly correlated with the degree of CpG methylation. The expression of CCNE1 was negatively correlated with CpG methylation in the promoter region (r = -0.760). The expression of CD79A, IL1R1, and TCF7 was negatively correlated with CpG methylation in the gene body region. In the PPI network constructed for all DEGs, CCNE1, CD79A, IL1R1, and TCF7 all had more than five nodes, suggesting that these genes play an important role in the age-related changes in the entire genetic network.

Combined analysis of transcriptome and methylation data. (A) Venn diagram of DEGs and DMGs in different groups. (B) Venn diagram of KEGG pathway significantly enriched by DEGs and DMGs. (C) A scatter plot and trend line (Pearson correlation) showing correlation between the log2 ratios of CCNE1, CD79A, IL1R1, and TCF7 expression from transcriptome and CpG methylation.

Figure 4. Combined analysis of transcriptome and methylation data. (A) Venn diagram of DEGs and DMGs in different groups. (B) Venn diagram of KEGG pathway significantly enriched by DEGs and DMGs. (C) A scatter plot and trend line (Pearson correlation) showing correlation between the log2 ratios of CCNE1, CD79A, IL1R1, and TCF7 expression from transcriptome and CpG methylation.


The immune-related indicators of giant pandas changed with age

Most of the immune-related physiological indicators in giant pandas indicated that immune function improved gradually from the young to the adult age group. WBCs, including NEU, monocytes (MONO), eosinophils (EOS), basophils (BAS), and lymphocytes (LY), are components of the innate and adaptive immune systems [17]. HGB is an iron-associated protein that is required to for oxygen transport in mammals [18]. Activation of the immune system leads to increased expression of transferrin receptors on NEUs, decreased serum iron levels (transferrin), and increased intracellular iron due to deposition by ferritin. T cell-dependent immune responses are reduced as a result of iron deficiency; thus, our findings suggest that HGB is also involved in the humoral immune response [19]. In our results, the number of WBC, NEU and HGB increased significantly from childhood to adulthood, but decreased significantly from adult to old age. This shows that the immune function of giant pandas improved gradually from young to adult, while decreased gradually from adult to old. Inflammation plays a crucial role in the control of pathogens and the formation of subsequent adaptive immune responses [20]. AST and GGT significantly increased as the pandas aged from adult to old, suggesting that the occurrence of inflammatory reaction gradually increased during the aging process [21]. This observation was supported by the changes in Cl- and Na+ concentrations. Cl- are the main ions in valval cells and are usually accompanied by Na+. Our results showed that the Cl- and Na+ concentrations decreased with age. The reduction of Cl- is thought to enhance the inflammatory response of endothelial cells [22].

Unique changes in metabolism-related indicators during the aging process in giant pandas

Intriguingly, aging is often accompanied by whole-body dysregulation of cholesterol metabolism [23]. A clinical manifestation of this process is an increase of low density lipoprotein cholesterol (LDL) and a decrease of high density lipoprotein (HDL) with age [24]. Compared with adult giant pandas in our study, plasma CHOL and LDL levels were obviously increased in old giant pandas. However, HDL levels increased slightly in the plasma of old giant pandas, which is in contrast to the results of studies in other mammals [24]. HDL is transported to the liver and then removed through the intestine and this is the mechanism of excess cholesterol removal in peripheral tissues; therefore, HDL levels are closely related to lifespan [25, 26]. Normally, HDL levels decreases by 1% per year, although favorable HDL characteristics are often observed in the offspring of centenarians [27], indicating that increased HDL levels confer exceptional longevity. In addition to the differences in HDL levels compared with the aging characteristics of other mammals, our results show that the GLU levels of giant pandas are also inconsistent with most mammalian aging characteristics. In aged mammals, the basal rate of gluconeogenesis is increased leading to increased blood glucose levels, while the hepatic association of glucose with glycogen is reduced [28]. In contrast to gluconeogenesis, the absorption of glucose by skeletal muscle, brain, and other energy-consuming tissues decreases with age, due to reduced insulin signaling [29], reduced insulin sensitivity [30], and reduced levels of glucose transporter [31]. This imbalance in supply and demand often leads to the occurrence of disorders of glucose metabolism in old age. The GLU levels of giant pandas in this study decreased significantly in old age, which is completely opposite to the pattern observed in other mammals. Multiple studies have shown that diet can reduce the effects of aging on cholesterol and glucose metabolism [32, 33]. Therefore, we speculate that the giant panda’s diet may be one of the reasons that its metabolic ability changes with age in a manner that is different from other mammals.

The susceptibility of giant pandas to different diseases changes with age

In order to further explore the molecular mechanism behind this change in blood physiological and biochemical indicators, we performed transcriptome sequencing on the blood of giant pandas of different age groups. However, the gene expression patterns of the three samples in the adult group are not completely similar, and the expression pattern of F08 is more similar to that of the two 28-year-old giant pandas. We are not clear what mechanism caused the specific expression pattern of F08, but such results suggest that there are differences in the expression patterns of different individuals at the same age. This may also be one of the reasons why there are fewer DEGs between the blood transcriptomes of adult and old giant pandas. Moreover, we used STEM software to divide the expression patterns of all age-related DEGs. In the KEGG pathway enriched by DEGs with different expression trends, most of belonged to the category of human diseases, indicating age-related variation in the risk of pandas suffering from diseases. Du et al.'s report [10] also showed that 27 of the 35 KEGG pathways enriched in DEGs between adult and old giant pandas are disease-related. In the human disease category, most pathways were enriched in cancer-related subcategories. In addition to the cancer pathway, the expression of DEGs in pathways in other cancer subcategories decreased continuously with age. And some recent research reports have shown the existence of ovarian cancer and testicular cancer in old giant pandas [34]. These results suggest that the risk of cancer in giant panda changes with age. In addition, the expression levels of DEGs enriched in the primary immunodeficiency pathway decreased significantly from young to adult and then stabilized, indicating that the innate immune system of giant pandas gradually improved with development from the young to the adult stage. The expression levels of DEGs associated with infection-related pathways were lowest in the young and usually peaked in old age. Unlike diseases such as cancer, infectious diseases are an important threat to the giant panda population [35, 36]. These results suggest that the prevention of infectious diseases may be important in extending its the lifespan of giant pandas in old age.

Transcriptome analysis indicates the immune and metabolic functions of giant pandas change with age

In accordance with the changes in blood physiological and biochemical indicators, DEGs with different trends were also significantly enriched in the pathways related to immune system and metabolism. Among the immune-related pathways, DEGs other than the B cell receptor signaling pathway showed the highest expression levels in old giant pandas. Multiple genes in the CD family, including CD19, CD22, CD72, CD79A, and CD79B, were enriched in B cell receptor signaling pathway. CD79A and CD19 were reported to inhibit the conversion of pro-B cells to pre-B cells [37]. CD19 [38] and CD22 [39] inhibit B cell proliferation under certain conditions, with significant downregulation of their expression in adults suggesting that this this inhibitory effect is weakened. These DEGs (CD79A, CD19, CD22) have a negative regulatory effect on B cell maturation, and their upregulation at an early age also indicates that the immune system of the young giant panda is not yet fully mature. Hematopoietic stem cells (HSC) refer to cells that have not yet matured and are the origin of all hematopoietic cells and immune cells. HSC can differentiate into red blood cells, WBC, and platelets. The function of HSC deteriorates with age with a concomitant effect on the function of the regenerative hematopoietic system [40]. However, in this study, hematopoietic cell profiles, platelet activation, complement and coagulation cascade pathways were all upregulated in old age. Blood physiological and biochemical indicator analyses also showed that the numbers of red blood cells and platelets were slightly reduced in the old giant pandas (P > 0.05). Collectively, these results suggest that there may be a mechanism that protects the function against the aging of HSC in pandas. Interestingly, we also found that the porphyrin and chlorophyll metabolism pathways were significantly upregulated in the old giant pandas. Chlorophyll plays role in reducing plasma levels of GLU [41] and CHOL [42], which may explain the difference in the trend of GLU and CHOL indicators between giant pandas and other mammals with age.

DNA methylation is involved in the changes of immune and metabolic functions of giant pandas with age

To further explore the regulatory mechanisms underlying these changes in immune and metabolic functions, we constructed the first dynamic map of blood DNA methylation in the giant panda. Previous studies have shown that CG methylation is the major type of methylation in mammals, including humans [43], mice [44], and poultry [45]. In accordance with this, we showed that the ratio of all methylated C sites and mCG in all methylated cytosine sites was 78.21%–81.01%, while less than 1% of the methylated cytosines were mCHG and mCHH. DNA methylation decreases globally with age. These results are similar to those found in other mammals [12, 46], indicating the DNA methylation patterns of mammals are similar and conserved. DNA methylation is thought to be involved in the regulation of gene expression, usually by inhibition [47]. However, only few of DEGs in each group are also differently methylated, this may be due to the fact that gene expression is also regulated in many ways other than methylation [48, 49]. In this study, we showed overlaps between the DEGs and DMGs in each group as well in the KEGG pathways enriched by DEGs and DMGs. CCNE1, CD79A, IL1R1, and TCF7 were identified as DEGs with multiple nodes in the PPI, and their methylation and expression showed a significant negative correlation. CCNE1 as a gene related to proliferation is increased in normal aging and carcinogenic people [50]. CD79A plays an important role in immunity. Studies in humans have shown that ILI1R1 increases with age [51], and its increase is associated with high blood pressure, diabetes, sudden heart failure, and a higher risk of mortality [52]. In particular, the age-related inhibition of TCF7 expression by methylation not only plays an important role in immune enhancement [53] and tumor suppression [54], but also in balancing blood glucose levels [55] and hematopoiesis [56]. More importantly, as a transcription factor, the methylation of TCF7 may affect the expression levels of other functional genes. These results indicate that DNA methylation is involved the process by which biological functions alter with age in the giant panda.

In general, our experimental system explains the age-related changes in blood physiological and biochemical indicators, transcriptome and methylation group map in giant pandas. Such results indicate that giant pandas have unique physiological characteristics at each age stage, so we need to differentiate the breeding of giant pandas of different ages according to their physiological characteristics. Moreover, Our results indicate that gene methylation participates in changes in immune and metabolic functions with age by regulating gene expression. These results provide an important reference in the development of strategies for the protection of giant pandas.

Materials and Methods

Ethics statement

All experiments involving animals in this study were conducted in strict accordance with the animal protection laws of the People’s Republic of China (October 26, 2018). Protocols for animal studies were approved by the Care and Use of Laboratory Animals of the Animal Ethics Committee of China Conservation and Research Center for the Giant Panda (Dujiangyan, China) (Approval No. 20180212) and the Key Laboratory of State Forestry and Grassland Administration on Conservation Biology of Rare Animals in the Giant Panda National Park.

Detection of physiological and biochemical indexes in blood

Thirty captive giant pandas aged 1–28 years were raised in China Conservation and Research Center for Giant Panda (Dujiangyan, China). These pandas were healthy (except for a few giant pandas with medical history) with typical spirit, appetite, and defecation patterns (Supplementary Table 8). During the period between 2017 and 2019, and 70 blood samples were collected as follows: 28 from the young group (1.5–5.5 years old), 16 from the adult group (7–20 years old), and 26 from the old group (21–28 years old). All samples were sent to Dujiangyan People’s Hospital for physiological and biochemical tests using the bc-5800 hematology analyzer (Mindray Technology Co., Ltd., China) and AU 2700 chemistry analyzer (Olympus, Japan), respectively.

Sample collection, library preparation, and sequencing

Blood (10 mL) were sampling from giant pandas in the young (n = 3), adult (n = 3), and old (n = 2) groups (Table 2). Total RNA and DNA were extracted from the blood samples using TRIzol Reagent and Blood tissue DNA extraction (Qiagen, Hilden, Germany), respectively. Qualified RNA and DNA were sent to Novogene (Tianjin, China) where libraries were generated. Subsequently, transcriptome and genome-wide methylation sequencing was conducted using the Illumina HiSeq 2000 and Illumina HiSeq 2500, respectively. All sequencing data have been submitted to the NCBI database (

Table 2. Information on sequencing samples used in our study.

Sample nameAge (year)GenderGenetic backgroundGroup
F011.5FemaleF4 of A × F1 of BYoung
M011.5MaleF3 of A × F1 of BYoung
M033MaleF3 of A × F2 of BYoung
M077MaleF1 of B × CAdult
F088FemaleF2 of A × DAdult
F1212FemaleF2 of A × EAdult
F2828FemaleA × FOld

RNA-seq data analysis

Clean reads were compared to the giant panda reference genome (GCF_000004335.2) using HISAT2 (version 2.1.0). The output SAM (sequencing alignment/mapping) file was converted to a BAM (binary alignment/mapping) file and sorted using SAMtools (0.1.19-44428cd). FeatureCounts (version 2.0.0) was used obtain the number of counts per gene. We divided the samples into Young vs. Adult, Young vs. Old and Adult vs. Old groups, and wrote a grouping file for each group. DEseq2 in R studio was used to read the data and group. After determination of the different genes in each group, the files were screened according to the FDR values. The DEGs were filtered based on FDR<0.05 and |log2FC|>1. String ( was used to construct a protein-protein interaction network for all the DEGs. The DEGs in each group were combined and their expression levels were input into STEM [57] software for analysis. A fold-change in expression >1.2 was considered to be a trend, and the genes were eventually divided into eight profiles. The KEGG functions of the DEGs in each profile were analyzed in KOBAS3.0 (

BS-seq data analysis

Clean reads were deduplicated and aligned against the bisulfite-converted reference sequence using Bismark Bisulfite Mapper (v0.15.0). Methylation status was determined using the Bismark Methylation Extractor Script and only uniquely aligned reads, the number of methylated and unmethylated CpG, and non-CpG (CHG and CHH, H representing A/C/T) sites were counted for each individual cytosine. ‘MethylKit’ in the R package was used to identify the differently methylated CpGs (DMCs) and DMRs. Genes that overlapped with DMR were defined as DMGs. DAVID software ( was used to analyze the function of DMGs, and the results were considered significant when P < 0.05. Genome-wide DNA methylation analysis was conducted according to the annotated protein-coding genes in giant panda, which was downloaded from the Ensembl ( database. The region from the transcription start site (TSS) to the transcription termination site (TTS) was defined as the gene body region. The genomic region 2.2 kb upstream of the TSS to 500 bp downstream of the TSS was considered as the proximal promoter region. CGI shores were defined as regions 2 kb in length adjacent to CGIs.

Combined analysis of transcriptome and WGBS data

The correlation coefficient and P-value between the log2 the value of Transcripts per million reads (TPM) and methylation of the gene body and promoter regions were calculated using Hmisc ( in the R package.

Author Contributions

Sampling: Xiaoyu Huang, Qingyuan Ouyang, Mingxia Ran, Linhua Deng, Guo Li, Tao Deng, Ming He, Ti Li, Haidi Yang, Guiquan Zhang; formal analysis: Xiaoyu Huang, Qingyuan Ouyang, Mingxia Ran; Funding acquision: Xiaoyu Huang; project administration: Jiwen Wang, Changjun Zeng, Heming Zhang; writing—original draft, Xiaoyu Huang, Qingyuan Ouyang, Minxia Ran; writing—review and editing, Bo Zeng, Shenqiang Hu, Mingyao Yang, Jiwen Wang; All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest

The authors declare no conflicts financial interest.


This work was supported by CCRCGP181910.


  • 1. Mainka SA, He T, Chen M, Dierenfeld ES. Hematologic and serum biochemical values for healthy captive giant pandas (Ailuropoda melanoleuca) at the Wolong reserve, Sichuan, China. J Zoo Wildl Med. 1995; 26:377–81.
  • 2. Nikolich-Žugich J. The twilight of immunity: emerging concepts in aging of the immune system. Nat Immunol. 2018; 19:10–19. [PubMed]
  • 3. Goronzy JJ, Weyand CM. Understanding immunosenescence to improve responses to vaccines. Nat Immunol. 2013; 14:428–36. [PubMed]
  • 4. Ford ES, Giles WH, Dietz WH. Prevalence of the metabolic syndrome among US adults: findings from the third National Health and Nutrition Examination Survey. JAMA. 2002; 287:356–9. [PubMed]
  • 5. Morley JE. Diabetes and aging: epidemiologic overview. Clin Geriatr Med. 2008; 24:395–405. [PubMed]
  • 6. Barzilai N, Huffman DM, Muzumdar RH, Bartke A. The critical role of metabolic pathways in aging. Diabetes. 2012; 61:1315–22. [PubMed]
  • 7. Yuan Y, Hakimi P, Kao C, Kao A, Liu R, Janocha A, Boyd-Tressler A, Hang X, Alhoraibi H, Slater E, Xia K, Cao P, Shue Q, et al. Reciprocal changes in phosphoenolpyruvate carboxykinase and pyruvate kinase with age are a determinant of aging in caenorhabditis elegans. J Biol Chem. 2016; 291:1307–19. [PubMed]
  • 8. Feng Z, Hanson RW, Berger NA, Trubitsyn A. Reprogramming of energy metabolism as a driver of aging. Oncotarget. 2016; 7:15410–20. [PubMed]
  • 9. Du L, Li W, Fan Z, Shen F, Yang M, Wang Z, Jian Z, Hou R, Yue B, Zhang X. First insights into the giant panda (Ailuropoda melanoleuca) blood transcriptome: a resource for novel gene loci and immunogenetics. Mol Ecol Resour. 2015; 15:1001–13. [PubMed]
  • 10. Du L, Liu Q, Shen F, Fan Z, Hou R, Yue B, Zhang X. Transcriptome analysis reveals immune-related gene expression changes with age in giant panda (Ailuropoda melanoleuca) blood. Aging (Albany NY). 2019; 11:249–262. [PubMed]
  • 11. Zhang X, Hu M, Lyu X, Li C, Thannickal VJ, Sanders YY. DNA methylation regulated gene expression in organ fibrosis. Biochim Biophys Acta Mol Basis Dis. 2017; 1863:2389–97. [PubMed]
  • 12. Kochmanski J, Marchlewicz EH, Cavalcante RG, Sartor MA, Dolinoy DC. Age-related epigenome-wide DNA methylation and hydroxymethylation in longitudinal mouse blood. Epigenetics. 2018; 13:779–92. [PubMed]
  • 13. Sen P, Shah PP, Nativio R, Berger SL. Epigenetic mechanisms of longevity and aging. Cell. 2016; 166:822–39. [PubMed]
  • 14. Ren J, Shen F, Zhang L, Sun J, Yang M, Yang M, Hou R, Yue B, Zhang X. Single-base-resolution methylome of giant panda’s brain, liver and pancreatic tissue. PeerJ. 2019; 7:e7847. [PubMed]
  • 15. You Y, Bai C, Liu X, Xia M, Jia T, Li X, Zhang C, Chen Y, Zhao S, Wang L, Wang W, Yin Y, Xiu Y, et al. Genome-wide analysis of methylation in giant pandas with cataract by methylation-dependent restriction-site associated DNA sequencing (MethylRAD). PLoS One. 2019; 14:e0222292. [PubMed]
  • 16. Li R, Fan W, Tian G, Zhu H, He L, Cai J, Huang Q, Cai Q, Li B, Bai Y, Zhang Z, Zhang Y, Wang W, et al. The sequence and de novo assembly of the giant panda genome. Nature. 2010; 463:311–17. [PubMed]
  • 17. Carrick JB, Begg AP. Peripheral blood leukocytes. Vet Clin North Am Equine Pract. 2008; 24:239–59. [PubMed]
  • 18. Gell DA. Structure and function of haemoglobins. Blood Cells Mol Dis. 2018; 70:13–42. [PubMed]
  • 19. Schrama JW, Schouten JM, Swinkels JW, Gentry JL, de Vries Reilingh G, Parmentier HK. Effect of hemoglobin status on humoral immune response of weanling pigs differing in coping styles. J Anim Sci. 1997; 75:2588–96. [PubMed]
  • 20. Cronkite DA, Strutt TM. The regulation of inflammation by innate and adaptive lymphocytes. J Immunol Res. 2018; 2018:1467538. [PubMed]
  • 21. Bronikowski AM, Carter PA, Morgan TJ, Garland T Jr, Ung N, Pugh TD, Weindruch R, Prolla TA. Lifelong voluntary exercise in the mouse prevents age-related alterations in gene expression in the heart. Physiol Genomics. 2003; 12:129–38. [PubMed]
  • 22. Yang H, Huang LY, Zeng DY, Huang EW, Liang SJ, Tang YB, Su YX, Tao J, Shang F, Wu QQ, Xiong LX, Lv XF, Liu J, et al. Decrease of intracellular chloride concentration promotes endothelial cell inflammation by activating nuclear factor-κB pathway. Hypertension. 2012; 60:1287–93. [PubMed]
  • 23. Mc Auley MT, Mooney KM. Lipid metabolism and hormonal interactions: impact on cardiovascular disease and healthy aging. Expert Rev Endocrinol Metab. 2014; 9:357–67. [PubMed]
  • 24. Wilson PW, Anderson KM, Harris T, Kannel WB, Castelli WP. Determinants of change in total cholesterol and HDL-C with age: the framingham study. J Gerontol. 1994; 49:M252–57. [PubMed]
  • 25. Ferrara A, Barrett-Connor E, Shan J. Total, LDL, and HDL cholesterol decrease with age in older men and women. The rancho bernardo study 1984-1994. Circulation. 1997; 96:37–43. [PubMed]
  • 26. Rader DJ, Schaefer JR, Lohse P, Ikewaki K, Thomas F, Harris WA, Zech LA, Dujovne CA, Brewer HB Jr. Increased production of apolipoprotein A-I associated with elevated plasma levels of high-density lipoproteins, apolipoprotein A-I, and lipoprotein A-I in a patient with familial hyperalphalipoproteinemia. Metabolism. 1993; 42:1429–34. [PubMed]
  • 27. Barzilai N, Gabriely I, Gabriely M, Iankowitz N, Sorkin JD. Offspring of centenarians have a favorable lipid profile. J Am Geriatr Soc. 2001; 49:76–79. [PubMed]
  • 28. Satrústegui J, Cuezva JM, Machado A. Increased basal gluconeogenesis in the aged rat. FEBS Lett. 1986; 197:159–63. [PubMed]
  • 29. Muzumdar R, Ma X, Atzmon G, Vuguin P, Yang X, Barzilai N. Decrease in glucose-stimulated insulin secretion with aging is independent of insulin action. Diabetes. 2004; 53:441–46. [PubMed]
  • 30. Karakelides H, Irving BA, Short KR, O’Brien P, Nair KS. Age, obesity, and sex effects on insulin sensitivity and skeletal muscle mitochondrial function. Diabetes. 2010; 59:89–97. [PubMed]
  • 31. dos Santos JM, Benite-Ribeiro SA, Queiroz G, Duarte JA. The effect of age on glucose uptake and GLUT1 and GLUT4 expression in rat skeletal muscle. Cell Biochem Funct. 2012; 30:191–97. [PubMed]
  • 32. Vartiainen E, Laatikainen T, Peltonen M, Juolevi A, Männistö S, Sundvall J, Jousilahti P, Salomaa V, Valsta L, Puska P. Thirty-five-year trends in cardiovascular risk factors in Finland. Int J Epidemiol. 2010; 39:504–18. [PubMed]
  • 33. Kalant N, Stewart J, Kaplan R. Effect of diet restriction on glucose metabolism and insulin responsiveness in aging rats. Mech Ageing Dev. 1988; 46:89–104. [PubMed]
  • 34. Gao Q, Wang C, Li D, Zhang H, Deng L, Li C, Chen Z. A case of giant panda ovarian cancer diagnosis and histopathology. BMC Vet Res. 2018; 14:311. [PubMed]
  • 35. Wang T, Xie Y, Zheng Y, Wang C, Li D, Koehler AV, Gasser RB. Parasites of the giant panda: a risk factor in the conservation of a species. Adv Parasitol. 2018; 99:1–33. [PubMed]
  • 36. Jin Y, Zhang X, Ma Y, Qiao Y, Liu X, Zhao K, Zhang C, Lin D, Fu X, Xu X, Wang Y, Wang H. Canine distemper viral infection threatens the giant panda population in China. Oncotarget. 2017; 8:113910–19. [PubMed]
  • 37. von Muenchow L, Tsapogas P, Albertí-Servera L, Capoferri G, Doelz M, Rolink H, Bosco N, Ceredig R, Rolink AG. pro-B cells propagated in stromal cell-free cultures reconstitute functional b-cell compartments in immunodeficient mice. Eur J Immunol. 2017; 47:394–405. [PubMed]
  • 38. Hatterer E, Barba L, Noraz N, Daubeuf B, Aubry-Lachainaye JP, von der Weid B, Richard F, Kosco-Vilbois M, Ferlin W, Shang L, Buatois V. Co-engaging CD47 and CD19 with a bispecific antibody abrogates B-cell receptor/CD19 association leading to impaired B-cell proliferation. MAbs. 2019; 11:322–334. [PubMed]
  • 39. Santos L, Draves KE, Boton M, Grewal PK, Marth JD, Clark EA. Dendritic cell-dependent inhibition of B cell proliferation requires CD22. J Immunol. 2008; 180:4561–69. [PubMed]
  • 40. Bartek J, Hodny Z. Ageing: old blood stem cells feel the stress. Nature. 2014; 512:140–41. [PubMed]
  • 41. Li Y, Cui Y, Hu X, Liao X, Zhang Y. Chlorophyll supplementation in early life prevents diet-induced obesity and modulates gut microbiota in mice. Mol Nutr Food Res. 2019; 63:e1801219. [PubMed]
  • 42. Takamoto H, Eguchi K, Kawabata T, Fujiwara Y, Takeya M, Tsukamoto S. Inhibitors for cholesterol ester accumulation in macrophages from Chinese cabbage. Biosci Biotechnol Biochem. 2015; 79:1315–19. [PubMed]
  • 43. Wijnands KP, Chen J, Liang L, Verbiest MM, Lin X, Helbing WA, Gittenberger-de Groot AC, van der Spek PJ, Uitterlinden AG, Steegers-Theunissen RP. Genome-wide methylation analysis identifies novel CpG loci for perimembranous ventricular septal defects in human. Epigenomics. 2017; 9:241–251. [PubMed]
  • 44. Tamai Y, Takemoto Y, Matsumoto M, Morita T, Matsushiro A, Nozaki M. Sequence of EndoA gene encoding mouse cytokeratin and its methylation state in the CpG-rich region. Gene. 1991; 104:169–76. [PubMed]
  • 45. Li J, Li R, Wang Y, Hu X, Zhao Y, Li L, Feng C, Gu X, Liang F, Lamont SJ, Hu S, Zhou H, Li N. Genome-wide DNA methylome variation in two genetically distinct chicken lines using MethylC-seq. BMC Genomics. 2015; 16:851. [PubMed]
  • 46. Johnson ND, Huang L, Li R, Li Y, Yang Y, Kim HR, Grant C, Wu H, Whitsel EA, Kiel DP, Baccarelli AA, Jin P, Murabito JM, Conneely KN. Age-related DNA hydroxymethylation is enriched for gene expression and immune system processes in human peripheral blood. Epigenetics. 2020; 15:294–306. [PubMed]
  • 47. Ehrlich M, Lacey M. DNA methylation and differentiation: silencing, upregulation and modulation of gene expression. Epigenomics. 2013; 5:553–68. [PubMed]
  • 48. Lawrence M, Daujat S, Schneider R. Lateral thinking: how histone modifications regulate gene expression. Trends Genet. 2016; 32:42–56. [PubMed]
  • 49. Zhao BS, Roundtree IA, He C. Post-transcriptional gene regulation by mRNA modifications. Nat Rev Mol Cell Biol. 2017; 18:31–42. [PubMed]
  • 50. Leiszter K, Galamb O, Sipos F, Krenács T, Veres G, Wichmann B, Kalmár A, Patai ÁV, Tóth K, Valcz G, Molnár B, Tulassay Z. Sporadic colorectal cancer development shows rejuvenescence regarding epithelial proliferation and apoptosis. PLoS One. 2013; 8:e74140. [PubMed]
  • 51. Coglianese EE, Larson MG, Vasan RS, Ho JE, Ghorbani A, McCabe EL, Cheng S, Fradley MG, Kretschman D, Gao W, O'Connor G, Wang TJ, Januzzi JL. Distribution and clinical correlates of the interleukin receptor family member soluble ST2 in the Framingham Heart Study. Clin Chem. 2012; 58:1673–81. [PubMed]
  • 52. Ho JE, Sritara P, deFilippi CR, Wang TJ. Soluble ST2 testing in the general population. Am J Cardiol. 2015; 115:22B–5B. [PubMed]
  • 53. Barra MM, Richards DM, Hansson J, Hofer AC, Delacher M, Hettinger J, Krijgsveld J, Feuerer M. Transcription factor 7 limits regulatory T cell generation in the thymus. J Immunol. 2015; 195:3058–70. [PubMed]
  • 54. Wu B, Chen M, Gao M, Cong Y, Jiang L, Wei J, Huang J. Down-regulation of lncTCF7 inhibits cell migration and invasion in colorectal cancer via inhibiting TCF7 expression. Hum Cell. 2019; 32:31–40. [PubMed]
  • 55. Kaur K, Vig S, Srivastava R, Mishra A, Singh VP, Srivastava AK, Datta M. Elevated Hepatic miR-22-3p Expression Impairs Gluconeogenesis by Silencing the Wnt-Responsive Transcription Factor Tcf7. Diabetes. 2015; 64:3659–69. [PubMed]
  • 56. Wu JQ, Seay M, Schulz VP, Hariharan M, Tuck D, Lian J, Du J, Shi M, Ye Z, Gerstein M, Snyder MP, Weissman S. Tcf7 is an important regulator of the switch of self-renewal and differentiation in a multipotential hematopoietic cell line. PLoS Genet. 2012; 8:e1002565. [PubMed]
  • 57. Ernst J, Bar-Joseph Z. STEM: a tool for the analysis of short time series gene expression data. BMC Bioinformatics. 2006; 7:191. [PubMed]