Genes & Genetic Systems
Online ISSN : 1880-5779
Print ISSN : 1341-7568
ISSN-L : 1341-7568
Full papers
DNA methylation plays an important role in iron-overloaded Tibetans
Qin ZhaoZhijing GeSuhong FuSha WanJing ShiYunhong WuYongqun Zhang
Author information
JOURNAL OPEN ACCESS FULL-TEXT HTML
Supplementary material

2022 Volume 97 Issue 2 Pages 55-66

Details
ABSTRACT

The prevalence of iron overload in Tibetans in Tibet is higher than that in Han. DNA methylation (DNAm) is closely related to iron metabolism and iron level. Nevertheless, the epigenetic status of Tibetans with iron overload is unknown, and we therefore aimed to explore whether the phenomenon observed in the Tibetan population is regulated by epigenetics. The results showed that 2.26% of cytosine was methylated in the whole genome, and that the rate of CG cytosine methylation was higher in individuals in the iron overload (TH) group than in those in the iron normal (TL) group. We analyzed differentially methylated genes (DMGs) in whole-genome bisulfite sequencing data from the TH and TL groups of high-altitude Tibetans. Protein–protein interaction and pathway analyses of candidate DMGs related to iron uptake and transport showed that epigenetic changes in 10 candidate genes (ACO1, CYBRD1, FLVCR1, HFE, HMOX2, IREB2, NEDD8, SLC11A2, SLC40A1 and TFRC) are likely to relate to iron overload. This work reveals, for the first time, changes of DNAm in Tibetan people with iron overload, which suggest that DNAm is a mechanism underlying differences in iron content between individuals in the high-altitude Tibetan population. Our findings should contribute to the study of iron metabolism and the overall health status of Tibetans.

INTRODUCTION

Iron is an essential element for growth and survival, and is also a cofactor for a large number of heme and non-heme ferritin proteins. Similarly, iron exists in heme or various non-heme forms in the diet (Nicolas et al., 2002). Non-heme iron is transported across the apical membrane of intestinal cells through divalent metal transporter 1 (DMT1) and transported to the circulatory system through the iron gate protein ferroportin 1 (FPN1). Newly absorbed iron is combined with plasma transferrin (TF) and distributed throughout the body, and is used for red bone marrow, which has a high requirement for iron. Iron-loaded TF binds to transferrin receptor (TFR) on the surface of most somatic cells. After the complex is endocytosed, iron enters the cytoplasm through DMT1, and can then be used for metabolic functions, stored in cytosolic ferritin, or exported from the cell via FPN1 (Anderson and Frazer, 2017). Notably, because ferritin is a soluble tissue protein that stores iron in the body, serum ferritin (SF) can correctly reflect the iron reserves of the liver (Rostoker et al., 2015). The main regulator of iron homeostasis in the body is hepcidin (HAMP), which is mainly produced in the liver. Liver cells produce ferritin due to excessive iron and inflammatory stimulation (such as interleukin-6), while the expression of HAMP is down-regulated due to hypoxia, anemia and erythropoiesis (Kawabata, 2017). In addition, the liver proteins human hemochromatosis protein (HFE), transferrin receptor 2 (TFR2), hemojuvelin (HJV) and bone morphogenetic proteins (BMPs) are necessary regulatory factors to activate the synthesis of ferritin (Borch-Iohnsen et al., 2009), and deficiency or mutation of their genes is associated with primary or secondary iron overload syndrome (Deugnier et al., 2017).

DNAm can affect the expression of these genes by regulating methylation levels in the promoter regions of genes for iron metabolism such as FPN1, TFR2, HAMP, HJV and BMP6 (Duan et al., 2020). When iron exceeds the body’s defensive ability, it will lead to oxidative damage and tissue dysfunction (de Lima et al., 2019), and iron-dependent cell ferroptosis. Iron overload may be asymptomatic or accompanied by major diseases of the liver, heart, endocrine glands, joints or other organs. Under the conditions of high altitude and low oxygen partial pressure, the body needs to absorb a large amount of iron, which participates in the synthesis of hemoglobin and then in the transport of oxygen. Furthermore, the iron absorption rate of in Tibetan men living at high altitudes is much higher than the average of Chinese people (Jia et al., 2005). Hence, the hypoxic environment causes a compensatory increase, in iron metabolic rate and in iron storage in the body, resulting in a significant increase in SF and TF in Tibetan people and in Han people when they emigrate to high altitude (Zhao et al., 2003). The saturation of ferritin and TF in Tibetan people (42.91%) is higher than that in Han people (31.96%), and adult males on the Qinghai–Tibet Plateau have a higher prevalence of iron overload at 4.2% (Suo et al., 2018). Moreover, studies have shown that the iron absorption of normal mice exposed to hypoxia for three days is increased 2–3-fold (Raja et al., 1990). Hypoxia is associated with a significant decrease of HAMP gene expression in the liver, which may be due to the positive response of the body to inflammation caused by hypoxia (Nicolas et al., 2002). Most importantly, DNAm may control the expression of iron-induced genes in liver, including HAMP, so the influence of epigenetics on iron homeostasis deserves further study (Sharp et al., 2018). Therefore, gene expression at the level of the individual can be controlled by many genetic or epigenetic factors. One possible mechanism is through DNAm. More importantly, DNAm is the main epigenetic mechanism of gene silencing. Increased DNAm leads to gene silencing, while decreased methylation activates gene expression (Adrian, 2002). In addition, DNAm may change the susceptibility to diseases or phenotypic traits of normal people. Iron metabolism is closely associated with DNAm, and the level and status of iron affect DNAm (Yamamoto et al., 2016).

