DNA methylation and transcriptome analysis reveal epigenomic differences among three macaque species

Abstract Macaques (genus Macaca) are the most widely distributed non‐human primates, and their evolutionary history, gene expression profiles, and genetic differences have been extensively studied. However, the DNA methylomes of macaque species are not available in public databases, which hampers understanding of epigenetic differences among macaque species. Epigenetic modifications can potentially affect development, physiology, behavior, and evolution. Here, we investigated the methylation patterns of the Tibetan macaque (M. thibetana; TM), Chinese rhesus macaque (M. mulatta lasiota; CR), and crab‐eating macaque (M. fascicularis; CE) through whole‐genome bisulfite sequencing from peripheral blood. We compared genome‐wide methylation site information for the three species. We identified 12,128 (CR vs. CE), 59,165 (CR vs. TM), and 39,751 (CE vs. TM) differentially methylated regions (DMRs) in the three macaques. Furthermore, we obtained the differentially expressed genes (DEGs) among the three macaque species. The differences between CR and CE were smaller at both the methylome and transcriptome levels than compared with TM (CR vs. TM and CE vs. TM). We also found a change in the density of single nucleotide mutations in DMRs relative to their flanking regions, indicating a potential mechanism through which genomic alterations may modulate methylation landscapes, thereby influencing the transcriptome. Functional enrichment analyses showed the DMR‐related genes were enriched in developmental processes and neurological functions, such as the growth hormone‐related pathway, insulin secretion pathway, thyroid hormone synthesis pathway, morphine addiction, and GABAergic synapses. These differences may be associated with variations in physiology and habitat among the macaques. Our study provides one of the first genome‐wide comparisons of genetic, gene expression, and epigenetic variations across different macaques. Our results should facilitate further research on comparative genomic and genetic differences in macaque species.


| INTRODUC TI ON
DNA methylation is the main mode of epigenetic modification, which can alter genetic expression without altering the DNA sequence and can silence genes, making them ineffective.As an essential gene regulatory process, DNA methylation facilitates phenotypic variability in response to environmental change, independent of genetic variation (Angers et al., 2010;Jeremias et al., 2018).Furthermore, DNA methylation plays a significant role in organismal development, gene expression profiles, genomic imprinting, disease susceptibility, and evolutionary dynamics (Mendizabal et al., 2016).Epigenetic mechanisms can easily be integrated into a model of phenotypic variation in multicellular organisms, offering insights into the phenotypic differences observed among genetically identical organisms (Wong et al., 2005).
While genetic disparities alone cannot fully account for observed phenotypic differences, research has demonstrated that DNA methylation can contribute to such variations.A classic example occurs in queen and worker bees, which share the same genome but exhibit different phenotypes.Notable: Lyko et al. (2010) reported significant differences in methylation levels in 500 genes between queen and worker bees, while Kucharski et al. (2008) used RNA interference to silence the DNMT3 gene, resulting in the development of young bees into queens.DNA methylation is also associated with phenotypic differences in human populations (Heyn et al., 2013), especially among twins (Fraga et al., 2005;Kaminsky et al., 2009).Hence, DNA methylation serves as a mechanism capable of elucidating phenotypic disparities and providing a theoretical foundation for such observations.Macaques, classified as endangered animals in their natural habitat, serve as important model animals in medical-related research (Gibbs et al., 2007;Xu et al., 2022).Despite their similar genetic backgrounds, phenotypic differences exist among macaque species (Gibbs et al., 2007).Tibetan macaques (Macaca thibetana), which are endemic to China, differ substantially from crab-eating macaques (M.fascicularis) in geographical range, body size, morphology, physiology, and behavior (Fan et al., 2014).Furthermore, compared with Chinese rhesus macaques (M.mulatta lasiota), TM exhibits a larger body size, a longer life span, a calmer temperament, and greater trainability (Solari & Baker, 2007).Previous analyses have identified the pattern of DNA differences among different macaques (Fan et al., 2014;Fan et al., 2018).For example, Fan et al. (2014) reported that Tibetan macaques possess more nonsynonymous variants in genes related to immune response and glucose metabolism than rhesus and crab-eating macaques.Studies have demonstrated that DNA methylation levels differ among species (Klughammer et al., 2023).However, no studies have focused on the DNA methylation levels between different macaques and their correlation with gene expression levels.
Thus, the small number of fixed genetic divergences at the DNA level may be insufficient for explaining all phenotypic differences.Consequently, there is an urgent need for further DNA methylation research.
To date, investigations into methylation patterns have included various primate species (Hernando-Herraez, Garcia-Perez, et al., 2015;Hodges et al., 2011;Molaro et al., 2011).However, the epigenomic differences among macaque species remain unclear.In this study, we used the latest high-quality rhesus macaque reference genome (M.mulatta, Mmul_10) (Warren Wesley et al., 2020) and performed whole-genome bisulfite sequencing (WGBS) on multiple M. mulatta lasiota, M. fascicularis, and M. thibetana individuals (with relationship (M.thibetana (M.mulatta lasiota, M. fascicularis))).We aimed to explore their epigenomic differences and investigate the relationship between epigenomic regulation and gene expression, as well as the association between genomic and epigenomic differences.This study should enhance our understanding of comparative genomics and genetic discrepancies within macaque species, thereby deepening our understanding of macaque health, disease, development, and more.

