Transcriptomics Analysis of Candida albicans Treated with Huanglian Jiedu Decoction Using RNA-seq

Candida albicans is the major invasive fungal pathogen of humans, causing diseases ranging from superficial mucosal infections to disseminated, systemic infections that are often life-threatening. Resistance of C. albicans to antifungal agents and limited antifungal agents has potentially serious implications for management of infections. As a famous multiherb prescription in China, Huanglian Jiedu Decoction (HLJJD, Orengedokuto in Japan) is efficient against Trichophyton mentagrophytes and C. albicans. But the antifungal mechanism of HLJDD remains unclear. In this study, by using RNA-seq technique, we performed a transcriptomics analysis of gene expression changes for C. albicans under the treatment of HLJDD. A total of 6057 predicted protein-encoding genes were identified. By gene expression analysis, we obtained a total of 735 differentially expressed genes (DEGs), including 700 upregulated genes and 35 downregulated genes. Genes encoding multidrug transporters such as ABC transporter and MFS transporter were identified to be significantly upregulated. Meanwhile, by pathway enrichment analysis, we identified 26 significant pathways, in which pathways of DNA replication and transporter activity were mainly involved. These results might provide insights for the inhibition mechanism of HLJDD against C. albicans.


Introduction
Candida albicans is the most prevalent opportunistic fungal pathogen implicated in superficial mucosal infections as well as invasive disseminated infections, especially in immunocompromised patients [1,2]. C. albicans infections are usually treated with antifungal agents, such as azoles, echinocandins, and polyene drugs. Limited by the number of available antifungal targets, the antifungal agents still remain restricted. The azoles are the most widely used drugs for treating pathogenic fungal infections. Sterol 14 -demethylase (ERG11) is an ancestral activity of the cytochrome P450 superfamily, which is required for ergosterol biosynthesis in fungi and cholesterol biosynthesis in mammals [3]. As a key enzyme of sterol biosynthesis, Erg11 is the main target for therapeutic azole antifungal drugs [4,5].
Widespread overuse of azole drugs for decades has led to the occurrence of drug-resistant isolates [6][7][8]. The prolonged and repeated treatment of OPC (oropharyngeal candidiasis) in AIDS patients has resulted in an increasing frequency of therapy failures caused by the emergence of fluconazoleresistant C. albicans strains. In one study, the levels of fluconazole resistance of a series of 17 clinical isolates taken from a single HIV-infected patient who was treated with azoles over 2 years increased over 200-fold [9]. In recent years, the incidence of azole-resistant strains of C. albicans has increased, especially the rapid emergence of fluconazoleresistant strains. In the vast majority of countries, far less than 10% of C. albicans strains isolated from 1997 to 2001 are resistant to fluconazole [10]. But recent study in China showed that the rate of fluconazole resistance in C. albicans was almost 14.1% [11]. In USA, compared with 2008, the proportion of cases identified from 2008 to 2013 from Georgia and Maryland with fluconazole resistance decreased (GA: 8.0% to 7.1%, −10%; MD: 6.6% to 4.9%, −25%), but the proportion of cases with an isolate resistant to an echinocandin 2 Evidence-Based Complementary and Alternative Medicine increased (GA: 1.2% to 2.9%, +147%; MD: 2.0% to 3.5%, +77%) [12]. So far, several resistance mechanisms of C. albicans have been well characterized: alterations in the sterol biosynthesis pathway, mutations in the ERG11 gene encoding the drug target enzyme, overexpression of the ERG11 gene, and overexpression of genes encoding efflux pumps [13]. Resistance of C. albicans to antifungal agents and limited antifungal agents has potentially serious implications for management of infections.
As a famous multiherb prescription in China, Huanglian Jiedu Decoction (HLJJD, Orengedokuto in Japan) is an aqueous extract of 4 herbal materials, Coptidis Rhizoma, Scutellariae Radix, Phellodendri Cortex, and Gardeniae Fructus with the ratio of 3 : 2 : 2 : 3. HLJJD was first mentioned in the book Wai-Tai-Mi-Yao compiled by Wang Tao in the Tang dynasty (about 752 AD), and it has been widely used in the clinical practice in China and officially listed in the Chinese Pharmacopoeia [14]. HLJJD has been widely used in the treatment of gastrointestinal disorders, inflammation and cardiovascular diseases, and Alzheimer's disease in China [15][16][17]. Modern pharmacological research also demonstrated multiple biological activities of HLJDD: decreasing levels of plasma glucose and blood lipid in type 2 diabetes mellitus [18,19]; increasing the cerebral blood flow, inhibiting the platelet aggregation; reducing hypertension and altering the gene expression profiles of spontaneous hypertensive rats [20][21][22][23]; reducing hepatic triglyceride accumulation, restraining the preadipocyte differentiation and lipid accumulation, inhibiting the lipid peroxidation, and preventing atherosclerosis [15,[24][25][26]. Anti-inflammatory effects of HLJDD were also investigated in some papers [18,[27][28][29]. Moreover, in Mugil cephalus, 1% modified HLJDD feeding for 28 days may be an optimal dose to prevent Lactococcus garvieae infection and could be used in aquaculture industries. In in vitro study, the modified HLJDD also activated the plasma bactericidal activities [30].
Each herb of HLJDD contains many chemical components. Some papers had reported the content determinations of components contained in HLJDD [31,32]. HLJDD contains multiple bioactive secondary metabolites, mainly including alkaloids from Coptidis Rhizoma and Phellodendri Cortex, flavonoids from Scutellariae Radix, and terpenes from Gardeniae Fructus. There are 4 typical compounds from HLJDD: geniposide, baicalin, berberine, and baicalein [27]. Further study showed that the combination of fluconazole and baicalein or berberine produced potently synergistic action in vitro, while baicalein and berberine showed weak antifungal activity when they were tested alone [33,34]. Our preliminary work showed that HLJDD is efficient against Trichophyton mentagrophytes and C. albicans [35]. HLJDD showed its impressive antifungal effect by multitarget and multichannel actions, due to the multiple components. For that reason, the use of HLJDD may be more beneficial to human health in fungal infection treatment as diverse mechanisms showed complementary effects between herbs. But the antifungal mechanism of HLJDD remains unclear.
RNA-seq (deep-sequencing of cDNA) has been used successfully to identify and quantify gene expression at a genome scale level under different conditions or in different cell types. Moreover, it is significantly more sensitive than microarray hybridization approaches [36]. This approach has already been used in C. albicans to generate a high-resolution map of the C. albicans transcriptome under several different environmental conditions [37]. The effect of berberine chloride on Microsporum canis infection was analyzed by the construction of a transcriptome of the M. canis cellular responses upon berberine treatment [38]. Therefore, in this study, by using RNA-seq technique, we performed a largescale analysis of gene expression changes when C. albicans was exposed to HLJDD, to better understand how HLJDD inhibits the growth of C. albicans.