In this study, we used whole-genome bisulfite sequencing (WGBS) for the first time to further study the epigenetics of iron homeostasis in Tibetan people and identified differentially methylated regions (DMRs) and related genes that may participate in the regulation of iron overload in the Tibetan population, which preliminarily revealed that iron overload in Tibetans may be regulated by the epigenetic mechanism of DNA methylation. This discovery illuminates the relationship between epigenetic status and iron metabolism in Tibet.

RESULTS

Genome-wide DNA methylation profiling in Tibetans

The raw data for the iron overload (TH) group and iron normal (TL) group (for details, see Materials and Methods) were 308.03 Gb and 303.85 Gb, respectively. The WGBS data information for group TH and group TL is recorded in Table 1. The clean reads of 285.75 Gb (average mapping rate 88.51%) (TH) and 281.52 Gb (88.42%) (TL) could be independently mapped to the human reference genome GRCh37.p13. Any fuzzy mapping and repeated reads were removed from further analysis.

Table 1. Summary of WGBS dataset
GroupSampleRaw readsClean readsMap ReadsRaw bases (Gb)Clean bases (Gb)Mapping rate (%)Q30 (%)
THTH1695203434662074430586457942104.2896.9188.5894.78
TH2691506272656501310582206976103.7395.9888.6894.64
TH3666776206634418148559918298100.0292.8688.2694.76
TLTL1679614888638837386559353664101.9493.4387.5694.47
TL2693907654662106988589101946104.0997.1388.9795.07
TL365215742262099648055095929097.8290.9688.7294.74

DNA methylation patterns in Tibetan genomes

The methylation level of cytosine for each sample was calculated. On average, 8,694,998 and 4,807,892 methylated cytosines (mCs) were detected in the TH and TL groups. The classification of mCs showed that the proportions in the TH group (average value 2.26%) and the TL group (average value 2.02%) were similar (Table 2). The proportion of methylated C sites (mC/CoveredC) in CG, CHG and CHH environments was 61.61%, 0.0197% and 0.0173%, respectively, in the TH group (Table 2). In summary, there are three classifications of methylation in Tibetans: mCG, mCHH (where H is A, C or T) and mCHG.

Table 2. Proportion of methylated cytosines in Tibetans of TH and TL groups
SamplemC (%)mCG (%)mCHG (%)mCHH (%)
TH12.6973.320.0230.020
TH21.8752.140.0170.015
TH32.2159.380.0190.017
TL11.4041.490.0150.012
TL22.4969.900.0210.019
TL32.1759.570.0190.017

To further compare the genome-wide distribution and methylation level of different functional genomic elements between the two groups, we analyzed the methylation status of six regions: the promoter, 5′ UTR (untranslated region), exons, introns, 3′ UTR and distal intergenic region. The highest level displayed by the CG context is far higher than that of the CHG and CHH contexts (Supplementary Fig. S1). Besides, the methylation levels of CG, CHG and CHH in the TL group and the TH group showed similar trends. Moreover, in the context of CG, the methylation level of DMRs located before the TSS was reduced, and, compared with the TL group, the TH group displayed further hypomethylation. In the environment of CHG and CHH, the trend was much milder. Compared with the TL group, the methylation level of the TH group was higher in most areas.

DMR information

We identified 201,591, 202,334 and 203,351 DMRs in the CG, CHG and CHH sequence contexts between TH and TL (methylation difference ≥ 10%, q-value < 0.05). In DMRs, 254,362 (CG: 84,609; CHG: 84,766; CHH: 84,987) sites were hypermethylated, while 352,914 (CG: 116,982; CHG: 117,568; CHH: 118,364) were hypomethylated. Therefore, DMRs are mainly located in the intron and exon regions. In the context of CG, CHG and CHH, respectively, 7,563, 7,587 and 7,625 DMRs were in the promoter.

GO and KEGG enrichment analysis of DMGs

To explore the gene function of DMGs in iron-overloaded Tibetans, GO and KEGG pathway analysis were used to identify DMGs. Meanwhile, given that CG context is the main type of DMR, the focus in functional enrichment analysis was CG DMGs. The CG-type DMRs cover a total of 12,274 DMGs. Thus, GO enrichment showed that DMGs were significantly enriched in a variety of biological processes. The top 10 most significantly enriched GO terms of molecular function, biological processes and cellular components were related to cell development, protein binding, signal regulation and cytoskeleton composition (Supplementary Fig. S2). Enrichment analysis of KEGG pathways showed that DMGs were significantly correlated with cancer, focal adhesion, platelet activation, cAMP signaling and other pathways (Supplementary Fig. S3).

