Potential Target Analysis of Triptolide Based on Transcriptome-Wide m6A Methylome in Rheumatoid Arthritis

Triptolide (TP), a major active component of the herb Tripterygium wilfordii Hook F (TwHF), has been shown to exert therapeutic potential against rheumatoid arthritis (RA). However, its molecular mechanism of action has not been fully elucidated. This study aimed to analyze the potential target of TP based on the discovery of differentially methylated and expressed genes (DMEGs) in RA using methylated RNA immunoprecipitation sequencing (MeRIP-seq) and RNA sequencing (RNA-seq). Five RA samples and ten control samples were obtained from China-Japan Friendship Hospital. The various levels of m⁶A methylation and genes expressed in the RA and control groups were compared by MeRIP-seq and RNA-seq. Bioinformatics explorations were also performed to explore the enriched biological roles and paths of the differentially expressed m⁶A methylation and genes. Molecular networks between TP target proteins and DMEGs were performed using Ingenuity Pathway Analysis (IPA) software. Potential target of TP was determined with Gene Expression Omnibus (GEO) database mining, molecular docking, and in vitro experiment validation. In total, 583 dysregulated m⁶A peaks, of which 295 were greatly upregulated and 288 were greatly downregulated, were identified. Similarly, 1,570 differentially expressed genes were identified by RNA-seq, including 539 upregulated and 1,031 downregulated genes. According to the deeper joint exploration, the m⁶A methylation and mRNA expression degrees of 35 genes varied greatly. Molecular networks between TP target proteins and DMEGs were constructed, and the results revealed that tubulin beta-2A chain (TUBB2A), insulin-like growth factor 2 mRNA-binding protein 3 (IGF2BP3), cytoplasmic dynein 1 intermediate chain 1 (DYNC1I1), and FOS-like 1 (FOSL1) were the most relevant genes that correlated with the target proteins of TP. The results of the GEO database showed that the gene expression of IGF2BP3 was increased in RA synovial tissue and consistent with the trend of our sequencing results of RA PBMCs. Molecular docking and in vitro experiment suggested that TP and IGF2BP3 had a high binding affinity and TP could decrease the mRNA expression of IGF2BP3 in PBMCs and MH7A.This research established a transcriptional map of m⁶A in RA PBMCs and displayed the hidden association between RNA methylation alterations and associated genes in RA. IGF2BP3 might be a potential therapeutic target of TP during RA treatment.

Introduction

As an autoimmune disease, rheumatoid arthritis (RA) has the characteristics such as chronic inflammatory variations in synovial tissue, bone, and joint cartilage, which lead to the substantial destruction of affected joints (Cross et al., 2014; Scherer et al., 2020). The etiology of RA has not been explained completely, but the pathogenesis of RA relies on the complicated interplay between genetic and environmental triggers (Doody et al., 2017). Modern genetic technologies used in large well-characterized clinical studies have determined the participation of genetic elements in the development of RA and have advanced our understanding of the genetics of the disease (Smolen et al., 2016). However, there is no genetic heterogeneity to interpret all the characteristics of RA. Previous studies have confirmed that RA is greatly dependent on age (Sparks et al., 2016), microbiome composition (Zhang et al., 2015), certain early life factors, such as birth weight (Jacobsson et al., 2003) and breastfeeding (Karlson et al., 2004), and environmental elements, including smoking (Pedersen et al., 2007; Sparks et al., 2016). An increasing amount of evidence suggests that epigenetics contributes to the pathogenesis of RA by bridging the environmental and genetic effects. Thus, epigenetic factors are another promising area for investigating the etiology of RA.

As an epigenetic regulatory mechanism, N⁶-methyladenosine (m⁶A) is the most abundant internal alteration of mRNA in eukaryotes. This alternation is reversible and regulated by “writers,” “erasers,” and “readers” (Fu et al., 2014). Recently, according to accumulating evidence, it has been reported that the m⁶A modification is indispensable for the biogenesis and roles of RNA and plays an important role in various disease pathological processes, including obesity (Loos and Yeo, 2014) and systemic lupus erythematosus (SLE) (Li et al., 2018; Luo et al. (2020b)). Many technologies for quantifying m⁶A-altered transcripts have been developed to determine the phenotypes induced by deleting proteins related to m⁶A modification and ways to affect m⁶A-altered transcripts at the molecular level, such as the transcriptome-wide scale of m⁶A modifications with methylated RNA immunoprecipitation sequencing (MeRIP-seq). Luo et al. (2020a) confirmed that, in RA, decreased alkB homolog 5 (ALKBH5), fat mass and obesity-associated protein (FTO), and YTH domain-containing family member 2 (YTHDF2) expression in peripheral blood mononuclear cells (PBMCs) is a risk factor for RA, and other research has shown that methyltransferase like 3 (METTL3) expression is greatly increased in the PBMCs of RA patients (Wang et al., 2019). Nevertheless, the transcriptome-wide m⁶A methylome in RA has not been fully characterized, and more in-depth research is rare.

Tripterygium wilfordii Hook F (TwHF), a traditional Chinese medicine, has been used for centuries to treat autoimmune diseases including RA in China. Previous studies have confirmed the effectiveness of TwHF in the treatment of RA (Lv et al., 2015; Zhou et al., 2018). Extracts of TwHF have also been tested in RA patients with good efficacy noted in the West (Tao et al., 2002; Goldbach-Mansky et al., 2009). As the main active ingredient in TwHF, increasing experimental evidence has verified triptolide’s (TP) anti-RA effect (Fan et al., 2016; Wang et al., 2018), and it has been considered as a promising anti-RA drug (Han et al., 2012). However, due to the complexity of disease, the underlying molecular mechanism of TP is not yet clear.

In this study, MeRIP-seq and RNA sequencing (RNA-seq) were performed, and differentially methylated peaks within mRNAs were found in RA and control samples. Meanwhile, differentially expressed mRNAs were identified using RNA-seq. In addition, further analysis of the combined MeRIP-seq and RNA-seq information was also performed for the identification of differentially methylated and expressed genes (DMEGs). Furthermore, the molecular network between TP target proteins and DMEGs was constructed using Ingenuity Pathway Analysis (IPA) software. A new potential target of TP was determined with Gene Expression Omnibus (GEO) database mining, molecular docking, and in vitro experiment validation.

Materials and Methods

Patients and Controls

Five patients suffering from RA were obtained from China-Japan Friendship Hospital. All patients had been diagnosed with RA according to the 2010 American College of Rheumatology/European League Against Rheumatism (ACR/EULAR) standards (Aletaha et al., 2010). In this study, the EULAR disease activity score (DAS28) (Prevoo et al., 1995) of five patients varied from 3.16 to 5.70. The patients were 52.00 (±9.08) years old on average. Ten subjects without RA or other autoimmune diseases were enrolled, in which five control samples (49.00 ± 7.62 years old on average) were used for MeRIP-seq, and the other five control samples (48.40 ± 5.46 years old on average) were used for the other experiments. The ethical committee of China-Japan Friendship Hospital approved the research with the ethical approval number 2020-133-K86.

PBMC Preparation and Total RNA Extraction

