Effects of aging on gene expression in blood of captive Tibetan macaques ( Macaca thibetana) and comparisons with expression in humans.

Changes in gene expression occur as animals, including primates, age. Macaques have long been used as a model species for primate evolution and biomedical studies. Here, to study gene expression in Tibetan macaques (Macaca thibetana, TMs) and its differences to humans, we applied RNA-Seq to obtain the blood transcriptomes of 24 TMs. In total, 2 523 age-associated differentially expressed genes (DEGs) were identified. Several pathways and processes that regulate aging, including the FoxO signaling pathway, autophagy, and platelet activation, were significantly enriched in the up-regulated DEGs. Two significantly age-related modules were identified by weighted gene co-expression network analysis (WGCNA). The TMs and humans shared 279 common DEGs, including 111 up-regulated and 141 down-regulated genes with advancing age in the same expression direction. However, 27 age-related DEGs presented the opposite expression direction in TMs as that in humans. For example, INPPL1, with inhibitory effects on the B cell receptor signaling pathway, was up-regulated in humans but down-regulated in TMs. In general, our study suggests that aging is a critical factor affecting gene expression in the captive TM population. The similarities and differences in gene expression patterns between TMs and humans could provide new insights into primate evolution and benefit TM model development.

DEAR EDITOR,

Changes in gene expression occur as animals, including primates, age. Macaques have long been used as a model species for primate evolution and biomedical studies. Here, to study gene expression in Tibetan macaques (Macaca thibetana, TMs) and its differences to humans, we applied RNA-Seq to obtain the blood transcriptomes of 24 TMs. In total, 2 523 age-associated differentially expressed genes (DEGs) were identified. Several pathways and processes that regulate aging, including the FoxO signaling pathway, autophagy, and platelet activation, were significantly enriched in the up-regulated DEGs. Two significantly age-related modules were identified by weighted gene co-expression network analysis (WGCNA). The TMs and humans shared 279 common DEGs, including 111 up-regulated and 141 down-regulated genes with advancing age in the same expression direction. However, 27 age-related DEGs presented the opposite expression direction in TMs as that in humans. For example, INPPL1, with inhibitory effects on the B cell receptor signaling pathway, was up-regulated in humans but down-regulated in TMs. In general, our study suggests that aging is a critical factor affecting gene expression in the captive TM population. The similarities and differences in gene expression patterns between TMs and humans could provide new insights into primate evolution and benefit TM model development.

Aging is an essential factor affecting metabolic activities and disease resistance. Thus, related changes in gene expression are an important area of research (Charruau et al., 2016; de Magalhães & Passos, 2018; Hong et al., 2008; Horvath et al., 2012; Peters et al., 2015; Reynolds et al., 2015; van den Akker et al., 2014). TMs share many physiological and genetic similarities with humans (Fan et al., 2014; Wei et al., 2006; Wu et al., 2014), and are considered a prospective model primate for several human diseases, such as glaucoma and diabetes, and for xenotransplantation (Fan et al., 2014; Yao et al., 2013; Zheng et al., 2019). In other species, multiple aspects of advancing age have been explored, including genetic and physiological changes (Charruau et al., 2016; Eghlidi et al., 2018; Lan et al., 2020; Peters et al., 2015; Stute et al., 2012; Wu et al., 2012, 2014). However, exclusive studies on the impacts of aging on gene expression profiles in non-human primates are limited, and no such research has been conducted on TMs to date.

To investigate changes in gene expression profiles with advancing age in TMs and to clarify their differences with humans, we used RNA-Seq to obtain the blood transcriptomes of 24 TMs. A total of 2 523 age-associated DEGs were then identified. Various pathways that regulate aging were significantly enriched in the up-regulated DEGs, including the FoxO signaling pathway, autophagy, and platelet activation. Humans and TMs shared 279 common DEGs, and TMs exhibited enhanced platelet activation, which has also been investigated in human evolution (Dannemann & Kelso, 2017; Gibbons, 2017). Based on weighted gene co-expression network analysis (WGCNA), the positive age-related module contained 1 375 core genes, and the gene enrichment results were similar to those obtained from the up-regulated DEGs. Thus, in general, our study suggests that aging is a critical factor affecting gene expression in the captive TM population. The similarities and differences in gene expression patterns between TMs and humans observed here could provide new insights into primate evolution and assist in TM model development.

