Effect of astragalosides on long non-coding RNA expression profiles in rats with adjuvant-induced arthritis.

Rheumatoid arthritis (RA) is an autoimmune disease of unknown etiology, which occurs in ~1.0% of the general population. Increasing studies have suggested that long non‑coding RNAs (lncRNAs) may serve important roles in various biological processes and may be associated with the pathogenesis of different types of disease, including RA. Astragalosides (AST) has been used as a traditional Chinese medicine for the treatment of RA. However, the mechanism underlying its therapeutic effect has remained unclear to date. Thus, there is an urgent need to elucidate the possible mechanism of AST in the treatment of RA from the perspective of lncRNAs. In the present study, the lncRNAs and mRNAs of a vehicle group, animal model group and AST treatment (control) group were determined by Arraystar Rat lncRNA/mRNA microarray. The differentially expressed genes with a fold change >1.5 and P<0.05 were selected and analyzed. Gene Ontology (GO) and pathway analysis was performed using the Database for Annotation, Visualization and Integration Discovery, and the coding‑non‑coding gene co‑expression network was drawn based on the correlation analysis between the differentially expressed lncRNAs and mRNAs. Based on node degree and the correlation between bioinformatics analysis and RA, the critical differentially expressed lncRNAs were selected, analyzed and verified by reverse transcription‑quantitative PCR (RT‑qPCR) analysis. The results showed that, following AST treatment, up to 75 lncRNAs and 247 mRNAs were found to be differentially expressed among the three groups. GO and pathway analysis manifested that 135 GO terms and 17 pathways were enriched by differentially expressed genes. Four lncRNAs (MRAK012530, MRAK132628, MRAK003448 and XR_006457) were selected as the critical lncRNAs and their trend in expression showed consistency between the RT‑qPCR and microarray data. In conclusion, AST had a regulatory effect on differentially expressed lncRNAs during the development of RA, and four lncRNAs could be selected as critical therapeutic targets of AST.

Introduction

Rheumatoid arthritis (RA), a systemic and autoimmune disease, is characterized by chronic inflammation, synovial hyperplasia, joint swelling and tenderness. This disease mainly occurs in individuals aged between 20 and 50 years old and its morbidity rate in women is 2-3 times higher than that in men, particularly in older women. The pathological changes of RA result in functional limitations, working disability and poor quality of life (1-4). One of the important therapeutic approaches for RA in western medicinal therapy is to alleviate pain. This control strategy is given priority mainly for symptomatic treatment. However, this control strategy is not satisfactory in the clinic, and side effects frequently occur (5). Therefore, there is an urgent need to develop novel treatment approaches for RA. Investigations involving traditional Chinese medicine (TCM) for RA have been developing for thousands of years. These ancient therapies remain in use at present as they are empirically effective and have stood the test of time (6,7).

Astragalus membranaceous (Fisch.), also known as Huang Qi, is one of the most commonly used herbal medicines in TCM (8). Astragalus and its compounds have been widely used for the treatment of RA for thousands of years in China (9-11). Astragalosides (AST) are the effective elements extracted from Astragalus, which mainly contains AST-I-VIII. Among these, AST-IV is its main active ingredient. Previous studies have demonstrated that AST has numerous biological effects, including anti-inflammatory, anti-oxidation, immunological regulation and anti-aging effects (12-15).

Accumulating evidence has suggested that, rather than being transcriptional noise, a number of non-coding RNAs (ncRNAs) serve as master regulators that affect the expression levels of dozens or even hundreds of target genes (16,17). Among them, long non-coding RNAs (lncRNAs), a recently discovered class of non-coding RNAs, are mRNA-like transcripts >200 nucleotides in length that have no or little protein-coding capacity and serve important roles in a variety of biological processes, including transcription, splicing and translation (18). In addition, lncRNAs are known to be associated with the pathogenesis of different diseases, including RA (19-21). However, to the best of our knowledge, there has been minimal focus on the differentially expressed lncRNAs in RA following AST treatment.

In the present study, in order to observe the effect of AST on the differentially expressed lncRNAs in adjuvant arthritis (AA) rats, and to evaluate the possible mechanism of AST, the Arraystar Rat lncRNA microarray was used. The results will advance current understanding of the possible mechanism underlying AST in treatment for RA from the perspective of lncRNAs and provide novel lncRNAs as potential treatment targets.

Materials and methods

Animal experiments

