Four RA-associated microarrays (GSE12021, GSE55235, GSE55457, and GSE55584) in the GPL96 platform were downloaded from Gene Expression Omnibus (GEO). Then RNA-seq data from RA and normal synovial tissues in the four microarrays were extracted and merged, and batch effect was reduced by the ComBat package in R software (www.datavis.ca/R/). Similarly, GSE48780 and GSE77298 were also retrieved from the GPL570 platform, and the transcriptome data were downloaded and merged before the batch effect was reduced. Moreover, RNA-seq and clinical data in GSE13026, GSE21537, and GSE45867 were also used to explore the relationships between ATF3 and RA status and the pharmacotherapy response.

The “edgeR” package in R software was used to screen out the differentially expressed genes with |log Fold-Change| ≥ 1 and false discovery rate (FDR) <0.05 in the GPL96 and GPL570 databases. After that, Gene Ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were performed to investigate enriched potential function and pathways in these differentially expressed genes with the “clusterProfiler” package in R software. Moreover, the differentially expressed genes were also uploaded to the STRING website (www.string-db.org) to visualize the function protein association network with a minimum required interaction score larger than 0.99.

WGCNA, a scale-free network algorithm based on gene expression microarray data, is widely used for finding clusters of highly correlated genes in bioinformatics applications. First, the genes were screened with SD larger than 0.5 among all samples and were imported to the co-expression network. In view of data integrity and little batch effect in the GPL96 and GPL570 databases, the sample trees were constructed with cut height = 2,000, and the selected samples were applied to build the sample dendrogram and clinical trait heatmap. Second, the pickSoftThreshold function in the “WGCNA” package in R software was used to estimate an optimum soft threshold power, which struck a balance between higher R ² and higher mean connectivity. Third, the adjacency values among the imported genes were calculated according to the optimum soft threshold power in both databases and were transformed to topological overlap matrix (TOM). Furthermore, the corresponding dissimilarity (1 − TOM) was also calculated, which was significant for hierarchical clustering of genes. Fourth, modules were identified with the dynamic tree cut method by hierarchical clustering of genes, which used dissimilarity (1 − TOM) as distance measurement with module size larger than 50 and deep split = 2. Fifth, the most clinically relevant module was selected, and the module membership (MM) and gene significance (GS) were calculated. After that, the two differentially expressed gene sets and two clinically relevant modules in the GPL96 and GPL570 databases were intersected, and it was found that ATF3 was a differentially expressed clinical-relevant gene, which might be a potential biomarker worthy of validation by in vitro cell experiments.

The transcriptome data and clinical data in GSE13026, GSE21537, and GSE45867 were downloaded to explore the relationships between ATF3 expression and RA status, lymphocyte aggregation, inflammation status, infliximab response, tocilizumab response, methotrexate response, joint location, age, and gender. In the light of the median expression of ATF3, all patients in the microarrays were divided into two groups, and then the clinical parameters of both the high and low ATF3 expression groups were compared. Moreover, ATF3 expression was also compared among groups classified by clinical parameters.

All transcriptome data in the GPL96 and GPL570 databases were normalized before being imported to the CIBERSORT tool to estimate the contents of 22 human immune cells in each RA patient. Then, the difference of infiltrating immune cells in the high and low ATF3 expression groups was observed and compared. A p < 0.05 was considered statistically significant.

Moreover, the RNA-seq data in both the high and low ATF3 expression groups were uploaded to Gene Set Enrichment Analysis (GSEA) to explore the differentially expressed gene-related signaling pathways. The enriched pathways were screened out based on an FDR < 0.25 and p < 0.05 after 1,000 permutations, and the top 10 enriched pathways in both groups were presented in multi-GSEA plots.

Human RA cell lines RA-FLS and MH7A (Beijing Beina Chuanglian Biotechnology Institute, Beijing, China) were cultured in Dulbecco’s Modified Eagle Medium/nutrient mixture F-12 (DMEM/F12; Gibco, Grand Island, NY, USA) supplemented with 10% heat-inactivated fetal bovine serum (FBS; Gibco, USA; 1% penicillin/streptomycin) in a humidified incubator containing 5% CO₂ at 37°C according to American Type Culture Collection suggestions.