Preparation of the Extract of HLJDD.
The herbal medicines of modified HLJDD were dried at 40 ∘ C for 24 h and then pulverized to powder using a mechanical blender. 0.5%, 1%, and 2% w/w of the powder was prepared and boiled for 30 min with 200 mL of deionized water, and the aqueous extracts were filtered through Whatman number 1 filter paper. The HLJDD residues were also boiled with another 200 mL of deionized water.

Determination of Sensitivity of the SC5314 Strain to
HLJDD. Antifungal susceptibility testing was performed by using the CLSI M27-A3 microbroth dilution method [40]. MICs were determined after growth at 30 ∘ C for 24 h for HLJDD. MICs were read as the lowest drug concentration producing a prominent decrease in turbidity translating to 100% growth reduction compared with the drug-free control.

Total RNA Extraction.
To identify genes in the early response of C. albicans to HLJDD, we treated the isolate with HLJDD at 20 mg/mL, the lowest drug concentration producing a prominent decrease in turbidity translating to 100% growth reduction compared with the drug-free control determined above. To extract total RNA, the cells of SC5314 were inoculated into YPD medium and cultured at 30 ∘ C overnight. Before SC5314 were harvested for RNA extraction, the culture was treated with HLJDD at 20 mg/mL for 3 h. The untreated culture was used as the control. Total RNA was isolated according to the protocol described by Alison et al. [41].

