Insights into the mechanisms of triptolide nephrotoxicity through network pharmacology-based analysis and RNA-seq

Introduction Triptolide (TPL) is a promising plant-derived compound for clinical therapy of multiple human diseases; however, its application was limited considering its toxicity. Methods To explore the underlying molecular mechanism of TPL nephrotoxicity, a network pharmacology based approach was utilized to predict candidate targets related with TPL toxicity, followed by deep RNA-seq analysis to characterize the features of three transcriptional elements include protein coding genes (PCGs), long noncoding RNAs (lncRNAs) and circular RNAs (circRNAs) as well as their associations with nephrotoxicity in rats with TPL treatment. Results & Discussion Although the deeper mechanisms of TPL nephrotoxcity remain further exploration, our results suggested that c-Jun is a potential target of TPL and Per1 related circadian rhythm signaling is involved in TPL induced renal toxicity.


Introduction
Over 3000 years of constant practice and optimization for the system of Traditional Chinese Medicine (TCM) have endowed its specific tradition that treasures in both scientific and medical fields (Li et al., 2007). A better understanding of therapeutic mechanisms of herb and herbal formulas from TCMs is of great significance for pharmacological study as they have played vital roles in clinical practice (Zuo et al., 2018). As one of the most renowned traditional Chinese medical herbs, Tripterygium wilfordii Hook f. (TWHF) has been applied in the treatment of multiple renal diseases such as membranous nephropathy (MN), nephrotic syndrome (NS) and refractory proteinuria since ancient China. Triptolide (TPL) is a major active component of TWHF as well as a promising compound for cancer therapy (Noel et al., 2019). Increasing evidence suggests that TPL can attenuate the progression of several types of tumor via varieties of approaches including target epigenetic networks (Noel et al., 2020), induce cancer cell apoptosis, enhance the effect of radiotherapy, inhibit metastasis and etc (Meng et al., 2014). TPL also shows potential immunosuppressive effect in autoimmune diseases treatment such as rheumatoid arthritis (Fan et al., 2018). However, the clinical application of TPL is restricted due to its hepatic, nephric, heart and gastrointestinal toxicity . The cytotoxic activities of TPL include introducing DNA damage and apoptosis, arresting cell cycle (Park and Kim, 2013), autophagy (You et al., 2018), and it involves in the production of reactive oxygen species (ROS), generation and depolarization of mitochondrial membrane potential (MMP) in different cell lines .
The advancement of bioinformatics as well as the booming development of compound/drug/diseases databases such as TCMSP (Ru et al., 2014), NIMS (Li et al., 2011) and comCIPHER (Zhao and Li, 2012) have facilitated network pharmacology as a feasible approach to explicate the material composition and molecular mechanism of drugs effectively since it seeks targets by constructing distinct networks and evaluating the molecular connections involved in the process of drug treatment . Network pharmacology has greatly enhanced the investigation of the molecular basis of herbal formula in the past decade (Li and Zhang, 2013). Through network pharmacology, Li et al. revealed the targets and pathways of niacin in the treatment of COVID-19 and colorectal cancer . Niu et al. found that IL6 is potentially regulated by phytochemicals in traditional Chinese medicine for COVID-19 treatment (Niu et al., 2021). On the other hand, RNA-seq has been widely used to affiliate the expression patterns of protein coding gene (PCG), long noncoding RNA (lncRNA) and circular RNA (circRNA). Increasing evidence has shown that lncRNAs and circRNAs are closely related with degenerative diseases (Bhatti et al., 2021), cancers (Anastasiadou et al., 2018) development (Fatica and Bozzoni, 2014;Di Agostino et al., 2020), aging (Jiang et al., 2021;Ge et al., 2022), and they have great potential to be utilized as drug targets in the near future (Matsui and Corey, 2017;He et al., 2021).
Although previous renal metabolic analysis revealed that Toll-like receptor signaling pathway and NF-kB signaling pathway played an important role in TPL-induced nephrotoxicity , the signatures of transcriptional elements are largely unexplored. In this study, we employed deep RNA-seq in female rat kidneys as well as network pharmacology-based analysis, to elucidate the principles of transcriptomic changes (include protein coding genes, lncRNAs and circRNAs) that associated with TPL and identify candidate targets for a better understanding of TPL renal toxicology.