Certain DMGs are involved in processes such as iron ion transport across the membrane, iron transport and cell–cell interaction through iron ions. The 10 most significantly enriched GO terms related to iron metabolism of DMGs between the TH group and the TL group are shown in Table 3. These results indicate that these DMGs may affect iron absorption, transport and metabolism, and thus lead to differences in iron content in the Tibetan population at high altitudes.

Table 3. GO terms related to iron metabolism enriched in DMGs
GO term*CategoryGene
iron ion transmembrane transporter activity (GO:0005381)molecular functionSLC25A37, TTYH1, SLC40A1
ferrous iron transmembrane transporter activity (GO:0015093)molecular functionSLC39A14, SLC11A2
ferrous iron transport (GO:0015684)biological processSLC11A2
iron-responsive element binding (GO:0030350)molecular functionFECH, IREB2, ACO1
iron ion transmembrane transport (GO:0034755)biological processTTYH1, SCARA5, SLC40A1
mitochondrial iron ion transport (GO:0048250)biological processSLC25A37
ferrous iron import (GO:0070627)biological processSLC11A2
cellular response to iron ion (GO:0071281)biological processSLC11A2, TFAP2A
iron ion import into cell (GO:0097459)biological processHFE, PICALM
response to iron ion (GO:0010039)biological processBCL2, CYBRD1, SLC11A2, ABAT, SLC6A3, BMP6, MDM2, CPOX, DRD2
*,  P < 0.05.

Iron metabolism candidate PPI network building and analysis of DMGs

To further explore the relationship between DNA methylation and iron overload in the Tibetan population, we screened candidate DMGs involved in iron metabolism, and established two limits for correlation analysis. First, 1,150 genes related to iron metabolism, high altitude hypoxia, erythropoiesis and cell damage were screened in NCBI. Using these four gene sets to make a Venn diagram (Supplementary Fig. S4), 62 candidate DMGs were extracted and related to iron metabolism in the Tibetan plateau population. Next, 20 DMGs involved in iron absorption, transport and metabolism appearing in Table 3 and 62 common genes related to iron metabolism obtained from the Venn diagram (Supplementary Fig. S4) were used to analyze the interaction between these DMGs, using the STRING website and Cytoscape to make a protein–protein interaction (PPI) network (Fig. 1). At the same time, the STRING database was used for GO, KEGG and Reactome pathway enrichment analysis of DMGs. The information related to DMGs participating in iron metabolism is shown in Supplementary Table S1. PPI network analysis predicted that DMGs involved in the steady state of iron ions in cells had strong correlations.

Fig. 1.

Construction of a PPI network map related to iron overload in the Tibetan population. According to the interaction index (confidence > 0.70), interactions between DMGs related to iron metabolism were analyzed via the STRING website using Cytoscape software. Individual nodes are hidden; the interaction between genes is represented by green arrows and dots, with dots representing the sources and arrows representing the targets of interactions. The red dots are key candidates for iron overload-related DMGs.

Regulatory effect of DMGs on the iron level in the Tibetan population at high altitudes

The transport of iron between cells is mediated by TF. Iron can enter and leave cells in the form of heme and iron carriers. When iron enters cells through the main pathway (TF endocytosis), heme is synthesized and iron–sulfur clusters gather in mitochondria (Kurz et al., 2008; Hower et al., 2009). Thus, the absorption and transport of iron in cells is closely related to the iron content in the body. To explore candidate genes related to DNA methylation and iron overload, 10 DMGs of iron metabolism (Table 4), ACO1, CYBRD1, FLVCR1, HFE, HMOX2, IREB2, NEDD8, SLC11A2, SLC40A1 and TFRC, were analyzed on the Reactome website and found to participate in the ‘Iron uptake and transport’ pathway (Fig. 2). CYBRD1 participated in the process in which heme reduces Fe3+ to Fe2+; SLC11A2 cotransports Fe2+ and H+ from the extracellular domain to the cytoplasm; SLC40A1 transports Fe2+ from the cytosol to the extracellular region, and participates in the oxidation of Fe2+ to Fe3+; the function of HMOX1 is to split heme; FLVCR1 transports heme from the cytosol to the extracellular matrix; ACO1 combines with 4Fe–4S to isomerize citrate to isocitrate; NEDD8 binds to CUL1 (in the SKP1–CUL1–FXBL5 complex) and then ubiquitinates IREB2.

