In silico analysis of serum miRNA profiles in seronegative and seropositive rheumatoid arthritis patients by small RNA sequencing

Ethical statement

All the experimental procedures were approved by the Institutional Review Board of Guangdong Provincial Hospital of Chinese Medicine (#BF2018-088-01) and conducted in accordance with the Declaration of Helsinki. Written informed consent had been obtained from all participants in this study.

Patients

In this study, nine healthy people (HC group), 12 patients with seronegative RA (N-RA group) and nine patients with seropositive RA (P-RA group) were recruited. All patients were recruited from the rheumatology clinic of the Second Affiliated Hospital of Guangzhou University of Chinese Medicine (Guangdong Provincial Hospital of Chinese Medicine, Guangzhou, China). RA patients were enrolled according to the following inclusion criteria: (1) RA was diagnosed according to the 2009 American College of Rheumatology (ACR) diagnostic criteria; (2) patients were classified as grade I, II or III according to the classification of American Rheumatic Association in 1987; (3) patients agreed to participate in the study and signed the informed consent. Then patients were diagnosed as seronegative or seropositive RA based on the presence of RF and ACPA.

Serum collection

Fasting venous blood were drawn and centrifugated to separate serum. Subsequently, the separated serum was stored at −80 °C for detection.

Construction, sequencing and miRNA identification of small RNA library

Total RNA was extracted from serum by TRIzol (Thermo-Fisher-Scientific, Waltham, MA, USA). Next, RNAs with the length of 18–30 nucleotides were enriched utilizing polyacrylamide gel electrophoresis. Subsequently, the 5′ and 3′ end adaptors were ligated to the RNA. Then the ligation products were reverse transcribed, and 140–160 base products were collected to construct small RNA library and sequenced by Illumina HiSeqTM 4000 in Beijing Genomics Institute (Shenzhen, Guangdong, China).

In order to get high-quality clean reads, raw reads were analyzed through in-house Perl scripts. After eliminating the low-quality reads with >10% poly-N sequences and <5% Phred score, all the clean reads were aligned with small RNAs using GeneBank database and Rfam database (11.0) to remove rRNA, scRNA, snoRNA, snRNA and tRNA. At the same time, all clean reads were also aligned with human reference genome (Grch37) by TopHat (v2.0.9) software.

Subsequently, all clean reads were blasted in the miRBase database (version 21) to identify the known miRNAs in the clean reads. Besides, all unannotated tags were aligned with the human reference genome (Grch37) to explore new candidate miRNAs according to the predicted genome location and hairpin structure by Mireap_v0.2 software.

MiRNA profiles

The levels of miRNAs were standardized to transcripts per million (TPM) using the following formula:

$$\mathrm{T}\mathrm{P}\mathrm{M}=\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{u}\mathrm{a}\mathrm{l}\mathrm{m}\mathrm{i}\mathrm{R}\mathrm{N}\mathrm{A}\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{n}\mathrm{t}\mathrm{i}\mathrm{t}\mathrm{y}/\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}\mathrm{c}\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{a}\mathrm{t}\mathrm{a}\mathrm{u}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{g}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{f}\mathrm{o}\mathrm{l}\mathrm{l}\mathrm{o}\mathrm{w}\mathrm{i}\mathrm{n}\mathrm{g}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{m}\mathrm{u}\mathrm{l}\mathrm{a}\times {10}^6$$

In addition, the total miRNA heat map was created to show the miRNA levels in different groups and to cluster miRNAs with similar patterns.

Analysis of differential miRNA

Based on the miRNA levels, the differential miRNAs among different groups were determined by the following formula:

$$p(x|y)={\left(\frac{N_2}{N_1}\right)}^y\frac{(x+y)!}{{x!y!(1+\frac{N_2}{N_1})}^{x+y+1}}\begin{array}{c}C(y\le y_{\mathrm{m}\mathrm{i}\mathrm{n}}|x)=\sum _{y=0}^{y\le y_{min}}p(y|x)\\ D(y\ge y_{\mathrm{m}\mathrm{a}\mathrm{x}}|x)=\sum _{y\ge y_{\mathrm{m}\mathrm{a}\mathrm{x}}}^{\mathrm{\infty }}p(y|x)\end{array}$$

MiRNAs with a fold change ≥2 and P value < 0.05 were considered as significantly differential miRNAs.

Validation of small RNA sequencing data by qPCR

