Transcriptome-Wide High-Throughput m6A Sequencing of Differential m6A Methylation Patterns in the Human Rheumatoid Arthritis Fibroblast-Like Synoviocytes Cell Line MH7A.

INTRODUCTION: N6-methyladenosine (m6A) is the most frequent internal modification in eukaryotic mRNAs and is closely related to the occurrence and development of many diseases, especially tumors. However, the relationship between m6A methylation and rheumatoid arthritis (RA) is still a mystery.
METHODS: Two high-throughput sequencing methods, namely, m6A modified RNA immunoprecipitation sequence (m6A-seq) and RNA sequence (RNA-seq) were performed to identify the differentially expressed m6A methylation in human rheumatoid arthritis fibroblast-like synoviocytes cell line MH7A after stimulation with TNF-α. Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were used to obtain enriched GO terms and significant KEGG pathways. Then, four candidate genes, Wilms tumor 1-associating protein (WTAP), receptor-interacting serine/threonine protein kinase 2 (RIPK2), Janus kinase 3 (JAK3) and tumor necrosis factor receptor SF10A (TNFRSF10A) were selected to further validate the m6A methylation, mRNA and protein expression levels in MH7A cells and synovial tissues of adjuvant arthritis (AA) rats by RT-qPCR and Western blot.
RESULTS: Using m6A-seq, we identified a total of 206 genes with differentially expressed m6A methylation, of which 118 were significantly upregulated and 88 genes were significantly downregulated. Likewise, 1207 differentially mRNA expressed mRNAs were obtained by RNA-seq, of which 793 were upregulated and 414 downregulated. Further joint analysis showed that the m6A methylation and mRNA expression levels of 88 genes changed significantly, of which 30 genes displayed increased m6A methylation and decreased mRNA expression, 57 genes displayed decreased m6A methylation and increased mRNA expression increased, and 1 gene displayed increased m6A methylation and increased mRNA expression. GO and KEGG analyses indicated that these unique genes were mainly enriched in inflammation-related pathways, cell proliferation and apoptosis. In addition, the validations of WTAP, RIPK2, JAK3 and TNFRSF10A were in accordance with the m6A and RNA sequencing results.
CONCLUSION: This study established the transcriptional map of m6A in MH7A cells and revealed the potential relationship between RNA methylation modification and RA related genes. The results suggested that m6A modification was associated with the occurrence and course of RA to some extent.

Introduction

N6-methyladenosine (m6A) is a common means of RNA modification, accounting for 0.1–0.4% of adenosine, and its biological significance has been gradually revealed.1,2 RNA m6A modification is widely found in the mRNAs of eukaryotes, specifically, m6A methylation means that RNA has been modified by methylation at the N6 site of the adenine base, and the modified adenosine is m6A adenine.3,4 Similar to DNA or histone methylation, RNA m6A methylation is also mainly catalyzed by methyltransferases and demethylases.5 This means that the m6A methylation modification is a dynamic reversible process.6 Various m6A methyltransferases, known as “writers,” including methyltransferase-like 3 (METTL3), METTL14 and Wilms tumor 1-associated protein (WTAP), catalyze the m6A methylation of mRNA.5,7 Various m6A demethylases are called “erasers”, including fat mass and obesity-associated protein (FTO) and ALKB family member 5 (ALKBH5), which function is to reduce the modified RNA to the original RNA.8 In addition, m6A binding protein acts as a “reader” and mainly includes the YT521-B homology (YTH) and heterogeneous nuclear ribonucleoprotein (HNRNP) families, such as HNRNPC, YTHDC1, YTHDC2, YTHDF1 and YTHDF2, which recognize m6A-modified RNA and thus regulate mRNA metabolism and function.9,10

For rheumatoid arthritis (RA), one of the most common autoimmune diseases, chronic inflammation, excessive synovial proliferation, joint swelling, and tenderness are the pathological characteristics.11,12 The differential proliferation of fibroblast-like synovial cells (FLSs) is an important factor in the progression of RA.13,14 With the deepening of research and the rapid development of technology, it has been found that there is an association between m6A methylation and RA. Luo et al15 confirmed that the expression of peripheral blood ALKBH5 and FTO was linked with autoantibody production and RA disease activity, which may be used as an indicator of activity and the inflammatory response. However, there are few studies about the effect of m6A methylation on RA FLSs.