Peripheral whole blood was collected from 12 male and 12 female TMs randomly selected from a healthy captive population from the Jianyang Monkey breeding farm, which is affiliated with the Zoological Institution of the Provincial Hospital of Sichuan, on 10 January 2017. Relatedness distance calculated from single-nucleotide polymorphisms (SNPs) obtained from RNA-Seq analysis confirmed that the TM individuals were not closely related. Blood sampling was carried out in strict adherence to the guidelines of the Animal Care and Use Committee of Sichuan University and the Animal Ethics Committee Guidelines of the Animal Facility of the West China Hospital. The selected TMs ranged from 3 to 18 years old according to their birth date or age-related morphological characteristics. The sex ratio in each growth stage was 1:1 (Supplementary Table S1). Whole blood was drawn and temporarily preserved in specialized blood collection tubes (PAXgene Blood RNA tubes), after which the white and red blood cells were separated.

All samples were sequenced using the Illumina Hiseq 2000 system at Novogene (China). The raw data were deposited in the NCBI Sequence Read Archive (SRA) under the project accession No. PRJNA516976. A total of 1.37×109 paired-end clean reads were obtained using the NGS QC Toolkit v2.3. The saturation test showed that sufficient clean reads were obtained for each sample (Supplementary Figure S1), with these reads then mapped to the rhesus macaque genome (MMUL_8.0.1.90) with Hisat2 (Kim et al., 2015). The average mapping efficiency was 88.68%. HTSeq-count v0.9.1 (Anders et al., 2015) was applied to qualify the expression level of each gene using union mode. We filtered insufficient reads (≤10 reads) in more than half the samples, and 12 907 genes from a total of 30 807 annotated genes in the rhesus macaque genome remained, which represented a typical whole-blood expression gene-set according to previous studies (Charruau et al., 2016; Peters et al., 2015; Tung et al., 2015). DESeq2 was then applied for DEG detection using negative binomial generalized linear models and Wald statistics (Love et al., 2014). A multi-factors method was implemented in model fitting with the formula, gene expression ~ sex+RIN (RNA integrity number)+rank (RNA quality rank)+age, to estimate significance of the continuous factor of age. All genes were also analyzed using the WGCNA package in R to reveal significant modules and core genes (Horvath et al., 2012). Both DEGs and core genes were further used to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses with Kobas3 (Xie et al., 2011) and the clusterProfiler package in R. In our study, a false discovery rate (FDR) of <0.05 was set as the significance threshold.

Principal component analysis (PCA) indicated that age significantly affected the expression dataset (Figure 1A, Supplementary Figure S2). Samples were segregated by age on PC1, which explained 29% of the variance. As an independent variable, age played an important role in shaping the whole-blood transcriptome patterns in TMs (Supplementary Figure S3), similar to that reported in previous study (Peters et al., 2015). In total, 2 523 genes met our significance cut-off threshold (FDR<0.05), including 1 358 up-regulated and 1 165 down-regulated genes with advancing age (Supplementary Table S2).

Figure 1

Gene expression changes in different Tibetan macaques (TMs)

A: Assessment of effects of variables on gene expression. Distribution of each sample in top two principal components. Color of dots indicates age, with cooler tones indicating younger individuals. F: Female; M: Male. B: Relationship between TM traits and modules identified by WGCNA. C: B cell receptor signaling pathway (mcc04662). Up-regulated DEGs with age are marked in purple, down-regulated DEGs with age are marked in light blue. Activation process is shown with red arrow, indirect activation process is shown with dotted arrow. D: FoxO signaling pathway (mcc04068). Expression of up-regulated DEGs with age in TMs are marked in purple (upstream of main pathway) and light brown (downstream of main pathway), down-regulated DEGs with age are marked in light blue. Activation process is shown with red arrow, indirect activation process is shown with dotted arrow.

Several previously identified age-related genes were also identified as age-related DEGs in this study, such as IGF-1 (P=0.016, FDR=0.069), IGF2R (P=0.001, FDR=0.012) (Hoopes et al., 2012), FOXO3 (P<0.001, FDR=0.006), andFOXO4 (P=0.004, FDR=0.03) (Göring et al., 2007) (Supplementary Table S2). Age-related DEGs identified in human blood (Peters et al., 2015) were compared with DEGs in our study. Among the 279 shared DEGs, 111 were up-regulated, 141 were down-regulated, and 27 showed different expression directions in TMs and humans with advancing age (Table 1).

Table 1

Age-related DEGs in Tibetan macaques with opposite expression to humans

