Disease Activity-Associated Alteration of mRNA m5 C Methylation in CD4+ T Cells of Systemic Lupus Erythematosus

Epigenetic processes including RNA methylation, post-translational modifications, and non-coding RNA expression have been associated with the heritable risks of systemic lupus erythematosus (SLE). In this study, we aimed to explore the dysregulated expression of 5-methylcytosine (m⁵C) in CD4⁺ T cells from patients with SLE and the potential function of affected mRNAs in SLE pathogenesis. mRNA methylation profiles were ascertained through chromatography-coupled triple quadrupole mass spectrometry in CD4⁺ T cells from two pools of patients with SLE exhibiting stable activity, two pools with moderate-to-major activity, and two pools of healthy controls (HCs). Simultaneously, mRNA methylation profiles and expression profiling were performed using RNA-Bis-Seq and RNA-Seq, respectively. Integrated mRNA methylation and mRNA expression bioinformatics analysis was comprehensively performed. mRNA methyltransferase NSUN2 expression was validated in CD4⁺ T cells from 27 patients with SLE and 28 HCs using real-time polymerase chain reaction and western blot analyses. Hypomethylated-mRNA profiles of NSUN2-knockdown HeLa cells and of CD4⁺ T cells of patients with SLE were jointly analyzed using bioinformatics. Eleven methylation modifications (including elevated Am, 3′OMeA, m¹A, and m⁶A and decreased Ψ, m³C, m¹G, m⁵U, and t⁶A levels) were detected in CD4⁺ T cells of patients with SLE. Additionally, decreased m⁵C levels, albeit increased number of m⁵C-containing mRNAs, were observed in CD4⁺ T cells of patients with SLE compared with that in CD4⁺ T cells of HCs. m⁵C site distribution in mRNA transcripts was highly conserved and enriched in mRNA translation initiation sites. In particular, hypermethylated m⁵C or/and significantly up-regulated genes in SLE were significantly involved in immune-related and inflammatory pathways, including immune system, cytokine signaling pathway, and interferon signaling. Compared to that in HCs, NSUN2 expression was significantly lower in SLE CD4⁺ T cells. Notably, hypomethylated m⁵C genes in SLE and in NSUN2-knockdown HeLa cells revealed linkage between eukaryotic translation elongation and termination, and mRNA metabolism. Our study identified novel aberrant m⁵C mRNAs relevant to critical immune pathways in CD4⁺ T cells from patients with SLE. These data provide valuable perspectives for future studies of the multifunctionality and post-transcriptional significance of mRNA m⁵C modification in SLE.

Introduction

Accumulating evidence has revealed that post-transcriptional RNA modifications in mRNA may serve as novel gene expression regulators. Such epigenetic markers include N6-methyladenosine (m⁶A), N1-methyladenosine (m¹A), 5-hydroxymethylcytosine (hm⁵C), 5-methylcytosine (m⁵C), pseudouridine (Ψ), and inosine (Roundtree et al., 2017). Specifically, as a novel modified form of cytosine, m⁵C methylation of mRNAs in eukaryotes mainly depends on catalysis via the RNA methyltransferase NSUN2 (known as an m⁵C “writer”) (Abbasi-Moheb et al., 2012; Auxilien et al., 2012; Flores et al., 2017). Moreover, m⁵C epigenetic modification plays a significant role in eukaryotic biological functions such as nuclear export regulation and protein transcription modulation (Squires et al., 2012; Flores et al., 2017; Trixl and Lusser, 2019). Additionally, several studies have demonstrated that mRNA modification by the epitranscriptome marker m⁵C constituted an integral part of gene regulatory networks in both mammals and plants (Amort et al., 2017; Cui et al., 2017; Dominissini and Rechavi, 2017).

Systemic lupus erythematosus (SLE) is an autoimmune disease with complicated clinical manifestations in which disease activity plays a key role in therapy and prognosis (Morel, 2017). Loss of the normal balance of activity in CD4⁺ T cells is associated with the development of SLE (Yin et al., 2015). In addition, epigenetics-based genetic predispositions and their interaction in T cells are thought to contribute to SLE pathogenesis and development (Lei et al., 2009; Yeung et al., 2019); these epigenetic modulators include DNA methylation, DNA 5-hmC, and altered long non-coding RNA expression (Coit et al., 2015; Zhao et al., 2016; Zhang et al., 2018; Guo et al., 2019; Ye et al., 2019). For example, the transcription factor RFX1 affects the epigenetic status of CD4⁺ T cells, which result in autoimmune responses in SLE (Zhao et al., 2010).

Compared to those in healthy controls (HCs) group, 5-hmC levels in genomic DNA of CD4⁺ T cells from patients with SLE present different profiles, which results in the aberrant regulation of gene transcription in SLE pathogenesis (Zhao et al., 2016). Similarly, differentially expressed DNA methylation in CD4⁺ T cells of patients with SLE contributes to the clinical phenotypes of this disorder and may reflect the potential mechanism (Zhao et al., 2014). Accordingly, epigenetic modification profiling in SLE has proven beneficial toward obtaining a better understanding of patients with different disease activity (Lawrence et al., 2016; Meroni and Penatti, 2016). However, few meaningful results regarding the association between CD4⁺ T cells of SLE and m⁵C modifications in mRNA have been reported.

In this study, we utilized mass spectrometry to identify that m⁵C mRNA methylation in SLE CD4⁺ T cells serves as a crucial form of epigenetic modification. Next, we analyzed m⁵C methylation distribution and variation in mRNA by bisulfite sequencing in CD4⁺ T cells from patients with SLE presenting different disease activity compared to those in HCs. Our data demonstrated significant differences of m⁵C patterns according to disease activity and that the affected mRNAs are involved in disease-associated pathways. Furthermore, our data also provided a novel viewpoint for identifying key m⁵C-marked genes that may participate in the pathogenesis of SLE.

Materials and Methods

Subjects and Information

For this study, we recruited 47 subjects with SLE, as detailed in Supplementary Table S1, from The First Affiliated Hospital of Wenzhou Medical University between December 2019 and February 2020. All patients fulfilled the diagnostic criteria of the American College of Rheumatology (Petri et al., 2012). In addition, disease activity was assessed according to the systemic lupus erythematosus disease activity index (SLEDAI) (Uribe et al., 2004) at the time of blood collection prior to administering glucocorticoids and immunosuppressive agents to the patients. A total of 46 age- and sex-matched HCs without heart failure, renal failure, arthralgia, or rheumatic disease, and free from inflammatory conditions were also recruited from same hospital. For mRNA liquid chromatography-coupled triple quadrupole tandem mass spectrometry (LC-MS/MS), RNA sequencing (RNA-Seq), and RNA bisulfite sequencing (RNA Bis-Seq) analyses, 20 patients with SLE and 18 HCs were divided into six groups [among these, patients with SLE were grouped according to their SLEDAI scores: patients with SLEDAI score < 5 were divided into SLE stable (SA) groups (SA1, SA2) and those with SLEDAI score ≥ 9 were allocated to SLE moderate/major active (SM-MA) groups (SM-MA1, SM-MA2)]. For each group, five samples of the same SLEDAI score were mixed and sequenced; similarly, each control group comprised nine mixed and sequenced samples of the age- and sex-matched HCs. In the study of the effect of NSUN2 expression levels on RNA aspects, 22 patients with SLE [10 inactive (SLEDAI score < 5) and 12 active (SLEDAI score > 5)] and 23 HCs were used for quantitative reverse-transcription polymerase chain reaction (qRT-PCR) detection. For evaluating the effects of NSUN2 expression levels on protein aspects, 5 patients with SLE [1 inactive and 4 active] and 5 HCs were used for western blot detection. This study was approved by the Medical Ethical Committees of The First Affiliated Hospital of Wenzhou Medical University. All participants in this research provided written informed consent.

CD4⁺ T Cells and RNA Isolation

Peripheral blood mononuclear cells were isolated from HCs and patients with SLE using human peripheral blood lymphocyte separation medium (Tianjin Hao Yang Biological Manufacture, Tianjin, China) within 6 h of sample collection. Peripheral blood CD4⁺ T cells were then purified following immunomagnetic separation using a human CD4⁺ T cell isolation kit (BD Biosciences, San Jose, CA, United States) according to the manufacturer’s instructions. The isolated human primary CD4⁺ T cells were lysed in TRIzol reagent (Invitrogen Life Technologies, Grand Island, NY, United States) for RNA purification. The isolated RNAs were digested by DNase I (Invitrogen, Waltham, MA, United States) to remove residual DNA, then collected in 25 μL of DNase/RNase-free water. Isolated RNA was stored at −80°C for use.

Analysis of mRNA Modifications via LC-MS/MS