Table 4. Methylation information for DMGs involved in the ‘Iron uptake and transport’ pathway
Gene symbolChromosomeMeth. diff.
(TH vs. TL)
P-valueAnnotation
ACO1chr9(32385001-32386000)−21.7278.16E-04intron
chr9(32425001-32426000)−11.6584.72E-05exon
chr9(32425001-32426000)−11.6584.72E-05intron
CYBRD1chr2(172395001-172396000)17.2439.15E-04intron
chr2(172410001-172411000)−15.0913.80E-04exon
chr2(172410001-172411000)−15.0913.80E-04intron
chr2(172414001-172415000)19.5408.47E-04exon
chr2(172414001-172415000)19.5408.47E-04intron
FLVCR1chr1(213069001-213070000)24.7086.20E-14exon
chr1(213069001-213070000)24.7086.20E-14intron
chr1(213069001-213070000)19.5236.48E-05exon
chr1(213069001-213070000)19.5236.48E-05intron
HFEchr6(26088001-26089000)20.1401.31E-06intron
chr6(26088001-26089000)22.5287.38E-06intron
HMOX2chr16(4523001-4524000)11.6452.03E-04TSS
IREB2chr15(78731001-78732000)−13.0272.65E-04exon
chr15(78731001-78732000)−13.0272.65E-04intron
NEDD8chr14(24688001-24689000)29.1002.02E-04intron
chr14(24694001-24695000)−11.1052.01E-07intron
SLC11A2chr12(51373001-51374000)−20.3744.12E-05exon
chr12(51373001-51374000)−20.3744.12E-05intron
chr12(51405001-51406000)26.9167.04E-04intron
SLC40A1chr2(190431001-190432000)−11.0058.46E-04intron
chr2(190444001-190445000)−16.9439.85E-04exon
chr2(190444001-190445000)−16.9439.85E-04intron
chr2(190446001-190447000)−12.2903.27E-05exon
chr2(190446001-190447000)−12.2903.27E-05intron
chr2(190449001-190450000)−19.8941.26E-04TSS
TFRCchr3(195805001-195806000)−13.9585.46E-05intron

Note: Meth. diff., difference in methylation levels between TH and TL. A positive number means that the methylation levels of this region in the TH group are higher than those in the TL group, and a negative number means that the methylation levels of this region in the TH group are lower than those in the TL group.

Fig. 2.

Annotation of the ‘Iron uptake and transport’ pathway of the DMGs using the Reactome tool. Yellow areas indicate where the candidate DMGs participate in the pathway.

Mutations in key proteins that transport, perceive and metabolize iron and promote its utilization can lead to the disorder of iron homeostasis, resulting in iron overload. These ten DMGs also play a role in the “cellular iron ion homeostasis” biological process (GO:0006879). Both the promoters and the gene body can regulate gene expression. Finally, our original data show that the difference in gene methylation is mainly distributed in the gene body. Therefore, among the ten genes, the methylation levels of FLVCR1, HMOX2 and HFE genes are up-regulated while ACO1, IREB2, SLC40A1 and TFRC are hypomethylated. SLC11A2, CYBRD1 and NEDD8 may contain more DMRs in the gene body (exon, intron), so DNA methylation is inconsistent (Fan et al., 2020). These 10 DMGs not only play an important role in the ‘Iron uptake and transport’ (R-HSA-917937.3) pathway, but also are enriched in iron metabolism or GO terms related to iron homeostasis, such as cellular iron homeostasis, cell response to iron ions and iron ion transmembrane transport. These results illuminate the mechanism of epigenetic disorder in iron-overloaded Tibetans.

Iron metabolism genes in biological information analysis based on the GEO database

To verify the relationship between these ten genes and iron overload, based on the GEO database, we used GEO2R to analyze four data sets, find changes of key genes’ expression during iron overload, and derive the relationship between the key genes obtained and iron overload using transcriptome sequencing data. The results are listed in Table 5. In these four data sets, ten genes have different expression in iron overload on the whole; some of these are opposite to methylated genes in our study. However, the CYBRD1, SLC11A2 and TFRC genes in some data sets have a highly significant logFC value (P < 0.01). Generally speaking, these ten key genes are expressed to different degrees during the final iron loading process of the body, and different expression trends in each data set may be related to different diseases and experimental conditions.

Table 5. Differential expression (logFC values) between iron overload and iron normal groups for our ten key DMGs in four GEO data sets
Gene symbolGSE160727GSE7357GSE9726GSE10421
ACO1−0.0700.2400.069−0.347
CYBRD10.3793.126*0.089−1.777*
HFEn.c.−1.474−0.0970.294
HMOX2−0.2640.0460.052−0.195
FLVCR1−0.206n.c.n.c.n.c.
IREB2−0.0170.3840.022−0.082
NEDD8−0.373−0.002−0.054−0.224
SLC11A2−0.2451.704−0.095−0.420*
SLC40A11.8210.2640.0930.162
TFRC0.219−0.870−0.485−1.457*
*,  P < 0.01; n.c., not calculated.

DISCUSSION

Epigenetic variation can regulate the gene expression and phenotypic traits of the population

