Blood transcriptome analysis revealing aging gene expression profiles in red panda

Sample collection

Peripheral blood samples of twelve red pandas were obtained during routine examinations at the Chengdu Research Base of Giant Panda Breeding, China (age and gender shown in Table S1). All red pandas were sampled in 2018. This study was approved by the Chengdu Institute of Biology Animal Use Ethics Committee (NO. 2018008). After collection, fresh blood samples were immediately stored in the PAXgene Blood RNA Tubes, and then stored at the −80 °C refrigerator until RNA extraction.

Library preparation and sequencing

Threefold red blood cell lysis buffer (TIANGEN, Beijing, China) was added to lyse erythrocytes and Sorvall Legend X1 (Thermo Scientific, Waltham, MA, USA) was used to centrifugate peripheral blood mononuclear cells (PBMCs). Total RNA was extracted by the M5 Universal RNA Mini Kit tissue/Cell RNA rapid Extraction Kit (Mei5 Biotechnology Co. Ltd, Beijing, China). After assessment and purification of RNA quality, samples with RIN (RNA integrity number) values higher than 7.5 were used for library construction and sequencing. Sequencing libraries were constructed by using NEBNext^® Ultra™ RNA Library Prep Kit for Illumina^® (NEB, Ipswich, MA, USA). All libraries were sequenced by using the Illumina NovaSeq 2000 platform with a paired-end sequencing length of 150 bp (PE150). Library construction and sequencing were performed by Novogene Bioinformatics Institute (Beijing, China).

Quality control and sequence alignment

The red panda genome assembly and reference annotation were downloaded from https://www.dnazoo.org/assemblies/Ailurus_fulgens. After sequencing, we obtained 12 sets of transcriptome data. For obtaining high-quality reads, raw data were refined by Trimmomatic (Bolger, Marc & Bjoern, 2014) to remove index, adapter and low-quality sequences. Parameter were set to SLIDINGWINDOW:5:20, LEADING:5, TRAILING:5, MINLEN:50, and others used default parameters. We used FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) to generate quality reports for evaluating the quality of the raw and clean reads. High-quality paired-end reads from 12 libraries were mapped on the red panda reference genome using HISAT2 (Kim, Langmead & Salzberg, 2015), separately. We used default parameters for scoring and the matching process. SAM files, generated by Hisat2, were sorted by SAMtools (Li et al., 2009) to generate BAM files.

Differentially expressed genes analysis

Raw read counts were the input in DESeq2 (Love, Huber & Anders, 2014) for each gene and were extracted by featureCounts (Liao, Smyth & Shi, 2014). P values adjusted (padj) ≤ 0.05 & —log2 Fold Change—≥ 1 was used to determine significant DEGs.

Functional annotation and gene enrich analysis

The gene sequences were searched against NCBI non-redundant (NR) database using BLASTX (Altschul et al., 1997) with a typical cutoff E-value of 1E-5. KEGG pathway annotation was performed on KAAS. Gene ontology (GO) was applied with the Blast2GO (Conesa et al., 2005) to obtain annotations of unknown genes. To classify functions of the DEGs, we performed GO enrichment analysis using GOstats (Falcon & Gentleman, 2007) and KEGG pathway enrichment analysis using clusterProfiler (Yu et al., 2012).

Cluster analysis of Mfuzz transcriptome expression pattern

The R package Mfuzz (Kumar & Futschik, 2007) clusters and visualize genes with similar expression patterns, and explores the dynamic changes of different genes over time. Mfuzz uses a new clustering algorithm of fuzzy C-means algorithm, called soft clustering. Clustering methods are often based on hard clustering, in which a gene is assigned to a cluster. In contrast, soft clustering methods can assign a gene to several clusters, and it is more noise robust and can avoid a priori pre-filtering of genes (Futschik & Carlisle, 2005). We divided all the genes into 15 clusters, the minimum score (membership) threshold was set to 0.25. We then performed GO enrichment analysis for each cluster.

protein–protein interaction network analysis

protein–protein interaction analysis can check the relationship between proteins expressed by DEGs, and search for hub genes. We first inputted immune-related DEGs into the STRING protein interaction database (Szklarczyk et al., 2019) and then generated protein–protein interaction networks. Next, we downloaded the TSV text file and imported it into Cytoscape (Smoot et al., 2011) software to visualize the annotated protein–protein interaction network.