| Genome-wide DNA methylation landscape of three macaque species
A total of 21 transcriptomes and 15 methylomes (from two different groups of macaques; see Table S1) were sequenced or downloaded from previous studies on three macaque species, including Tibetan macaques (M.thibetana; TM), Chinese rhesus macaques (M.mulatta lasiota; CR), and crab-eating macaques (M.fascicularis; CE).The bisulfite conversion rate of all samples was greater than 99% (Table S2).All methylome data were mapped to the M. mulatta genome (Mmul_10), with a mean efficiency of 75.75% (Table S2).
After using deduplicate_bismark, the remaining 86.65%-97.60% of reads were used for subsequent analysis (Table S3).At the genome-wide level, 2.92%-4.24% of cytosine sites were found to be methylated.Of these methylated sites, 82.45%-90.52% were in the mCpG context, 2.25%-4.44% in the mCHG context, and 7.23%-13.11% in the mCHH context.The overall DNA methylation patterns showed comparable genome-wide methylation ratios among the different species, with 68.80%-75.50% of cytosine sites in the CpG context being methylated and less than 1% of cytosine sites in the CHG and CHH contexts being methylated (Figure S1a; Table S3).Methylation level distribution analysis demonstrated that the methylation levels at the CpG site were relatively high, with relatively low methylation levels at the CHH and CHG sites (Figure S1b).
We delineated trends in the methylation levels of gene bodies and flanking regions (Figure S1c).Notably, the region proximal to the transcription start site (dotted line on the left) showed significant fluctuations in methylation levels, gradually decreasing toward the transcription start site.In particular, CR showed a significant reduction in methylation levels near the transcription start site.Circos plots of mean methylation levels (Figure S2) showed both hypermethylated and hypomethylated regions in each chromosome of the | 3 of 14 WANG et al.
three species.In addition, centromeres and their flanking regions exhibited relatively low levels of methylation.
To examine the factors shaping the distribution patterns of CpG methylation, we conducted principal component analysis (PCA) using the methylation data from the samples.As shown in Figure 1a, the two main components affecting methylation distribution accounted for 18.9% and 11% of the variance, respectively.The samples were also clustered by species.Thus, to clarify the methylation relationship among the three species, we performed cluster analysis of the CpG information of the samples (Figure 1b).The results were aligned with the PCA findings, confirming that samples from the same species formed distinct clusters, with a particularly close correlation between CR and CE.hypermethylated DMRs n = 30,338; hypomethylated DMRs n = 9413).In each group, hypermethylation and hypomethylation refer to the former species relative to the latter.The CR versus CE group exhibited a higher number of hypomethylated DMRs, whereas the CR versus TM and CE versus TM groups had a greater number of hypermethylated DMRs.However, for most DMRs, the mean methylation level was higher in CR than in CE and TM, and higher in CE than in TM (Figure 2a-f).We also excluded the effect of age-and sex-related DMR on inter-species DMR and found minimal overlap between age-and sex-related DMR and inter-species DMR, including only 0.1513% (168/111044) age-and 0.3899% (433/111044) sex-related common DMRs with inter-species DMRs (see Materials and Methods, Tables S4-S8).

