Conprehensive Analysis of circRNA-microRNA-mRNA Network Revealed Novel Regulatory Mechanism in Trichosposon Asahii Infection

Background: Trichosporon asahii (T. asahii) invasive infection frequently occurs in immunodecient hosts with high mortality, but the pathogenesis of T. asahii infection remains elusive. Circular RNAs (circRNAs) are a type of endogenous noncoding RNAs that participant various disease processes. However, the mechanism of circRNAs in T. asahii infection are still completely unknown. Methods: RNA sequencing (RNA-seq) was performed to analyse the expression proles of circRNA, microRNA (miRNA), and mRNA in THP-1 cells infected with T. asahii or uninfected samples. Part of the RNA-seq results were veried by RT-qPCR. Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were used to analyse the differentially expressed (DE) mRNA. The circRNA-miRNA-mRNA network was constructed and veried with dual-luciferase reporter assay and overexpression experiment. Results: A total of 46 circRNAs , 412 mRNAs and 47 miRNAs were DE after T. asahii infection at 12 h. GO and KEGG analysis showed that the DE mRNAs were primary linked to the leukocyte migration involved in inammatory response, Toll-like receptor signaling pathway, and TNF signaling pathway. A competing endogenous RNA (ceRNA) network was constructed with ve DE circRNAs, ve DE miRNAs and 42 DE mRNAs. Among them, we veried that hsa_circ_0065336 indirectly regulate PTPN11 expression by sponging miR-505-3p. Conclusions: These data revealed a comprehensive circRNA-associated ceRNA network during T. asahii infection, thus providing new insights to clarify the pathogenesis between T. asahii-host interation.


Background
Trichosporon asahii (T. asahii) is the most common pathogenic species in the genus trichosporon that widely distributed in tropical and temperate regions, usually through built-in catheters, intestinal mucosa translocation of microorganisms or respiratory inhalation to invade the human body [1]. T. asahii are usually associated with super cial mycosis in immunocompetent hosts but could cause invasive infections in immunosuppressed patients. With the abuse of immunosuppressants and broad-spectrum antibiotics and the extensive development of traumatic surgery, the incidence of invasive trichosporonosis has increased over the past years [2]. Despite that various antifungal drugs have been used in the treatment of invasive trichosporonosis, such as amphotericin B, ucytosine, caspofungin and uconazole, the mortality of invasive trichosporonosis is still high and requires further analysis [3]. Circular RNAs (circRNAs) were a type of endogenous noncoding RNAs that are widely distributed in eukaryotic cells and have the function of regulating gene expression. They are constituted by a covalently closed cyclic structure with high abundance, relative conserved, tissue-speci c and cell-speci c expression patterns [4]. Although circRNAs have been discovered for more than 40 years, they were mainly regarded as 'junk' generated by mis-splicing of transcripts [5,6]. In recent years, with the wide application of RNA sequencing (RNA-seq) technology and the rapid development of circRNA-speci c computational tools, various circRNAs have been reported to be involved in the pathogenesis of cancer, neurological disorders, cardiovascular diseases, diabetes mellitus and autoimmune diseases [7][8][9][10][11].
However, circRNAs associated with fungal infections have not been reported up to date.
MicroRNAs (miRNAs) are small noncoding RNAs that can negatively control their target gene expression post-transcriptionally [12,13]. Recent studies have characterized the miRNAs expression pro les and identi ed several critical miRNAs that participate immune and in ammatory responses following fungal exposure [14]. However, the upstream regulation mechanism of miRNAs are poorly understood. It has been reported that circRNAs can regulate gene expression by binding miRNA response elements (MREs) called competing endogenous RNAs (ceRNAs) mechanism [15]. But whether circRNAs could serve as ceRNAs and involve in the interactions between T. asahii and its host remain unknown.
In this study, we systematically investigated the alteration of circRNAs, miRNAs, and mRNAs expression pro les in macrophages during T. asahii infection by using RNA-sEq. A large number of dysregulated circRNAs, miRNAs and mRNA were identi ed after T. asahii infection. CeRNAs network revealed the putative regulatory roles of circRNAs in cellular response to T. asahii infection. Furthermore, we identi ed that hsa_circ_0065336 could indirectly regulate PTPN11 expression by sponging miR-505-3p. In summary, this is the rst report of the expression pro le and function analyses of circRNAs after fungal exposure in mammalian cell.