Animals, pathological measurements and ethic statements
Female Sprague-Dawley (SD) rats, weighing 170-190g, were purchased from Guangdong medical laboratory animal center (Guangzhou, China) and housed in the animal facility of our institute under a pathogen-free condition. Rats were fed in an ad arbitrium diet and with free access to water. TPL was purchased from MedChemExpress (New Jersey, USA). Rats were randomly divided into control (Ctrl, n = 3), low dosage of TPL (L-TPL, n = 6) and high dosage (H-TPL, n = 6) groups. The L-TPL and H-TPL rats were separately administrated by oral gavage at a dose of 0.2 mg/kg and 0.4mg/kg for 28 days. Blood samples were collected for testing blood urea nitrogen (BUN) and serum creatinine (Scr) using one-way anova method among three groups. Coronal renal tissue was sectioned for H & E staining following standard protocols. Renal parenchyma was dissected for RNA-seq. The animal protocol of this study was approved by the institutional ethics review board of Shenzhen PKU-HKUST Medical Center (No. 2020252) and the authors declare that all the procedures have carefully followed the animal protocol. This study was in accordance with ARRIVE guidelines (https://arriveguidelines.org).

RNA isolation and sequencing
Three rats per group were randomly selected from Ctrl and H-TPL groups for RNA-seq. Total RNA was extracted from kidney using Trizol reagent (Invitrogen Cat#15596026) following standard protocols and subjected to the preparation of ribosome depletion RNA sequencing library by illumina platform.

Data availability
The annotation files of novel lncRNAs and circRNAs and the raw data were submitted to the Genome Sequence Archive in BIG Data Center, (Beijing Institute of Genomics (BIG), Chinese Academy of Sciences (http://bigd.big.ac.cn/gsa) , under the bioproject PRJCA010363 with accession No. CRA008544.

Target prediction by network pharmacology-based analysis
A step-wise workflow was utilized to predict candidate targets related with TPL nephrotoxicity. Firstly, TPL related targets (TPL-RT) were collected from TCMSP database (http://tcmspw.com/ tcmsp.php) and published studies from PubMed (https:// pubmed.ncbi.nlm.nih.gov/) database that related with TPL. Then, we used "nephrotoxicity" as the keyword to acquire the known nephrotoxicity related targets (nephrotoxicity-RT) from GeneCard, OMIN and DRUGBANK databases, respectively. The overlapped targets (OT) between TPL-RT and nephrotoxicity-RT were retained and subjected to STRING database to construct their protein-protein interaction (PPI) networks. Protein pairs with correlation r-value > 0.9 were regarded as high-quality networks (hq-ntw) and were visualized by cytoscape (Shannon et al., 2003). Gene enrichment analysis was employed to classify proteins within hq-ntw.

Molecular docking analysis
The 2D structure of TPL was downloaded from PubChem database (https://pubchem.ncbi.nlm.nih.gov/), then the structure was subjected to optimization by Chem3D software (https:// library.bath.ac.uk/chemistry-software/chem3d). PyMOL (https:// pymol.org/2/) was utilized to remove the water molecues and small ligands from the protein structures of targets downloaded from PDB database (http://www.rcsb.org/) for subsequent step. The molecular docking was finally performed and visualized using the hydrogen bonded protein structure and optimized TPL structure via AutoDockTools software (https://www.scripps.edu/sanner/software/ adt/Tutorial/index.html).

Western-blot
Western-blot assay was performed as we have previously described (Jiang et al., 2021), blots were cut according to the sizes of target proteins prior to hybridisation with antibodies during blotting and exposed by Bio-rad imaging system. Antibody information see Supplementary Table 1.

Novel circRNA identification
Candidate circRNA were identified using two tools include find_circ (Memczak et al., 2013) and CIRCexplorer  following the standard tutorials with default parameters. Only transcripts that recognized as circRNA by both find_circ and CIRCexplorer were subjected for further analysis. Spliced reads per billion mapping (SRPBM) value for circRNA was calculated through: SRPBM = Circular reads Ã 10 9 Total mapped reads Ã read length

Prediction of interactions between circRNAs and miRNAs
The potential interactions between circRNAs and miRNAs were predicted by miRanda (John et al., 2004) using default parameters.

Identification of differentially expressed PCGs and lncRNAs
RPKM (Reads per Kilobase per Million Reads) was calculated via formula: RPKM = Total exon reads mapped reads(millions) Ã exon length(kb) . By comparing the RPKM values, thresholds of |log(Fold change)| > 1 and p-value< 0.05 were set to define significantly differentially expressed genes. False discovery rate (FDR) was used for adjusting p-value. Unsupervised clustering was employed to uncover unknown relationships among genes and biological samples.
Weighted gene co-expression network analysis analysis WGCNA analysis was performed following its official tutorial . Briefly, the normalized FPKM values of PCGs and lncRNAs were pooled and generated to adjacency matrix and subjected to "dynamicTreeCut" package  to filter out outliner samples. Then we used "pickSoftThreshold" function to calculate soft power values for predicting block-wise modules.

Gene enrichment analysis
Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway mapping are useful strategies for gene functional classification (Kanehisa, 2019;Gene Ontology Consortium, 2021). Genes were classified to Gene Ontology and KEGG terms via online tool DAVID (https://david.ncifcrf.gov/) with default parameters.

Renal pathological changes
The serum creatinine (Scr) and blood urea nitrogen (BUN) were tested and we found that H-TPL rats showed significant higher Scr (p-value = 0.0189) than Ctrl and L-TPL, indicated that treatment of high dosage of TPL induced declined renal function ( Figures 1A, B). H & E staining found that both L-TPL and H-TPL rats exhibited tubular atrophy and renal blood vessel congestion comparing with Ctrl group (Figures 1C, D),. A total of 10 representative views (3-4 views/rat, n=3) were selected for each group for renal tubular injury scores evaluations, non-significant differences was observed between L-TPL and H-TPL( Figure S1, S2). However, wider intercellular gap between renal tubules was observed in H-TPL than L-TPL ( Figures 1C, D).

Potential targets of TPL nephrotoxicity predicted by network pharmacology-based analysis
Network pharmacology-based analysis was performed to predict potential targets related with TPL renal toxicity. A total of 537 targets were predicted by GeneCards, OMIM, and DRUGBANK databases and 31 targets by the TCMSP database, we found that 17 were overlapped. The visualization of the interplay among TPL, nephrotoxicity and the 17 overlapped targets (OTs) was shown in Figures 2A, B by cytoscape3 (Shannon et al., 2003). These OTs were ported to STRING database to acquire their protein-protein interaction (PPI) pairs ( Figure 2B). Contribution score (CS) was used to assess the importance of interested genes in contributing TPL renal toxicity, it is the number of gene nodes that correlated with each overlapped target (OT) within the PPI network. The contribution scores of 17 OTs were shown in Figure 2C, among which STAT3, TNF and JUN rank the top 3 genes with ≥ 11 gene nodes. By ranking the CSs, genes with CSs ≥ 4 were regarded as candidate targeted genes related to TPL renal toxicity (CTGs-TPL) ( Figure 2C). Potential targets were subjected gene enrichment analysis and the top KEGG terms were shown in Figure 2D. Western-blot (WB) assay was employed to validate the association between the top-ranked CTGs-TPL include Stat3 and c-Jun, it showed that c-Jun were decreased in both L-TPL and H-TPL rats ( Figure 2E), suggesting that c-Jun is a potential target of TPL. In addition, the activation status of c-Jun, the phosphorylated c-Jun (pc-Jun) was also inhibited slightly ( Figure 2E), Considered the vital importance of Stat3 in renal function, two key factors that play essential roles in the upstream and/or downstream of Stat3 signaling include Jak2 and IL17 were selected to investigate their expression levels by WB assays although Stat3 and phosphorylated-Stat3 (pStat3) showed non-significant changes with TPL treatment ( Figure S3). Nevertheless, neither Jak2/phosphorylated-Jak2 (p-Jak2) nor IL17 showed response to TPL treatment ( Figure S2), demonstrating that TPL may not serve as a potential ligand for either Stat3/IL17 or Stat3/Jak2 signaling in the process of TPL nephrotoxicity.
To discovery the structure-based associations between identified targets and TPL, a molecular docking based strategy was applied to predict the ligand-target interactions between TPL and interested targets include CD86, IL4, CXCL8, STAT3 and CD40. Their interactions were visualized in Figure 3.

Dysregulated protein coding genes
As H-TPL rats showed aggravated renal pathological changes with wider intercellular gap between renal tubules, we selected H-TPL renal tissues instead of L-TPL for deep RNA-seq to elucidate the underlying molecular mechanisms of TPL induced renal toxicity. A total of 178 up-regulated and 152 downexpressed PCGs were obtained (Table S2). The gene enrichment analysis (GEA) indicates that up-regulated PCGs are classified (rich factor >  (hsa04977), Nitrogen metabolism (hsa00910), glycine, serine and threonine metabolism (hsa00260), steroid biosynthesis (hsa00100), citrate cycle (TCA cycle) (hsa00020) and proximal tubule bicarbonate reclamation (hsa04964) ( Figure 4A). The downexpressed PCGs were significantly enriched in circadian rhythm (hsa04710) signaling (rich factor > 10) ( Figure 4B). Co-expression correlation of protein-protein pairs was calculated by WGCNA. The networks of protein-protein interaction (PPI) were shown in Figure S4, the core genes include Ptcd3 and Cdk1.

10) to vitamin digestion and absorption
Next to liver, kidney exerts the second most robust rhythms of circadian gene expression Myung et al., 2019). Target genes predicted to be associated with triptolide nephrotoxicity by network pharmacology-based analysis. (A) Visulization of interpalys between TPL and its predicted targets; (B) PPI network of TPL renal toxicity candidate target genes; (C) List of top 15 target genes ranked by contribution scores and they were regarded as candidate targeted genes related to TPL renal toxicity (CTGs-TPL); (D) KEGG terms of potential targets by gene enrichment analysis; (E) Western-blot assay of c-Jun and pc-Jun in Ctrl, L-TPL and H-TPL groups (from left to right, n = 3/group).
Among these dysregulated PCGs in H-TPL rats, we noticed that multiple circadian genes such as Per1, Per2, Per3 and Cry were significantly downexpressed. Previous study found that Per1 in kidney is important for renal sodium handling and necessary for maintaining homeostasis (Douma et al., 2022); therefore, Per1 was selected and validated by western blot assay, it showed that Per1 was decreased in TPL treated kidney ( Figure 4C). Although the roles of circadian genes are unknown in the mechanisims of TPL renal toxicity, our results suggested that TPL may relate with the circadian pace of kidney function. It is interesting that we noticed that c-Jun was not among the significant dysregulated genes by TPL treatment, suggesting that c-Jun may involve in TPL toxicity in renal tissues via posttranscriptional regulation.

Significantly differentially expressed lncRNAs
After a strict filtering pipeline, a total of 6061 novel lncRNAs were identified and combined with the annotated lncRNAs of rat genome (ALRG) for next-step analysis. The length distribution and exon number density plots were shown in Figure 5A, B and Figure  S5. The majority of novel lncRNAs and ALRG own 2 exons. Unlike ALRG that generally enriched in 200 -500bp, the length of novel lncRNAs are mostly distributed in 200 -500, 500 -1000 and > 3500 bp ( Figures 5A, B). A total of 131 up-and 119 down-expressed lncRNAs were identified as significantly differentially expressed lncRNAs (SDElncs) in H-TPL (p-value< 0.05). Increasing studies have demonstrated that lncRNAs usually owns the capacity to regulate their nearby genes (Statello et al., 2021). To illustrate the potential roles of SDElncs, genes that locate within 100kb of SDElncs were acquired and own strong correlations with SDElncRNAs (weight value > 0.8) were defined as target genes ( Figure 5C). A total of 26 genes were identified and subjected for GEA. We surprisingly found target genes, alike with the dysregulated PCGs, were also enriched in vitamin digestion and absorption (hsa04977) and metabolic related pathways such as Alanine, aspartate and glutamate metabolism (hsa00250) ( Figure 5D), which suggesting that abnormal metabolism of amino acids was potentially related with TPL nephrotoxicity.

CircRNA signatures
A total of 1529 high-quality novel circRNAs were identified. By calculating the SRPBM values, only 7 circRNAs were found deferentially expressed in H-TPL rats with p-value< 0.05 (Table  S1). As circRNAs can be served as miRNA sponges, we utilized miRanda tool (Enright et al., 2003) to predict the potential connections between dysregulated circRNAs and miRNAs. Ranking by tot scores, the top 10 circRNA-miRNA pairs were shown in Table 1. Previous studies reveled that miR-207 was upregulated in renal and urine of rats with renal fibrosis and decreased in ischemia-reperfusion injury (IRI) model of mouse (Wei et al., 2010;Shi et al., 2018). We found that both circ: Chr 6:124934385- 124981303 and circ: Chr11: 66774701-66795896 are strongly targeted with miR-207 (Tot scores > 1000), suggesting these two circRNAs may be involved with mi-207 related renal function regulation although the deep mechanisms is unknown.

Discussion
TPL has been applied as an useful compound for treatment of multiple renal diseases for decades; however its toxicity largely limited its clinical practice. Metabolites has been studied in a broad field such as screening for new therapeutic targets, discovery and validation of disease biomarkers. Multitude studies have applied metabonomics technology to investigate TPL, the regulating mechanisms and the toxicities (Du et al., 2014;; however, the transcriptional changes of TPL nephrotoxicity were rarely reported. In this study, a combined approach of network pharmacology method and RNA-seq was used to elucidate the molecular mechanisms of TPL nephrotoxicity. RNA-seq analysis found that a series of circadian genes, such as Per1-3, were significantly dysregulated in renal tissues along with H-TPL treatment. Per1-3 are closely related with renal rhythm. Per1 acts as a circadian clock transcription factor and was regulated by aldosterone, a steroid hormone increases blood pressure via elevating blood volume and Na+ retention (Douma et al., 2022). Myung et al. demonstrated that the mouse kidney of adenine diet induced chronic kidney disease (CDK) model displayed disorganization of Per2 expression (Myung et al., 2019). Per3  exerts dynamic expression patterns in pan renal carcinoma . Through western blot assay, we validated that Per1 was decreased along with TPL treatment, suggesting that Per1 is involved in the regulation of TPL nephrotoxicity.
For identifying candidate targets of TPL nephrotoxicity by network pharmacology based analysis, Huang et al. employed GeneMANIA database and screened out 39 direct-targets in male rats . Although there is no evidence suggested that TPL toxicity has sexual difference, our study utilized female rats to identify TPL regulated proteins and found that female rats seem to tolerant the renal toxicty under both L-TPL and H-TPL treatment with low renal injure rate. Comparing with Huang's study, a more strict filtering standard and different databases were used and we gained highly consistent targets. Our further investigations by western-blots validated that c-Jun protein is a potential target of TPL. c-Jun protein is a widely expressed transcription factor associated with a variety of diseases include human renal diseases (Blau et al., 2012). In glomerular and tubular cells, c-Jun was activated and its activation involves in the regulation of renal inflammation and/or fibrosis (De Borst et al., 2007).
RNA-seq and network pharmacology are different techniques to elucidate relevant candidate molecular targets from two perspectives. Network pharmacology is a novel approach that widely applied for discovering the targets involved in the process of TCM compounds or modern drugs treatment in a specific disease via integrating biomedical, pharmacological and computational approaches while RNA-seq can gain us a cohort of genes with differential expression directly. A combination of these two techniques definitely provide a more comprehensive knowledge of molecular mechanisms in the process of TPL induced renal toxicity. In this study, Pe1 and c-Jun are two candidates related with TPL nephrotoxicity identified by these two analyses, respectively. Although little evidence has been implied on the connections between Per1 and c-Jun and the mechanisms among c-Jun, Per1, TPL and nephrotoxicity are remain explored, our study suggested that c-Jun protein and Per1 are possibly to be involved in TPL induced renal toxicity via two independent pathways.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics statement
The animal study was reviewed and approved by Shenzhen PKU-HKUST Medical Center.