Integrative Analysis of m6A Regulator-Mediated RNA Methylation Modification Patterns and Immune Characteristics in Lupus Nephritis

Introduction

Lupus nephritis (LN) is the most common and most serious manifestation of systemic lupus erythematosus (SLE). It is also a major cause of morbidity and mortality in patients with SLE (Anders et al., 2020). Current treatments for LN are often ineffective and have strong adverse effects. In the last 50 years, only one drug has been developed for the treatment of SLE and LN, and other well-designed clinical trials have been unsuccessful (Thanou and Merrill, 2014; Narain and Furie, 2016). It is widely accepted that LN is caused by autoimmune and inflammatory responses owing to the loss of tolerance to endogenous nuclear material, which activates complement, pro-inflammatory pathways, and resident renal cells (Anders et al., 2020). Previously, the immune response to LN was mainly determined from the analysis of blood samples, which does not effectively reflect the immune state of the kidneys. Therefore, further investigation on the immune characteristics, including immune cell infiltration in LN kidney tissue, may be key in revealing its pathological mechanism and providing insight for the development of new immunotherapies for LN.

Genetic susceptibility can partially explain immune dysregulation in LN. Single-egg twins with the same gene only show a disease consistency of around 20–40%, suggesting that in addition to genetic susceptibility, epigenetics influenced by environmental factors also play an important role in SLE (Javierre et al., 2010) RNA methylation has been widely studied in epigenetic research. N6-methyladenosine (m⁶A) methylation is the most common RNA post-transcriptional modification that regulates gene expression outside of DNA sequences in eukaryotes and plays a key role in diseases progression (Meyer and Jaffrey, 2014; Huang et al., 2018). It is a reversible process mediated by an expanding list of m⁶A binding proteins (“readers”), adenosine methyltransferases (“writers”), and potential m⁶A demethylating enzymes (“erasers”) (Zaccara et al., 2019).

Current studies have demonstrated that m⁶A methylation is involved in immune regulation. For example, Han et al. (2019) discovered that the m⁶A binding protein YTHDF1 prolongs neoantigen-specific immunity through m⁶A methylation modification of mRNA. YTHDF1 is also involved in antigen cross-presentation and cross-priming of CD8⁺ T cells. Li et al. (2017) demonstrated that the m⁶A “writer” protein METTL3 regulates the homeostasis and differentiation of mouse T cells. However, no study has attempted to explore how m⁶A modification plays a role in LN, and the association between m⁶A modification and immune characteristics remains to be elucidated. The aim of this study was to clarify the role of m⁶A RNA methylation modification in LN and explore how m⁶A affects the immune status of LN.

Materials and Methods

Collection and Preprocessing of Data

The research strategy is presented in Figure 1. We collected gene expression data of patients with LN and healthy living donors from the Gene Expression Omnibus (GEO) database.¹ Four datasets (GSE32591 (Berthier et al., 2012), GSE69438 (Ju et al., 2015), GSE127797 (Almaani et al., 2019), and GSE112943) were selected in our study. GSE32591 contained 93 samples, which included 47 tubulointerstitium samples (32 LN samples, 15 living healthy donor samples), whereas 46 glomeruli samples (32 LN samples, 14 living healthy donor samples) were obtained from GPL14663 (Affymetrix Genechip HG-U133A). The platform for GSE69438 was GPL11670 (Affymetrix Human Genome U133 Plus 2.0 Array). It contained 42 tubulointerstitium samples, including 16 LN samples. The platform for GSE127797 was GPL24299 (Affymetrix Human Transcriptome Array 2.0), which contained 47 LN tubulointerstitium samples and 41 LN glomeruli samples. GSE127797 was the only dataset that included the pathological classifications of patients with LN. GSE112943 contained 21 whole kidney tissue samples, including 14 LN and 7 control, and sequencing was performed on GPL10558 (Illumina HumanHT-12 V4.0 expression beadchip). We then divided the samples into the glomeruli, tubulointerstitium, and whole kidney tissue for subsequent analysis. In total, 87 glomeruli samples (73 LN, 14 living healthy donors), 110 tubulointerstitium samples (95 LN, 15 living healthy donors), and 21 whole kidney tissues (14 LN, 7 controls) were included in our study to evaluate the expression levels of m⁶A regulators. All probes were converted into gene symbols, and median gene expression was used to represent the average expression level when multiple probes corresponded to the same gene symbol. We normalized the expression data from different datasets using the robust multi-array average, merged them together, and used the sva library for ComBat Batch correction to remove batch effects (Leek et al., 2012).

FIGURE 1

Figure 1. Study flow diagram. GEO, Gene Expression Omnibus; LN, lupus nephritis; HLA, human leukocyte antigen; LASSO, least absolute shrinkage and selection operator; GSVA, gene set variation analysis; ssGSEA, single sample gene set enrichment analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes; GO-BP, Gene Ontology Biological Processes; WGCNA, weighted gene co-expression network analysis; GO, Gene Ontology; PPI, protein-protein interaction.

m⁶A RNA Methylation Regulator Detection

The m⁶A RNA methylation regulator list we used was based on previous publications (Zhao et al., 2017; He et al., 2019; Zaccara et al., 2019). Then, the R package “limma” was applied to determine expression differences of m⁶A regulators between LN samples and healthy samples (including the glomeruli, tubulointerstitium, and kidney tissues) (Ritchie et al., 2015).

Development and Validation of m⁶A Regulator Gene Signature for LN