Transcriptome sequencing and sequence alignment

Twelve peripheral blood samples (PBMCs) of red pandas were sequenced by Illumina HiSeq™ 2000 using pair-end sequencing. 311 million raw pair-end reads and 296 million clean reads were obtained after quality control by removing adapter and low-quality reads. HISAT2 results showed that high-quality reads had a good alignment with the reference genome of the red panda, with an average alignment rate of 87.35% (Table S1). Our data have been deposited into the CNGB Sequence Archive (CNSA) (Guo et al., 2020) of the China National GeneBank DataBase (CNGBdb) (Chen et al., 2020a) with an accession number CNP0002175.

Identification of age-related differentially expressed genes

The lifespan of wild red pandas is about 12 years and up to 18 years in captivity (Ali Khan & Kamal Naidu, 1981). The samples were categorized into juvenile, adults and geriatric by age according to a previous study (Delaski, Ramsay & Gamble, 2015). Principal component analysis (PCA) results showed that juvenile, adult, and geriatric individuals could be divided into three groups in different dimensions (Fig. 1C), indicating that there were differences in gene expression patterns between the three age groups. In total, 2219 genes were significantly differentially expressed in red panda blood samples between geriatric and adult individuals. Among the 2219 DEGs, 1422 genes were up-regulated while 797 genes were down-regulated in geriatric compared to adult individuals (Fig. 1A). We found that 1690 genes were significantly differentially expressed between adult and juvenile individuals. Among the 1690 DEGs, 725 genes were up-regulated while 965 genes were down-regulated in adults compared to juveniles (Fig. 1B). We found that 563 genes were differentially expressed in the two comparisons (Fig. 1D). All DEGs are shown in Supplementary 2.

Gene ontology enrichment analysis of differentially expressed genes

To comprehensively understand the biological role of DEGs, we conducted GO enrichment analysis. Results showed that for the up-regulated genes between geriatric and adult individuals, almost all the significantly enriched GO terms in the biological process category were associated with the innate immune system (Supplementary 1 and Supplementary 2), such as neutrophil activation (GO:0042119), leukocyte activation involved in immune response (GO:0002366), leukocyte degranulation (GO:0043299), and defense response (GO:0006952). Obsolete cytoplasmic vesicle part (GO:0044433) and signaling receptor binding (GO:0005102) were the most significantly enriched GO term in the cellular component and molecular function categories, respectively. For down-regulated genes, the most significantly enriched biological process GO term was chromosome segregation (GO:0007059). Regulation of CD40 signaling pathway (GO:2000348) and B cell activation (GO:0042113) were also enriched. Condensed chromosome (GO:0000793) and molecular adaptor activity (GO:0060090) were the most significantly enriched GO term in the cellular component and molecular function categories, respectively (Supplementary 1 and Supplementary 2).

For the up-regulated genes between adults and juveniles, DNA replication (GO:0006260), chromosomal region (GO:0005654) and ATP-dependent microtubule motor activity (GO:0008574) were the most significantly enriched GO term in the biological process, cellular component and molecular function categories, respectively (Supplementary 1 and Supplementary 2). For down-regulated genes, the significantly enriched GO terms in the biological process category were mainly associated with organ development, such as animal organ development (GO:0048513), and multicellular organism development (GO:0007275). Cell periphery (GO:0071944) and sequence-specific DNA binding (GO:0043565) were the most significantly enriched GO term in the cellular component and molecular function categories, respectively (Supplementary 1 and Supplementary 2).

Pathway enrichment analysis of differentially expressed genes

Hypergeometric test with a P-value cutoff of 0.05 was used as the criteria for pathway detection. Down-regulated genes between adults and juveniles were significantly enriched in four pathways, including Parathyroid hormone synthesis, secretion and action (ko04928), Axon guidance (ko04360) and Cholinergic synapse (ko04725). Down-regulated genes between geriatric and adult individuals were significantly enriched in six pathways, including Intestinal immune network for IgA production (ko04672) and B cell receptor signaling pathway (ko04662); more details can be found in Supplementary 1.