| Differentially methylated regions (DMRs) between different macaque species
To further analyze the distribution of methylation differences across the genome, we visualized the lengths of the DMRs (Figure S3).Interestingly, the DMRs in the three groups exhib- However, the highest DMR distribution density in the CE versus TM group was on chromosome 4 (0.4122%).The lowest DMR density in the CR versus CE and CR versus TM groups was on chromosome 13 (CR vs. CE = 0.0877%; CR vs. TM = 0.5533%).
Interestingly, the lowest DMR density distribution in the CE versus TM group was on chromosome 19 (0.3082%).We also plotted the density distribution of DMRs and corresponding regions at the chromosomal level (Figure 3a).First, the distribution of DMRs in each chromosome was not uniform but fluctuated, with some F I G U R E 1 Display of overall methylation pattern and cluster analysis of CpG locus information.(a) PCA of CpG site methylation information.(b) Cluster analysis performed based on the methylation information of CpG sites.Distance method: "correlation"; Clustering method: "ward".chromosomal regions showing high DMR density.For example, the start of chromosomes 14 and 19 exhibited a higher DMR density, suggesting that differential methylation between species may prefer to occur in specific regions.
Next, we annotated functional regions of DMRs according to annotation files and counted the number of different genomic elements that overlapped with the DMRs (Figure 3b-d).Among the genomic elements, the proportion of introns was highest in all three groups: CR versus CE n = 11,244 (37.7%),CR versus TM groups n = 27,001 (40.17%), and CE versus TM n = 19,215 (58.65%).The CE versus TM group contained fewer overlapping promoters with DMRs (2342, 7.15%) compared to the other two groups.In addition, the proportion of CpG islands and flanking-related areas covered by DMRs was small (0.39%-1.3%).Analysis of the coverage of repeating elements and DMRs indicated a high proportion of coverage over long and short-interspersed nuclear elements (LINEs and SINEs) (Figure S6).
Finally, we obtained the DMR-related genes (DMGS, overlapping with gene body) and DMR-related promoters (DMPs, overlapping with promoters) in the three groups and found more common DMGs (2947) than common DMPs (244) among the groups (Figure S7).To better understand the biological functions of DMGs and DMPs, we conducted Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses.Based on GO enrichment analysis of DMGs (Figure S8a-c; Table S9), negative regulation of cell communication (GO:0010648), cell projection morphogenesis (GO:0048858), and peptidyl-amino acid modification (GO:0018193) were enriched in the CR versus CE group, cellular macromolecule catabolic processes (GO:0044265) and cellular protein catabolic processes (GO:0044257) were enriched in the CR versus TM group, and small GTPase-mediated signal transduction (GO:0007264) was enriched in the CE versus TM group.Furthermore, overlapping DMGs (common in three macaques) were enriched in categories related to neuro development, including neuron development (GO:0048666), axon development (GO:0061564), and dendrite development (GO:0016358), and categories related to cell development, including cell morphogenesis (GO:0000902), cell projection morphogenesis (GO:0048858), and cellular component morphogenesis (GO:0032989).In the KEGG enrichment analysis, we plotted the 20 most significant pathways in each group (Figure S9a-c; Table S10).Several enriched KEGG pathways were related to essential hormones in organisms, including, but not limited to, growth hormone synthesis, secretion, and action (KEGG:04935), the insulin signaling pathway (KEGG:04910), and the thyroid hormone signaling pathway (KEGG:04919).We selected the SOCS1 gene  Analysis of DMPs showed that only the CR vs CE and CR vs TM groups were enriched in GO categories and KEGG pathways.
Based on GO enrichment analysis, both groups were enriched in actin filament-based processes (GO:0030029), the Wnt signaling pathway (GO:0016055), and actin cytoskeleton organization (GO:0030036) (Figure S10a-b; Table S11).Based on KEGG pathway  S12).DMPs in the CE vs TM group and shared among all three groups showed no significant enrichment.

| Gene expression analysis
Several studies have investigated gene expression differences at different physiological stages within specific Macaca species (Lan et al., 2020;Yan et al., 2020).However, research on gene expression differences between macaque species is scarce.Given the potential regulatory role of DNA methylation in gene expression, we examined the gene expression profiles of the three macaque species.Cluster analysis showed that samples were primarily clustered according to species (Figure 4a).Similarly, PCA demonstrated that the samples were mainly distributed according to species, with PC1 accounting for 47% of the variance, representing a major factor in species variation (Figure 4b).The distance between CR and CE was closer than that between CR and TM (Figure 4b).We next analyzed and identi-  S13).
The quantity distribution exhibited a pattern similar to that observed for the DMGs, that is CR vs TM > CE versus TM > CR versus CE.
To elucidate the biological functions of the DEGs, we conducted GO and KEGG enrichment analyses (Figure S13a-c; Table S14).For  S15).
Among the three groups, we identified 339 shared DEGs (Figure S12).Based on further GO and KEGG enrichment analysis (Tables S14 and S15), these DEGs were enriched in inclusion body (GO:0016234) and porphyrin metabolism (KEGG:00860).Porphyrin metabolism is recognized as a vital process involved in heme synthesis (Bonkovsky et al., 2013).