QPCR assays were performed to confirm the reliability of the small RNA sequencing data according to previous studies (Xu et al., 2018; Zhong et al., 2014). First, reverse transcription for miRNAs were performed by the miRcute miRNA First-Strand cDNA Synthesis Kit (Tiangen, Beijing, China). Next, miRNA levels were normalized to the level of 5sRNA using the ΔΔCT method. The sequences of primers used in this study were listed as follows: 5sRNA forward: 5′-AAGCCTACAGCACCCGGTAT-3′, reverse: 5′-GTAACGCCCGATCTCGTCT-3′.

Fold change (log2) values of miRNAs were presented as mean $\pm$ standard error of the mean (SEM). Moreover, statistical analyses were performed on fold change values using ordinary one-way ANOVA, marked as *P < 0.05 compared with other groups.

Prediction of miRNA target genes

Based on the sequences of miRNAs, candidate target genes of miRNAs were predicted by RNAhybrid (v2.1.2), PITA, Miranda (v3.3a) and TargetScan. Furthermore, the intersection of target genes predicted by the above software and the website were more reliable. In addition, the miRNA-target gene network was constructed by using Cytoscape software (v3.6.0).

Enrichment analysis of GO and KEGG pathways for miRNA target genes

Based on the Gene Ontology (http://www.geneontology.org/) database, all the target genes of differential miRNAs were matched to GO terms, and the number of genes matching each GO term was also counted. In addition, GO terms with significant enrichment compared to human genome background were defined via hypergeometric test analysis. Besides, KEGG was used for molecular pathway enrichment analysis to identify pathways in which target genes of differential miRNAs are involved.

Identification of miRNAs in serum

Thirty small RNA libraries including 12 N-RA samples, nine HC samples and nine RA samples were constructed and sequenced to demonstrate the miRNA profiles. Clean reads were obtained after filtering low quality reads. Next, rRNA, scRNA, snoRNA, snRNA and tRNA were removed from the clean reads through comparing with GenBank and Rfam (11.0). Then the clean reads were mapped to the human reference genome (Grch37) by using TopHat. Results indicated that 180, 201 and 216 known miRNAs (Table S1) and 15, 11 and 18 new miRNAs (Table S2) were found in N-RA, HC and P-RA samples, respectively.

Prolife of differential miRNAs in serum

Levels of serum miRNAs were determined by reading count and TPM analysis. Colors from blue to red stand for z-score got through the dimensionality reduction of FPKM value and reveal decreasing miRNA levels in each group. Compared with HC group, five differential miRNAs (one up-regulated and four down-regulated) were found in N-RA group, including hsa-miR-362-5p, hsa-miR-4429, hsa-miR-378e, hsa-miR-302c-3p and hsa-miR-378g (Fig. 1A and Table 1). Compared with P-RA group, there were six differential miRNAs (two up-regulated and four down-regulated) in N-RA group, including hsa-miR-362-5p, hsa-miR-708-3p, hsa-miR-6741-5p, hsa-miR-3127-5p, hsa-miR-6855-5p and hsa-miR-187-3p (Fig. 1B and Table 2). Compared with HC group, seven differential miRNAs (three up-regulated and four down-regulated) were identified in P-RA group, including hsa-miR-6855-5p, hsa-miR-187-3p, hsa-miR-371a-5p, hsa-miR-302a-3p, hsa-miR-320e, hsa-miR-218-5p and hsa-miR-378e (Fig. 1C and Table 3). Among them, the level of hsa-miR-362-5p in N-RA group was upregulated compared with that in HC group and P-RA group, the level of hsa-miR-378e in N-RA group and P-RA group was downregulated compared with that in HC group, and the level of hsa-miR-6855-5p and hsa-miR-187-3p in P-RA group was upregulated compared with those in N-RA group and HC group (P < 0.05). These results suggest that serum hsa-miR-362-5p might be used as a diagnostic marker for seronegative RA, while serum hsa-miR-6855-5p and hsa-miR-187-3p might be promising diagnostic markers for seropositive RA.

Validation of small RNA sequencing data

Subsequently, levels of serum hsa-miR-362-5p, hsa-miR-6855-5p and hsa-miR-187-3p were validated in 11 HC samples, 10 N-RA samples and 10 P-RA samples by qPCR. There was a 165.3% increase in serum hsa-miR-362-5p level in N-RA group compared to that in HC group, while there was a 178.4% increase in serum hsa-miR-362-5p level in N-RA compared to that in P-RA group (Fig. 2A). Besides, there was a 575.6% increase in serum hsa-miR-6855-5p level in N-RA compared to that in HC group, while there was a 585.1% increase in serum hsa-miR-6855-5p level in N-RA compared to that in P-RA group (Fig. 2B). Moreover, there was no difference of serum hsa-miR-187-3p level among three groups (Fig. 2C). Above results showed that the level of serum hsa-miR-362-5p was consistent with that determined by small RNA sequencing (Fig. 2A). However, levels of serum hsa-miR-6855-5p and hsa-miR-187-3p were opposite to results of small RNA sequencing (Figs. 2B and 2C). Thus, these data suggested that serum hsa-miR-362-5p should be the promising diagnostic marker for seronegative RA.

Prediction of differential miRNA target genes

MiRNA often plays a vital role in biological process by regulating the expression of target genes. To explore the effects of miRNAs on seronegative and seropositive RA, target genes of differential miRNAs were predicted to understand potential functions of differential miRNAs. Results showed that there were 4,933 recognized target genes of differential miRNAs between N-RA and HC group (Fig. 3A), 7,496 recognized target genes of differential miRNAs between N-RA and P-RA group (Fig. 3B), and 4,778 recognized target genes of differential miRNAs between P-RA and HC group (Fig. 3C). Among them, target genes of differential miRNAs between N-RA and HC groups included dysferlin (DYSF), myosin light chain 2 (MYL2), glucosylceramidase beta (GBA), tumor necrosis factor (TNF), CD4 molecule (CD4), mannosidase endo-alpha like (MANEAL), triggering receptor expressed on myeloid cells like 2 (TREML2) (Fig.4A and Table S3). Besides, target genes of differential miRNAs between N-RA and P-RA group included phosphatase, orphan 2 (PHOSPHO2), GBA, TNF, CD79b molecule (CD79B), CD4, calcium and integrin binding family member 2 (CIB2), progestin and adipoQ receptor family member 6 (PAQR6) (Fig.4B and Table S3). Moreover, target genes of differential miRNAs between P-RA and HC group included steroidogenic acute regulatory protein (STAR), MYL2, GBA, TNF, CD4, mitochondrial ribosomal protein L11 (MRPL11), oxytocin receptor (OXTR) (Fig.4C and Table S3). Results showed that most target genes of differential miRNAs among three groups were consistent. In addition, miRNA-mRNA networks revealed that target genes could be regulated by both upregulated miRNAs and downregulated miRNAs (Figs. 4A–4C and Table S3).

Functional analysis of differential miRNAs

To further clarify the function of differential miRNAs, GO and KEGG pathway enrichment were utilized to analyze the physiological processes and molecular pathways contributed to the regulation of target genes. GO enrichment analysis showed that target genes of differential miRNAs between N-RA and HC group were mainly involved in binding, catalytic activity, molecular transduction activity, transport activity, receptor activity (Fig. 5A). Besides, target genes of differential miRNAs between N-RA and P-RA group mainly contributed to catalytic activity, molecular transduction activity, receptor activity, transport activity and nucleic acid transcription factor activity (Fig. 5B). Furthermore, target genes of differential miRNAs between P-RA and HC group were mainly involved in binding, catalytic activity, molecular transduction activity, transport activity and nucleic acid transcription factor activity (Fig. 5C). These results indicated that molecular functions of differential miRNAs among three groups were almost identical.

In addition, KEGG enrichment analysis showed that target genes of differential miRNAs between N-RA and HC group were mainly involved in VEGF signaling pathway, natural killer cell mediated cytotoxicity, leukocyte endothelial migration, N-glycan biosynthesis, sphingolipid metabolism, and antigen processing presentation (Fig. 6A). Besides, target genes of differential miRNAs between N-RA and P-RA group were mainly involved in vitamin B6 metabolism, natural killer cell mediated cytotoxicity, N-glycan biosynthesis, B-cell receptor signaling pathway, and antigen processing and presentation (Fig. 6B), while target genes of differential miRNAs between P-RA and HC group mainly contributed to transcriptional misregulation in cancer, vitamin B6 metabolism, natural killer cell mediated cytotoxicity, sphingolipid metabolism, and N-glycan biosynthesis (Fig. 6C). These results revealed that differential miRNAs among the three groups regulated similar molecular pathways.