Differential expression and correlation of immunoregulation related piRNA in rheumatoid arthritis

Abstract

Background:

Although PIWI-interacting RNAs (piRNAs) have recently been associated with germline development and many human diseases, their expression pattern and relationship in autoimmune diseases remain indistinct. This study aimed to investigate the presence and correlation of piRNAs in rheumatoid arthritis (RA).

Methods:

We first analyzed the expression profile of piRNAs using small RNA sequencing in peripheral leukocytes of three new-onset untreated RA patients and three healthy controls (HCs). We then selected piRNAs related to immunoregulation by bioinformatics analysis and verified them in 42 new-onset RA patients and 81 HCs by RT-qPCR. Furthermore, a receiver operating characteristic curve was generated to quantify the diagnostic performance of these piRNAs. A correlation analysis was conducted to observe the link between piRNA expression and RA clinical characteristics.

Results:

A total of 15 upregulated and 9 downregulated piRNAs among 1,565 known piRNAs were identified in peripheral leukocytes of RA patients. Dysregulated piRNAs were enriched in numerous pathways related to immunity. After selection and validation, two immunoregulation piRNAs (piR-hsa-27620 and piR-hsa-27124) were significantly elevated in RA patients and have good abilities to distinguish patients from controls, which have the potential to serve as biomarkers. PIWI and other proteins implicated in the piRNA pathway were also associated with RA.

1 Introduction

Rheumatoid arthritis (RA) is a chronic inflammatory autoimmune disease that will bring about accumulating joint damage and irreversible disability (1); therefore, timely and precise diagnosis is of high importance in RA treatment to prevent adverse clinical endpoints. Typically, RA is often diagnosed by elevated anti-cyclic citrullinated peptide antibodies (ACPAs), rheumatoid factor (RF), and other specific autoantibody markers (2). Although the specific pathogenesis remains unclear, it is believed that RA has a strong genetic component (3). For this reason, the identification of novel genetic biomarkers and the integration of multiple indicators may provide better diagnostic insights for clinicians.

Recently, a great number of studies have shown that non-coding RNA (ncRNA), as an important gene regulatory factor, participates in the occurrence and development of RA through diverse pathways (4). PIWI-interacting RNA (piRNA, 24–32 nt), compared with microRNA (miRNA), long non-coding RNA (lncRNA), etc., is a relatively unexplored ncRNA that is characterized by specific interaction with PIWI protein, a member of the Argonaute family, to form an RNA-induced silencing complex to protect genome integrity by suppressing transposable elements (TEs) (5). In fact, in addition to negative regulation, the piRNA/PIWI system can also stabilize the expression of target genes and induce translational activation according to different binding sites (6). Diversified functional regulation modes also show great research potential. Since piRNA was initially mainly found in mammalian germ lines and embryonic stem cells (7, 8), most of the early studies on piRNA were focused on this (9, 10). Nowadays, with the expansion of the research scope, piRNA has been gradually found in the heart, islet, nervous system, and other somatic cells (11–13) and is closely associated with various diseases, such as breast cancer, hepatocellular carcinoma, bladder cancer, multiple myeloma, and so on (14–20). Elena et al. identified that many piRNAs were expressed differentially after lineage commitment towards osteo- and chondrogenesis of human bone marrow mesenchymal stromal cells (hMSCs) (21). Zhong et al. found that piR30840 interacts with Piwil4 and Ago4 and targets the intron of the IL-4 gene, resulting in decreased IL-4 secretion from CD4+T cells and inhibiting Th2 cell differentiation (22). Furthermore, similar to miRNA, piRNA can be detected in easily accessible samples (such as serum, plasma, saliva, and urine) (23) and remain stable regardless of repeated freeze-thawing or long-term incubation, which has the potential to be a disease biomarker (24). Pleštilová et al. reported that the piRNA/PIWI system exists in synovial tissues and synovial fibroblasts of RA patients and might be involved in signaling pathways for host defense and inflammation (25). However, whether piRNAs are indeed associated with RA remains unclear. We also do not know if piRNA dysregulation occurs in patients’ peripheral blood (which is an easier sample to obtain than synovial tissue), except for synovial tissue.

Accordingly, this study aims to reveal the expression of piRNA/PIWI in peripheral leukocytes of RA patients and their correlation with clinical indexes, and preliminarily explore potential mechanisms based on bioinformatics analysis, which will provide preliminary research data for further exploration of the role of the piRNA/PIWI system in the RA pathogenesis.

2 Materials and methods

2.1 Patients and samples

All new-onset untreated RA patients and healthy controls (HCs) were recruited from the West China Hospital of Sichuan University between October 2020 and March 2023. The inclusion criteria for new-onset untreated patients were fulfilled with the ACR/EULAR 2010 Classification Criteria for Rheumatoid Arthritis (26) and diagnosed at West China Hospital. The exclusion criteria were as follows: treated with immunosuppressive agents, glucocorticoids, and DMARDs before diagnosis; long-term smoking or passive smoking history; and suffering from other autoimmune diseases, acute co-infections, cancers, and other systemic chronic diseases simultaneously. During the same period, gender- and age-matched healthy volunteers were enrolled in the control group. All blood samples were collected before drug treatment. DAS28-CRP, anti-cyclic citrullinated peptide antibody (ACPA), C-reactive protein (CRP), rheumatoid factor (RF), and many other clinical data were collected when patients were recruited. This study was approved by the Ethics Committee of the West China Hospital of Sichuan University.