| Combined methylation and transcriptome analysis
To further explore the methylation of genes with different expression intensities near transcription start sites, we grouped the genes according to transcripts per million (TPM) intensity (TPM = 0, 0 < TPM ≤1, 1 < TPM ≤10, and TPM ≥10) to obtain methylation trends.From the perspective of gene expression intensity, methylation levels in the gene body and flanking regions changed with increasing TPM (Figure 5a-d; Figure S15a-h).In general, gene expression intensity increased with the decreasing methylation level in the 2-kb region upstream of the gene body.In contrast, gene expression intensity increased slowly with methylation level of the gene body.Regions near the transcription start sites were associated with gene expression intensity and methylation levels (Jones, 2012).The methylation level in the 2-kb region downstream of the gene body fluctuated around 0.7 with no apparent regularity.In the gene body region, we compared differences in methylation levels between different species under the same expression intensity (Figure 5a-d; Figure S15a-h).Results showed that the methylation level near the transcription start site decreased with increasing gene expression intensity, with the greatest decrease in CR, followed by CE and TM.
These results suggest potential differences between species in how DNA methylation regulates gene expression.
To explore the association between the methylome and transcriptome, we conducted multi-omics and enrichment analyses.
We defined DEGs with DMRs on the promoter or gene body as DMR-related DEGs (Table S16).According to analysis, the CR versus CE, CR versus TM, and CE versus TM groups contained 545, 3314, and 1143 DMR-related DEGs, respectively (Figure 5e-g).
Enrichment analysis was performed as above (Figure S16a-c; Tables S17 and S18).Due to the limited number of genes, the DMR-related DEGs in the CR versus CE group showed no GO and KEGG enrichment.However, the DMR-related DEGs in the CR vs TM group were enriched in several GO:MF categories, including protein kinase activity (GO:0004672) (Figure S16a) and KEGG pathways related to phenotypic differences, including growth hormone synthesis, secretion, and action (KEGG:04935) (Figure 5h).These results suggest that the phenotypic differences corresponding to KEGG pathways may be related to variation in gene expression levels caused by methylation differences.The DMR-related DEGs in the CE versus TM group were enriched in several GO:CC categories (Figure S16b), and KEGG pathways (e.g., platelet activation (KEGG:04611)) (Figure S16c) also enriched in the CR versus TM group (Figure 5h).

| Genomic variation of DMRs
Previous research has indicated that human hypermethylated DMRs exhibit more base mutations than other nearby regions (Hernando-Herraez, Heyn, et al., 2015).To confirm the association between DMRs and base mutations among macaque species, we calculated the base mutation density within DMRs in different groups.Results revealed that hypermethylated DMRs displayed a higher base mutation density than their flanking regions, with methylation changes and an accelerated base mutation rate (Figure 5i-k).Conversely, the analysis showed a decrease in base mutation density in the CR versus CE and CR versus TM hypomethylated DMRs relative to the lateral region (Figure S17a,b), but an increase in base mutation density in the CE versus TM hypomethylated DMRs (Figure S17c).These findings suggest potential genetic conservation of base mutations in these regions.Moreover, both hypermethylated and hypomethylated DMRs displayed distinct base mutation patterns in the flanking sequence.

| DISCUSS ION
Despite the close relationship among macaque species, variations in their phenotypes, such as body size, exist.Whole-genome sequencing studies have shed light on the evolutionary history and genetic distinctions among macaque species (Fan et al., 2014;Higashino et al., 2012;Osada et al., 2015;Yan et al., 2011).In our prior genomic study, we observed similarities in the genetic backgrounds of macaque species but also identified significant differences in specific genes, such as immune response genes and insulin-related genes, among the different species (Fan et al., 2014).
However, genetic variations at the DNA level alone cannot fully account for the observed variations in phenotypes among macaque species.Previous epigenetic studies on macaques have been limited, focusing on either a fraction of the genome or only on one species (e.g., rhesus macaques) without integrating genome and transcriptome data (Hata et al., 2013;Howard et al., 2011;Massart et al., 2014;Provençal et al., 2012).Here, we employed WGBS to obtain single-base resolution methylation profiles of three macaque species and applied multi-omics analysis of related genes affecting phenotypic differences.This study provides one of the first genome-wide comparisons of genetic, gene expression, and epigenetic variation among different macaque species.

Despite the close genotypic relationship between CR and TM,
which diverged approximately three million years ago, they also exhibit significant phenotypic differences (Fan et al., 2014;Perelman et al., 2011).Various studies have investigated distinctions among macaque species from different perspectives.For instance, Li et al. (2020) used copy number variation (CNV) to construct a neighbor-joining tree of nine Macaca species.Song et al. (2021) revealed that complex introgressive gene flow played an important role in the evolutionary history of Macaca and suggested that macaque evolution can be more accurately understood as phylogenetic networks rather than trees.Cui et al. (2019) reported similarities in the gut microbiota of CR and CE.They also found significantly higher alpha diversity in the rectal microbiota of CE compared to CR, despite a similar number of operational taxonomic units (OTUs), as well as significantly different relative abundances of dominant taxa, despite similar taxa at both the phylum and genus levels (Cui et al., 2019).