Total RNA was extracted from each sample using TRIzol reagent/RNeasy Mini Kit (Qiagen, Hilden, Germany). Total RNA was quantified and qualified using an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, United States), NanoDrop ND-1000 (Thermo Fisher Scientific Inc., Waltham, MA, United States), and 1% agarose gel. mRNA was isolated from total RNA using the NEBNext Poly(A) mRNA Magnetic Isolation Module (E7490; NEB, Ipswich, MA, United States) according to the manufacturer’s protocol. Purified mRNA was hydrolyzed to single nucleosides which were then dephosphorylated by enzyme treatment. Pretreated nucleosides solution was deproteinized using a Sartorius 10,000 Da molecular weight cut-off spin filter (Göttingen, Germany). Analysis of nucleoside mixtures was performed using an Agilent 6460 QQQ MS with ESI Jetstream ionization operated in positive ion mode, with SB-Aq 3.5 μm 2.1 × 150 mm high-performance LC column (Agilent Technologies) conditions according to the following solvent gradient [Solution A, high-performance LC-grade water with the relevant amount of formic acid to obtain a final formic acid concentration of 0.1% (vol/vol); Solution B, 100% acetonitrile with the relevant amount of formic acid to achieve a final formic acid concentration of 0.1% (vol/vol)]. The voltages and source gas parameters were as follows: gas temperature, 350°C; gas flow, 7 l min^–1; nebulizer, 40 psi; sheath gas temperature, 350°C; sheath gas flow, 11 l min^–1; capillary voltage, 3,500 V; and VCharging, 500 V. The molecular transition ions were quantified in multiple reaction monitoring mode (Supplementary Table S2). Multiple reaction monitoring peak information of modified nucleosides for each sample was extracted using Agilent Qualitative Analysis software. Peaks with signal-to-noise ratio ≥ 5 were considered as detectable nucleosides. Peak areas were then normalized to the quantity of purified mRNA of each sample.

Next Generation Sequencing for mRNA Gene Expression Profiling

Total RNA of each sample was extracted using TRIzol reagent/RNeasy Mini Kit/other kits. Total RNA (1 μg) with RNA integrity value > 6.5 was used for subsequent library preparation. Next generation sequencing library preparations were constructed according to the manufacturer’s protocol (NEBNext Ultra Directional RNA Library Prep Kit for Illumina). Differentially expressed transcripts were defined with a foldchange ≥ 1.2, and FDR ≤ 0.05, similar to previous studies (Chen X. et al., 2019). The detailed experimental procedure and analytic strategy can be found in our previous study (Guo et al., 2018).

Sequencing of Bisulfite-Converted RNAs (RNA-Bis-Seq) and Bioinformatics Analyses

Briefly, mRNA was purified from total RNA using the Dynabeads^® mRNA Purification Kit (Thermo Fisher Scientific Inc.). mRNA was bisulfite-converted and purified using the EZ RNA Methylation Kit (Zymo Research, Los Angeles, CA, United States). RNA libraries were then constructed using the NEBNext^® Ultra II Directional RNA Library Prep Kit (NEB) according to the manufacturer’s instructions. The library quality was evaluated using the BioAnalyzer 2100 system (Agilent Technologies, Inc.). Library sequencing was performed on an Illumina HiSeq instrument with 150 bp paired-end reads. Paired-end reads were harvested from the Illumina HiSeq sequencer and quality controlled based on a Q score of Q30, followed by 3′ adaptor-trimming and removal of poor-quality reads using cutadapt software (v1.9.3). Next, clean reads of busulfite-treated libraries were aligned to the reference genome (UCSC HG19) using meRanGs (a component of meRanTK) software with default parameters. The methylation status of each C within the genome was extracted using meRanCall (a component of meRankTK) software. meRanCompare (a component of meRanTK) software was used to identify differentially methylated sites (DMSs). The Ensembl genome features were used to annotate the methylated sites and DMSs. Only sites with coverage depth (methylated C number + non-methylated C number) ≥ 10, m⁵C methylation level ≥ 0.1, and methylated cytosine depth ≥ 5 were considered as credible m⁵C sites. Each group contained two sample replicates, and only overlapping m⁵C sites between two sample replicates were used to the following analyses. The m⁵C sites were annotated by applying BEDTools intersectBed (Quinlan and Hall, 2010). The distribution of m⁵C sites was analyzed as previously described (Yang et al., 2017). To acquire the sequence preference proximal to m⁵C sites, 21 nt sequences centered on each m⁵C site were extracted using BEDTools (Quinlan and Hall, 2010); and logo plots were generated using WebLogo¹ and motif enrichment analysis was done using MEME.², MetaPlotR software (Olarerin-George and Jaffrey, 2017) was used to map the distribution of methylation peaks in each sample on metagene. Moreover, meRanTK software (Rieder et al., 2016) was used to analyze the conversion rate of RNA-Bis-Seq; this could count the number of effectively covered C (and mutated C), the coverage of each C, and the number of Cs converted to Ts.

For the hypermethylated m⁵C transcripts in SLE samples, transcripts with differences in mean m⁵C level ≥ 0.05 between the groups (SA and HC; SM-MA and HC) were considered to be statistically significant. To further explore the critical role of m⁵C modification in SLE, transcripts were separated into four groups based on whether the modifications resulted in hyper- or hypomethylation of m⁵C, and based on the up- or down-regulation of gene expression. Only the intersection of significant hyper- or hypomethylated transcripts, as well as significant DEGs (| FC| > 1.2, FDR < 0.05), was subjected to downstream analysis (Chen X. et al., 2019). Metascape³ is an effective and efficient tool for experimental biologists to comprehensively analyze and interpret OMICs-based studies in the big data era (Zhou et al., 2019). For the list of hyper-/hypomethylated genes, pathway and process enrichment analysis was carried out using the following ontology sources: Kyoto Encyclopedia of Genes and Genomes Pathway, Gene Ontology Biological Processes, Reactome Gene Sets, Canonical Pathways, and CORUM. All genes in the genome were used as the enrichment background. Terms with p-value < 0.01, minimum count of 3, and enrichment factor > 1.5 (ratio between the observed counts and the counts expected by chance) were collected and grouped into clusters based on their membership similarities. To further capture the relationships between the terms, a subset of enriched terms was selected and rendered as a network plot, where terms with a similarity > 0.3 were connected by edges. We selected the terms with the best p-values from each of the 20 clusters. Moreover, for the hypermethylated genes, protein–protein interaction enrichment analysis was also carried out using the following databases: BioGrid, InWeb_IM, OmniPath. The resultant network contained the subset of proteins that form physical interactions with at least one other member in the list. If the network contained between 3 and 500 proteins, the Molecular Complex Detection (MCODE) algorithm was applied to identify densely connected network components. Reactome Pathway significant enrichment analysis was performed to analyze the Hyper-Up related gene list or potential target genes for NSUN2 using the STRING database⁴ (Szklarczyk et al., 2017). Significance was determined based on false discovery rate < 0.05.

RT-qPCR

NSUN2 expression levels in the CD4⁺ T cells from 45 subjects (22 patients with SLE and 23 HCs) were detected by qRT-PCR using the QuantiNova SYBR Green PCR Kit (Qiagen). All qRT-PCR reactions were carried out using an Applied Biosystems QuantStudio^TM 3 Real-Time PCR Instrument (Thermo Fisher Scientific Inc.). For each reaction, 1 μL of diluted cDNA was mixed with 10 μL of 2 × SYBR Green PCR Master Mix. A final volume of 20 μL was achieved by the addition of 1.4 μL forward and reverse primers (10 μmol). The conditions for PCR amplification were as follows: 95°C for 2 min, followed by 40 cycles of 95°C for 5 s and 60°C for 10 s. The specificity of the primer amplicons was tested by melting curve analysis. All samples were tested in triplicate. The data were analyzed using the comparative threshold cycle (Ct) method. GAPDH was used as a control, and the relative quantification of NSUN2 in CD4⁺ T cells was calculated using the following equation: amount of target = 2^–Δct, where ΔCt = Ct_NSUN2 – Ct_GAPDH. The following gene-specific primers were used for qRT-PCR analysis: NSUN2: 5′-GAACTTGCCTGGCACACAAAT-3′ and 5′-TGCTAACAGCTTCTTGACGA CTA-3′, and GAPDH: 5′-CAGGGCTGCTTTTAACTCTGGTAA-3′ and 5′-GGGTG GAATCA TATTGGAACATGT-3′.

Cell Culture and Western Blot