The protocol was approved by the Committee on the Ethics of Animal Experiments of Anhui University of Chinese Medicine (Hefei, China). A total of 40 male Sprague Dawley rats (age, 6-8 weeks; SPF-grade; weight, 200±20 g) were obtained from the Experimental Animal Center of Anhui Medical University (Hefei, China). All rats were healthy and had never received drugs previously. The animals were housed in the animal facilities of Anhui University of Chinese Medicine. The animals were housed individually under a set temperature (18-22°C) and humidity (40-60%) with a 12-h light/dark cycle, free access to standard laboratory food and water, and environmental noise kept to a minimum range. All efforts were made to minimize animal suffering and the number of animals used.

Drugs and reagents

AST was purchased from Nanjing Zhi Bai Cui Biology Technology Co. Ltd. (Nanjing, China; purity >98.5%), mainly containing AST-I-VIII. AST-IV was purchased from Nanjing Zelang Medical Technology Co. Ltd. (Nanjing, China, purity >98%). Freund's complete adjuvant (FCA) was purchased from Sigma; Merck KGaA (Darmstadt, Germany). Sodium pentobarbital was obtained from Shanghai Chemical Reagent Company (Shanghai, China). TRIzol reagent was obtained from Invitrogen; Thermo Fisher Scientific, Inc. (Waltham, MA, USA). The C-reactive protein (CRP) kit was obtained from Beijing Lidman Co., Ltd. (Beijing, China). The white blood cell (WBC), red blood cell (RBC), hemoglobin (HGB) and hematocrit (HCT) kits were obtained from Sysmex Corporation (Shenhu, Japan). The ELISA kits of interleukin (IL)-1β (cat. no. 20180622), IL-4 (cat. no. 20180526), IL-10 (cat. no. 20180615) and tumor necrosis factor (TNF)-α (cat. no. 20180702) were obtained from Shanghai Elisa Biotech Co., Ltd. (Shanghai, China).

The rats were adaptively fed for 1 week. All rats were divided randomly into a vehicle group, animal model group, AST treatment (control) group and AST-IV group (n=10). The animal model group and control group were established with a single intracutaneous injection of 0.1 ml of FCA into the right hind metatarsal footpad (22). The vehicle group was administered with liquid paraffin for comparison at the same time. From day 12 after injection of FCA, AST (200 mg/kg) and AST-IV (80 mg/kg) were administered by gavage once a day for 12 days in the control group. The rats in the vehicle and animal model groups were given 0.9% saline.

Hind paw swelling and polyarthritis index

From day 16, swelling of the left hind paw was measured and the poly-arthritis index was assessed once every 4 days until the animals were sacrificed. Swelling of the paw was determined by water volume (23): Swelling (∆ml)=paw thickness (ml) tested value-paw thickness (ml) original value. The severity of arthritis was graded on a scale of 0-4 with the following criteria (24): 0=no edema or swelling, 1=slight edema and limited erythema, 2=slight edema and erythema from the ankle to the tarsal bone, 3=moderate edema and erythema from the ankle to the tarsal bone, and 4=edema and erythema from the ankle to the entire leg. The arthritis score for each rat was the sum of the severity in all four limbs (maximum 16 points per rat).

Synovial tissue samples

On day 28 following immunization, the rats were anesthetized by intraperitoneal injection with sodium pentobarbital (60 mg/kg of body weight). Three rats were randomly selected from each group. Synovial tissue samples were harvested from all rats and stored at −80°C until use.

Histological examination

The ankle joints and synovial tissue were removed and fixed in 10% neutral formalin and 2.5% glutaraldehyde stationary liquid, respectively. Histological changes were subsequently observed by H&E staining and transmission electron microscopy (TEM).

Measurements of CRP, WBC, RBC, HGB and HCT

The content of CRP in serum was measured using an automatic biochemical analytic instrument; the levels of CRP, WBCs, RBCs, HGB and HCT in the whole blood were measured using a complete blood cell counter.

Measurements of IL-1β, TNF-α, IL-4 and IL-10

The serum contents of IL-1β, TNF-α, IL-4 and IL-10 were determined using the respective ELISA kits, according to the manufacturer's instructions.

Microarray and data analysis

RNA extraction

Total RNAs were extracted from the three groups using TRIzol® reagent (Thermo Fisher Scientific, Inc.), according to the manufacturer's protocol. RNA quantity and quality were measured using a NanoDrop ND-1000 instrument. RNA integrity was assessed by 2% standard denaturing agarose gel electrophoresis.

Microarray