Gene name	log2FC	P-value	Gene description
NET1	0.100402	2.30E–05	Neuroepithelial cell transforming 1
MAML2	0.064429	0.000227	Mastermind like transcriptional coactivator 2
TBC1D4	0.053905	0.000237	TBC1 domain family, member 4
LPGAT1	0.073952	0.000393	Lysophosphatidylglycerol acyltransferase 1
CEP135	0.083235	0.00043	Centrosomal protein 135
TGFBR2	0.050934	0.000482	TGF-beta receptor type-2 precursor
LARP4	0.046432	0.001619	La ribonucleoprotein domain family member 4
NIN	0.059529	0.001998	Ninein
RORC	0.065932	0.002403	RAR related orphan receptor C
RNF144A	0.037814	0.005089	Ring finger protein 144A
PGRMC2	0.023864	0.006044	Progesterone receptor membrane component 2
TIGD7	0.040531	0.006794	Tigger transposable element derived 7
CRLF3	0.030661	0.006891	Cytokine receptor like factor 3
FNIP1	0.044619	0.007115	FNIP1
PDE3B	0.034851	0.007153	Phosphodiesterase 3B
OSBPL5	–0.10469	1.58E–08	Oxysterol binding protein like 5
CDKN1C	–0.15178	5.59E–05	Cyclin dependent kinase inhibitor 1C
IL18BP	–0.04766	0.00071	Interleukin 18 binding protein
PALLD	–0.09838	0.000947	Palladin, cytoskeletal associated protein
EPB41L4A	–0.08794	0.001517	Erythrocyte membrane protein band 4.1 like 4A
INPPL1	–0.03325	0.002177	Inositol polyphosphate phosphatase like 1
IL12RB1	–0.05667	0.002345	Interleukin 12 receptor subunit beta 1
TIMP1	–0.11489	0.003691	Metalloproteinase inhibitor 1
GNLY	–0.07513	0.004472	Granulysin
FAM46A	–0.02272	0.005602	Family with sequence similarity 46 member A
TAGLN	–0.07185	0.008341	Transgelin
AGPAT4	–0.04314	0.00835	1-acylglycerol-3-phosphate O-acyltransferase 4

Based on WGCNA, we explored the correlations among gene modules and all possible variables, with particular interest in those related to age. We identified 33 modules, including two gene modules (blue and royal-blue) significantly related to age (Figure 1B). The positive age-related blue module (R2=0.64, P=0.0008) contained 1 375 core genes, 906 of which were age-related DEGs. These core genes were also significantly enriched in the FoxO signaling pathway, B cell receptor signaling pathway, platelet activation, and innate immunity (Supplementary Table S3). Multiple genes showed very strong connections (weight>0.42) to FOXO3 and FOXO4 by WGCNA (Supplementary Figure S4). The royal-blue module contained 384 genes that were negatively associated (R2=–0.58, P=0.003) with advancing age, including 144 DEGs related to age. The genes in the royal-blue module were mainly related to ribosome (KEGG:03010) and spliceosome (KEGG:03040) (Supplementary Table S4).

Figure 1

In the TMs, the 1 165 down-regulated DEGs with advancing age were enriched in basal amide and peptide metabolism from rRNA processing and ribosome biogenesis to the translation process (Supplementary Tables S5, S6). Suppressing translation is an essential part of the mechanism of aging (López-Otín et al., 2013; Peters et al., 2015). For cell components, the down-regulated DEGs were mainly located in mitochondrial parts, such as the mitochondrial envelope and membrane. Two biological processes (rRNA assembly and protein translation) were impacted by age in the TMs, with similar results also confirmed in human lymphocytes (Hong et al., 2008). Furthermore, several ribosomal proteins (RNPs) were also significantly down-regulated with age based on WGCNA, i.e., small ribosomal protein L12 (RPL12) and large subunit ribosomal protein L37A (RPL37A). These RNPs are essential components of the ribosome, and their inhibition can decrease ribosome biogenesis. Ribosome biogenesis and translation are among the most energy-consuming cellular processes and may reflect energy metabolism changes in cells (Granneman & Tollervey, 2007; Lindqvist et al., 2018). Although it is still difficult to resolve the initial step of aging, the rRNA-protein-DNA downward spiral may be an important link in aging processes in TMs, as found in other organisms (Reynolds et al., 2015; Yin et al., 2017).

Based on our results at the gene expression level, the B cell receptor signaling pathway showed increased expression with age (Figure 1C; Supplementary Tables S7, S8). The B cell receptor signaling pathway is supposed to improve with advancing age (Simon et al., 2015), and we found 22 age-related DEGs enriched in this pathway, including 13 up-regulated DEGs and nine down-regulated DEGs. The greater number of up-regulated than down-regulated genes suggests that the up-regulated DEGs with advancing age may play crucial roles in the B cell receptor signaling pathway, and many genes function near the nucleus through the reaction pathway, like AP1 (Figure 1C).