Univariate logistic regression was used to preliminarily screen variables in the identified m⁶A regulator list, and LASSO regression was used to select the best predictive features while fitting a generalized linear model and avoiding overfitting (Friedman et al., 2010). m⁶A regulators with non-zero LASSO regression coefficients were included in the multivariate logistic regression (MLR) analysis. The p-value in the MLR was based on the Wald test, and statistical significance was set at p < 0.05. Forest plots were drawn using the R package “ggplot2” to visually describe the results of the logistic regression. The receiver operating characteristic curve (ROC) and the average optimism of the area under the curve (AUC) quantified the predicted probabilities of the model. The risk score for each sample was calculated as follows:

$$Riskscore=\sum _{i=1}^{n}Coe{f}_i*x_i$$

where Coef_i indicates the coefficients of MLR and x_i is the gene expression value of each m⁶A regulator.

Correlation Between m⁶A Regulators and Immune Characteristics

The CIBERSORTx with 1,000 permutations was used to evaluate the abundance of infiltrating immunocytes.² The inclusion criterion was as follows: CIBERSORT, p < 0.05. We conducted single-sample gene-set enrichment analysis (ssGSEA) to assess immune response activity. We downloaded these gene sets from the ImmPort database (Bhattacharya et al., 2014).³ Lastly, Spearman correlation analysis was used to determine the correlation between m⁶A regulators and immune characteristics.

Unsupervised Cluster Analysis of m⁶A Modification Patterns in LN

Based on 18 identified m⁶A regulators, unsupervised cluster analysis was performed to determine distinct m⁶A subtypes using the R package “ConsensusClusterPlus,” and the consensus clustering algorithm ran 1,000 times to guarantee the robustness of clustering (Wilkerson and Hayes, 2010). The Kruskal test was used to compare the differences in m⁶A regulator expression and immune characteristics between subtypes. Principal component analysis was performed with the R package “PCA.”

Pathway Enrichment Analysis of the Two m⁶A Patterns

We downloaded the gene sets “h.all.v7.4.symbols” and “c2.cp.kegg.v7.4.symbols” from the MSigDB database. The gene set variation analysis (GSVA) algorithm was used to calculate the pathway activation score, which was conducted using the R package “GSVA” (Hänzelmann et al., 2013). The R package “limma” was used to compare the differences in pathway activation score between two subtypes, and a p-value < 0.01 was the cut-off criterion (Ritchie et al., 2015).

Identification of m⁶A Modification Pattern Markers

m⁶A modification subtypes-related differentially expressed genes (DEGs) between two distinct m⁶A subtypes (p < 0.0001) were defined as m⁶A related genes. m⁶A related genes were enriched in biological processes (BP), cellular component (CC), and molecular function (MF) terms in Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and were visualized with a bubble plot. We performed enrichment analysis with the cut-off criterion of the Q-value at < 0.05, and it was conducted using the “clusterProfiler” package (Yu et al., 2012).

Weighted gene co-expression network analysis (WGCNA) was conducted to identify m⁶A subtype-related genes and gene modules that characterize the pathways or functions of subtypes based on gene profiles using the WGCNA R package (Langfelder and Horvath, 2008). Correlations between different modules and subgroups were analyzed using Pearson’s correlation.

We used the STRING database⁴ to construct a protein-protein interaction (PPI) network for genes from the key module of the WGCNA. Visualization was performed using Cytoscape (Szklarczyk et al., 2019).⁵

Clinical Correlation With m⁶A Pattern Markers

The Nephroseq V5 tool⁶ was used to determine the correlation between m⁶A markers and renal function. We downloaded the expression data of markers and used “ggplot2” to replot the scatter plots.

Results

Landscape of m⁶A RNA Methylation Regulators in LN

Currently, there are 23 m⁶A RNA methylation regulators that have been widely studied, including 8 writers, 13 readers, and 2 erasers. Figure 2A shows the m⁶A regulators with functions and crosstalk between regulators and the immune microenvironment. The regulatory interactions of these 23 m⁶A regulators are shown in Figure 2B. First, the m⁶A regulator gene expression values of glomeruli (GSE32591 and GSE127797), tubulointerstitium (GSE32591, GSE69438, and GSE127797), and whole kidney tissues (GSE112943) of LN and healthy samples were evaluated. The expression of m⁶A regulators was the most considerably different in the glomeruli between LN and healthy samples. In total, 18 m⁶A RNA methylation regulators were identified in the glomeruli (Figure 2C). Significant expression differences in the 13 regulators (p < 0.05) were observed between LN and healthy samples, including WTAP, RBM15B, LRPPRC, and FTO (p < 0.001). Differences in the expression of m⁶A regulators between LN and healthy samples in the tubulointerstitium were not significant. As shown in Figure 2D, only six expressions altered m⁶A regulators in 17 identified m⁶A regulators. Significant expression differences in the 13 m⁶A regulators (p < 0.05) were observed among 21 m⁶A identified regulators in whole kidney tissue (Figure 2E). Taken together, the most significant differences in the expression of m⁶A regulators between LN and healthy samples were observed in the glomeruli. Thus, we selected the glomeruli samples for further detailed analysis.

FIGURE 2

Figure 2. Landscape of m⁶A RNA methylation regulators in LN. (A) m⁶A RNA methylation modification regulated by m⁶A “writer,” “reader,” and “eraser,” which is involved in the immune microenvironment of LN. (B) Protein-protein interaction (PPI) network composed of 23 m⁶A regulators. (C) Violin plot demonstrating the expression level of 18 m⁶A regulators in glomeruli between living donor and LN. (D) Violin plot demonstrating the expression level of 17 m⁶A regulators in tubulointerstitium between living donor and LN. (E) Violin plot demonstrating the expression level of 21 m⁶A regulators in kidney whole tissue between living donor and of LN. (F) Volcano plot showing a summary of the expression differences of 18 m⁶A regulators between the healthy and LN patients’ glomerular samples. (G) Correlations between 18 m⁶A regulators in LN glomeruli samples. The two respective scatterplots show the two pairs of m⁶A regulators with the highest correlation, HNRNPA2B1 and RBM15B with the most negative correlation, and YTHDC1 and FMR1 with the most positive correlation.