RNA extraction, library construction and sequencing
Total RNA was extracted from both mock-infected and T. asahii-infected cells at 12 h time point by TRIzol reagent (Invitrogen, USA) according to the manufacturer's protocol. The concentration and purity of total RNA was analyzed by NanoDrop ND-2000 (Thermo Fisher Scienti c, USA) and the integrity was evaluated by Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA). We constructed two types of libraries from six samples for sequencing as follows: (1) rRNA-depleted RNA library: The rRNA was removed in total of RNA (5 µg) using a Ribo-off rRNA Depletion kit (Vazyme, China). After rRNA-depleted RNA puri cation, cDNA libraries were constructed using VAHTSTM Stranded mRNA-seq Library Prep Kit for Illumina® kits (Vazyme, China) according to the manufacturer's instructions. The Illumina Hiseq XTen platform (Sangon Biotech, China) was used to perform sequencing analysis. (2) miRNA library: Total RNA was collected and quanti ed by Qubit2.0 RNA detection kit, The 3′ adapter and 5′ adapter were ligated by using T4 RNA Ligase 2 (New England Biolabs, USA) and T4 RNA Ligase 1 (New England Biolabs, USA), respectively.
Reverse transcription and PCR ampli cation of ligation products. The 140 ~ 150 bp PCR products were enriched and quanti ed with Qubit 2.0 DNA kit. Finally the cDNA library was sequenced on the Illumina Hiseq XTen platform (Sangon Biotech, China).

circRNA identi cation and characterization analysis
The raw sequencing data were ltered by Trimmomatic software to obtain high quality clean data [16]. the quality control sequences were mapped to the human reference genome hg38 using BWA-MEM [17], Then subjected to CIRI2 to identify circRNAs [18]. The BEDtools was used to determine the source of circRNAs based on circRNAs position and gene annotation information [19].
Abundance quantify and DE analysis mRNAs abundance were detected by StringTie software and normalized to transcripts per million (TPM) [20]. The formula for calculating TPM is as follows: . The DESeq2 was used for DE analysis of mRNAs [21]. A q-value < 0.01, |log 2 FoldChange (FC)| > 1 and Mean TPM > 5 in at least one group were set as thresholds for signi cant DE mRNAs. circRNAs expression were quanti ed based on the number of back-spliced junction read pairs and circRNAs length were calculated using Reads Per Kilobase of transcript per Million mapped reads (RPKM ). The formula is shown as follows: . |FC| > 1.5 and P < 0.05 were considered signi cantly DE. The miRNA level was analyzed and normalized counts to reads per million (RPM). The RPM was calculated as below: . The miRNAs were considered DE only when the | log 2 FC| > 1 and P < 0.05.

RT-qPCR validation
To verify the results of RNA-seq, we performed RT-qPCR analysis on some of the DE circRNAs, miRNAs and mRNAs. The total RNA was extracted from T. asahii infected and uninfected cells, and reverse transcribed to cDNA using a PrimeScript™ RT reagent Kit with gDNA Eraser (Perfect Real Time) (Takara, Japan) according to the manufacturer's instructions. qPCR reactions were performed on ABI StepOnePlusTM Real-Time PCR System (Applied Biosystems, USA). The divergent primers encompassing the back-splicing junctions for circRNAs and the convergent primers for mRNAs and miRNAs were synthesized by Sangon Biotech (China), and the sequences were shown in (Additional le 1: Table S1). βactin was used as an internal control for circRNAs and mRNAs, and U6 served as the endogenous control for miRNAs. We used the 2 −△△CT method to analyze the data. Each experiment was repeated in triplicate.
GO and KEGG pathway analysis GO enrichment and KEGG pathway analysis were performed to determine the role of the DE mRNAs and ceRNA network-associated target mRNAs. P-value < 0.05 is a threshold indicate signi cant difference.
ceRNA network analysis The circRNA-miRNA-mRNA network was constructed as the following standards: (1) DE circRNAs, miRNAs and mRNAs; (2) The expression of circRNA/mRNA and miRNA must have the opposite trend; (3) Pearson correlation coe cient (r) > 0.8 was set as a threshold between circRNA and mRNA; (4) the coexpressed circRNA-mRNA pairs should have the same putative miRNA binding sites based on miRanda [22]. Cytoscape software (v3.8.0) was used to construct and visualize the ceRNA network [23].