The rat lncRNA array was designed for profiling the lncRNAs and mRNAs. Approximately 9,000 lncRNAs were selected from the databases of NCBI RefSeq (https://www.ncbi.nlm.nih.gov/refseq/), UCSC (http://genome.ucsc.edu/), all_mrna records and orthologs of mouse lncRNAs. The probes for coding genes were printed once, whereas the probes for lncRNAs were printed three times. The probes for housekeeping genes and negative probes were printed multiple times to ensure hybridization quality.

RNA labeling and array hybridization

Sample labeling and array hybridization were performed according to the Agilent One-Color Microarray-Based Gene Expression Analysis protocol (Agilent Technology) with minor modifications. Briefly, mRNA was purified from total RNA following the removal of rRNA (mRNA-ONLY™ Eukaryotic mRNA Isolation kit, Epicentre, Illumina, Inc.). Each sample was then amplified in a circulating water bath at 40°C for 2 h with a Quick Amp Labeling kit, One-Color (Agilent Technologies, Inc.; cat. no. 5190-0442), according to the manufacturer's protocol. Subsequently, each sample was transcribed into fluorescent cRNA along the entire length of the transcripts without 3′ bias utilizing a mixture of oligo (dT) and random primers (Arraystar Flash RNA Labeling kit, Arraystar). The labeled cRNAs were purified using an RNeasy Mini kit (Qiagen, Inc.). The concentration and specific activity of the labeled cRNAs (pmol Cy3/µg cRNA) were measured using the NanoDrop ND-1000 instrument. Subsequently, 1 µg of each labeled cRNA was fragmented by adding 5 µl 10X blocking agent and 1 µl 25X fragmentation buffer, heating the mixture at 60°C for 30 min, and finally adding 25 µl 2X GE hybridization buffer to dilute the labeled cRNA. The hybridization solution (50 µl) was dispensed into the gasket slide and assembled to the lncRNA expression microarray slide. The slides were incubated for 17 h at 65°C in an Agilent hybridization oven. The hybrid-ized arrays were washed The hybridized arrays were washed with Gene Expression Wash Buffer 1 (cat. no. 5188-5325; Agilent Technologies) and Gene Expression Wash Buffer 2 (cat. no. 5188-5326; Agilent Technologies), and scanned using the Agilent DNA Microarray Scanner (model G2505C) (25).

Data analysis

Agilent Feature Extraction software (version 11.0.1.1) was used to analyze the acquired array images. Quantitative normalization and subsequent data processing were performed using the GeneSpring GX v12.1 software package (Agilent Technologies). Following quantitative normalization of the raw data, lncRNAs and mRNAs that had at least three of nine samples have flags in Present or Marginal ('All Targets Value') were selected for further data analysis. Differentially expressed lncRNAs and mRNAs were identified through P-value and fold change filtering. Hierarchical clustering was performed using homemade scripts.

Functional group analysis

To assess functional enrichment, Gene Ontology (GO) biological process (BP) term and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses of the differentially expressed mRNAs were performed using the Database for Annotation, Visualization and Integration Discovery (https://david.ncifcrf.gov/; v6.8).

Co-expression network

It is well known that lncRNAs perform their function by positively or negatively regulating gene expression at the transcriptional and post-transcriptional level. Therefore, it was inferred that the function of each lncRNA was based on the function of connected mRNAs. Correlation between the expression of lncRNAs and mRNAs was evaluated using Pearson's correlation coefficient (PCC). The lncRNA-mRNA pairs with PCC ≥0.94 and P<0.05 were selected to construct the co-expression network and was visualized using Cytoscape 3.2.8 software (http://cytoscape.github.io/).

Key lncRNAs

The selection of key lncRNAs was based on two aspects: i) node degree: It is well known that hub nodes serve critical roles in the biological network, therefore all node degrees of the lncRNA-mRNA network were calculated. ii) Correlation between bioinformatics analysis and RA: GO and pathway annotations were performed for those lncRNAs with the top 10 node degrees. According to the node degree and bioinformatics analysis, four lncRNAs were selected as potential key lncRNAs.

Key lncRNA-mRNA sub-networks

Every key lncRNA, its linked lncRNAs and mRNAs in the co-expression network were extracted and used to reconstruct new sub-networks in Cytoscape software. For further analysis, GO and pathway annotations were performed for the key lncRNAs by using their first mRNA neighbors in the key lncRNA-mRNA sub-networks.

Reverse transcription-quantitative PCR (RT-qPCR) analysis

The total RNA was extracted using TRIzol reagent and then cDNA was synthesized with PrimeScript™ RT reagent kit (Takara Bio, Inc.) at 65°C for 5 min and 42°C for 60 min. RT-qPCR analysis was performed using a LightCycler 480 II (Roche, Germany) with TB Green™ Advantage® qPCR Premix (Takara Bio, Inc.). The thermocycling conditions were as follows: Pre-denaturation at 95°C for 5 min, 40 cycles of 95°C for 15 sec and 60°C for 60 sec. The expression of each lncRNA was determined as the fold change using the 2−ΔΔCq method (26). The primers were designed and shown in Table I.

Table I

Primer sequences.

Gene name	Primer sequence (5′-3′)		Product length, bp
	Forward	Reverse