Interestingly, the expression of all 13 altered m⁶A regulators was downregulated in LN compared with healthy samples in the glomeruli (Figure 2F). The decrease in fold change of RBM15 was the largest among these genes, whereas the decrease in eraser protein FTO levels was the most statistically significant. Note that there was no significant difference in the expression of the well-studied writer METTL13 between LN and healthy samples. In the correlation analysis for 18 m⁶A regulators, we observed that close relationships among m⁶A regulators which means they function together (Figure 2G). It contributes to explore the specific mechanism of aberrant m⁶A modification in LN. Figure 2G also shows the two pairs of m⁶A regulators with the highest positive/negative correlation. The reader HNRNPA2B1 and writer RBM15B were the most negatively correlated, whereas the readers YTHDC1 in the nuclei and FMR1 in the cytoplasm were the most positively correlated.

m⁶A Regulators Have the Potential to Distinguish Between Healthy and LN Samples

To better understand the contribution of m⁶A regulators to the progression of LN, we established an m⁶A regulator gene signature. First, 16 regulators were selected from 18 identified m⁶A regulators by univariate logistic regression analysis (Figure 3A). LASSO regression was performed to further screen the m⁶A candidates, in which 13 regulators with non-zero coefficients were included in the multivariate logistic regression (Figures 3B,C). Finally, multivariate logistic regression analysis demonstrated that METTL3, WTAP, YTHDC2, YTHDF1, FMR1, and FTO were independently associated with LN (Figure 3D). It should be noted that the well-studied writer protein METTL3 is an independent risk factor for LN, although its expression does not differ between LN and in healthy samples. How METTL3 plays a role in LN might be an interesting topic for further exploration. For the gene signature model, the AUC for the derivation sets was 0.949, indicating that this model performed well in classifying healthy and LN samples (Figure 3E). In the independent external validation set (GSE112943), the AUC was 0.962, which suggests its ability to classify the samples (Figure 3F). In addition, Figure 3G shows that there was a significant difference in m⁶A risk scores between LN and healthy samples. The risk scores for LN were noticeably higher than those for healthy samples (Figure 3G). The distribution of risk scores and gene profiles based on the six selected m⁶A regulators is shown in Figure 3H. We observed that the expression of WTAP, YTHDC2, and FTO decreased in the high-risk group.

FIGURE 3

Figure 3. m⁶A regulators have the potential to distinguish between healthy and LN individuals. (A) Univariate logistic regression revealed 16 LN-related m⁶A regulators (P < 0.05). (B,C) Feature selection by LASSO regression model. (B) By verifying the optimal parameter (lambda) in the LASSO model, the partial likelihood deviance (binomial deviance) curve was plotted vs. log (lambda). Dotted vertical lines were drawn based on 1 SE of the minimum criteria (the 1-SE criteria). (C) Thirteen features with non-zero coefficients were selected by optimal lambda. A coefficient profile plot was produced against the log (lambda) sequence in (B). (D) Multivariate logistic analysis distinguished six independent risk factors and risk scores for LN were calculated using the LASSO Logistic regression algorithm. (E,F) The predictive value of the m⁶A regulator gene signature in the derivation (E) and validation (F) sets by calculating the pooled AUC. 0.9 < AUC ≤ 1 indicates that the gene signature has high accuracy. (G) Distribution of risk scores in healthy and LN samples. (H) Risk score distribution based on the 6 m⁶A RNA modification regulator signature and gene expression profiles between our study groups. Patients were divided into high-risk and low-risk groups by the black dotted line, which indicates the median cut-off value.

m⁶A Regulators Are Related to Immune Microenvironment in LN

To further elucidate the relationship between m⁶A regulators and immune characteristics, we analyzed the correlations between them. Immune characteristics include immune cell infiltration, immune response activity, and Human Leukocyte Antigen (HLA) genes. The abundance of 22 infiltrating immunocytes in the glomeruli of LN was evaluated using the CIBERSORTx algorithm (Supplementary Figure 1A). Eosinophils were excluded from the correlation analysis because of the lack of expression in all samples. Several infiltrating immunocytes were correlated with m⁶A regulators but were mostly weakly correlated (Figure 4A). Among all immunocytes, activated NK cells were closely correlated with m⁶A regulators, and these were most positively correlated with HNRNPA2B1 and most negatively correlated with RBM15B. This indicates that NK-activated cell infiltration in LN is regulated by HNRNPA2B1 and RBM15B. Supplementary Figure 1B shows the expression differences of each immune response between LN and healthy samples. Correlation analysis demonstrated that most immune reaction pathways were closely related to m⁶A regulators (Figure 4B). Cytokinesis, inflammation pathway, interferon receptor activity, interferon-mediated signaling pathway, and TGF-β pathway were correlated with most of the 18 m⁶A regulators, indicating that immune dysregulation in LN is affected by m⁶A RNA methylation. Both the most positive and negative correlations between regulators and immune reactions were related to the reader protein YTHDC1, indicating that YTHDC1 exerts important functions in the cytokinesis and inflammation pathways in LN. Similarly, HLA genes were closely correlated with m⁶A regulators (Figure 4C). HLA-DRA was most positively correlated with HNRNPA2B1, with a correlation coefficient of 0.67. HLA-F was most negatively correlated with LRPPRC, with a correlation coefficient of –0.61. These indicate that HLA gene expression in LN was affected by m⁶A regulators. The differences in HLA gene expression between LN and healthy samples were shown in Supplementary Figure 1C.

FIGURE 4

