Zheng Zhou1,2,3†
Jinjin Liu1,2,3†Yun Chen1,2,3Bingxuan Ren1,2,3Siyuan Wan4Yao Chen1,2,3Yanhong He1,2,3
Qiuyang Wei5,6Haiyan Gao1,2,3,7Lixiang Liu1,2,3*
Hongmei Shen1,2,3*- 1Disorders Control, Centre for Endemic Disease Control, Chinese Centre for Disease Control and Prevention, Harbin Medical University, Harbin, Heilongjiang, China
- 2Key Laboratory of Etiology and Epidemiology, National Health Commission & Education Bureau of Heilongjiang Province, Harbin Medical University, Harbin, Heilongjiang, China
- 3Heilongjiang Provincial Key Laboratory of Trace Elements and Human Health, Harbin Medical University, Harbin, Heilongjiang, China
- 4Department of Preventive Medicine, Qiqihar Medical University, Qiqihar, Heilongjiang, China
- 5First Clinical Medical Department, Heilongjiang University of Chinese Medicine, Harbin, Heilongjiang, China
- 6Second Department of Endocrinology, The First Affiliated Hospital of Heilongjiang University of Chinese Medicine, Harbin, China
- 7Department of Clinical Laboratory, The Sixth Affiliated Hospital of Harbin Medical University, Harbin, China
Background: Hashimoto thyroiditis (HT), a prevalent autoimmune disorder, is not yet thoroughly understood, especially when it comes to the influence of epigenetics in its pathogenesis. The primary goal of this research was to probe the DNAm profile across the genome in the whole blood derived from patients suffering from HT.
Method: Using the Illumina 850K BeadChip, we conducted a genome-wide DNAm assessment on 10 matched pairs of HT sufferers and healthy individuals. Genes with differential methylation (DMGs) were identified and underwent functional annotation via the databases of Gene Ontology and Kyoto Encyclopedia of Genes and Genomes. The transcriptional significance of potential epigenetic biomarker genes was corroborated through qRT-PCR.
Results: The DNAm profiling across the genome indicated an overall reduction in methylation in HT subjects in comparison with their healthy counterparts. We detected 283 DMPs (adjusted P < 0.05 and |Δβ| > 0.1), among which 152 exhibited hypomethylation and 131 demonstrated hypermethylation. Further analysis exposed a noteworthy concentration of hypermethylated DMPs in the 3´UTR, North Shore, and CpG islands, while there was a significant decrease in the Open Sea (all P < 0.001). The 283 DMPs were broadly distributed from chromosome 1 to 22, with chromosome 6 harboring the most DMPs (n = 51) and chromosome 12 carrying the most DMGs (n = 15). The SLFN12 gene, which presented with extreme hypomethylation in its promoter DMPs among HT patients, was identified as the epigenetic marker gene. Consequently, the SLFN12 mRNA expression was markedly upregulated in HT, displaying a negative relationship with its methylation levels. The area under curve (AUC) value for the SLFN12 gene among HT patients was 0.85 (sensitivity: 0.7, specificity: 0.7), a significant difference compared with healthy controls. The methylation levels of all DMPs in SLFN12 gene were negatively correlated with TSH and one CpG site (cg24470734) was positively assocciated with FT4.
Conclusion: This investigation presents an initial comprehensive DNAm blueprint for individuals with HT, which permits clear differentiation between HT subjects and normal controls through an epigenetic lens. The SLFN12 gene plays a pivotal role in the onset of HT, suggesting that the methylation status of this gene could serve as a potential epigenetic indicator for HT.
1 Introduction
Currently, Hashimoto thyroiditis (HT) is recognized as the predominant autoimmune disorder of the endocrine system, standing as the leading cause of hypothyroidism. It’s distinguished by the presence of antibodies that target thyroid antigens and the lymphocytic penetration of the thyroid gland (1). The onset of HT can be instigated by both genetic predisposing factors and influences from the environment. Twin studies have allowed us to estimate that genetic elements account for 70%-80% of the susceptibility to HT (2, 3). As a result, it’s of utmost importance to pinpoint the genes that possess a genetic susceptibility impacting the initiation and progression of HT through whole-genome sequencing.
In recent times, the potential implications of epigenetics in thyroid-related disorders have begun to gain significant interest. Epigenetics, defined as inheritable modifications that impact gene expression without inducing alterations to the DNA sequence, serves to elucidate the possible correlations among genetic elements, environmental influencers, and the development of HT (4). Of the various epigenetic mechanisms, DNA methylation (DNAm) has been the subject of extensive investigation. DNAm entails the integration of a methyl group into a cytosine unit in a cytosine-phosphate-guanine (CpG) dinucleotide sequence, enabling further regulation of gene expression by impeding the attachment of transcription factors to DNA (5). Multiple research efforts have delved into DNAm segments of HT susceptibility genes. For instance, it was discovered that DNAm in the ICAM-1 promoter area was diminished in HT patients and exhibited a negative correlation with thyrotropic immunoglobulin concentration, which ultimately influenced the adhesion of cytotoxic T cells to their targets (6). In a separate study, significantly elevated levels of DNAm were identified in the promoter region of PTPN22 in young patients with HT, relative to controls; this verified the role of PTPN22 as an epigenetically determined susceptibility gene for thyroid autoimmunity (7).
To conclude, it’s evident that alterations in DNAm hold a significant correlation with HT. While many studies have explored the potential changes in methylation levels of susceptibility genes amongst HT patients, to our understanding, no prior research has published a genome-wide DNAm study specifically involving HT patients sourced from the general population. Hence, in the present research, we conducted a comprehensive DNAm investigation using whole blood specimens from both HT patients and healthy individuals. Our study’s objective was to identify the crucial differentially methylated genes (DMGs) through the means of bioinformatics analysis and subsequently pinpoint the epigenetic biomarker gene linked with HT, thereby furnishing pivotal molecular insights into HT’s pathogenesis.
2 Materials and methods
2.1 Study participants
In the Shandong Province of China, June 2019 saw the involvement of a collective group of 30 paired matches, comprising of HT-affected subjects and healthy control individuals. The diagnosis of HT was effectively established by an experienced endocrinologist, using two key indicators: 1. The existence of serum thyroid peroxidase antibodies (TPOAb), thyroglobulin antibodies (TgAb), or a combination of both; 2. Apparent signs of goiter, changes in echogenicity, or multiple hypoechoic regions as determined by thyroid ultrasound tests. On the other hand, the criteria for selection of control participants included: 1. Healthy individuals of similar sex, age, locale, and BMI to the HT-affected subjects; 2. Absence of personal or familial histories of autoimmune disorders or other thyroid complications, no chronic or acute illnesses, no record of long-term thyroid-related medication or hormone use, and non-pregnant status; 3. A clean bill of health in terms of goiter; negative for both TgAb and TPOAb; normal thyroid function indicators, and no irregularities spotted in thyroid ultrasound examinations. Blood samples were drawn from all participants and kept in a -80°C storage for later examinations. As shown in Figure 1A, out of the 30 pairs, 10 were selected for the genome-wide DNAm analysis while the entire group was used for mRNA expression examination of the epigenetic biomarker gene. The research upheld the standards set in the Declaration of Helsinki and received approval from the Ethics Review Committee of Harbin Medical University (No. hrbmuecdc20200320). The participants provided their consent by signing the necessary forms.
Figure 1 A series of depictions pertaining to the 850K BeadChip, the density plot, and the PCA plot. (A) Presents the experimental procedure. (B) Illustrates the quality assurance and data preprocessing measures adopted for the 850K BeadChip. (C) Depicts a density plot. Beta values express the levels of methylation in both HT patients and healthy controls. (D) Showcases a PCA plot grounded on the 719,455 CpG sites present in the 850K BeadChip, with no clear demarcation discernible between HT patients and the controls. (E) Reveals a PCA plot grounded on 283 differentially methylated positions (DMPs). The plot shows clear clustering differences between HT patients and healthy controls.
2.2 Assessment of thyroid function and detection of thyroid antibody
Measurements of free triiodothyronine (FT3), free thyroxine (FT4), thyroid stimulating hormone (TSH), TgAb, and TPOAb were obtained through the application of chemiluminescence immunoassay techniques, a service facilitated by Siemens Inc., Tarrytown, NY, USA. The defined standard values for these measures were set as follows: for FT3, an acceptable range was considered to be 3.1 pmol/L to 6.8 pmol/L; FT4 standard value resided between 11.5 pmol/L and 22.7 pmol/L; TSH was designated as normal if it fell between 0.27 mIU/L and 4.20 mIU/L. To detect thyroid antibodies, benchmarks were set with the normal threshold for TgAb and TPOAb each maintained within the bounds of 0–60 U/mL.
2.3 Genome-wide DNAm measurement
The TIANGEN Extraction Kit (TIANGEN, Beijing, China) was employed for the procedure of genome-wide DNA extraction. The purity and concentration of the harvested DNA were subsequently confirmed with the assistance of the Nanodrop 2000 instrument (Thermo Fisher, Waltham, USA). Following this, we processed 500 ng of DNA from each specimen for bisulfite conversion, a procedure facilitated by the EZ DNA Methylation Kits (Zymo Research, California, USA). Post conversion, the resulting products were processed further using the Illumina Infinium Human Methylation 850K BeadChip (Illumina Inc, California, USA), in strict accordance with the manufacturer’s directives.
2.4 Genome-wide DNAm BeadChip analysis
The BeadChip was scanned using an Illumina iScan instrument. The raw intensity data, also known as IDAT, was loaded into Rstudio for further analysis. Using the Chip Analysis Methylation Pipeline (ChAMP) package (2.14.0), the data underwent preprocessing, normalization, and group comparison processes (8). To identify confounding factors, singular value decomposition analysis was executed. To elaborate, the minfi package (1.30.0) was used for preprocessing the raw data, and technical variations were normalized with SWAN (9). Certain probes were excluded from the analysis based on the following specific criteria: 1. Any CpG probes that had a detection P ≥ 0.01 in one or more samples; 2. probes with fewer than three beads in 5% of the samples; 3. non-CpG probes; 4. cross-reactive probes and probes associated with single nucleotide polymorphisms (SNPs); 5. probes located on the X or Y chromosomes.
2.5 Differential methylation position analysis and the distribution analysis of DMPs
A β value was utilized to signify the methylation extent for each CpG site, with the scale ranging from 0 (0% methylation) to 1 (100% methylation of a specific CpG dinucleotide). The Δβ was computed by deducting the β value of the control group from that of the HT cases. A CpG site with Δβ > 0 was deemed as hypermethylated, whereas hypomethylation was assigned when Δβ < 0. Differential methylated positions (DMPs) were discerned by contrasting mean β values in the HT group with those in the normal group for a given CpG site through ChAMP analysis (10). We established DMP criteria by setting a P value of less than 0.05 and a β value discrepancy between groups greater than 0.1 (|Δβ| > 0.1). The DMPs were then segregated into distinct groups based on various gene regions: TSS1500, TSS200, 5’UTR, 1stExon, gene body, 3’UTR, and intergenic region (IGR). In accordance with different CpG island regions, the DMPs were also divided into: North Shelf, North Shore, CpG Island, South Shore, South Shelf, and Open Sea, wherein shores were 0–2 kb away from a CpG island and shelves were 2–4 kb distant. The entire array data processing and analysis were executed in Rstudio.
2.6 Gene ontology and Kyoto encyclopedia of genes and genomes enrichment
Every identified DMP was matched to its respective DMGs. To delve into the prospective biological roles of these DMGs, we employed the R package clusterProfiler for conducting GO analysis and KEGG pathway examination. To identify the significantly enriched GO terms and pathways, we set a P value of less than 0.05 as the cut-off criterion.
2.7 Selection of epigenetic biomarker gene
Subsequently, we embarked on a mission to further identify potential epigenetic biomarker genes linked to HT. We adopted a strategic set of screening criteria for DMGs previously screened through Illumina 850K BeadChip. The rules for potential epigenetic biomarker genes included: 1. DMGs situated within promoter regions (namely, TSS1500, TSS200, 5′UTR, and 1stExon); 2. DMGs with strikingly high or low methylation levels (with a β value exceeding 0.7 or under 0.3); 3. DMGs possessing at least three corresponding DMPs.
2.8 Quantitative real-time PCR
RNA extraction was performed from total blood samples utilizing RNAiso Plus (Takara, Dalian, China). The concentration of RNA was ascertained using the NanoDrop 2000 (Thermo Fisher, Waltham, USA). RNA quality was deemed satisfactory with an OD 260/280 ratio ranging between 1.8 and 2.0. Each specimen was further processed using a QuantStudio™5 Real-Time PCR system (Thermo Fisher, Waltham, USA). The amplification reaction was set up with the following constituents: 5.0 μL SYBR Green (Roche, Basel, Switzerland), 1.0 μL cDNA, 0.5 μL of each primer (upstream and downstream), and 3.0 μL ddH2O. The PCR protocol was as follows: 1. initial denaturation at 95°C for 10 min; 2. forty cycles of 95°C for 15s and 60°C for 1 min; 3. melt curve stage at 95°C for 15s, 60°C for 1 min, and 95°C for 15s. The primers were designed as such: SLFN12-F: 5′-TGG CCT CTT TTG GAA TGG CA-3′; SLFN12-R: 5′-CAA GCA GCC CAG ATC CAC AGA CC-3′; β-actin-F: 5′-CCT TCC TGG GCA TGG AGT CCT G-3′; β-actin-R: 5′-GGA GCA ATG ATC TTG ATC TTC-3′. The 2-ΔΔCt method facilitated the analysis of mRNA levels of the epigenetic biomarker gene.
2.9 Statistical analysis
Data from the 850K BeadChip, as well as GO and KEGG, underwent statistical evaluation in the RStudio software (Version 1.2.1335) (http://www.rstudio.com/) within an R coding environment (version 3.6.0) (https://www.R-project.org). Participant demographic details and thyroid functionality data were inputted into Microsoft Office Excel 2016 and analyzed statistically using the SPSS 23.0 software. The Kolmogorov–Smirnov approach was used to test for normal data distribution. Data adhering to normal distribution are displayed as mean ± standard deviation (mean ± SD), and mean discrepancies were examined using student’s t-test. Non-normally distributed data are expressed through the median along with the 25th and 75th percentiles. Comparisons among various groups were executed using the Mann–Whitney U test. The chi-square (χ2) test was applied for comparing rates among groups with categorical results. The associations between relative mRNA expression and the methylation levels of SLFN12 were scrutinized via Pearson correlation analysis. The diagnostic value and the cutoff point were ascertained through receiver operating characteristic (ROC) curve analysis. Correlation between the DNAm levels of 5 DMPs in SLFN12 and variables (age, FT3, FT4, and TSH) was assessed through Spearman’s rank or Pearson correlation analyses. A p-value of less than 0.05 was taken to be statistically significant.
3 Results
3.1 Global DNAm profile of HT patients
Table 1 delineates the demographic attributes and thyroid functions of 10 pairs consisting of Hashimoto’s thyroiditis (HT) patients and healthy control subjects analyzed using the 850K BeadChip. All the subjects classified under the HT group were in compliance with the HT diagnostic criteria. Supplementary Table 1 provides the demographic details for 30 pairs of HT patients and control subjects and their thyroid functions for validation. The BeadChip detected a total of 865,918 CpG sites for future manipulation. Applying the selection process as described in our materials and methods section, a total of 719,455 CpG sites were retained for subsequent scrutiny (Figure 1B). Following this, normalization of data and quality control was performed on these 719,455 CpG sites. As observable from the density plot (Figure 1C), the methylation level distribution of HT patients and their healthy counterparts showed considerable similarity after the aforementioned procedure. We then executed principal component analysis (PCA) on the 719,455 CpG sites (Figure 1D). The outcome demonstrated no clear demarcation between the HT patients and healthy control subjects.
Table 1 The basic information of HT patients and healthy controls in 850K.
Next, we scrutinized the DNAm levels in various genomic locations in both HT patients and their healthy counterparts (Figure 2). First and foremost, the patterns of genome-wide DNAm spanning 719,455 CpG sites exhibited a bimodal distribution. Most CpG sites, in both HT patients and control groups, showed high (β > 0.8) or low (β < 0.2) levels of DNAm. In a general sense, the methylation levels in HT were considerably lower compared to the healthy controls, and this difference was statistically significant (P < 0.001). In the second place, a substantial number of unmethylated CpG sites (β < 0.1) were found distributed within the promoter and 5’UTR region. The promoter region in HT patients displayed higher methylation levels compared to the healthy control group, and this discrepancy was statistically significant (P = 0.006). Last but not least, both IGR and gene body displayed a downward shift in median methylation in HT patients when juxtaposed with healthy controls (both P < 0.001). These outcomes point to the conclusion that the DNAm levels within different gene regions undergo significant changes in HT patients.
Figure 2 Comparison of DNA methylation levels across various genomic locations in HT patients versus healthy controls. The violin plots reflect DNA methylation at a genomic scale (n = 719,455), within promoters (n = 142,552), spanning intergenic regions (n = 199,991), across gene body regions (n = 270,613), within 3′UTRs (n = 18,023), and 5′UTRs (n = 61,970). The embedded box plots in each violin plot signify the interquartile range, and the red lines delineate a significant statistical deviation between HT patients and the healthy controls. For all the cases, the y-axis denotes the methylation level (β value) on a scale from 0 to 1.
3.2 DMPs in HT patients
In contrast to the 10 healthy controls, HT patients revealed a total of 283 differential methylation positions (DMPs) (adjusted P < 0.05, |Δβ| > 0.1). A volcano plot was used to illustrate the DMPs, where it was shown that of the 283 DMPs, 152 (53.7%) demonstrated hypomethylation, and 131 (46.3%) exhibited hypermethylation (Figure 3A). The top 10 DMPs showing hypermethylation and hypomethylation were ranked according to |Δβ| and are presented in Supplementary Table 2. Utilizing the 283 DMPs, the PCA distinguished two unique clusters separating the HT patients from the healthy controls (Figure 1E). The unsupervised hierarchical cluster heatmap of 283 DMPs distinctly segregated HT patients and healthy controls (Figure 3B). Each row within the heatmap denotes a CpG probe while each column represents an independent sample. Given the known relationship between the specific genomic location of DMPs and their effect, we then endeavored to determine the genomic distribution of DMPs. The distribution of hypomethylated DMPs in the 5’UTR was significantly lower (2.63%) compared to the probes in 850K BeadChip, whereas the distribution of hypermethylated DMPs in the 3’UTR was significantly increased (6.10%) in comparison to hypomethylated DMPs (both P < 0.001) (Figure 3C left panel). As for the CpG island regions, the hypermethylated DMPs were significantly overrepresented in North Shore (16.03%) but notably lower in Open Sea (47.32%) when compared to the probes in 850K BeadChip and hypomethylated DMPs (all P < 0.001). Furthermore, in comparison to hypomethylated DMPs, hypermethylated DMPs were significantly overrepresented in CpG islands (23.7%) (P < 0.001) (Figure 3C right panel). The 283 DMPs in HT patients were broadly distributed across chromosomes 1 through 22 (Figure 4). Out of the 22 chromosomes, chromosome 6 had the most DMPs (n = 30) and hypomethylated DMPs (n = 22); chromosome 11 had the highest number of hypermethylated DMPs (n = 15). Additionally, chromosome 12 housed the largest number of differentially methylated genes (DMGs) (15 genes). All hypermethylated and hypomethylated DMGs on chromosomes 1 through 22 are listed in Supplementary Table 3.
Figure 3 The dispersion of differentially methylated probes (DMPs) in HT patients. (A) A Volcano plot visualizing probe-level methylation in HT patients compared to healthy controls. The diagram presents the relationship between Δβ values (along the x-axis) and adjusted P values (-log10 transformed adjusted P values) on the y-axis. Each dot is representative of an individual probe. The demarcations for adjusted P = 0.05 and |Δβ| = 0.1 are indicated by horizontal and vertical dashed lines, respectively. Red and blue dots symbolize hyper- and hypomethylated DMPs, correspondingly. (B) Heatmap of the 283 DMPs distinguishing HT patients from healthy controls. Each row stands for a unique DMP and each column symbolizes a participant. The diagram consists of 283 rows (indicative of 283 DMPs) and 20 columns (denoting 10 HT patients versus 10 healthy controls). On the top, red and blue bars represent HT patients and healthy controls, respectively. Color transition from red to blue signifies a shift from hypermethylation to hypomethylation levels. (C) Bar charts illustrating the distribution of DMPs with respect to gene region (left panel) and CpG island region (right panel). The total probes in the 850K BeadChip available for analysis (grey; n = 719,455), hypermethylated DMPs (red; n = 131), and hypomethylated DMPs (blue; n = 152) as per gene region and CpG island region. P values were computed by Chi squared tests. ‘a’ suggests the difference is statistically significant compared with probes in the 850K BeadChip; ‘b’ indicates that the difference is statistically significant compared to hypomethylated DMPs.
Figure 4 Circos diagram of differentially methylated probes (DMPs). The circos plot delineates the distribution of 283 DMPs across 22 chromosomes; the exterior layer depicts the physical location of each chromosome; the intermediate layer illustrates chromosome segments; the innermost layer signifies the methylation site. Here, red color corresponds to hypermethylation DMPs, while blue represents hypomethylation DMPs.
3.3 GO functional analysis and KEGG pathway analysis
Among the 283 identified DMPs, 147 had biological data available and corresponded to differentially methylated genes (DMGs). A total of 246 biological process (BP) terms were enriched by these DMGs, of which 59 demonstrated significant enrichment (all P < 0.010) (Supplementary Table 4). For cellular component (CC) terms, 42 were enriched by all DMGs with 10 being significant (all P < 0.010) (Supplementary Table 5). When considering molecular function (MF), 37 terms were enriched by all DMGs and 13 showed significant enrichment (all P < 0.010) (Supplementary Table 6). The top 10 terms enriched by BP, CC, and MF are demonstrated in Figures 5A, C, E respectively. Concurrently, Figures 5B, D, F highlight all the DMGs enriched in the top 5 terms of BP, CC, and MF, respectively. A subset of 18 DMGs emerged in the top 5 BP terms, consisting of 8 hypermethylated and 10 hypomethylated genes. In the top 5 CC terms, a total of 20 DMGs were enriched, composed of 7 hypermethylated and 13 hypomethylated genes. Similarly, in the top 5 MF terms, 16 DMGs were discovered including 3 hypermethylated and 13 hypomethylated genes. A total of 85 pathways were enriched by all genes, with 17 pathways showing significant DMGs enrichment (all P < 0.050) (Table 2). Among these pathways, cell adhesion molecules (hsa04514), Th17 cell differentiation (hsa04659), and calcium signaling pathway (hsa04020) were noted for their relevance to autoimmune disease. The top 10 pathways significantly enriched by KEGG of the 147 DMGs are shown in Figure 6A, while Figure 6B highlights the 20 DMGs enriched in these top 10 signaling pathways.
Figure 5 HT-related differentially methylated genes (DMGs) explored via Gene Ontology (GO) analysis. Panels (A, C, E) present bubble charts, each showcasing the top 10 enriched terms under biological process (BP), cellular component (CC), and molecular function (MF) respectively. The vertical axis represents the specific GO pathway’s name, while the P-value is plotted on the horizontal axis. The bubble size corresponds to the count of differentially methylated genes participating in that pathway. The color intensity of each bubble relates to the –log10 P-value. Panels (B, D, F) feature the GO Chord plots illustrating the association between DMGs and their relevant BP, CC, and MF terms. Various color blocks denote distinct GO terms with the symbols of targeted genes displayed on the chart’s left. The colored lines on the right establish the gene involvement under each GO term. Red blocks denote DMGs with hypermethylation, and blue blocks correspond to hypomethylated DMGs.
Table 2 17 significant pathways of KEGG analysis.
Figure 6 Pathway analysis of HT-related differentially methylated genes (DMGs) based on Kyoto Encyclopedia of Genes and Genomes (KEGG). Panel (A) depicts a bubble chart, illuminating the top 10 KEGG-enriched pathways. Panel (B) shows that among the top 10 signaling pathways, 20 DMGs are notably enriched. The dot size is indicative of the count of DMGs enriched.
3.4 Hypomethylation of SLFN12 gene, an epigenetic biomarker gene for HT
The procedure for screening DMGs potentially linked to the onset of HT, with the aim of identifying epigenetic biomarker genes, is presented in Figure 7A. Initially, our focus was confined to the DMGs found on the promoter. Sixty-five DMPs resided in the promoter region, corresponding to a total of 46 DMGs (Refer Supplementary Table 7). To further refine the pool of DMGs, we took into account the methylation level of each DMP, selecting only those DMPs with an absolute high methylation level (β > 0.7) or an absolute low methylation level (β < 0.3). This strategy resulted in 30 DMPs that corresponded to 22 DMGs (Refer Supplementary Table 8). Our final criterion was to choose a DMG that encompasses at least 3 DMPs, which led us to select the SLFN12 gene as an epigenetic biomarker. SLFN12 hosted 5 DMPs, distributed across TSS1500, TSS200, and the 1st Exon, and all exhibited methylation levels below 0.2 in HT cases (Figure 7B). Post identification of the epigenetic biomarker gene, we assessed the mRNA levels of the SLFN12 gene using whole blood from 30 matched pairs of HT patients and healthy controls. The results revealed a significant elevation of SLFN12 in HT patients (3.882 ± 1.355 vs. 1.026 ± 0.083, P < 0.001) (Figure 7C). Figures 7D–H illustrate the correlation between DNAm and SLFN12 gene mRNA expression in 10 matched pairs of HT patients and healthy controls. The data suggest a negative correlation between the DNAm levels of SLFN12 and its expression across all 5 DMPs in the 850K BeadChip (all P < 0.05). The area under the curve (AUC) value of 0.85 (sensitivity: 0.7, specificity: 0.7) for the 5 DMPs mapped to SLFN12, a potential biomarker, among HT patients compared to healthy controls, indicates the excellence of our model (Figure 7I). ROC curves were also constructed for each DMP and the AUC for all DMPs exceeded 0.85 (Refer Supplementary Table 9).
Figure 7 Unveiling the epigenetic marker gene in HT. Panel (A) represents a schematic demonstrating the number of DMPs discerned in each stage during the screening progression of the epigenetic marker gene, with SLFN12 singled out for verification, as denoted by the red dotted rectangle. Panel (B) illustrates methylation states at five distinct locations of the SLFN12 gene among HT patients versus healthy controls (mean ± SD). Panel (C) presents SLFN12 mRNA expression levels, where ** denotes P < 0.001. Panels (D–H) exhibit the correlation established between DNA methylation and relative mRNA expression across five specific SLFN12 gene sites. Panel (I) presents the ROC curve tailored for SLFN12.
3.5 Correlations between SLFN12 DNA methylation levels and age, as well as thyroid function
Figure 8 displays negative correlations between the DNAm levels of two DMPs situated in the SLFN12 gene (cg03251655 and cg24470734) and age (both P < 0.050). Furthermore, all DMPs in the SLFN12 gene show negative correlations with TSH levels (all P < 0.050). Additionally, one DMP (cg24470734) in the SLFN12 gene demonstrates a positive correlation with FT4 levels (P < 0.050).
Figure 8 Heat map of correlations among SLFN12 DNA methylation levels, age, and thyroid function. * denotes P < 0.050; ** denotes P < 0.010.
4 Discussion
Employing the Epigenome-wide association study (EWAS) is a commonly adopted approach for unearthing biomarkers within populations and for elucidating the molecular underpinnings of disease susceptibility. The 850K BeadChip serves as one of the instruments used in this methodology (11). In recent years, epigenetics has provided a potential link between genetics and autoimmune thyroid disease phenotype with some studies reporting a relationship between changes in DNAm patterns of some genes and HT (12). However, to the best of our knowledge, this is the first comprehensive study providing distinct DNAm profiles of HT patients screened from the general population based on EWAS. The 850K BeadChip we used was able to cover more than 850,000 CpGs, which greatly has advanced our understanding of the association between DNAm and HT, and demonstrated the importance of DNAm modification in the development of HT.
Prior research on a genome-wide scale has suggested alterations in the DNAm profiles in autoimmune disease patients (13). Several investigations have demonstrated that a characteristic trait of autoimmune diseases, including systemic lupus erythematosus, multiple sclerosis, and chronic spontaneous urticaria, is global hypomethylation (14). Our study findings align with these conclusions. Specifically, we found that in relation to the genome-wide DNAm levels (719,455 CpGs), the levels in patients with HT were significantly lower than in healthy individuals. Additionally, our findings indicated that despite a decrease in DNAm in other gene regions, the DNAm was heightened in promoter regions in HT patients, a result that corroborates many previous research conclusions (15). From these observations, it can be concluded that DNAm anomalies occur in the entire genome of HT patients. This is a key consideration for future HT research and may have practical implications in a clinical setting.
Initially, at the level of DMPs, our study observed distinct methylation patterns in HT, consisting of 283 DMPs that significantly deviated from healthy controls. Heatmap investigations and PCA results suggest that these DMPs serve as distinguishing features between HT patients and healthy controls, offering valuable theoretical basis for the continued exploration of HT and epigenetics. Subsequently, past research has demonstrated that the distribution of DMPs within the promoter region varies across diseases. In autoimmune diseases, a predominant characteristic is hypomethylated sites in the promoter region DMPs (16, 17), a finding that is reinforced by our research. Contrastingly, in carcinoma, hypermethylated sites are more prevalent within promoter region DMPs (18, 19). This indicates that the role of hypomethylation sites in the promoter region in autoimmune diseases’ pathogenesis warrants further investigation. Moreover, existing literature reveals a connection between chromosomal abnormalities and DNAm (19), implying that defining the distribution and properties of DMPs across various chromosomes can enhance disease comprehension. Our study discovered that the 283 DMPs in HT patients were extensively dispersed across autosomal chromosomes, with their count varying greatly across different chromosomes. The maximum number of DMPs (30 sites) were located on chromosome 6, and the most DMGs (15 genes) were found on chromosome 12. Interestingly, both chromosome 6 and chromosome 12 contained the DMG related to autoimmune diseases. On chromosome 6, the HLA-DPB1 gene, a crucial gene in T lymphocyte activation and antigen presentation, is listed among the HT susceptibility genes (20). Concurrently, the C1R gene on chromosome 12 has been linked to the development of systemic autoimmune diseases (21, 22). Therefore, the distribution of DMPs across the chromosomes of HT patients is neither random nor uniform, suggesting that chromosomes with multiple DMPs or DMGs should be the focus in subsequent HT studies.
In our KEGG pathway analysis, we found that pathways involving cell adhesion molecules (hsa04514), Th17 cell differentiation (hsa04659), and the calcium signaling pathway (hsa04020) warrant careful consideration, given their established links with autoimmune diseases. Firstly, pathways associated with cell adhesion molecules are key to the adherence, rolling, and migration of leukocytes in inflamed tissues. Cell adhesion molecules specifically influence the migration of Th17 to inflammation sites, which effectively modulates immune response (23, 24). Secondly, recent findings underscore the critical role of Th17 cells, and their cellular and secretory components in the pathogenesis and progression of HT. Numerous studies have reported an elevated proportion of Th17 cells in HT patients (25, 26). In this context, the enrichment of the SMAD3 gene in this pathway is noteworthy. Known for its vital role in autoimmunity (27, 28), SMAD3 can control the activation of the TGF-β receptor, thereby influencing the differentiation of T cells into Th17 cells, and ultimately impacting the immune response. Thirdly, recent research has highlighted the crucial role of calcium signaling pathways in autoimmune diseases and inflammatory reactions. Kaufmann et al. demonstrated that calcium influx is a critical modulator of mitochondrial function and oxidative stress in Th17 cell-mediated multiorgan autoimmune disease and inflammation (29). Further, it has been shown that in several autoimmune diseases, B cell receptor calcium signaling exhibits defects, rendering B cells pathogenic and leading to the development of autoimmune diseases (30). While these pathways have been linked to autoimmune diseases, previous studies have not explored potential abnormalities in the methylation patterns of these pathways. Our research suggests that alterations in the methylation levels of key genes within these pivotal pathways might be implicated in HT pathogenesis.
Utilizing the state-of-the-art high-methylation BeadChip, we were able to filter a significant number of DMGs. The primary objective of our research was to pinpoint the crucial epigenetic biomarker gene, utilizing the data from the 850K BeadChip. Given the vast number of DMGs and the lack of a universally agreed approach to select epigenetic biomarker genes, we first excluded all DMPs that were not within the promoter region, given the ongoing debate over the influence of methylation sites outside the promoter region on gene expression (31). Further, DMPs with DNAm levels exceeding 0.7 or falling below 0.3 were selected, as such extreme levels are more likely indicative of lost expression or abundant expression. This approach aids in the discovery of potential biomarkers with functional relevance to disease phenotypes (16). Lastly, prior research has established that methylation patterns of multiple DMPs on promoters constitute key mechanisms influencing gene expression (15). Through these stringent criteria, we identified SLFN12 as a potential epigenetic biomarker gene for HT. The SLFN family, discovered by Schwarz in 1988, encompasses genes that induce T cell cycle arrest, thus deactivating T cells (32). Fortunately, numerous studies have explored the relationship between SLFN12 methylation and autoimmune diseases and inflammatory responses. Evidence suggests elevated SLFN12 methylation levels in CD4+ and CD8+ T cells in the peripheral blood of patients with multiple sclerosis, a chronic inflammatory disease of the central nervous system. Such changes to SLFN12 methylation pattern are thought to affect the activation and inflammatory response of T cells (33). In another study, researchers observed an increase in SLFN12 gene methylation levels in the peripheral blood of patients with allergic rhinitis, with methylation levels correlating positively with inflammation severity, suggesting SLFN12’s association with immune response and inflammation (34). In this research, we identified five significant hypomethylated DMPs in the SLFN12 gene promoter region in HT patients. Although this diverges from the hypermethylation pattern demonstrated in SLFN12 in other inflammatory diseases, the shift to a hypomethylation pattern in HT patients substantially affected the expression of the SLFN12 gene, showing a significant negative correlation. Additionally, a biomarker panel encompassing all discovered DMPs (especially cg11346248) on genes involving SLFN12 exhibited outstanding diagnostic potential, accurately differentiating between HT patients and healthy controls. Since thyroid function is a crucial HT disease indicator, this article also seeks to establish the precise correlation between SLFN12 methylation levels and thyroid function. Relevant research has suggested that genetic factors contribute to interindividual variations in FT4, FT3 and TSH levels (35). In 2021, a cohort study identified six CpG sites linked to elevated levels of FT3, and two CpG sites associated with elevated levels of TSH (36). Our findings indicated that methylation levels of all DMPs in the SLFN12 gene were negatively correlated with TSH. Notably, one CpG site (cg24470734) had the largest |Δβ| among all CpG sites, and its methylation levels were positively correlated with FT4. Therefore, this finding implies that cg24470734 could have a significant impact on the link between thyroid function indicators. Consequently, we surmised that SLFN12 gene methylation levels could affect thyroid function and predict changes in HT’s thyroid function. As a result, SLFN12 may have a vital role in HT pathogenesis and its methylation pattern might serve as a potential epigenetic marker for HT. Nevertheless, these conclusions necessitate further validation through expansion of the sample size and more rigorous experimental evidence.
The strengths of our investigation are multifaceted. Ours is the pioneering study to undertake a genome-wide methylation analysis in HT, thus laying a robust theoretical groundwork for the involvement of epigenetics in HT. Additionally, the novelty of our research lies in the fact that for the first time, biological processes and pathways related to HT were evaluated from a methylation standpoint using bioinformatic analysis. This enabled the identification of potential methylation biomarkers for HT, and their transcriptional levels were subsequently verified. Nonetheless, it is important to recognize the limitations in our current study. A primary shortcoming of our study is that we only performed DNA methylation profiling on leukocytes present in the whole blood. To gain a more detailed understanding of the epigenetic changes associated with HT, conducting cell type-specific analysis after whole blood samples sequencing would be beneficial. A secondary limitation of our study pertains to the relatively small sample size. In future studies, a broader participant representation is required to corroborate our findings. Additionally, we will further investigate the evolution of SLFN12 gene methylation levels in HT patients over time and disease progression, as well as explore the relationship between SLFN12 gene methylation levels and HT treatment.
5 Conclusion
In essence, this investigation lays the foundation for an initial genome-wide DNAm profile for patients with HT, facilitating a clear distinction between HT subjects and healthy individuals from an epigenetic standpoint. Modifications in the methylation landscape of key genes, orchestrating processes such as cell adhesion molecules, calcium signaling pathways, and Th17 cell differentiation, may hold implications for HT pathogenesis. The gene SLFN12 emerges as a potential epigenetic marker for HT and could bear substantial relevance to the disease’s progression. The comprehensive data yielded from this investigation offer significant insights, thereby enhancing our comprehension of aberrant DNAm contributing to HT’s pathogenesis.
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: BioProject accession number: PRJNA1005791.
Ethics statement
The studies involving humans were approved by Harbin Medical University Ethics Review Committee (No. hrbmuecdc20200320). The studies were conducted in accordance with the local legislation and institutional requirements. The participants provided their written informed consent to participate in this study.
Author contributions
ZZ: Conceptualization, Data curation, Investigation, Methodology, Software, Writing – original draft. JL: Writing – review & editing, Conceptualization, Data curation, Investigation, Software. YC: Writing – review & editing, Data curation, Investigation, Software, Supervision, Validation. BR: Data curation, Investigation, Methodology, Software, Writing – review & editing. SW: Conceptualization, Investigation, Writing – review & editing. YC: Writing – review & editing, Methodology, Supervision. YH: Validation, Writing – review & editing. QW: Software, Validation, Writing – review & editing. HG: Validation, Writing – review & editing. LL: Supervision, Validation, Visualization, Writing – review & editing. HS: Funding acquisition, Project administration, Resources, Supervision, Validation, Writing – review & editing.
Funding
The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This study was supported by grants from the National Natural Science Foundation of China (82073490).
Acknowledgments
We would like to thank the Institute of Endemic Disease Control of Shandong for their assistance in collecting samples and data. Thanks to all the volunteers who participated and donated samples to this study.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fendo.2023.1259903/full#supplementary-material
References
1. Caturegli P, De Remigis A, Rose NR. Hashimoto thyroiditis: clinical and diagnostic criteria. Autoimmun Rev (2014) 13(4-5):391–7. doi: 10.1016/j.autrev.2014.01.007
2. Skov J, Calissendorff J, Eriksson D, Magnusson P, Kämp O, Bensing S, et al. Limited genetic overlap between overt Hashimoto's thyroiditis and graves' Disease in twins: A population-based study. J Clin Endocrinol Metab (2021) 106(4):1101–10. doi: 10.1210/clinem/dgaa956
3. Hansen PS, Brix TH, Iachine I, Kyvik KO, Hegedüs L. The relative importance of genetic and environmental effects for the early stages of thyroid autoimmunity: a study of healthy Danish twins. Eur J Endocrinol (2006) 154(1):29–38. doi: 10.1530/eje.1.02060
4. Simmonds MJ. GWAS in autoimmune thyroid disease: redefining our understanding of pathogenesis. Nat Rev Endocrinol (2013) 9(5):277–87. doi: 10.1038/nrendo.2013.56
5. Moore LD, Le T, Fan G. DNA methylation and its basic function. Neuropsychopharmacology (2013) 38(1):23–38. doi: 10.1038/npp.2012.112
6. Shalaby SM, Mackawy AMH, Atef DM, Atef RM, Saeed J. Promoter methylation and expression of intercellular adhesion molecule 1 gene in blood of autoimmune thyroiditis patients. Mol Biol Rep (2019) 46(5):5345–53. doi: 10.1007/s11033-019-04990-6
7. Kyrgios I, Giza S, Fragou A, Tzimagiorgis G, Galli-Tsinopoulou A. DNA hypermethylation of PTPN22 gene promoter in children and adolescents with Hashimoto thyroiditis. J Endocrinol Invest (2021) 44(10):2131–8. doi: 10.1007/s40618-020-01463-7
8. Morris TJ, Butcher LM, Feber A, Teschendorff AE, Chakravarthy AR, Wojdacz TK, et al. ChAMP: 450k chip analysis methylation pipeline. Bioinformatics (2014) 30(3):428–30. doi: 10.1093/bioinformatics/btt684
9. Maksimovic J, Gordon L, Oshlack A. SWAN: Subset-quantile within array normalization for illumina infinium HumanMethylation450 BeadChips. Genome Biol (2012) 13(6):R44. doi: 10.1186/gb-2012-13-6-r44
10. Hochberg Y, Benjamini Y. More powerful procedures for multiple significance testing. Stat Med (1990) 9(7):811–8. doi: 10.1002/sim.4780090710
11. Wei S, Tao J, Xu J, Chen X, Wang Z, Zhang N, et al. Ten years of EWAS. Adv Sci (2021) 8(20):e2100727. doi: 10.1002/advs.202100727
12. Lafontaine N, Wilson SG, Walsh JP. DNA methylation in autoimmune thyroid disease. J Clin Endocrinol Metab (2023) 108(3):604–13. doi: 10.1210/clinem/dgac664
13. Zhang R, Chang C, Jin Y, Xu L, Jiang P, Wei K, et al. Identification of DNA methylation-regulated differentially expressed genes in RA by integrated analysis of DNA methylation and RNA-Seq data. J Transl Med (2022) 20(1):481. doi: 10.1186/s12967-022-03664-5
14. Wilson AS, Power BE, Molloy PL. DNA hypomethylation and human diseases. Biochim Biophys Acta (2007) 1775(1):138–62. doi: 10.1016/j.bbcan.2006.08.007
15. Lin R, Li C, Liu Z, Wu R, Lu J. Genome-wide DNA methylation profiling identifies epigenetic signatures of gastric cardiac intestinal metaplasia. J Transl Med (2020) 18(1):292. doi: 10.1186/s12967-020-02453-2
16. Qi Y, Zhang L, Yang X, Tang B, Xiao T. Genome-wide DNA methylation profile in whole blood of patients with chronic spontaneous urticaria. Front Immunol (2021) 12:681714. doi: 10.3389/fimmu.2021.681714
17. Joseph S, George NI, Green-Knox B, Treadwell EL, Word B, Yim S, et al. Epigenome-wide association study of peripheral blood mononuclear cells in systemic lupus erythematosus: Identifying DNA methylation signatures associated with interferon-related genes based on ethnicity and SLEDAI. J Autoimmun (2019) 96:147–57. doi: 10.1016/j.jaut.2018.09.007
18. Vymetalkova V, Vodicka P, Pardini B, Rosa F, Levy M, Schneiderova M, et al. Epigenome-wide analysis of DNA methylation reveals a rectal cancer-specific epigenomic signature. Epigenomics (2016) 8(9):1193–207. doi: 10.2217/epi-2016-0044
19. Li J, Liang Y, Fan J, Xu C, Guan B, Zhang J, et al. DNA methylation subtypes guiding prognostic assessment and linking to responses the DNA methyltransferase inhibitor SGI-110 in urothelial carcinoma. BMC Med (2022) 20(1):222. doi: 10.1186/s12916-022-02426-w
20. Lee HJ, Li CW, Hammerstad SS, Stefan M, Tomer Y. Immunogenetics of autoimmune thyroid diseases: A comprehensive review. J Autoimmun (2015) 64:82–90. doi: 10.1016/j.jaut.2015.07.009
21. Demirkaya E, Zhou Q, Smith CK, Ombrello MJ, Deuitch N, Tsai WL, et al. Brief report: deficiency of complement 1r subcomponent in early-onset systemic lupus erythematosus: the role of disease-modifying alleles in a monogenic disease. Arthritis Rheumatol (2017) 69(9):1832–9. doi: 10.1002/art.40158
22. Macedo AC, Isaac L. Systemic lupus erythematosus and deficiencies of early components of the complement classical pathway. Front Immunol (2016) 7:55. doi: 10.3389/fimmu.2016.00055
23. Velázquez FE, Anastasiou M, Carrillo-Salinas FJ, Ngwenyama N, Salvador AM, Nevers T, et al. Sialomucin CD43 regulates T helper type 17 cell intercellular adhesion molecule 1 dependent adhesion, apical migration and transendothelial migration. Immunology (2019) 157(1):52–69. doi: 10.1111/imm.13047
24. Beumer W, Effraimidis G, Drexhage RC, Wiersinga WM, Drexhage HA. Changes in serum adhesion molecules, chemokines, cytokines, and tissue remodeling factors in euthyroid women without thyroid antibodies who are at risk for autoimmune thyroid disease: a hypothesis on the early phases of the endocrine autoimmune reaction. J Clin Endocrinol Metab (2013) 98(6):2460–8. doi: 10.1210/jc.2012-4122
25. Nanba T, Watanabe M, Inoue N, Iwatani Y. Increases of the Th1/Th2 cell ratio in severe Hashimoto's disease and in the proportion of Th17 cells in intractable Graves' disease. Thyroid (2009) 19(5):495–501. doi: 10.1089/thy.2008.0423
26. Li D, Cai W, Gu R, Zhang Y, Zhang H, Tang K, et al. Th17 cell plays a role in the pathogenesis of Hashimoto's thyroiditis in patients. Clin Immunol (2013) 149(3):411–20. doi: 10.1016/j.clim.2013.10.001
27. Ichiyama K, Sekiya T, Inoue N, Tamiya T, Kashiwagi I, Kimura A, et al. Transcription factor Smad-independent T helper 17 cell induction by transforming-growth factor-β is mediated by suppression of eomesodermin. Immunity (2011) 34(5):741–54. doi: 10.1016/j.immuni.2011.02.021
28. Zhang Q, Cui F, Fang L, Hong J, Zheng B, Zhang JZ. TNF-α impairs differentiation and function of TGF-β-induced Treg cells in autoimmune diseases through Akt and Smad3 signaling pathway. J Mol Cell Biol (2013) 5(2):85–98. doi: 10.1093/jmcb/mjs063
29. Kaufmann U, Kahlfuss S, Yang J, Ivanova E, Koralov SB, Feske S. Calcium signaling controls pathogenic Th17 cell-mediated inflammation by regulating mitochondrial function. Cell Metab (2019) 29(5):1104–18.e1106. doi: 10.1016/j.cmet.2019.01.019
30. Hemon P, Renaudineau Y, Debant M, Le Goux N, Mukherjee S, Brooks W, et al. Calcium signaling: from normal B cell development to tolerance breakdown and autoimmunity. Clin Rev Allergy Immunol (2017) 53(2):141–65. doi: 10.1007/s12016-017-8607-6
31. Guo H, Zhu P, Yan L, Li R, Hu B, Lian Y, et al. The DNA methylation landscape of human early embryos. Nature (2014) 511(7511):606–10. doi: 10.1038/nature13544
32. Xu J, Chen S, Liang J, Hao T, Wang H, Liu G, et al. Schlafen family is a prognostic biomarker and corresponds with immune infiltration in gastric cancer. Front Immunol (2022) 13:922138. doi: 10.3389/fimmu.2022.922138
33. Rhead B, Brorson IS, Berge T, Adams C, Quach H, Moen SM, et al. Increased DNA methylation of SLFN12 in CD4+ and CD8+ T cells from multiple sclerosis patients. PLoS One (2018) 13(10):e0206511. doi: 10.1371/journal.pone.0206511
34. North ML, Jones MJ, MacIsaac JL, Morin AM, Steacy LM, Gregor A, et al. Blood and nasal epigenetics correlate with allergic rhinitis symptom development in the environmental exposure unit. Allergy (2018) 73(1):196–205. doi: 10.1111/all.132
35. Samollow PB, Perez G, Kammerer CM, David F, Patrick Z, Lorena MH, et al. Genetic and environmental influences on thyroid hormone variation in Mexican Americans. J Clin Endocrinol Metab (2004) 89(7):3276–84. doi: 10.1210/jc.2003-031706
Keywords: Hashimoto thyroiditis, DNA methylation, epigenetics, autoimmune endocrinopathies, SLFN12
Citation: Zhou Z, Liu J, Chen Y, Ren B, Wan S, Chen Y, He Y, Wei Q, Gao H, Liu L and Shen H (2023) Genome-wide DNA methylation pattern in whole blood of patients with Hashimoto thyroiditis. Front. Endocrinol. 14:1259903. doi: 10.3389/fendo.2023.1259903
Received: 17 July 2023; Accepted: 16 November 2023;
Published: 24 November 2023.
Edited by:
Zhiguo Xie, Central South University, ChinaReviewed by:
Xiaofeng Wu, Pfizer (United States), United StatesNelson George, Eunice Kennedy Shriver National Institute of Child Health and Human Development (NIH), United States
Copyright © 2023 Zhou, Liu, Chen, Ren, Wan, Chen, He, Wei, Gao, Liu and Shen. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Hongmei Shen, c2hlbmhtMTE5QGhyYm11LmVkdS5jbg==; Lixiang Liu, bGl1bGl4aWFuZ0BocmJtdS5lZHUuY24=
†These authors share first authorship