β-actin	5′-CGAGTACAACCTTCTTGCAGC-3′	5′-ACCCATACCCACCATCACAC-3′	202
MRAK012530	5′-TGATAACAGCATGAGACCCAA-3′	5′-AAGCCAAGAACAAGAGCACAA-3′	160
MRAK132628	5′-TTCCAGAGGACCTTCTCAATG-3′	5′-GACAAATCCCACCCTTACCTAC-3′	171
MRAK003448	5′-TTTTGTGAGTAAGGGAGGGTGA-3′	5′-GGGGATTAAGAAGCCATGCTAT-3′	157
XR_006457	5′-CAGCAAGCAGGAGTACGATGA-3′	5′-AGTGACAGTCCGCCTAGAAGC-3′	71

Statistical analysis

All statistical analyses were performed using the SPSS 17.0 statistical software package (SPSS, Inc., Chicago, IL, USA). One-way analysis of variance was used to perform statistical analysis with Tukey's post hoc test for multiple comparisons. Numerical data are presented as the mean ± standard error. P<0.05 was considered to indicate a statistically significant difference.

Results

Secondary inflammation was observed in the right hind paws from day 12, which indicated successful establishment of the AA animal model. The results showed that the hind paw swelling in the AA rats was markedly increased compared with that in the vehicle group and, between days 20 and 28, swelling of the right hind paws was significantly attenuated by AST treatment (Table II). Compared with that in the vehicle rats, the polyarthritis index of the AA rats was significantly increased. Administration with AST resulted in a marked decrease in the polyarthritis index compared with that in the AA rats between days 20 and 28 (Table III).

Table II

Effect of astragalosides on hind paw swelling in adjuvant arthritis rats.

Group	Dose, mg/kg	Day 16, ∆ml	Day 20, ∆ml	Day 24, ∆ml	Day 28, Δml
Vehicle	-	0.203±0.070	0.235±0.077	0.258±0.066	0.288±0.068
Animal model	-	0.573±0.172a	0.750±0.172a	0.883±0.118a	0.782±0.104a
Control	200	0.526±0.138	0.605±0.129b	0.705±0.138c	0.632±0.112c

Analysis of variance and Tukey's test were used to perform statistical analysis in multiple comparisons. Data are expressed as the mean ± standard error (n=10/group).

P<0.01 vs. vehicle group;

P<0.05 and

P<0.01 vs. animal model group.

Table III

Effect of astragalosides on the polyarthritis index in adjuvant arthritis rats.

Group	Dose, mg/kg	Day 16	Day 20	Day 24	Day 28
Vehicle	-	0.40±0.52	0.30±0.67	0.20±0.42	0.30±0.48
Animal model	-	6.90±2.47a	8.50±2.12a	9.80±1.93a	8.00±1.76a
Control	200	5.50±1.58	6.40±2.59b	7.60±2.64b	5.8±1.81c

Analysis of variance and Tukey's test were used to perform statistical analysis in multiple comparisons. Data are expressed as the mean ± standard error (n=10/group).

P<0.01 vs. vehicle group;

P<0.05 and cP<0.01 vs. animal model group.

Levels of CRP, WBCs, RBCs, HGB and HGB

As shown in Table IV, compared with the vehicle group, the levels of CRP and WBCs were significantly elevated in the animal model group. Compared with those in the animal model group, the levels of CRP and WBCs were decreased following treatment with AST. However, there were no significant changes in the levels of RBCs, HGB or HCT among the three groups.

Table IV

Effect of astragalosides on biomarkers of inflammation in adjuvant arthritis rats.

Group	Dose, mg/kg	CRP, ng/l	WBC, 109/l	RBC, 1012/l	HGB, g/l	HCT, %
Vehicle	-	4.99±1.23	8.73±1.44	7.16±0.75	151.48±12.48	44.69±4.65
Animal model	-	13.45±2.58a	12.85±1.47a	7.53±0.66	142.2±12.25	43.04±4.15
Control	200	10.27±2.86b	10.69±1.69b	7.26±0.77	145.23±14.93	44.07±4.89

Analysis of variance and Tukey's test were used to perform statistical analysis in multiple comparisons. Data are expressed as the mean ± standard error (n=10/group).

P<0.01 vs. vehicle group;

P<0.05 vs. animal model group. CRP, C-reactive protein; WBC, white blood cell; RBC, red blood cell; HGB, hemoglobin; HCT, hematocrit.

Levels of IL-1β, TNF-α, IL-4 and IL-10

Compared with the vehicle group, the animal model group exhibited higher levels of IL-1β and TNF-α, and lower levels of IL-4 and IL-10. Following AST intervention, compared with those in the animal model group, the levels of IL-1β and TNF-α were markedly decreased, and those of IL-4 and IL-10 were significantly increased in the control group (Table V).