HEK293T [American Type Culture Collection (ATCC^§), Manassas, VA, United States] and Jurcat cells were, respectively, cultured in Dulbecco’s modified Eagle’s medium and Roswell Park Memorial Institute 1640 supplemented with 10% fetal bovine serum at 37°C in a humidified atmosphere containing 5% CO₂ (SANYO, MCO-175, Osaka prefecture, Japan). HEK293T and Jurcat cells and CD4⁺ T cells from five HCs and five patients with SLE were lysed using protein lysis buffer (Beyotime Institute of Biotechnology, Beijing, China) supplemented with protease inhibitor cocktail (Pierce, Rockford, IL, United States) at 4°C for 20 min. Protein samples were separated using 10% sodium dodecyl sulfate-polyacrylamide gel electrophoresis and then electrophoretically transferred to polyvinylidene difluoride membranes (Millipore, Billerica, MA, United States). Anti-NSUN2 antibody (20854-1-AP; Proteintech, Rosemont, IL, United States) was diluted with primary antibody dilution buffer (Beyotime Institute of Biotechnology) to 1:5,000, and anti-GAPDH antibody (GOOD HERE, Hangzhou, China) was also diluted to 1:1,000. The membranes were then washed with TBST buffer five times for 5 min each and incubated with horseradish peroxidase-conjugated goat anti-rabbit IgG secondary antibody (1:5,000 dilution) (Cell Signaling Technology, Danvers, MA, United States) for 1.5 h at 37°C. Bands were detected using enhanced chemiluminescence and visualized with a Gel Doc 2,000 (BioRad, Hercules, CA, United States).

Statistical Analysis

Statistical analysis was performed using SPSS 22.0 software (IBM, Armonk, NY, United States) and GraphPad Prism version 8.0.1 software (GraphPad Software, La Jolla, CA, United States). An independent sample t-test was performed to analyze the difference of the mRNA expression levels in CD4⁺ T cells of patients with SLE and HCs; the data are presented as the means ± standard deviation. A non-parametric Mann–Whitney U-test was performed to analyze the difference of the mRNA m⁵C methylated levels or protein expression levels in CD4⁺ T cells of patients with SLE and HCs; the data are presented as 50% quantile (25% quantile, 75% quantile). A two-sided p < 0.05 was considered to represent a statistically significant difference.

Results

mRNA Methylation Profiling of HCs and Patients With SLE Presenting Diverse Disease Activity

We isolated mRNA from the CD4⁺ T cells of 10 patients with SLE exhibiting stable activity (SA group), 10 patients with moderate/major activity (SM-MA group), and 18 HCs (HC group), then combined equal amounts of mRNA from 5 or 9 individuals, respectively, into one pool for each group. Finally, each group consisted of two separate pools for analysis. We generated mRNA methylomes for six separate pools using LC-MS/MS and identified that the mRNA levels were differently modified between HCs and patients with SLE exhibiting diverse disease activity (Figures 1A,B and Supplementary Table S2). A total of 11 modifications (including m⁵C, Ψ, m⁶A, and m¹A) previously identified in mRNA were detected in our study among these groups. Compared with those of HCs, the Am, 3′OMeA, m¹A, and m⁶A levels in CD4⁺ T cells of SLE were elevated, whereas those of m⁵C, Ψ, m³C, m¹G, m⁵U, and t⁶A were decreased (Figure 1A and Supplementary Figure S1). As it has been reported that m⁵C is a newly discovered internal mRNA modification in eukaryotes (Amort et al., 2017) that regulates immune response including oncogene activation (Chen X. et al., 2019), in this study, we further focused on the m⁵C level in overall mRNA. Compared to those in HCs, the m⁵C/C levels in CD4⁺ T cells were markedly lower in both SA and SM-MA groups (Figure 1C). Moreover, the m⁵C/C levels in CD4⁺ T cells were decreased in patients with SLE exhibiting increasing disease activity.

FIGURE 1

Figure 1. Detection of mRNA modifications by LC-MS/MS among healthy controls (HCs) and systemic lupus erythematosus (SLE) patients with different disease activity. (A) Heatmap of normalized abundance (modification/canonical nucleotide) of 11 mRNA modifications detected by LC-MS/MS between HCs and SLE patients. Red indicates a high z-score, whereas blue indicates a low z-score. (B) LC-MS/MS extracted ion chromatograms of modified nucleotides analyzed in CD4⁺ T cells mRNA of HCs and SLE patients. S₄U, 4-thiouridine; m³C, 3-methylcytosine; m⁵C, 5-methylcytosine; m⁵U, 5-methyluridine; m¹A, N¹-methyladenosine; Ψ, pseudouridine; m⁶A, N6-methyladenosine; 3′-OMeA, 3’-O-methyladenosine; Am, 2’-O-methyladenosine; m¹G, 1-methylguanosine; t₆A, N6-threonylcarbamoyladenosine. (C) Global dynamics of calibrated m⁵C/C levels in HCs and SLE patients with different disease activity (SA, SM-MA). Bars show the mean of individual biological replicates (n = 2).

Distribution Profiling of m⁵C in mRNA of Patients With SLE Exhibiting Different Disease Activity and HCs

To obtain a transcriptome-wide landscape of m⁵C profiling, we further performed mRNA Bis-Seq analysis on mRNA samples purified from CD4⁺ T cells of patients contributing to the SA, SM-MA, and HCs pools according to a recently described study (Yang et al., 2017). The overlapping m⁵C sites in two independent pools from each group were selected for follow-up analysis. For example, a total of 233 m⁵C sites identified in both SM-MA patient replicates (“high-confidence” set) were used in subsequent bioinformatics analyses (Figure 2A and Supplementary Table S3). Overall, the m⁵C levels (approximately 62.8%) in mRNA of CD4⁺ T cells of HCs were considerably higher compared with those from both SA and SM-MA groups (Figure 2B), as determined by LC-MS/MS analyses. Furthermore, the overall m⁵C level in mRNA of CD4⁺ T cells from the SM-MA group (19.6%) was relatively lower than that in the SA group (25.3%). Notably, the number of m⁵C-modified mRNA molecules exhibited opposite changes to the number of m⁵C-modified sites with increasing disease activity (Figure 2C). Among the m⁵C sites/mRNAs identified in CD4⁺ T cells of SA and SM-MA groups, more m⁵C-containing gene transcripts were observed with fewer m⁵C methylation sites (297/158 and 233/186, respectively) within mRNAs (Figure 2C), which is contrary to the results in CD4⁺ T cells of HCs (2436/81).

FIGURE 2

Figure 2. Distribution profiles of m⁵C in mRNAs from systemic lupus erythematosus CD4⁺ T cells. (A) Venn diagram showing overlap of m⁵C sites within mRNAs between two SLE moderate/major active (SM-MA) pool replicates. (B) Bar chart showing level of m⁵C sites among healthy controls (HCs) and SLE stable (SA) and SLE moderate/major active (SM-MA) patients (***p < 0.001). (C) Bar chart showing number of m⁵C sites and m⁵C-containing mRNA among HCs and SA and SM-MA patients. (D) Venn diagram showing overlap of genes within m⁵C-containing mRNAs in HCs and SA and SM-MA patients. (E) Proportion of genes harboring different numbers of m⁵C sites in three groups. The majority of genes harboring only one m⁵C site except that in the HC group. (F) Bar chart showing percentage of mRNA m⁵C sites identified in each sequence context (CG, CHG, and CHH, where H = A, C, or U) in CD4⁺ T cells of HCs and SA and SM-MA patients. (G) Distribution of m⁵C sites along mRNA transcripts. The moving averages of percentages of mRNA m⁵C sites are shown (raw data; only HC2 and SM-MA2 group are shown here) (left). The averages of percentages of mRNA m⁵C sites within 5′UTR, StartC, CDS, stopC, and 3′UTR regions in transcriptomes are shown. Distribution of m⁵C sites peaks along mRNA transcripts are shown (right). (H) Sequence frequency logo for the sequences proximal to mRNA m⁵C sites among HCs and SLE patients. (I) Top differential enrichment mode motifs (SA vs. HC; SM-MA vs. HC) enriched across the sequences proximal to mRNA m⁵C sites identified from HCs and SA and SM-MA patients (only E-value < 0.05 was shown here).

A Venn diagram was used to analyze the overlap of genes within m⁵C-containing mRNAs in HC, SA, and SM-MA groups. An overlap of 13 m⁵C mRNAs between the groups was found as shown in Figure 2D. Additionally, 46 m⁵C-containing mRNAs were detected in only HCs, 93 only in the SA group, and 129 only in the SM-MA group. Next, the proportion of genes harboring different numbers of m⁵C sites in the three groups was further analyzed. Compared with that of HCs, the percentage of m⁵C-containing gene transcripts harboring one m⁵C site was obviously higher in CD4⁺ T cells of SA and SM-MA groups (44.4% vs. 79.1% and 84.4%, respectively), whereas the percentage of transcripts harboring three m⁵C sites was significantly lower (50.6% vs. 4.4% and 1.1%) (Figure 2E). In turn, the m⁵C-containing gene transcript percentage containing two or three m⁵C sites in SLE was slightly higher than that in HCs. Moreover, the percentage of mRNA m⁵C sites identified in each sequence context (CG, CHG, and CHH, where H = A, C, or U) in HCs and patients with SLE exhibiting different disease activity markedly differed, with an increase in CG and CHG regions, respectively, from 10.2% and 18.3% in HCs to 24.2% and 31.3% in the SA group, further increasing to 28.3% and 46.4% in the SM-MA group; conversely, a sharp decrease was observed in CHH regions from 71.4% in HCs to 44.4% in SA, then to 25.3% in SM-MA groups (Figure 2F).