In this study, using m6A-seq and RNA-seq, we identified differential m6A methylation and mRNA expression in MH7A cells, a type of human RA FLSs cell line that has been widely used in the study of RA. The sequencing and analysis process are shown in Figure 1.

Figure 1

A schematic diagram of m6A-seq and RNA-seq analyses of MH7A cells. MH7A cells were stimulated with TNF-α, and total RNA was extracted from MH7A cells. Then, the RNA was fragmented, and m6A RNA was isolated as described by the m6A-RNA immunoprecipitation (Me-RIP) assay. The m6A-seq library and RNA-seq library of these mRNAs were constructed, and then sequenced for analysis.

Materials and Methods

Cell Lines and Reagents

MH7A cells were provided by JENNIO Biological Technology (Guangzhou, CHN). MH7A cells were cultured in DMEM medium (Gibco), which 10% FBS (Gemini) and antibiotics (penicillin and streptomycin) were added and maintained at 37°C, 5% CO2 atmosphere. MH7A cells were stimulated with TNF-α (10 ng/mL), which was defined as the TNF-α group, and unstimulated MH7A cells were defined as the control group. After culturing for 24 h, MH7A cells were collected for the subsequent experiment.16

Animal Model

Male Sprague-Dawley rats (8 weeks old, 180–200 g) were purchased from the Anhui Experimental Animal Center. All rats were housed in the animal facility of the First Affiliated Hospital of Anhui University of Chinese Medicine and provided adequate water and food. Animal experiments in this study had been approved by the Animal Ethics Committee of Anhui University of Chinese Medicine, Ethical Approval Number: AHUCM-rats-004. The rats were adaptively fed for 1 week and then subcutaneous injection of complete Freund’s adjuvant into the left hindfoot toe to establish an adjuvant arthritis (AA) rat model. After induction of the immune response for 20 days, all rats were anesthetized by the intraperitoneal injection of 1% sodium pentobarbital (60 mg/kg) and sacrificed by exsanguination of the abdominal aorta. Then, the knee joints were removed, and the synovial tissues were separated. Experimental methods were performed according to the study of Hui Jiang et al.17 The synovial tissues of AA rats were defined as a model group. Conversely, the synovial tissues of normal rats were defined as the control group.

Construction of Library RNA-Seq and m6A-Seq

MH7A cells were digested with trypsin and centrifuged to collect cell precipitates, and the total RNA was extracted by the Total RNA Rapid Extraction Kit (TR205-200, Tianmo, China). The RNA after processing was determined by an Agilent 2100 Bioanalyzer (Agilent Technologies, United States), and the RNA concentration was detected by a NanoDrop one (Thermo Fisher Scientific, United States). VAHTS mRNA Capture Beads was used for capturing poly (A) RNA, which was fragmented by VAHTS 2X Frag/Prime Buffer. Some of the fragments were used for constructing the RNA-seq library with the VAHTS mRNA-seq V2 Library Prep Kit for Illumina, and others were used for the MeRIP-seq library following the instructions of the EpiMark N6-Methyladenosine Enrichment Kit. The obtained libraries were quantified by Qubit 3.0 Fluorometer (Life Technologies) and then sequenced on an Illumina NovaSeq 6000 (Illumina), and the sequencing process was entrusted to Origin-Biotech Inc (Ao-Ji Biotech, Shanghai, China).

Differential Gene Expression Analysis

The quality visualization of RNA-seq reads was accomplished by FastQC (version. 0.11.3).18 Trimming was performed by fastp for known Illumina TruSeq adapter sequences, poor reads, and ribosomal RNA reads.19 The trimmed reads were then mapped to the human reference genome (hg38) using Hisat2.20,21 And p-value < 0.05 and an absolute value of a fold change >2 were used as criteria for screening differentially expressed genes.

Differential m6A Peaks Analysis

Similar to the RNA-seq analysis process, the raw reads were subjected to QC, and trimmed, and m6A-enriched peaks were obtained by mapped to hg38. StringTie (version: 1.3.0) was performed for each gene count from trimmed reads.21 Gene counts were normalized by TMM,22 and FPKM was performed using Perl script.23 EdgeR was used to determine differential gene expression,24 with significance determined by a p-value <0.05 and an absolute value of a fold change >2.25