Previous studies have shown that excessive or insufficient iron is harmful to the body, and it is therefore very important to control the dynamic balance of iron by regulating its input, storage and discharge (Fibach and Rachmilewitz, 2017). The effect of iron on some diseases is directly related to iron metabolism, or secondary to basic diseases (Reichert et al., 2017), whereas genetic cofactors, ambient environments and an excessively rich diet might trigger a multifactorial and dynamic process to elevate the iron store in the body. Accordingly, in liver cancer tissues, hepcidin (HAMP) gene transcription is inhibited, and the promoter is hypermethylated. Similarly, long-term exposure of colon cells to iron can result in significant hypomethylation, which may be one of the pathogenic factors in iron-mediated intestinal diseases (Tao et al., 2014). Furthermore, DNA methylation may control the expression of HAMP and other liver iron-sensing genes, and may lead to gene silencing. The influence of epigenetics on iron homeostasis may have a similar relationship with the risk of iron metabolism disorders (Sharp et al., 2018). One of the mechanisms of DNA methylation in regulating gene expression is the combination and interaction between transcription factors and DNA (Hu et al., 2013). A recent study has shown that both promoter and gene region can regulate gene expression (Tang et al., 2015). However, DNA methylation in the promoter is negatively correlated with transcription, and methylation in the gene body has a positive or negative effect on transcription. CpG islands, genomic regions rich in CG, are mostly hypomethylated in normal living conditions (Tang et al., 2015). It has been proved that epigenetic variation has an effect on regulating gene expression and phenotypic traits in a population (Fibach and Rachmilewitz, 2017). Epigenetic DNA methylation affects many biological processes, including hypoxia adaptation (Zhang et al., 2017a), and the Tibetan population shows unique genetic adaptability to high-altitude environments. Methylation modification plays an important role in iron overload. A 2-fold increase in the brain ferritin level of H67D HFE mutant mice is related to inhibition of the activity and expression of DNA methyltransferases (DNMTs) and a decrease in the overall methylation level of the brain (Ye et al., 2019). The significant decrease in SF levels observed in rats after bariatric surgery may be due to an increased DNA methylation level in the HAMP promoter region in liver tissue and the suppression of HAMP gene expression (Huang et al., 2020). By inhibiting the expression of BDH2 in CD4+ T cells of lupus erythematosus patients, the iron concentration in CD4+ T cells of these patients can be increased and the overall DNA methylation level can be reduced (Zhao et al., 2018). Therefore, we hypothesized that DNA methylation is involved in the regulation of iron homeostasis in a hypoxic environment.

Iron-overloaded Tibetans have higher methylation levels

Single-base methylation of iron-overloaded Tibetans was deciphered by WGBS. The methylation level of the whole genome in the TH group was about 2.26%, and the proportion of mCG, mCHH and mCHG in the TH group was slightly higher than that in the TL group, which was consistent with the results observed in iron overload mice, where iron overload can enhance genomic DNA methylation (Yamamoto et al., 2016). The ratio and level of CG methylation types in the genome is the highest, which is consistent with previous findings in humans (Rademacher et al., 2014). DNA methylation in transcribed genomic regions is related to the regulation of alternative splicing. This study also obtained the sequence preference of mCGs, mCHGs and mCHHs. CAG is the most common sequence motif at mCHG sites, while CAT is the most common sequence motif at mCHH sites. There was no significant difference between the two groups in sequence motifs. In addition, a significant decrease in methylation levels near TSSs was observed in the TH group, which was similar to previous results (Zhang et al., 2017b). The methylation level of the TH group was lower than that of the TL group in the CG environment. We proved that CGIs were enriched in many gene bodies and regulatory regions. DNMTs play an essential role in DNA methylation and transcriptional regulation in the genome (Loaeza-Loaeza et al., 2020). We found that DNMT1, DNMT3a and DNMT3b have significant DMRs in this study, indicating that methylation of DNMTs is related to this phenomenon. These results provide conclusive evidence that iron overload leads to epigenetic effects in the whole genome, and provide a comprehensive analysis of genes and pathways related to iron overload.

Epigenetic regulation is involved in the molecular mechanism of iron homeostasis in the Tibetan plateau population

We generated an interactive network of candidate DMGs (Fig. 1) to determine whether DMGs play a decisive role in iron metabolism. Our analysis shows that DMGs reflect not only the influence of iron overload on Tibetans but also the specific signals related to ‘iron absorption and transport’. These DMGs showed methylation changes in genes and regulatory regions. Overall, a comprehensive analysis of affected genes and pathways provides evidence that iron overload might occur in response to genome-wide epigenetic effects. We analyzed the Reactome pathway of iron metabolism-related DMGs, and found 10 DMGs (Fig. 2); remarkably, HMOX2 is closely related to high-altitude hypoxia, HFE and NEDD8 are highly related to cell damage, and six genes (ACO1, FLVCR1, HFE, SLC11A2, SLC40A1 and TFRC) are closely related to the biosynthesis of red blood cells. Iron overload can lead to oxidative stress and DNA damage, which may be related to compensatory polycythemia in the Tibetan population.

