- 1Department of Clinical Laboratory, Shanghai Tenth People’s Hospital, Tongji University, Shanghai, China
- 2Department of Pathology, Shanghai Tenth People’s Hospital, Tongji University, Shanghai, China
- 3Department of Pathology, Universal College of Medical Sciences, Bhairahawa, Nepal
- 4Anhui Medical University Shanghai Clinical College, Hefei, China
Endometriosis (EM) is a chronic neuroinflammatory disorder that is associated with pain and infertility that affects ∼10% of reproductive-age women. The pathophysiology and etiology of EM remain poorly understood, and diagnostic delays are common. Exploration of the underlying molecular mechanism, as well as novel diagnostic biomarkers and therapeutic targets, is urgently needed. Inflammation is known to play a key role in the development of lesions, which are a defining feature of the disorder. In our research, the CIBERSORT and WGCNA algorithms were used to establish a weighted gene co-expression network and to identify macrophage-related hub genes using data downloaded from the GEO database (GSE11691, 7305). The analysis identified 1,157 differentially expressed genes (DEGs) in EM lesions, of which five were identified as being related to M2 macrophages and were validated as differentially expressed by qRT-PCR and immunohistochemistry (IHC). Of these putative novel biomarker genes, bridging integrator 2 (BIN2), chemokine receptor 5 (CCR5), and macrophage mannose receptor 1 (MRC1) were upregulated, while spleen tyrosine kinase (SYK) and metalloproteinase 12 (ADAM12) were downregulated in ectopic endometria vs. normal endometria. Meanwhile, 23 potentially therapeutic small molecules for EM were obtained from the cMAP database, among which topiramate, isoflupredone, adiphenine, dexverapamil, MS-275, and celastrol were the top six molecules with the highest absolute enrichment values. This is our first attempt to use the CIBERSORT and WGCNA algorithms for the identification of novel Mϕ2 macrophage-related biomarkers of EM. Our findings provide novel insights into the impact of immune cells on the etiology of EM; nevertheless, further investigation of these key genes and therapeutic drugs is needed to validate their effects on EM.
Introduction
Endometriosis (EM) is a chronic, estrogen-dependent inflammatory disease that is characterized by abnormal growth of endometrial tissue outside the uterine cavity (Zannoni et al., 2016). High prevalence rates (10–15%) of endometriosis are reported in women between 20 and 40 years of age (Hickey et al., 2014). The causes and pathophysiology of endometriosis are still unclear. EM is mainly diagnosed by laparoscopy (Rolla, 2019). Patient treatments involve both drug therapy and laparoscopy surgery. The application of therapeutic agents is usually suppressive rather than curative, and the most commonly used progestogens that exhibit demonstrated efficacy are not recommended for long-term use because of their extensive side effects (Rolla, 2019). Disease recurrence is very common soon after the discontinuation of treatment with medicines (Falcone and Flyckt, 2018). Therefore, complete laparoscopic surgical resection of ectopic endometrial tissues is the most common treatment (Falcone and Flyckt, 2018). However, the progression of disease is very slow, and often takes 7–10 years before the onset of symptoms, which results in delays in diagnosis and optimal treatment timing (Bakhtiarizadeh et al., 2018). Thus, additional serious research work must be carried out to explore the underlying molecular mechanisms of EM. Similarly, novel diagnostic biomarkers and therapeutic targets should be identified to improve the early diagnosis and treatment of EM, which can improve the quality of life and cure rates of patients (Fassbender et al., 2015).
Various studies have established that the immune system plays a major role in the pathophysiology and symptomatology of EM, and immune cells, including natural killer (NK) cells, macrophages, neutrophils, and the CD4 T helper cell system, are dysregulated in women with EM (Beste et al., 2014; Jeung et al., 2016; Patel et al., 2018; Zondervan et al., 2018). Similarly, a significant increase in macrophages (Mϕs) allows cytotoxic T helper cells to release excessive amounts of inflammatory cytokines, thus facilitating the creation of a pro-inflammatory environment in the endometrium and eventually the development of endometritis (Vallvé-Juanico et al., 2019a). In addition, the polarization/distribution of pro-inflammatory Mϕ type 1 and anti-inflammatory Mϕ type 2 macrophages also has a significant role in the development of EM (Hudson et al., 2020). Since the polarization state of Mϕs is influenced by the microenvironment, there is a significant reduction in pro-inflammatory Mϕ type 1 cells and a significant increase in anti-inflammatory Mϕ type 2 cells in the EM microenvironment. Aberrant increase and activation of anti-inflammatory Mϕ type 2 could stimulate the abnormal gene expression that is associated with the ectopic endometrium (Vallvé-Juanico et al., 2019b). Therefore, recognition of potential biomarkers associated with Mϕ infiltration can help to understand and explore their roles in immune pathogenesis in EM, which can contribute to the management and treatment of endometriosis patients.
Over the last few years, with the advancement and availability of numerous online integrated bioinformatics tools, the identification and establishment of distinct molecular markers and signaling pathways for different diseases have become easier (Yang et al., 2018; Fachal et al., 2020). Additionally, several bioinformatic studies on EM have identified various genes involved in EM, such as CXCL12, MALAT1, CC2D2A, AEBP1, HOXB6, and IER3 (Bakhtiarizadeh et al., 2018; Chen et al., 2020; Hu et al., 2020). Weighted gene co-expression network analysis (WGCNA) is the most valuable and widely used tool for establishing co-expression gene networks and identifying important hub gene modules for cancers (Langfelder and Horvath, 2008). In this study, we used the normal endometrium and ectopic endometrium-related gene expression datasets, GSE11691 and GSE7305, which were extracted from the Gene Omnibus (GEO) database, to identify potential macrophage-related biomarkers in EM using WGCNA. Similarly, CIBERSORT, which has been successfully applied for estimating the infiltration of immune cells in prostate cancer (Zhao et al., 2019) and kidney cancer (Zhang et al., 2019), was used to determine the extent of distinct immune cells, including macrophage infiltration in EM. In addition, the most significant modules related to macrophage infiltration were identified, and the prognostic genes in these modules were further analyzed and verified. Moreover, with the application of the cMAP database, we identified the small molecules that may have an effect on EM by using differentially expressed genes (DEGs). This is the first time that WGCNA was used to explore macrophage-related genes in EM, which could provide novel insights for the early diagnosis at the molecular level and treatment of patients with EM.
Results
Data Processing Strategy and mRNA Expression Profiles
Figure 1 illustrates the study protocol. The mRNA expression profiles (Hull et al., 2008; Falcone and Flyckt, 2018) of 19 pairs of endometriosis and healthy endometrial samples in the GSE11691 and GSE7305 datasets were downloaded from the GEO database, which were then merged and batch-normalized for the following analysis.
FIGURE 1. Workflow for identification of the novel macrophage M2-related biomarkers through the CIBERSORT and WGCNA algorithms, and exploration of potential therapeutic agents in endometriosis by using data downloaded from the GEO database (CV: distinct variant genes expressed in 19 endometriosis samples with coefficient values > 0.1; DEGs: differentially expressed genes, ROC curve: receiver operating characteristic curve, IHC: immunohistochemistry, and RT-PCR: reverse transcription polymerase chain reaction).
Estimation of the Immune-Infiltration Level in Endometriosis
CIBERSORT algorithm analyzes the abundance of the 22 distinct types of immune cell infiltration in 19 pairs of samples. Three distinct types of macrophages in EM samples were selected as trait data for WGCNA (Supplementary Table S1).
Construction of the Co-Expression Network and Identification of Hub Modules Related to Macrophages Using Weighted Gene Co-Expression Network Analysis
A total of 2,395 genes with correlation coefficients higher than 0.1 (CV genes) (Supplementary Table S2) were identified and used for construction of gene co-expression network analysis using the WGCNA R package. A scale-independent topological network (soft-thresholding power nine scale-free R2 = 0.9) and mean connectivity network were established (Figures 2A,B). Dynamic hybrid cutting was enforced for the construction of a hierarchical clustering tree by splitting the dendrogram at relevant transition points. Single genes are presented as the leaves of trees, and multiple genes with analogical expression data are presented as branches of the dendrogram tree. Branches that contained analogously expressed genes were considered to be gene modules. A total of 16 gene modules were generated with the fusion of similar modules (Figure 2C).
FIGURE 2. Hierarchical cluster formation based on the soft threshold power (β). (A) Estimation of the scale independence index of the 1–20 soft threshold power (β = 9). (B) Determination of the mean connectivity of the 1–20 soft threshold power. (C) Identification and establishment of gene co-expression modules that are shown as different colors in hierarchical clustering. (D) Relationships of consensus module eigengenes and Mϕ. The rows in the figure correspond to consensus modules, and the columns correspond to macrophage subtypes. The numbers in each module indicated the correlation coefficients to show the association between the corresponding module and Mϕ, along with the p values shown below in parentheses (red indicates a positive correlation and green indicates a negative correlation).
Among these 16 modules, the correlations (R2) of the yellow and red modules with those of Mϕ2 were 0.67 and 0.59, respectively. Meanwhile, the correlation (R2) between the purple module and Mϕ1 was 0.54. However, the correlations of other modules with Mϕ were <0.5 (Figure 2D). Therefore, the yellow module that showed the highest connectivity with Mϕ2 (R2 = 0.67, p = 2.2e−32) was identified as a hub module (Figures 2D, 3A). A total of 36 hub genes with module membership values > 0.8 and gene significance values > 0.5 were discerned from the yellow hub module (Supplementary Table S3).
FIGURE 3. Enrichment analysis and PPI construction of hub genes in key modules. (A) Module membership in yellow modules and gene significance for Mϕ2. Each yellow dot indicates a gene. Genes with module membership values >0.8 and gene significances >0.5 inside the green box are candidate hub genes. (B) GO analysis of hub genes: biological process (BP), cellular composition (CC), and molecular function (MF) of hub genes. (C) KEGG analysis of hub genes identified the top 18 pathways, dot sizes reflect the number of genes associated with relative pathways, and dot colors indicate the p-values. (D) Protein–protein interaction (PPI) network of hub genes was constructed using the STRING database. (E) Histogram based on the number of protein–protein acting nodes shown in the network (only the top 30 are shown).
Enrichment and Protein–Protein Interaction Analysis of the Hub Gene
The biological processes (BP), cellular components (CC), and molecular functions (MF) of the hub genes were identified and explored using Gene Ontology (GO) analysis. In Figure 3B, the top 10 enrichment terminologies related to BF, CC, and MF are shown. The top three enriched terms of BP were related to the immune response: T-cell activation, positive regulation of cell activation, and regulation of tumor necrosis factor (TNF) superfamily cytokine production. Similarly, the external side of the plasma membrane, secretory granule membrane, and endocytic vesicle were the top three highly enriched CC terminology fits on GO analysis. Meanwhile, immune receptor activity, immunoglobulin binding, and MHC protein binding were the top three highly enriched terminologies related to MF. Details of the GO analysis reports are available in the supplementary materials (Supplementary Table S4).
Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of hub genes identified 18 immune-related pathways, including the chemokine signaling pathway, NK cell-mediated cytotoxicity, phagocytosis, B-cell receptor signaling pathway, and primary immunodeficiency (Figure 3C). Details of the KEGG analysis are included in the supplementary materials (Supplementary Table S5).
Protein–protein interaction (PPI) network analysis of EM-related hub genes using Cytoscape STRING software constructed a PPI network that consisted of 36 genes, 36 nodes, and 208 edges (Figures 3D,E). All 36 hub genes were selected for subsequent analysis (Supplementary Table S3).
Identification and Exploration of Diagnostic Biomarkers
To identify novel EM-related diagnostic biomarkers, 1,157 differentially expressed genes (DEGs) with log2FC > 0.5 and false discovery rate–adjusted p values (FDR) < 0.05 were screened. Of the 1,157 DEGs, 568 genes were upregulated (Supplementary Table S6), and 589 genes were downregulated (Supplementary Table S7). The 1,157 DEGs are shown in the volcano plots and heat maps (Figures 4A,B). GO analysis of 1,157 DEGs showed that they could regulate distinct biological processes (BPs), which were mainly involved in connective tissue development, embryonic organ development, endothelial cell proliferation, epithelial cell proliferation, and so on (Figure 4C). KEGG pathway analysis of DEGs showed their involvement in the progression of EM with regulation of distinct PI3K−Akt and Rap1 signaling pathways (Figure 4D), which was consistent with other studies (Shigetomi et al., 2012; Tsai et al., 2018; Madanes et al., 2020).
FIGURE 4. DEGs enrichment analysis. (A) Volcano plots of DEGs between endometrium and endometriosis, where red dots indicate upregulated genes, green dots indicate downregulated genes, and gray dots indicate no substantial difference (|log2 (fold change)|>0.5 and corrected p-values <0.05). (B) Heat map illustrating the DEGs in both the endometrium and endometriosis samples, where red indicates significantly upregulated genes, green indicates significantly downregulated genes, and white indicates genes with no substantial difference. (C) GO analysis of DEGs. The center dots indicate the different important biological processes, and the sizes of the center dots reflect the number of genes involved in biological processes. The small dots that radiate from the central dot connected by lines with the same color represent those genes with the same biological function, and the colors of the small dots represent the fold change of the genes that they reflect. (D) KEGG pathway analysis of DEGs revealed the top 11 relevant pathways. The dot sizes refer to the number of genes involved in the pathway, and the dot colors depict the p-values.
Subsequent analysis of 1,157 DEGs and hub genes using a Venn diagram identified five overlapping genes (e.g., SYK, BIN2, ADAM12, CCR5, and MRC1) (Figure 5A) as potential key genes. The differential expressions of SYK, BIN2, ADAM12, CCR5, and MRC1 between the endometrium and endometriosis are shown in Figure 5B, with significant p-values of 0.005, <0.001, <0.001, 0.002, and 0.001, respectively. Among the five key genes, SYK and ADAM12 were downregulated, while the three other key genes were upregulated in endometriosis.
FIGURE 5. Identification and validation of the diagnostic values of the key genes. (A) Venn diagram of DEGs and hub genes. Commonly intersecting genes were identified as key genes, which included SYK, BIN2, ADAM12, CCR5, and MRC1. (B) Differential expression of five key genes between the eutopic endometrium and ectopic endometrium. (C) ROC curves for SYK, BIN2, ADAM12, CCR5, and MRC1.
Similarly, to validate the diagnostic significance of all five key genes, ROC analysis was carried out (Figure 5C). The AUC values of all genes were greater than 0.75, of which MRC1 (AUC = 0.8864, SE = 0.06, p < 0.0001) had the maximum AUC value and SYK (AUC = 0.7645, SE = 0.08, p = 0.005) had the minimum AUC value.
qRT-PCR and Immunohistochemical Analysis
qRT-PCR and IHC were carried out on 10 pairs of matched eutopic and ectopic endometrial tissue samples to validate the expression levels of five potential diagnostic biomarkers. qRT-PCR was carried out using the total RNA that was extracted from 10 pairs of normal endometrium and EM tissues and showed higher expressions of BIN2, CCR5, and MRC1 and lower expressions of SYK and ADAM12 in endometriosis tissue than in non-endometriosis tissue (Figure 6A).
FIGURE 6. Results of qRT-PCR and IHC analysis. (A) Relative expressions of five key genes in the eutopic endometrium and ectopic endometrium, as determined by qRT-PCR analysis. (*p < 0.05, **p < 0.01, ***p < 0.005, ****p < 0.001). (B) Relative protein expressions of five key genes in the eutopic endometrium and ectopic endometrium, as determined by IHC techniques.
Meanwhile, IHC staining of EM tissues using anti-SYK, anti-BIN2, anti-ADAM12, anti-CCR5, and anti-MRC1 showed higher expressions of BIN2, CCR5, and MRC1 and lower expressions of SYK and ADAM12 relative to normal endometrium tissues (Figure 6B).
In addition, the Turku Endometriosis Database was used to validate the expressions of five key genes in various tissues of EM patients. The expression trend of the five key genes was consistent with the findings in this study (Supplementary Figure S1). Moreover, the expression of these genes in the peritoneal fluid and other lesions also presented discrepancies (Supplementary Figure S1). Therefore, to verify the above results, these five key genes have the potential to be diagnostic biomarkers for EM, although they need to be further validated.
Identification of Small Molecular Therapeutic Agents
With the application of the cMAP database, 23 potentially important small molecules targeting EM were screened (Supplementary Table S8). As shown in Table 1, the six small molecules with the highest absolute enrichment values were chosen, which indicated significant correlations with EM and were as follows: topiramate (enrichment = −0.972, p = 0.0000), isoflupredone (enrichment = −0.948, p = 0.00026), adiphenine (enrichment = −0.936, p = 0.0000), dexverapamil (enrichment = 0.947, p = 0.0054), MS-275 (enrichment = 0.961, p = 0.00262), and celastrol (enrichment = 0.993, p = 0.0000). Figure 7 shows the 3D conformers of the above six chemicals that presumably can reverse or induce specific gene expression, hence affecting the state of EM.
FIGURE 7. 3D conformers of the top six candidate therapeutic agents with the highest absolute enrichments. For the identification of small molecules, the DEGs and Broad Institutes Connectivity Map (cMAP) database (https://portals.broadinstitute.org/cmap) were used, and the connectivity cutoff value used was >0.8. Molecules with the highest absolute enrichment values, namely, those that are strongly associated with the DEGs, indicate their potential therapeutic effects on EM.
Materials and Methods
Acquisition and Processing Gene Expression Data
The GSE7305 and GSE11691 microarray data (raw data CEL files) were downloaded from the GEO database (http://www.ncbi.nlm.nih.gov/geo/). The GSEA7305 dataset, which was extracted on the GPL570 platform (Affymetrix U133 Plus 2.0 Array), comprises 10 ovarian endometriosis samples and 10 matched control endometrium samples from the same patients (follicular phase, n = 2; luteal phase, n = 8) (Hever et al., 2007). Meanwhile, the GSE11691 dataset, which was derived on the GPL96 platform (Affymetrix Human Genome U133A Array), consists of nine eutopic endometria and nine matched ectopic lesion endometrium from either the broad ligament (visceral peritoneum) or parietal peritoneum of nine 20- to 46-year-old women with endometriosis (AFS-r stages 2–4; proliferative phase, n = 5; secretory phase, n = 4) (Hull et al., 2008). By using the RMA (robust multichip averaging) algorithm (Irizarry et al., 2003), the GSEA7305 and GSE11691 data were initially modified, and they were combined into an integrated dataset by “Perl.” Integrated datasets were eventually batch-normalized using the R packages “sva” and “limma” (Ritchie et al., 2015).
Estimation of Extent of Immune Cell Infiltration in Endometriosis
In this study, we applied the CIBERSORT algorithm to determine the immune cell subsets of GSE705 and GSE1169 (Newman et al., 2015; Lin et al., 2020). The CIBERSORT algorithm is an excellent method in contrast to traditional deconvolution methods for estimating infiltrating immunity since it can analyze unspecified data and noise, which makes it an excellent tool for calculating the abundances of specialized cells within the mixed matrix.
Establishment of Co-Expression Networks and Identification of Macrophage-Related Hub Modules Using Weighted Gene Co-Expression Network Analysis
To establish the co-expression network, we initially identified the distinct variant genes that were expressed in 19 endometriosis samples. A total of 2,395 genes with coefficient values > 0.1 were selected (Lin et al., 2020). Similarly, we applied the “WGCNA” R package to establish a weighted gene co-expression network of 2395 CV genes, which not only reduced the size of computation of the whole network but also maintained a scale-free topological network (Langfelder and Horvath, 2008). Then, Pearson’s correlation matrices were determined with the similarity matrix converted from the expression of individual transcripts. The similarity matrix was later converted to an adjacency matrix on the basis of this equation: amn = |cmn|β, where amn is the adjacency between paired genes, cmn is Pearson’s correlation coefficient between the paired genes, and β is the soft-power threshold. Next, a topological overlap matrix (TOM) was established from the adjacency matrix when β = 9. The value of β was also used to evaluate the connectivity characteristics in the co-expression network. In addition, the average linkage hierarchical clustering and dynamic hybrid cutting method were used for the construction of a dendrogram of the TOM matrix and to categorize the genes with identical expression patterns into distinct modules. A bottom-up algorithm with a module minimum size cutoff value of 30 was set, and for merging identical modules, a threshold of 0.25 was set. Component analysis of each module was carried out using module eigengenes. For identification and exploration of important modules, the correlations among the infiltration levels of Mϕ and the modules were calculated by Pearson’s tests. A distinct module with p < 0.05 was considered to be significantly correlated with Mϕ infiltration. The Mϕ subset and modules with the maximum correlation coefficients were considered to be hub modules.
Enrichment and Protein–Protein Interaction Analysis of Hub Genes
In the hub module, genes with gene significance values > 0.5 and module membership values > 0.8 were described as hub genes. All hub genes were subjected to GO, KEGG, and PPI analyses.
GO annotation contained three parts: BP, CC, and MF . KEGG analysis was carried out to explore the interactive networks of hub genes, which thereby facilitated the retrieval of their genetic information. For the GO and KEGG analyses, the hub gene symbols were initially converted to gene IDs using the R package “org.Hs.eg.db.” Similarly, the “clusterProfiler,” “org.Hs.eg.db,” “enrichplot,” and “ggplot2” R packages were used to determine the GO biological processes and KEGG pathway analysis of the target genes using gene ID. Eventually, the gene IDs used in the KEGG pathway were retrieved and again converted into the Symbol gene by “Perl.” A threshold of FDR <0.05 was utilized as the cutoff criterion.
We performed PPI network analysis using the STRING database (https://string-db.org/) to explore the interactions among proteins that were encoded by the hub genes (Szklarczyk et al., 2019). For the PPI network analysis, all hub genes were subjected to the STRING database. The minimum interaction score was set >0.4 to establish the PPI network. The disconnected nodes in the PPI network were hidden. R software was used to construct a histogram of the top 30 genes with the highest number of nodes.
Identification of Novel Diagnostic Markers for Endometriosis
Initially, the merged and batch-normalized gene expression data transformed with log2 were used to screen the DEGs with log (fold-change) values ≥ 0.5 and p-values< 0.05 by the “limma” R package. A total of 1,157 DEGs were identified and subjected to GO and KEGG analyses in the same way as the hub genes.
Next, the DEGs and hub genes were subjected to an online Web-based Venn diagram tool (http://bioinformatics.psb.ugent.be/webtools/Venn/) for the identification of key genes. In this study, intersecting DEGs and hub genes were considered to be key genes. In addition, the differential expressions of key genes in the endometrium and endometriosis were analyzed using the R packages “limma” and “beeswarm.”
Next, ROC analysis was conducted with the key genes using GraphPad Prism 6.01 to validate their diagnostic value.
qRT-PCR and Immunohistochemical Staining Analysis
qRT-PCR and IHC analysis were carried out using 10 ectopic endometrial samples from five broad ligament endometriosis (follicular phase, n = 2; luteal phase, n = 3) and five ovarian endometriosis (follicular phase, n = 1; luteal phase, n = 4) patients who underwent intrauterine surgery at Shanghai Tenth People’s Hospital from May 2019 to July 2020. All endometriosis patients were diagnosed with stage IV endometriosis according to the staging criteria of endometriosis that is defined by American Fertility Society revised (AFS-r) classification. Meanwhile, 10 eutopic endometrial tissues from the above patients were selected as the control group. The ages of the patients who contributed to the tissues were between 20 and 40 years. The study was permitted by the Ethics Committee of Shanghai Tenth People’s Hospital.
Total RNA from the ectopic and eutopic endometrial tissues was extracted using the TRIzol reagent (Invitrogen, Thermo Scientific, United States) and was then reverse-transcribed into cDNA using the Primescript reverse transcription reagent kit (Takara, Dalian, China). The relative expressions of five key genes were determined by qRT-PCR using the SYBR Green reagent kit (Takara, Otsu, Shiga, Japan) in an ABI 7500 PCR system (Applied Biosystems). The primer sequences of SYK, BIN2, ADAM12, CCR5, and MRC1 were as follows: SYK-F: 5′-CCTGGCGCAGGTGGAC-3′ and SYK-R: 5′-GTAGTGGTGTGCCTTCCTCC-3′, BIN2-F: 5′-TACCTGGAGCAAGAAACGCA-3′ and BIN2-R: 5′- AGGTCCTTGTACAGCTTGTGG-3, ADAM12-F: 5′-CGCTCGAAATTACACGGGTC-3′ and ADAM12-R: 5′-CAGCGAGGTTTGGTGTGTTG-3′, CCR5-F: 5′-TCCAGTGAGAAAAGCCCGTA-3′ and CCR5-R: 5′- GGAACGGATGTCTCAGCTCT-3′, and MRC1-F: 5′-ACCTGCGACAGTAAACGAGG-3′ and MRC1-R: 5′- TGTCCGCTTCATCAGCCATT-3′. GAPDH was used as the endogenous control.
For IHC staining, the paraffin-embedded endometrial tissues were cut into 3-µm-thick and 3-mm-diameter thin sliced sections. Paraffin-embedded thin sections of tissue were deparaffinized in xylene and hydrated through increasing alcohol concentrations. Antigen retrieval from the tissue sections was carried out with the subsequent treatment of tissue sections with sodium citrate (pH = 6) and hydrogen peroxide for 15 min, and washing was performed with PBS. The tissue sections were incubated in 5% BSA for 2 h. Subsequently, the tissue sections were treated with primary anti-SYK, anti-BIN2, anti-ADAM12, anti-CCR5, and anti-MRC1 antibodies (all antibodies were purchased from Abcam) and incubated overnight at 4°C. The sections were incubated with biotinylated secondary rabbit antibodies that were conjugated with streptavidin-HRP (1:200) for 1 h at room temperature, stained with 3,3′-diaminobenzidine solution (DAB Sigma, United States), and finally counterstained with hematoxylin stain. The stained tissue sections were visualized and photographed using a Lecia microscope (Japan). The protein expressions in the stained tissue sections were analyzed using Image-Pro Plus 6.0 software (Media Cybernetics, United States).
Eventually, the expressions of key genes were verified using the Turku Endometriosis Database (https://endometdb.utu.fi/), which contained gene expression data from 1811 different EM tissue samples, 242 patients, and controls.
Identification of Small Molecular Therapeutic Agents
The Broad Institutes Connectivity Map (cMAP) database (https://portals.broadinstitute.org/cmap) was used for the identification of small candidate molecules related to EM (Lamb et al., 2006). For the identification of small candidate chemical molecules, the DEGs were introduced into the cMAP database for gene set enrichment analysis, and the connectivity cutoff value was set >0.8. Small molecules with high absolute enrichment values can possibly reverse or induce biological states that are encoded in specific gene expression markers to have potential therapeutic effects on EM. Clue command (https://clue.io/command) and PubChem (https://pubchem.ncbi.nml.gov) were used for the extraction of detailed information and 3D confirmation of the established small molecules.
Statistical Analysis
All bioinformatics analyses were carried out using R 4.0.2 software or Perl. GraphPad Prism 6.01 was used for the ROC analyses. Image-Pro Plus 6.0 (Media Cybernetics, United States) was used for the analysis and interpretation of ICH findings. The Kruskal–Wallis test (SPSS 18.0) was used for the assessment of multiple independent variables. The findings of the study with p values < 0.05 were considered to be statistically significant.
Discussion
EM is a chronic, estrogen-dependent inflammatory disease that primarily affects women of reproductive age. The diagnosis and treatment of the disease primarily depends upon laparoscopy. The progression and diagnosis of EM is very slow, which results in delayed golden treatment periods (Rolla, 2019). Several studies have reported the dysregulation of distinct immune cells, including Mϕs, in EM, which illustrated the crucial role of Mϕs in disease progression (Hudson et al., 2020; Jeung et al., 2016; Vallvé-Juanico et al., 2019a; Vallvé-Juanico et al., 2019b). In this study, 36 hub genes that were correlated with Mϕ2 infiltration were identified, and the potential roles of these genes in EM progression were further explored. Among the 36 hub genes, only the five hub genes, namely, SYK, BIN2, ADAM12, CCR5, and MRC1, have a potential diagnostic value and are recognized as potential novel prognostic biomarkers in EM. In addition, we identified 23 small molecules that could affect EM progression.
For constructing distinct modules using WGCNA, 2395 CV genes and Mϕ infiltration levels in ectopic endometria were selected. The correlation coefficients were determined for the identification of highly significant Mϕ-associated gene modules. Genes in the hub modules with MM values > 0.8 and gene significance values > 0.5 were considered to be hub genes. The GO and KEGG analyses of the hub genes suggested that they are strongly immune-related genes. Venn diagram analysis of hub genes and DEGs identified the five overlapping genes that were recognized as potential key genes. Of the five key genes, SYK and ADAM12 were downregulated, whereas BIN2, CCR5, and MRC1 were upregulated in ectopic endometria compared with normal endometria. Subsequently, ROC analysis of five individual key genes was carried out, and the AUC area was estimated. All of the AUCs of the five key genes were >0.75, which suggested that they possess the diagnostic value.
Moreover, the relative expressions of five distinct key genes and proteins in endometriosis and non-endometriosis tissues were determined by qRT-PCR and IHC, respectively. In addition, we used the Turku Endometriosis Database for the validation of the expression genes and noted that the expression trends of the five genes between the patient endometria and control endometria were consistent with our results, although it was not clear whether these differences were statistically significant. Considering the above findings, these five key genes have a significant role in EM progression and diagnosis.
SYK is a potent signaling mediator and plays a crucial role in signal transduction in distinct cells. Accumulating evidence on SYK suggests its potential role in the development of various allergic conditions, autoimmune diseases, and cancer (Polak et al., 2020). However, the roles of SYK in EM have not yet been explored. ADAM12 is a multi-domain transmembrane and secreted protein. It releases and activates numerous biologically significant ligands, including tumor necrosis factor (TNF), epidermal growth factor (EGF), transforming growth factor-a (TGF-a), amphiregulin, and HB-EGF, which were found to be associated with the progression of distinct diseases (e.g., liver fibrosis, obesity, and asthma) and cancers (Kveiborg et al., 2005; Rocks et al., 2008; Roy and Moses, 2012; Nyren-Erickson et al., 2013). In addition, numerous studies on ADAM12 have stated its significant role in the pathogenesis of EM (Burney et al., 2007; Painter et al., 2011). ADM12 might be a potential target for the treatment of various diseases and malignancies.
BIN2 gene is a cytoplasmic protein–encoding gene that can interact with the cellular membrane and cytoskeleton, and thus affects podosome formation, motility, and phagocytosis (Sánchez-Barrena et al., 2012). Its role in cancer progression is unknown. However, researchers have used TCGA to establish and explore its involvement in the progression of diseases, including endometrial, cervical, and ovarian cancers (Cancer Genome Atlas, 2012).
CCR5 is a reported receptor for RANTES (regulation-activation, normal T-cell-expressed and T-cell-secreted). It not only regulates leucocyte activation and migration but also plays a significant role in endometrial, stromal, and glandular cellular apoptosis (Antinolo, 2004). Moreover, the expression of CCR5 has been found to be significantly increased in women with EM (Hornung et al., 2001).
MRC1 is highly expressed in macrophages such as Mϕ2. The MRC1 expression on macrophages significantly increased in distinct inflammatory conditions, such as EM, suggesting its significant role in endometriosis development (Woo et al., 2017).
In addition, as there are no effective therapeutic drugs for EM (Falcone and Flyckt, 2018; Rolla, 2019), distinct small molecules that could be effective in EM were identified with DEGs using the cAMP database in this study. Celastrol, MS-275, and dexverapamil are the top three small molecules with higher enrichment correlations.
Celastrol, a naturally extracted compound with anti-inflammatory and antioxidant properties, also acts as an HSP90 (Trott et al., 2008), NF-κB pathway (Kang et al., 2013), and topoisomerase inhibitor. Therefore, it can be used in the treatment of inflammation and cancers (Lee et al., 2012). MS-27 is an HDAC inhibitor and cell cycle inhibitor that has also shown potent antitumor activity (Tang et al., 2020). Cancer cells can develop resistance to multiple therapeutic drugs, which is the major challenge for the effective treatment of cancer patients. However, dexverapamil, a calcium channel blocker, modulates multidrug resistance in tumors (Keilhauer et al., 1994). Although EM is a mild slowly progressive endometrial disease, it exhibits tumor-identical distinguishing characteristics, including lump formation, metastasis, invasion, and relapse. Therefore, the discovered small molecules might have a significant effect on the treatment of EM patients.
Conclusion
This is our first attempt to use the CIBERSORT and WGCNA algorithms for identifying the novel Mϕ2 macrophage-related biomarkers of endometriosis. Five potential diagnostic biomarkers and therapeutic target genes were identified through various validations and experiments, among which SYK and ADAM12 were downregulated genes, while BIN2, CCR5, and MRC1 were upregulated genes in EM. In addition, celastrol, MS-275, and dexverapamil were analyzed as potential therapeutic medicines in the treatment of endometriosis patients. Our findings provide novel insight into EM at the immune and molecular levels; nevertheless, further investigations of these key genes and therapeutic drugs are needed to validate their effect on EM.
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: GSE7305 and GSE11691 microarray data were downloaded from the GEO database (http://www.ncbi.nlm.nih.gov/geo/).
Ethics Statement
The studies involving human participants were reviewed and approved by the Ethical Institutions of Pathology of Shanghai Tenth People’s Hospital. Written informed consent for participation was not required for this study in accordance with the national legislation and the institutional requirements.
Author Contributions
ZC and RB: data analysis and original draft writing. QL and ML: verification experiment implementation. LZ and MZ: helped to collect endometrium samples and proofread the references. FS: edited the manuscript and provided funding support. LF and SZ: study design, supervision, and funding support. All authors contributed to the article and agreed with the submission.
Funding
The authors’ research was funded by the Shanghai Municipal Health Commission (No. 20204Y0088), Cultivation Foundation of Shanghai Tenth People’s Hospital (No. 04.03.19.155), and the National Natural Science Foundation of China (Nos. 81930066 and 81772941).
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.
Acknowledgments
We thank the ethical institutions and department of pathology of Shanghai Tenth People’s Hospital for their cooperation, as well as the financial support institutions and all authors who contributed to this research.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmolb.2021.656145/full#supplementary-material
References
Antinolo, G. (2004). Analysis of the Involvement of CCR5- 32 and CCR2-V64I Variants in the Development of Endometriosis. Mol. Hum. Reprod. 10, 155–157. doi:10.1093/molehr/gah026
Bakhtiarizadeh, M. R., Hosseinpour, B., Shahhoseini, M., Korte, A., and Gifani, P. (2018). Weighted Gene Co-expression Network Analysis of Endometriosis and Identification of Functional Modules Associated With its Main Hallmarks. Front. Genet. 9, 453. doi:10.3389/fgene.2018.00453
Beste, M. T., Pfäffle-Doyle, N., Prentice, E. A., Morris, S. N., Lauffenburger, D. A., Isaacson, K. B., et al. (2014). Molecular Network Analysis of Endometriosis Reveals a Role for C-Jun-Regulated Macrophage Activation. Sci. translational Med. 6, 222ra16. doi:10.1126/scitranslmed.3007988
Burney, R. O., Talbi, S., Hamilton, A. E., Vo, K. C., Nyegaard, M., Nezhat, C. R., et al. (2007). Gene Expression Analysis of Endometrium Reveals Progesterone Resistance and Candidate Susceptibility Genes in Women with Endometriosis. Endocrinology 148, 3814–3826. doi:10.1210/en.2006-1692
Cancer Genome Atlas (2012). Comprehensive Molecular Portraits of Human Breast Tumours. Nature 490, 61–70. doi:10.1038/nature11412
Chen, M., Zhou, Y., Xu, H., Hill, C., Ewing, R. M., He, D., et al. (2020). Bioinformatic Analysis Reveals the Importance of Epithelial-Mesenchymal Transition in the Development of Endometriosis. Sci. Rep. 10, 8442. doi:10.1038/s41598-020-65606-9
Fachal, L., Aschard, H., Beesley, J., Barnes, D. R., Allen, J., Kar, S., et al. (2020). Fine-mapping of 150 Breast Cancer Risk Regions Identifies 191 Likely Target Genes. Nat. Genet. 52, 56–73. doi:10.1038/s41588-019-0537-1
Falcone, T., and Flyckt, R. (2018). Clinical Management of Endometriosis. Obstet. Gynecol. 131, 557–571. doi:10.1097/aog.0000000000002469
Fassbender, A., Burney, R. O., O, D. F., D’HoogheF, T., and Giudice, L. (2015). Update on Biomarkers for the Detection of Endometriosis. Biomed. Res. Int. 2015, 1130854–1130914. doi:10.1155/2015/130854
Hever, A., Roth, R. B., Hevezi, P., Marin, M. E., Acosta, J. A., Acosta, H., et al. (2007). Human Endometriosis Is Associated with Plasma Cells and Overexpression of B Lymphocyte Stimulator. Proc. Natl. Acad. Sci. 104, 12451–12456. doi:10.1073/pnas.0703451104
Hickey, M., Ballard, K., and Farquhar, C. (2014). Endometriosis. Bmj 348, g1752. doi:10.1136/bmj.g1752
Hornung, D., Bentzien, F., Wallwiener, D., Kiesel, L., and Taylor, R. N. (2001). Chemokine Bioactivity of RANTES in Endometriotic and Normal Endometrial Stromal Cells and Peritoneal Fluid. Mol. Hum. Reprod. 7, 163–168. doi:10.1093/molehr/7.2.163
Hu, W., Xie, Q., Xu, Y., Tang, X., and Zhao, H. (2020). Integrated Bioinformatics Analysis Reveals Function and Regulatory Network of miR-200b-3p in Endometriosis. Biomed. Res. Int. 2020, 1–7. doi:10.1155/2020/3962953
Hudson, Q. J., Ashjaei, K., Perricos, A., Kuessel, L., Husslein, H., Wenzl, R., et al. (2020). Endometriosis Patients Show an Increased M2 Response in the Peritoneal CD14+low/CD68+low Macrophage Subpopulation Coupled with an Increase in the T-Helper 2 and T-Regulatory Cells. Reprod. Sci. 27, 1920–1931. doi:10.1007/s43032-020-00211-9
Hull, M. L., Escareno, C. R., Godsland, J. M., Doig, J. R., Johnson, C. M., Phillips, S. C., et al. (2008). Endometrial-peritoneal Interactions during Endometriotic Lesion Establishment. Am. J. Pathol. 173, 700–715. doi:10.2353/ajpath.2008.071128
Irizarry, R. A., Hobbs, B., Collin, F., Beazer-Barclay, Y. D., Antonellis, K. J., Scherf, U., et al. (2003). Exploration, Normalization, and Summaries of High Density Oligonucleotide Array Probe Level Data. Biostatistics (Oxford, England) 4, 249–264. doi:10.1093/biostatistics/4.2.249
Jeung, I., Cheon, K., and Kim, M.-R. (2016). Decreased Cytotoxicity of Peripheral and Peritoneal Natural Killer Cell in Endometriosis. Biomed. Research International 2016, 1–6. doi:10.1155/2016/2916070
Kang, H., Lee, M., and Jang, S.-W. (2013). Celastrol Inhibits TGF-β1-Induced Epithelial-Mesenchymal Transition by Inhibiting Snail and Regulating E-Cadherin Expression. Biochem. Biophysical Res. Commun. 437, 550–556. doi:10.1016/j.bbrc.2013.06.113
Keilhauer, G., Romerdahl, C., Kupper, H., and Schlick, E. (1994). Modulation of Multidrug-Resistant Tumors by Dexverapamil. Oncol. Res. Treat. 17, 572–576. doi:10.1159/000218480
Kveiborg, M., Fröhlich, C., Albrechtsen, R., Tischler, V., Dietrich, N., Holck, P., et al. (2005). A Role for ADAM12 in Breast Tumor Progression and Stromal Cell Apoptosis. Cancer Res. 65, 4754–4761. doi:10.1158/0008-5472.Can-05-0262
Lamb, J., Crawford, E. D., Peck, D., Modell, J. W., Blat, I. C., Wrobel, M. J., et al. (2006). The Connectivity Map: Using Gene-Expression Signatures to Connect Small Molecules, Genes, and Disease. Science 313, 1929–1935. doi:10.1126/science.1132939
Langfelder, P., and Horvath, S. (2008). WGCNA: an R Package for Weighted Correlation Network Analysis. BMC Bioinformatics 9, 559. doi:10.1186/1471-2105-9-559
Lee, J.-H., Won, Y.-S., Park, K.-H., Lee, M.-K., Tachibana, H., Yamada, K., et al. (2012). Celastrol Inhibits Growth and Induces Apoptotic Cell Death in Melanoma Cells via the Activation ROS-dependent Mitochondrial Pathway and the Suppression of PI3K/AKT Signaling. Apoptosis 17, 1275–1286. doi:10.1007/s10495-012-0767-5
Lin, J., Yu, M., Xu, X., Wang, Y., Xing, H., An, J., et al. (2020). Identification of Biomarkers Related to CD8+ T Cell Infiltration with Gene Co-expression Network in clear Cell Renal Cell Carcinoma. Aging 12, 3694–3712. doi:10.18632/aging.102841
Madanes, D., Bilotas, M. A., Bastón, J. I., Singla, J. J., Meresman, G. F., Barañao, R. I., et al. (2020). PI3K/AKT Pathway Is Altered in the Endometriosis Patient's Endometrium and Presents Differences According to Severity Stage. Gynecol. Endocrinol. 36, 436–440. doi:10.1080/09513590.2019.1680627
Newman, A. M., Liu, C. L., Green, M. R., Gentles, A. J., Feng, W., Xu, Y., et al. (2015). Robust Enumeration of Cell Subsets from Tissue Expression Profiles. Nat. Methods 12, 453–457. doi:10.1038/nmeth.3337
Nyren-Erickson, E. K., Jones, J. M., Srivastava, D. K., and Mallik, S. (2013). A Disintegrin and Metalloproteinase-12 (ADAM12): Function, Roles in Disease Progression, and Clinical Implications. Biochim. Biophys. Acta (BBA) - Gen. Subjects 1830, 4445–4455. doi:10.1016/j.bbagen.2013.05.011
Painter, J. N., Nyholt, D. R., Morris, A., Zhao, Z. Z., Henders, A. K., Lambert, A., et al. (2011). High-density fine-mapping of a Chromosome 10q26 Linkage Peak Suggests Association between Endometriosis and Variants Close to CYP2C19. Fertil. Sterility 95, 2236–2240. doi:10.1016/j.fertnstert.2011.03.062
Patel, B. G., Lenk, E. E., Lebovic, D. I., Shu, Y., Yu, J., and Taylor, R. N. (2018). Pathogenesis of Endometriosis: Interaction between Endocrine and Inflammatory Pathways. Best Pract. Res. Clin. Obstet. Gynaecol. 50, 50–60. doi:10.1016/j.bpobgyn.2018.01.006
Polak, A., Bialopiotrowicz, E., Krzymieniewska, B., Wozniak, J., Stojak, M., Cybulska, M., et al. (2020). SYK Inhibition Targets Acute Myeloid Leukemia Stem Cells by Blocking Their Oxidative Metabolism. Cell Death Dis 11, 956. doi:10.1038/s41419-020-03156-8
Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). Limma powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies. Nucleic Acids Res. 43, e47. doi:10.1093/nar/gkv007
Rocks, N., Estrella, C., Paulissen, G., Quesada‐Calvo, F., Gilles, C., Guéders, M. M., et al. (2008). The Metalloproteinase ADAM‐12 Regulates Bronchial Epithelial Cell Proliferation and Apoptosis. Cel Prolif. 41, 988–1001. doi:10.1111/j.1365-2184.2008.00557.x
Rolla, E. (2019). Endometriosis: Advances and Controversies in Classification, Pathogenesis, Diagnosis, and Treatment. F1000Res 8, 529. doi:10.12688/f1000research.14817.1
Roy, R., and Moses, M. A. (2012). ADAM12 Induces Estrogen-independence in Breast Cancer Cells. Breast Cancer Res. Treat. 131, 731–741. doi:10.1007/s10549-011-1431-4
Sánchez-Barrena, M. J., Vallis, Y., Clatworthy, M. R., Doherty, G. J., Veprintsev, D. B., Evans, P. R., et al. (2012). Bin2 Is a Membrane Sculpting N-BAR Protein that Influences Leucocyte Podosomes, Motility and Phagocytosis. PLoS One 7, e52401. doi:10.1371/journal.pone.0052401
Shigetomi, H., Higashiura, Y., Kajihara, H., and Kobayashi, H. (2012). A Potential Link of Oxidative Stress and Cell Cycle Regulation for Development of Endometriosis. Gynecol. Endocrinol. 28, 897–902. doi:10.3109/09513590.2012.683071
Szklarczyk, D., Gable, A. L., Lyon, D., Junge, A., Wyder, S., Huerta-Cepas, J., et al. (2019). STRING V11: Protein-Protein Association Networks with Increased Coverage, Supporting Functional Discovery in Genome-wide Experimental Datasets. Nucleic Acids Res. 47, D607–d613. doi:10.1093/nar/gky1131
Tang, S., Lian, X., Jiang, J., Cheng, H., Guo, J., Huang, C., et al. (2020). Tumor Suppressive Maspin-Sensitized Prostate Cancer to Drug Treatment Through Negative Regulating Androgen Receptor Expression. Front. Cel Dev. Biol. 8, 573820. doi:10.3389/fcell.2020.573820
Trott, A., West, J. D., Klaić, L., Westerheide, S. D., Silverman, R. B., Morimoto, R. I., et al. (2008). Activation of Heat Shock and Antioxidant Responses by the Natural Product Celastrol: Transcriptional Signatures of a Thiol-Targeted Molecule. MBoC 19, 1104–1112. doi:10.1091/mbc.e07-10-1004
Tsai, H.-W., Huang, M.-T., Wang, P.-H., Huang, B.-S., Chen, Y.-J., and Hsieh, S.-L. (2018). Decoy Receptor 3 Promotes Cell Adhesion and Enhances Endometriosis Development. J. Pathol. 244, 189–202. doi:10.1002/path.5000
Vallvé-Juanico, J., Houshdaran, S., and Giudice, L. C. (2019a). The Endometrial Immune Environment of Women with Endometriosis. Hum. Reprod. Update 25, 565–592. doi:10.1093/humupd/dmz018
Vallvé-Juanico, J., Santamaria, X., Vo, K. C., Houshdaran, S., and Giudice, L. C. (2019b). Macrophages Display Proinflammatory Phenotypes in the Eutopic Endometrium of Women with Endometriosis with Relevance to an Infectious Etiology of the Disease. Fertil. Sterility 112, 1118–1128. doi:10.1016/j.fertnstert.2019.08.060
Woo, J.-H., Yang, Y.-I., Ahn, J.-H., Choi, Y. S., and Choi, J.-H. (2017). Interleukin 6 Secretion from Alternatively Activated Macrophages Promotes the Migration of Endometriotic Epithelial Cells†. Biol. Reprod. 97, 660–670. doi:10.1093/biolre/iox118
Yang, X., Zhu, S., Li, L., Zhang, L., Xian, S., Wang, Y., et al. (2018). Identification of Differentially Expressed Genes and Signaling Pathways in Ovarian Cancer by Integrated Bioinformatics Analysis. Ott 11, 1457–1474. doi:10.2147/ott.S152238
Zannoni, L., Forno, S. D., Paradisi, R., and Seracchioli, R. (2016). Endometriosis in Adolescence: Practical Rules for an Earlier Diagnosis. Pediatr. Ann. 45, e332–335. doi:10.3928/19382359-20160727-03
Zhang, S., Zhang, E., Long, J., Hu, Z., Peng, J., Liu, L., et al. (2019). Immune Infiltration in Renal Cell Carcinoma. Cancer Sci. 110, 1564–1572. doi:10.1111/cas.13996
Zhao, S. G., Lehrer, J., Chang, S. L., Das, R., Erho, N., Liu, Y., et al. (2019). The Immune Landscape of Prostate Cancer and Nomination of PD-L2 as a Potential Therapeutic Target. J. Natl. Cancer Inst. 111, 301–310. doi:10.1093/jnci/djy141
Keywords: endometriosis, CIBERSORT, weighted gene co-expression network analysis, M2 macrophages, diagnostic biomarker genes, therapeutic agents
Citation: Cui Z, Bhandari R, Lei Q, Lu M, Zhang L, Zhang M, Sun F, Feng L and Zhao S (2021) Identification and Exploration of Novel Macrophage M2-Related Biomarkers and Potential Therapeutic Agents in Endometriosis. Front. Mol. Biosci. 8:656145. doi: 10.3389/fmolb.2021.656145
Received: 20 January 2021; Accepted: 09 June 2021;
Published: 06 July 2021.
Edited by:
Mario Hiroyuki Hirata, University of São Paulo, BrazilReviewed by:
Philippa Saunders, University of Edinburgh, United KingdomChaoyan Yue, Fudan University, China
Copyright © 2021 Cui, Bhandari, Lei, Lu, Zhang, Zhang, Sun, Feng and Zhao. 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: Lijin Feng, ZmVuZ2xqaW5fMjA0OUB0b25namkuZWR1LmNu; Shasha Zhao, emhhb3NoYTIwMTJAMTI2LmNvbQ==
†These authors have contributed equally to this work