(5R)-5-Hydroxytriptolide (LLDT-8) induces substantial epigenetic mediated immune response network changes in fibroblast-like synoviocytes from rheumatoid arthritis patients

Tripterygium is a traditional Chinese medicine that has widely been used in the treatment of rheumatic disease. (5R)-5-hydroxytriptolide (LLDT-8) is an extracted compound from Tripterygium, which has been shown to have lower cytotoxicity and relatively higher immunosuppressive activity when compared to Tripterygium. However, our understanding of LLDT-8-induced epigenomic impact and overall regulatory changes in key cell types remains limited. Doing so will provide critically important mechanistic information about how LLDT-8 wields its immunosuppressive activity. The purpose of this study was to assess the effects of LLDT-8 on transcriptome including mRNAs and long non-coding RNA (lncRNAs) in rheumatoid arthritis (RA) fibroblast-like synoviocytes (FLS) by a custom genome-wide microarray assay. Significant differential expressed genes were validated by QPCR. Our work shows that 394 genes (281 down- and 113 up-regulated) were significantly differentially expressed in FLS responding to the treatment of LLDT-8. KEGG pathway analysis showed 20 pathways were significantly enriched and the most significantly enriched pathways were relevant to Immune reaction, including cytokine-cytokine receptor interaction (P = 4.61 × 10−13), chemokine signaling pathway (P = 1.01 × 10−5) and TNF signaling pathway (P = 2.79 × 10−4). Furthermore, we identified 618 highly negatively correlated lncRNA-mRNA pairs from the selected significantly differential lncRNA and mRNA including 27 cis-regulated and 591 trans-regulated lncRNA-mRNAs modules. KEGG and GO based function analysis to differential lncRNA also shown the enrichment of immune response. Finally, lncRNA-transcription factor (TF) and lncRNA-TF-mRNA co-expression network were constructed with high specific network characteristics, indicating LLDT-8 would influence the expression network within the whole FLS cells. The results indicated that the LLDT-8 would mainly influence the FLS cells systemically and specially in the process of immune related pathways.