PBMCs from each subject were separated using Lymphoprep™ with a density of 1.077 ± 0.001 g/ml, according to the manufacturer’s instructions (Stemcell Technologies, Canada). Total RNA in PBMCs was extracted, purified with TRIzol (Invitrogen, Carlsbad, United States) according to the manufacturer’s instructions, and stored at −80°C.

MeRIP-Seq Assays

A NanoDrop ND-1000 (NanoDrop, Wilmington, United States) was used to quantify the total RNA amount and purity of every sample. A Bioanalyzer 2100 (Agilent, CA, United States) with an RIN > 7.0 was employed to assess the RNA integrity, followed by confirmation by electrophoresis in a denaturing agarose gel. Ribosomal RNA was depleted from approximately more than 25 μg of total RNA from a specific adipose type, based on the instructions of the Epicenter Ribo-Zero Gold Kit (Illumina, San Diego, United States). After purification, the 7-min fragmentation of ribosomal-depleted RNA into small pieces was performed with a magnesium RNA fragmentation module (NEB, United States) at 86°C. Next, the 2 h incubation of the cleaved RNA fragments was performed at 4°C with the provided m⁶A antibody (Synaptic Systems, Germany) in IP buffer (50 mM Tris-HCl, 750 mM NaCl, and 0.5% Igepal CA-630). Next, SuperScript™ II reverse transcriptase (Invitrogen) was used to reverse-transcribe the IP RNA for the creation of cDNA, which was then used to synthesize U-labeled second-stranded DNAs with E. coli DNA polymerase I (NEB, United States), RNase H (NEB, United States), and dUTP solution (Thermo Fisher, United States). An A-base was placed on the blunt ends of every strand for ligation to the indexed adapters. There was a T-base overhang for ligating the adapter to the A-tailed-fragmented DNA on each adapter. During the ligation of single- or dual-index adapters to the fragments, size was selected with AMPureXP beads. After treatment of the U-labeled second-stranded DNAs with the heat-labile UDG enzyme (NEB, United States), amplification of the ligated products was performed with PCR with 3-min initial denaturation at 95°C, eight cycles of denaturation at 98°C for 15 s, 15 s annealing at 60°C, 30 s extension at 72°C, and 5 min final extension at 72°C. The final cDNA library showed an average insert size of 300 ± 50 bp. Eventually, 2 × 150 bp paired-end sequencing (PE150) was performed on an Illumina Novaseq™ 6000 (LC-Bio Technology Co., Ltd., Hangzhou, China), according to the vendor’s suggested instructions (Dominissini et al., 2013).

Data Analysis

The reads with adapter contamination, low quality bases, and uncalled bases with default coefficients were removed with fastp software (Chen et al., 2018). Then, fastp was used to verify the sequence quality of the IP and input samples. Reads were mapped to the reference genome Homo sapiens (Version: v101) with HISAT2 (Kim et al., 2015). The mapped reads of IP and input libraries were assessed by R package exomePeak (Meng et al., 2014) for identifying m⁶A peaks with a bed or bigwig format that can be changed for visualization on IGV software. De novo and known motif finding used MEME and HOMER according to localization of the motif in terms of peak summit (Bailey et al., 2009; Heinz et al., 2010). R package ChIPseeker was employed to annotate the called peaks by intersection with the gene architecture (Yu et al., 2015). Then, the expression level was determined with StringTie for all mRNAs from input libraries by using the calculation of FPKM [total exon fragments/mapped reads (millions) × exon length (kB)] (Pertea et al., 2015). A log2 (fold change) > 1 or log2 (fold change) <-1 and p value <0.05 by R package edgeR were used to select the differentially expressed mRNAs (Robinson et al., 2010).

Gene Ontology and Kyoto Encyclopedia of Genes and Genomes Analysis

GO and KEGG pathway analyses were conducted on m⁶A peaks and mRNAs. GO exploration was performed, which included biological processes, cellular components, and molecular functions. On the basis of the KEGG database, the potential functions were analyzed by pathway analysis.

Venn Analysis

Overlapped dysregulated mRNAs between the RA and control (Con) groups detected by RNA-Seq were characterized by Venn analysis. Various colors in the pie diagram show the number of mRNAs with or without overlapping.

Analysis of Molecular Networks Between TP and DMEGs