We next determined the enrichment feature of m⁵C sites in mRNA transcripts. m⁵C peaks were divided into 5′-untranslated region, StartC, coding sequence (CDS), StopC, and 3′untranslated regions according to their locations in RNA transcripts (Figure 2G). Notably, m⁵C sites were also enriched in startC and CDS regions immediately downstream of translation initiation sites (Figure 2G), which is similar to previous observation in HeLa cells (Yang et al., 2017). Similarly, a sequence frequency logo identified that m⁵C sites were embedded in CG-rich environments (Figure 2H). The most significantly enriched motifs in patients with SLE and HCs were further analyzed and presented in Supplementary Figure S2. Moreover, the top differential enrichment mode motifs (SA vs. HC; SM-MA vs. HC) enriched across the sequences proximal to mRNA m⁵C sites were further identified (Figure 2I). These were BYHBCWGSCCRGRRA (p = 1.0e−013) and GGVGMGGG (p = 3.1e−002) between SA and HC groups, which were present in approximately 49 and 5% of the methylated sequences proximal to mRNA m⁵C sites, respectively (Figure 2I). Between SM-MA and HC groups, CAGGRRRRNDR (> 75% of the methylated sequences proximal to mRNA m⁵C sites) was most significantly enriched (p = 3.7e−018) (Figure 2I), which indicated that CCRGRRA and CAGGRR might constitute the specific motif of m⁵C-modified genes in SA and SM-MA groups, respectively.

Pathway Analysis of the Genes Containing Significantly Altered m⁵C Sites (Differentially Methylated Genes) Between Patients With SLE and HCs

Accumulating evidence has shown that numerous mRNAs with m⁵C methylation promote the pathogenesis of human diseases such as urothelial carcinoma (Chen X. et al., 2019). To provide functional insights regarding whether mRNAs carrying m⁵C methylation in CD4⁺ T cells are linked to SLE, we explored the overall changes of transcript m⁵C methylation by mapping the m⁵C distribution in CD4⁺ T cells through mRNA-Bis-Seq data between SA/SM-MA and HC groups. We found 78 m⁵C hypomethylated transcripts and 131 m⁵C hypermethylated transcripts in CD4⁺ T cells of the SA group compared to those in HCs (Figure 3A and Supplementary Table S4). Moreover, compared to HCs, 80 transcripts with decreased m⁵C levels and 166 transcripts with increased m⁵C levels were found in CD4⁺ T cells of the SM-MA group (Figure 3B and Supplementary Table S4).

FIGURE 3

Figure 3. Pathway analysis of genes containing significantly altered m⁵C sites (differentially methylated genes, DMGs) between healthy controls (HCs) and systemic lupus erythematosus (SLE) patients. Correlation of m⁵C methylation between SLE stable (SA) (A) or SLE moderate/major active (SM-MA) (B) group with HC group. Regression lines and 95% confidence intervals (shaded regions) are shown. Spearman’s correlation coefficients (R) were calculated using the “cor.test” function in the statistical language R. Hypermethylation (red points) indicates that the m⁵C transcripts with altered methylation level in SLE patients is higher than that in HCs. Hypomethylation (blue points) indicates that the m⁵C transcripts with altered methylation level in SLE patients are lower than that in HCs. (C,D) Network of enriched terms colored by cluster ID for differentially hypermethylated genes (C, SA vs. HC; D, SM-MA vs. HC). Each term is represented by a circle node, where its size is proportional to the number of input genes that are related to that term, and its color represents its cluster identity. Terms with a similarity score > 0.3 are linked by an edge (thickness of the edge represents the similarity score). (E,F) The molecular complex detection (MCODE) components identified for the differentially hypermethylated genes (E, SA vs. HC; F, SM-MA vs. HC). Pathway and process enrichment analysis has been applied to each MCODE component independently, and the three best-scoring terms according to p-value have been retained as the functional description of the corresponding components. (G,H) Bar graph of enriched terms across the differentially hypomethylated genes (E, SA vs. HC; F, SM-MA vs. HC), colored according to p-values.

Additionally, we further investigated the functions of differentially methylated genes in patients with SLE using Metascape. m⁵C hypermethylated genes in both SA and SM-MA groups were found to be significantly enriched in immune-related pathways, including cytokine-mediated signaling pathway and homeostasis of number of cells (Figures 3C,D and Supplementary Table S5), which consistent with the function of m⁵C-modified transcripts in HEK293 cells (Sun et al., 2019). Among these, MCODEs analyses of hypermethylated genes in SA/SM-MA groups for the top pathways and process enrichment result identified mRNA splicing in the protein–protein interaction networks of both SA versus HC and SM-MA versus HC (Figures 3E,F), which is similar to the results of a previous study (Yang L. et al., 2019). Gene ontology enrichment analysis was also performed for m⁵C hypomethylated transcripts between SA/SM-MA and HC groups. Overall, m⁵C hypomethylated transcripts in SA or SM-MA groups were found to be significantly enriched in eukaryotic translation elongation and protein methylation (Figures 3G,H and Supplementary Table S6).

Functional Enrichment of Differentially Expressed Genes (DEGs) Between Patients With SLE and HCs

As shown in Supplementary Figure S3A, principal component analysis results of mRNA profiles indicated that the gene expression patterns were similar within the same group, however, significantly different between HC, SA and SM-MA groups. In SA versus HCs and SM-MA versus HCs, gene ontology analysis showed that the DEGs in both groups were significantly enriched in some essential molecular functions (e.g., protein binding and chemokine activity), cellular components (e.g., plasma membrane, cytosol and extracellular region and space), and biological processes (e.g., immune response, innate immune response, cytokine-mediated signaling pathway, and type I interferon signaling pathway) (Supplementary Figure S3B), which was similar to our previous study on peripheral blood mononuclear cells of SLE (Guo et al., 2018). Moreover, pathway analysis of DEGs in the SA versus HC and SM-MA versus HC groups showed that DEGs in both comparisons were linked to inflammatory pathways and immune-related pathways in organismal systems (e.g., chemokine signaling, Toll-like receptor signaling, and IL-17 signaling pathway), human disease (e.g., rheumatoid arthritis and SLE), and environmental information processing (e.g., cytokine–cytokine receptor interaction and TNF-signaling pathway) (Supplementary Figure S3C).

Integration Analyses of m⁵C-Containing mRNA Methylation and mRNA Transcript Expression

As previously reported, the up-regulation of numerous oncogene RNAs with hypermethylated m⁵C sites has been causally related to human bladder cancer (Chen X. et al., 2019). Here, we further explored the association between RNAs with m⁵C modification and mRNA transcript expression by integrated RNA-Seq and RNA-Bis-Seq data. We found no strong correlation between m⁵C-containing mRNA methylation and mRNA transcript expression levels in patients with SLE (Figure 4A), which was consistent with previous report (Wang et al., 2019). Moreover, mRNA m⁵C modification has been reported to up-regulate mRNA expression by stabilizing mRNA (Chen X. et al., 2019). Toward this end, we also identified 27 and 52 up-regulated expressed mRNAs exhibiting hypermethylated m⁵C in the CD4⁺ T cells of SA versus HCs and SM-MA versus HC groups, respectively (Figure 4A and Supplementary Table S7). Notably, ADAR (Roth et al., 2018), an SLE-related gene, was identified in these data as exhibiting both m⁵C hypermethylation and mRNA up-regulation in patients with SLE (Supplementary Table S7). Moreover, we further investigated the enrichment of signaling pathways of up-regulated expressed RNAs exhibiting hypermethylated m⁵C levels using Reactome. In particular, these transcripts were enriched prominently in categories of immune system, neutrophil degranulation, cytokine signaling in immune system, and interferon signaling (Figure 4B and Supplementary Table S8), which indicated that up-regulated genes with high m⁵C-levels modification participate in the flares and remission of patients with SLE through disruption of the patient’s immune system.

FIGURE 4

Figure 4. Conjoint analysis of differentially methylated genes and differentially expressed genes. (A) Distribution of mRNAs with a significant change in both m⁵C modification and gene expression levels in healthy controls (HCs) and SLE stable (SA)/SLE moderate/major active (SM-MA) patients, respectively. (B) Reactome pathways analysis of hypermethylated mRNAs with up-regulated genes in SA versus HC (left) or SM-MA versus HC (right) groups, respectively.

Expression Level and Regulation of NSUN2 in CD4⁺ T Cells From Patients With SLE

As NSUN2 is the only known mRNA m⁵C methyltransferase (Li et al., 2017), we further explored the expression level of NSUN2 in CD4⁺ T cells. Compared with that in HCs, its mRNA and protein expression in CD4⁺ T cells from patients with SLE was significantly decreased (p < 0.01 and p < 0.05, respectively) (Figures 5A,B). Similar to the down-regulated m⁵C methylation levels in transcripts observed in siNSUN2-HeLa cells (Yang et al., 2017), we also identified 15 m⁵C hypomethylated transcripts in both SA versus HC and SM-MA versus HC groups (Figure 5C). Reactome pathways analysis revealed that these 15 transcripts were mainly enriched in representative transcription-related pathways including eukaryotic translation elongation and termination, peptide chain elongation, and mRNA translation and metabolism (Figure 5D), which indicted that these abnormally hypomethylated transcripts likely contribute to the risk of developing SLE.