Table V

Effect of astragalosides on the levels of IL-1β, TNF-α, IL-4 and IL-10 in adjuvant arthritis rats.

Group	Dose, mg/kg	IL-1β, ng/l	TNF-α, ng/l	IL-4, ng/l	IL-10, ng/l
Vehicle	-	35.39±10.68	162.65±25.4	94.23±15.11	197.66±21.1
Animal model	-	73.99±13.18a	284.2±41.28a	61.11±13.51a	125.8±28.96a
Control	200	57.90±12.05b	223.78±36.47c	78.47±9.38b	157.39±14.14b

Analysis of variance and Tukey's test were used to perform statistical analysis in multiple comparisons. Data are expressed as the mean ± standard error (n=10/group).

P<0.01 vs. vehicle group;

P<0.05 and

P<0.01 vs. animal model group. IL, interleukin; TNF, tumor necrosis factor.

Histopathology in rats with AA

H&E staining of samples from the vehicle group showed normal arrangement of syno-vial cells, without synovial hyperplasia or inflammatory cell infiltration. However, samples from the animal model group exhibited local thickening of the surrounding synovium tissue and inflammatory cell infiltration accompanied by spalling pannus formation. In the control group, less pathological tissue injury was present compared with that in the animal model group (Fig. 1A).

Figure 1

Effect of AST on the histopathology of ankle joints and synovial tissues in the vehicle, animal model and control groups. (A) Effect of AST on the histopathology of ankle joints by H&E staining (magnification, ×200). TEM images show the effects of AST on the histopathology of (B) type A synoviocytes and (C) type B synoviocytes in synovial tissues, as shown by TEM (magnification, ×20,000). a, normal synoviocytes; b, articular cavity; c, incrassation of synoviocytes; d, formation of pannus tissue; e, inflammatory cell infiltration; f, nucleus; g, mitochondrion; h, lysosome; i, golgi complex; j, mitochondrial swelling; k, alleviation of swollen mitochondria; l, rough endoplasmic reticulum expansion failure; m, mitochondrial vacuolar and degeneration; n, rough endoplasmic reticulum. AST, astragalosides; TEM, transmission electron microscopy.

TEM revealed that two basic types of synoviocyte in the vehicle group. Type A synoviocytes were similar to macrophages with numerous mitochondria, which possessed the main characteristics of the phase of phagocytosis with multitudinous vacuoles, lysosomes, parts of rough endoplasmic reticulum and ribosomes dispersed in the cytoplasm. Type B synoviocytes were similar to fibroblasts with the main characteristics of the phase of protein secretion. Chromatin was loose and was distributed around the nucleus with a fine granular appearance. As shown in Fig. 1B and C, the vehicle group exhibited abundant rough endoplasmic reticulum, which were ranged in stratiform and without vacuoles. In the animal model group, mitochondria and vacuoles were significantly increased, and mitochondria were swollen in the type A synoviocytes; rough endoplasmic reticulum and dense bodies were significantly increased. In addition, the endoplasmic reticulum cistern was swollen and there was destruction of the structure in the type B synoviocytes. Following AST intervention, the vacuoles were reduced and there was alleviation of mitochondrial swelling in the type A synoviocytes. Furthermore, the rough endoplasmic reticulum was reduced and chromatin in the cell nucleus exhibited pyknosis. The volume of the cell was also reduced in the type B synoviocytes.

Analysis of differentially expressed lncRNAs

Differentially expressed lncRNAs were distinguished among the vehicle, animal model and control groups according to the filtering criteria (P<0.05, fold change >1.5). It was found that there were 260 significantly differentially expressed lncRNAs between the animal model group and the vehicle group, and the expression levels of 915 lncRNAs were significantly altered in the control group compared with those in the animal model group. There were 75 differentially expressed lncRNAs, which were co-expressed differentially in the three groups. A hierarchical clustering analysis was applied to the differentially expressed lncRNAs based on their expression levels among samples, the results are shown in Fig. 2A.

Figure 2

Hierarchical cluster analysis of differentially expressed lncRNAs and mRNAs among the vehicle, animal model and control groups. (A) Hierarchical cluster analysis of 75 differentially expressed lncRNAs. (B) Hierarchical cluster analysis of 247 differentially expressed mRNAs. Hierarchical cluster analysis showed that the differentially expressed genes ultimately clustered into two major branches, including upregulated genes (labeled in red) and downregulated genes (labeled in green). V, vehicle group, M, animal model group, C, control group; lncRNAs, long non-coding RNAs.

Analysis of differentially expressed mRNAs