Cells were added to the 6-well plate at a density of 1 × 10⁶/well and incubated for 24 h. Then, a fresh medium containing lentivirus and 8 μg/ml of polybrene [OBiO Technology (Shanghai) Corp., Ltd., Shanghai, China] was added. The addition amount of lentivirus was calculated according to the required multiplicity of infection (MOI). Preliminary experiments determined that the optimal MOI for lentivirus of short hairpin (Sh)-ATF3 (Sh-ATF3), the lentivirus of control (Sh-Ctrl), ATF3-overexpression (ATF3-OE), and control lentivirus (Ctrl-OE) was 100. The culture medium was changed within 24 h after infection, and 1‰ of puromycin (Sigma, St. Louis, MO, USA) was added to the complete medium to kill the wild-type RA-FLS and MH7A. About 14 days later, the majority of the wild-type cells died. The efficiency of infection was determined by quantitative real-time PCR.

The gene expression levels of Sh-ATF3 and ATF3-OE were evaluated by qRT-PCR. Total RNA was isolated using Invitrogen TRIzol reagent (Thermo Fisher Scientific, Waltham, MA, USA) according to the manufacturer’s instructions. The isolated total RNA was then reverse transcribed into cDNA using the QuantScript RT cDNA Synthesis Kit (Bio-Rad, Hercules, CA, USA). The Maxima SYBR Green qPCR kit (Thermo Fisher Scientific, USA) was used for qRT-PCR analysis. The oligonucleotide primer sequences including forward and reverse used in these experiments were presented as follows: ATF3: 5′-CCT​CTG​CGC​TGG​AAT​CAG​TC-3′ and 5′-TTC​TTT​CTC​GTC​GTC​GCC​TCT​TTT​T-3′, IL-1β: 5′-GTC​GGA​GAT​TCG​TAG​CTG​GA-3′ and 5′-GTC​GGA​GAT​TCG​TAG​CTG​GA-3′, IL-6: 5′-ACT​CAC​CTC​TTC​AGA​ACG​AAT​TG-3′ and 5′-CCA​TCT​TTG​GAA​GGT​TCA​GGT​TG-3′, IL-8: 5′-ACT​GAG​AGT​GAT​TGA​GAG​TGG​AC-3′ and 5′-AAC​CCT​CTG​CAC​CCA​GTT​TTC-3′, and β-actin: 5′-CAT​GTA​CGT​TGC​TAT​CCA​GGC-3′ and 5′-CTC​CTT​AAT​GTC​ACG​CAC​GAT-3′. The relative mRNA level was calculated by the 2^−ΔΔCt method.

Stable cell lines were constructed using the ATF3 lentivirus. The cells were plated into 96-well plates at 1 × 10³ cells/well. Cell viability was measured at 1, 2, 3, and 4 days using Cell Counting Kit-8 (CCK-8) reagent (96992; Sigma, St. Louis, MO, USA) according to the manufacturer’s instructions. Optical density (OD) was measured at 450 nm using an ELx800 Absorbance Microplate Reader (BioTek, Winooski, VT, USA).

RA-FLS and MH7A were detached as single-cell suspension after 48-h transfection and fixed in 75% ethanol overnight at 4°C. Fixed cells were washed and stained with 25 mg/ml of propidium iodide (PI; Sigma, USA) in phosphate-buffered saline (PBS) containing 0.1% Triton and 10 mg/ml of RNase (Thermo Fisher Scientific, USA) on ice for 0.5 h in the dark. The cell cycle was analyzed by flow cytometry on a FACSCalibur system (BD Biosciences, San Jose, CA, USA). The percentage of cells in different cell cycle phases was calculated by ModFit 3.0 (Verity Software House, Inc., Topsham, ME, USA).

The Annexin V-APC/7-AAD kit (Multisciences, Hangzhou, China) was used to visualize apoptotic cells according to the manufacturer’s instructions. RA-FLS and MH7A were infected with Sh-ATF3, Sh-Ctrl, ATF3-OE, and Ctrl-OE lentiviruses. Cells were washed with PBS and recombined in 1× binding buffer, then treated with Annexin VAPC and 7-AAD for 15 min at room temperature (RT) in the dark, and then analyzed by flow cytometry.

Cell migration and invasion (Matrigel added) were evaluated using the transwell assay (Corning, New York, NY, USA). Cells (4 × 10⁴ cells/well for migration and invasion assay of ATF3 knockdown and 2 × 10⁴ cells/well for migration and invasion assay of ATF3 overexpression) in serum-free DMEM were seeded into the upper chamber (pre-coated with Matrigel (R&D Systems, Minneapolis, MN, USA) for invasion assay), and media containing 10% FBS were added into the lower chamber. After being incubated for 24 h, the cells were fixed with 70% ethanol, stained with 0.1% crystal violet, and photographed under an inverted microscope. The number of cells in five random fields per well was calculated, and the average value was taken as the migration and invasion numbers.

