- Department of Neurology, The Second Affiliated Hospital of Harbin Medical University, Harbin, China
Myasthenia gravis (MG) is an autoimmune disease associated with autoantibody production that leads to skeletal muscle weakness. The molecular mechanisms underlying MG are not fully understood. We analyzed the gene expression profile (GSE85452) and methylation profile (GSE85647) of MG samples from the GEO database to identify aberrantly methylated-differentially expressed genes. By integrating the datasets, we identified 143 hypermethylation-low expression genes and 91 hypomethylation-high expression genes. Then we constructed PPI network and ceRNA networks by these genes. Phosphatase and tensin homolog (PTEN) and Abelson tyrosine-protein kinase (ABL)1 were critical genes in both PPI networks and ceRNA networks. And potential MG associated lncRNAs were selected by comprehensive analysis of the critical genes and ceRNA networks. In the hypermethylation-low expression genes associated ceRNA network, sirtuin (SIRT)1 was the most important gene and the lncRNA HLA complex (HC) P5 had the highest connection degree. Meanwhile, PTEN was the most important gene and the lncRNA LINC00173 had the highest connection degree in the hypomethylation-high expression genes associated ceRNA network. LINC00173 was validated to be upregulated in MG patients by qRT-PCR (P = 0.005), which indicated LINC00173 might be a potential biomarker for MG. These results provide a basis for future studies on the molecular pathogenesis of MG.
Introduction
Myasthenia gravis (MG) is an autoimmune disease that involves cellular immunity and complement activation and is characterized by impaired transmission at the neuromuscular junction, leading to skeletal muscle weakness and autoantibody production. The disease is mainly caused by the loss of immune tolerance, which induces the production of antibodies mainly against acetylcholine receptor (AChR), triggering an abnormal immune response (Gilhus and Verschuuren, 2015). The underlying cause of MG is unclear, but genetic factors may determine an individual’s disease risk. For example, the association between human leukocyte antigen (HLA) class I and II genes and MG is well established (Vandiedonck et al., 2005). MG-associated susceptibility genes such as those encoding cytotoxic T lymphocyte-associated protein (CTLA)-4, tumor necrosis factor (TNF)-α, protein tyrosine phosphatase non-receptor type (PTPN)22, and interleukin (IL)-10 have been implicated in other autoimmune diseases (Avidan et al., 2014).
In addition to the abovementioned genes, environmental factors such as pollutants, drugs, and pathogens are thought to increase the risk of MG (Cavalcante et al., 2013). Epigenetic modifications play an important role in gene regulation; they include DNA methylation, chromatin remodeling, non-coding RNA, and histone modifications. DNA methylation involves the addition of a methyl group to the 5th carbon position of cytosine at CpG dinucleotides, which are typically found in gene promoters; this prevents transcription factors from binding, thereby inhibiting gene expression (Schultz et al., 2015). CTLA4 methylation may promote the development of MG by increasing cytokine expression through upregulation of autoantibodies against AChR and erythrocyte acetylcholinesterase (Fang et al., 2018). Growth hormone secretagogue receptor (GHSR) hypermethylation was observed in thymoma samples of MG patients, especially in late stages of disease (Coppede et al., 2020). It was also reported that G0/G1 switch (G0S)2 was upregulated in B cells and CD8+ T cells of AChR-positive MG patients, but this was significantly attenuated by G0S2 methylation (Xu et al., 2020). These findings suggest that gene methylation contributes to the pathogenesis of MG.
Non-coding (nc) RNAs such as micro (mi) RNAs and long (l) ncRNAs modulate the expression of protein-coding genes in a variety of biological processes and play a key role in the occurrence of MG (Jonas and Izaurralde, 2015; Chen et al., 2017). For example, members of the lethal (let)-7 lncRNA family were downregulated in peripheral blood mononuclear cells (PBMCs) of MG patients, and the level of let-7c was negatively correlated with that of IL-10 (Jiang et al., 2012). Moreover, the lncRNA interferon gamma (IFN-γ) antisense RNA (IFNG-AS)1 was found to be differentially expressed in MG patients and negatively regulated the expression of HLA-DRB and HLA-DOB (Luo et al., 2017). The competing endogenous (ce) RNA theory states that transcripts sharing common miRNA binding sites compete to bind identical miRNAs for mutual regulation of expression levels (Salmena et al., 2011). Moreover, mRNAs, lncRNAs, and other types of RNA can adsorb miRNAs through miRNA binding elements and thereby affect their function. LncRNAs are endogenous sponges that regulate miRNA activity (Wang et al., 2015). The lncRNA small nucleolar RNA host gene (SNHG)16 was shown to regulate the expression of IL-10 by sponging let-7c-5p as a ceRNA in MG (Wang J. et al., 2020); and metastasis-associated lung adenocarcinoma transcript (MALAT)1, acting as a ceRNA for the miRNA miR-338-3p, directly induced male-specific lethal complex subunit (MSL)2 expression in MG (Kong et al., 2019). However, it remains unclear whether aberrantly methylated genes are related to ceRNAs in MG.
To answer this question, in this study, we analyzed the gene expression and methylation profiles of MG samples from the Gene Expression Omnibus (GEO) database. Hypermethylated genes with low expression (hypermethylation-low expression genes) and hypomethylated genes with high expression (hypomethylation-high expression genes) were identified by integrating differentially expressed genes (DEGs) and differentially methylated genes (DMGs). Functional enrichment analysis and protein–protein interactions (PPI) were investigated by bioinformatic methods. Based on the results, we constructed an aberrantly methylated-differentially expressed genes associated ceRNA network (AMCEN) for MG and selected potential biomarkers by comprehensive analysis of critical genes and ceRNA networks. This study may provide molecular-level insight into the pathogenesis of MG.
Materials and Methods
Data Sources
Gene expression (GSE85452) and gene methylation (GSE85647) datasets for MG were obtained from the GEO.1 GSE85452 (platform: GPL10558 HumanHT-12 v4.0 Expression BeadChip; Illumina, San Diego, CA, United States) was based on monocytes from 10 MG patients and 9 controls; and GSE85647 (platform: GPL13534 HumanMethylation450 BeadChip; Illumina) was based on monocytes from 9 MG patients and 7 controls. In order to avoid the influence of genomic similarity of monozygotic twin in the datasets on the experimental results, if monozygotic twins in the datasets were not all MG patients, they were excluded from the analysis. We manually searched PubMed2 for MG risk miRNAs in articles published before December 1st, 2020 using the keywords “microRNA” or “miRNA” or “miR” and “MG,” with the species limited to “Homo sapiens.” We then manually selected miRNAs that met the following criteria: (i) differentially expressed with statistical significance between MG patients and healthy controls; and (ii) detected by reliable experimental methods (microarray, real-time PCR, etc.). Experimentally verified mRNA–miRNA interactions were downloaded from miRTarBase (release 7.0) (Chou et al., 2018); the data were supported by luciferase reporter assay or western blotting. miRNA–lncRNA interactions were obtained from starBase v3.0 (Li et al., 2014), DIANA-LncBase v2.0 (Paraskevopoulou et al., 2013), and LncACTdb v2.0 (Wang et al., 2019) which contain interactions that have been high-throughput experimentally validated (e.g., HITS-CLIP, PAR-CLIP, iCLIP, CHIP, and CLASH).
Data Processing
We used GEO2R online software3 to analyze expression data for MG patients and healthy controls. For the GES85452 and GSE85647 datasets, P < 0.05 and | t| > 2 was used as the cutoff criteria to identify DEGs and DMGs. Hypermethylation-low expression genes were obtained by overlapping hypermethylated and downregulated genes, and hypomethylation-high expression genes were obtained by overlapping hypomethylated and upregulated genes. DMG–miRNA interactions were determined by matching hypermethylation-low expression genes or hypomethylation-high expression genes with mRNA–miRNA interaction data from miRTarBase (release 7.0). We then identified mRNA–miRNA–lncRNA ceRNA triples by overlapping miRNAs from DMG–miRNA and miRNA–lncRNA interaction datasets. To increase the reliability of the data, the manually selected MG-associated miRNAs were also matched with the mRNA–miRNA–lncRNA ceRNA triples to obtain MG-associated aberrantly methylated-differentially expressed involved ceRNA interactions.
Functional Enrichment Analysis
We performed functional enrichment analysis of hypermethylation-low expression and hypomethylation-high expression genes with the Database for Annotation, Visualization and Integrated Discovery (DAVID) platform4 (Dennis et al., 2003). Gene Ontology (GO) analysis included biological process, cellular component and molecular function terms (Gene Ontology, 2006); and we also conducted a Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis (Kanehisa and Goto, 2000). P < 0.05 was set as the cutoff for statistical significance.
Analysis of GSEA Enrichment
We next performed GSEA analysis to avoid missing the genes which actually play an important part during the process of identifying differential expressed genes among the data. Gene set permutations were performed 1,000 times for each analysis. The normalized enrichment score (NES) > 1 and nominal P-value < 0.05 were regarded as significant.
PPI Network Construction and Module Analysis
The construction of PPI networks is important for elucidating the molecular mechanisms underlying cellular activities. We used the Search Tool for the Retrieval of Interacting Genes (STRING)5 database and Cytoscape v3.8.16 to construct PPI networks of hypermethylation-low expression and hypomethylation-high expression genes, respectively. An interaction score of 0.4 was taken as the cutoff. The Molecular Complex Detection (MCODE) plugin of Cytoscape software was used to define modules in the PPI network with degree cutoff = 2, node score cutoff = 0.2, k-core = 2, and maximum depth = 100. Hub genes were selected with cytoHubba based on a connection degree > 10.
Cumulative Hypergeometric Test
Aberrantly methylated mRNA–lncRNA ceRNA interaction pairs sharing the same miRNA were identified by cumulative hypergeometric testing (Zhang G. et al., 2018). The P-value for each interaction pair was computed using the following formula:
where N is the total number of miRNAs from DMG–miRNA interactions; n and m are the number of miRNAs associated with 1 mRNA and 1 lncRNA, respectively, and x is the number of miRNAs shared by the mRNA and lncRNA. Interactions of aberrantly methylated mRNA–lncRNA ceRNA pairs with a P-value < 0.05 were considered significant.
LncRNA–mRNA Co-expression Analysis
We performed a co-expression correlation analysis of lncRNA–mRNA interaction pairs by calculating the Pearson correlation coefficient (PCC) of lncRNA and mRNA expression data downloaded from the dbGaP database (Tryka et al., 2014). The P-values of the co-expression analysis were subjected to false discovery rate (FDR) adjustment. The cutoff values were PCC > 0.5 and FDR < 0.05.
Construction of the AMCEN
We constructed the aberrantly methylated-differentially expressed genes associated ceRNA network (AMCEN) based on aberrantly methylated-differentially expressed mRNA–miRNA–lncRNA ceRNA triples using Cytoscape v3.8.1. The AMCEN for MG was constructed based on the theory that lncRNAs share common miRNA binding sites with mRNAs and function as ceRNAs to regulate mRNAs. For a given lncRNA–miRNA–mRNA interaction, both mRNA and lncRNA shared common miRNAs and were co-expressed for merging into a competing triplet. After identifying and assembling all lncRNA–miRNA–mRNA competing triplets by using cumulative hypergeometric test and PCC analysis, we constructed the AMCEN. The network was visualized using Cytoscape software, in which nodes represent miRNAs, genes and lncRNAs, and edges represent their interactions. Topologic features including degree and betweenness distribution were analyzed for all nodes in the AMCEN.
Clinical Samples
A total of 27 patients with myasthenia gravis (MG group) admitted to the Second Affiliated Hospital of Harbin Medical University were selected as the experimental subjects. A total of 27 sex- and age- matched healthy subjects without autoimmune diseases and without any drug therapies were selected as the control group (Control group). According to the Chinese Guidelines for the Diagnosis and Treatment of Myasthenia Gravis (2015 Edition), MG patients should meet the following diagnostic criteria: first, patients should have the typical clinical characteristics of fluctuating fatigue of muscles, and in addition, it should meet any of the following three criteria: (i) pharmacological examination: positive for neostigmine test; (ii) Electrophysiological characteristics: repeated electrical nerve stimulation recorded complex muscle action and single fiber EMG recorded the neuromuscular junction for longer transmission; (iii) Serum antibody was positive. All the selected MG patients met the above diagnostic criteria and did not receive immunosuppressant treatment, and obtained informed consent of the selected patients. Peripheral blood samples were collected from each subject in tubes containing ethylenediaminetetraacetic acid, and PBMCs were isolated using lymphocyte separation medium. This study was approved by the Ethics Committee of The Second Affiliated Hospital of Harbin Medical University.
Real-Time PCR Analysis
Total RNA was extracted from PBMCs using Trizol reagent (Sigma Life Science, Darmstadt, Germany) following the manufacturer’s instructions. Reverse transcription of total RNA into complementary DNA (cDNA) was performed using a Transcriptor First Strand cDNA Synthesis Kit (Roche, Switzerland) according to the manufacturer’s instructions. LINC00173 expression level was detected by quantitative real-time polymerase chain reaction (qRT-PCR) using FastStart Universal SYBR Green Master (Roche, Switzerland). The relative expression level was calculated using the 2–ΔΔCt method using glyceraldehyde 3-phosphate dehydrogenase for normalization. The primer used for qRT-PCR was illustrated in Supplementary Table 1.
Results
Identification of DEGs and DMGs in MG
The GSE85452 and GSE85647 datasets were separately screened with GEO2R online software for DEGs or DMGs. In GSE85452, 987 DEGs were upregulated and 771 were downregulated and in GSE85647, 3336 genes were hypermethylated and 1810 were hypomethylated. A total of 143 hypermethylation-low expression genes were obtained by overlapping 3336 hypermethylated and 771 downregulated genes; and 91 hypomethylation-high expression genes were obtained by overlapping 1,810 hypomethylated and 987 upregulated genes (Figure 1).
Figure 1. Analysis of differentially expressed genes and differentially methylated genes. (A) Representative heatmap of the top 40 differentially expressed genes in dataset GSE85452 (20 up-regulated genes and 20 down-regulated genes). (B) Representative volcanic map of differentially expressed genes distribution in GES85452. (C) Identification of aberrantly methylated-differentially expressed genes in gene expression dataset (GSE85452) and gene methylation dataset (GSE85647).
GO and KEGG Pathway Analyses of Aberrantly Methylated DEGs
To determine the functions of the identified hypermethylation-low expression and hypomethylation-high expression genes, GO and KEGG pathway analyses were carried out with DAVID. The top five significant GO terms and KEGG pathways are shown in Figure 2 and Table 1. Hypermethylation-low expression and hypomethylation-high expression genes were enriched in the biological processes of cellular response to DNA damage stimulus; covalent chromatin modification; antigen processing and presentation of exogenous peptide antigen via major histocompatibility complex (MHC) class I, TAP-dependent; IFN-γ-mediated signaling pathway; and type I IFN signaling pathway. For cellular component, the genes showed enrichment in nucleoplasm, membrane, cytosol, nucleus, and cytoplasm; and for molecular function, the genes were mostly enriched in protein binding. These results suggest that the identified genes play an important role in the pathogenesis of MG. Additionally, the KEGG pathway enrichment analysis showed that the genes were significantly enriched in viral myocarditis and antigen processing and presentation, indicating that the pathogenesis of MG is closely related to autoimmunity.
Figure 2. GO functional enrichment and KEGG pathway analysis of hypermethylation-low expression genes and hypomethylation-high expression genes. (A) Bar graph of BP, CC, MF and KEGG. (B) Distribution of hypermethylation-low expression genes and hypomethylation-high expression genes in each GO term of BP, CC, MF. Red represents hypermethylation-low expression genes; Blue represents hypomethylation-high expression genes. (C) KEGG pathway analysis of hypermethylation-low expression genes and hypomethylation-high expression genes. Different colors represent different pathways. (D) Genes distribution in different KEGG pathways.
Table 1. GO and KEGG analysis of hypermethylation-low expression genes and hypomethylation-high expression genes.
GSEA Enrichment in GSE85452
GSEA enrichment analysis showed that lipid droplet organization, chaperone binding, actin filaments, regulation of cellular amino acid metabolic process, negative regulation of BMP signaling pathway, endopeptidase complex, viral release from host cell, nucleobase containing compound kinase activity, protein demethylase activity had the largest NES, and the details were reported in Supplementary Figure 1 and Supplementary Table 2.
Construction of a PPI Network, Module Analysis, and Hub Gene Selection
To investigate the interactions of hypermethylation-low expression or hypomethylation-high expression genes, we used the STRING database along with MCODE and cytoHubba to construct a PPI network and define the modules and hub genes (Figures 3A,B). The top 5 hub genes for the hypermethylation-low expression gene PPI network were Von Hippel–Lindau tumor suppressor (VHL), zinc finger and BTB domain-containing protein (ZBTB)16, WD repeat and SOCS box-containing protein (WSB)1, tripartite motif containing (TRIM)4, and ring finger protein (RNF)144B; and the top 5 hub genes in the hypomethylation-high expression gene PPI network were phosphatase and tensin homolog (PTEN), ATP synthase F1 subunit delta (ATP5D), Abelson tyrosine-protein kinase (ABL)1, poly(rC)-binding protein (PCBP)2, and human leukocyte antigen (HLA) class II histocompatibility antigen, DP alpha 1 chain (HLA-DPA1) (Figures 3C,D).
Figure 3. PPI networks and modules of hypermethylation-low expression genes and hypomethylation-high expression genes. (A) PPI network of hypermethylation-low expression genes. The bigger and darker the dot is, the higher degree the dot has. (B) Top two modules of hypermethylation-low expression genes PPI network. (C) PPI network of hypomethylation-high expression genes. The bigger and darker the dot is, the higher degree the dot has. (D) The top module of hypomethylation-high expression genes PPI network.
AMCEN and Topologic Analysis
We constructed aberrantly methylated gene-associated ceRNA networks by matching hypermethylation-low expression or hypomethylation-high expression genes with mRNA–miRNA and miRNA–lncRNA interactions and performing cumulative hypergeometric testing and lncRNA–mRNA co-expression analysis to screen triples for AMCEN construction (Figures 4, 5). We calculated the betweenness of the nodes in the AMCEN, with higher values reflecting a more important role in maintaining tight connectivity in the network. There were 744 mRNA–miRNA–lncRNA ceRNA triples of hypermethylation-low expression genes with 22 mRNAs, 46 miRNAs, and 151 lncRNAs (Figure 4); and 197 mRNA–miRNA–lncRNA ceRNA triples of hypomethylation-high expression genes with 9 mRNAs, 36 miRNAs, and 81 lncRNAs (Figure 5). In the hypermethylation-low expression gene-associated ceRNA network, sirtuin (SIRT)1 was the most important gene and the lncRNA HCP5 had the highest degree. In the hypomethylation-high expression gene-associated ceRNA network, PTEN was the most important gene and the lncRNA LINC00173 had the highest degree.
Figure 4. Construction of AMCEN by hypermethylation-low expression genes and topological property. (A) The ceRNA network. Purple circles represent mRNAs, red triangles represent miRNAs, green diamonds represent lncRNAs and lines represent their regulatory interactions. (B) The node degree distribution of the global ceRNA network. (C) The node betweenness distribution of the network. (D) Degree distribution of hypermethylation-low expression genes. (E) Degree distribution of lncRNAs.
Figure 5. Construction of AMCEN by hypomethylation-high expression genes and topological property. (A) The ceRNA network. Purple circles represent mRNAs, red triangles represent miRNAs, green diamonds represent lncRNAs and lines represent their regulatory interactions. (B) The node degree distribution of the ceRNA network. (C) The node betweenness distribution of the network. (D) Degree distribution of lncRNAs.
miRNAs are a type of ceRNA that play an important role in ceRNA networks. We manually selected 192 MG risk miRNAs to increase the reliability of our analysis and overlapped these with aberrantly methylated ceRNA interactions. This yielded a hypermethylation-low expression gene-associated ceRNA network comprising 333 mRNA–miRNA–lncRNA ceRNA triples including 14 genes, 17 miRNAs, and 122 lncRNAs (Figure 6A); and a hypomethylation-high expression gene-associated ceRNA network comprising 97 mRNA–miRNA–lncRNA ceRNA triples including 6 genes, 14 miRNAs, and 49 lncRNAs (Figure 6B). These aberrantly methylated DEGs and ceRNA networks are presumed to be involved in the pathogenesis of MG.
Figure 6. AMCENs after matching miRNAs. (A) Hypermethylation-low expression genes mediated ceRNA network. Purple circles represent mRNAs, red triangles represent miRNAs, green diamonds represent lncRNAs and lines represent their regulatory interactions. (B) Hypomethylation-high expression genes mediated ceRNA network. Purple circles represent mRNAs, red triangles represent miRNAs, green diamonds represent lncRNAs and lines represent their regulatory interactions.
Identification of Biomarkers Based on Aberrantly Methylated-Differentially Expressed Genes and ceRNA Networks
By matching genes in ceRNA networks and hub genes of the MG PPI network, we identified miRNAs that interact with PTEN including hsa-miR-106b, hsa-miR-10b, hsa-miR-181a, hsa-miR-181b, hsa-miR-181c, hsa-miR-193a, hsa-miR-221, hsa-miR-23a, hsa-miR-25, hsa-miR-26a, hsa-miR-425, and hsa-miR-93. The lncRNAs co-expressed with PTEN in the ceRNA networks were LINC00173, family with sequence similarity 13, antisense RNA A1 (FAM13A-AS1), ankyrin repeat and SOCS box-containing antisense RNA (ASB16-AS)1, Opa-interacting protein 5 antisense RNA (OIP5-AS)1, transducer of ERBB2 antisense RNA (TOB1-AS)1, LINC01126, nucleolar protein interacting with the FHA domain of MKI67 antisense RNA (NIFK-AS)1, LINC01001, muscleblind-like splicing regulator antisense RNA (MBNL1-AS)1, nuclear transcription factor Y subunit gamma antisense RNA (NFYC-AS)1, taurine upregulated gene (TUG)1, LINC00593, DLG-associated protein 1 antisense RNA (DLGAP1-AS)2, HLA complex group (HCG)11, LINC00467, and WD repeat and FYVE domain-containing 3 antisense RNA (WDFY3-AS)2. The miRNA hsa-miR-30a combined with ABL1 in the ceRNA networks; the co-expressed lncRNAs were erythrocyte membrane protein band 4.1-like 4A antisense RNA (EPB41L4A-AS)1, small nucleolar RNA host gene (SNHG)1, sorting nexin 29 pseudogene (SNX29P)2, C terminal-binding protein antisense RNA (CTBP1-AS)2, plasmacytoma variant translocation (PVT)1, metastasis-associated lung adenocarcinoma transcript (MALAT)1, miR-22 host gene (MIR22HG), U-bx domain-containing 5 antisense RNA (UBOX5-AS)1, StAR-related lipid transfer domain-containing 7 antisense 1 (STARD7-AS1), LINC00115, keratin-associated protein 5-5 antisense RNA (KRTAP5-AS)1, tumor protein 73 antisense RNA (TP73-AS)1, LINC00469, LINC00926, and LINC00667. Based on the degree of lncRNAs in ceRNA networks, the lncRNAs LINC00173, FAM13A-AS1, and OIP5-AS1 were found to be associated with PTEN in ceRNA networks (Table 2). Thus, PTEN and ABL1 are critical genes in the pathogenesis of MG and the lncRNAs LINC00173, FAM13A-AS1, and OIP5-AS1 are potential disease biomarkers.
Validation of LINC00173 Expression
The expression level of LINC00173 was evaluated by qRT-PCR in PBMCs of MG patients and control subjects. LINC00173 was upregulated in MG patients compared with controls (P = 0.005), which indicated LINC00173 might be a potential biomarker for MG (Figure 7).
Figure 7. LINC00173 expression in MG. LINC00173 was up regulated in PBMCs from MG patients compared with healthy controls.
Discussion
MG is an autoantibody-mediated disease of the nervous system with an ill-defined etiology. Because of the complex pathogenesis and heterogeneity of MG, there is no single treatment that can be applied to all MG patients (Sanders et al., 2016). DNA methylation marks can serve as a biomarker for MG (Fang et al., 2018). The discovery of ceRNA networks has provided a basis for investigating the molecular mechanisms of MG, which could aid diagnosis, treatment, and prognostic assessment. In this study, we identified 143 hypermethylation-low expression and 91 hypomethylation-high expression genes by analyzing microarray-based gene expression and gene methylation data for MG. DEGs and DMGs reflected the genes and pathways modulated by aberrant methylation in MG. We then constructed an AMCEN to determine the main factors contributing to MG pathogenesis.
Hypermethylation-low expression and hypomethylation-high expression genes in MG were enriched in biological processes of cellular response to DNA damage stimulus; covalent chromatin modification; antigen processing and presentation of exogenous peptide antigen via MHC class I, TAP-dependent; IFN-γ–mediated signaling pathway; and type I IFN signaling pathway. As an autoimmune disease, MG is caused by both genetic and external factors (Meriggioli and Sanders, 2009). Genetics can influence disease susceptibility, which can be modulated by external factors such as infection or diet (Gilhus and Verschuuren, 2015) that lead to chromatin modification. Infections can induce the production of antibodies that act as autoantibodies against AChR; in fact, exacerbation of MG is often associated with infection (Hehir and Silvestri, 2018). Thus, aberrant antigen processing and presentation may contribute to the pathogenesis of MG. IFN-γ level was shown to be elevated in MG patients (Link et al., 1994). Viral infection and an increased level of type I IFN have been linked to thymic changes similar to those observed in MG patients with anti-AChR antibodies (Berrih-Aknin and Le Panse, 2014). Consistent with a viral etiology for MG and the role of immune responses, the results of the KEGG pathway analysis revealed significant enrichment of the viral myocarditis and antigen processing and presentation pathways.
The top 5 hub genes in the PPI network of hypermethylation-low expression genes were VHL, ZBTB16, WSB1, TRIM4, and RNF144b. VHL disease is caused by mutation of the VHL gene on chromosome 3p25-26 (Kaelin, 2007). A case was reported of an antibody-positive MG patient with thymoma concurrent with VHL (Sheth et al., 2005). Although this was a unique case, it suggests a potential genetic link between MG and VHL. ZBTB16, also known as promyelocytic leukemia zinc finger (PLZF), is a transcription factor that promotes the recruitment of effector T helper cells during the development of innate lymphocyte lineages and is also essential for the development of osteoblasts and spermatogonia (Mao et al., 2017). Another study found that DNA methylation in the promoter region of the ZBTB16 gene may regulate gene expression during osteoblast differentiation (Marofi et al., 2019). Thus, ZBTB16 may contribute to abnormal T cell development in MG, and ZBTB16 methylation is a potential mechanism of MG pathogenesis. WSB1, TRIM4, and RNF144b are associated with multiple cancers including neuroblastoma (Shichrur et al., 2014), chordoma (Wang C. B. et al., 2020) and hepatocellular carcinoma (Dong et al., 2020). Abnormal methylation of TRIM4 in immune pathways may be associated with neural tube defects (Zhang et al., 2019). RNF144B is thought to be involved in the development of chordoma through a ceRNA network (Wang C. B. et al., 2020). In the PPI network of hypomethylation-high expression genes, the top 5 hub genes were PTEN, ATP5D, ABL1, PCBP2, and HLA-DPA1. PTEN is a tumor suppressor and regulator of metabolism (Chen et al., 2018) that is expressed in some types of thymoma, although methylation of the gene has not been detected in thymoma samples (Masunaga et al., 2017). ATP5D was shown to be differentially expressed in amyotrophic lateral sclerosis patients and was linked to mitochondrial dysfunction in synaptic clefts (Engelen-Lee et al., 2017). ABL1 and PCBP2 are associated with immunity (Shi et al., 2014; Rafei et al., 2019), and HLA-DPA1 has been implicated in juvenile-onset MG (Feng et al., 2019). Elucidating the functions of these genes may provide new insights into MG pathogenesis.
PTEN and ABL1 were found to be components of hypermethylation-low expression and hypomethylation-high expression gene-associated ceRNA networks. PTEN was involved in 33 mRNA–miRNA–lncRNA triples. Previous study had found that LINC00173 promotes the proliferation and migration of vascular endothelial cells in lung squamous cell carcinoma by acting as ceRNA to bind with miR-511-5p and regulate VEGFA expression (Chen et al., 2020). LINC00173, which had the highest node degree in the hypomethylation-high expression gene-associated ceRNA network, has been investigated as part of the ceRNA networks of lung cancer (Chen et al., 2020), breast cancer (Fan et al., 2020), and glioma (Du et al., 2020), but its molecular mechanism in autoimmune diseases of the nervous system was still unclear. As a constituent of one of the 15 ABL1-associated ceRNA triples, the lncRNA MALAT1 was shown to be downregulated in MG patients and to act as a ceRNA for miR-338-3p (Kong et al., 2019). In the hypermethylation-low expression gene-associated ceRNA network, SIRT1 had the highest node degree. SIRT1 has been linked to immune-related diseases such as tumors and autoimmune diseases (Yu et al., 2018) and was decreased in the PBMCs of multiple sclerosis patients during relapses (Martin et al., 2015). The lncRNAs HCP5, OIP5-AS1, and LINC00894 were the most important lncRNAs in the ceRNA network. HCP5 is associated with several autoimmune diseases including psoriatic arthritis (Liu et al., 2008), systemic lupus erythematosus (Ciccacci et al., 2014), and Graves’ disease (Lane et al., 2020); and OIP5-AS1 was found to play a critical role in the ceRNA network of multiple sclerosis (Ding et al., 2021). LINC00894 may be related to drug resistance in breast cancer (Zhang X. et al., 2018). NR4A3, which is associated with MG (Bernasconi et al., 2003), was involved in 3 ceRNA triples; one of these includes zinc finger 674 antisense RNA (ZNF674-AS)1, which has been investigated in relation to cancers. Little is known about the functions of VAC14-AS1 and SNX29P2.
This study had several limitations. Firstly, lncRNAs and the ceRNA hypothesis have not been extensively investigated in MG, so most lncRNAs that have not been previously reported in this context and require validation. Secondly, the DEGs and aberrantly methylated genes were screened from an online database, and the relatively small number of samples reduced the reliability of the data. Finally, we constructed the AMCEN and identified several critical ceRNA triples, but the actual relationships between their components need to be verified experimentally.
Conclusion
In Conclusions, we identified 143 hypermethylation-low expression and 91 hypomethylation-high expression genes by overlapping DEGs and DMGs. These genes were mainly enriched in chromatin modification and immune process. Phosphatase and tensin homolog (PTEN) and Abelson tyrosine-protein kinase (ABL)1 were critical genes in both PPI networks and ceRNA networks. This revealed that PTEN and ABL1 were critical genes while the lncRNAs LINC00173, FAM13A-AS1, and OIP5-AS1 were potential biomarkers in MG. LINC00173 was validated to be upregulated in MG patients compared with controls. These findings provide insight into the molecular pathogenesis of MG as well as novel therapeutic targets for disease treatment.
Data Availability Statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.
Ethics Statement
The studies involving human participants were reviewed and approved by the Ethics Committee of The Second Affiliated Hospital of Harbin Medical University. The patients/participants provided their written informed consent to participate in this study.
Author Contributions
SX, JW, and LW conceived, designed the study, and revised the manuscript. SX, TW, XL, and HZ collected the data. LL, XK, and SL were performed the bioinformatics analysis. SX, XW, and HG analyzed and visualized the results. SX and LL conducted the experiments. SX performed the statistical analyses. SX and TW wrote the manuscript. All authors read, edited and approved the final manuscript.
Funding
This work was supported by the National Natural Science Foundation of China (Grant Numbers 81820108014, 82071407, 81771361, and 81901277), the National Key Research and Development Project (Grant Number 2018YFE0114400), and the Postdoctoral Foundation of Heilongjiang Province (Grant Number LBH-TZ1019).
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/fgene.2021.726751/full#supplementary-material
Supplementary Figure 1 | Gene Set Enrichment Analysis (GSEA) of mRNAs in GSE85452.
Footnotes
- ^ https://www.ncbi.nlm.nih.gov/geo/
- ^ http://www.ncbi.nlm.nih.gov/pubmed
- ^ https://www.ncbi.nlm.nih.gov/geo/geo2r/
- ^ https://david.ncifcrf.gov/
- ^ http://string.embl.de/
- ^ http://www.cytoscape.org/
References
Avidan, N., Le Panse, R., Berrih-Aknin, S., and Miller, A. (2014). Genetic basis of myasthenia gravis - a comprehensive review. J. Autoimmun. 52, 146–153. doi: 10.1016/j.jaut.2013.12.001
Bernasconi, P., Passerini, L., Annoni, A., Ubiali, F., Marcozzi, C., Confalonieri, P., et al. (2003). Expression of transforming growth factor-beta1 in thymus of myasthenia gravis patients: correlation with pathological abnormalities. Ann. N. Y. Acad. Sci. 998, 278–283. doi: 10.1196/annals.1254.031
Berrih-Aknin, S., and Le Panse, R. (2014). Myasthenia gravis: a comprehensive review of immune dysregulation and etiological mechanisms. J. Autoimmun. 52, 90–100. doi: 10.1016/j.jaut.2013.12.011
Cavalcante, P., Cufi, P., Mantegazza, R., Berrih-Aknin, S., Bernasconi, P., and Le Panse, R. (2013). Etiology of myasthenia gravis: innate immunity signature in pathological thymus. Autoimmun. Rev. 12, 863–874. doi: 10.1016/j.autrev.2013.03.010
Chen, C. Y., Chen, J., He, L., and Stiles, B. L. (2018). PTEN: tumor suppressor and metabolic regulator. Front. Endocrinol. (Lausanne) 9:338. doi: 10.3389/fendo.2018.00338
Chen, J., Liu, A., Wang, Z., Wang, B., Chai, X., Lu, W., et al. (2020). LINC00173.v1 promotes angiogenesis and progression of lung squamous cell carcinoma by sponging miR-511-5p to regulate VEGFA expression. Mol. Cancer 19:98. doi: 10.1186/s12943-020-01217-2
Chen, Y. G., Satpathy, A. T., and Chang, H. Y. (2017). Gene regulation in the immune system by long noncoding RNAs. Nat. Immunol. 18, 962–972. doi: 10.1038/ni.3771
Chou, C. H., Shrestha, S., Yang, C. D., Chang, N. W., Lin, Y. L., Liao, K. W., et al. (2018). miRTarBase update 2018: a resource for experimentally validated microRNA-target interactions. Nucleic Acids Res. 46(D1), D296–D302. doi: 10.1093/nar/gkx1067
Ciccacci, C., Perricone, C., Ceccarelli, F., Rufini, S., Di Fusco, D., Alessandri, C., et al. (2014). A multilocus genetic study in a cohort of Italian SLE patients confirms the association with STAT4 gene and describes a new association with HCP5 gene. PLoS One 9:e111991. doi: 10.1371/journal.pone.0111991
Coppede, F., Stoccoro, A., Nicoli, V., Gallo, R., De Rosa, A., Guida, M., et al. (2020). Investigation of GHSR methylation levels in thymomas from patients with Myasthenia Gravis. Gene 752:144774. doi: 10.1016/j.gene.2020.144774
Dennis, G. Jr., Sherman, B. T., Hosack, D. A., Yang, J., Gao, W., Lane, H. C., et al. (2003). DAVID: database for annotation, visualization, and integrated discovery. Genome Biol. 4:3.
Ding, Y., Li, T., Yan, X., Cui, M., Wang, C., Wang, S., et al. (2021). Identification of hub lncRNA ceRNAs in multiple sclerosis based on ceRNA mechanisms. Mol. Genet. Genomics 296, 423–435. doi: 10.1007/s00438-020-01750-1
Dong, Z. R., Zhou, W., Sun, D., Yan, Y. C., Yang, C. C., Yang, Y. F., et al. (2020). Role of the E3 Ubiquitin Ligase TRIM4 in predicting the prognosis of hepatocellular carcinoma. J. Cancer 11, 4007–4014. doi: 10.7150/jca.37164
Du, Q., Liu, J., Tian, D., Zhang, X., Zhu, J., Qiu, W., et al. (2020). Long noncoding RNA LINC00173 promotes NUTF2 expression through sponging miR-765 and facilitates Tumorigenesis in Glioma. Cancer Manag. Res. 12, 7211–7217. doi: 10.2147/CMAR.S262279
Engelen-Lee, J., Blokhuis, A. M., Spliet, W. G. M., Pasterkamp, R. J., Aronica, E., Demmers, J. A. A., et al. (2017). Proteomic profiling of the spinal cord in ALS: decreased ATP5D levels suggest synaptic dysfunction in ALS pathogenesis. Amyotroph. Lateral Scler. Frontotemporal Degener. 18, 210–220. doi: 10.1080/21678421.2016.1245757
Fan, H., Yuan, J., Li, X., Ma, Y., Wang, X., Xu, B., et al. (2020). LncRNA LINC00173 enhances triple-negative breast cancer progression by suppressing miR-490-3p expression. Biomed. Pharmacother. 125:109987. doi: 10.1016/j.biopha.2020.109987
Fang, T. K., Yan, C. J., and Du, J. (2018). CTLA-4 methylation regulates the pathogenesis of myasthenia gravis and the expression of related cytokines. Medicine (Baltimore) 97:e0620. doi: 10.1097/MD.0000000000010620
Feng, X., Li, W., Song, J., Liu, X., Gu, Y., Yan, C., et al. (2019). HLA typing using next-generation sequencing for Chinese juvenile- and adult-onset myasthenia gravis patients. J. Clin. Neurosci. 59, 179–184. doi: 10.1016/j.jocn.2018.10.077
Gene Ontology, C. (2006). The gene ontology (GO) project in 2006. Nucleic Acids Res. 34, D322–D326. doi: 10.1093/nar/gkj021
Gilhus, N. E., and Verschuuren, J. J. (2015). Myasthenia gravis: subgroup classification and therapeutic strategies. Lancet Neurol. 14, 1023–1036. doi: 10.1016/S1474-4422(15)00145-3
Hehir, M. K., and Silvestri, N. J. (2018). Generalized myasthenia gravis: classification, clinical presentation, natural history, and epidemiology. Neurol. Clin. 36, 253–260. doi: 10.1016/j.ncl.2018.01.002
Jiang, L., Cheng, Z., Qiu, S., Que, Z., Bao, W., Jiang, C., et al. (2012). Altered let-7 expression in Myasthenia gravis and let-7c mediated regulation of IL-10 by directly targeting IL-10 in Jurkat cells. Int. Immunopharmacol. 14, 217–223. doi: 10.1016/j.intimp.2012.07.003
Jonas, S., and Izaurralde, E. (2015). Towards a molecular understanding of microRNA-mediated gene silencing. Nat. Rev. Genet. 16, 421–433. doi: 10.1038/nrg3965
Kaelin, W. G. (2007). Von Hippel-Lindau disease. Annu. Rev. Pathol. 2, 145–173. doi: 10.1146/annurev.pathol.2.010506.092049
Kanehisa, M., and Goto, S. (2000). KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 28, 27–30. doi: 10.1093/nar/28.1.27
Kong, X., Wang, J., Cao, Y., Zhang, H., Lu, X., Wang, Y., et al. (2019). The long noncoding RNA MALAT-1 functions as a competing endogenous RNA to regulate MSL2 expression by sponging miR-338-3p in myasthenia gravis. J. Cell Biochem. 120, 5542–5550. doi: 10.1002/jcb.27838
Lane, L. C., Kus, A., Bednarczuk, T., Bossowski, A., Daroszewski, J., Jurecka-Lubieniecka, B., et al. (2020). An intronic HCP5 variant is associated with age of onset and susceptibility to graves disease in UK and Polish Cohorts. J. Clin. Endocrinol. Metab. 105, e3277–e3284. doi: 10.1210/clinem/dgaa347
Li, J. H., Liu, S., Zhou, H., Qu, L. H., and Yang, J. H. (2014). starBase v2.0: decoding miRNA-ceRNA, miRNA-ncRNA and protein-RNA interaction networks from large-scale CLIP-Seq data. Nucleic Acids Res. 42, D92–D97. doi: 10.1093/nar/gkt1248
Link, J., Navikas, V., Yu, M., Fredrikson, S., Osterman, P. O., and Link, H. (1994). Augmented interferon-gamma, interleukin-4 and transforming growth factor-beta mRNA expression in blood mononuclear cells in myasthenia gravis. J. Neuroimmunol. 51, 185–192. doi: 10.1016/0165-5728(94)90080-9
Liu, Y., Helms, C., Liao, W., Zaba, L. C., Duan, S., Gardner, J., et al. (2008). A genome-wide association study of psoriasis and psoriatic arthritis identifies new disease loci. PLoS Genet. 4:e1000041. doi: 10.1371/journal.pgen.1000041
Luo, M., Liu, X., Meng, H., Xu, L., Li, Y., Li, Z., et al. (2017). IFNA-AS1 regulates CD4(+) T cell activation in myasthenia gravis though HLA-DRB1. Clin. Immunol. 183, 121–131. doi: 10.1016/j.clim.2017.08.008
Mao, A. P., Ishizuka, I. E., Kasal, D. N., Mandal, M., and Bendelac, A. (2017). A shared Runx1-bound Zbtb16 enhancer directs innate and innate-like lymphoid lineage development. Nat. Commun. 8:863. doi: 10.1038/s41467-017-00882-0
Marofi, F., Vahedi, G., Solali, S., Alivand, M., Salarinasab, S., Zadi Heydarabad, M., et al. (2019). Gene expression of TWIST1 and ZBTB16 is regulated by methylation modifications during the osteoblastic differentiation of mesenchymal stem cells. J. Cell Physiol. 234, 6230–6243. doi: 10.1002/jcp.27352
Martin, A., Tegla, C. A., Cudrici, C. D., Kruszewski, A. M., Azimzadeh, P., Boodhoo, D., et al. (2015). Role of SIRT1 in autoimmune demyelination and neurodegeneration. Immunol. Res. 61, 187–197. doi: 10.1007/s12026-014-8557-5
Masunaga, A., Omatsu, M., Kunimura, T., Uematsu, S., Kamio, Y., Kitami, A., et al. (2017). Expression of PTEN and its pseudogene PTENP1, and promoter methylation of PTEN in non-tumourous thymus and thymic tumours. J. Clin. Pathol. 70, 690–696. doi: 10.1136/jclinpath-2016-204220
Meriggioli, M. N., and Sanders, D. B. (2009). Autoimmune myasthenia gravis: emerging clinical and biological heterogeneity. Lancet Neurol. 8, 475–490. doi: 10.1016/S1474-4422(09)70063-8
Paraskevopoulou, M. D., Georgakilas, G., Kostoulas, N., Reczko, M., Maragkakis, M., Dalamagas, T. M., et al. (2013). DIANA-LncBase: experimentally verified and computationally predicted microRNA targets on long non-coding RNAs. Nucleic Acids Res. 41, D239–D245. doi: 10.1093/nar/gks1246
Rafei, H., Kantarjian, H. M., and Jabbour, E. J. (2019). Recent advances in the treatment of acute lymphoblastic leukemia. Leuk. Lymphoma 60, 2606–2621. doi: 10.1080/10428194.2019.1605071
Salmena, L., Poliseno, L., Tay, Y., Kats, L., and Pandolfi, P. P. (2011). A ceRNA hypothesis: the Rosetta Stone of a hidden RNA language? Cell 146, 353–358. doi: 10.1016/j.cell.2011.07.014
Sanders, D. B., Wolfe, G. I., Benatar, M., Evoli, A., Gilhus, N. E., Illa, I., et al. (2016). International consensus guidance for management of myasthenia gravis: executive summary. Neurology 87, 419–425. doi: 10.1212/WNL.0000000000002790
Schultz, M. D., He, Y., Whitaker, J. W., Hariharan, M., Mukamel, E. A., Leung, D., et al. (2015). Human body epigenome maps reveal noncanonical DNA methylation variation. Nature 523, 212–216. doi: 10.1038/nature14465
Sheth, M. N., Nations, S. P., Wolfe, G. I., and Trivedi, J. R. (2005). Von hippel-lindau disease associated with thymoma and myasthenia gravis. J. Clin. Neuromuscul. Dis. 7, 59–61. doi: 10.1097/01.cnd.0000185585.00247.cd
Shi, C. S., Qi, H. Y., Boularan, C., Huang, N. N., Abu-Asab, M., Shelhamer, J. H., et al. (2014). SARS-coronavirus open reading frame-9b suppresses innate immunity by targeting mitochondria and the MAVS/TRAF3/TRAF6 signalosome. J. Immunol. 193, 3080–3089. doi: 10.4049/jimmunol.1303196
Shichrur, K., Feinberg-Gorenshtein, G., Luria, D., Ash, S., Yaniv, I., and Avigad, S. (2014). Potential role of WSB1 isoforms in growth and survival of neuroblastoma cells. Pediatr. Res. 75, 482–486. doi: 10.1038/pr.2014.2
Tryka, K. A., Hao, L., Sturcke, A., Jin, Y., Wang, Z. Y., Ziyabari, L., et al. (2014). NCBI’s database of genotypes and phenotypes: dbGaP. Nucleic Acids Res. 42, D975–D979. doi: 10.1093/nar/gkt1211
Vandiedonck, C., Giraud, M., and Garchon, H. J. (2005). Genetics of autoimmune myasthenia gravis: the multifaceted contribution of the HLA complex. J. Autoimmun. 25, 6–11. doi: 10.1016/j.jaut.2005.09.010
Wang, C. B., Wang, Y., Wang, J. J., and Guo, X. L. (2020). LINC00662 triggers malignant progression of chordoma by the activation of RNF144B via targeting miR-16-5p. Eur. Rev. Med. Pharmacol. Sci. 24, 1007–1022. doi: 10.26355/eurrev_202002_20151
Wang, J., Cao, Y., Lu, X., Wang, X., Kong, X., Bo, C., et al. (2020). Identification of the regulatory role of lncRNA SNHG16 in myasthenia gravis by constructing a competing endogenous RNA network. Mol. Ther. Nucleic Acids 19, 1123–1133. doi: 10.1016/j.omtn.2020.01.005
Wang, P., Li, X., Gao, Y., Guo, Q., Wang, Y., Fang, Y., et al. (2019). LncACTdb 2.0: an updated database of experimentally supported ceRNA interactions curated from low- and high-throughput experiments. Nucleic Acids Res. 47(D1), D121–D127. doi: 10.1093/nar/gky1144
Wang, P., Zhi, H., Zhang, Y., Liu, Y., Zhang, J., Gao, Y., et al. (2015). miRSponge: a manually curated database for experimentally supported miRNA sponges and ceRNAs. Database (Oxford) 2015:bav098. doi: 10.1093/database/bav098
Xu, L., Li, Z., Li, Y., Luo, Z., Luo, Y., Xiao, B., et al. (2020). The expression pattern and regulatory mechanism of the G0/G1 switch gene 2 (G0S2) in the pathogenesis and treatment of AChR Myasthenia Gravis (MG). Mediators Inflamm. 2020:4286047. doi: 10.1155/2020/4286047
Yu, Q., Dong, L., Li, Y., and Liu, G. (2018). SIRT1 and HIF1alpha signaling in metabolism and immune responses. Cancer Lett. 418, 20–26. doi: 10.1016/j.canlet.2017.12.035
Zhang, G., Sun, H., Zhang, Y., Zhao, H., Fan, W., Li, J., et al. (2018). Characterization of dysregulated lncRNA-mRNA network based on ceRNA hypothesis to reveal the occurrence and recurrence of myocardial infarction. Cell Death Discov. 4:35. doi: 10.1038/s41420-018-0036-7
Zhang, H., Guo, Y., Gu, H., Wei, X., Ma, W., Liu, D., et al. (2019). TRIM4 is associated with neural tube defects based on genome-wide DNA methylation analysis. Clin. Epigenetics 11:17. doi: 10.1186/s13148-018-0603-z
Keywords: lncRNA, ceRNA, myasthenia gravis, biomarker, network
Citation: Xu S, Wang T, Lu X, Zhang H, Liu L, Kong X, Li S, Wang X, Gao H, Wang J and Wang L (2021) Identification of LINC00173 in Myasthenia Gravis by Integration Analysis of Aberrantly Methylated- Differentially Expressed Genes and ceRNA Networks. Front. Genet. 12:726751. doi: 10.3389/fgene.2021.726751
Received: 17 June 2021; Accepted: 20 August 2021;
Published: 16 September 2021.
Edited by:
Xiaojian Shao, National Research Council Canada (NRC-CNRC), CanadaReviewed by:
Yaqoub Ashhab, Palestine Polytechnic University, PalestineFabricio Alves Barbosa da Silva, Oswaldo Cruz Foundation (Fiocruz), Brazil
Copyright © 2021 Xu, Wang, Lu, Zhang, Liu, Kong, Li, Wang, Gao, Wang and Wang. 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: Jianjian Wang, d2FuZ2ppYW5fNDI3QDE2My5jb20=; Lihua Wang, d2FuZ2xoMjExQDE2My5jb20=