Functional Enrichment Analysis

Based on the crucial differentially expressed genes obtained by m6A-seq and RNA-seq, Gene Ontology (GO)26 and Kyoto Encyclopedia of Genes and Genomes (KEGG)26 analyses were performed by the clusterProfifiler of R package. Then, the top 30 GO terms and pathways were selected according to the P-value and the degree of enrichment.

PPI Network

The differentially expressed genes (DEGs) were imported into the STRING database, which contains comprehensive information about interactions between proteins;27 thus, the interaction relationship between DEGs was obtained. The protein–protein interaction (PPI) network was constructed based on importing the data into Cytoscape 3.5.1 software, and then, the network was analyzed by Network Analyzer. The DEGs with interactions with combined scores greater than 0.4 were selected to construct a protein–protein interaction network diagram.28

RT-qPCR Validation

RT-qPCR was used to detect candidate gene m6A methylation and mRNA expression levels. Total RNA from MH7A cells and synovial tissues of AA rats was extracted by Mireasy Mini Kit (cat 217004, Qiagen). A Nanodrop nd-2000 spectrophotometer (Thermo Fisher Science, Walsham, Massachusetts, USA) was used to determine the concentration and purity of RNA. Then, cDNA reverse transcription and RT-qPCR reactions were performed using StepOne SYBR Green Real-time PCR Master Mix (ABclonal, RK20404). Reactions proceeded using the following conditions: 95°C for 1 min, followed by 40 cycles of 95°C for 5 s, 60°C for 30 s. Melting curves ranging from 60°C to 95°C were included in each run.

Western Blot

Total protein of MH7A cells and synovial tissues of AA rats were extracted. SDS loading buffer was added to the extracted protein and the mixture was heated in boiling water for 10 min to denature the protein. A 10% polyacrylamide gel was prepared by an SDS-PAGE gel preparation kit (cat P0012A, Beyotime), and then, the protein samples were added to each well for electrophoretic separation. The protein was transferred to a polyvinylidene difluoride membrane and placed in 5% skim milk to seal at room temperature for 2 h. The membranes were incubated with primary antibody against WTAP (cat ab195380, Abcam, 1:2000), RIPK2 (cat bs-3546R, Bioss, 1:1000), JAK3 (cat ab45141, Abcam, 1:1000) and TNFRSF10A (cat bs-0591R, Bioss, 1:1000) overnight at 4°C with slow shaking. Then, the membranes were transferred to horseradish peroxidase-labeled secondary antibodies (1:20000) and incubated continuously at room temperature for 1 h. The protein was developed by a multifunction chemiluminescence imaging system and quantitatively analyzed by ImageJ software.

Statistical Analysis

The experimental data are presented as the mean ± standard deviation (SD). Statistical analysis was performed by using SPSS 23.0 software. Paired Student’s t-tests were used to detect the differences between the two groups. When the p-value < 0.05, the results were considered to be statistically significant.

Results

Transcriptome-Wide m6A-Seq Revealed Global m6A Modification Patterns in MH7A Cells

We compared m6A methylation peaks at each site in MH7A cells after stimulation with TNF-α. The differences and overlaps in m6A methylation between the individuals are shown by the Venn diagram in Figure 2A. We found 13,763 m6A peaks in the control group and 13,621 m6A peaks in the TNF-α group, of which 12,896 m6A peaks were common between the two groups. Compared with the control group, 725 peaks appeared and 867 peaks disappeared in the TNF-α group, indicating that there was a significant difference in the m6A modification pattern of MH7A cells after stimulation with TNF-α.

Figure 2

Transcriptome-wide m6A-Seq revealed global m6A modification patterns in MH7A cells after stimulation with TNF-α. (A) Venn diagram of m6A modified genes in the control group and the TNF-α group. (B) The distribution patterns of m6A peaks in different chromosomes. (C) The number of m6A peaks in per chromosome. (D) The number of m6A peaks in each group. (E), Violin plot of the relative abundance of m6A peaks in each group.