Material and Methods
Patients, cell culture and LLDT-8 treatment. FLS cells were derived from the synovial tissues of patients with RA in Guanghua hospital who had undergone total joint replacement surgery during June to August 2015. Written informed consent to collect FLS cells for LLDT-8 studies was obtained from all participating RA patients. This study was approved by the academic advisory board of Guanghua Hospital (No. 2015-SRA-01) and all methods were performed in accordance with the relevant guidelines and regulations. All patients were random enrolled and all of them fulfilled the American College of Rheumatology classification criteria for RA 28 . Clinical data were collected at the time of sample collection. The detail clinical information is shown in Table 1.
An overview of the cellular processing is shown in Fig. 1. Briefly, joint tissues were minced into 1 × 1 × 1 mm pieces and treated for 2 hours (h) with 2-4 mg/ml of collagenase (Serva, GERMANY) in DMEM at 37 °C in 5% CO2. Dissociated cells were then centrifuged at 500 × g, re-suspended in DMEM supplemented with 10% FCS (Gibco, USA), and plated in 75 cm 2 flasks. The cultures were kept at 37 °C in 5% CO2 and the culture medium was replaced every 2-3 days. Mycoplasma assay was performed on synovial fibroblasts to avoid mycoplasma contamination. When cells approached confluence, they were passed after diluting 1:3 with fresh medium until used. The purity of the cells was verified by flow cytometric analysis. The FLS cells of passages 4th were seeded in 6-well plates. As the control group, the FLS cells were cultured in DMEM including 10 ng/ml TNF-α and 10 ng/ml IL-17  Fig. 1) which is consistent with a previous report 11,29 . In the process to validate the microarray data, TNF-α treatment was applied to simulate a pro-inflammatory reaction and then LLDT8 treatment was applied to evaluate the effect of LLDT-8 on the inflammatory reaction as measured by real-time PCR (RT-PCR).
RNA extraction and microarray hybridization. The total RNA was extracted by TRIzol reagent and quantified by the NanoDrop ND-2000 (Thermo Scientific, Waltham, MA, USA). The integrity of RNA was assessed using the Agilent Bioanalyzer 2100 (Agilent Technologies Inc, Santa Clara, CA, USA). Customized Agilent Human lncRNA microarray (4 × 180 K, Design ID: 062918) was used in this study. 30,656 probes located in mRNA transcripts from Entrez Gene and 78,243 probes in lncRNA. the catalog of the lncRNA were collected from the integration of Broad Institute, Human Body Map lncRNA, TUCP catalog, UCSC lncRNA Transcripts, GENCODE 18, NONCODE V4.0, Ensembl, RefSeq, Ultra-conserved region encoding LncRNA (UCR), lncRNAdb and ncRNA database. The total RNA was transcribed to double-stranded cDNA, synthesized into cRNA (Low Input Quick-Amp Labeling Kit, one-color, Agilent), and labeled with Cyanine-3-CTP. The labeled cRNAs were then hybridized onto the microarray. After washing, the arrays were scanned by the Agilent Scanner G2505C (Agilent Technologies Inc., Santa Clara, CA, USA). Other parameters in ncRNA labeling, microarray hybridization, and washing procedure were performed according to the manufacturer's standard protocols. RT-PCR was applied to validate parts of significant differential genes within same samples with > 3 times technical repeats.
Image scanning and data analysis. The array images were analyzed using Feature Extraction software (version 10.7.1.1, Agilent Technologies Inc., Santa Clara, CA, USA), and the raw data were obtained and basically analyzed with Genespring. The raw data were initially normalized with the quantile algorithm. The probes that at least in one out of two conditions had flags in "P" were chosen for further data analysis. Differentially expressed genes or lncRNAs were identified through a combination of fold change and p-values calculated from a t-test (details in the Statistical analysis section). The threshold set for upregulated and downregulated genes was a fold change ≥ 2.0 and a P ≤ 0.05.
Identification of cis-regulated mRNAs of the differential lncRNAs. Based on these 5 biological replications with or without LLDT-8 treatment, differential expressed lncRNA, cis-regulatory genes/mRNAs were identified. The mRNAs were identified as "cis-regulated mRNAs" when (1) the mRNA loci are within 100 k windows up-and down-stream of the given lncRNA, and (2) the Pearson correlation of lncRNA-mRNA expression is significant (P ≤ 0.01). GO analysis and KEGG analysis were applied to determine the roles of these differentially expressed mRNAs. Finally, Hierarchical Clustering was performed to display the distinguishable genes' expression pattern among samples.
Functional prediction of selected differential lncRNAs. The overall gene function distribution (co-expression, ontology and pathway) of the differential lncRNAs obtained in the experiment was identified as follows. For each differential lncRNA, the Pearson correlation of its expression value with the expression value of each mRNA was calculated, and a P < 0.05 was selected. The enrichment of functional terms in annotation of Figure 1. Flowchart of the study. In the study, we collected 5 RA joint tissues from joint replacement surgery. FLS cells were purified with cell culture and validated by flow cytometry assay. When cells approached confluence, they were passed after diluting 1:3 with fresh medium until used. The purity of the cells was verified by flow cytometric analysis. The FLS, from passages 4, were seeded in 6-well plates. As the control group, the FLS cells were cultured in DMEM including 10 ng/ml TNF-α and 10 ng/ml IL-17 (PeproTech, USA) for 12 h. In the treatment group, LLDT-8 (Shanghai Pharma, Shanghai) was dissolved with 2% DMSO and diluted with DMEM (with TNF-α and IL-17) to 100 nM for 24 hours. RT-PCR, Western blot and microarray were then conducted to different groups. (2019) 9:11155 | https://doi.org/10.1038/s41598-019-47411-1 www.nature.com/scientificreports www.nature.com/scientificreports/ differential expressing gene or co-expressed mRNAs was statistically evaluated using the hypergeometric cumulative distribution function 30,31 .