Cell transfection
The sequence of hsa_circ_0065336 was cloned into overexpression vector pCD5-ciR (Geneseed,

Statistical analysis
The measured data were presented as the means ± standard deviation (SD). Statistical signi cance of different groups were examined using Student's t-test. P-value < 0.05 was de ned as statistically signi cant.

Identi cation and characterization of circRNAs
Six rRNA-depleted RNA libraries from T. asahii-infected and uninfected samples at 12 h were sequenced to analyse the pro les of circRNAs using RNA-sEq. The results showed that a total of 10539 unique circRNAs were identi ed, among which 4595, 3315 and 4456 circRNAs in three control samples, while 3422, 4423 and 4091 circRNAs in three T. asahii-infected samples (Fig. 1a), respectively. According to the position of circRNA on the genome, we found exonic circRNAs were the most abundant circRNAs with a number of 10388 (98.57%), whereas intronic and intergenic circRNAs accounted only for 151 (1.43%) (Fig. 1b). 82.69% of circRNAs were less than 2 k bp and the average length was about 1.54 kb (Fig. 1c).

Identi cation of DE circRNAs
To investigate the DE pro le of circRNAs in response to T. asahii infection, we employed hierarchical clustering analysis, the results indicated that the circRNA expression pro les were signi cantly changed after T. asahii infection. A total of 46 circRNAs were signi cantly dysregulated in THP-1 at 12 h post infection, with 19 circRNAs up-regulated and 27 circRNAs down-regulated (Fig. 2a, 2b and Additional le 1: Table S2). To verify the credibility of RNA-seq results, we randomly selected four signi cantly changed circRNAs (hsa_circ_0077736, hsa_circ_0065336, hsa_circ_0108123 and hsa_circ_0029624) for RT-qPCR veri cation. The RT-qPCR results exhibited an expression trend consisted with the RNA-seq results (Fig. 2c).

Identi cation of DE mRNAs
We also examined the change of mRNAs expression in the same condition. Cluster heatmap revealed that the differential expression patterns of mRNA in THP-1 after T. asahii infection compared with controls. A total of 95139 mRNA were identi ed and 412 mRNAs were DE between the two groups, of which, 276 mRNAs up-regulated and 136 mRNAs down-regulated (Fig. 3a, 3b and Additional le 1: Table S3). To further con rm the reliability of RNA-seq, four differentially expressed mRNAs (NECTIN2, MAPK14, CCL1 and CCL2) were randomly selected and validated with RT-qPCR. The results were highly consistent with RNA-seq data (Fig. 3c).

Identi cation of DE miRNAs
We further analyzed the differentially expressed miRNAs from the same six samples. A total of 3218 miRNAs were detected, among which 47 miRNAs were signi cantly dysregulated between the two groups, including 32 miRNAs up-regulated and 15 miRNAs down-regulated (Fig. 4a, 4b and Additional le 1: Table  S4). Four differentially expressed miRNAs (miR-4792, miR-193b-3p, miR-23b-3p and miR-133a-3p) were randomly selected and veri ed by RT-qPCR. As expected, the RT-qPCR results were consistent with that of RNA-seq (Fig. 4c).

GO enrichment and KEGG analysis of DE mRNAs
We performed GO enrichment and KEGG pathway analysis of the DE mRNAs. The GO enrichment analysis results showed that the DE mRNAs were highly enriched in some biological processes, such as leukocyte migration involved in in ammatory response, innate immune response in mucosa, leukocyte migration and phagocytosis (Fig. 5a). In addition, KEGG pathway analysis indicated that the plenty of mRNAs were primarily involved in Toll-like receptor signaling pathway, phagosome, endocytosis and TNF signaling pathway (Fig. 5b). These data implied that the cellular antifungal response was activated by T. ashiii infection.

Analysis of circRNA-miRNA-mRNA regulatory network
Recent studies showed that circRNAs can indirectly regulate gene expression by act as a miRNA sponge through a ceRNA mechanism. So we constructed a circRNA-miRNA-mRNA network according to the ceRNA hypothesis. The screening criterias are listed in the Methods. 44 cirRNA-miRNA-mRNA pathways were established, including ve circRNAs, ve miRNAs, and 42 mRNAs ( Fig. 5c and Additional le 1: Table  S5). Owing to there are so many mRNAs involved in the ceRNA network, we performed GO enrichment and KEGG pathway analysis of these above mRNAs. The GO enrichment analysis results showed that the DE mRNAs in the ceRNA network were mainly involved in leukocyte migration involved in in ammatory response, transmembrane receptor protein tyrosine kinase signaling pathway and phagocytosis (Fig. 5d). Furthermore, KEGG pathway analysis indicated that neurotrophin signaling pathway, lysosome, and natural killer cell mediated cytotoxicity were mostly involved (Fig. 5e).

Veri cation of ceRNA network
In order to demonstrate the authenticity of ceRNA network, we selected the predicted hsa_circ_0065336/miR-505-3p/PTPN11 pathway based on the literature and GO enrichment and KEGG pathway analysis in the previous step. We rst measured the expression level of hsa_circ_0065336, miR-505-3p and PTPN11 in THP-1 post infected with T. asahii at 0 h, 3 h, 6 h, 9 h, and 12 h time points, respectively. The expression of hsa_circ_0065336 and PTPN11 were up-regulated in a time-dependent manner, and reached a peak at 12 h. Whereas miR-505-3p was down-regulated starting at 3 h after infection, then maintained declining level at 9 h post infection, which is highly similar to the results of RNA-seq.
Bioinformatics prediction analysis showed that miR-505-3p targeted both hsa_circ_0065336 and PTPN11 based on miRanda database. In order to con rm that, the full-length of hsa_circ_0065336-WT and MUT were constructed and subcloned into dual-luciferase report vector pmir-GLO, then co-transfected with miR-505-3p mimics or NC into 293T cells. The results showed that hsa_circ_0065336-WT with miR-505-3p mimics could signi cantly decrease the luciferase activity, whereas hsa_circ_0065336-MUT was not affected. As well, the luciferase activity of PTPN11 3'UTR-WT was signi cantly decreased by miR-505-3p mimics compared with control groups.
To further veri ed the hsa_circ_0065336/miR-505-3p/PTPN11 pathway in THP-1 cells after T. asahii infection, we transfected the hsa_circ_0065336 over expression vector into THP-1 cells and subsequently infected with T. asahii for 12 h. The results indicated that overexpression of hsa_circ_0065336 could signi cant elevated PTPN11 level compared with control. In addition, the increase of PTPN11 induced by hsa_circ_0065336 overexpression could be reversed by miR-744-3p mimics. Taken together, these data identi ed the hsa_circ_0065336/miR-505-3p/PTPN11 pathway in THP-1 upon T. asahii infection.

Discussion
Recently, several studies have reported the cricRNA expression pro les and ceRNA network in viral and bacteria infection diseases [24][25][26][27]. Nevertheless, the expression pro les and potential roles of circRNAs in fungi infections have not been explored.
In this study, we systematically investigated circRNA, miRNA, and mRNA expression pro les in T. asahiiinfected THP-1 cells compared with uninfected samples. The characteristic of circRNA expression, type, length and host gene number were similar with previous reports [28,29]. We found a total of 46 DE circRNAs, 412 DE mRNAs and 47 DE miRNAs upon T. asahii infection after 12 h, of these mRNAs, GO enrichment and KEGG pathway analysis showed that they were primarily enriched in immune-related pathways, suggesting that T. ashii infection could activate a strong antifungal response, contributing to fungi pathogenic progression.
The ceRNA hypothesis was rst proposed by Salmena et al [30] in 2011 described how mRNAs, transcribed pseudogenes, and lncRNAs "talk" to each other using MREs as letters of a new language. Many studies have identi ed that circRNAs could act as ceRNAs to sponge miRNAs, indirect regulating the expression of genes and participate in the occurrence and development of various diseases, especially cancer. Since no data about circRNA-associated ceRNA have been reported in fungal infection, we discovered a systematic circRNA-miRNA-mRNA regulatory network after T. asahii infection based on the RNA-seq and ceRNA rules. Five circRNAs, ve miRNAs, and 42 mRNAs were screened and constructed the ceRNA regulatory network. Further more, we performed GO and KEGG analysis on the 42 predicted targeting mRNAs among the ceRNA network to explore the function of the ve circRNAs. The results showed that many host immune response pathways were involved, including in ammatory response, phagocytosis and neurotrophin signaling pathway, indicating that circRNAs may be the principal regulators of host response during T. asahii infection.
Moreover, We selected hsa_circ_0065336/miR-505-3p/PTPN11 pathway to veri ed the reliability of ceRNA network using RT-qPCR, dual luciferase reporter assay and overexpression experiment. PTPN11 is a cytoplasmic tyrosine phosphatase that participate the signaling of growth factors, cytokines and hormones. Recent study found that PTPN11 is crucial for the induction of pro-in ammatory cytokines and control candida albicans infection through Dectin-1/PTPN11/Syk/NF-κB pathway [31]. miR-505-3p has been showed to be down-regulated in several autoimmune disease, including primary biliary cirrhosis and in ammatory bowel disease [32,33]. Rafael et al. [34] reported that miR-505-3p was decreased in familial hypercholesterolemia, which could up-regulated several chemokine receptor through target RUNX1 transcription factor in macrophages. hsa_circ_0065336 has been reported expression in brain and several cell types, but no study has assessed its function in any disease or physiological process. Here, we identi ed the hsa_circ_0065336 and PTPN11 expression were simultaneously up-regulated, whereas miR-505-3p expression was down-regulated in T. asahii infected THP-1 cells. Dual-luciferase reporter assay demonstrated that miR-505-3p can direct target hsa_circ_0065336 and PTPN11. Moreover, function assay revealed that PTPN11expression was up-regulated by hsa_circ_0065336 while inhibited by miR-505-3p. Therefore, the hsa_circ_0065336/miR-505-3p/PTPN11 regulatory axis was found and veri ed for the rst upon T. asahii infection.

Conclusions
As far as we know, this is the rst report to conjoint analyze the expression pro les of circRNA, miRNA and mRNA upon fungal infection in mammal cells. We construct a circRNA-associated ceRNA regulatory network, which mainly involved in fungal-host immune response. In conclusion, this study could provide novel insight for clarifying the pathogenesis of T. asahii infection.