The analysis of the m6A methylation distribution at different chromosome loci found that the chromosomes with the most m6A methylation were chromosome 1 with 1365 m6A methylation peaks, chromosome 2 with 922 m6A methylation peaks and chromosome 19 with 976 m6A methylation peaks, respectively (Figure 2B and C). Then, by comparing the m6A methylation levels, we found an average of 12,658 peaks in the control group and 12,599 peaks in the TNF‐α group, respectively (Figure 2D). The analysis of the degree of enrichment of m6A methylation showed that the average logarithmic fold-enrichment of the control group was 4.91, while that of the TNF-α group was 5.19 (Figure 2E).

Differential m6A Methylation of Genes in MH7A Cells After Stimulation with TNF-α

Using the filtering criteria of fold change > 2 and p-value < 0.05, 206 genes with differential m6A methylation genes in MH7A cells were identified after stimulation with TNF-α, of which 118 m6A hypermethylated genes and 88 m6A hypomethylated genes were identified (Figure 3A and B). The top 10 hypermethylated genes and the top 10 hypomethylated genes are shown in Table 1.

Figure 3

Genes with differential m6A methylation modification genes in MH7A cells after stimulation with TNF-α. (A) Volcano plot representation of microarray data on the differentially expressed m6A methylation genes between the control group and TNF-α group. The blue and red dots to the left and to the right of the two vertical lines indicate more than a 1.5 fold change and represent the differentially expressed m6A methylation genes with statistical significance. Statistical significance was defined as a fold change of 1.5 and a p-value of 0.05 between the control and TNF-α groups. (B) Hierarchical cluster analysis of differentially expressed m6A methylation genes between the control group and TNF-α group. Hierarchical clustering showed that the differentially expressed m6A methylation genes ultimately clustered into two major branches, including hypermethylated genes, which are labeled in red, and hypomethylated genes, which are labeled in green. The darker the color, the more significant the difference. (C) The sector graph shows the proportion of m6A methylation distribution in the indicated regions. (D) The number of genes with differential m6A methylation in per chromosome. (E) PPI of differentially expressed m6A methylation genes. (F) The top 30 enriched GO terms of differentially expressed m6A methylated genes in MH7A cells. (G) The top 30 pathway enrichments of differentially expressed m6A methylated genes in MH7A cells.

Table 1

List of the Top 10 Hypermethylated Genes and the Top 10 Hypomethylated Genes

Genes	Gene Description	Chromosome	Peak Start	Peak End	Fold Change	P-value	Up/Down
TNFAIP3	TNF alpha induced protein 3	6	137,867,214	137,883,312	0.016	0.012	Down
IFNAR2	Interferon alpha and beta receptor subunit 2	21	33,229,901	33,265,675	0.043	0.042	Down
CD83	CD83 molecule	6	14,117,256	14,136,918	0.070	0.049	Down
NUAK2	NUAK family kinase 2	1	205,302,063	205,321,745	0.101	0.002	Down
OAS1	2ʹ-5ʹ-oligoadenylate synthetase 1	12	112,906,783	112,933,222	0.105	0.012	Down
DCAF16	DDB1 and CUL4 associated factor 16	4	17,800,655	17,810,758	0.132	0.047	Down
RIPK2	Receptor interacting serine/threonine kinase 2	8	89,757,806	89,791,064	0.149	0.001	Down
SMPD2	Sphingomyelin phosphodiesterase 2	6	109,440,724	109,443,919	0.149	0.022	Down
RGS20	Regulator of G protein signaling 20	8	53,851,795	53,959,303	0.155	0.030	Down
MCOLN2	Mucolipin 2	1	84,925,583	84,997,113	0.156	0.024	Down
FGF1	Fibroblast growth factor 1	5	142,592,178	142,698,070	34.908	0.046	Up
PRDM16	PR/SET domain 16	1	3,069,168	3,438,621	10.500	0.045	Up
PMF1-BGLAP	PMF1-BGLAP readthrough	1	156,212,982	156,243,332	9.404	0.021	Up
SPCS2	Signal peptidase complex subunit 2	11	74,949,247	74,979,033	9.080	0.049	Up
TRIB2	Tribbles pseudokinase 2	2	12,716,889	12,742,734	7.610	0.008	Up
TGFBR3	Transforming growth factor beta receptor 3	1	91,680,343	91,906,335	7.566	0.033	Up
PIR	Pirin	X	15,384,799	15,493,564	7.201	0.036	Up
TMEM8B	Transmembrane protein 8B	9	35,814,451	35,865,518	6.825	0.033	Up
TGM2	Transglutaminase 2	20	38,127,385	38,166,578	5.890	0.035	Up
DEPTOR	DEP domain containing mTOR interacting protein	8	119,873,717	120,050,918	5.833	0.022	Up