Identification of transcription factors (TF) associated to differential lncRNAs. The transcription
factor/chromatin regulation complexes that may possibly play a co-regulatory role with lncRNAs were identified 32,33 . In brief, the set between the lncRNA co-expression coding genes and the target genes of transcription factors/chromatin regulation complex was collected respectively. The enrichment level of the set was determined using hypergeometric distribution, thus the transcription factors significantly associated with differential lncR-NAs were finally screened. Finally, the co-expression networks among lncRNA, TF, and target genes were built with Cytoscape 3.4.0 34 .
Western Blot and Immunofluorescence. Western Blot and Immunofluorescence were applied to detect the protein expression change and nuclear localization to show the effect of LLDT-8 treatment to RF cells in our validation stage and the procedures are same as the previous study 2 . Briefly, FLS cells were seeded at a density of 1 × 10 7 cells per well (6-well plates) in DMEM supplemented with 10% FBS and 1% penicillin. After cell attachment, the culture medium was replaced by DMEM supplemented with A panel: DMSO (2%), B panel: TNF-α (10 ng/ml) and IL-17 (10 ng/ml) or C panel: TNF-α (10 ng/ml), IL-17 (10 ng/ml) and LLDT-8 (0, 50 nM, 100 nm, respectively) for 12 hours. The cells were collected and lysed by 1 nM PMSF SDS on ice for 30 min. The cellular lysates were loaded, and proteins were separated on SDS-PAGE and transferred to nitrocellulose filter. The blots were blocked with 5% BSA TBST for 1 h at room temperature, then probed with rabbit anti-mouse antibodies against p-P65(1:1000) and GAPDH (1:6000) overnight at 4 °C. After five washes, the blots were subsequently incubated with a HRP-linked secondary antibody (1:2500) for 90 min at room temperature. The blots were visualized by ECLTM Prime Western Blotting Detection Reagent. FLS with different treatments were fixed using 4% (v/v) formaldehyde in PBS for 15 min and mounted onto slides using VECTASHIELD Antifade Mounting Medium containing DAPI (Vector Laboratories, H-1200) and visualized using an inverted SP5 confocal microscope (Leica).
Statistical analysis. One sample t-test was applied to compare the mean difference of the expressed level before and after the treatment of LLDT-8. All the statistical analysis was made using R version 3.2.2. We selected the top signals using the threshold on both the fold change > 2 and the P-value < 0.05 from mRNA and lncRNA. In addition, we calculated the Pearson correlation for the expression values on every different combination of lncRNA and mRNA. The pair of IncRNA-mRNA will be selected if their correlation test reaches the threshold (P-value < 0.05). Also, in order to ensure the selected pairs are meaningful, we only selected the correlated pairs in which lncRNA plays as the inhibition of mRNA expression (correlation coefficient, r < −0.9). False discover rate (FDR) correction was performed by Benjamini-Hochberg procedure with default R function.
Genome-wide lncRNA profile of RA FLS before and after LLDT-8 treatment. The lncRNAs with both fold change ≥ 2.0 and P-value ≤ 0.05 from the t-test were identified as differentially expressed lncRNAs. Our results showed that 360 lncRNAs of RA FLS were significantly changed with LLDT-8 treatment. Those P-values ranged from 4.88 × 10 −2 to 2.92 × 10 −7 , and the fold change spanned from 2.00 to 18.44. Among these lncRNAs, 56% (203) were downregulated and 44% (157) were upregulated (Supp. Table 2). Our result indicates not only mRNAs but also lncRNA transcripts would be widely changed during the treatment of LLDT-8. In addition, we found the target genes of the differential lncRNAs that could also be a strong indication of the change caused by LLDT-8 treatment. In addition, as Fig. 2E shows, cluster analysis based on the 30 co-expression genes of differential lncRNA (ENST00000584923) could separate the samples quite accurately with or without LLDT-8 treatment, indicating that the lncRNA network and the mRNA network in the cells were highly interacted and coordinately changed by LLDT-8 treatment.  Volcano plot for the differential express genes. B. Heatmap plot for differential expressed genes to show the effect of LLDT-8 to FLS cells. C. qPCR validation to significant differential expressed genes from microarray. D. Gene ontology analysis to differential expressed genes (N = 3 technical repeats in control, TNF-α and LLDT8 treatment group). TNF-α treatment is applied to simulate inflammatory reaction and then LLDT8 treatment is applied to check the effect of LLDT8 to inflammatory reaction. E. Heatmap plot to the co-expressed mRNA of differential lncRNA have the power to separate samples with or without LLDT-8 treatment.
mRNA was selected if their correlation test reached the threshold (P-value < 0.05). Also, in order to ensure the selected pairs is meaningful, we only kept those correlated pairs which show strong negative correlation (correlation coefficient < −0.9). Our result indicated that there were 9,666 pairs of lncRNA-mRNA with strong negative correlation. Moreover, 618 of them were also identified as the significantly differential mRNA and lncRNA from the t-test in which 70 pairs of them were located in the same chromosome. Furthermore, 13 out of those 70 pairs were located in chromosome-6, and 10 out of those them were located in chromosome-17. Such kind of chromosome-distribution disequilibrium suggests the influence of the LLDT-8 is not a random effect but with specific target and influence to FLS. 27 cis-regulated (Table 2 and Supp. Table 3) and 591 trans-regulated lncRNA-mRNAs modules were identified (Table 3 and Supp. Table 4).
Functional prediction of selected differential lncRNAs analysis were conducted both with GO and KEGG strategies. Gene Ontology analysis showed that the term of innate immune responses was significantly enriched and 17 differential lncRNAs were identified to be involved (Supp. Table 5). KEGG analysis identified term of rheumatoid arthritis was significant enriched and 4 differential lncRNAs (NONHSAG028996, NONHSAT142637, ENST00000584934 and NR_028330.1) were identified to be involved (Supp. Table 6).
We also identified large number of co-expressed lncRNA and TFs. The co-expression network between TF and lncRNA was shown in Fig. 4A. We discovered one large sub-network and large number isolated interaction between TFs and lncRNA. This network indicated lncRNA and TF would regulate the gene expression together in FLS cells. When we build lncRNA, TF, target gene network together (Fig. 4B), the result is quite similar as our expectation, and large numbers of genes were identified regulated by lncRNA and TF simultaneously.