Immune related genes

Combined with the ImmPort database, we detected 163 immunity related DEGs between geriatrics and adults, and 116 immune related DEGs between adults and juveniles. Genes (PRKCB, BLNK, IGLV2-33, BTK, CD72, IGKV5-2, IGHE, RASGRP3) related to the BCR signaling pathway were down-regulated when comparing geriatrics to adults, and almost all genes belonging to natural killer cell were up-regulated (Fig. 2). Similarly, almost all genes belonging to natural killer cell were up-regulated when comparing adults to juveniles and genes (HSPA1A, THBS1) involved in antigen processing and presentation were down-regulated (Fig. 3). The expression level of CXCR3, BMP3, and KLRC1 continued to increase with increasing age, and the expression level of BLNK, CD72, PLXNA2 continued to decrease with increasing age.

Cluster analysis of Mfuzz transcriptome expression pattern

Cluster1, Cluster3, and Cluster10 were related to leukocyte activation, and all genes experienced a small decline and then a significant rise. Cluster6 was related to humoral immunity, where genes (MEF2C, PRKCB, ATM, TFRC, BANK1, MSH2, PELI1, EXO1, STAP1) experienced a slight rise followed by a significant decline. Cluster8 was related to cellular immunity, where genes (CD8A, CD3G, IL12RB1, XCL1, CBLB, CD6, KLRK1, ZNF683, EOMES, TBX21, CCL5, PRLR, IL7R, PTPN22) experienced a significant rise followed by a slight decline. Cluster5, Cluster11, and Cluster13 were mainly involved in animal organ development, where there was a significant decline from juveniles to adults and little difference between adults and geriatrics. Cluster14 was related to response to external stimulus where there was little difference from juveniles to adults, but a sharp increase between adult and geriatric individuals (Fig. 4). More details of each cluster are shown in Supplementary 2.

Protein–protein interaction network of differentially expressed genes

We performed a protein–protein interaction analysis of 33 immune-related DEGs that were common in the two comparisons (Table 1). CXCR3, BLNK, and CCR4 were at the key position of the interaction network (Supplementary 1).

Innate immune changes

Innate immunity includes neutrophils, macrophages, dendritic cells, natural killer cells, antimicrobial agents, opsonins, and cytokines (Rosenzweig & Holland, 2011). Neutrophils are the most abundant leukocytes, defending against pathogens by degranulation and the release of neutrophil extracellular traps (Pruchniak, Arazna & Demkow, 2013). Neutrophil degranulation pathway genes, such as ITGB2, CD14, ITGAD, CXCR1, LCN2, S100A9, and SLC11A1, were at first down-regulated in juveniles and then up-regulated in geriatrics. TNFR superfamily members (TNFRSF1B, TNFRSF4, TNFRSF9, TNFRSF14) play important roles in both innate and adaptive immunity (Collette et al., 2003; He et al., 2017) and were at first down-regulated in juveniles and then up-regulated as individuals aged. TNFRSF4 contributes to the activation and survival of neutrophils (Jin et al., 2019) and TNFRSF9 is indispensable for innate immunity against Candida albicans infection (Tran, Nguyen & Kwon, 2021). TNFRSF14 is involved in innate mucosal defense as its activation leads to activation of neutrophil effector functions, including respiratory bursts, degranulation, and interleukin-8 release (Haselmayer et al., 2006; Shui & Kronenberg, 2013).