FIGURE 5

Figure 5. Down-regulation of NSUN2 is involved in the pathogenesis of systemic lupus erythematosus. (A) NSUN2 mRNA expression level in CD4⁺ T cells from SLE patients (n = 22) and healthy controls (HCs) (n = 23). Data are presented as 2^–Δct relative to GAPDH expression (mean ± standard deviation; **p < 0.01). (B) NSUN2 protein expression levels in CD4⁺ T cells from five SLE patients, five HCs, and two positive control cell lines (Jurkat and HEK293T). The gray scale of WB was measured by ImageQuant software and compared between SLE patients and HCs (*p < 0.05). NSUN2 protein expression was normalized to that of GAPDH. (C) Overlap of differentially hypomethylated genes in siNSUN2-HeLa cells with SLE stable (SA) versus HC and SLE moderate/major active (SM-MA) versus HC hypomethylated mRNAs. (D) Reactome pathways analysis of potential methylated target genes of methylase NSUN2.

As suggested in the model shown in Figure 6, NSUN2-down-regulated expression in SLE contributes to the down-regulation of m⁵C levels. Among the affected transcripts, hypermethylated m⁵C transcripts in CD4⁺ T cells of patients with SLE mainly participated in the cytokine-related signaling pathway, and mRNA splicing, stabilization, and translation, whereas hypomethylated m⁵C transcripts were mainly involved in translation elongation in the pathogenesis of SLE.

FIGURE 6

Figure 6. Proposed model for the landscape of NSUN2-mediated regulation of m⁵C methylation in the pathogenesis of systemic lupus erythematosus. m⁵C formation is promoted by NSUN2. Genes with dysregulated m⁵C methylation participate in the flares and remission course of SLE.

Discussion

Recent studies have identified that DNA hypomethylation might contribute to the immune response and over-reactivity in SLE CD4⁺ T cells (Zhang et al., 2013; Zhao et al., 2014; Yeung et al., 2019). In addition, biochemical modifications to mRNA, especially m⁶A and m⁵C, have also been shown to be of considerable importance with regard to the biological functions of developing disease (Li et al., 2018; Chen X. et al., 2019). As a newly discovered methylation modification in eukaryotic mRNAs, m⁵C methylation was reported as an unprecedented mechanism of regulating oncogene activation in bladder carcinoma (Chen X. et al., 2019). Nevertheless, the regulation and biological role of m⁵C modification in SLE pathogenesis remains unclear. Here, we displayed transcriptome-wide m⁵C modification profiles of SLE CD4⁺ T cells with different disease activity, investigating gene expression and disease-related pathways. As the number of CD4⁺ T cells able to be isolated from the peripheral blood is limited and patients with SLE exhibit heterogeneity, we pooled RNA of CD4⁺ T cells from 5/9 subjects equally into a single pool. As reported previously, for epigenetic signs, RNA pooling allowed the effective focus on actual epigenetic patterns under the background of the numerous random inter-individual differences (Docherty et al., 2009; Toperoff et al., 2012; Zhao et al., 2014). Additionally, sample pooling ensures the ability to obtain multi-omics information from the same samples and precise detection of average RNA methylation and expression levels in large groups of mRNAs from individual patients. And PCA results of mRNA profiles indicated that pooled samples are valid for the analysis of methylation profiling.

In this study, we revealed the modifications in global methylated mRNA in SLE CD4⁺ T cells for the first time, illuminating multiple patterns of mRNA modification in SLE. We also uncovered the primary differences of m⁵C modifications including prevalence, number, and unique distribution along transcripts in HCs and patients with SLE. Notably, we identified that m⁵C levels were decreased in CD4⁺ T cells of patients with SLE compared with those of HCs along with lower m⁵C levels associated with severe disease activity in patients with SLE. Consistent with DNA hypomethylation in SLE, the results also indicated that m⁵C hypomethylation in SLE CD4⁺ T cells might increase T cell self-reactivity in patients with SLE. In comparison, a recent study reported that the levels of DNA hydroxymethylation, transited from DNA m⁵C, were increased and some genes exhibited promoter region hypomethylation in CD4⁺ T cells from patients with SLE (Zhao et al., 2016), which may explain the potential mechanism of m⁵C in SLE. Notably, in SLE patients, we observed increased m⁵C-containing transcripts with low modification levels in CD4⁺ T cells, primarily containing one or two m⁵C sites, among which certain mRNAs were determined to likely be involved in SLE pathogenesis. It should be noted that the overall modification level of one sample is primarily determined by the number of m⁵C modified transcripts in that sample, as well as the level of m⁵C modification of each transcript. Therefore, an increase in the number of m⁵C modified transcripts in one sample does not necessarily indicate an overall higher level of modified m⁵C.

Notably, our analysis of the sequences proximal to mRNA m⁵C sites revealed several hitherto unknown features of m⁵C distribution in SLE mRNA. First, m⁵C sites are located in the vicinity of the translation start sites of mRNAs, and these patterns of distribution are highly conserved in patients with SLE. A similar pattern of distribution of m⁵C sites has been reported in HeLa cells (Yang et al., 2017), cancers (Chen X. et al., 2019), and mouse embryonic stem cells and brain tissues (Amort et al., 2017), whereas the patterns differ in Arabidopsis (Cui et al., 2017). Second, m⁵C sites were observed to be enriched in the GC-rich region with the highest enrichment motifs of m⁵C sites being CCRGRRA and CAGGRR. This may be relevant as the regulatory effects of RNA modification are usually interpreted by their reader proteins (Roundtree et al., 2017). However, to date, only three known m⁵C binding proteins are known to recognize the m⁵C modified mRNAs: the nuclear ALYREF promotes mRNA export (Yang et al., 2017), and the cytoplasmic YBX1 and TRM4B maintain mRNA stability in mammals and plants (Cui et al., 2017; Chen X. et al., 2019). Moreover, the most enriched motifs for m⁵C methylated peak summits in Arabidopsis thaliana were HACCR and CTYCTYC (Cui et al., 2017), which different from those in patients with SLE. The specific motifs in these patients implied the presence of other as-yet unknown potential reader(s), which needs to be validated though additional study. Nevertheless, these distribution features, together with the observation of low m⁵C levels in mRNA with high disease activity, imply that m⁵C modification in CDSs may play a crucial role in mediating the flares and remission of SLE.

In additional, compared with the HC group, we identified 209 (including 131 hypermethylated and 78 hypomethylated transcripts) and 246 (including 131 hypermethylated and 78 hypomethylated transcripts) dysregulated m⁵C-methylated transcripts in the SA and SM-MA group, respectively. These abnormally methylated transcripts in SLE CD4⁺ T cells might contribute to SLE pathogenesis. CD4⁺ T cells constitute the principle components of the adaptive immune system that can secrete cytokines; moreover, naive CD4⁺ T cells can differentiate into a variety of different cell subsets. However, dysregulation of CD4⁺ T cell differentiation is crucial to the pathogenesis of SLE. Increased TGF-β and IL-6 serve as major contributors to the imbalance of Th17 and regulatory T cell differentiation (Raphael et al., 2015), which promotes SLE development. Consistent with this, our transcription-wide analysis identified transcripts with increased m⁵C in patients with both stable and active SLE as being involved in several important inflammatory biological process and immune-related pathways including the cytokine-mediated signaling pathway, homeostasis of number of cells, and interferon response signaling, which suggested that transcripts exhibiting m⁵C hypermethylation may play a very vital role in the pathogenesis of SLE. For example, cytokines (such as interleukin-6, interleukin-10, and tumor necrosis factor) in patients with SLE can contribute to SLE susceptibility through the promotion of autoantibody production and inflammation (Tsokos et al., 2016). Moreover, in this study, although the methylated cytosine depth of both transforming growth factor B1 and interleukin-6 receptor transcripts were < 5, both were m⁵C hypermethylated in CD4⁺ T cells of SA/SM-MA groups, which might stabilize their mRNA expression and lead to the imbalance of Th17/regulatory T cell differentiation and, thus, likely contribute to SLE pathology.

Interferon signaling has also been reported to constitute a crucial driver of the development of autoimmunity in lupus mice (Ikeda et al., 2017; Zeng et al., 2019) and demonstrates close association with SLE disease activity in human (Ronnblom and Elkon, 2010). Furthermore, homeostasis of the number of cells and the immune system is also responsible for the development of autoimmune diseases (Saferding and Bluml, 2019). In addition, a common major pathway of hypermethylated m⁵C transcripts is mRNA splicing. Recently, growing evidence has supported differential splicing of mRNAs as a mechanism in T cells (Schaub and Glasmacher, 2017) that is linked to presentation of immune-related diseases, considering that some transcripts with decreased m⁵C methylation level were significantly enriched in eukaryotic translation elongation. As inflammation is controlled by transcriptional regulation including the pathogen-associated molecular process of mRNA translation elongation (Bianco et al., 2019), these findings suggest that m⁵C modification may serve to remodel the transcriptional regulation in patients with SLE.