Cellular regulation effect of LLDT-8 treatment to FLS involved in NF-κB pathway.
In order to validate our discoveries that LLDT-8 could provide pharmacological effect to RA therapy via the interaction with rheumatology related pathway (Fig. 2D), we validate the effect of LLDT-8 to one of most important www.nature.com/scientificreports www.nature.com/scientificreports/ rheumatology related pathway, NF-kappa B signaling pathway. Western-blot shown that the protein level of p-p65 are strongly inhibited with increasing concentration of LLDT-8 treatment to FLS cells stimulated by TNF-α and IL-17 (Fig. 5A). In addition, we also demonstrate that LLDT-8 would inhibit the nuclear translocation of the p65 (Fig. 5B) which is the significant effect of TNF and IL-17, indicating LLDT-8 has significant inhibition to TNF and IL-17 effects and therefore have potential pharmacological effect to RA therapy.

Discussion
Tripterygium is a traditional Chinese medicine, which has widely been used in the treatment of rheumatic disease. (5R)-5-hydroxytriptolide (LLDT-8) is an extracted compound from Tripterygium and has been showed lower cytotoxicity and relatively higher immunosuppressive activity. However, how LLDT-8 influences RA FLS cells is still unknown, especially in the level of lncRNAs. In this study of genome-wide microarray assay, we identified large number of lncRNA and mRNAs responding to the LLDT-8 treatment. In addition, the first LLDT-8  Table 3. Most significant cis-lncRNA-mRNA regulation pairs for the significant differential lncRNA Footnote: Coordination is based on GRCh37. P-value A indicate the significance of the differential expression while P-value B indicate the significance of the correlation between lncRNA and gene expression. Abbreviation: FC, Fold change; R, Correlation coefficient. www.nature.com/scientificreports www.nature.com/scientificreports/ regions 23,35,36 . Considering that our result showed the significant change of the lncRNA during LLDT-8 treatment, lncRNA would be play important role in the pathogenesis and therapy of RA.
Our microarray data are consistent with protein results in our previous study in which we found LLDT-8 inhibited IL-1β, IL-6, IL-21 and IL-23 while LLDT-8 increased the IL-10 protein level. Meanwhile, our research were supported by numerous of previous research, such as Dr. Karouzakis and colleagues found synovial fibroblasts from early rheumatoid arthritis shown DNA methylation changes 37 . Our research demonstrate the immunosuppressive function of LLDT-8 was related to the interaction with ncRNA network. Our results were also consistent with n'sh (PMC2833606) in which they found CCL8 are abundantly-expressed in RA synovial tissue while we found LLDT-8 could significant decrease CCL8 (Table 3). In another study 38 , p38 MAPK inhibition treatment showed increased TNFRSF18 which is similar with LLDT-8 treatment in our study which also highlights the anti-inflammatory function of LLDT-8 38 . What's more, FAM129A was demonstrated low-expression in RSK2 deficiency mice which have earlier and exacerbated inflammation of bone destruction and in our study, we found FAM129A was significantly up-regulated after LLDT-8 treatment 39 . All these evidence support our conclusion that LLDT-8 own immune-suppressive function and epigenetic regulation mediated by lncRNA might play important roles.
We conducted the power analysis to estimate the minimum samples size to identify significant (q-value < 0.5) mRNA changes based on our dataset GSE84074 with bootstrap resampling. In the bootstrap resampling, we selected several important rheumatoid arthritis genes such as MMP genes (MMP1, MMP3), CCL genes (CCL3, CCL4, and CCL26), and IL genes (IL6, IL8, IL19, and IL24). We found when the sample size is larger than 4, majority of these rheumatoid arthritis genes can be detected with power > 0.8 ( Supplementary Fig. 2), therefore, we believe our sample size is sufficient to identify the most interesting inflammation and autoimmune associated signals.
Our study also provided a probability to compare the lncRNA profiles with RA, primary Sjögren's syndrome 40 , gastric carcinogenesis 41 , endometrial carcinoma 42 , bladder cancer 43 , esophageal squamous cell carcinoma 44 , non-small-cell lung cancer 45 and cervical cancer 46 , since our data has been deposited to GEO database. These diseases might share some common aberrant lncRNAs relevant with inflammation and immune response and therefore provide some special insight to the therapy of RA.
In summary, our results provided a good prospect in further clinical tests of LLDT-8 for its therapeutic potential in the treatment of RA. We also suggest that more studies need to be done in order to further investigate the biological functions of IncRNAs in more patient samples. However, those initial results suggest that lncRNA would be a good biomarker and drug target in the future drug discovery, especially for LLDT-8. Our study also demonstrate epigenetics including DNA methylation [47][48][49] and lncRNA 50,51 are playing important roles in etiology and therapy of RA.

Conclusion
Genome-wide lncRNA and transcriptome analysis of response to treatment with (5R)-5-Hydroxytriptolide (LLDT-8) in cells shown LLDT-8 would mainly influence the RA cells systemically and especially in the process of immune network regulation.

Data Availability
The microarray data was deposited in the Gene Expression Omnibus (GEO accession: GSE84074).