According to the filtering criteria P<0.05 and fold change >1.5, 675 differentially expressed mRNAs were found between the vehicle and AA animal model groups, and the expression of 1,717 mRNAs were significantly altered between the animal model and control groups. Subsequently, 247 differentially expressed mRNAs were found among the three groups. Their distinct expression patterns are also shown in Fig. 2B through hierarchical clustering analysis.

GO and pathway analysis

A total of 135 GO terms were found to be significantly altered, including acetyl-CoA metabolic process, cofactor metabolic process and coenzyme metabolic process. The top 10 significantly regulated BPs are shown in Fig. 3A. Pathway analysis demonstrated that 17 pathways were significantly enriched among the upregulated and downregulated genes, which were mainly associated with the TCA cycle, pyruvate metabolism, and oxidative phosphorylation. The top 10 significant pathways are shown in Fig. 3B. According to the functional annotation and fold change, 15 upregulated and five downregulated mRNAs were selected (Table VI) to clarify the functions of the differentially expressed lncRNAs. The results revealed that the upregulated genes were mainly associated with regulation of the inflammatory response, the immune response and various metabolic processeses, and the downregulated genes were mainly associated with various types of signaling pathways.

Figure 3

BP and pathway analysis of differentially expressed genes. (A) Top 10 significant GO terms in BPs. (B) Top 10 significant pathways. GO, Gene Ontology, BP, biological process; DE, differentially expressed.

Table VI

Functional annotation of 20 differentially expressed genes.

Number	Gene	Trend	Fold change	P-value	Gene annotation
1	Dio3	Up	5.77	9.93×10−4	Oxidation-reduction process
2	Mlph	Up	7.31	5.38×10−5	Intracellular protein transport
3	Tm6sf2	Up	2.99	2.47×10−3	Regulation of lipid metabolic process
4	Gys2	Up	6.44	5.75×10−3	Glycogen biosynthetic process, response to glucose
5	Mogat2	Up	7.34	0.011	Lipid metabolic process
6	Hgd	Up	5.32	2.14×10−3	Cellular amino acid metabolic process, oxidation-reduction process
7	Dlat	Up	3.96	1.73×10−3	Acetyl-coa biosynthetic process from pyruvate, glucose metabolic
8	Fmo1	Up	1.95	0.013	Organic acid metabolic process, NADPH oxidation, drug metabolic
9	Gcgr	Up	6.81	8.30×10−5	G-protein coupled receptor signaling pathway, regulation of glycogen metabolic process
10	Me1	Up	5.71	2.33×10−3	Oxidation-reduction process, regulation of NADP metabolic process
11	L2hgdh	Up	2.66	0.021	Cellular protein metabolic process, oxidation-reduction process
12	Pdhb	Up	2.38	0.024	Glucose metabolic process, acetyl-coa biosynthetic process
13	Aco1	Up	2.71	2.55×10−3	Tricarboxylic acid cycle
14	Acly	Up	5.43	1.03×10−4	Fatty acid biosynthetic process
15	Acy1	Up	2.02	3.65×10−3	Cellular amino acid metabolic process
16	Angptl4 D	own	4.76	4.25×10−5	Cell differentiation, negative regulation of apoptotic process
17	C1qtnf3 D	own	3.20	3.81×10−4	Fat cell differentiation, positive regulation of cytokine secretion
18	Dhx58 D	own	1.82	4.27×10−4	Immune system process, regulation of innate immune response
19	Cblc	Down	2.62	0.017	Negative regulation of MAP kinase activity, cell surface receptor signaling pathway
20	Npr3 D	own	5.60	2.13×10−3	Osteoclast proliferation, phospholipase C-activating G-protein coupled receptor signaling pathway

lncRNA-mRNA co-expression network

The lncRNA-mRNA co-expression network among the vehicle, animal model and control groups with differentially expressed lncRNAs and mRNAs, which had PCC values ≥0.94, were constructed using the Cytoscape program and then visualized. As shown in Fig. 4, the lncRNA-mRNA co-expression network was composed of 58 lncRNA nodes, 229 mRNA nodes and 1,617 edges. Based on the node degree and correlation between bioinformatics analysis and RA, four significant lncRNAs were selected, two of which were upregulated (MRAK012530 and MRAK132628) and two of which were downregulated (MRAK003448 and XR_006457) in the animal model group.

Figure 4

Co-expression network of differentially expressed lncRNAs and mRNAs following AST treatment. Squares represent lncRNAs, circles represent mRNAs. Upregulated genes are labeled in red, downregulated genes are labeled in green. lncRNAs, long non-coding RNAs.

Key lncRNAs-mRNA sub-networks