Figure 4. Correlation between m⁶A regulator expression and immune characteristics in LN. (A) Heatmap of the correlations between 18 m⁶A regulators and 21 immunocytes (eosinophils with no expression were removed in all samples). The two respective scatterplots show the m⁶A regulator and immunocyte with the highest positive or negative correlation. (B) Heatmap of the correlations between 18 m⁶A regulators and immune response gene sets. The two respective scatterplots show m⁶A regulators and immune response gene sets with the highest positive or negative correlation. (C) Heatmap of the correlations between 18 m⁶A regulators and 18 HLA genes. The two respective scatterplots show m⁶A regulators and HLA genes with the highest positive or negative correlation.

Identification of m⁶A RNA Methylation Subtypes Based on 18 m⁶A Regulators in LN and Clinical Correlation of 2 Subtypes

To identify m⁶A RNA methylation modification patterns of LN, we conducted unsupervised clustering based on the expression similarity of m⁶A regulators in LN and k = 2 seemed to be an adequate selection resulted in 2 distinct subtypes (Figures 5A–C). Two m⁶A subtypes had significantly different populations in PCA (Figure 5D). To investigate the relationship between clinical characteristics and m⁶A subtypes, we used data from GSE127797, including the pathological stages of LN, to create a correlation heatmap. There were 38 LN samples, consisting of subtype 1 with 12 samples and subtype 2 with 26 samples. Most patients with mixed proliferative and membranous LN (class III+V and IV+V) belong to subtype 1, whereas most with pure proliferative LN (class III and IV) or pure membranous LN (class V) belong to subtype 2. Significant differences in m⁶A regulator gene profiles were observed between the two subtypes (Figures 5E,F). FTO, HNRNPC, HNRNPA2B1, YTHDC2, ZC3H13, YTHDC1, YTHDF3, FMR1, and LRPPRC were highly expressed in m⁶A subtype 1, whereas the other regulators were highly expressed in m⁶A subtype 2.

FIGURE 5

Figure 5. Identification of two distinct m⁶A modification subtypes in LN and clinical correlation of two subtypes. (A) Consensus clustering of cumulative distribution function (CDF) for k = 2–9. (B) Elbow plot shows relative change in area under CDF curve. (C) Consensus clustering matrix for k = 2. (D) Principal component analysis (PCA) of two m⁶A subtypes in LN. (E) Heatmap of the clinical features of two clusters comparing the stages of LN and gene profiles between m⁶A subtypes in GSE127797. (F) The two m⁶A subtypes exhibit distinct expression statuses of 18 m⁶A RNA methylation regulators.

Immune Characteristics and Biological Functional Characteristics of Two Distinct m⁶A Subtypes

To further determine the characteristics of the two m⁶A subtypes, we compared the abundance of infiltrating immune cells, activity of immune responses and HLA gene expression value between the two distinct m⁶A subtypes. More infiltrating activated NK cells (p = 0.001), memory resting CD4 T cells (p = 0.009), and activated dendritic cells (p = 0.03) were observed in subtype 1, whereas more plasma cells (p = 0.009), naïve CD4 T cells (p = 0.03), and macrophages M0 (p = 0.003) were observed in subtype 2 (Figure 6A). For immune reactions, m⁶A subtype 1 had a stronger immune response than subtype 2. There were 11 immune reactions, including MHC-I-mediated antigen processing, cytokinesis, and interferon-mediated signaling pathways, which were more active in subtype 1, whereas BCR, CTL, inflammation, and IL-12 pathways were more active in subtype 2 (Figure 6B). Different expression levels of HLA genes between the two m⁶A subtypes were also observed. For example, subtype 1 had a higher expression of HLA-C, HLA-DPA1, and HLA-DRA, whereas subtype 2 had a higher expression of HLA-DOB and HLA-DQB2 (Figure 6C). Taken together, m⁶A modification patterns were shaped with different immune characteristics, suggesting that m⁶A RNA methylation regulators might play an important role in immune microenvironment regulation in LN. To investigate the biological functional pathways that m⁶A may affect, we conducted GSVA to assess the enrichment of biological pathways. Figure 6D shows the enrichment difference of the HALLMARKS pathways between the two subtypes, indicating that protein secretion and UV-response pathways are more enriched in subtype 1, whereas myogenesis and KRAS signaling pathways are more enriched in subtype 2. Some KEGG pathways, including regulation of autophagy, TGF-β signaling pathway, and antigen processing and presentation were more enriched in subtype 1, whereas other pathways such as cytokine-cytokine receptor interaction, intestinal immune network for IgA production, and the JAK-STAT signaling pathway were more enriched in subtype 2 (Figure 6E). It should be noted that the highly enriched pathways in both subtypes both included immune-related pathways.

FIGURE 6

Figure 6. Differences in immune characteristics between m⁶A subtypes and functional enrichment analysis in two m⁶A subtypes. (A) Differences in abundance of 22 infiltrating immunocytes. (B) Differences in the activity of 22 immune response gene sets in two m⁶A subtypes. (C) Expression differences of 18 HLA genes in two m⁶A subtypes. (D) Differences in HALLMARKS pathway enrichment between m⁶A subtypes. (E) KEGG pathways with significant differences in enrichment between m⁶A subtypes.

Identification of m⁶A Methylation Modification Markers and Clinical Correlation of Markers With Renal Function