2.2 RNA extraction

Total RNA was extracted from peripheral leukocytes using QIAzol^® Lysis Reagent (Qiagen, Hilden, Germany) and purified by the miRNeasy mini-Kit (Qiagen, Hilden, Germany). The concentration and purity of the RNA were assessed by a NanoVue Plus microspectrophotometer (GE Healthcare, Madison, WI, USA).

2.3 Small RNA sequencing and bioinformatic analysis

Three RA patients and three HCs were selected for detecting small RNAs with small RNA sequencing. The RNA quality was assured with an Agilent 2100 chip (Agilent Technologies, USA). Small RNA library preparations and sequencing were conducted at Shanghai Oebiotech Co. (Shanghai, China).

The basic reads were converted into raw data by base calling. Low-quality reads were filtered, and the reads with 5’ primer contaminants and poly (A) were removed. The reads without the 3’ adapter and insert tag, the reads shorter than 15 nt and longer than 41 nt from the raw data, were filtered, and the clean reads were obtained. The clean reads from other known types of small RNAs (rRNAs, tRNAs, snRNAs, snoRNAs, and miRNAs) were removed by aligning them to Rfam v.10.1 (http://www.sanger.ac.uk/software/Rfam) (27), miRBase databases (https://www.mirbase.org/) (28), and GenBank databases (http://www.ncbi.nlm.nih.gov/genbank/). Through the above several methods to filter, obtained clean reads were aligned with the known piRNA from the piRBase (http://www.regulatoryrna.org/database/piRNA/) (29) with Bowtie software without mismatches.

The piRNA expression levels of each sample were normalized using TPM. Differentially expressed (DE) piRNAs were identified with the threshold of fold change >2 and P <0.05, calculated with the DEG algorithm (30) in the R package for experiments with biological replicates. The targets of DE piRNAs were predicted by Miranda software and piRNAQuest V.2 (31), with the parameters as follows: S ≥ 150, ΔG ≤ −30 kcal/mol. Gene ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of DE piRNA-target- gene were respectively performed using R based on the hypergeometric distribution.

2.4 Analysis of gene profile targeted by immunoregulation-related piRNA

The GEO database was used to obtain transcriptome sequencing data. Screening criteria included the following: (1) Homo sapiens expression profiling by high throughput sequencing; (2) peripheral blood of RA patients or synovial tissue from joint synovial biopsies; (3) datasets contain no less than five samples. Finally, two datasets, GSE193193 and GSE 77298, were selected for further analysis (Table Supplementary 1). The statistical programming language R (version 4.0.3) (http://www.r-project.org/) was used for log2 transformation of the data. The “SVA” package was used for batch correction. The R package “limma” was used to identify differentially expressed genes (DEGs), which were defined as the absolute value of log2 fold change |log2FC| >1 and adjusted P <0.05. To better visualize these genes, the pheatmap package was used to make heatmaps.

2.5 Quantitative reverse transcriptase-polymerase chain reaction (RT-qPCR)

To verify the expression levels of immunoregulation-related piRNA identified from RA patients by RNA sequencing, five upregulated piRNAs (piR-hsa-27620, piR-hsa-27124, piR-hsa-23940, piR-hsa-27400, and piR-hsa-30353) and two downregulated piRNAs (piR-hsa-25672 and piR-hsa-4946) were selected and validated with the All-in-One™ miRNA qRT-PCR Detection Kit (GeneCopoeia, Rockville, MD, United States). All-In-One 5× RT MasterMix and BlasTaq™ 2× qPCR MasterMix (ABM, Vancouver, BC, Canada) were used to estimate the expression of PIWI and other proteins implicated in the piRNA pathway.

Template RNA was reverse transcribed into cDNA in a Veriti 96-Well Thermal Cycler (Applied Biosystems, Forster City, CA, USA) and measured by an ultraviolet spectrophotometer (Eppendorf, Hamburg, Germany). Real-time PCR was performed on the QuantStudio^®Q5 Real-Time PCR system (Applied Biosystems). NRT (no reverse transcriptase) and NTC (no template control) served as negative controls. U6 and GAPDH served as the reference genes. The relative expression level was calculated by the 2^−ΔΔCT method. All primers (Table S2) were synthesized by Sangon Biotech (Shanghai) Co., Ltd.

2.6 Statistical analysis

All data were analyzed using SPSS (version 26.0; SPSS Inc., Chicago, IL, USA), SAS (version 9.4; SAS Institute, Cary, NC), and GraphPad Prism v9.0 (GraphPad Software, San Diego, CA, USA) software. When the clinical data were lower than the limit of detection, the detection limit was used instead of the true value. For quantitative data, the Shapiro–Wilk normality test was used for normality analysis. The normal distribution data was presented by mean ± standard deviation and analyzed by the Student-t test. Non-normally distributed data were represented by the median (interquartile spacing) [M(IQR)] and analyzed by the Mann–Whitney test. Categorical data were evaluated using the Chi-square test. Logistic regression analysis and receiver operating characteristic (ROC) curves were generated to determine the value of piRNAs as diagnostic biomarkers. Spearman’s correlation (r_s) was used to analyze the correlation of candidate piRNAs with the levels of clinical data. A difference of P <0.05 was considered statistically significant. G*Power was used to estimate the sample size and compute statistical power (32). A power greater than 0.8 was considered acceptable.

3 Results

3.1 Characteristics of the study population

The clinical information of the studied patients is presented in Table 1. There were no differences in age (P = 0.136) or gender (P = 0.360). The RA laboratory diagnostic index, ACPA, RF, and AKA of patients were significantly higher than those of healthy individuals (all P <0.001). Patients had elevated CRP and C3 levels, indicating active inflammation. Furthermore, the blood routine of RA patients showed significantly different patterns compared with that of healthy individuals.

3.2 The piRNA dysregulation profile was revealed by small RNA sequencing in peripheral leukocytes of RA patients

To comprehend the patterns of piRNA expression, small RNA sequencing was used to explore the small RNA expression profiles in the peripheral leukocytes of healthy people and RA patients. Sequence alignment annotation results of six blood samples showed that a total of 1,565 known piRNAs were identified. The length distribution displayed two different peaks, one at the 21–23 nt position, corresponding to miRNA in size, and another concentrated in the 26–28 nt position, corresponding to piRNA in size (Figure 1A). Furthermore, to explore the differential expression of piRNA in RA patients and healthy controls, we conducted expression analysis with the DEG algorithm. The box– whisker plot indicated no abnormal expression in the six samples (Figure S1). Then, after the piRNA expression levels of the two groups were standardized, we used fold change ≥ 2 and P < 0.05 as the standard to screen DE piRNAs. The overall picture of piRNA expression between the two groups was displayed in the volcano plot (Figure 1B). Hierarchical clustering revealed the expression of DE piRNAs in each sample (Figure 1C). A total of 24 DE piRNAs were recognized, including nine downregulated and 15 upregulated piRNAs (Table S3).

Figure 1

3.3 GO term and KEGG pathway enrichment analysis of DE piRNAs

Target genes of the DE piRNAs were calculated by the Miranda algorithm and piRNAQuest V.2, and 48,908 genes were predicted. GO and KEGG pathway analyses were performed on target genes to investigate the role of differentially expressed piRNA in RA pathogenesis.

Based on the enriched dysregulated piRNA targets originated from the gene annotation, the top 10 GO classifications of biological process (BP), cell components (CC), and molecular function (MF) were shown according to routine GO classification algorithms (Figure 1D). In BP, target genes of DE piRNAs were significantly enriched in transcription, regulation of transcription, signal transduction, regulation of transcription from the RNA polymerase II promoter, cell differentiation, the G- protein-coupled receptor signaling pathway, and cell adhesion. In CC, they were mainly concentrated on the cytosol, nucleus, cytoplasm, nucleoplasm, extracellular exosome, and mitochondrion. In MF, genes were primarily related to metal ion binding, ATP binding, DNA binding, RNA binding, identical protein binding, zinc ion binding, and so on. The top 10 GO processes of upregulated/downregulated piRNA targets were also enriched in the regulation of transcription, signal transduction, cell differentiation, and cell adhesion in the BP subgroup (Figures S2A, B).

KEGG pathway analysis assists us in better understanding the biological function of genes. Figure 1E illustrates the top 20 ranked significant pathways for DE piRNA targets in KEGG. The targets of upregulated piRNAs were primarily enriched in the PI3K–Akt signaling pathway, Rap1 signaling pathway, cAMP signaling pathway, axon guidance, regulation of the actin cytoskeleton, calcium signaling pathway, proteoglycans in cancer, and so on (Figure S2C). Downregulated piRNA targets were also enriched in many important pathways, such as the PI3K– Akt signaling pathway, MAPK signaling pathway, Th17 cell differentiation, cGMP–PKG signaling pathway, and ECM–receptor interaction (Figure S2D).

Subsequently, to screen piRNAs that may be related to immune regulation, GO and KEGG enrichment analyses of target genes of the 24 DE piRNAs were performed, and the target genes enriched in the immune system, immune signaling pathways, and immune diseases were screened. Based on the target gene enrichment analysis, piRNAs that regulate those target genes were screened. A total of seven piRNAs highly correlated with RA immune regulation were selected for subsequent RT-qPCR verification, and the screening results are shown in Table S4.

3.4 The expression of gene targets of immunoregulation-related piRNAs

To better understand the expression of genes that were potentially regulated by immunoregulation-related piRNA, two datasets (GSE193193 and GSE77298) were selected from the GEO database to analyze their expression in the peripheral blood and synovial tissue of RA patients, respectively. GSE193193, including five peripheral blood mononuclear cells (PBMCs) from RA patients and five PBMCs from HCs, was used to analyze the expression profile of genes targeted by immunoregulation-related piRNAs in the peripheral blood of RA patients (Figure 1F, Table S5). KLRD1, TBX21, and IL23R, the potential targets of piR-hsa-25672, piR-hsa-27620, and piR-hsa-23940, respectively, were all significantly downregulated, while the level of C4B, the potential target of piR-hsa-25672, was significantly higher in RA patients than that in HCs. The GSE77298 dataset was used to explore these gene profiles in synovial tissue (Figure 1G). Differential expression analysis showed that 18 genes (STAT1, TNFSF13, MMP3, FCGR2A, CCL5, etc.) were significantly differentially expressed in synovial tissue (Table S6).

3.5 Validation of immunoregulation-related piRNAs with RT-qPCR

RT-qPCR was performed to validate the expression levels of seven candidate piRNAs in 42 new-onset RA patients and 81 HCs (Table 2). According to the qPCR results, the levels of two piRNAs [piR-hsa-27620 (P < 0.001, power = 0.98) and piR-hsa-27124 (P < 0.001, power = 0.95)] were significantly higher in RA patients than those in HCs. Nevertheless, no significant difference was observed in the other five piRNAs (Figures 2A–G).

Figure 2

3.6 Evaluation of the diagnostic value of immunoregulation-related piRNAs

With the purpose of understanding the ability of immunoregulation-related piRNAs to distinguish between RA and HC, ROC curves were generated to identify the diagnostic value of those two piRNAs (Figures 2H–J). The AUC values of piR-hsa-27620 and piR-hsa-27124 were respectively 0.7851 [95% CI (0.6982, 0.8721)] and 0.7390 [95% CI (0.6416, 0.8363)], which both had high sensitivity (79.00% and 81.50%, respectively) and specificity (69.00% and 59.50%, respectively) when the largest Youden’s index was the optimal operating point (Table 3). Logistic regression analysis with ROC curves was used to construct the combined indicators. Compared with the discriminative ability of clinical indicators, piRNA alone has insufficient diagnostic ability; however, when combining piRNAs with other indicators, the discriminative ability of the combination index can exceed that of any single clinical indicator.

3.7 PIWIL protein profiles in the peripheral blood of RA patients and HC

The PIWI proteins belong to the Argonaute family, which has crucial functions in piRNA biogenesis and functional complex formation. We first measured the expression of the PIWIL1-4 mRNA. PIWIL2 and PIWIL4 were expressed similarly in RA and HC, while the mRNAs of PIWIL1 and PIWIL3 were not detectable (Figure 3A). In addition, we also examined genes (DDX39B, MYBL1, HENMT1, and PLD6) related to piRNA metabolism and found that HENMT1 was significantly overexpressed in RA patients (Figure 3B).

Figure 3

3.8 Correlation between piRNAs/PIWI system and RA clinical characteristics

Our data have suggested that the expression levels of piR-hsa-27620, piR-hsa-27124, and HENMT1 were elevated in the peripheral blood of RA patients. To further study their association with disease characteristics, we performed correlation analysis and presented the groups with significant correlations in the form of heatmap and scatter plots (Figure 3C, Figure S3). The results showed that piR-hsa-27620 was positively correlated with IgA (r_s = 0.4616, P = 0.0176), while negatively correlated with IgM (r_s = −0.5133, P = 0.0087). Furthermore, piR-hsa-27124 was positively correlated with neutrophil count (r_s = 0.3069, P = 0.0481), and CRP levels (r_s = 0.3190, P = 0.0395). Furthermore, despite the lack of direct correlation between disease duration, DAS28-CRP, and piRNAs expression, we discovered through stratified analysis that piR-hsa-27124 dramatically decreased in patients with high disease activity and disease duration of more than 1 year (Figures 3D, E).

HENMT1, a methyltransferase that adds a 2’-O-methyl group at the 3’-end of piRNAs, increased gradually over the course of the disease (r_s = 0.6119, P = 0.0200) and was negatively correlated with CRP, C3, and anti-CCP (r_s = −0.5308, P = 0.0440; r_s = −0.7720, P = 0.0420; r_s = −0.5542, P = 0.0397, respectively) (Figures 3F–I).

4 Discussion

In this study, 1,565 piRNAs were identified in the peripheral leukocytes of RA patients from southwest China by small RNA sequencing, and 24 piRNAs were differentially expressed significantly. Then the GO term and KEGG pathway enrichment analysis were conducted, and seven immunoregulation-related piRNAs were selected; among them, piR-hsa-27620 and piR-hsa-27124 were not only significantly upregulated in patients with rheumatoid arthritis but also may be involved in immune regulation, which has potential research value as biological markers of RA.

RA is a chronic inflammatory rheumatoid disease mediated by immunity, which will eventually cause disabling joint destruction and is closely related to genetic background. The piRNA/PIWI system has long been reported to be involved in epigenetic regulation (5, 6, 33), which is also associated with RA (34–36). Transcriptome sequencing data show that piRNA clusters exist in a variety of immune cells related to inflammatory reactions, such as neutrophils, macrophages, dendritic cells, B cells, CD4+T cells, and NK cells (37–39). Therefore, we analyzed and validated that the piRNA/PIWI system is dysregulated in the peripheral leukocytes of RA.

Studies have found that piRNAs may be the largest class of non-coding RNAs and have more regulation modes compared with other ncRNAs (6, 40), which means that the biological processes affected by dysregulated piRNAs will be more extensive in the pathogenesis of diseases. Our sequencing results also showed that the number of fragments in the piRNA region at 26–28 nt was higher than in other regions. Bioinformatics analysis showed that DE piRNAs might participate in the transcriptional regulation of various RA regulatory molecules, such as STAT3, FCGR2A, CXCL12, etc. GO analysis results showed that piRNA might intervene in the disease by participating in transcriptional regulation, intracellular signal transduction, RNA polymerase II promoter transcriptional regulation, cell differentiation, and other processes, indicating that the whole functional system of piRNAs exists in RA and may change in response to disease status. KEGG pathway analysis results showed that piRNA may be involved in many RA-related signaling pathways, including the Rap1 signaling pathway, the PI3K–Akt signaling pathway, the MAPK signaling pathway, Th17 cell differentiation, apoptosis mediated by natural killer cells, and osteoclast differentiation. These are all available in Table S4. Therefore, we can speculate that piRNA may be involved in the pathogenesis of RA through transcription and post-transcriptional regulation, affecting immune cell differentiation and the secretion of inflammatory factors. Additionally, according to the GEO database, we analyzed the expression of genes targeted by immunoregulation-related piRNAs in the peripheral blood and synovial tissues of RA patients. The results showed that many gene targets are dysregulated in RA patients, which suggests indirectly the potential function of dysregulated piRNAs. Furthermore, it has been reported that the piRNA/PIWI system can also stabilize the expression of target genes and induce translational activation according to the different binding sites, in addition to negative regulation (6). That is to say, the expression of gene targets regulated by piRNAs does not always exhibit an inverse trend to that of the piRNAs themselves. Similarly, we also observed that upregulated piRNAs may target upregulated genes, for instance, piR-hsa-27620, STAT1, piR-hsa-27400, and TNFSF13. Of course, it is still imperative to conduct clinical studies on the same patient population and molecular mechanisms to obtain direct evidence, which will help illustrate the position of immunoregulation-related piRNAs in the inflammatory pathway.

Since discovering that the expression of piRNA is not limited to germ cell lines and could be used for auxiliary diagnosis and monitoring of diseases (23, 24), more and more studies have demonstrated this point of view. Mai et al. demonstrated that piR-hsa-54265 is an oncogenic RNA in the development of colorectal cancer (CRC) and a valuable biomarker better than other serum tumor markers routinely used in clinic, which would decrease dramatically after surgical resection of CRC but increase again when tumor relapses (41). Lipps et al. identified the piRNA DQ593039 as a promising biomarker reflecting the severity of pulmonary hypertension (42). The results of transcriptome-wide piRNA profiling revealed that up to 103 piRNAs nominally correlated with genome-wide significant risk SNPs for Alzheimer’s disease (AD) (43). With these studies, our research demonstrated that piRNA expression is indeed not limited to germlines and might have a potential correlation with rheumatoid arthritis. The ROC curve built with 42 RA patients and 81 controls yielded that piR-hsa-27620 had the best diagnostic accuracy [AUC (95% CI): 0.7851 (0.6982, 0.8721), SEN: 79.00%, SPE: 69.00%], while piR-hsa-27124 had the highest sensitivity [AUC (95% CI): 0.7390 (0.6416, 0.8363), SEN: 81.50%, SPE: 59.50%], which have good abilities to distinguish patients with rheumatoid arthritis from healthy individuals.

To better demonstrate the landscape of the piRNA/PIWI system in RA peripheral blood, we assessed the expression of PIWI protein and piRNA biogenesis-related enzymes subsequently. Similar to the synovial tissue of RA (25), only PIWIL2 and PIWIL4 expressions were detected in the peripheral blood of RA. Due to the large number of enzymes involved in piRNA biogenesis (5), we first observed their expression in peripheral blood based on preliminary sequence data from public databases, and then selected four representative enzymes for further verification based on their functional roles in piRNA processing. Among them, MYBL1 can promote transcription of piRNA clusters, DDX39B is important for enabling nuclear export of piRNA precursor transcripts to sites of piRNA production, HENMT1 is required to modify the 2’ hydroxyl at the piRNA 3’ end, and PLD6 participates in establishing the 3’ ends of responder pre-piRNAs. Our findings suggested that except for PLD6, the other three enzymes were all expressed in RA peripheral blood, and HENMT1 was significantly elevated. Interestingly, studies have already found that DDX39B is associated with RA susceptibility independently of HLA-DRB1 (44), further strengthening the evidence of the relationship between piRNAs and RA. Moreover, though PLD6 was not detectable in RA, it is reported that carrying HLA-B*27 was associated with robust hypomethylation of PLD6 in ankylosing spondylitis patients (45). Afterward, correlation analysis was conducted to better align these results with clinical information. According to our results, neutrophil, an essential immune cell in modulating autoimmunity (46), may have more piR-hsa-27124 and piR-hsa-25672. Positively correlated with IgA level, piR-hsa-27620 might be involved in RA mucosal immunity. Besides, piR-hsa-27124 was associated with inflammation, and patients with moderate symptoms in the early stages of RA exhibited higher expression levels than those with more severe symptoms and longer disease courses. Meanwhile, HENMT1, necessary for piRNA stability, presented an increasing trend with the duration of RA, implying that the piRNA/PIWI system might initiate the stabilization mechanism in response in the later stages of the disease.

The piRNA/PIWI system expression profile in the peripheral leukocytes of RA patients was discussed and verified for the first time in our project, and it was found that the piRNA/PIWI system was dysregulated and that two immunoregulation-related piRNAs, i.e., piR-hsa-27620 and piR-hsa-27124, had value as potential biomarkers to help improve the diagnostic ability of current indicators. This project also had some limitations. First, we only explored the landscape of the piRNA/PIWI system in new-onset untreated RA patients; therefore, further exploration is required to determine the impact of clinical treatment on this system. Second, due to the limitations of current studies on the function of piRNA, this study only made superficial preliminary assumptions through bioinformatics analysis and comparison with clinical indicators and could not clarify the role of piRNAs in the pathogenesis of RA in detail. The deficiencies of the current study point to the direction of our next experiment. It is believed that with the expansion of the breadth and depth of the study, the function of immunoregulation-related piRNAs will be revealed gradually.

References

• 1

  SmolenJSAletahaDBartonABurmesterGREmeryPFiresteinGSet al. Rheumatoid arthritis. Nat Rev Dis Primers (2018) 4:18001. doi: 10.1038/nrdp.2018.1

  • CrossRef
  • Google Scholar
  • View reference in article

• 2

  AletahaDSmolenJS. Diagnosis and management of rheumatoid arthritis: a review. JAMA (2018) 320(13):1360–72. doi: 10.1001/jama.2018.13103

  • CrossRef
  • Google Scholar
  • View reference in article

• 3

  OkadaYEyreSSuzukiAKochiYYamamotoK. Genetics of rheumatoid arthritis: 2018 status. Ann Rheum Dis (2019) 78(4):446–53. doi: 10.1136/annrheumdis-2018-213678

  • CrossRef
  • Google Scholar

• 4

  WangJYanSYangJLuHXuDWangZ. Non-coding RNAs in rheumatoid arthritis: from bench to bedside. Front Immunol (2019) 10:3129. doi: 10.3389/fimmu.2019.03129

  • CrossRef
  • Google Scholar
  • View reference in article

• 5

  OzataDMGainetdinovIZochAO'CarrollDZamorePD. PIWI-interacting RNAs: small RNAs with big functions. Nat Rev Genet (2019) 20(2):89–108. doi: 10.1038/s41576-018-0073-3

  • CrossRef
  • Google Scholar
  • View reference in article

• 6

  RamatASimoneligM. Functions of PIWI proteins in gene regulation: new arrows added to the piRNA quiver. Trends Genet (2021) 37(2):188–200. doi: 10.1016/j.tig.2020.08.011

  • CrossRef
  • Google Scholar
  • View reference in article

• 7

  AravinAGaidatzisDPfefferSLagos-QuintanaMLandgrafPIovinoNet al. A novel class of small RNAs bind to MILI protein in mouse testes. Nature (2006) 442(7099):203–7. doi: 10.1038/nature04916

  • CrossRef
  • Google Scholar
  • View reference in article

• 8

  WatanabeTTakedaATsukiyamaTMiseKOkunoTSasakiHet al. Identification and characterization of two novel classes of small RNAs in the mouse germline: retrotransposon-derived siRNAs in oocytes and germline small RNAs in testes. Genes Dev (2006) 20(13):1732–43. doi: 10.1101/gad.1425706

  • CrossRef
  • Google Scholar
  • View reference in article

• 9

  RooversEFRosenkranzDMahdipourMHanCTHeNChuva de Sousa LopesSMet al. Piwi proteins and piRNAs in mammalian oocytes and early embryos. Cell Rep (2015) 10(12):2069–82. doi: 10.1016/j.celrep.2015.02.062

  • CrossRef
  • Google Scholar
  • View reference in article

• 10

  HuangXAYinHSweeneySRahaDSnyderMLinH. A major epigenetic programming mechanism guided by piRNAs. Dev Cell (2013) 24(5):502–16. doi: 10.1016/j.devcel.2013.01.023

  • CrossRef
  • Google Scholar
  • View reference in article

• 11

  La GrecaAScarafiaMAHernandez CanasMCPerezNCastanedaSColliCet al. PIWI-interacting RNAs are differentially expressed during cardiac differentiation of human pluripotent stem cells. PloS One (2020) 15(5):e0232715. doi: 10.1371/journal.pone.0232715

  • CrossRef
  • Google Scholar
  • View reference in article

• 12

  HenaouiISJacovettiCGuerra MolletIGuayCSobelJEliassonLet al. PIWI-interacting RNAs as novel regulators of pancreatic beta cell function. Diabetologia (2017) 60(10):1977–86. doi: 10.1007/s00125-017-4368-2

  • CrossRef
  • Google Scholar
  • View reference in article

• 13

  RajasethupathyPAntonovISheridanRFreySSanderCTuschlTet al. A role for neuronal piRNAs in the epigenetic control of memory-related synaptic plasticity. Cell (2012) 149(3):693–707. doi: 10.1016/j.cell.2012.02.057

  • CrossRef
  • Google Scholar
  • View reference in article

• 14

  HashimARizzoFMarcheseGRavoMTaralloRNassaGet al. RNA Sequencing identifies specific PIWI-interacting small non-coding RNA expression patterns in breast cancer. Oncotarget (2014) 5(20):9901–10. doi: 10.18632/oncotarget.2476

  • CrossRef
  • Google Scholar

• 15

  LawPTQinHChingAKLaiKPCoNNHeMet al. Deep sequencing of small RNA transcriptome reveals novel non-coding RNAs in hepatocellular carcinoma. J Hepatol (2013) 58(6):1165–73. doi: 10.1016/j.jhep.2013.01.032

  • CrossRef
  • Google Scholar
  • View reference in article

• 16

  EcksteinMJungRWeigeltKSikicDStohrRGeppertCet al. Piwi-like 1 and -2 protein expression levels are prognostic factors for muscle invasive urothelial bladder cancer patients. Sci Rep (2018) 8(1):17693. doi: 10.1038/s41598-018-35637-4

  • CrossRef
  • Google Scholar
  • View reference in article

• 17

  AiLMuSSunCFanFYanHQinYet al. Myeloid-derived suppressor cells endow stem-like qualities to multiple myeloma cells by inducing piRNA-823 expression and DNMT3B activation. Mol Canc (2019) 18(1):88. doi: 10.1186/s12943-019-1011-5

  • CrossRef
  • Google Scholar
  • View reference in article

• 18

  RajanKSVelmuruganGGopalPRamprasathTBabuDDKrithikaSet al. Abundant and altered expression of PIWI-interacting RNAs during cardiac hypertrophy. Heart Lung Circ (2016) 25(10):1013–20. doi: 10.1016/j.hlc.2016.02.015

  • CrossRef
  • Google Scholar
  • View reference in article

• 19

  ZhangXLiangZZhangYZhuMZhuYLiSet al. Specific PIWI-interacting small noncoding RNA expression patterns in pulmonary tuberculosis patients. Epigenomics (2019) 11(16):1779–94. doi: 10.2217/epi-2018-0142

  • CrossRef
  • Google Scholar
  • View reference in article

• 20

  KimKW. PIWI proteins and piRNAs in the nervous system. Mol Cells (2019) 42(12):828–35. doi: 10.14348/molcells.2019.0241

  • CrossRef
  • Google Scholar
  • View reference in article

• 21

  Della BellaEMenzelUBasoliVTourbierCAliniMStoddartMJ. Differential regulation of circRNA, miRNA, and piRNA during early osteogenic and chondrogenic differentiation of human mesenchymal stromal cells. Cells (2020) 9(2):398. doi: 10.3390/cells9020398

  • CrossRef
  • Google Scholar
  • View reference in article

• 22

  ZhongFZhouNWuKGuoYTanWZhangHet al. A SnoRNA-derived piRNA interacts with human interleukin-4 pre-mRNA and induces its decay in nuclear exosomes. Nucleic Acids Res (2015) 43(21):10474–91. doi: 10.1093/nar/gkv954

  • CrossRef
  • Google Scholar
  • View reference in article

• 23

  MeiYClarkDMaoL. Novel dimensions of piRNAs in cancer. Cancer Lett (2013) 336(1):46–52. doi: 10.1016/j.canlet.2013.04.008

  • CrossRef
  • Google Scholar
  • View reference in article

• 24

  YangXChengYLuQWeiJYangHGuM. Detection of stably expressed piRNAs in human blood. Int J Clin Exp Med (2015) 8(8):13353–8.

  • Google Scholar
  • View reference in article

• 25

  PlestilovaLNeidhartMRussoGFrank-BertonceljMOspeltCCiureaAet al. Expression and regulation of PIWIL-proteins and PIWI-interacting RNAs in rheumatoid arthritis. PloS One (2016) 11(11):e0166920. doi: 10.1371/journal.pone.0166920

  • CrossRef
  • Google Scholar
  • View reference in article

• 26

  AletahaDNeogiTSilmanAJFunovitsJFelsonDTBinghamCO3rdet al. 2010 Rheumatoid arthritis classification criteria: an American college of Rheumatology/European league against rheumatism collaborative initiative. Arthritis Rheum (2010) 62(9):2569–81. doi: 10.1002/art.27584

  • CrossRef
  • Google Scholar
  • View reference in article

• 27

  Griffiths-JonesSBatemanAMarshallMKhannaAEddySR. Rfam: an RNA family database. Nucleic Acids Res (2003) 31(1):439–41. doi: 10.1093/nar/gkg006

  • CrossRef
  • Google Scholar
  • View reference in article

• 28

  Griffiths-JonesSSainiHKvan DongenSEnrightAJ. miRBase: tools for microRNA genomics. Nucleic Acids Res (2008) 36(Database issue):D154–8. doi: 10.1093/nar/gkm952

  • CrossRef
  • Google Scholar
  • View reference in article

• 29

  WangJZhangPLuYLiYZhengYKanYet al. piRBase: a comprehensive database of piRNA sequences. Nucleic Acids Res (2019) 47(D1):D175–D80. doi: 10.1093/nar/gky1043

  • CrossRef
  • Google Scholar
  • View reference in article

• 30

  JinYHammellM. Analysis of RNA-seq data using TEtranscripts. Methods Mol Biol (2018) 1751:153–67. doi: 10.1007/978-1-4939-7710-9_11

  • CrossRef
  • Google Scholar
  • View reference in article

• 31

  GhoshBSarkarAMondalSBhattacharyaNKhatuaSGhoshZ. piRNAQuest V.2: an updated resource for searching through the piRNAome of multiple species. RNA Biol (2022) 19(1):12–25. doi: 10.1080/15476286.2021.2010960

  • CrossRef
  • Google Scholar
  • View reference in article

• 32

  FaulFErdfelderEBuchnerALangAG. Statistical power analyses using G*Power 3.1: tests for correlation and regression analyses. Behav Res Methods (2009) 41(4):1149–60. doi: 10.3758/BRM.41.4.1149

  • CrossRef
  • Google Scholar
  • View reference in article

• 33

  RayfordKJCooleyARumphJTArunARachakondaGVillaltaFet al. piRNAs as modulators of disease pathogenesis. Int J Mol Sci (2021) 22(5):2373. doi: 10.3390/ijms22052373

  • CrossRef
  • Google Scholar
  • View reference in article

• 34

  HammakerDFiresteinGS. Epigenetics of inflammatory arthritis. Curr Opin Rheumatol (2018) 30(2):188–96. doi: 10.1097/BOR.0000000000000471

  • CrossRef
  • Google Scholar
  • View reference in article

• 35

  AiRLaragioneTHammakerDBoyleDLWildbergAMaeshimaKet al. Comprehensive epigenetic landscape of rheumatoid arthritis fibroblast-like synoviocytes. Nat Commun (2018) 9(1):1921. doi: 10.1038/s41467-018-04310-9

  • CrossRef
  • Google Scholar
  • View reference in article

• 36

  NairNWilsonAGBartonA. DNA Methylation as a marker of response in rheumatoid arthritis. Pharmacogenomics (2017) 18(14):1323–32. doi: 10.2217/pgs-2016-0195

  • CrossRef
  • Google Scholar
  • View reference in article

• 37

  ZhangXHeXLiuCLiuJHuQPanTet al. IL-4 inhibits the biogenesis of an epigenetically suppressive PIWI-interacting RNA to upregulate CD1a molecules on Monocytes/Dendritic cells. J Immunol (2016) 196(4):1591–603. doi: 10.4049/jimmunol.1500805

  • CrossRef
  • Google Scholar
  • View reference in article

• 38

  BahnJHZhangQLiFChanTMLinXKimYet al. The landscape of microRNA, piwi-interacting RNA, and circular RNA in human saliva. Clin Chem (2015) 61(1):221–30. doi: 10.1373/clinchem.2014.230433

  • CrossRef
  • Google Scholar

• 39

  YanZHuHYJiangXMaierhoferVNebEHeLet al. Widespread expression of piRNA-like molecules in somatic tissues. Nucleic Acids Res (2011) 39(15):6596–607. doi: 10.1093/nar/gkr298

  • CrossRef
  • Google Scholar

• 40

  LinXXiaYHuDMaoQYuZZhangHet al. Transcriptome−wide piRNA profiling in human gastric cancer. Oncol Rep (2019) 41(5):3089–99. doi: 10.3892/or.2019.7073

  • CrossRef
  • Google Scholar

• 41

  MaiDZhengYGuoHDingPBaiRLiMet al. Serum piRNA-54265 is a new biomarker for early detection and clinical surveillance of human colorectal cancer. Theranostics (2020) 10(19):8468–78. doi: 10.7150/thno.46241

  • CrossRef
  • Google Scholar

• 42

  LippsCNorthePFigueiredoRRohdeMBrahmerAKramer-AlbersEMet al. Non-invasive approach for evaluation of pulmonary hypertension using extracellular vesicle-associated small non-coding RNA. Biomolecules (2019) 9(11):666. doi: 10.3390/biom9110666

  • CrossRef
  • Google Scholar

• 43

  QiuWGuoXLinXYangQZhangWZhangYet al. Transcriptome-wide piRNA profiling in human brains of alzheimer's disease. Neurobiol Aging (2017) 57:170–7. doi: 10.1016/j.neurobiolaging.2017.05.020

  • CrossRef
  • Google Scholar

• 44

  Quinones-LombranaALopez-SotoABallina-GarciaFJAlperi-LopezMQueiro-SilvaRLopez-VazquezAet al. BAT1 promoter polymorphism is associated with rheumatoid arthritis susceptibility. J Rheumatol (2008) 35(5):741–4.

  • Google Scholar

• 45

  CoitPKaushikPCaplanLKerrGSWalshJADubreuilMet al. Genome-wide DNA methylation analysis in ankylosing spondylitis identifies HLA-B*27 dependent and independent DNA methylation changes in whole blood. J Autoimmun (2019) 102:126–32. doi: 10.1016/j.jaut.2019.04.022

  • CrossRef
  • Google Scholar

• 46

  O'NeilLJKaplanMJ. Neutrophils in rheumatoid arthritis: breaking immune tolerance and fueling disease. Trends Mol Med (2019) 25(3):215–27. doi: 10.1016/j.molmed.2018.12.008

  • CrossRef
  • Google Scholar