Then, we analyzed the number of m6A methylation modifications in different gene structures. As shown in Figure 3C the differentially expressed genes were mainly distributed in 3ʹ-UTR (31%), with similar percentages in CDS (18%), promoter (18%), 5ʹ-UTR (17%) and terminator (16%). In addition, we also analyzed the distribution of genes with differential m6A methylation genes in chromosomes (Figure 3D). We found that most chromosomes had more than one differentially methylated gene, including 27 on chromosome 1, 17 genes were on chromosome 5 and 16 genes were on chromosomes 6 and 8. Based on interactions with combined scores ≥ 0.4,27 the PPI network analysis constructed interaction networks for genes with differential m6A modification, as shown in Figure 3E.

Simultaneously, the results of GO analysis showed that 854 enriched GO terms in the biological process, cellular component and molecular function categories, particularly in the tumor necrosis factor-mediated signaling pathway, tube formation, regulation of cytokine-mediated signaling pathway, and more. Similarly, KEGG analysis found that 93 pathways were enriched, especially RNA degradation, the NOD-like receptor signaling pathway and apoptosis. The top 30 enriched GO terms and pathways are shown in Figure 3F and G.

Differences in mRNA Expression of MH7A Cells After Stimulation with TNF-α

As shown in Figure 4A and B, a fold change >2 and p-value <0.05 were the screening conditions, and 1207 differentially expressed mRNAs in MH7A cells after stimulation with TNF-α were found, of which 793 were upregulated genes and 414 were downregulated. The top 10 upregulated mRNAs and the top 10 down-regulated mRNAs are shown in Table 2. Based on interactions with combined scores ≥ 0.4,27 the PPI network analysis constructed interaction networks with differentially expressed mRNA genes, as shown in Figure 4C.

Figure 4

Differences in the mRNA expression genes in MH7A cells after stimulation with TNF-α. (A) Volcano plot representation of microarray data on the differentially expressed mRNA genes between the control group and TNF-α group. The blue and red dots to the left and to the right of the two vertical lines indicate more than a 1.5-fold change and represent the differentially expressed genes with statistical significance. Statistical significance was defined as a fold change of 1.5 and a p-value of 0.05 between the control and TNF-α groups. (B) Hierarchical cluster analysis of differentially expressed mRNA genes between the control group and TNF-α group. Hierarchical clustering showed that the differentially expressed genes ultimately clustered into two major branches, including upregulated genes, which are labeled in red, and downregulated genes, which are labeled in green. The darker the color, the more significant the difference. (C) PPI of differentially expressed mRNA genes. (D) The top 30 enriched GO terms of differentially expressed mRNA genes in MH7A cells. (E) The top 30 pathway enrichments of differentially expressed mRNA genes in MH7A cells.

Table 2

List of the Top 10 Up-Regulated mRNAs and the Top 10 Down-Regulated mRNAs

Genes	Gene Description	Fold Change	P-value	Up/Down
MSMP	Microseminoprotein, prostate associated	0.001460	2.31E-07	Down
PNLDC1	PARN like, ribonuclease domain containing 1	0.004044	0.001036	Down
AC018523.2	Novel protein	0.006028	9.21E-05	Down
HOXC10	Homeobox C10	0.009453	6.44E-13	Down
AL096711.2	Novel protein	0.010732	0.000204	Down
TRIM6-TRIM34	TRIM6-TRIM34 readthrough	0.010957	0.005221	Down
AC019257.8	None	0.011574	0.005931	Down
CCKAR	Cholecystokinin A receptor	0.011842	0.001152	Down
NPIPA3	Nuclear pore complex interacting protein family member A3	0.012015	0.006723	Down
SLC26A4	Solute carrier family 26 member 4	0.013341	0.000298	Down
CXCL10	C-X-C motif chemokine ligand 10	155.6325	2.96E-05	Up
ELOVL2	ELOVL fatty acid elongase 2	169.7185	0.000111	Up
OASL	2ʹ-5ʹ-oligoadenylate synthetase like	179.0229	1.29E-06	Up
AL645922.1	Novel complement component 2 (C2) and complement factor B (CFB) protein	197.2662	2.59E-07	Up
CCL20	C-C motif chemokine ligand 20	211.8648	8.73E-07	Up
IL34	Interleukin 34	245.5716	1.26E-08	Up
MMP1	Matrix metallopeptidase 1	304.0153	1.94E-24	Up
CCL5	C-C motif chemokine ligand 5	540.4435	6.78E-61	Up
SERPINB2	Serpin family B member 2	621.2370	1.16E-11	Up
CXCL8	C-X-C motif chemokine ligand 8	1394.952	3.11E-67	Up