As reported previously, m⁵C-containing transcript modification and mRNA expression levels show no strong correlation in bladder cancers (Chen X. et al., 2019), as was also observed in patients with SLE in this study. Through integration with our previous data (Guo et al., 2018), we found that some dysregulated mRNAs likely participate in the immune system in either peripheral blood mononuclear cells or CD4⁺ T cells, whereas some transcripts significantly enriched in chemokine activity processes were detected only in CD4⁺ T cells. This may be due to the distinctions in disease phenotype in the sample pool, sample size, and sample type. As mRNA m⁵C modifications may lead to up-regulated mRNA expression by stabilizing mRNA in the immune system (Chen X. et al., 2019), we focused on up-regulated transcripts with increased m⁵C levels. We identified that significantly up-regulated transcripts exhibiting hypermethylated m⁵C levels in patients with both stable and active SLE were associated with SLE-related biological process and immune-related pathways, such as the immune system, cytokine signaling pathway, and interferon signaling. These pathways are the same as the major pathways in transcripts with increased m⁵C modification are involved, indicating that up-regulated transcripts with increased m⁵C levels may play a central part in the flares of patients with SLE. However, further functional research will be needed to clarify the functional connection between RNA m⁵C modification and gene expression in the development of SLE.

Recently, m⁵C methylation in eukaryotic cells has been identified to largely depend on the function of the known “writer” NSUN2, which functions as an mRNA methyltransferase (Bohnsack et al., 2019; Yang Y. et al., 2019). Moreover, it was reported that NSUN2 affects cell biological processes including cell proliferation, migration, and mRNA metabolism by regulating mRNA m⁵C modification (Wang et al., 2011; Courtney et al., 2019). Compared to that in HCs, the NSUN2 expression (including mRNA and protein) in CD4⁺ T cells of patients with SLE was significantly lower, which further confirmed the Bis-Seq data of overall m⁵C levels in SLE. Consistent with our results, Han et al. (2020) found NSUN2 up-regulation, as well as aberrant mRNA m⁵C gain and loss in PM2.5-induced pulmonary fibrosis mouse models. This phenomenon has also been described in bladder cancer and Zebrafish (Chen X. et al., 2019; Yang Y. et al., 2019). In summary, the results for aberrant methylated modifications, as determined by RNA-Bis-Seq, are in part due to the down-regulation of NSUN2; however, we cannot rule out the presence or role of other m⁵C methyltransferases in SLE patients resulting in the generation of more m⁵C-containing transcripts. Additionally, this provided indirect evidence that the abnormal NSUN2 may constitute a good target associated with m⁵C modification in SLE CD4⁺ T cells. Through comparison with the RNA-Bis-Seq data of NSUN2-knockdown HeLa cells (Yang et al., 2017), we found that 15 transcripts were similarly m⁵C hypomethylated in CD4⁺ T cells of SLE, the majority of which primarily participated in eukaryotic translation elongation and termination and RNA metabolism. This observation was consistent with previous studies (Dominissini and Rechavi, 2017; Yang et al., 2017; Chen X. et al., 2019), indicating that the majority of core hypomethylated m⁵C genes in SLE are influenced by NSUN2 expression. However, to some extent, CD4⁺ T cells and HeLa cells differ in the aspects of molecular function, biological process, and diseases models. The relationship between NSUN2-mediated RNA methylation and gene expression in patients with SLE is still unclear. Further studies on immune cells and animal models will be required to obtain a common understanding of the role of NSUN2 in patients with SLE.

Moreover, based on the integrated analysis of m⁵C-containing mRNA methylation and mRNA transcript expression, 5-methylcytosine is expected to be confirmed as a potential and pivotal modification in the pathogenesis of SLE. Consequent to NSUN2 dysregulated expression, m⁵C methylation levels were down-regulated in SLE. In these patients, RNAs exhibiting m⁵C hypomethylation, enriched in translation initiation regions, participated in translation elongation, whereas hypermethylated mRNAs were involved in SLE-related and cytokine-mediated signaling pathways, and mRNA splicing, stabilization, and translation. It has been found that several mRNA m⁶A methyltransferases (including METTL3/14, WTAP, RBM15/15B, etc.) are relevant to installation of m⁶A modification in mRNA (Chen X.Y. et al., 2019). In addition to NSUN2, we speculated that other writer(s) might exist and regulate m⁵C modification in SLE, resulting in the generation of more m⁵C-containing transcripts. However, the small sample size and heterogeneous sample pools constituted limitations in this study. For example, although a large proportion of the modifications exhibited no distinct difference in the global methylation profile of the HC group, several modifications (Ψ, m⁵U, m¹G and s⁴U) were found to differ from each other. Subsequent studies using larger sample sizes will be needed to validate these overall modification levels in SLE. In addition, considering that the most abundant modification was that of m⁶A, its mRNA should be jointly analyzed with m⁵C in the context of SLE pathogenesis.

Conclusion

In summary, this study presented the first transcriptome-wide m⁵C methylated map of CD4⁺ T cells from patients with SLE and HCs. Our results offered a potential link between abnormal m⁵C RNA modifications and immune-related pathways. This novel epigenetic modification may provide a different perspective for better understanding the pathogenesis of SLE. Moreover, targeting m⁵C modifications is likely to become a therapeutic strategy in SLE in the future.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation, to any qualified researcher.

Ethics Statement

The studies involving human participants were reviewed and approved by the Medical Ethical Committees of The First Affiliated Hospital of Wenzhou Medical University. The patients/participants provided their written informed consent to participate in this study.

Author Contributions

GG and HW performed the experiments, analyzed and interpreted the data, and drafted the manuscript. XS performed the experiments and statistical analysis. LY, KY, and ZC acquired the data and provided material support. HZ, ZJ, and XX contributed to the conception and design of the study, and analyzed and interpreted the data, supervised the study, provided the project funding, revised the manuscript, and finally approved the version of the manuscript for publication. All authors read and approved the final manuscript.

Funding

This study was supported by the Zhejiang Provincial Natural Science Foundation of China (No. LY20H100003), Medical Health Science and Technology Project of Zhejiang Provincial Health Commission (No. 2019KY453), and Zhejiang Science and Technology Bureau (No. 2012C33126).

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 patients and healthy volunteers who participated in this study.

Supplementary Material

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

FIGURE S1 | Global dynamics of mRNA modifications (calculated as ratio of modified/unmodified nucleotide). The maximal ratio of each modification was set to 1. Points represent individual biological replicates (n = 2) with bars showing the mean.

FIGURE S2 | Top classic mode motifs enriched across the sequences proximal to mRNA m⁵C sites identified from healthy controls and SLE stable and SLE moderate/major active patients, respectively (only E-value < 0.05 is shown here).

FIGURE S3 | Pathway analysis of differentially expressed mRNAs in the CD4⁺ T cells from healthy controls (HCs) and systemic erythematosus lupus patients (SLE). (A) Principal component analysis (PCA) of the mRNA expression profiles in CD4⁺ T cells from HCs and SLE patients. This is another way to visualize sample-to-sample distances. In this ordination method, the data points are projected onto the 2D plane such that they spread out in the two directions that explain most of the differences. The x-axis separates the data points the most. The values of the samples in this direction are referred to as PC1. The y-axis (it must be orthogonal to the first direction) separates the data points the second most. The values of the samples in this direction are referred to as PC2. The percent of the total variance that is contained in the direction is shown on the axis label. Note that these percentages do not add to 100%, because there are more dimensions that contain the remaining variance (although each of these remaining dimensions will explain less than the two that we see). This analysis was performed using R language. (B) Gene ontology analysis of differentially expressed mRNAs. Red indicates molecular function (MF), green indicates cellular component (CC), and blue indicates biological process (BP). (C) Pathway analysis of differentially expressed mRNAs. Red indicates organismal systems, green indicates human diseases, and blue indicates environmental information processing.

TABLE S1 | Clinical features for systemic lupus erythematosus (SLE) patients included in our study.

TABLE S2 | Multiple-reaction monitoring table of mRNA modifications.

TABLE S3 | Dataset of all identified m5C sites in CD4+ T cells of HCs and SLE patients.

TABLE S4 | Genes with the differential m5 C levels in SLE CD4+ T cells compared to healthy controls.

TABLE S5 | The enriched go terms in genes with differential RNA hypermethylation in SLE CD4+ T cells compared with healthy controls according to metascape analysis.

TABLE S6 | The enriched go terms in genes with differential RNA hypomethylation in SLE CD4+ T cells compared with healthy controls according to metascape analysis.

TABLE S7 | Increased m5C levels and up-regulated expression in SLE CD4+ T cells compared with healthy controls.