The new key lncRNA-mRNA sub-networks were constructed by extracting the key lncRNAs and their linked mRNAs from the co-expression network. As shown in Fig. 5, there were six lncRNA nodes, 13 mRNA nodes and 19 edges (Fig. 5A-a), seven lncRNA nodes, 14 mRNA nodes and 21 edges (Fig. 5B-a), eight lncRNA nodes, 30 mRNA nodes and 38 edges (Fig. 5C-a), and 12 lncRNA nodes, 39 mRNA nodes and 51 edges (Fig. 5D-a) in the lncRNA MRAK012530, MRAK132628, MRAK003448 and XR_006457-mRNA sub-networks, respectively. The top 10 significantly regulated BP and pathways for the key lncRNAs, determined using their first mRNAs neighbors in the key lncRNA-mRNA sub-networks, are shown in Fig. 5A-b and -c, B-b and -c, C-b and -c and D-b and -c.

Figure 5

Sub-network of the key lncRNAs and informatics analysis. (A) lncRNA MRAK012530 (a) key lncRNA-mRNA sub-network, (b) top 10 significant GO biological processes and (c) top 10 significant pathways. (B) lncRNA MRAK132628 (a) key lncRNA-mRNA sub-network, (b) top 10 significant GO biological processes and (c) top 10 significant pathways. (C) lncRNA MRAK003448 (a) key lncRNA-mRNA sub-network, (b) top 10 significant GO biological processes and (c) top 10 significant pathways. (D) lncRNA XR_006457 (a) key lncRNA-mRNA sub-network, (b) top 10 significant GO biological processes and (c) top 10 significant pathways. In the sub-networks, squares represent lncRNAs, circles represent mRNAs. Upregulated genes are labeled in red, downregulated genes are labeled in green. lncRNAs, long non-coding RNAs; GO, Gene Ontology.

RT-qPCR validation

In order to verify the reliability of the microarray data, a number of notable candidate lncRNAs (MRAK012530, MRAK132628, MRAK003448 and XR_006457) were selected and analyzed using RT-qPCR analysis. As shown in Fig. 6, the results demonstrated that lncRNAs MRAK012530 and MRAK132628 were upregu-lated and lncRNAs MRAK003448 and XR_006457 were downregulated in the animal model relative to their levels in the vehicle group. Following AST and AST-IV intervention, the expression levels of lncRNAs MRAK012530 and MRAK132628 were significantly downregulated and those of lncRNAs MRAK003448 and XR_006457 were significantly upregulated. The RT-qPCR results were consistent with those of the microarray analysis, indicating that the false positive results of the gene microarray caused by individual differences and tissue heterogeneity could be excluded, that these four key lncRNAs may be used as potential critical therapeutic targets of AST, and that AST-IV was the main active ingredient of AST, which serves a major role in regulating differentially expressed lncRNAs.

Figure 6

Validation of differentially expressed long non-coding RNA by reverse transcription-quantitative PCR. ##P<0.01 between the animal model group and vehicle group; **P<0.01 between the animal model group and AST group or AST-IV group. AST, astragaloside.

Discussion

RA is a progressive and debilitating autoimmune disease (27). Synovial tissue is principally composed of fibroblast-like synovial cells, which are closely associated with the pathogenesis and development of RA, particularly joint damage (28). lncRNAs, a recently identified type of ncRNA, are involved in numerous biological processes, including cellular differentiation, cell growth and cell cycle regulation (29). Evidence from previous research indicates the aberrant expression of lncRNAs in the pathogenesis of RA (30,31). In our previous study, it was also confirmed that up to 260 lncRNAs (170 upregulated and 90 downregulated lncRNAs) were found to be differentially expressed in the synovium between an AA animal model and normal rats (25). However, fewer studies have examined the expression profile of lncRNAs in RA following AST treatment. In the present study, 75 differentially expressed lncRNAs and 247 differentially expressed mRNAs were found in synovial tissues of the vehicle, animal model and control groups. Using the node degree and correlation between bioinformatics analysis and RA, the results suggested that four lncRNAs (MRAK012530, MRAK132628, MRAK003448 and XR_006457) could be selected as critical therapeutic targets of AST. In order to exclude the false positive results of the gene microarray caused by individual differences and tissue heterogeneity, RT-qPCR was used to verify the four lncRNAs. The results showed consistency between the RT-qPCR and microarray data, which indicated that these four lncRNAs have increased reliability and repeatability and can be used as potential treatment targets of AST.

To investigate the potential regulatory roles of lncRNAs, GO and pathway analyses were used to reveal the biological functions enriched among the differentially expressed coding-genes (32). In the present study, the results of the GO analysis revealed that differentially expressed genes were mainly involved in acetyl-CoA metabolic process, cofactor metabolic process and coenzyme metabolic process. The pathway analysis showed that 17 pathways were enriched and primarily involved in the citrate cycle (TCA cycle), glyoxylate and dicarboxylate metabolism, pyruvate metabolism, oxidative phosphorylation and PPAR signaling pathway, which may have essential roles in the occurrence and development of RA (33-37).