Cellular iron homeostasis is maintained through coordinated post-transcriptional regulation of genes associated with iron metabolism. The absorption of iron in the diet is mainly in the duodenum. It was previously found that ferrous ions were absorbed from the intestinal lumen through the apical membrane of intestinal cells, and released into the portal circulation through the basilar lateral membrane. SLC11A2 encodes the divalent cation transporter DMT1, which is the only known transmembrane iron transporter involved in iron uptake in cells, and mediates the absorption of ferrous ions into the cells (Tandy et al., 2000) and the cellular response to iron ions. Diseases associated with SLC11A2 include hypochromic microcytic anemia with iron overload 1 and hypochromic microcytic anemia. SLC40A1 encodes the metal transporter FPN1, localized to the basolateral membrane of intestinal epithelial cells, which mediates the entry of ferrous iron entering the portal vein (Schimanski et al., 2005), and is considered as the key coordinator of the iron balance between intracellular and systematic iron homeostasis (Le Gac et al., 2013). FPN1 is highly expressed in macrophages (iron outflow-mediated heme decomposition), which mediates iron efflux during hemoglobin rupture (Schimanski et al., 2005). Correspondingly, CYBRD1 is an iron-regulated protein, a functional ferric reductase that can reduce Fe2+ in the brush border membrane and the airway (McKie et al., 2001). Our results showed that GO annotations related to this gene include ferric-chelate reductase activity and bound metal ions. In the case of increased erythropoiesis, SLC11A2, CYBRD1 and SLC40A1 were up-regulated (Zaahl et al., 2004; Viatte et al., 2005). HFE binds to TFR and reduces its affinity for iron-loaded transferrin. The HFE gene is also involved in the perception of extracellular and intracellular iron levels, while inhibiting the expression of HAMP under conditions of intracellular and extracellular iron overload (Mehta et al., 2017). Indeed, the iron level of HFE knockout mice is higher than that of normal mice (Coppin et al., 2007). Hence, SLC40A1 gene silencing induces iron retention and enhances ferritin synthesis in human macrophages (Gallí et al., 2015). Ultimately, the SLC40A1, HFE and TFR2 genes are related to the phenotypes of hyperferremia and iron overload (Wang et al., 2017).

Under hypoxia, the specific expression of the HMOX2 gene can regulate hemoglobin metabolism downstream of the hypoxic metabolic pathway, leading to more effective hemoglobin degradation, which is helpful to maintain a relatively low hemoglobin level in high-altitude areas and to reduce cell damage caused by hypoxia (Yang et al., 2016). Moreover, heme transporter 1 (FLVCR1) is located on the plasma membrane and plays a key role in erythropoiesis by protecting developing erythrocytes from heme toxicity (Quigley et al., 2004).

Because Tibetans have lived on the Qinghai–Tibet Plateau for thousands of years and have developed mechanisms to adapt to the hypoxic environment, and because iron is an important part of hemoglobin and erythrocytes, we speculate that Tibetans need to consume more iron to produce more erythrocytes to adapt to this environment. Iron transport and absorption are closely related to iron overload, and any disruption of this link may lead to iron overload in the blood. According to the function of DMGs that participate in iron absorption and transmission pathways (Fig. 2), we propose a molecular mechanism for iron overload. Changes in the methylation levels of the SLC11A2, CYBRD1, HFE and TFRC genes in people with iron overload make the body constantly absorb iron. CYBRD1 combines with heme to transform Fe3+ to Fe2+, and DMT1 continues to deliver extracellular Fe2+ and H+ into the cell. HFE binds to TFRC dimer and reduces its affinity for iron-loaded TF (Feder et al., 1998). Holo-transferrin (holoTF) and TFRC specifically bind to form a dimer, which then transports iron from the plasma membrane to endosomes (Willingham et al., 1984). When the endosomal pH reaches 6.0, more Fe3+ dissociates from the holoTF:TFRC dimer (Hémadi et al., 2006), and the iron ion level in the blood increases. Free TF and TFR1 circulate to the cell surface for periodic iron binding and uptake (Willingham et al., 1984). The intestinal absorption of dietary iron is unbalanced, the iron released by macrophages increases, and the cells are insensitive to the high intracellular iron content, which could lead to overload of iron in the body. At the same time, iron overload, heme deposition and oxidative damage in cells are increased. The methylation level of SCL40A1, HMOX2, TFRC, IREB2 and ACO1 is decreased, while the methylation level of FLVCR1 is increased, which may lead to the maintenance of intracellular iron homeostasis. FPN1 colocalizes with hephaestin and transports more Fe2+ from the cytosol to the extracellular environment (Chen et al., 2010), which reduces the level of iron in cells. HMOX2 cleaves a small amount of heme to form biliverdin and Fe2+ (Kutty et al., 1994). FLVCR1 weakens the ability of heme to move from the cytosol to the extracellular environment (Rey et al., 2008), which may lead to the ineffective decomposition of heme, as a result of which the body is attacked by heme toxicity, and heme deposition and oxidative damage are aggravated. ACO1 combines with a 4Fe–4S cluster as an aconitase (Philpott et al., 1994), isomerizing citrate to isocitrate (Dupuy et al., 2006) and reducing the iron content in the cell, and may promote the organism’s metabolism and reduce the damage caused by iron overload. NEDD8 is covalently attached to CUL, and forms an E3 ubiquitin ligase complex with FBXL5 and SKP1 (Salahudeen et al., 2009), and then activates its related E3 ubiquitin ligase activity. The protein complex of SKP1–CUL1–FBXL5 ubiquitin ligase is degraded by iron-related ubiquitination and IREB2 to decrease iron storage (Vashisht et al., 2009). These reactions are the protective responses to iron overload.