TABLE S8 | Reactome pathways enriched in genes with increased m5C levels and up-regulated expression in SLE CD4+ T cells compared with healthy controls.

Footnotes

1. ^ http://weblogo.threeplusone.com/
2. ^ http://meme-suite.org/
3. ^ http://metascape.org/gp/index.html
4. ^ https://string-db.org/

References

Abbasi-Moheb, L., Mertel, S., Gonsior, M., Nouri-Vahid, L., Kahrizi, K., Cirak, S., et al. (2012). Mutations in NSUN2 cause autosomal-recessive intellectual disability. Am. J. Hum. Genet. 90, 847–855. doi: 10.1016/j.ajhg.2012.03.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Amort, T., Rieder, D., Wille, A., Khokhlova-Cubberley, D., Riml, C., Trixl, L., et al. (2017). Distinct 5-methylcytosine profiles in poly(A) RNA from mouse embryonic stem cells and brain. Genome Biol. 18:1. doi: 10.1186/s13059-016-1139-1131

PubMed Abstract | CrossRef Full Text | Google Scholar

Auxilien, S., Guérineau, V., Szweykowska-Kulińska, Z., and Golinelli-Pimpaneau, B. (2012). The human tRNA m (5) C methyltransferase Misu is multisite-specific. RNA Biol. 9, 1331–1338. doi: 10.4161/rna.22180

PubMed Abstract | CrossRef Full Text | Google Scholar

Bianco, C., Thompson, L., and Mohr, I. (2019). Repression of eEF2K transcription by NF-κB tunes translation elongation to inflammation and dsDNA-sensing. Proc. Natl. Acad. Scie. U.S.A. 116, 22583–22590. doi: 10.1073/pnas.1909143116

PubMed Abstract | CrossRef Full Text | Google Scholar

Bohnsack, K. E., Hobartner, C., and Bohnsack, M. T. (2019). Eukaryotic 5-methylcytosine (m(5)C) RNA methyltransferases: mechanisms, cellular functions, and links to disease. Genes 10:102. doi: 10.3390/genes10020102

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, X., Li, A., Sun, B. F., Yang, Y., Han, Y. N., Yuan, X., et al. (2019). 5-methylcytosine promotes pathogenesis of bladder cancer through stabilizing mRNAs. Nat Cell Biol. 21, 978–990. doi: 10.1038/s41556-019-0361-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, X. Y., Zhang, J., and Zhu, J. S. (2019). The role of m(6)A RNA methylation in human cancer. Mol. Cancer 18:103. doi: 10.1186/s12943-019-1033-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Coit, P., Yalavarthi, S., Ognenovski, M., Zhao, W., Hasni, S., Wren, J. D., et al. (2015). Epigenome profiling reveals significant DNA demethylation of interferon signature genes in lupus neutrophils. J. Autoimmun. 58, 59–66. doi: 10.1016/j.jaut.2015.01.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Courtney, D. G., Tsai, K., Bogerd, H. P., Kennedy, E. M., Law, B. A., Emery, A., et al. (2019). Epitranscriptomic Addition of m(5)C to HIV-1 transcripts regulates viral gene expression. Cell Host Microbe 26, 217.e6–227e6. doi: 10.1016/j.chom.2019.07.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Cui, X., Liang, Z., Shen, L., Zhang, Q., Bao, S., Geng, Y., et al. (2017). 5-Methylcytosine RNA methylation in Arabidopsis Thaliana. Mol. Plant 10, 1387–1399. doi: 10.1016/j.molp.2017.09.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Docherty, S. J., Davis, O. S., Haworth, C. M., Plomin, R., and Mill, J. (2009). Bisulfite-based epityping on pooled genomic DNA provides an accurate estimate of average group DNA methylation. Epigenetics Chromatin 2:3. doi: 10.1186/1756-8935-2-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Dominissini, D., and Rechavi, G. (2017). 5-methylcytosine mediates nuclear export of mRNA. Cell Res. 27, 717–719. doi: 10.1038/cr.2017.73

PubMed Abstract | CrossRef Full Text | Google Scholar

Flores, J. V., Cordero-Espinoza, L., Oeztuerk-Winder, F., Andersson-Rolf, A., Selmi, T., Blanco, S., et al. (2017). Cytosine-5 RNA methylation regulates neural stem cell differentiation and motility. Stem Cell Rep. 8, 112–124. doi: 10.1016/j.stemcr.2016.11.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, G., Wang, H., Shi, X., Ye, L., Wu, K., Lin, K., et al. (2018). NovelmiRNA-25 inhibits AMPD2 in peripheral blood mononuclear cells of patients with systemic lupus erythematosus and represents a promising novel biomarker. J. Transl. Med. 16:370. doi: 10.1186/s12967-018-1739-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, G., Wang, H., Ye, L., Shi, X., Yan, K., Lin, K., et al. (2019). Hsa_circ_0000479 as a Novel diagnostic biomarker of systemic lupus erythematosus. Front. Immunol. 10:2281. doi: 10.3389/fimmu.2019.02281

PubMed Abstract | CrossRef Full Text | Google Scholar

Han, X., Liu, H., Zhang, Z., Yang, W., Wu, C., Liu, X., et al. (2020). Epitranscriptomic 5-methylcytosine profile in PM2.5-induced mouse pulmonary fibrosis. Genomics Proteomics Bioinformatics S1672-0229, 30020–30026. doi: 10.1016/j.gpb.2019.11.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Ikeda, K., Hayakawa, K., Fujishiro, M., Kawasaki, M., Hirai, T., Tsushima, H., et al. (2017). JAK inhibitor has the amelioration effect in lupus-prone mice: the involvement of IFN signature gene downregulation. BMC Immunol. 18:41. doi: 10.1186/s12865-017-0225-229

PubMed Abstract | CrossRef Full Text | Google Scholar

Lawrence, M., Daujat, S., and Schneider, R. (2016). Lateral thinking: how histone modifications regulate gene expression. Trends Genet. 32, 42–56. doi: 10.1016/j.tig.2015.10.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Lei, W., Luo, Y., Lei, W., Luo, Y., Yan, K., Zhao, S., et al. (2009). Abnormal DNA methylation in CD4+ T cells from patients with systemic lupus erythematosus, systemic sclerosis, and dermatomyositis. Scand. J. Rheumatol. 38, 369–374. doi: 10.1080/03009740902758875

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, L. J., Fan, Y. G., Leng, R. X., Pan, H. F., and Ye, D. Q. (2018). Potential link between m(6)A modification and systemic lupus erythematosus. Mol. Immunol. 93, 55–63. doi: 10.1016/j.molimm.2017.11.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Q., Li, X., Tang, H., Jiang, B., Dou, Y., Gorospe, M., et al. (2017). NSUN2-Mediated m5C methylation and METTL3/METTL14-Mediated m6A methylation cooperatively enhance p21 translation. J. Cell Biochem. 118, 2587–2598. doi: 10.1002/jcb.25957

PubMed Abstract | CrossRef Full Text | Google Scholar

Meroni, P. L., and Penatti, A. E. (2016). Epigenetics and systemic lupus erythematosus: unmet needs. Clin. Rev. Allergy Immunol. 50, 367–376. doi: 10.1007/s12016-015-8497-8494

PubMed Abstract | CrossRef Full Text | Google Scholar

Morel, L. (2017). Immunometabolism in systemic lupus erythematosus. Nat. Rev. Rheumatol. 13, 280–290. doi: 10.1038/nrrheum.2017.43

PubMed Abstract | CrossRef Full Text | Google Scholar

Olarerin-George, A. O., and Jaffrey, S. R. (2017). MetaPlotR: a Perl/R pipeline for plotting metagenes of nucleotide modifications and other transcriptomic sites. Bioinformatics 33, 1563–1564. doi: 10.1093/bioinformatics/btx002

PubMed Abstract | CrossRef Full Text | Google Scholar

Petri, M., Orbai, A. M., Alarcon, G. S., Gordon, C., Merrill, J. T., Fortin, P. R., et al. (2012). Derivation and validation of the systemic lupus international collaborating clinics classification criteria for systemic lupus erythematosus. Arthritis Rheum 64, 2677–2686. doi: 10.1002/art.34473

PubMed Abstract | CrossRef Full Text | Google Scholar

Quinlan, A. R., and Hall, I. M. (2010). BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics 26, 841–842. doi: 10.1093/bioinformatics/btq033

PubMed Abstract | CrossRef Full Text | Google Scholar

Raphael, I., Nalawade, S., Eagar, T. N., and Forsthuber, T. G. (2015). T cell subsets and their signature cytokines in autoimmune and inflammatory diseases. Cytokine 74, 5–17. doi: 10.1016/j.cyto.2014.09.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Rieder, D., Amort, T., Kugler, E., Lusser, A., and Trajanoski, Z. (2016). meRanTK: methylated RNA analysis ToolKit. Bioinformatics 32, 782–785. doi: 10.1093/bioinformatics/btv647

