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.

   Keywords: DNA methylation, epigenetics, gene expression, Macaca, neuro
   development

1. INTRODUCTION

   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., [36]2010; Jeremias et al., [37]2018). Furthermore, DNA methylation
   plays a significant role in organismal development, gene expression
   profiles, genomic imprinting, disease susceptibility, and evolutionary
   dynamics (Mendizabal et al., [38]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., [39]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. ([40]2010) reported
   significant differences in methylation levels in 500 genes between
   queen and worker bees, while Kucharski et al. ([41]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., [42]2013),
   especially among twins (Fraga et al., [43]2005; Kaminsky et
   al., [44]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., [45]2007; Xu et al., [46]2022). Despite their similar genetic
   backgrounds, phenotypic differences exist among macaque species (Gibbs
   et al., [47]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., [48]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, [49]2007). Previous analyses have identified the
   pattern of DNA differences among different macaques (Fan et
   al., [50]2014; Fan et al., [51]2018). For example, Fan et
   al. ([52]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., [53]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., [54]2015;
   Hodges et al., [55]2011; Molaro et al., [56]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., [57]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.

2. RESULTS

2.1. 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 [58]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 [59]S2). All
   methylome data were mapped to the M. mulatta genome (Mmul_10), with a
   mean efficiency of 75.75% (Table [60]S2). After using
   deduplicate_bismark, the remaining 86.65%–97.60% of reads were used for
   subsequent analysis (Table [61]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 [62]S1a;
   Table [63]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 [64]S1b).

   We delineated trends in the methylation levels of gene bodies and
   flanking regions (Figure [65]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 [66]S2) showed both
   hypermethylated and hypomethylated regions in each chromosome of the
   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 [67]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 [68]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.

FIGURE 1.

   FIGURE 1
   [69]Open in a new tab

   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”.

2.2. Differentially methylated regions (DMRs) between different macaque
species

   Three groups of DMRs were identified from the 20 auto‐chromosomes: that
   is, CR versus CE (total n = 12,128; hypermethylated DMRs n = 5887;
   hypomethylated DMRs n = 6241; Table [70]S4), CR versus TM (total
   n = 59,165; hypermethylated DMRs n = 42,834; hypomethylated DMRs
   n = 16,332), and CE versus TM (total n = 39,751; 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 [71]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 [72]S4–S8).

FIGURE 2.

   FIGURE 2
   [73]Open in a new tab

   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.

   To further analyze the distribution of methylation differences across
   the genome, we visualized the lengths of the DMRs (Figure [74]S3).
   Interestingly, the DMRs in the three groups exhibited 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 [75]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 [76]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%). 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 [77]3a). First, the distribution of DMRs in
   each chromosome was not uniform but fluctuated, with some 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.

FIGURE 3.

   FIGURE 3
   [78]Open in a new tab

   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.

   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 [79]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 [80]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 [81]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 [82]S8a–c;
   Table [83]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 [84]S9a–c; Table [85]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
   (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 (chr20: 11220711–11,221,349,
   diff.Methy = −0.317753274, p < 0.01) also covered most regions of 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 [86]S10a–b;
   Table [87]S11). Based on KEGG pathway 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 [88]S11a,b; Table [89]S12). DMPs in the
   CE vs TM group and shared among all three groups showed no significant
   enrichment.

2.3. Gene expression analysis

   Several studies have investigated gene expression differences at
   different physiological stages within specific Macaca species (Lan et
   al., [90]2020; Yan et al., [91]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 [92]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 [93]4b). The distance between CR and CE was
   closer than that between CR and TM (Figure [94]4b). We next analyzed
   and identified 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 [95]4c–e;
   Table [96]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.

FIGURE 4.

   FIGURE 4
   [97]Open in a new tab

   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.

   To elucidate the biological functions of the DEGs, we conducted GO and
   KEGG enrichment analyses (Figure [98]S13a–c; Table [99]S14). For 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 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 [100]S14a–c; Table [101]S15).

   Among the three groups, we identified 339 shared DEGs
   (Figure [102]S12). Based on further GO and KEGG enrichment analysis
   (Tables [103]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., [104]2013).

2.4. 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 [105]5a–d; Figure [106]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, [107]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 [108]5a–d;
   Figure [109]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.