There is a complex relationship between DNA methylation and gene function during the iron absorption and transport of iron-overloaded Tibetans. ACO1, CYBRD1, FLVCR1, HFE, HMOX2, IREB2, NEDD8, SLC11A2, SLC40A1 and TFRC are key regulatory genes for iron overload in Tibetans, and their DNA methylation status may be the key functional factor for iron overload. Similarly, the expression of these genes is also different under iron overload. Our research is currently limited to the analysis of the methylation levels of the genome, and there is no measurement of actual gene expression levels. Targeted methylation detection of the population should also be expanded. The methylation of these genes may contribute to the observed significant changes in ferritin content of Tibetan people at high altitudes. The epigenetic mechanisms involved in the regulation of these genes and the genetic regions involved in iron metabolism need further study.

In this study, a high-resolution single-base methylation profile of a Tibetan iron overload genome was constructed by WGBS. Methylation differences between TH and TL groups at the genome level were detected, and 10 DMGs that may be related to the change of iron level in the Tibetan human body were identified. Genes involved in iron absorption and transportation will undergo methylation changes in the iron-overloaded human body. DNAm can influence the iron status of the human body by affecting the function of genes. Some changes in the level of gene methylation observed in our research may be the cause of iron overload or the result of iron overload. Therefore, the epigenetic molecular mechanism of iron homeostasis and metabolism needs to be further studied. We believe that this work not only provides new clues for elucidating the epigenetic mechanism of iron overload in the Tibetan population at high altitudes but also provides a new idea for the pathogenesis of iron metabolism.

MATERIALS AND METHODS

Sample preparation and DNA isolation

Six Tibetan men were divided into a high serum ferritin group (TH; n = 3) and a normal serum ferritin group (TL; n = 3) according to their serum ferritin content. Their age range was 30–72 years, and the average SF levels of the TH and TL groups were 1,372.3 and 180.9 ng/ml. Iron overload diagnostic criteria used European and American standards; after the elimination of active inflammation, liver disease, tumor, hemolysis, alcohol and other factors, serum ferritin > 1,000 ng/ml can be diagnosed as iron overload (Anonymous, 2011). The three criteria for inclusion in the group of six samples were as follows: 1) They are living at high altitudes above 2,500 meters. 2) They are all male (their average age is 49.65 ± 13.85). 3) Iron overload: serum ferritin > 1,000 ng/ml and iron normal: serum ferritin < 200 ng/ml. Nine criteria for exclusion were: 1) Hematological diseases such as anemia, polycythemia vera and myelodysplastic syndrome. 2) Patients with secondary hypertension. 3) Definitive patients with atherosclerotic cardiovascular disease, such as coronary heart disease, coronary stenting, coronary bypass surgery, stroke, peripheral vascular disease, revascularization, carotid endarterectomy and carotid stenting. 4) Patients with heart failure with reduced systolic function (EF < 50%). 5) Creatinine > 3 times normal, glomerulonephritis, renal insufficiency, ALT > 3 times normal and patients with liver dysfunction. 6) Chronic obstructive pulmonary disease and sleep apnea syndrome. 7) Congenital heart disease and heart surgery. 8) History of venous thromboembolism. 9) Use of statins. All six subjects have no kinship.

Two milliliters of peripheral venous blood was collected from each person, and transported on ice. The blood samples were stored at –80 ℃. Total genomic DNA was extracted from the blood using a TIANamp Genomic DNA Kit (Tiangen Biotech, Beijing, China) according to the manufacturer’s instructions. Agarose gel electrophoresis was used to monitor the integrity and purity of DNA.

All procedures involving human participants were performed in accordance with the ethical standards of the Ethics Committee of the Hospital of Chengdu Office of People’s Government of Tibetan Autonomous Region in Chengdu (Approval ID: 2018, 64). Informed consent was obtained from all individual participants included in the study.

Library construction