RNA Sequencing and Assembling.
Three independent experiments were performed for the study of either control C. albicans or C. albicans with HLJDD treatment. Shear cDNA into 300-500 bp fragments using ultrasonic apparatus (Fisher) and purify it with Ampure beads (Agencourt, America). Library of all the samples was constructed according to the procedure of NEBNext5 UltraTM RNA Library Prep Kit for Illumina (NEB, America). Sequencing library was checked with Onedrop quantitation, 2% agarose gel electrophoresis detection, and high sensitivity of DNA chip detection. Paired-end sequencing of cDNA was carried out with Illumina Hiseq TM2000. Raw data was filtered by removing reads with adaptor sequences, as well as low quality reads. Then, clean reads were obtained and mapped to reference sequences using SOAP (2.21) [42].

Gene Prediction and Annotation.
Trinity software was used to assemble the clean reads into contigs and BLAST (2.2.23) was used to do gene prediction. Predicted sequences ( -value < 1.0 − 05) were annotated with information from GenBank NR, GO, and KEGG using BLAST2GO (2.2.5). GO classification was conducted using WEGO [43].

Analysis of Differential Expressed
Genes. The expression level for each gene is determined by the numbers of reads uniquely mapped to the specific gene and the total number of uniquely mapped reads in the sample. The gene expression level is calculated by using RPKM (Reads Per kb per Million reads) method [36]. Then, NOI seq method was applied to screen differentially expressed genes between two groups, with the threshold of significance as fold change of RPKM ≥ 3 and Probability ≥ 0.8 [44].

Enrichment Analysis of GO and KEGG Pathways.
Enrichment analysis was performed by hypergeometric test to find significantly enriched GO terms and KEGG pathways in DEGs. False discovery rate (FDR) of pathways was calculated. The threshold of significance of pathways was set as FDR < 0.05.

Real-Time Quantitative Reverse Transcription-(qRT-) PCR.
To evaluate the validation of RNA-seq results, we conducted quantitative real-time (RT) PCR assays for determination of expression of 8 genes. Gene expression levels were calculated using the 2 −ΔΔCt method [45]. For each sample, PCR amplifications with primer pair actin-F and actin-R for the quantification of expression of actin gene were performed as a reference. The experiment was repeated 3 times.

RNA Sequencing and Gene Prediction.
Approximately 12,000,000 raw reads were obtained from each sample. After filtering by quality, about 96% clean reads were mapped. Summary of mapping result was shown in Table 1. Using the longest sequence of a subgroup as the unigene as the reference sequence, we got 6057 predicted protein-encoding genes totally. The data have been submitted to NCBI under BioProject accession number PRJNA314910.