Several important genes involved in antiviral response (OASL, ISG20, STAT4) were significantly up-regulated in adults. OASL is an interferon-inducible antiviral protein that can boost innate host defense by enhancing the sensitivity of RIG-I activation (Zhu, Ghosh & Sarkar, 2015). ISG20, as an interferon-inducible 3′–5′exonuclease, inhibits replication of several human and animal RNA viruses and plays an important role in the host’s defenses against pathogens (Degols, Eldin & Mechti, 2007; Zhou et al., 2011). STAT4 promotes activation of RIG-I signaling to initiate more type I IFN production in antiviral innate immunity (Zhao et al., 2016). Natural killer cells have been known to be an essential part of the innate immune system (Averdam et al., 2009). Our results indicated that killer cell lectin-like receptors (KLRC1, KLRD1, KLRK1, KLRG1, KLRF1) were all up-regulated in adults. However, IRF5, MX1, and TLR5, playing crucial roles in innate immune were down-regulated in adults. Interferon regulatory factor 5 (IRF5), a member of the interferon regulatory factor (IRF) gene family, is a key transcription factor for the activation of innate immune responses in the Toll-like receptor signaling pathway (Takaoka et al., 2005). Myxovirus resistance 1 (MX1), an interferon-induced gene that encodes a GTPase, is considered to be an important part of innate immune responses for antiviral host defense (Spitaels et al., 2019). Toll-like receptors (TLRs) are a key component of the innate immune system that recognize microbial infection and trigger antimicrobial host defense responses (Akira & Takeda, 2004). TLR5 binding to bacterial flagellin activates signaling through the transcription factor NF-kappa B and triggers an innate immune response to the invading pathogen (Yoon et al., 2012). Taken together, we consider that the innate immunity of red pandas increases significantly in geriatric individuals, whereas its changes from the juvenile phase to adulthood remain obscure, which are probably influenced by many factors during the rapid physiological development period. Likewise, the immune changes are also complex.

Adaptive immune changes

The adaptive immune system consists of T cells, B cells and their antigen-specific receptors (TCR and BCR) (Martin, 2014). T cell receptor signaling pathway genes (IL12RB1, XCL1, TBX21, CD8A, IL7R, JAG1, CD3G) experienced an up-regulation and then a down-regulation from the juvenile phase to geriatric phase. The T cell-specific T-box transcription factor (TBX21) is vital to the regulation of the immune system because this factor induces the development and differentiation of T (H)1 as well as contributing to cell-mediated immunity against intracellular pathogens (Chung et al., 2003; Suttner et al., 2009). CD8A that encodes a critical coreceptor of cytotoxic T cells—the CD8 alpha chain of the dimeric CD8 protein—is critical for cell-mediated immune defense and T-cell development (Xu et al., 2014). Signaling mediated by interleukin 7 receptor (IL7R) is essential for normal T-cell development and homeostasis. IL7R inactivation causes severe combined immunodeficiency in humans (Jiang et al., 2005; Puel et al., 1998). Furthermore, genes (THEMIS, PTPRC, BMP4, PRDM1, CD6, CD28, EOMES, PRLR, IFNG) related to T cell activation were up-regulated in adults. Protein tyrosine phosphatase receptor type C (PTPRC), also known as CD45, is an essential regulator of T and B cell antigen receptor-mediated activation (Al Barashdi et al., 2021). CD6 facilitates adhesion between T cells and antigen-presenting cells and modulates thymocyte selection and T lymphocyte activation (Santos, Oliveira & Carmo, 2016). We further found Interferon-regulatory factor 4 (IRF4) was down-regulated in geriatric individuals. IRF4 is induced after T cell receptor (TCR) stimulation, and the level of its expression correlates with the strength of TCR signal (Harberts et al., 2021). The transcription factor IRF4 controls activation and differentiation of CD8+ T cells and is essential for the function of CD8+ memory T cells (Yao et al., 2013). Consequently, our results showed that the cellular immunity of red pandas at first increased in juveniles and then decreased with increasing age