To gain further insight into which genes are involved in the biological processes affected by m⁶A regulators, we identified m⁶A-related genes, and enrichment analysis of these genes was performed. The top 10 pathways in the BP and GO were mainly RNA or protein modification pathways and immune-related pathways, such as neutrophil degranulation and MHC-I-mediated antigen processing (Figure 7A). This confirmed that m⁶A methylation modification was associated with the immune microenvironment in LN. The KEGG enrichment analysis revealed that m⁶A modification patterns-related genes were more enriched in protein processing in the endoplasmic reticulum and Salmonella infection pathways (Figure 7B). Based on m⁶A related genes, we performed WGCNA to identify module hub genes (Figures 7C–E). Three gene modules were established, including the nonsense gray module, based on their similar expression spectrum (Figure 7F). The MEturquoise module genes were most positively correlated with m⁶A subtype 1 (R² = 0.78) (Figure 7G), indicating that MEturquoise is a key module. Then, genes in MEturquoise were used to construct the PPI network (Figure 7H). If the module membership (MM) of genes in the turquoise module was > 0.8, and their gene significance (GS) was > 0.6, these genes were considered the hub genes of the turquoise module. Finally, we overlapped the central nodes in the PPI and hub genes of the turquoise module, and seven m⁶A RNA methylation modification markers (CDC5L, CDC40, HNRNPU, NUDT21, PAPOLA, POLR2B, and WBP4) were identified (Figure 7I). To further elucidate the roles of these m⁶A markers in LN, correlation analysis between markers and GFR was carried out in the Nephroseq database (Figures 8A–G). Among the seven markers, only CDC40 was positively correlated with GFR, thus, a higher expression of CDC40 indicates better renal function in patients with LN have and may play a protective role against LN. The other six markers, CDC5L, HNRNPU, NUDT21, PAPOLA, POLR2B, and WBP4, were all negatively correlated with GFR, indicating that these genes may aggravate kidney damage in patients with LN.

FIGURE 7

Figure 7. Pathway enrichment analysis of m⁶A regulator related genes (A,B) and identification of m⁶A methylation pattern markers in LN (C–I). (A) Enrichment analysis of GO biological process, cellular component, and molecular function. (B) Bubble plot of KEGG enrichment pathways. (C) Clustering dendrogram of two m⁶A modification subtypes in LN. (D) Scale-free fitting index analysis and mean connectivity of soft threshold power from 1 to 20. (E) Clustering dendrograms for m⁶A regulator-related genes. According to dynamic tree cutting, the genes were clustered into different modules through hierarchical clustering and merged when the correlation of the modules is > 0.8. Each color represents each module. (F) Correlation heatmap between module eigen genes and m⁶A subtypes. (G) Scatter plot of m⁶A subtype 1 in the turquoise module. In the turquoise module, GS and MM show a very significant correlation, indicating that the genes of the turquoise module are highly related to the m⁶A modification subtype. The dots in the red box indicate that the module membership of these genes is > 0.8, and their gene significance > 0.6, meaning that these dots are the hub genes of the turquoise module. (H) PPI analysis network of m⁶A methylation-related genes from the turquoise module 10, the central nodes in PPI are marked in red, orange, and yellow. (I) Venn diagram of seven m⁶A modification markers. The central nodes of PPI (green set) were overlapped with the hub genes in the turquoise module (blue set) by weighted correlation network analysis.

FIGURE 8

Figure 8. Relationship between seven m⁶A methylation pattern markers and renal function (glomerular filtration rate).

Discussion

LN is an autoimmune disease characterized by symptoms of inflammation. Immune response dysregulation mediated by genetic and environmental factors leads to the occurrence and development of LN (Anders et al., 2020). Many studies have confirmed that m⁶A methylation modification exerts critical functions in the development of diseases, especially malignancies. However, little research has been conducted on m⁶A methylation in LN. Our study is the first to investigate the roles of m⁶A regulators in LN and reveal the association with m⁶A methylation modification and immune characteristics. Firstly, significant differences in the expression of most m⁶A regulators between healthy individuals and LN were observed in the glomeruli. This is mainly because LN is a form of glomerulonephritis. We also identified an m⁶A regulator gene signature that included METTL3, WTAP, YTHDC2, YTHDF1, and FMR1 after LASSO-logistic regression. LN and healthy samples were easily distinguished, which highlights that m⁶A methylation modification patterns differ between LN and healthy samples.

Then, we demonstrated a correlation between m⁶A regulators and immune characteristics. A series of immune reactions were considerably correlated with m⁶A regulators, especially MHC-I-mediated antigen processing, cytokinesis, inflammation pathway, and interferon-mediated signaling pathway. Most m⁶A regulators were found to be strongly correlated with the IFN signaling pathway. Recent studies have shown that type I interferon (IFN-I) is an important risk factor for the occurrence and progression of LN (Ding et al., 2021), indicating that m⁶A methylation modification plays a key role in the development of LN. Most HLA genes were closely correlated with m⁶A regulators. Studies have identified that HLA-DR3, HLA-DR4, HLA-DR11, and HLA-DR15 can promote or improve kidney damage in LN (Munroe and James, 2015; Iwamoto and Niewold, 2017). For immune cell infiltration, activated NK cells were most strongly correlated with m⁶A regulators. Activated NK cells were most positively correlated with HNRNPA2B1 and most negatively correlated with RBM15B. NK cells are an important link between the innate and adaptive immune systems. Postól et al. (2008) found that the onset of glomerulonephritis in NZBxNZW (F1) mice (SLE model) can be delayed by long-term depletion of NK cells, indicating that functional defects in NK cells may induce the development of LN (Spada et al., 2015; Segerberg et al., 2019).

However, the overall correlation between various immunocytes and m⁶A regulators was found to be generally weak. One possible reason for this is the limitations of previous technical tools. The samples collected for RNA sequencing contain very limited immune cells, which might not precisely reflect the abundance of infiltrating immunocytes (Stewart et al., 2020).