PubMed Abstract | CrossRef Full Text | Google Scholar

Ronnblom, L., and Elkon, K. B. (2010). Cytokines as therapeutic targets in SLE. Nat. Rev. Rheumatol. 6, 339–347. doi: 10.1038/nrrheum.2010.64

PubMed Abstract | CrossRef Full Text | Google Scholar

Roth, S. H., Danan-Gotthold, M., Ben-Izhak, M., Rechavi, G., Cohen, C. J., Louzoun, Y., et al. (2018). Increased RNA editing may provide a source for autoantigens in systemic lupus erythematosus. Cell Rep. 23, 50–57. doi: 10.1016/j.celrep.2018.03.036

PubMed Abstract | CrossRef Full Text | Google Scholar

Roundtree, I. A., Evans, M. E., Pan, T., and He, C. (2017). Dynamic RNA modifications in gene expression regulation. Cell 169, 1187–1200. doi: 10.1016/j.cell.2017.05.045

PubMed Abstract | CrossRef Full Text | Google Scholar

Saferding, V., and Bluml, S. (2019). Innate immunity as the trigger of systemic autoimmune diseases. J. Autoimmun. 110:102382. doi: 10.1016/j.jaut.2019.102382

PubMed Abstract | CrossRef Full Text | Google Scholar

Schaub, A., and Glasmacher, E. (2017). Splicing in immune cells-mechanistic insights and emerging topics. Int. Immunol. 29, 173–181. doi: 10.1093/intimm/dxx026

PubMed Abstract | CrossRef Full Text | Google Scholar

Squires, J. E., Patel, H. R., Nousch, M., Sibbritt, T., Humphreys, D. T., Parker, B. J., et al. (2012). Widespread occurrence of 5-methylcytosine in human coding and non-coding RNA. Nucleic Acids Res. 40, 5023–5033. doi: 10.1093/nar/gks144

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, Z., Xue, S., Xu, H., Hu, X., Chen, S., Yang, Z., et al. (2019). Effects of NSUN2 deficiency on the mRNA 5-methylcytosine modification and gene expression profile in HEK293 cells. Epigenomics 11, 439–453. doi: 10.2217/epi-2018-2169

PubMed Abstract | CrossRef Full Text | Google Scholar

Szklarczyk, D., Morris, J. H., Cook, H., Kuhn, M., Wyder, S., Simonovic, M., et al. (2017). The STRING database in 2017: quality-controlled protein-protein association networks, made broadly accessible. Nucleic Acids Res. 45, D362–D368. doi: 10.1093/nar/gkw937

PubMed Abstract | CrossRef Full Text | Google Scholar

Toperoff, G., Aran, D., Kark, J. D., Rosenberg, M., Dubnikov, T., Nissan, B., et al. (2012). Genome-wide survey reveals predisposing diabetes type 2-related DNA methylation variations in human peripheral blood. Hum. Mol. Genet. 21, 371–383. doi: 10.1093/hmg/ddr472

PubMed Abstract | CrossRef Full Text | Google Scholar

Trixl, L., and Lusser, A. (2019). The dynamic RNA modification 5-methylcytosine and its emerging role as an epitranscriptomic mark. Wiley Interdiscip. Rev. RNA 10:e1510. doi: 10.1002/wrna.1510

PubMed Abstract | CrossRef Full Text | Google Scholar

Tsokos, G. C., Lo, M. S., Costa Reis, P., and Sullivan, K. E. (2016). New insights into the immunopathogenesis of systemic lupus erythematosus. Nat. Rev. Rheumatol. 12, 716–730. doi: 10.1038/nrrheum.2016.186

PubMed Abstract | CrossRef Full Text | Google Scholar

Uribe, A. G., Vila, L. M., McGwin, G. Jr., Sanchez, M. L., Reveille, J. D., and Alarcon, G. S. (2004). The systemic lupus activity measure-revised, the mexican systemic lupus erythematosus disease activity Index (SLEDAI), and a modified SLEDAI-2K are adequate instruments to measure disease activity in systemic lupus erythematosus. J. Rheumatol. 31, 1934–1940.

PubMed Abstract | Google Scholar

Wang, B., Hu, J. J., Yan, C. F., Su, H. H., Ding, J. C., Guo, Y. Y., et al. (2011). Characterization of human cytomegalovirus UL145 and UL136 genes in low-passage clinical isolates from infected Chinese infants. Med. Sci. Monit. 17, CR423–CR431.

PubMed Abstract | Google Scholar

Wang, Y., Zeng, L., Liang, C., Zan, R., Ji, W., Zhang, Z., et al. (2019). Integrated analysis of transcriptome-wide m(6)A methylome of osteosarcoma stem cells enriched by chemotherapy. Epigenomics 11, 1693–1715. doi: 10.2217/epi-2019-2262

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, L., Perrera, V., Saplaoura, E., Apelt, F., Bahin, M., Kramdi, A., et al. (2019). m(5)C methylation guides systemic transport of messenger RNA over graft junctions in plants. Curr. Biol. 29:e2465. doi: 10.1016/j.cub.2019.06.042

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, Y., Wang, L., Han, X., Yang, W. L., Zhang, M., Ma, H. L., et al. (2019). RNA 5-methylcytosine facilitates the maternal-to-zygotic transition by preventing maternal mRNA decay. Mol. Cell 75:e1111. doi: 10.1016/j.molcel.2019.06.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, X., Yang, Y., Sun, B. F., Chen, Y. S., Xu, J. W., Lai, W. Y., et al. (2017). 5-methylcytosine promotes mRNA export - NSUN2 as the methyltransferase and ALYREF as an m(5)C reader. Cell Res. 27, 606–625. doi: 10.1038/cr.2017.55

PubMed Abstract | CrossRef Full Text | Google Scholar

Ye, H., Wang, X., Wang, L., Chu, X., Hu, X., Sun, L., et al. (2019). Full high-throughput sequencing analysis of differences in expression profiles of long noncoding RNAs and their mechanisms of action in systemic lupus erythematosus. Arthritis Res Ther. 21:70. doi: 10.1186/s13075-019-1853-1857

PubMed Abstract | CrossRef Full Text | Google Scholar

Yeung, K. S., Lee, T. L., Mok, M. Y., Mak, C. C. Y., Yang, W., Chong, P. C. Y., et al. (2019). Cell lineage-specific genome-wide DNA methylation analysis of patients with paediatric-onset systemic lupus erythematosus. Epigenetics 14, 341–351. doi: 10.1080/15592294.2019.1585176

PubMed Abstract | CrossRef Full Text | Google Scholar

Yin, Y., Choi, S.-C., Xu, Z., Perry, D. J., Seay, H., Croker, B. P., et al. (2015). Normalization of CD4+ T cell metabolism reverses lupus. Sci. Trans. Med. 7, 274ra218. doi: 10.1126/scitranslmed.aaa0835

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeng, J., Meng, X., Zhou, P., Yin, Z., Xie, Q., Zou, H., et al. (2019). Interferon-alpha exacerbates neuropsychiatric phenotypes in lupus-prone mice. Arthritis Res. Ther. 21:205. doi: 10.1186/s13075-019-1985-1989

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Zhao, M., Sawalha, A. H., Richardson, B., and Lu, Q. (2013). Impaired DNA methylation and its mechanisms in CD4(+)T cells of systemic lupus erythematosus. J. Autoimmun. 41, 92–99. doi: 10.1016/j.jaut.2013.01.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Z., Shi, L., Song, L., Maurer, K., Petri, M. A., and Sullivan, K. E. (2018). Overall downregulation of mRNAs and enrichment of H3K4me3 change near genome-wide association study signals in systemic lupus erythematosus: cell-specific effects. Front. Immunol. 9:497. doi: 10.3389/fimmu.2018.00497

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, M., Liu, S., Luo, S., Wu, H., Tang, M., Cheng, W., et al. (2014). DNA methylation and mRNA and microRNA expression of SLE CD4+ T cells correlate with disease phenotype. J. Autoimmun. 54, 127–136. doi: 10.1016/j.jaut.2014.07.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, M., Sun, Y., Gao, F., Wu, X., Tang, J., Yin, H., et al. (2010). Epigenetics and SLE: RFX1 downregulation causes CD11a and CD70 overexpression by altering epigenetic modifications in lupus CD4+ T cells. J. Autoimmun. 35, 58–69. doi: 10.1016/j.jaut.2010.02.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, M., Wang, J., Liao, W., Li, D., Li, M., Wu, H., et al. (2016). Increased 5-hydroxymethylcytosine in CD4(+) T cells in systemic lupus erythematosus. J. Autoimmun. 69, 64–73. doi: 10.1016/j.jaut.2016.03.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, Y., Zhou, B., Pache, L., Chang, M., Khodabakhshi, A. H., Tanaseichuk, O., et al. (2019). Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat. Commun. 10:1523. doi: 10.1038/s41467-019-09234-9236

PubMed Abstract | CrossRef Full Text | Google Scholar