Meanwhile, the results of GO analysis showed that 2980 enriched GO terms in the biological process, cellular component and molecular function categories, particularly in response to interferon-alpha, I-kappa B/NF-kappa B complex and CXCR chemokine receptor binding, and more. Similarly, KEGG analysis found that 191 pathways were enriched, particularly the Toll-like receptor signaling pathway, rheumatoid arthritis and cell apoptosis, and more. The top 30 enriched GO terms and pathways are shown in Figure 4D and E.

Conjoint Analysis of m6A-Seq and RNA-Seq Data of MH7A Cells After Stimulation with TNF-α

A conjoint analysis was conducted for m6A-seq and RNA-seq data. We found 88 genes that commonly had differential m6A methylation levels and differentially expressed mRNA levels. Among them, there were 57 genes with m6A hypomethylation and downregulated mRNA expression, 30 genes with m6A hypermethylation and downregulated mRNA expression and 1 gene with m6A hypermethylation and upregulated mRNA expression (Figure 5A and B). Based on interactions with combined scores ≥ 0.4,27 the PPI network analysis constructed interaction networks for these common genes, as shown in Figure 5C.

Figure 5

Joint analysis of m6A methylation and mRNA expression in MH7A cells after stimulation with TNF-α. (A) Four quadrant graph of genes with differential m6A methylation and differentially expressed mRNA levels. (B) Venn diagram of genes with differential m6A methylation and differentially expressed mRNA levels. (C) PPI of genes with differentially expressed m6A methylation and differential expressed mRNA. (D) The top 30 GO enrichments of genes with differentially expressed m6A methylation and differentially expressed mRNA in MH7A cells. (E) The top 30 pathway enrichments of genes with differentially expressed m6A methylation and differentially expressed mRNA in MH7A cells.

Then, we analyzed these common genes by GO and KEGG assays. The GO analysis showed that there were 532 enriched GO terms and these genes were particularly enriched in the tumor necrosis factor-mediated signaling pathway, I-kappaB kinase/NF-kappaB signaling, and cellular response to tumor necrosis factor (Figure 5C and E). Similarly, KEGG analysis found that 59 pathways were enriched, and these genes were closely related to inflammation-related pathways such as the Toll-like receptor signaling pathway, NOD-like receptor signaling pathway, Jak-STAT signaling pathway and cell apoptosis. The top 30 enriched GO terms and pathways are shown in Figure 5D and E.

Correlation Analysis, RT-qPCR and Western Blot Verification of Genes with Differential m6A Modification and Differential mRNA Expression

As shown in Figure 6A and B, joint m6A modification and mRNA expression joint analysis showed changes in m6A methylation abundance and mRNA abundance in WTAP, RIPK2, JAK3 and TNFRSF10A. Compared with the control group, the m6A methylation abundance of WTAP, RIPK2 and JAK3 in MH7A cells after TNF-α stimulation was significantly decreased, while the mRNA abundance were increased. In contrast, the m6A methylation abundance of TNFRSF10A was significantly increased, while the mRNA abundance was significantly decreased. To validate the sequencing results, RT-qPCR and Western blot were applied to detect the mRNA and protein expression levels of WTAP, RIPK2, JAK3 and TNFRSF10A in MH7A cells and synovial tissues of AA rats. The results showed that the mRNA and protein expression trends of these genes were in accordance with the sequencing results (Figure 6C–G).