In our study, two distinct m⁶A RNA methylation modification subtypes were identified based on the m⁶A regulator gene expression profiles using unsupervised clustering. The differences between the 2 subtypes included the following aspects. As for immune characteristics, m⁶A subtype 1 of LN was characterized by increased immune response activation and a higher HLA gene expression profile. A higher abundance of infiltrating immune cells was observed in m⁶A subtype 2, including that of plasma cells, naïve CD4 T cells, M0 macrophages, and dendritic cells. Plasma cells play a key role in the development of SLE and LN (Crickx et al., 2021). A higher abundance of infiltrating T cells, memory resting CD4 T cells, and activated NK cells was found in subtype 1. The deposition of immunoglobulins produced by plasma cells in the glomeruli is the initial trigger for LN. Now, the focus of treatment for LN is targeted B-cell therapy. The two m⁶A subtypes of LN have the potential to be used to develop targeted immunotherapy.

Additionally, the pathological stages of LN in the two subtypes were also considerably different. Most patients with mixed proliferative and membranous LN (class III+V and IV+V) belong to subtype 1, whereas most with pure proliferative LN (class III and IV) or pure membranous LN (class V) belong to subtype 2. No current research has been conducted to illustrate the relationship between LN classification and immune status. Our results initially suggest that a greater activation of immune reaction pathways occurs in class III+V and IV+V, and a higher abundance of infiltrating plasma cells occurs in class III, IV, and V. In general, class III+V and IV+V have more complicated pathological changes than class III, IV, or V because their lesions are mixed proliferative lesions in class V (Parikh et al., 2020). LN classification is established according to differences in prognosis and is the gold standard for guiding treatment (Levey and Coresh, 2012; Mackay et al., 2019; Mageau et al., 2019). In our study, a strong correlation was observed between LN classification and m⁶A subtypes. As subtype1 is associated with mixed, complex, and more severe lesions compared with subtype2, uncovering the key differences between the subtypes will contribute to preventing the aggravation of LN. Molecular subtyping is a widely used strategy in malignancies, and targeted treatment plans can be formulated based on different molecular types to improve patient prognosis (Teo et al., 2019). The two m⁶A modification subtypes of LN have the potential to be considered as an alternative classification of LN. Furthermore, from a functional pathway perspective, genes of m⁶A subtype 1 are more enriched in the TGF-β signaling pathway, MTOR signaling pathway, and autophagy regulation.

Finally, seven m⁶A methylation modification markers were identified. CDC5L, HNRNPU, NUDT21, PAPOLA, POLR2B, and WBP4 were negatively correlated with GFR (an indicator of kidney function), whereas CDC40 was positively correlated with GFR. The protein encoded by CDC5L has been shown to be a positive regulator of the G2/M stage of the cell cycle. Zhou et al. (2020) found that CDC5L also regulates cell proliferation and metastasis in lung adenocarcinoma through promoter methylation. POLR2B encodes the second largest subunit of RNA polymerase II (Pol II), which is involved in RNA splicing and modification (Wang et al., 2021). CDC40, HNRNPU, WBP4, NUDT21, and PAPOLA are all involved in precursor mRNA splicing. S-adenosylmethionine (SAM) is a methyl donor for almost all cell methylation events. Scarborough et al. (2021) reported that NUDT21 regulates intracellular SAM levels. As shown, m⁶A RNA methylation modification plays a key role in mRNA splicing, suggesting that the m⁶A markers identified in this study are closely related to the m⁶A modification process. In addition, Fuyuno et al. (2016) proposed that the mutation of the WBP4 locus may lead to the occurrence of inflammatory bowel disease. The nuclear matrix protein HNRNPU is also considered as a nuclear virus dsRNA sensor for DNA and RNA viruses (Lin et al., 2010). m⁶A markers may be related to immune disorders and inflammatory responses, which again highlights that m⁶A-regulators can regulate immune characteristics. At present, research on LN mainly focuses on genetics and clinical advances, whereas epigenetic research is rare. There is also almost no research on m⁶A RNA methylation modification. We took the lead in identifying the role of m⁶A regulators in LN and exploring their relationship with immune characteristics. The various results in our study indicate that m⁶A methylation modification is a new direction for research regarding the pathogenesis of LN.

Our study has certain limitations. First, we were unable to obtain more clinical data for each patient, such as sex, age, treatment, and prognosis, for the longitudinal analysis. We could not perform a correlation analysis between m⁶A patterns, pathological stages, and other clinical characteristics of all samples. Second, we included as many samples as possible in the GEO database that met our requirements, but the sample size was still small (73 LN, 14 living healthy donors). Studies with larger sample sizes are required in the future. Additionally, the expression changes of some identified m⁶A regulators between LN and healthy samples were small and our findings were mainly obtained through bioinformatics analysis, which is needed to be further verified experimentally. However, the good predictive performance of our identified m⁶A regulator gene signature was verified using an external data set. The correlation between seven m⁶A markers identified from GEO data and kidney function was verified using data from the Nephroseq database. We have sufficient reasons to believe that m⁶A methylation plays an important role in the development of LN.

In summary, we comprehensively assessed the role of m⁶A methylation in the glomeruli of patients with LN, established an m⁶A regulators signature that can easily distinguish LN and healthy individuals, and identified two distinct m⁶A subtypes based on 18 m⁶A regulators. The two distinct m⁶A subtypes in LN were determined with significant differences in m⁶A regulators expression, immune microenvironment, biological functional pathways, and clinical characteristics. We uncovered an association between m⁶A subtypes and immune characteristics, which can be used to develop targeted immunotherapy. Moreover, seven m⁶A subtype markers were identified and all of them demonstrated a meaningful correlation with GFR, indicating that they are potential prognostic biomarkers.

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 in the article/Supplementary Material.

Author Contributions

HZ, ZL, and DL conceived the project. HZ analyzed the data, explained the results, drafted the original manuscript, and drew the figures. SP and JD discussed the draft manuscript. FL and GL revised the manuscript. All the authors approved the final version of the manuscript.