FIGURE 5.

   FIGURE 5
   [110]Open in a new tab

   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.

   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 [111]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 [112]5e–g). Enrichment analysis was performed as above
   (Figure [113]S16a–c; Tables [114]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 [115]S16a) and KEGG pathways
   related to phenotypic differences, including growth hormone synthesis,
   secretion, and action (KEGG:04935) (Figure [116]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 [117]S16b), and KEGG
   pathways (e.g., platelet activation (KEGG:04611)) (Figure [118]S16c)
   also enriched in the CR versus TM group (Figure [119]5h).

2.5. 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., [120]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 [121]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 [122]S17a,b), but an increase in base mutation density in the
   CE versus TM hypomethylated DMRs (Figure [123]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.

3. DISCUSSION

   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., [124]2014; Higashino et
   al., [125]2012; Osada et al., [126]2015; Yan et al., [127]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., [128]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., [129]2013; Howard et al., [130]2011;
   Massart et al., [131]2014; Provençal et al., [132]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.

3.1. Genetic differences and DNA methylation profiles

   Macaques are the most widely distributed nonhuman primates and serve as
   valuable animal models in biomedical research (Gibbs et
   al., [133]2007). Previous studies have reported the whole‐genome
   sequences of M. mulatta lasiota, M. fascicularis, and M. thibetana (Fan
   et al., [134]2014; Fang et al., [135]2011; Osada et al., [136]2015),
   revealing high similarity at the DNA level (with genome divergences
   around 0.28%–0.40%) (Fan et al., [137]2014; Osada et al., [138]2008;
   Osada et al., [139]2010). 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., [140]2014; Perelman et al., [141]2011).

   Various studies have investigated distinctions among macaque species
   from different perspectives. For instance, Li et al. ([142]2020) used
   copy number variation (CNV) to construct a neighbor‐joining tree of
   nine Macaca species. Song et al. ([143]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. ([144]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., [145]2019).

   Fan et al. ([146]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., [147]2005). The three macaque species
   studied here exhibit variations in their habitats and behaviors
   (Fleagle, [148]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 [149]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., [150]2016).

3.2. 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, [151]1996; Shively et al., [152]1997; Shively &
   Wallace, [153]2001), as well as hypothalamic‐pituitary‐adrenal axis
   function at an endocrine level (Kaufman et al., [154]2007). The growth
   and development of organisms are influenced by the endocrine signaling
   system in conjunction with signal transduction pathways
   (Neave, [155]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 [156]S10 and S12).

   Notably, the growth hormone synthesis, secretion, and action pathway
   (KEGG:04935) was enriched in all DMG groups (Table [157]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., [158]2006; Krebs & Hilton, [159]2001).
   Overexpression of SOCS1 can completely inactivate the growth hormone
   (GH) signaling pathway (Ram & Waxman, [160]1999). Furthermore, SOCS1
   can inhibit GH receptor‐dependent JAK2 tyrosine phosphorylation and
   decrease GH activity (Ram & Waxman, [161]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., [162]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 down‐regulation 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., [163]2010). Studies have reported significant
   differences in blood glucose between CR and TM, with TM exhibiting
   relatively higher blood sugar levels (Yang et al., [164]2010). Fan et
   al. ([165]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 [166]S10 and S12), as well as enrichment in insulin receptors
   (INSR; ENSMMUG00000028907), which are essential components of
   insulin‐related pathways (Tables [167]S10 and S12). In mammals, changes
   in insulin receptor levels are of physiological importance (Payankaulam
   et al., [168]2019). In animals, the physiological function of INSR is
   closely related to tyrosine kinase signal transduction (Ward et
   al., [169]2007). Heterozygous INSR mutant mice generally exhibit normal
   growth and fertility but show deficiencies in downstream signaling
   (Batista et al., [170]2019). In humans, heterozygous carriers of null
   mutations of INSR exhibit abnormal glucose tolerance, highlighting the
   importance of INSR gene expression (Psiachou et al., [171]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 [172]S10and S12), which are
   essential in the regulation of growth, development, and metabolism
   (Cheng et al., [173]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.

3.3. 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 [174]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, [175]2011; Pinard et al., [176]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 receptor‐related 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., [177]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.

4. CONCLUSIONS

   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.

5. MATERIALS AND METHODS

5.1. Samples of BS‐seq

   Peripheral blood samples were collected from M. mulatta lasiota (CR;
   age: 1.5–18; male n = 3; female n = 2), M. fascicularis (CE; age:
   2.3–11; male n = 2; female n = 3), and M. thibetana (TM; age: 3–15;
   male n = 3; female n = 2). Specific sample information is provided in
   the [178]Supplementary Materials. We used whole blood from all species
   to represent general epigenomes of peripheral blood tissues
   (Hernando‐Herraez et al., [179]2013; Reinius et al., [180]2012).
   Interspecies DNA methylation variation may also result from genetic
   differences and/or environmental factors and other stochastic events
   (Hernando‐Herraez, Heyn, et al., [181]2015). Therefore, we used captive
   samples from all species to ensure homogenous environments.

5.2. BS‐seq data processing and analysis

   Quality control and adapter trimming were performed using the
   NGSQCToolkit_v2.3 (Patel & Jain, [182]2012). Bismark was used to map
   the sequences after quality control with the reference genome (M.
   mulatta, Mmul_10), using the following parameters: ‐N 1‐L 20 ‐‐
   Score_min L,0, −0.2 (Krueger & Andrews, [183]2011; Zhou et
   al., [184]2020). To obtain sample methylation coverage, we extracted
   methylation information from previously generated BAM files using
   Bismark (Krueger & Andrews, [185]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.

5.3. 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., [186]2018). We used the R
   package circlize to determine methylation distribution in the genome
   (Gu et al., [187]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., [188]2012).

5.4. Identification and analysis of DMRs

   To identify DMRs between species, we analyzed CpG methylation site
   information using the R package DSS (Feng et al., [189]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., [190]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., [191]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, [192]2004).
   RepeatMasker uses default parameters to analyze the reference genome
   (Smit, [193]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., [194]2019), while GO
   analysis was performed using the R package clusterProfiler with
   FDR <0.05 (Wu et al., [195]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., [196]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, [197]2010).

5.5. 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., [198]2021).
   Transcriptome data of M. fascicularis (CE; age: 15–17; male n = 4;
   female n = 6) were sequenced and analyzed using the same procedures as
   previous studies (Wang et al., [199]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, [200]2012). The processed M. fascicularis reads were mapped to
   the reference genome (M. mulatta, Mmul_10) using HISAT2 v2.1.0 (Kim et
   al., [201]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., [202]2015). The genome sequence
   and annotations were downloaded from Ensembl release‐98. For M.
   fascicularis, specific Ensembl GTF files (Mmul_10.98) were used as the
   reference annotation file to guide the assembly process to obtain the
   TPM expression value and raw read counts for each gene and transcript.

   We used the generalized linear model (GLM) in the R package DESeq2
   (Fitting formula: ~species + sex) to analyze pairwise differential
   expression in transcriptome data of the three species (Love et
   al., [203]2014). The Wald test was used for pairwise comparisons of the
   differentially expressed transcriptomes, and the GLM was fit to
   estimate size factors and dispersity of all samples from each species.
   The significance threshold of DEGs was set as adjusted p < 0.05 and abs
   (log2Fold‐Change) > 1. To eliminate interference, we excluded genes
   with count values <10 and loci that were not autosomal. The obtained
   gene expression information was used for subsequent correlation
   analysis.

5.6. Analysis of genetic differences in DMRs

   We analyzed the base mutation information in the DMRs and flanking
   regions to investigate the relationship between genomic differences and
   DNA methylation. Genomic data for the three species were obtained from
   unpublished laboratory work. Previous studies excluded insertions and
   deletions from the genome when calculating base differences between
   pairs of species (Hernando‐Herraez, Heyn, et al., [204]2015).

FUNDING INFORMATION

   This work was supported by the National Natural Science Foundation of
   China (32070413).

CONFLICT OF INTEREST STATEMENT

   The authors declare that they have no competing interests.

Supporting information

   Appendix S1.
   [205]Click here for additional data file.^ (9.4MB, pdf)

   Appendix S2.
   [206]Click here for additional data file.^ (13.3MB, xlsx)

   Appendix S3.
   [207]Click here for additional data file.^ (132.6KB, pdf)

ACKNOWLEDGMENTS