Stable infection of RA-FLS and MH7A was treated with or without 10 ng/ml of tumor necrosis factor-α (TNF-α) (Propotec, R&D Systems, Minnesota, USA) for 24 h, and the supernatants were collected to determine the levels of IL-1β, IL-6, and IL-8 with ELISA kits.

All statistical analyses were performed in the R software and GraphPad Prism 9. Quantitative data of two groups were compared by Student’s t-test, and quantitative data of three or more groups were compared by one-way ANOVA or Welch’s test. A p < 0.05 was regarded as statistically significant.

The flowchart shows all essential and significant procedures in our study (Figure 1). The transcriptome data in GSE12021, GSE55235, GSE55457, and GSE55584 were downloaded and merged as the GPL96 database. Similarly, the RNA-seq data in GSE48780 and GSE77298 were also acquired and merged as the GPL570 database. The ComBat function was used to reduce batch effect among different microarrays in the two merged databases, and the normalizeBetweenArrays function was used to normalize gene expression among different samples in the GPL570 and GPL96 databases (Figures 2A,B). After data normalization, the heatmaps showed the top 50 differentially expressed genes between synovial tissues and RA tissues (Figures 2C,D). Meanwhile, the volcano plots showed that there were 163 and 707 differentially expressed genes with FDR < 0.05 and |log Fold-Change| ≥ 1 in the GPL570 and GPL96 databases, respectively (Figures 2E,F).

We used the “clusterProfiler” package in R software to explore the potential GO and KEGG function enrichment in those differentially expressed genes. The GO enrichment analysis showed that those differential genes in the GPL570 database were significantly enriched in the muscular system process, muscle contraction in the subset of biological process (BP); enriched in contractile fibers, myofibrils in the subset of cell compartment (CC); and enriched in the structural constituent of muscle, actin binding in the subset of molecular function (MF) (Figure 2G). Meanwhile, the GO enrichment analysis showed that those differential genes in the GPL96 database were significantly enriched in immune- and inflammation-associated functions, such as lymphocyte differentiation, T-cell differentiation in the subset of BP; external side of the plasma membrane, immunological synapse in the subset of CC; and cytokine receptor binding, chemokine activity in the subset of MF (Figure 2H).

Moreover, KEGG enrichment analysis showed that those differentially expressed genes in the two databases were all significantly enriched in some inflammation-associated pathways, such as cytokine–cytokine receptor interaction, chemokine signaling pathway, viral protein interaction with cytokine and cytokine receptor, and IL-17 signaling pathway (Figures 2I,J), which were in accordance with inflammatory attributes of RA.

Furthermore, the differentially expressed genes in the two databases were also uploaded to STRING website to explore the protein–protein interaction (PPI). Interestingly, we found that ATF3 was incorporated into the two PPI networks with a minimum required interaction score larger than 0.99 and had close relationships with JUN and ATF4 (Figures 3A,B).

One hundred six samples, including seven normal synovial tissues and 99 RA synovia tissues, in the GPL570 database were presented in the sample dendrogram and trait heatmap with cut height = 2,000 in order to ensure data integrity (Figure 4A). Similarly, the sample dendrogram and trait heatmap showed that 29 normal synovial tissues and 45 RA synovia tissues were included in the GPL96 database (Figure 4B). The pickSoftThreshold function in the “WGCNA” package was applied to estimate optimum soft threshold powers, which could strike a balance between R ² and mean connectivity and corresponded to five and six in the GPL570 and GPL96 databases, respectively (Figures 4C–F). After the soft threshold power was chosen, the adjacency value was calculated among the selected genes with SD > 0.5 and was transformed to TOM. The gene clustering dendrogram was constructed according to TOM-based dissimilarity, and the gene modules were identified and merged with the dynamic tree cut method (cut height = 0.25) by hierarchical clustering of genes in the GPL570 and GPL96 databases (Figures 4G–I and J–L).

The module–trait correlation heatmap showed that the ME-greenyellow module was the most positively correlated with normal synovial tissues and was the most negatively correlated with RA synovial tissues (Figure 5A). Meanwhile, the module–trait correlation heatmap of GPL96 indicated that ME-salmon module had the most positive correlation with normal samples and the most negative correlation with RA samples (Figure 5B). Moreover, the correlation plots between MM and GS showed that the genes in the ME-greenyellow module and ME-salmon module had a significant correlation with normal and RA samples and could be regarded as a potential set of biomarkers (Figures 5C–F).