Identification and Verification of Differentially Expressed
Genes. By using the threshold of significance as fold change of RPKM ≥ 3 and Probability ≥ 0.8, we obtained a total of 735 differentially expressed genes (DEGs), including 700 upregulated genes and 35 downregulated genes (Supporting Information

Effects of HLJDD Treatment on the Genes Involved in
Sterol Biosynthesis. As the most widely used antifungal drugs, azoles can block fungal sterol biosynthesis pathway. Thus, effects of HLJDD on the genes involved in sterol biosynthesis were analyzed in detail. Expression of 23 genes involved in sterol biosynthesis was detected in the RNA-seq analysis, and expression of 8 genes showed a more than 2-fold increase after C. albicans was treated with HLJDD; only the genes encoding sterol 24-C-methyltransferase (ERG6) and C-8 sterol isomerase (ERG2) were upregulated by more than 3 times ( Table 3). None of these 23 genes was downregulated significantly (Probability > 0.8) by HLJJD.

Effects of HLJDD Treatment on the Genes Encoding Multidrug Transporters.
In C. albicans, upregulation of multidrug transporter genes is one of the well-documented mechanisms of resistance to azole antifungal agents [9,[46][47][48]. Two families of multidrug transporters, the ABC (ATP-binding cassette) transporter family (Cdr1p and Cdr2p) and the major facilitator superfamily (MFS, CaMdr1p), have been shown to be involved in resistance to azole antifungal agents [47,48]. Thus, we also paid attention to the multidrug transporter genes. In genome sequences of C. albicans, a total of 36 genes are annotated as multidrug transporters. In this study, expression of 32 genes was detected by the RNA-seq, and 7 genes were identified to be significantly upregulated more than 3 times by HLJDD treatment (Table 4), including CDR2 (Candida Drug Resistance) from the family of ABC transporters. Cdr2 has been shown as the principal mediators of resistance to azoles due to transport phenomena [47,48].

Enrichment Analysis of GO and KEGG Pathways.
GO and KEGG assignments were used to classify the genes in the response of C. albicans to HLJDD. By GO classification analysis, the percentage and distribution of top-level GO terms were portrayed in the 3 categories: (A) cellular component; (B) molecular function, and (C) biological process  Figure 1). A high percentage of genes were assigned to "cell," "cell part," "binding," "catalytic," "cellular process," and "metabolic process" (Figure 1). By enrichment analysis, with FDR < 0.05, 23 significant GO terms and 3 significant KEGG pathways were identified (Supporting Information Table S3). These significant pathways were mainly associated with DNA replication and transporter activity. The maps with highest unigene representation were meiosis (cal04113; 23 unigenes), followed by cell cycle (cal04111; 23 unigenes), and DNA replication (cal03030; 11 unigenes).

Discussion
C. albicans is the most prevalent opportunistic fungal pathogen causing superficial to systemic infections in immunocompromised individuals [1,2]. The concomitant use of drugs and the lack of available drugs frequently result in the occurrence of drug-resistant isolates and strains display multidrug resistance (MDR). In search of novel fungicides, efficiency of medicinal plants against fungi has been reported, but studies on their underlying mechanisms are very few [49]. In this study, we explored a famous multiherb prescription in China, Huanglian Jiedu Decoction (HLJJD, Orengedokuto in Japan), for its antifungal potential. Our preliminary work showed that HLJDD is efficient against Trichophyton mentagrophytes and C. albicans [35]. HLJDD showed its impressive antifungal effect by multitarget and multichannel actions, but studies on the underlying mechanisms are very few. To determine the antifungal mechanism of HLJDD against C. albicans, we performed a large-scale analysis of gene expression changes when C. albicans was exposed to HLJDD, to better understand how HLJDD inhibits the growth of C. albicans. Due to the multiple components of HLJDD, it is most likely that the antifungal effect is multitarget and multichannel actions. KEGG analysis suggested that 3 cellular functions were affected in C. albicans upon HLJDD treatment, including meiosis, cell cycle, and DNA replication. Most genes (56 genes) involved in the 3 cellular functions were upregulated excepted for 1 gene, potential hexose transporter (XP 719596.1). Among these genes, Spo22 (also called Zip4) (XP 718811.1) was upregulated obviously upon HLJDD treatment. Zip4/Spo22 was shown to be a central protein of the SICs (synapsis initiation complexes), from which the polymerization of the transverse filament proceeds. In S. cerevisiae, Zip4/Spo22 was identified as a member of the ZMM group of proteins that also includes Zip1, Zip2, Zip3, Msh4, Msh5, and Mer3 which together control the formation of class I COs [50][51][52]. In Arabidopsis thaliana, Zip4/Spo22 function in class I CO formation is conserved with budding yeast. However, mutation in AtZIP4 does not prevent synapsis, showing that both aspects of the Zip4 function (i.e., class I CO maturation and synapsis) can be uncoupled [51].
Azoles are the most widely used antifungal drugs, which target on cytochrome P450 lanosterol 14 -demethylase encoded by the ERG11 gene. In Fusarium graminearum, using a deep serial analysis of gene expression (DeepSAGE) sequencing approach, the transcriptional response of F. graminearum to tebuconazole (a widely used azole fungicide) was profiled. Expression of 23 genes involved in sterol biosynthesis was detected in the DeepSAGE analysis, and expression of 9 genes showed a more than 5-fold increase after the fungus was treated with tebuconazole. None of these 23 genes was downregulated by more than 5 times by tebuconazole [53]. Thus, effects of HLJDD on the genes involved in sterol Evidence-Based Complementary and Alternative Medicine 5   biosynthesis were analyzed in detail. Expression of 23 genes involved in sterol biosynthesis was detected in the RNA-seq analysis, and expression of 8 genes showed a more than 2-fold increase after the fungus was treated with HLJDD, only the genes encoding sterol 24-C-methyltransferase (ERG6) and C-8 sterol isomerase (ERG2) were upregulated by more than 3 times ( Table 3). None of these 23 genes was downregulated significantly (Probability > 0.8) by HLJJD. These results indicate that HLJDD might also affect sterol biosynthesis of C. albicans.
Overexpression of multidrug resistance efflux transporter genes in several fungi was found to be correlated with azole resistance [54]. In C. albicans, a number of efflux transporter genes have been cloned and characterized. Two families of multidrug transporters, the ABC (ATP-binding cassette) transporter family (Cdr1p and Cdr2p) and the major facilitator superfamily (MFS, CaMdr1p), have been shown to be involved in resistance to azole antifungal agents [47,48]. Expression of 32 genes out of 36 genes annotated as multidrug transporters in genome sequences of C. albicans was detected by the RNA-seq sequencing. Expression of 13 genes was upregulated by more than 2 times by HLJDD; meanwhile, only 2 genes were significantly downregulated including CDR4. In addition, expression of only 4 genes was upregulated by more than 3 times by HLJDD, including CDR2, which plays an important role in azole resistance ( Table 4). The upregulated expression of these genes may be related to efflux of HLJJD, which provides supporting evidence to previous studies on expression level.
Previous study examined changes in the gene expression profile of C. albicans following exposure to representatives of the 4 currently available classes of antifungal agents, the azoles (ketoconazole), polyenes (amphotericin B), echinocandins (caspofungin), and nucleotide analogs (5flucytosine). And the data showed that none of the differentially regulated genes found exhibited similar changes in expression for all 4 classes of drugs. Thus, the response of C. albicans to different drugs seems to be highly specific [55]. Ketoconazole exposure increased the expression of genes involved in lipid, fatty acid, sterol metabolism, and several genes associated with azole resistance, including CDR1 and CDR2 [56]. It is surprising that HLJDD increased the expression of genes involved in sterol metabolism and azole resistance (CDR1 and CDR2). Considering the similarity of expression changing pattern, it is possible that HLJDD affects sterol metabolism. And further experiments are required to confirm this hypothesis.

Conclusions
In conclusion, we performed a transcriptomics analysis of gene expression changes for C. albicans under treatment of HLJDD using RNA-seq technique. Overall, a total of 6057 predicted protein-encoding genes were identified. Further gene expression analysis revealed a total of 735 differentially expressed genes (DEGs), including 700 upregulated genes and 35 downregulated genes. Intensive bioinformatics analysis identified 26 significant pathways, and DNA replication and transporter activity were mainly involved.
In addition, genes encoding multidrug transporters such as ABC transporter and MFS transporter were identified to be significantly upregulated. Overall, the results from this study might provide insights in understanding of the mechanisms for the response of C. albicans to HLJJD. Furthermore, this work demonstrates the potential utility of the RNA-seq technique in antifungal studies.