Genomic DNA was extracted from the samples and used to construct a library. After the sample was quantified, a certain proportion of negative control (lambda DNA) was added and the genomic DNA was randomly sheared to 200–300 bp and used to construct a library as shown in Supplementary Fig. S5. Following bisulfite treatment using an EZ DNA Methylation-Gold Kit (Zymo Research), the unmethylated C becomes U (it becomes T after PCR amplification), while the methylated C remains unchanged. PCR amplification was then performed to obtain the final DNA library, which was quantified and paired-end-sequenced on the Illumina HiSeq 2500 platform (OE Biotech, Shanghai, China) and called the original read.

Quality control

To eliminate the influence of sequencing errors on the results, it is necessary to control the quality of the original sequencing data so as to obtain clear readings. The preprocessing software was NGSQCtoolkit (Dai et al., 2010), and the quality filtering standards were as follows: (1) Filter out readings with high-quality bases (Q20) less than 70%, that is, more than 70% of the bases should be greater than or equal to 20. (2) Remove the low-quality base at the 3′ end, and the sequence of the base having a mass value of less than 20 will be removed from the 3′ end. (3) Reads containing n (non-ATCG) are removed.

Bismark software (Hansen et al., 2012) and the reference genome GRCh37.p13 (https://www.ncbi.nlm.nih.gov/assembly/GCF_000001405.25/) were used to compare and analyze the methylation data. Bismark uses the sodium bisulfite conversion algorithm to convert all C bases in the reads and the reference genome into T bases, and then compares them to locate the reads correctly and identify the methylation sites. The process of data analysis is shown in Supplementary Fig. S6. The sequencing data were deposited in the National Center for Biotechnology Information (NCBI), with the accession number PRJNA559755.

Estimating methylation levels

To ensure accurate detection of methylation sites, sites with ≥ 5 covered reads were regarded as truly detected C sites (CoveredC). In addition, because each site includes methylated reads and unmethylated reads, and the detection probability of methylated reads should follow a binomial distribution, a binomial distribution test was performed on the number of methylated and unmethylated C reads at each site to identify whether the site was a real methylated site. Assuming that the number of methylated Cs at a certain site is x, the read coverage of this site is n, and the methylation conversion rate is p, then it is necessary to verify whether the occurrence of x methylated C is reliable under the conditions that the probability is p and the sequencing depth is n. After FDR detection, the sites with q-value ≤ 0.05 are regarded as true methylation sites. Methylation level for a C site represents the segment of methylated C sites, and the ratio of methylated cytosines is equal to the number of methylated cytosine sites divided by the detected cytosine sites. In this study, Cm was defined as the number of supporting methylation reads at a single C site, while Cu was defined as the number of supporting unmethylation reads at a single C site. The methylation level of a site and region is shown in the equation below. Additionally, we calculated DNA methylation levels at three different scales: chromosome, gene region and genomic feature.   

Methylation   Level   of   C   site= C m C m + C u ×100%

Identification of differentially methylated regions (DMRs)

DMR detection used methylKit (Akalin et al., 2012) software. First, the methylation level data were entered in the sample genome. The sliding window and sliding step length were set to 1,000 bp. Multi-sample methylation data were combined for DMR analysis. The logistic regression model analyzes the DMR and predicted the corresponding DMR differences between the groups. Differences in methylation levels and statistical significance of all windows between all groups were determined, using cutoff values of difference ≥ 10% and q-value ≤ 0.05, respectively.

Functional enrichment analysis

If there is a DMR in a gene and its upstream 2 kb, the gene is regarded as a differentially methylated gene (DMG). All DMGs were annotated by Gene Ontology (GO, http://www.geneontology.org/) and Kyoto Encyclopedia of Genes and Genomes (KEGG). GO analysis was performed with the GOseq R software package and corrected by the Benjamini–Hochberg method. The Kobas website (Xie et al., 2011) was used to test the statistically significant enrichment of DMGs in the KEGG pathways. The STRING database (Szklarczyk et al., 2019) (http://string-db.org/) was used to analyze the interaction networks and functional enrichment of the selected DMGs, and the network data were imported into Cytoscape (Shannon et al., 2003) software to output the results. The Reactome website (Fabregat et al., 2018) (https://reactome.org/) was used to examine the enrichment of selected DMGs involved in iron-related pathways.

Biological information analysis

Iron overload-related data sets GSE160727 (Votavova et al., 2021), GSE7357 (Coppin et al., 2007), GSE9726 (Rodriguez et al., 2007) and GSE10421 (Kautz et al., 2008) were downloaded from the Gene Expression Omnibus (GEO) database, and DEGs were analyzed by GEO2R. The specific samples and groups used in each data set are listed in Supplementary Table S2.

ACKNOWLEDGMENTS

This study was funded by the Sichuan Science and Technology Program (Grant No. 2019YJ0642).

REFERENCES
 
© 2022 The Author(s).

This is an open access article distributed under the terms of the Creative Commons BY 4.0 International (Attribution) License (https://creativecommons.org/licenses/by/4.0/legalcode), which permits the unrestricted distribution, reproduction and use of the article provided the original source and authors are credited.
https://creativecommons.org/licenses/by/4.0/legalcode
feedback
Top