Although thousands of lncRNAs have been identified and recorded in several public databases, the functional characterization of lncRNAs remains in its infancy. Until now, only a small number of lncRNAs have been well functionally annotated (38,39). Improved knowledge has suggested that lncRNAs, owing to their increased length, can regulate microRNA (miRNA) abundance by binding and sequestering them, acting as so-called miRNA sponges. In this manner, lncRNAs can regulate the expression of target mRNAs (40). Therefore, it has been shown that an efficient way to infer the potential function of lncRNAs is by examining related mRNAs, whose functions have been annotated (41-43). In the present study, according to the fold change and functional annotation, 20 mRNAs were selected to examine the possible mechanism underlying the effect of AST treatment in RA.

Deiodinase, iodothyronine type III (Dio3), is involved in the oxidation-reduction process and positive regulation of multicellular organism growth (44). The oxidation-reduction reaction is an integrated aspect of cell metabolism and dynamic equilibrium. This reaction occurs between oxygen metabolites and certain amino acids, proteins or enzymes, which have important physiological functions of reducing its activity and macromolecular substances. Those products mainly mediate endothelial cells and tissue injury and aggravate the injury of inflammation-stimulating factors (45). In this study, AST not only inhibited the expression of Dio3, but also downregulated several other oxidation-reduction reaction-related genes, including homogentisate 1,2-dioxygenase, flavin containing dimethylaniline monoxygenase 1, malic enzyme 1 and L-2-hydroxyglutarate dehydrogenase, which indicated that AST regulated the oxidation-reduction reaction.

In our previous study, using gas chromatography time-of-flight mass spectrometry and a variety of multivariate statistical methods, including PCA, PLS-DA and OPLS-DA, it was confirmed that L-dopa, 1,4-dihydroxy-2-naphthoic acid and tartronic acid, which are primarily associated with amino acid, glycogen metabolism and the TCA cycle, were markedly altered (46,47). In the present study, it was found that certain genes were significantly altered in the animal model group, including transmembrane 6 superfamily member 2, monoacylglycerol O-acyltransferase, ATP citrate lyase and angiopoietin-like 4, which are related to lipid metabolism (48,49), glycogen synthase 2, dihydrolipoamide S-acetyltransferase, glucagon receptor and pyruvate dehydrogenase E1 beta subunit, which are important in several biological functions including glycogen metabolic process and response to glucose (50,51), and aconitase 1, which is involved in TAC (52,53), indicating there were some metabolic disorders during the development of RA, which is consistent with our previous studies (46,47). Following the administration of AST, the expression of these genes returned to normal levels, indicating that AST can regulate disorders of metabolism.

C1q and tumor necrosis factor-related protein 3 (C1qtnf3) is involved in the negative regulation of interleukin-6 (IL-6) secretion (54). IL-6 is involved in inflammatory immune cell multiplication and aggravates the development of inflammation in RA (55,56). DEXH (Asp-Glu-X-His) box polypeptide 58 (Dhx58) is involved in the RIG-I-like receptor signaling pathway, which can negatively regulate the innate immune response (57). Casitas B-lineage lymphoma c (Cblc) is involved in the JAK-STAT signaling pathway and is an important negative regulator in the T-cell and B-cell antigen receptor signaling pathway (58,59). Several experiments have demonstrated that the T cell receptor signaling pathway and JAK-STAT signaling pathway are involved in RA (60-63). Following AST treatment, the expression levels of C1qtnf3, Dhx58 and Cblc were increased. The levels of CRP, WBCs, IL-1β and TNF-α were decreased and the levels of IL-4 and IL-10 were significantly increased following treatment with AST. These results suggested that AST inhibited inflammation in RA.

The hyperplasia of synoviocytes is one of the major morphologic features of RA. It has been identified that synoviocytes can cause severe joint destruction by secreting various cytokines, inflammatory mediators and proteinases (64,65). Natriuretic peptide receptor 3 (Npr3) negatively regulates cell proliferation (66). In the present study, the level of Npr3 was significantly increased following AST treatment compared with that in the AA animal model group. These results suggested that AST inhibited hyperplasia of the synoviocytes.

In conclusion, the present study advances current understanding of the possible mechanism of AST in the treatment of RA from the perspective of lncRNAs, and provides novel lncRNAs as potential critical therapeutic targets of AST. In addition, AST-IV was the main active ingredient of AST, which serves a major role in regulating differentially expressed lncRNAs.