Fan et al. (2018) used whole-genome sequences of several
Macaca species to generate and analyze single nucleotide variant (SNV) data and found that the genetic distance between TM and CE was closer than that between TM and CR, with CR and CE exhibiting the closest genetic relationship.Our study revealed a relationship between DNA methylation patterns, differential gene expression, and species-relatedness.Specifically, the greater the genetic distance between species, the more pronounced the differences in methylation and gene expression profiles.Similarly, we found that the methylation difference between CR and CE was the smallest, while that between CR and TM was the largest.
Moreover, the number of DEGs in the CR versus CE group was the lowest, whereas that in the CR versus TM group was the highest.
In addition, the distribution of DMR density in the CR versus CE and CR versus TM groups was more similar, indicating that methylation differences between CR and other species may occur at specific locations.These findings provide evidence of the correlation of the genome with the methylation profile, and the methylation profiles with the transcriptome.

Species have developed distinct adaptation patterns spanning
from the genome to methylation groups and transcriptomes to adapt to different environments (Wong et al., 2005).The three macaque species studied here exhibit variations in their habitats and behaviors (Fleagle, 2013), and their adaptability to the environment could not solely be determined by genetics but also influenced by DNA methylation.Our functional enrichment analyses revealed that the DMR-related genes (DMGs and DMPs) were enriched in developmental processes and several signaling pathways, including insulin secretion, thyroid hormone, morphine addiction, and the GABAergic synapse (Tables S10 and S12).Notably, the Wnt signaling pathway was enriched in multiple analyses.The Wnt signaling pathway plays a vital role in early embryonic development, organogenesis, tissue regeneration, and other physiological processes (Arkell et al., 2013).These patterns suggest that the DMRs are located in genomic regions involved in brain development and other traits, such as energy metabolism, which are particularly relevant to recent evolutionary history, as epigenetic involvement in human and primate brain evolution has been demonstrated (Mendizabal et al., 2016).

| DMRs in metabolism-related pathways
Considerable phenotypic differences, including variations in body size, exist among macaque species.TM exhibits the largest body size within the genus Macaca, while CE macaques are notably smaller than CR macaques.Macaque body size is influenced by several factors, such as social status and foraging behavior at an ethological level (Saito, 1996;Shively et al., 1997;Shively & Wallace, 2001), as well as hypothalamic-pituitary-adrenal axis function at an endocrine level (Kaufman et al., 2007).The growth and development of organisms are influenced by the endocrine signaling system in conjunction with signal transduction pathways (Neave, 2008).However, the specific impact of different levels of methylation on body size in the three macaque species has not yet been revealed.In our analysis of differential methylation among macaques, we identified several DMGs and DMPs that were enriched in pathways related to body size, growth, and development, including growth hormone synthesis, secretion, and action (KEGG:04935), insulin signaling pathway (KEGG:04910), thyroid hormone signaling pathway (KEGG:04919), and thyroid hormone synthesis (KEGG:04918) (Tables S10 and S12).
Notably, the growth hormone synthesis, secretion, and action pathway (KEGG:04935) was enriched in all DMG groups (Table S10).Among the identified DMGs, SOCS1 is involved in the cell growth and cytoskeletal recombination stages of the growth hormone synthesis, secretion, and action pathway (Keren et al., 2006;Krebs & Hilton, 2001).Overexpression of SOCS1 can completely inactivate the growth hormone (GH) signaling pathway (Ram & Waxman, 1999).Furthermore, SOCS1 can inhibit GH receptor-dependent JAK2 tyrosine phosphorylation and decrease GH activity (Ram & Waxman, 1999).Methylation of SOCS1 can eliminate the inhibition of the JAK/STAT pathway at the methylation level, whereas demethylation of SOCS1 can suppress cell growth and promote anchorage-independent growth by inducing apoptosis (Yoshikawa et al., 2001).Therefore, we hypothesize that the up-regulation of SOCS1 gene expression in CR and CE compared to TM may be associated with the larger TM body size.
DNA methylation can affect gene transcription near the transcription start site or in the gene body region.Therefore, the downregulation of SOCS1 gene expression in TM may be attributed to hypermethylation.As a key downstream gene of the growth hormone pathway, SOCS1 plays a crucial role in the growth and development of organisms.Hence, the down-regulation of SOCS1 gene expression in TM may be associated with its relatively large body size.
Blood glucose level is an essential factor affecting animal growth and development.Glucose-induced insulin secretion and its potentiation constitute the principal mechanism of insulin release, while glucose metabolism generates adenosine triphosphate (ATP) to support daily activities (Seino et al., 2010).
Studies have reported significant differences in blood glucose between CR and TM, with TM exhibiting relatively higher blood sugar levels (Yang et al., 2010).Fan et al. (2014) reported that Tibetan macaques contain more nonsynonymous mutations in genes involved in glucose metabolism, insulin secretion, and insulin receptor signaling pathways than rhesus and crab-eating macaques.However, genomic differences alone may not fully explain the variations in genes related to insulin pathways, necessitating the examination of epigenetic differences.Our analysis revealed enrichment in the insulin signaling pathway (KEGG:04910) and insulin secretion (KEGG:04911) in all comparison groups (Tables S10 and S12), as well as enrichment in insulin receptors (INSR; ENSMMUG00000028907), which are essential components of insulin-related pathways (Tables S10 and S12).In mammals, changes in insulin receptor levels are of physiological importance (Payankaulam et al., 2019).In animals, the physiological function of INSR is closely related to tyrosine kinase signal transduction (Ward et al., 2007).Heterozygous INSR mutant mice generally exhibit normal growth and fertility but show deficiencies in downstream signaling (Batista et al., 2019).In humans, heterozygous carriers of null mutations of INSR exhibit abnormal glucose tolerance, highlighting the importance of INSR gene expression (Psiachou et al., 1993).We identified INSR as a DMR-related DEG in the CR versus CE and CR versus TM groups, with CR showing lower INSR gene expression compared to CE and TM.Therefore, both genomic and epigenetic differences were detected in identical/similar pathways, suggesting that these pathways may play crucial roles in the variation of blood glucose levels among macaque species.
Our results further showed the enrichment of DMGs in the thyroid hormone signaling pathway (KEGG:04919) in all comparison groups and thyroid hormone synthesis (KEGG:04918) in CR versus CE, CR versus TM, and the intersection of DMG groups (Tables S10and   S12), which are essential in the regulation of growth, development, and metabolism (Cheng et al., 2010).KEGG analysis identified 58 (CR vs. CE), 95 (CR vs. TM), and 59 (CE vs. TM) DMGs enriched in the thyroid hormone signaling pathway (KEGG:04919).Interestingly, seven DMGs in the three groups (ADCY1, ADCY2, ADCY3, ADCY5, ADCY7, ADCY8, and ADCY9) were enriched in the same enzyme, adenylate cyclase (AC: EC 4.6.1.1),with ADCY5 and ADCY6 also identified as DMR-related DEGs in the CR versus TM group.The enzyme protein kinase C (PKC: EC 2.7.11.13), related to PRKCA (ENSMMUG00000018378) and PRKCB (ENSMMUG00000014271), was found in all DMG groups.Furthermore, PRKCA was identified as a DMR-related DEG in the CR versus TM group.We found that a certain number of thyroid hormone pathway-related genes were differentially methylated, with some also showing differential expression.
As the thyroid hormone-related pathways regulate physiological functions, we consider that differential methylation may lead to corresponding physiological differences in the three macaque species.