B cell receptor signaling pathway genes (STAP1, PRKCB, MEF2C, BANK1) experienced an up-regulation then down-regulation from juvenile to geriatric individuals. PRKCB plays a central role in propagating NF-Kappa B signaling and cell proliferation downstream of the BCR and is vital to the formation of germinal centers and plasma cells (Tsui et al., 2018). MEF2C has a critical function in mature B cell survival and proliferation after BCR stimulation (Wilker et al., 2008). B cells can differentiate into plasma cells and secrete immunoglobulins against the pathogen (Geng et al., 2020). We found genes (ATAD5, IGJ, IGHE) related to immunoglobulin were up-regulated in adults. ATAD5 deficiency has been shown to reduce B cell division and Igh recombination (Zanotti et al., 2015). Immunoglobulin J chain (IgJ) is a small polypeptide that serves to regulate polymer formation and secretory activity of Immunoglobulin A (IgA) and Immunoglobulin M (IgM) (Brandtzaeg, 1983). IgE is produced primarily in response to parasitic infections and can also be involved in allergic reactions (Larsson et al., 2020). IGHE encodes a constant region of the IgE heavy chain (Wagner, 2009). Genes (BLNK, MYB, BTK, FCRL1, FCRL5, MSH2, PAX5, BCL11A) related to B cell activation were down-regulated in geriatric individuals. PAX5 is essential for many aspects of B cell development including immunoglobulin rearrangement, pre-B cell receptor signaling and maintaining cell identity in mature B cells (Holmes, Pridans & Nutt, 2008). We found that BLNK was a hub gene in our protein–protein interaction analysis. BLNK is a central linker protein connecting B cell receptor-associated kinases with multiple signaling pathways, which regulates the pre-BI to pre-BII transition and Ig light chain gene rearrangement, playing a non-redundant key function in B cell development (Lagresle-Peyrou et al., 2014). This suggests that BLNK may play a critical role in the progression of age-related immune changes in the red panda. These results show that the humoral immunity of red pandas increases then decreases with age. Taken together, we conclude that the adaptive immunity of red pandas increases then decreases with age.

Gene expression changes associated with disease

The main cause of death in adult and geriatric red pandas was cardiovascular disease, and the proportion of death from cardiovascular disease in juveniles (31–365 days), adults (1–10 years) and geriatrics (>10 years) was 2.94%, 20.65% and 24.62%, respectively (Delaski, Ramsay & Gamble, 2015). Age is the biggest risk factor for cardiovascular disease, and the prevalence increases sharply with age. In geriatric individuals, we observed an up-regulation of PTX3, a prototypic long pentraxin, which is produced by somatic cells and innate immune cells after pro-inflammatory stimulation. In peripheral blood tissue, PTX3 was found to be significantly elevated in human patients with heart failure and acute myocardial infarction, contributing to the occurrence of cardiovascular events and cardiovascular risk factors (Inforzato et al., 2015; Liu et al., 2011; Peri et al., 2000). Therefore, PTX3 may be an effective target point for the prevention and treatment of cardiovascular disease in red pandas. The aging process and its related diseases often cause energy metabolism disorder and oxidative damage (Mattson, Maudsley & Martin, 2004). The most significant up-regulated gene, CP, between geriatric and adult red pandas, encodes a serum ferroxidase ceruloplasmin that is related to iron metabolism (Hellman & Gitlin, 2002). Clinical and epidemiological studies have shown that iron metabolism plays an important role in multiple aging-associated disorders, especially in cardiovascular disease (Roijers et al., 2011). Higher levels of serum ceruloplasmin, as an independent risk factor for cardiovascular disease, can be a marker for metabolic stresses associated with metabolic syndrome (Kim et al., 2002).

Vaccines appear as the most efficient and economical medical interventions against infectious diseases and have also been used to defend against pathogen infections in red pandas (Dumpa et al., 2019; Qin et al., 2007). However, some vaccines for red pandas have not produce consistent antibody titers, thus the protective effect did not meet expectations (Loeffler et al., 2007; Zhang et al., 2017) especially in geriatric individuals. One reason for this may be a decline in immune system function associated with aging and leads to a low vaccine response and vaccine longevity as well as reducing the preventive effectiveness of vaccination (Giefing-Kroll et al., 2015; van den Berg et al., 2018). Reduced adaptive immunity leads to reduced vaccine responses and vaccine longevity (Lord, 2013). Flagellin can trigger an innate immune response by interacting with a cell surface receptor TLR5 (Uematsu et al., 2006), subsequently stimulating T cell responses and initiating adaptive immunity (Applequist et al., 2005). In both murine models and clinical trials, incorporating flagellin into vaccines has been shown to enhance vaccination efficiency in the old (Leng et al., 2011; Taylor et al., 2011). Our study has indicated that TLR5 was significantly up-regulated in geriatric red pandas, suggesting the possibility of enhancing the immune response to vaccination by incorporating flagellin, via the TLR5 innate immune pathway to address the problem of decreased vaccine response. This research provides an orientation for us to develop effective vaccines for endangered animals such as red pandas in the future. This work provides an insight into gene expression changes associated with aging, and paves the way to develop effective disease prevention and treatment strategies for red pandas in the future.