Figure 6

Correlation analysis, RT-qPCR and Western blot verification of genes with differential m6A modification and differential mRNA expression. (A) Integrative genomics viewer (IGV) plots of m6A methylation abundances and expression abundances in the control group and the TNF-α group. The x-axis represents the genomic position and the y-axis represents the sequence read number. a, the m6A methylation abundances and expression abundances of WTAP. b, the m6A methylation abundances and expression abundances of RIPK2. c, the m6A methylation abundances and expression abundances of JAK3. d, the m6A methylation abundances and expression abundances of TNFRSF10A. (B) The m6A methylation levels were measured by RT-qPCR in the control group and the TNF-α group. a, the m6A methylation level of WTAP. b, the m6A methylation level of RIPK2. c, the m6A methylation level of JAK3. d, the m6A methylation level of TNFRSF10A. (C) The mRNA expression levels in MH7A cells were measured by RT-qPCR. a, the mRNA expression level of WTAP. b, the mRNA expression level of RIPK2. c, the mRNA expression level of JAK3. d, the mRNA expression level of TNFRSF10A. #p < 0.05, ##p < 0.01 compared with the control group. (D) The mRNA expression levels were measured by RT-qPCR in synovial tissues of AA rats. a, the mRNA expression level of WTAP. b, the mRNA expression level of RIPK2. c, the mRNA expression level of JAK3. d, the mRNA expression level of TNFRSF10A. #p < 0.05, ##p < 0.01 compared with the control group. (E) The protein expression levels were measured by Western blot in MH7A cells and synovial tissues of AA rats. a, the protein expression level of WTAP in MH7A cells. b, the protein expression level of WTAP in synovial tissues of AA rats. c, the protein expression level of RIPK2 in MH7A cells. d, the protein expression level of RIPK2 in synovial tissues of AA rats. e, the protein expression level of JAK3 in MH7A cells. f, the protein expression level of JAK3 in synovial tissues of AA rats. g, the protein expression level of TNFRSF10A in MH7A cells. h, the protein expression level of TNFRSF10A in synovial tissues of AA rats. (F), Semiquantitative analysis of protein expression in MH7A cells. a, semiquantitative analysis of WTAP protein. b, semiquantitative analysis of RIPK2 protein. c, semiquantitative analysis of JAK3 protein. d, semiquantitative analysis of TNFRSF10A protein. #p < 0.05, ##p < 0.01 compared with the control group. (G) Semiquantitative analysis of protein expression in synovial tissues of AA rats. a, semiquantitative analysis of WTAP protein. b, semiquantitative analysis of RIPK2 protein. c, semiquantitative analysis of JAK3 protein. d, semiquantitative analysis of TNFRSF10A protein. #p < 0.05, ##p < 0.01 compared with the control group.

Discussion

As an important method of RNA modification, m6A methylation reversibly changes the local structure of RNA and affects the biological function such as the transcription and translation of RNA, thus affecting the expression of downstream target genes.29,30 Therefore, it has an intimate connection to the occurrence and development of many diseases. In recent years, an increasing number of studies on m6A and diseases have been carried out.31,32 However, the transcriptome-wide distributions of m6A modifications in MH7A cells and the specific function of m6A methylation differences in the pathogenesis of RA are not clear.

The m6A and RNA sequence approaches were utilized to construct the m6A methylation library and RNA library to investigate differentially m6A methylation and gene expression in MH7A cells. Our results showed that there were 206 differentially expressed m6A methylation genes and 1207 differentially mRNAs in MH7A cells after stimulation with TNF-α. By analyzing the transcriptional region with the m6A peaks, the differentially expressed m6A methylated genes were mainly concentrated in the 3ʹ-UTR, which was quite consistent with the previous findings of Dominissini et al.33 On the other hand, according to the joint analysis of genes with differentially expressed m6A methylation and mRNA, there were significant differences between the m6A methylation and mRNA expression levels of 88 genes. We suspected that these potential genes, through m6A methylation, are involved in the occurrence and development of RA.