Funding

This work was supported by the National Science and Technology Major Project of China (Grant No. 2020ZX09201-009), the Innovation Scientists and Technicians Troop Construction Projects of Henan Province (No. 182101510002), and the Excellent Youth Fund Project of Henan Province (No. 202300410363).

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.

Acknowledgments

We would like to thank all the contributors of the transcriptome sequencing data for this study.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2021.724837/full#supplementary-material

Footnotes

1. ^ http://www.ncbi.nlm.nih.gov/geo/
2. ^ https://cibersort.stanford.edu/
3. ^ https://string-db.org/
4. ^ https://cytoscape.org/
5. ^ http://v5.nephroseq.org/

References

Almaani, S., Prokopec, S. D., Zhang, J., Yu, L., Avila-Casado, C., Wither, J., et al. (2019). Rethinking lupus nephritis classification on a molecular level. J. Clin. Med. 8:1524. doi: 10.3390/jcm8101524

CrossRef Full Text | Google Scholar

Anders, H.-J., Saxena, R., Zhao, M.-H., Parodis, I., Salmon, J. E., and Mohan, C. (2020). Lupus nephritis. Nat. Rev. Dis. Primers 6:7.

Google Scholar

Berthier, C. C., Bethunaickan, R., Gonzalez-Rivera, T., Nair, V., Ramanujam, M., Zhang, W., et al. (2012). Cross-species transcriptional network analysis defines shared inflammatory responses in murine and human lupus nephritis. J. Immunol. 189, 988–1001. doi: 10.4049/jimmunol.1103031

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhattacharya, S., Andorf, S., Gomes, L., Dunn, P., Schaefer, H., Pontius, J., et al. (2014). ImmPort: disseminating data to the public for the future of immunology. Immunol. Res. 58, 234–239. doi: 10.1007/s12026-014-8516-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Crickx, E., Tamirou, F., Huscenot, T., Costedoat-Chalumeau, N., Rabant, M., Karras, A., et al. (2021). Evolution of kidney antibody secreting cells molecular signature in lupus patients with active nephritis upon immunosuppressive therapy. Arthritis Rheumatol. (Hoboken, NJ) 73, 1461–1466. doi: 10.1002/art.41703

PubMed Abstract | CrossRef Full Text | Google Scholar

Ding, X., Ren, Y., and He, X. I. F. N.-I. (2021). Mediates lupus nephritis from the beginning to renal fibrosis. Front. Immunol. 12:676082.

Google Scholar

Friedman, J., Hastie, T., and Tibshirani, R. (2010). Regularization paths for generalized linear models via coordinate descent. J. Stat. Softw. 33, 1–22.

Google Scholar

Fuyuno, Y., Yamazaki, K., Takahashi, A., Esaki, M., Kawaguchi, T., Takazoe, M., et al. (2016). Genetic characteristics of inflammatory bowel disease in a Japanese population. J. Gastroenterol. 51, 672–681.

Google Scholar

Han, D., Liu, J., Chen, C., Dong, L., Liu, Y., Chang, R., et al. (2019). Anti-tumour immunity controlled through mRNA m 6 A methylation and YTHDF1 in dendritic cells. Nature 566, 270–274. doi: 10.1038/s41586-019-0916-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Hänzelmann, S., Castelo, R., and Guinney, J. (2013). GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 14:7. doi: 10.1186/1471-2105-14-7

PubMed Abstract | CrossRef Full Text | Google Scholar

He, L., Li, H., Wu, A., Peng, Y., Shu, G., and Yin, G. (2019). Functions of N6-methyladenosine and its role in cancer. Mol. Cancer 18:176.

Google Scholar

Huang, H., Weng, H., Sun, W., Qin, X., Shi, H., Wu, H., et al. (2018). Recognition of RNA N(6)-methyladenosine by IGF2BP proteins enhances mRNA stability and translation. Nat. Cell Biol. 20, 285–295. doi: 10.1038/s41556-018-0045-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Iwamoto, T., and Niewold, T. B. (2017). Genetics of human lupus nephritis. Clin. Immunol. 185, 32–39. doi: 10.1016/j.clim.2016.09.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Javierre, B. M., Fernandez, A. F., Richter, J., Al-Shahrour, F., Martin-Subero, J. I., Rodriguez-Ubreva, J., et al. (2010). Changes in the pattern of DNA methylation associate with twin discordance in systemic lupus erythematosus. Genome Res. 20, 170–179. doi: 10.1101/gr.100289.109

PubMed Abstract | CrossRef Full Text | Google Scholar

Ju, W., Nair, V., Smith, S., Zhu, L., Shedden, K., Song, P. X. K., et al. (2015). Tissue transcriptome-driven identification of epidermal growth factor as a chronic kidney disease biomarker. Sci. Transl. Med. 7:316ra193.

Google Scholar

Langfelder, P., and Horvath, S. (2008). WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 9:559.

Google Scholar

Leek, J. T., Johnson, W. E., Parker, H. S., Jaffe, A. E., and Storey, J. D. (2012). The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics 28, 882–883. doi: 10.1093/bioinformatics/bts034

PubMed Abstract | CrossRef Full Text | Google Scholar

Levey, A. S., and Coresh, J. (2012). Chronic kidney disease. Lancet 379, 165–180.

Google Scholar

Li, H. B., Tong, J., Zhu, S., Batista, P. J., Duffy, E. E., Zhao, J., et al. (2017). m(6)A mRNA methylation controls T cell homeostasis by targeting the IL-7/STAT5/SOCS pathways. Nature 548, 338–342. doi: 10.1038/nature23450

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, R. K., Hsieh, Y. S., Lin, P., Hsu, H. S., Chen, C. Y., Tang, Y. A., et al. (2010). The tobacco-specific carcinogen NNK induces DNA methyltransferase 1 accumulation and tumor suppressor gene hypermethylation in mice and lung cancer patients. J. Clin. Invest. 120, 521–532. doi: 10.1172/jci40706