| DMRs in the neurosystem-related pathway
In total, 34 (CR vs. CE), 69 (CR vs. TM), and 51 (CE vs. TM) DMGs were enriched in the GABAergic synapse pathway (KEGG:04727) (Tables S10 and S12).Gamma-aminobutyric acid (GABA) is the most abundant inhibitory neurotransmitter in the mammalian central nervous system and binds to three major classes of receptors, that is, GABAA, GABAB, and GABAC (Chalifoux & Carter, 2011;Pinard et al., 2010).GABAA and GABAC receptors are ionotropic and mediate rapid GABA responses by triggering chloride channel opening, while GABAB receptors are metabotropic and mediate slower GABA responses by activating G-proteins and influencing second messenger systems.We found that DMGs enriched in the GABAergic synapse pathway were related to protein subunits in the GABAA and GABAB receptors.Furthermore, GABBR2 was identified as a DMG in the CR versus CE, CR vs TM, and CE versus TM groups.At the transcriptome level, GABBR2 gene expression was down-regulated in CR and CE compared with TM.In addition, the GABAA receptorrelated GABRA2 gene (ENSMMUG00000003614) was identified as a DMR-related DEG in the CR versus TM group.
The morphine addiction pathway (KEGG:05032) is closely associated with the GABAergic synapse pathway as morphine can indirectly excite ventral tegmental area dopamine neurons by reducing inhibitory synaptic transmission mediated by GABAergic neurons (Yan et al., 2008).Here, we identified 44 (CR vs. CE), 70 (CR vs. TM), and 64 (CE vs. TM) DMGs enriched in this pathway.
Our findings revealed differential methylation and gene expression patterns in neurosystem-related genes among the three Macaca species, suggesting potential neurological differences among them.
The observed variations in methylation and gene expression in neurosystem-related genes may contribute to the differences in body size observed among the Macaca species due to endocrine regulation.

