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.

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.

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).

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).

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.

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.

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.

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.

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).

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.

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.

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.

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.

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.

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).

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.

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).

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).

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).

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.