The key KO (KEGG Orthology) in the FoxO signaling pathway is FOXO. Here, two age-related DEGs (FOXO3 and FOXO4) were identified in our study (Figure 1D). FOXO3 and FOXO4 are regulated by several energy metabolism genes, including IGF1R and AMPK. In subsequent processes, the two genes also regulate the expression of a series of genes involved in aging-related functions such as cell cycle regulation, apoptosis, autophagy, oxidative stress resistance and DNA repair, glycolysis/gluconeogenesis metabolism, and immune-regulation (de Magalhães & Passos, 2018; Kenyon, 2010). Furthermore, the FoxO signaling pathway was also identified as a positive age-related pathway based on WGCNA (Supplementary Figure S4 and Table S3). The up-regulation of FOXO3 and FOXO4 may enhance these functions and help organisms in resistance to aging (van der Horst & Burgering, 2007).

Several studies have investigated the relationships between age and gene expression in blood, and the results have been commonly used in clinical disease detection (Aziz et al., 2007; Lapointe et al., 2004; Mesko et al., 2011). These studies have emphasized the tendency that age-related DEGs are mainly enriched in metabolic level aspects, consistent with our study. Hence, we argue that metabolic responses are also essential in dealing with changes during aging in TMs.

TMs are important animal models in biomedical research due to their close genetic relatedness to humans, sharing a common ancestor ~25 million years ago (Mya) (Fan et al., 2014; Yao et al., 2013; Zheng et al., 2019). Therefore, it is important to compare gene expression patterns with advancing age between TMs and humans. Here, we found that TMs and humans shared 279 age-related DEGs, as well as many common enriched functions. In brief, 90.32% (252) of age-related DEGs were regulated in the same direction between humans and TMs. For example, the conserved nucleophosmin/nucleoplasmin 3 (NPM3) gene was down-regulated with advancing age. The nucleophosmin/nucleoplasmin protein participates in various significant cellular activities like sperm chromatin remodeling, nucleosome assembly, genome stability, ribosome biogenesis, DNA duplication, and transcriptional regulation (Frehlick et al., 2007; Shu & Zhang, 2007). Thus, NPM3 expression could be potentially used as a molecular marker for age evaluation.

However, TMs and humans also had genes showing different expression directions with advancing age. Among the 279 commonly shared age-related DEGs, 27 showed the opposite expression direction in the TM and human blood transcriptomes (Table 1). These 27 genes were enriched in biological process terms, such as wound healing involved in the inflammatory response (GO:0002246), response to stimulus (GO:0050896), and regulation of cellular process (GO: 0050794), suggesting possible physiological differences in responding to disease or environmental change in aging processes. For example, INPPL1, which functions in the B cell receptor signaling pathway, could affect B cells and subsequent humoral immunity (Figure 1C). INPPL1 with inhibitory effects was up-regulated in humans but down-regulated in TMs.

Another significant enhancement in older TMs was platelet activation, which is associated with coagulation and wound healing. Blood coagulation is considered a risk factor for arteriosclerosis (Favaloro et al., 2014; Mari et al., 2008), and increases during human aging. Why humans retain this molecular mechanism with aging remains unknown. We detected a similar enhancement of the coagulation function in the TMs. Of the 19 up-regulated DEGs enriched in coagulation in humans, three (F2R, MAPK1, PIK3CG) were shared with the TMs. We assume that this mechanism may be a complement to the TM immune response and helpful in accelerating the healing of infected wounds. Coagulation also played an important role in human evolution, as shown in studies of the Neanderthal fossil genome (Dannemann & Kelso, 2017; Gibbons, 2017). However, the less common formation of thrombosis in TM individuals suggests differences in the coagulation and anticoagulation mechanisms between TMs and humans (Lyman et al., 2018; Sheffield et al., 1981). This finding may offer new insight into medicinal development for arteriosclerosis.

In conclusion, we investigated gene expression profiles in TMs through RNA-Seq and identified age-related genes and their functional pathways. In general, our study suggests that aging is a critical factor affecting gene expression in the TM population. The identified age-related genes should provide a rich trove of data for future aging studies. The similarities and differences in gene expression profiles between TMs and humans are also discussed in our study, which could provide new insights into primate evolution and benefit TM model development.

DATA AVAILABILITY

The raw data have been deposited in the NCBI Sequence Read Archive (SRA) database under SRA project accession No. PRJNA516976.

Supplementary data to this article can be found online.

Click here for additional data file.

COMPETING INTERESTS

The authors declare that they have no competing interests.

AUTHORS’ CONTRIBUTIONS

C.Y., X.Z., L.Z., Q.Y., and M.Z. performed the bioinformatics analyses; L.Z., Z.Y., and L.Z. collected the samples; J.X., M.P., and J.L. revised the manuscript; C.Y., Q.Y., and Z.F. wrote the manuscript; B.S. and Z.X.F. designed and supervised the study. All authors read and approved the final version of the manuscript.