| CON CLUS IONS
This study represents the first investigation into the genome-wide DNA methylation patterns of rhesus, Tibetan, and crab-eating macaques, providing valuable insights into their epigenetics and evolutionary biology.Furthermore, the transcriptomic analysis allowed us to investigate the relationship between differential methylation and gene expression.The enrichment of DMRs in pathways related to the nervous system and hormone regulation suggests an influence on energy metabolism and daily activities, as well as differentiation among species at the epigenetic level, which may partially explain the differences in body size and evolutionary history among macaques.However, our study only focuses on the differences between species, and the impact of other important factors, such as sex and age, on epigenetic inheritance also needs more in-depth research.Besides, it is important to note that our study was limited to whole blood, and future research encompassing a wider range of tissues and samples within the genus Macaca will be necessary to comprehensively elucidate epigenomic differences.Including more species and samples in future studies will enhance our understanding of epigenetic variations among macaques.
To obtain sample methylation coverage, we extracted methylation information from previously generated BAM files using Bismark (Krueger & Andrews, 2011).As DNA methylation primarily occurs in CpG dinucleotides in most animals, we only studied methylation at the CpG sites.To avoid the influence of sex differences and mitochondrial genomes, we removed methylation information from the X, Y, and MT chromosomes in subsequent analysis.

| DNA methylation landscape
To explore the distribution pattern of DNA methylation in the three species, we used the R package ViewBS to analyze the methylation information for each sample (Huang et al., 2018).We used the R package circlize to determine methylation distribution in the genome (Gu et al., 2014).The methylation-level distribution map was generated using a window size of 50 kb, while the DMR density statistics were calculated using a window size of 1e7 bp.The CpG sites with a methylation coverage ≥5 in the samples were extracted, and PCA and cluster analysis were performed using the R package methylKit, filters covering less than 10% of bases and more than 99.9% of bases, and a PCA plot with 95% confidence ellipses (Akalin et al., 2012).

| Identification and analysis of DMRs
To identify DMRs between species, we analyzed CpG methylation site information using the R package DSS (Feng et al., 2014).The DMR extraction parameters were δ > 0.1 and p < 0.01, and the other parameters were the default values: minlen = 50, minCpG = 3, dis.merge = 50, pct.sig = 0.5.The obtained DMR information was compared with the information of various components of the reference genome using the R package GenomicRanges to search for overlapping areas (Lawrence et al., 2013).The gene annotation files were downloaded from the Ensembl website.The CpG_island information files were downloaded from the UCSC database.Promoter regions were defined as the sequence from 2-kb before the transcription start site to the transcription start site.We used the online website g:Profiler to convert the ENST ID of the promoter into the ENSG ID of the corresponding gene (Raudvere et al., 2019).Finally, we obtained coverage information between DMRs and different components (UTR3, UTR5, promoter, gene body, exon, intron, CpG_island, and CpG_island-related areas).Gene body regions with DMR coverage were defined as DMGs, while promoter regions with DMR coverage were defined as DMPs.To analyze the relationship between DMRs and transposable elements (TEs), we utilized RepeatMasker software to identify the TE regions (Smit, 2004).RepeatMasker uses default parameters to analyze the reference genome (Smit, 2004).
We then searched for overlaps between TEs and DMRs.Enrichment analysis was also performed to explore the biological functions of DMR-related genes.KEGG analysis was performed using g:Profiler with FDR <0.05 (Raudvere et al., 2019), while GO analysis was performed using the R package clusterProfiler with FDR <0.05 (Wu et al., 2021).
We also identified DMRs between male and female individuals using DSS same as the above method.Furthermore, we used the age-related DMRs in rhesus macaques from our previous research (Zhou et al., 2020).Then, to indicate the effects of the age and sex to the species-related DMRs, we identified the overlap DMRs between species-related and age-related, and between species-related and sex-related DMRs, using bedtools v2.25.0 (Quinlan & Hall, 2010).

| Gene expression analysis
Another group of macaques, sequenced differently from the methylation group, was used for gene expression studies.Raw counts and TPM values of M. thibetana (TM; age: 15-18; male n = 2; female n = 3) and M. mulatta lasiota (CR; age: 15; male n = 3; female n = 3) were obtained from published laboratory studies (Wang et al., 2021).
For each sample with strict standards (high-quality paired-end reads with more than 90% of bases with Q-value ≥30), we used NGSQCToolkit v2.3.3 to obtain high-quality reads (clean reads) and remove low-quality paired-end reads or reads containing adaptors (Patel & Jain, 2012).The processed M. fascicularis reads were mapped to the reference genome (M.mulatta, Mmul_10) using HISAT2 v2.1.0(Kim et al., 2015).Each alignment output file was assembled into a separate transcriptome using StringTie v1.3.6 to generate a transcriptome GTF file (Pertea et al., 2015).The genome sequence and annotations ited similar length distribution patterns (CR vs. CE: length range 51-8920 bp, average length 272.86 bp; CR vs. TM: length range 51-4995 bp, average length 278.40 bp; CR vs. CE: length range 51-4994 bp, average length 253.53 bp).At the chromosome level, Chromosome 1 had the largest number of DMRs in all three groups (1119 in the CR versus CE; 5054 in the CR versus TM; 3302 in the CE versus TM groups) (Figure S4).To determine the distribution of DMR density in each chromosome, we divided the length of the DMR on each chromosome by the total length of the chromosome (Figure S5).Results showed that the DMRs in the CR versus CE and CR versus TM groups reached maximum distribution density on chromosome 19 (CR vs. CE = 0.2696%; CR vs. TM = 0.897%).