We conducted GO and KEGG analyses of differentially expressed m6A methylated genes. Biological processes and pathways related to inflammation, cell characteristics and cell fate, especially the NOD-like receptor single pathway, NF-κB pathway and cell apoptosis were enriched. NOD‐like receptor family, pyrin domain containing 3 (NLRP3) gene polymorphisms have a direct relationship with inflammation in the course of RA. Guo et al34 found that NLRP3 inflammatory bodies were highly activated in the synovium of RA patients and CIA rats, and the inhibition of NLRP3 could significantly improve the pathological symptoms. The NF-κB pathway plays an important role in the regulation of immune and inflammatory responses. Roman JA et al35 found that NF-κB is differentially activated in RA synovium, and the NF-κB pathway regulates the expression of many inflammatory factors, such as cytokines and cell adhesion molecules, in the synovium and synovial fluid of RA patients. Apoptosis is an active programmed death of cells that maintains the stability of the body.36 Both apoptosis and cell proliferation are precisely regulated by genes to ensure the dynamic balance of the number of cells in the body.37 It has been found that the damage to articular cartilage and bone tissue in the process of RA is closely related to the excessive proliferation of synovial cells and insufficient apoptosis.38,39

Then, according to biological function, we selected WTAP, RIPK2, JAK3 and TNFRSF10A for further analysis. The sequencing, RT-qPCR and Western blot results in MH7A cells after stimulation with TNF-α all clearly indicated that the m6A methylation levels of WTAP, RIPK2 and JAK3 were significantly reduced, while the mRNA and protein expression levels were significantly increased. Interestingly, the m6A methylation as well as mRNA and protein expression levels of TNFRSF10A presented the completely opposite trends compared with the above three genes. To further validate the above experimental results, we established and detected the mRNA and protein expression levels of the above genes in synovial tissues of AA model rats. The expression levels of WTAP, RIPK2, JAK3 and TNFRSF10A showed the same trends between the cellular level and animal level.

WTAP, as an important component of the m6A methyltransferase complex, does not have N6-methyladenine methyltransferase activity but is required for efficient RNA methylation in vivo and for the localization of METTL3 and METTL14 to nuclear speckles.40,41 WTAP has been found to be strongly related to the development of many diseases, including autoimmune diseases.42,43 Tong et al44 found that mice with the specific deletion of m6A methyltransferases in regulatory T cells would develop severe autoimmune diseases such as enlarged lymph nodes and an enlarged spleen, indicating that WTAP and other methyltransferases are closely related to T cell function and involved in the occurrence of immune diseases. TNFRSF10A, also known as TRAIL-R1 or DR4, is a member of the tumor necrosis factor receptor family, which can regulate cell function and so on.45,46 TNFRSF10A can not only induce apoptosis but also inhibit inflammation, and its ligands can activate the NF-κB pathway and promote cell proliferation and migration.47,48 Li et al49 found that TNFRSF10A was involved in apoptosis induced by endoplasmic reticulum stress, and TNFRSF10A knockout reduced the expression level of TNFRSF10B, thereby protecting cells from apoptosis. RIPK2 is the downstream target of the pattern recognition receptors NOD1 and NOD2 and is closely related to inflammatory diseases.50,51 Usluoglu N et al52 found that RIPK2 knockout could reduce the phosphorylation of p38MAPK, ERK and Ikappa B, leading to a reduction in the inflammatory factor IL-12. In addition, Krieg A et al53 also found that the inhibition of endogenous RIPK2 significantly reduced cell apoptosis. JAK3 belongs to the intracellular nonreceptor tyrosine kinase family, which mediates the signal produced by cytokines and is transmitted through the JAK-STAT signal pathway.54 Some studies have reported that JAK3 plays an extensive role in the pathogenesis of many autoimmune diseases.55,56 Byung-Hak Kim et al57 found that the upregulation of JAK3/STATs is closely associated with arthritic inflammation, and JAK3 antagonists inhibit JAK3 activity to reduce inflammation in the body. From our verification results, there were divergences in m6A methylation and mRNA expression levels in these genes. This suggests that the pathogenesis of RA is related to changes in the RNA methylation levels of these genes.

Conclusion

In summary, our study revealed the differential expression of m6A methylation and mRNA genes in MH7A cells, provided the m6A transcriptome map of MH7A cells, and revealed the potential relationship between RNA methylation modification and RA related genes. Nevertheless, the specific and mechanism of m6A methylation in RA remains to be further studied in the future.