PubMed Abstract | CrossRef Full Text | Google Scholar

Mackay, M., Dall’Era, M., Fishbein, J., Kalunian, K., Lesser, M., Sanchez-Guerrero, J., et al. (2019). Establishing surrogate kidney end points for lupus nephritis clinical trials: development and validation of a novel approach to predict future kidney outcomes. Arthritis Rheumatol. (Hoboken, NJ) 71, 411–419. doi: 10.1002/art.40724

PubMed Abstract | CrossRef Full Text | Google Scholar

Mageau, A., Timsit, J. F., Perrozziello, A., Ruckly, S., Dupuis, C., Bouadma, L., et al. (2019). The burden of chronic kidney disease in systemic lupus erythematosus: a nationwide epidemiologic study. Autoimmun. Rev. 18, 733–737. doi: 10.1016/j.autrev.2019.05.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Meyer, K. D., and Jaffrey, S. R. (2014). The dynamic epitranscriptome: N6-methyladenosine and gene expression control. Nat. Rev. Mol. Cell Biol. 15, 313–326. doi: 10.1038/nrm3785

PubMed Abstract | CrossRef Full Text | Google Scholar

Munroe, M. E., and James, J. A. (2015). Genetics of lupus nephritis: clinical implications. Semin. Nephrol. 35, 396–409. doi: 10.1016/j.semnephrol.2015.08.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Narain, S., and Furie, R. (2016). Update on clinical trials in systemic lupus erythematosus. Curr. Opin. Rheumatol. 28, 477–487. doi: 10.1097/bor.0000000000000311

PubMed Abstract | CrossRef Full Text | Google Scholar

Parikh, S. V., Almaani, S., Brodsky, S., and Rovin, B. H. (2020). Update on lupus nephritis: core curriculum 2020. Am. J Kidney Dis. 76, 265–281. doi: 10.1053/j.ajkd.2019.10.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Postól, E., Meyer, A., Cardillo, F., de Alencar, R., Pessina, D., Nihei, J., et al. (2008). Long-term administration of IgG2a anti-NK1.1 monoclonal antibody ameliorates lupus-like disease in NZB/W mice in spite of an early worsening induced by an IgG2a-dependent BAFF/BLyS production. Immunology 125, 184–196. doi: 10.1111/j.1365-2567.2008.02835.x

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Scarborough, A. M., Flaherty, J. N., Hunter, O. V., Liu, K., Kumar, A., Xing, C., et al. (2021). SAM homeostasis is regulated by CFI(m)-mediated splicing of MAT2A. Elife 10:e64930.

Google Scholar

Segerberg, F., Lundtoft, C., Reid, S., Hjorton, K., Leonard, D., Nordmark, G., et al. (2019). Autoantibodies to killer cell immunoglobulin-like receptors in patients with systemic lupus erythematosus induce natural killer cell hyporesponsiveness. Front. Immunol. 10:2164.

Google Scholar

Spada, R., Rojas, J. M., and Barber, D. F. (2015). Recent findings on the role of natural killer cells in the pathogenesis of systemic lupus erythematosus. J. Leukoc. Biol. 98, 479–487. doi: 10.1189/jlb.4ru0315-081rr

PubMed Abstract | CrossRef Full Text | Google Scholar

Stewart, B. J., Ferdinand, J. R., and Clatworthy, M. R. (2020). Using single-cell technologies to map the human immune system - implications for nephrology. Nat. Rev. Nephrol. 16, 112–128. doi: 10.1038/s41581-019-0227-3

PubMed Abstract | CrossRef Full Text | Google Scholar

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.

Google Scholar

Teo, W. Y., Sekar, K., Seshachalam, P., Shen, J., Chow, W. Y., Lau, C. C., et al. (2019). Relevance of a TCGA-derived glioblastoma subtype gene-classifier among patient populations. Sci. Rep. 9:7442.

Google Scholar

Thanou, A., and Merrill, J. T. (2014). Treatment of systemic lupus erythematosus: new therapeutic avenues and blind alleys. Nat. Rev. Rheumatol. 10, 23–34. doi: 10.1038/nrrheum.2013.145

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J., Ahimaz, P. R., Hashemifar, S., Khlevner, J., Picoraro, J. A., Middlesworth, W., et al. (2021). Novel candidate genes in esophageal atresia/tracheoesophageal fistula identified by exome sequencing. Eur. J. Hum. Genet. 29, 122–130. doi: 10.1038/s41431-020-0680-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilkerson, M. D., and Hayes, D. N. (2010). ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics 26, 1572–1573. doi: 10.1093/bioinformatics/btq170

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, G., Wang, L. G., Han, Y., and He, Q. Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 16, 284–287. doi: 10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

Zaccara, S., Ries, R. J., and Jaffrey, S. R. (2019). Reading, writing and erasing mRNA methylation. Nat. Rev. Mol. Cell Biol. 20, 608–624. doi: 10.1038/s41580-019-0168-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, B. S., Roundtree, I. A., and He, C. (2017). Post-transcriptional gene regulation by mRNA modifications. Nat. Rev. Mol. Cell Biol. 18, 31–42. doi: 10.1038/nrm.2016.132

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, R. H., Zhang, J. T., Chen, C., Xu, Z. H., Lv, X. B., Ye, L., et al. (2020). Identification of CDC5L as bridge gene between chronic obstructive pulmonary disease and lung adenocarcinoma. Epigenomics 12, 1515–1529. doi: 10.2217/epi-2020-0112

PubMed Abstract | CrossRef Full Text | Google Scholar