F
Methylation levels of differential methylation regions (DMRs).(a-c) The percent methylation levels of DMR in groups (a) CR versus CE, (b) CR versus TM, and (c) CE versus TM were plotted."t-test, p < 0.001" means that methylation levels differ significantly between any two macaque species.(d-f)The DMR heatmap of group (d) CR versus CE, (e) CR versus TM, and (f) CE versus TM.Each heatmap shows 2000 DMRs randomly selected from each group of DMRs.(ENSMMUG00000010956), which is involved in the regulation of physiological functions, to analyze the coverage of DMRs on genes.We analyzed the methylation level of SOCS1 and found a significantly higher methylation level in one region of the gene in the TM group compared to the CR and CE groups.Moreover, the related DMR in the CR versus TM group (chr20: 11220602-11,221,998, diff.Methy = −0.349173575,p < 0.01) covered SOCS1 (chr20: 11220610-11,221,257, reverse chain).The related DMR in the CE vs TM group

F
Distribution pattern of DMRs.(a) Circos plot of DMR density distribution.Different color areas in the figure represent the density distribution of DMRs in different groups.The red line in the outermost box represents the centromere.(b-d) The percentage of genomic regions overlapped with DMR.The different DMR groups from (b) to (d) are CR versus CE, CR versus TM, and CE versus TM. enrichment analysis, both groups were enriched in the Rap1 signaling pathway (KEGG:04015), growth hormone synthesis, secretion, and action (KEGG:04935), and axon guidance (KEGG:04360) (Figure S11a,b; Table fied three groups of differentially expressed genes (DEGs, |log2Fold-Change| > 1 and FDR <0.05), namely, CR versus CE (total n = 1677, up-regulated DEGs n = 820, down-regulated DEGs n = 857), CR versus TM (total n = 5162, up-regulated DEGs n = 2371, down-regulated DEGs n = 2791), and CE versus TM (total n = 2895, up-regulated DEGs n = 1245, down-regulated DEGs n = 1650) (Figure 4c-e; Table GO analysis, three GO:BP categories associated with pathogen infection were enriched in the CR versus CE group, including defense response to virus (GO:0051607), defense response to symbiont (GO:0140546), and response to virus (GO:0009615).Several categories related to organismal responses to compounds were F I G U R E 4 Gene expression profiles of three rhesus species.(a) Results of cluster analysis of gene expression profiles in transcriptome samples.(b) Results of PCA of gene expression profiles in transcriptome samples.(c-e) The species-specific DEGs of (c) CR versus CE, (d) CR versus TM, and (e) CEvsTM.The genes marked in the figure are DEGs with LogFC >10.enriched in the CE versus TM group, including cellular response to toxic substances (GO:0097237) and response to oxygen-containing compounds (GO:1901700).KEGG analysis indicated that the CR versus TM group was enriched in pathways previously associated with DMGs, particularly growth hormone synthesis, secretion, and action (KEGG:04935).Notably, DEGs in the CE versus TM group were significantly enriched in serotonergic synapse (KEGG:04726) and platelet activation (KEGG:04611), along with pathways linked to cancer and disease (Figure S14a-c; Table

F
Combined analysis of multi-omics.(a-d) Methylation trends of 4 groups of CR genes with different expression intensities in the gene body and flanking regions.The corresponding expression intensity groups from (a) to (d) were TPM = 0, 0 < TPM ≤1, 1 < TPM ≤10, and TPM ≥10, respectively.(e-g) The multi-omics differential gene. of (e) CR versus CE, (f) CR versus TM, and (g) CE versus TM.Part of the genes marked with asterisks were DMR-related DEGs.(h) KEGG enrichment results of DMR-related DEGs of CR versus TM. (i-k) SNV density of DMR and its flanking regions.In the figure, the abscissa origin DMR* represents DMR and ± 1 kb flanking sequence.±1 ~ ±4 on the abscissa represents the 1 kb region from DMR ± 1 ~ ± 4 kb.The different SNV density groups from (i) to (k) are the hypermethylated DMGs of CR versus CE, CR versus TM, and CE versus TM.