By using “triptolide” as a key word, we searched for the potential human target proteins of TP in the PubChem platform (http://pubchem.ncbi.hlm.nih.gov). Then, the human target proteins and DMEGs were imported into the Ingenuity Pathway Analysis (IPA) platform, and molecular networks were built. In the molecular networks, nodes represent molecules, and an edge (line) shows the biological association between two nodes. At least, one reference from a textbook, literature, or canonical data saved in the IPKB supported all edges. Nodes with various shapes represent the functional class of the gene goods.

Microarray Data

The gene expression data of GSE55235, GSE55457, and GSE55584 based on the GPL96 platform, and of GSE36700 based on the GPL570 platform were downloaded from the GEO database. The datasets contained varying number of synovial tissue samples and were selected for subsequent analysis.

Molecular Docking

The TP’s structure was obtained from the PubChem platform (Kim et al., 2021). The PDB file for the insulin-like growth factor 2 mRNA-binding protein 3 (IGF2BP3) protein’s crystal structure was downloaded from the RCSB Protein Data Bank (PDB) database (Burley et al., 2021). Preparing and preprocessing of TP and IGF2BP3 protein were performed using AtuoDock Tools (Tanchuk et al., 2016) according to the tutorial and manual (http://vina.scripps.edu/manual.html). TP was docked with IGF2BP3 protein by AutoDock Vina (Trott and Olson, 2010). Generally, the docking score less than -5 kcal/mol was considered to have a strong binding affinity of a compound-target pair. PyMol was used for visualizing the results of AutoDock Vina (Alexander et al., 2011).

Cell Culture

The isolated human PBMCs of the control group were cultured in RPMI 1640 medium (GIBCO, Gaithersburg, United States) and fibroblast-like synoviocytes (FLSs) cell line MH7A (JENNIO Biological Technology, Guangzhou, China) were cultured in DMEM medium (GIBCO) with 10% fetal bovine serum (GIBCO) and 1% penicillin–streptomycin in an incubator with 5% CO₂ at 37°C.

Treatment of PBMCs and MH7A With TP

PBMCs or MH7A were placed into 6-well plates at a concentration of 1.0 × 10⁶ cells/ml and 1 × 10⁵ cells/ml, respectively, followed by 2 h of pretreatment with or without 100 ng/ml of lipopolysaccharides (LPS) (Sigma, St. Louis, MO, United States) or 20 ng/ml recombinant human tumor necrosis factor (TNF)-α (PeproTech, New Jersey, United States). Then, PBMCs were incubated with 6.25 nM or 12.5 nM TP (purity: 99.8%, National Institutes for Food and Drug Control, Beijing, China) for another 24-h, MH7A was incubated with 12.5 nM or 25 nM TP for another 24-h. Finally, the cells were harvested for real-time PCR analysis.

Quantitative Real-Time PCR

The PrimeScriptTM RT Reagent Kit with gDNA Eraser (Takara, Dalian, China) was used to reverse-transcribe the extracted RNA into cDNA using the instructions provided by the manufacturer. Quantitative real-time PCR analysis was performed with TB GreenTM Premix Ex TaqTM (Takara, Dalian, China) and the QuantStudio 5 Real-Time PCR System (Applied Biosystems, Bedford, MA, United States). The relative RNA expression was quantified with the 2^−ΔΔCt method and β-actin was used as an endogenous control. Supplementary Table S2 shows the primers used.

Statistical Analysis

The data are shown as the mean ± SEM. The statistical significance for differences between two groups was determined with two-tailed, unpaired Student’s t-test or one-way analysis of variance (ANOVA) according to Tukey’s post hoc analysis. All explorations were performed with GraphPad Prism 8.0 (La Jolla, CA, United States). A p < 0.05 was regarded as statistically significant.

Results

Overview of the m⁶A Methylation Map in Human PBMCs

To reveal the global m⁶A modification patterns in RA PBMCs, MeRIP-seq analysis was performed. According to Figure 1A, the m⁶A consensus motif of GGACU was presented with highly enriched in human PBMCs, which conforms to the reported consensus motif RRACH (R = G or A; H = A, C or U) (Dominissini et al., 2012). In total, 41,026 and 41,184 m⁶A peaks from 22,080 to 22,206 m⁶A-modified transcripts in Con and RA PBMCs were identified, respectively (Figures 1B,C). In contrast to the unique m⁶A peaks of the RA group (20,589) and the Con group (20,435), there were 20,595 and 20,591 common m⁶A modification peaks in the RA and Con groups, respectively (Figure 1B). In addition, 10,839 and 10,832 common m⁶A-modified genes were identified in Con and RA PBMCs, respectively, and 11,241 and 11,374 unique m⁶A-modified genes were found in Con and RA PBMCs, respectively (Figure 1C). According to the exploration of the m⁶A methylation distribution at various chromosome loci, the chromosomes with the greatest m⁶A methylation were chromosome 1 with 7,998 m⁶A methylation peaks, chromosome 17 with 5,244 m⁶A methylation peaks and chromosome 19 with 6,527 m⁶A methylation peaks (Figure 1D). The m⁶A distribution patterns within the entire transcripts were also analyzed, and the results showed a similar pattern of total m⁶A distribution in the Con and RA groups. For the dysregulated peaks, 40.71% were in the 3′-UTR area, 20.58% were in the 5′-UTR area, 12.69% were in the first exon, and 26.01% were in the residual exons (Figures 1E,F). Furthermore, the m⁶A peak abundance of common m⁶A-altered genes between PBMCs of the Con and RA patients was compared. We discovered 295 hypermethylated m⁶A peaks and 288 hypomethylated m⁶A peaks in the RA PBMCs relative to the Con PBMCs (log2 |fold change| > 1 and p value <0.05; Figure 1G). Table 1 shows the top 20 differentially methylated genes (DMGs), and Supplementary Table S3 displays the list of all DMGs.

FIGURE 1

FIGURE 1. Characteristics of m⁶A methylation in RA PBMCs. (A) Recognition of the top consensus motif from m⁶A peaks determined in PBMCs from the RA and Con groups. (B) Number of m⁶A peaks recognized by m⁶A-seq in the RA and Con groups. (C) Summary of the m⁶A-altered genes identified using m⁶A-seq. (D) Distribution patterns of m⁶A peaks in various chromosomes. (E) m⁶A modification distribution in different gene contexts. (F) Accumulation of m⁶A peaks across transcripts fallen into5′-UTR, CDS, and 3′-UTR. (G) Recognition of the abundance of hypermethylated and hypomethylated m⁶A peaks in the RA PBMCs relative to that of the Con group. m⁶A, N⁶-methyladenosine; RA, rheumatoid arthritis; PBMCs, peripheral blood mononuclear cells; and Con, control group.

TABLE 1

TABLE 1. Top 20 differentially methylated m⁶A peaks (RA/Con).

Differentially Methylated Genes Participated in Important Biological Functions and Pathways

To study the biological importance of m⁶A modification in RA, GO enrichment and KEGG pathway analyses were performed. GO analysis results showed that both hypermethylated and hypomethylated genes in the PBMCs participated in the following functions and pathways: “regulation of transcription, DNA-templated,” “signal transduction,” and “cell differentiation” (ontology: biological process); “cytoplasm,” “membrane,” and “nucleus” (ontology: cellular component); and “protein binding,” “nucleotide binding,” “metal ion binding,” and “DNA binding” (ontology: molecular function) (Figures 2A,D). In KEGG subclass analyses, hypermethylated and hypomethylated genes were involved in the following functions and pathways: “cell motility,” “cellular community-eukaryotes,” and “cell growth and death” (cellular processes); and “signal transduction” and “signaling molecules and interaction” (environmental information processing). In addition, hypermethylated and hypomethylated genes were associated with “immune disease” (human diseases) and “immune system” (organismal systems) (Figures 2B,E). For KEGG pathway enrichment, we found that the hypermethylated genes were greatly associated with the “calcium signaling pathway” and the “B-cell receptor signaling pathway” (Figure 2C). However, the hypomethylated genes were significantly related to “sphingolipid metabolism” and “ECM-receptor interaction” (Figure 2F).

FIGURE 2

FIGURE 2. Biological function and pathway analysis of differentially methylated genes between the RA and Con groups. (A) GO enrichment exploration of hypermethylated peaks. (B) KEGG enrichment subclass analysis of hypermethylated peaks. (C) Top 20 KEGG pathways exploration of hypermethylated peaks. (D) GO enrichment exploration of hypomethylated peaks. (E) KEGG enrichment subclass analysis of hypomethylated peaks. (F) Top 20 KEGG pathways in the hypomethylated peaks. RA, rheumatoid arthritis; Con, control group; GO, Gene Ontology; and KEGG, Kyoto Encyclopedia of Genes and Genomes.

Overview of the Transcriptome Profiles

To clarify the altered gene expression and the corresponding signaling pathways in RA PBMCs, RNA-seq was performed to identify the differentially expressed genes (DEGs). In Figure 3A, the differences and overlaps of genes in the RA and Con groups are shown using a Venn diagram. We found that 2,252 genes in the Con group, 2,477 genes in the RA group, and 29,979 genes were common between the two groups. Next, applying the standard of p value <0.05 and fold change (FC) ≥ 2, 1,570 genes were shown to be greatly dysregulated compared with the Con group, including 539 upregulated genes and 1,031 downregulated genes (Figure 3B). Figure 3C shows the hierarchical clustering of the DEGs. The full list of upregulated and downregulated genes is contained in Supplementary Table S4.

FIGURE 3

FIGURE 3. RNA-seq data analysis and the differentially expressed genes between the RA and Con groups. (A) Venn diagram of RNA-seq data displaying the common and unique genes between the RA and Con groups. (B) Volcano plots showing the differentially expressed genes in RA PBMCs compared with those in the Con group. (C) Heatmap plots displaying the differentially expressed genes in RA PBMCs vs. those in the Con group. RA, rheumatoid arthritis; Con, control group; and PBMCs, peripheral blood mononuclear cells.

Differentially Expressed Genes Participated in Important Biological Functions and Pathways

To analyze the biological functions and signaling pathways involved in DEGs, GO and KEGG analyses were performed. The outcomes of the GO analysis suggested the enrichment of these genes in the following functions and pathways: “signal transduction” and “cell adhesion” (ontology: biological process); “membrane,” “cytoplasm,” and “integral component of membrane” (ontology: cellular component); and “protein binding” and “DNA binding” (ontology: molecular function). Upregulated and downregulated genes were also involved in “neutrophil degranulation” and “regulation of transcription, DNA-templated” (ontology: biological process), respectively (Figures 4A,D).

FIGURE 4

FIGURE 4. Biological function and pathway analysis of differentially regulated genes between the RA and Con groups. (A) GO enrichment analysis of the upregulated genes with great diversity. (B) KEGG enrichment subclass exploration of the upregulated genes with significant differences. (C) Top 20 KEGG pathway enrichments of the upregulated genes with great diversity. (D) GO enrichment exploration of the significantly downregulated genes. (E) KEGG enrichment subclass exploration of the highly downregulated genes. (F) Top 20 KEGG pathway enrichments of the significantly downregulated genes. RA, rheumatoid arthritis; Con, control group; GO, Gene Ontology; and KEGG, Kyoto Encyclopedia of Genes and Genomes.

In KEGG subclass analyses, upregulated and downregulated genes were involved in the following functions and pathways: “cell motility,” “cellular community-eukaryotes,” “cell growth and death,” and “transport and catabolism” (cellular processes); and “membrane transport,” “signaling molecules and interaction,” and “signal transduction” (environmental information processing). In genetic information processing, upregulated and downregulated genes were involved in “folding, sorting, and degradation.” In metabolism, “glycan biosynthesis and metabolism” was associated with the upregulated and downregulated genes. In addition, upregulated and downregulated genes were associated with “immune disease” (human diseases) and “immune system” (organismal systems) (Figures 4B,E). According to KEGG pathway exploration, the upregulated genes were strongly related to “complement and coagulation cascades,” “Jak-STAT signaling pathway” and “IL-17 signaling pathway” (Figure 4C). However, the downregulated genes were significantly related to “natural killer cell mediated cytotoxicity,” “antigen progressing and representation,” and “protein digestion and absorption” (Figure 4F).

Conjoint Analyses of MeRIP-Seq and RNA-Seq Data

Furthermore, by the conjoint analysis of MeRIP-seq and RNA-seq data, 35 DMEGs changed significantly, including 13 hypermethylated and upregulated genes (hyper-up), 12 hypermethylated and downregulated genes (hyper-down), five hypomethylated and upregulated genes (hypo-up), and five hypomethylated and downregulated genes (hypo-down) (Figure 5A). Table 2 shows the list of 35 DMEGs. The GO analysis revealed that the DMEGs were associated with the following functions and pathways: “protein transport,” “intracellular protein transport,” and “multicellular organism development,” (biological process); “membrane,” “nucleus,” and “cytosol” (cellular component); and “protein binding,” “nucleotide binding,” and “GTP binding” (molecular function) (Figure 5B). In KEGG subclass analysis, the DMEGs were involved in “cellular community” and “transport and catabolism” (cellular processes) and “signal transduction” (environment information processing). In human diseases, DMEGs were related to “immune disease” (Figure 5C). The KEGG pathway analysis displayed the primary enrichment of these DMEGs in “lysosome,” “phagosome,” and “component and coagulation cascades” (Figure 5D).

FIGURE 5

FIGURE 5. Conjoint analysis of the MeRIP-seq and RNA-seq data. (A) Distribution of genes with a great variation in both the m⁶A modification and mRNA levels in the PBMCs of RA relative to that of Con. (B) GO enrichment exploration of genes with great variation in both the m⁶A modification and mRNA levels. (C) KEGG subclass of genes with great variation in both the m⁶A modification and mRNA levels. (D) KEGG pathway enrichments of differentially expressed genes at both the m⁶A modification and mRNA levels. PBMCs, peripheral blood mononuclear cells; RA, rheumatoid arthritis; Con, control group; GO, Gene Ontology; m⁶A, N⁶-methyladenosine; and KEGG, Kyoto Encyclopedia of Genes and Genomes.

TABLE 2

TABLE 2. Differentially expressed genes with differentially methylated peaks.

Potential Target Analysis of TP and Hub Gene Validation

To construct the molecular network between TP target proteins and candidate DMEGs, we first searched and found 169 target proteins of TP (Supplementary Table S5) in the PubChem platform. Then, the molecular network of TP target proteins and DMEGs was built. The result revealed that the tubulin beta-2A chain (TUBB2A), IGF2BP3, cytoplasmic dynein 1 intermediate chain 1 (DYNC1I1), and FOS-like 1 (FOSL1) were the most relevant genes that correlated with the target proteins of TP (Figure 6A). In addition, the gene expression of IGF2BP3 in RA synovial tissue of different GEO datasets was increased (Figures 6B–E) and consistent with the trend of our sequencing results of RA PBMCs. So, we chose IGF2BP3 as a candidate gene for further analysis. To verify the binding between IGF2BP3 and TP, we performed molecular docking. The crystal structure of IGF2BP3 (PDB, ID: 6fq1) was obtained from RCSB Protein Data Bank and the chemical structure of TP was downloaded from the PubChem platform. The docking scores between IGF2BP3 and TP is -8.6 kcal/mol (Figure 6F), suggesting that TP and IGF2BP3 have a higher binding affinity. The m⁶A peak distributions of IGF2BP3 between the Con and RA groups are shown in Supplementary Figure S1A. Then, the mRNA expression of IGF2BP3 in PBMCs and MH7A were validated by RT-PCR and found to be increased, similar to our sequence data, while after treatment with different concentrations of TP, IGF2BP3 mRNA level was decreased (Figures 6G,H). These results suggested that IGF2BP3 might be a new potential target of TP during the treatment of RA.

FIGURE 6

FIGURE 6. Potential target analysis of TP. (A) Network analysis between genes with great variation in both the m⁶A modification and mRNA levels and target proteins of TP. (B–E) Gene expression of IGF2BP3 in synovial tissue samples of GEO datasets. (F) Molecular docking of IGF2BP3 and TP. In the picture on the right, thick sticks represent the TP molecule, and thin lines represent residues in the protein binding site. (G,H) mRNA expression of IGF2BP3 in PBMCs and MH7A. *p < 0.05, **p < 0.01. TP, triptolide; IGF2BP3, insulin-like growth factor 2 mRNA-binding protein 3; GEO, Gene Expression Omnibus; and PBMCs, peripheral blood mononuclear cells.

Discussion

RA is an autoimmune disease that is characterized by chronic inflammation and progressive bone destruction of joints (Smolen et al., 2016). Multiple studies have suggested that the interplay among immunity, lifestyle, and genetic and environmental factors plays an important role in RA pathogenesis. Recently, an increasing body of evidence has revealed that epigenetics contributes to the pathogenesis of RA, probably by integrating genetic and environmental effects (Gray, 2014; Smolen et al., 2016). m⁶A, an epigenetic regulatory mechanism, is the most abundant internal modification of mRNA in eukaryotes and modulates the gene expression program at the posttranscriptional level. Variant m⁶A levels exert various biological functions such as cell differentiation, neuronal signaling, and immune tolerance. m⁶A modification is involved in various diseases, such as obesity, diabetes, and SLE (Church et al., 2010; Shen et al., 2015; Luo et al., 2020a). A recent study showed that METTL3 can promote activation and inflammation in FLSs and the adjuvant-induced arthritis animal model using nuclear factor (NF)-κB signaling pathway (Shi et al., 2021). Another study revealed that the m⁶A modification of TGM2 mRNA contributes to the inhibitory activity of sarsasapogenin in RA FLS (Lin et al., 2022). Nevertheless, as a complex disease, the effect of m⁶A modification on RA progression remains to be clarified.

Numerous studies have shown that peripheral blood represents the main highway for the RA immune system. PBMCs are immune cells and can initiate the autoimmune inflammatory process directed against target organs. Therefore, the gene expression signatures of PBMCs can shed light on the molecular features of the immune cells in the targeted organs in RA progression (Wang et al., 2022). Our research obtained the first MeRIP-seq profile in RA PBMCs. MeRIP-seq revealed 41,026 and 41,184 m⁶A peaks from 22,080 to 22,206 m⁶A-modified transcripts in Con and RA PBMCs, respectively. We found the enrichment of the m⁶A peaks of RA and Con PBMCs in the 3′-UTR area and near the stop codon, which were conformed with the results from past studies (Dominissini et al., 2012; Meyer et al., 2012; He and He, 2021). These sites were m⁶A-specific. Compared with the Con group, 295 hypermethylated m⁶A peaks and 288 hypomethylated m⁶A peaks in the RA PBMCs were identified. These differentially methylated peaks were greatly enriched in “signal transduction,” “cell differentiation,” “regulation of transcription,” “multicellular organism development,” “cell motility,” and “cellular community-eukaryotes,” suggesting the conserve and basic effects of m⁶A in controlling growth and cell fate specification (Zaccara et al., 2019). In KEGG pathway enrichment, it was observed that upregulated m⁶A modification genes are associated with the “calcium signaling pathway” and the “B-cell receptor signaling pathway.” Calcium (Ca²⁺) signaling is a common hub for many signaling pathways that take part in tolerance mechanisms, and B cells express a large panel of ion channels and transporters that contribute to the production and regulation of Ca²⁺ signals, which exert a significant effect on B-cell growth and fate (Hemon et al., 2017). Signaling by intracellular Ca²⁺ has been shown to be involved in RA pathogenesis. Autoreactive B cells also exert a key effect on the pathogenesis of RA. Defects in Ca²⁺ influx kinetics and enhancement of basal [Ca²⁺]_cyt have been found in T lymphocytes isolated from RA patients (Toldi et al., 2013). These studies suggested that m⁶A modification may boost the appearance and progression of RA by controlling the calcium signaling pathway or the B-cell receptor signaling pathway. Interestingly, the hypomethylated genes were significantly related to “sphingolipid metabolism.” A previous study reported that sphingolipids are a category of complicated lipids with a backbone of sphingoid bases and that several bioactive sphingolipids are involved in the pathological processes of RA (Gomez-Larrauri et al., 2020), which indicates that m⁶A modification may affect RA progression by regulating sphingolipid metabolism. However, the precise functions of sphingolipid metabolism in RA still need to be studied in-depth, and the m⁶A-altered genes associated with this pathway in RA must be further clarified.

Next, through RNA-seq analyses, 1,570 significantly dysregulated genes were identified, including 539 upregulated genes and 1,031 downregulated genes. These dysregulated genes were significantly involved in “signal transduction,” “cell adhesion” and “neutrophil degranulation.” In KEGG pathway analyses, several canonical pathways associated with RA inflammation were significantly enriched, including “complement and coagulation cascades,” “Jak-STAT signaling pathway” and “IL-17 signaling pathway.” However, the downregulated genes were significantly related to “natural killer cell mediated cytotoxicity.”

To elucidate the mechanism of m⁶A modification in influencing RA development, conjoint explorations of MeRIP-seq and RNA-seq data were performed, and 35 DMEGs were found. The enrichment of genes in “protein transport” and “protein binding” was observed, and KEGG pathway analysiss howed that these DMEGs were enriched primarily in “lysosome.” Lysosomes are membrane-bound organelles involved in the processes of degrading and recycling cellular waste, cellular signaling, and energy metabolism (Bonam et al., 2019). Increasing evidence implicates the effects of lysosomal dysfunction in more common diseases such as inflammatory and autoimmune disorders, including RA. Weissmann G concluded that various cathepsins, acid phosphatases, and other enzymes in humans participate in inflammation and joint damage (Weissmann, 1966). The earlier results indicated that the “lysosome” function in which DMEGs are involved might play an important role in RA pathogenesis and might be a therapeutic target.

TP is the main active ingredient extracted from TwHF, and it has been considered a promising anti-RA drug (Fan et al., 2018), but the underlying molecular mechanism of TP on RA has still not been fully elucidated. To better discover potential targets of TP and elucidate the molecular mechanism of TP in our further study, we chose TP as a therapeutic candidate and constructed the molecular network of TP target proteins and candidate DMEGs obtained from our MeRIP-seq and RNA-seq results. Our results showed that there was a close relationship between the target proteins of TP and the DMEGs of RA. The TUBB2A, IGF2BP3, DYNC1I1, and FOSL1 had high-degree nodes with the target proteins of TP in the network, suggesting that they might exert key effects on the response to TP. Considering the critical role of synovial and FLSs in RA (Nygaard and Firestein, 2020), we searched the GEO database and found four different datasets of RA synovial tissue in which the gene expression of IGF2BP3 was increased. This result was consistent with the trend of our sequencing results of RA PBMCs and suggested that IGF2BP3 might be an important virulence gene in RA.

IGF2BP3 is a highly conserved paralog of IGF2BPs and primarily exerts oncogenic effects in cancer (Mancarella and Scotlandi, 2019). In recent years, some studies have increasingly recorded the effects of IGF2BP3 on basic processes in cancer biology, and its overexpression has been extensively related to adverse patient results in various tumors, such as colon cancer and gastric carcinogenesis (Zhou et al., 2017; Yang et al., 2020). However, studies on IGF2BP3 in RA are rare. In our results, the docking scores between IGF2BP3 and TP is -8.6 kcal/mol, indicating that TP and IGF2BP3 have a higher binding affinity. The mRNA expression of IGF2BP3 in PBMCs and MH7A was increased, similar to our sequence data, while after treatment with different concentrations of TP, IGF2BP3 mRNA level was decreased, demonstrating that IGF2BP3 might be a new target of TP, which provide a possible strategy for the treatment of RA. Of course, elucidating the fundamental mechanisms by which m⁶A and TP regulate gene expression is an important direction of our future studies.

Conclusion

This study established a transcriptional map of m⁶A in RA PBMCs and displayed the hidden association between m⁶A methylation modification and associated genes in RA. The results highlight the importance of m⁶A modification as a gene regulatory system in RA. IGF2BP3 might be a potential therapeutic target of TP during the treatment of RA. More research based on these epigenetic clues should be conducted in the future to deepen the understanding of the mechanisms of TP and m⁶A modification in RA.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: National Center for Biotechnology Information (NCBI) BioProject database under accession number GSE193193.

Ethics Statement

The studies involving human participants were reviewed and approved by the ethical committee of China-Japan Friendship Hospital. The patients/participants provided their written informed consent to participate in this study.

Author Contributions

CX and YXu designed the framework of the article. DF wrote the manuscript and analyzed the data. BL contributed to the data analyses. XG, QZ, QY, and ZW took part in discussions related to the manuscript. XX, YXi, and QW performed the experiments. The manuscript was revised and approved by CX and BW. All authors made contributions to the article and approved the submitted version.

Funding

The National Natural Science Foundation of China (Grant number 82073677), Elite Medical Professionals Project of China-Japan Friendship Hospital (No. ZRJY2021-TD01), and Beijing Chinese Medicine Science and Technology Development Fund Project (Youth Planning Project): QN-2020-32 financially supported this research.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

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

References

Aletaha, D., Neogi, T., Silman, A. J., Funovits, J., Felson, D. T., Bingham, C. O., et al. (2010). 2010 Rheumatoid Arthritis Classification Criteria: an American College of Rheumatology/European League against Rheumatism Collaborative Initiative. Ann. Rheum. Dis. 69 (9), 1580–1588. doi:10.1136/ard.2010.138461

PubMed Abstract | CrossRef Full Text | Google Scholar

Alexander, N., Woetzel, N., and Meiler, J. (2011). bcl::Cluster : A Method for Clustering Biological Molecules Coupled with Visualization in the Pymol Molecular Graphics System. IEEE Int. Conf. Comput. Adv. Bio Med. Sci. 2011, 13–18. doi:10.1109/ICCABS.2011.5729867

PubMed Abstract | CrossRef Full Text | Google Scholar

Bailey, T. L., Boden, M., Buske, F. A., Frith, M., Grant, C. E., Clementi, L., et al. (2009). MEME SUITE: Tools for Motif Discovery and Searching. Nucleic Acids Res. 37, W202–W208. doi:10.1093/nar/gkp335

PubMed Abstract | CrossRef Full Text | Google Scholar

Bonam, S. R., Wang, F., and Muller, S. (2019). Lysosomes as a Therapeutic Target. Nat. Rev. Drug Discov. 18 (12), 923–948. doi:10.1038/s41573-019-0036-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Burley, S. K., Bhikadiya, C., Bi, C., Bittrich, S., Chen, L., Crichlow, G. V., et al. (2021). RCSB Protein Data Bank: Powerful New Tools for Exploring 3D Structures of Biological Macromolecules for Basic and Applied Research and Education in Fundamental Biology, Biomedicine, Biotechnology, Bioengineering and Energy Sciences. Nucleic Acids Res. 49 (D1), D437–D451. doi:10.1093/nar/gkaa1038

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, S., Zhou, Y., Chen, Y., and Gu, J. (2018). Fastp: an Ultra-fast All-In-One FASTQ Preprocessor. Bioinformatics 34 (17), i884–i890. doi:10.1093/bioinformatics/bty560

PubMed Abstract | CrossRef Full Text | Google Scholar

Church, C., Moir, L., Mcmurray, F., Girard, C., Banks, G. T., Teboul, L., et al. (2010). Overexpression of Fto Leads to Increased Food Intake and Results in Obesity. Nat. Genet. 42 (12), 1086–1092. doi:10.1038/ng.713

PubMed Abstract | CrossRef Full Text | Google Scholar

Cross, M., Smith, E., Hoy, D., Carmona, L., Wolfe, F., Vos, T., et al. (2014). The Global burden of Rheumatoid Arthritis: Estimates from the Global burden of Disease 2010 Study. Ann. Rheum. Dis. 73 (7), 1316–1322. doi:10.1136/annrheumdis-2013-204627

PubMed Abstract | CrossRef Full Text | Google Scholar

Dominissini, D., Moshitch-Moshkovitz, S., Salmon-Divon, M., Amariglio, N., and Rechavi, G. (2013). Transcriptome-wide Mapping of N(6)-methyladenosine by m(6)A-Seq Based on Immunocapturing and Massively Parallel Sequencing. Nat. Protoc. 8 (1), 176–189. doi:10.1038/nprot.2012.148

PubMed Abstract | CrossRef Full Text | Google Scholar

Dominissini, D., Moshitch-Moshkovitz, S., Schwartz, S., Salmon-Divon, M., Ungar, L., Osenberg, S., et al. (2012). Topology of the Human and Mouse m6A RNA Methylomes Revealed by m6A-Seq. Nature 485 (7397), 201–206. doi:10.1038/nature11112

PubMed Abstract | CrossRef Full Text | Google Scholar

Doody, K. M., Bottini, N., and Firestein, G. S. (2017). Epigenetic Alterations in Rheumatoid Arthritis Fibroblast-like Synoviocytes. Epigenomics 9 (4), 479–492. doi:10.2217/epi-2016-0151

PubMed Abstract | CrossRef Full Text | Google Scholar

Fan, D., Guo, Q., Shen, J., Zheng, K., Lu, C., Zhang, G., et al. (2018). The Effect of Triptolide in Rheumatoid Arthritis: From Basic Research towards Clinical Translation. Int. J. Mol. Sci. 19 (2). doi:10.3390/ijms19020376

PubMed Abstract | CrossRef Full Text | Google Scholar

Fan, D., He, X., Bian, Y., Guo, Q., Zheng, K., Zhao, Y., et al. (2016). Triptolide Modulates TREM-1 Signal Pathway to Inhibit the Inflammatory Response in Rheumatoid Arthritis. Int. J. Mol. Sci. 17 (4), 498. doi:10.3390/ijms17040498

PubMed Abstract | CrossRef Full Text | Google Scholar

Fu, Y., Dominissini, D., Rechavi, G., and He, C. (2014). Gene Expression Regulation Mediated through Reversible m⁶A RNA Methylation. Nat. Rev. Genet. 15 (5), 293–306. doi:10.1038/nrg3724

PubMed Abstract | CrossRef Full Text | Google Scholar

Goldbach-Mansky, R., Wilson, M., Fleischmann, R., Olsen, N., Silverfield, J., Kempf, P., et al. (2009). Comparison of Tripterygium Wilfordii Hook F versus Sulfasalazine in the Treatment of Rheumatoid Arthritis: a Randomized Trial. Ann. Intern. Med. 151 (4), 229W49–4051. doi:10.7326/0003-4819-151-4-200908180-00005

CrossRef Full Text | Google Scholar

Gomez-Larrauri, A., Presa, N., Dominguez-Herrera, A., Ouro, A., Trueba, M., and Gomez-Muñoz, A. (2020). Role of Bioactive Sphingolipids in Physiology and Pathology. Essays Biochem. 64 (3), 579–589. doi:10.1042/EBC20190091

PubMed Abstract | CrossRef Full Text | Google Scholar

Gray, S. G. (2014). Epigenetic-based Immune Intervention for Rheumatic Diseases. Epigenomics 6 (2), 253–271. doi:10.2217/epi.13.87

PubMed Abstract | CrossRef Full Text | Google Scholar

Han, R., Rostami-Yazdi, M., Gerdes, S., and Mrowietz, U. (2012). Triptolide in the Treatment of Psoriasis and Other Immune-Mediated Inflammatory Diseases. Br. J. Clin. Pharmacol. 74 (3), 424–436. doi:10.1111/j.1365-2125.2012.04221.x

CrossRef Full Text | Google Scholar

He, P. C., and He, C. (2021). m6 A RNA Methylation: from Mechanisms to Therapeutic Potential. EMBO J. 40 (3), e105977. doi:10.15252/embj.2020105977

PubMed Abstract | CrossRef Full Text | Google Scholar

Heinz, S., Benner, C., Spann, N., Bertolino, E., Lin, Y. C., Laslo, P., et al. (2010). Simple Combinations of Lineage-Determining Transcription Factors Prime Cis-Regulatory Elements Required for Macrophage and B Cell Identities. Mol. Cel. 38 (4), 576–589. doi:10.1016/j.molcel.2010.05.004

CrossRef Full Text | Google Scholar

Hemon, P., Renaudineau, Y., Debant, M., Le Goux, N., Mukherjee, S., Brooks, W., et al. (2017). Calcium Signaling: From Normal B Cell Development to Tolerance Breakdown and Autoimmunity. Clin. Rev. Allergy Immunol. 53 (2), 141–165. doi:10.1007/s12016-017-8607-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Jacobsson, L. T., Jacobsson, M. E., Askling, J., and Knowler, W. C. (2003). Perinatal Characteristics and Risk of Rheumatoid Arthritis. BMJ 326 (7398), 1068–1069. doi:10.1136/bmj.326.7398.1068

PubMed Abstract | CrossRef Full Text | Google Scholar

Karlson, E. W., Mandl, L. A., Hankinson, S. E., and Grodstein, F. (2004). Do breast-feeding and Other Reproductive Factors Influence Future Risk of Rheumatoid Arthritis? Results from the Nurses' Health Study. Arthritis Rheum. 50 (11), 3458–3467. doi:10.1002/art.20621

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, D., Langmead, B., and Salzberg, S. L. (2015). HISAT: a Fast Spliced Aligner with Low Memory Requirements. Nat. Methods 12 (4), 357–360. doi:10.1038/nmeth.3317

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, S., Chen, J., Cheng, T., Gindulyte, A., He, J., He, S., et al. (2021). PubChem in 2021: New Data Content and Improved Web Interfaces. Nucleic Acids Res. 49 (D1), D1388–D1395. doi:10.1093/nar/gkaa971

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 m6A 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

Lin, X., Tao, C., Zhang, R., Zhang, M., Wang, Q., and Chen, J. (2022). N6-methyladenosine Modification of TGM2 mRNA Contributes to the Inhibitory Activity of Sarsasapogenin in Rheumatoid Arthritis Fibroblast-like Synoviocytes. Phytomedicine 95, 153871. doi:10.1016/j.phymed.2021.153871

PubMed Abstract | CrossRef Full Text | Google Scholar

Loos, R. J., and Yeo, G. S. (2014). The Bigger Picture of FTO: the First GWAS-Identified Obesity Gene. Nat. Rev. Endocrinol. 10 (1), 51–61. doi:10.1038/nrendo.2013.227

PubMed Abstract | CrossRef Full Text | Google Scholar

Luo, Q., Gao, Y., Zhang, L., Rao, J., Guo, Y., Huang, Z., et al. (2020a). Decreased ALKBH5, FTO, and YTHDF2 in Peripheral Blood Are as Risk Factors for Rheumatoid Arthritis. Biomed. Res. Int. 2020, 5735279. doi:10.1155/2020/5735279

PubMed Abstract | CrossRef Full Text | Google Scholar

Luo, Q., Rao, J., Zhang, L., Fu, B., Guo, Y., Huang, Z., et al. (2020b). The Study of METTL14, ALKBH5, and YTHDF2 in Peripheral Blood Mononuclear Cells from Systemic Lupus Erythematosus. Mol. Genet. Genomic Med. 8 (9), e1298. doi:10.1002/mgg3.1298

PubMed Abstract | CrossRef Full Text | Google Scholar

Lv, Q. W., Zhang, W., Shi, Q., Zheng, W. J., Li, X., Chen, H., et al. (2015). Comparison of Tripterygium Wilfordii Hook F with Methotrexate in the Treatment of Active Rheumatoid Arthritis (TRIFRA): a Randomised, Controlled Clinical Trial. Ann. Rheum. Dis. 74 (6), 1078–1086. doi:10.1136/annrheumdis-2013-204807

PubMed Abstract | CrossRef Full Text | Google Scholar

Mancarella, C., and Scotlandi, K. (2019). IGF2BP3 from Physiology to Cancer: Novel Discoveries, Unsolved Issues, and Future Perspectives. Front. Cel Dev Biol 7, 363. doi:10.3389/fcell.2019.00363

PubMed Abstract | CrossRef Full Text | Google Scholar

Meng, J., Lu, Z., Liu, H., Zhang, L., Zhang, S., Chen, Y., et al. (2014). A Protocol for RNA Methylation Differential Analysis with MeRIP-Seq Data and exomePeak R/Bioconductor Package. Methods 69 (3), 274–281. doi:10.1016/j.ymeth.2014.06.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Meyer, K. D., Saletore, Y., Zumbo, P., Elemento, O., Mason, C. E., and Jaffrey, S. R. (2012). Comprehensive Analysis of mRNA Methylation Reveals Enrichment in 3' UTRs and Near Stop Codons. Cell 149 (7), 1635–1646. doi:10.1016/j.cell.2012.05.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Nygaard, G., and Firestein, G. S. (2020). Restoring Synovial Homeostasis in Rheumatoid Arthritis by Targeting Fibroblast-like Synoviocytes. Nat. Rev. Rheumatol. 16 (6), 316–333. doi:10.1038/s41584-020-0413-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Pedersen, M., Jacobsen, S., Garred, P., Madsen, H. O., Klarlund, M., Svejgaard, A., et al. (2007). Strong Combined Gene-Environment Effects in Anti-cyclic Citrullinated Peptide-Positive Rheumatoid Arthritis: a Nationwide Case-Control Study in Denmark. Arthritis Rheum. 56 (5), 1446–1453. doi:10.1002/art.22597

PubMed Abstract | CrossRef Full Text | Google Scholar

Pertea, M., Pertea, G. M., Antonescu, C. M., Chang, T. C., Mendell, J. T., and Salzberg, S. L. (2015). StringTie Enables Improved Reconstruction of a Transcriptome from RNA-Seq Reads. Nat. Biotechnol. 33 (3), 290–295. doi:10.1038/nbt.3122

PubMed Abstract | CrossRef Full Text | Google Scholar

Prevoo, M. L., van 't Hof, M. A., Kuper, H. H., van Leeuwen, M. A., van de Putte, L. B., and van Riel, P. L. (1995). Modified Disease Activity Scores that Include Twenty-Eight-Joint Counts. Development and Validation in a Prospective Longitudinal Study of Patients with Rheumatoid Arthritis. Arthritis Rheum. 38 (1), 44–48. doi:10.1002/art.1780380107

PubMed Abstract | CrossRef Full Text | Google Scholar

Robinson, M. D., Mccarthy, D. J., and Smyth, G. K. (2010). edgeR: a Bioconductor Package for Differential Expression Analysis of Digital Gene Expression Data. Bioinformatics 26 (1), 139–140. doi:10.1093/bioinformatics/btp616

PubMed Abstract | CrossRef Full Text | Google Scholar

Scherer, H. U., Häupl, T., and Burmester, G. R. (2020). The Etiology of Rheumatoid Arthritis. J. Autoimmun. 110, 102400. doi:10.1016/j.jaut.2019.102400

CrossRef Full Text | Google Scholar

Shen, F., Huang, W., Huang, J. T., Xiong, J., Yang, Y., Wu, K., et al. (2015). Decreased N(6)-methyladenosine in Peripheral Blood RNA from Diabetic Patients Is Associated with FTO Expression rather Than ALKBH5. J. Clin. Endocrinol. Metab. 100 (1), E148–E154. doi:10.1210/jc.2014-1893

PubMed Abstract | CrossRef Full Text | Google Scholar

Shi, W., Zheng, Y., Luo, S., Li, X., Zhang, Y., Meng, X., et al. (2021). METTL3 Promotes Activation and Inflammation of FLSs through the NF-Κb Signaling Pathway in Rheumatoid Arthritis. Front. Med. (Lausanne) 8, 607585. doi:10.3389/fmed.2021.607585

PubMed Abstract | CrossRef Full Text | Google Scholar

Smolen, J. S., Aletaha, D., and Mcinnes, I. B. (2016). Rheumatoid Arthritis. Lancet 388 (10055), 2023–2038. doi:10.1016/S0140-6736(16)30173-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Sparks, J. A., Chang, S. C., Deane, K. D., Gan, R. W., Kristen Demoruelle, M., Feser, M. L., et al. (2016). Associations of Smoking and Age with Inflammatory Joint Signs Among Unaffected First-Degree Relatives of Rheumatoid Arthritis Patients: Results from Studies of the Etiology of Rheumatoid Arthritis. Arthritis Rheumatol. 68 (8), 1828–1838. doi:10.1002/art.39630

PubMed Abstract | CrossRef Full Text | Google Scholar

Tanchuk, V. Y., Tanin, V. O., Vovk, A. I., and Poda, G. (2016). A New, Improved Hybrid Scoring Function for Molecular Docking and Scoring Based on AutoDock and AutoDock Vina. Chem. Biol. Drug Des. 87 (4), 618–625. doi:10.1111/cbdd.12697

PubMed Abstract | CrossRef Full Text | Google Scholar

Tao, X., Younger, J., Fan, F. Z., Wang, B., and Lipsky, P. E. (2002). Benefit of an Extract of Tripterygium Wilfordii Hook F in Patients with Rheumatoid Arthritis: a Double-Blind, Placebo-Controlled Study. Arthritis Rheum. 46 (7), 1735–1743. doi:10.1002/art.10411

PubMed Abstract | CrossRef Full Text | Google Scholar

Toldi, G., Bajnok, A., Dobi, D., Kaposi, A., Kovács, L., Vásárhelyi, B., et al. (2013). The Effects of Kv1.3 and IKCa1 Potassium Channel Inhibition on Calcium Influx of Human Peripheral T Lymphocytes in Rheumatoid Arthritis. Immunobiology 218 (3), 311–316. doi:10.1016/j.imbio.2012.05.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Trott, O., and Olson, A. J. (2010). AutoDock Vina: Improving the Speed and Accuracy of Docking with a New Scoring Function, Efficient Optimization, and Multithreading. J. Comput. Chem. 31 (2), 455–461. doi:10.1002/jcc.21334

CrossRef Full Text | Google Scholar

Wang, J., Yan, S., Lu, H., Wang, S., and Xu, D. (20192019). METTL3 Attenuates LPS-Induced Inflammatory Response in Macrophages via NF-Κb Signaling Pathway. Mediators Inflamm. 2019, 3120391. doi:10.1155/2019/3120391

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, S., Zuo, S., Liu, Z., Ji, X., Yao, Z., and Wang, X. (2018). Study on the Efficacy and Mechanism of Triptolide on Treating TNF Transgenic Mice with Rheumatoid Arthritis. Biomed. Pharmacother. 106, 813–820. doi:10.1016/j.biopha.2018.07.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., Xie, X., Zhang, C., Su, M., Gao, S., Wang, J., et al. (2022). Rheumatoid Arthritis, Systemic Lupus Erythematosus and Primary Sjögren's Syndrome Shared Megakaryocyte Expansion in Peripheral Blood. Ann. Rheum. Dis. 81 (3), 379–385. doi:10.1136/annrheumdis-2021-220066

PubMed Abstract | CrossRef Full Text | Google Scholar

Weissmann, G. (1966). Lysosomes and Joint Disease. Arthritis Rheum. 9 (6), 834–840. doi:10.1002/art.1780090611

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, Z., Wang, T., Wu, D., Min, Z., Tan, J., and Yu, B. (2020). RNA N6-Methyladenosine Reader IGF2BP3 Regulates Cell Cycle and Angiogenesis in colon Cancer. J. Exp. Clin. Cancer Res. 39 (1), 203. doi:10.1186/s13046-020-01714-8

CrossRef Full Text | Google Scholar

Yu, G., Wang, L. G., and He, Q. Y. (2015). ChIPseeker: an R/Bioconductor Package for ChIP Peak Annotation, Comparison and Visualization. Bioinformatics 31 (14), 2382–2383. doi:10.1093/bioinformatics/btv145

PubMed Abstract | CrossRef Full Text | Google Scholar

Zaccara, S., Ries, R. J., and Jaffrey, S. R. (2019). Reading, Writing and Erasing mRNA Methylation. Nat. Rev. Mol. Cel Biol 20 (10), 608–624. doi:10.1038/s41580-019-0168-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, X., Zhang, D., Jia, H., Feng, Q., Wang, D., Liang, D., et al. (2015). The Oral and Gut Microbiomes Are Perturbed in Rheumatoid Arthritis and Partly Normalized after Treatment. Nat. Med. 21 (8), 895–905. doi:10.1038/nm.3914

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, Y., Huang, T., Siu, H. L., Wong, C. C., Dong, Y., Wu, F., et al. (2017). IGF2BP3 Functions as a Potential Oncogene and Is a Crucial Target of miR-34a in Gastric Carcinogenesis. Mol. Cancer 16 (1), 77. doi:10.1186/s12943-017-0647-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, Y. Z., Zhao, L. D., Chen, H., Zhang, Y., Wang, D. F., Huang, L. F., et al. (2018). Comparison of the Impact of Tripterygium Wilfordii Hook F and Methotrexate Treatment on Radiological Progression in Active Rheumatoid Arthritis: 2-year Follow up of a Randomized, Non-blinded, Controlled Study. Arthritis Res. Ther. 20 (1), 70. doi:10.1186/s13075-018-1563-6

PubMed Abstract | CrossRef Full Text | Google Scholar