Finally, we intersected the ME-greenyellow module genes and differentially expressed genes in the GPL570 database and the ME-salmon module genes and differentially expressed genes in the GPL96 database, and we found that only one gene, ATF3, was screened out (Figure 5G).

In addition to the transcriptome and associated clinical information of the GPL570 and GPL96 databases, the RNA-seq and clinical data in GSE13026, GSE21537, and GSE45867 were also downloaded to investigate the relationships between ATF3 expression and clinical parameters, such as age, gender, joint location, RA status, lymphocyte aggregation status in the synovium, inflammation status in the synovium, infliximab response, tocilizumab response, and methotrexate response. All patients in these microarrays were divided into the high and low ATF3 groups in the light of the median expression of ATF3.

In GSE13026, we found that ATF3 expression in long-standing RA patients is significantly lower than that in normal individuals, and its expression in long-standing RA patients was lower than that in early RA patients. Similarly, there was a higher proportion of normal synovium in the high ATF3 expression group (Figures 6A,B).

Figures 6C–F show that the lymphocyte aggregation status in the synovium seemed to have no relationships with ATF3 expression in GSE21537. Moreover, the clinical data in GSE21537 indicated that cases with good infliximab response accounted for a higher percentage in the high ATF3 group than in the low ATF3 group, although patients with good infliximab response just had a slightly higher ATF3 expression than those with no or moderate response (Figures 6G,H).

GSE48780 in the GPL570 database indicated that ATF3 expression was significantly different among different joint synovium and its expression in hip synovium was significantly higher than that in the knees and hands (Figures 6I,J). Furthermore, ATF3 expression in inflammation-infiltrated synovium was slightly higher than that in non-inflammation synovium, and patients with inflammation-infiltrated synovium accounted for a higher percentage in the high ATF3 group (Figures 6K,L). After that, in mining the clinical data of GSE48780, we also found that ATF3 expression seemed higher in female patients, although the difference between genders was not significant (Figures 6M,N).

In GSE45867, we explored the relationship between ATF3 and age, and we found that ATF3 had a slight decrease along with the increment of age (Figure 6O). Meanwhile, we discovered a positive correlation between ATF3 and DAS-28 score (Figure 6P), which was a significant index to evaluate RA activity and is widely used in clinical practice. Moreover, we also discovered that ATF3 expression in RA patients with a good response was significantly higher than that in patients with a limited response before treatment with tocilizumab or methotrexate (Figures 6Q,R). More interestingly, ATF3 expression in RA synovium after treatment with tocilizumab or methotrexate was significantly lower in patients with good response compared with those with a limited response (Figures 6S,T).

Taken together, in the light of the correlation between ATF3 and DAS-28 score, the relationship between ATF3 expression and methotrexate/tocilizumab response, and the change of ATF3 expression after tocilizumab treatment, we assumed that ATF3 expression in the synovium might promote RA progression.

According to the median expression of ATF3, patients in the GPL570 and GPL96 databases can be divided into the high and low ATF3 groups. GSEA showed that these pathways enriched in the high ATF3 group were associated with metabolisms, such as glycine serine and threonine metabolism, O-glycan biosynthesis, primary bile acid biosynthesis in the GPL570 database, and fatty acid metabolism, fructose and mannose metabolism, glycerophospholipid metabolism, and glycolysis gluconeogenesis in the GPL96 database (Figures 7A,C). Similarly, GSEA indicated that KEGG pathways enriched in the low ATF3 group were related to inflammation and immunoregulation, such as antigen processing and presentation, natural killer cell-mediated cytotoxicity, primary immunodeficiency in the GPL570 database, allograft rejection, antigen processing and presentation, and graft versus host disease in the GPL96 database (Figures 7B,D).

We applied the CIBERSORT tool to estimate the contents of 22 human immune cells in each RA patient. Figures 7E,F show the relative percentage of different infiltrated immune cells in RA samples in the GPL570 and GPL96 databases. Moreover, the immune cell correlation matrix was constructed, and we noticed that activated NK cells were the most positively correlated with resting mast cells, while resting mast cells were the most negatively correlated with activated mast cells in the GPL570 database. Similarly, monocytes were the most positively correlated with activated mast cells, and memory B cells were the most negatively correlated with naive B cells in the GPL96 database (Figures 7G,H).

Finally, we analyzed and compared the contents of infiltrated immune cells between the high and low ATF3 groups in the GPL570 and GPL96 databases (Figures 7I,J). Interestingly, we discovered that the contents of activated mast cells in the high ATF3 group were widely and significantly higher than those in the low ATF3 group in both databases.