Digital Gene Expression Analysis of Microsporum canis Exposed to Berberine Chloride

Berberine, a natural isoquinoline alkaloid of many medicinal herbs, has an active function against a variety of microbial infections including Microsporum canis (M. canis). However, the underlying mechanisms are poorly understood. To study the effect of berberine chloride on M. canis infection, a Digital Gene Expression (DGE) tag profiling was constructed and a transcriptome analysis of the M. canis cellular responses upon berberine treatment was performed. Illimina/Hisseq sequencing technique was used to generate the data of gene expression profile, and the following enrichment analysis of Gene Ontology (GO) and Pathway function were conducted based on the data of transcriptome. The results of DGE showed that there were 8476945, 14256722, 7708575, 5669955, 6565513 and 9303468 tags respectively, which was obtained from M. canis incubated with berberine or control DMSO. 8,783 genes were totally mapped, and 1,890 genes have shown significant changes between the two groups. 1,030 genes were up-regulated and 860 genes were down-regulated (P<0.05) in berberine treated group compared to the control group. Besides, twenty-three GO terms were identified by Gene Ontology functional enrichment analysis, such as calcium-transporting ATPase activity, 2-oxoglutarate metabolic process, valine catabolic process, peroxisome and unfolded protein binding. Pathway significant enrichment analysis indicated 6 signaling pathways that are significant, including steroid biosynthesis, steroid hormone biosynthesis, Parkinson’s disease, 2,4-Dichlorobenzoate degradation, and tropane, piperidine and Isoquinoline alkaloid biosynthesis. Among these, eleven selected genes were further verified by qRT-PCR. Our findings provide a comprehensive view on the gene expression profile of M. canis upon berberine treatment, and shed light on its complicated effects on M. canis.


Introduction
Dermatophytes are pathogenic fungi that have the ability to invade the keratinized structure and infect skin, hair and nails of human and animals [1]. The genera Microsporum canis (M. canis) was found in both man and animal and the species M. canis is zoonotic in nature. Previous studies revealed that M. canis is the major cause of dermatophytosis in pets and rabbits. In total, 21 isolates of M. canis were collected from rabbits with or without skin lesions [1]. Although there are a variety of drugs, including azoles, propylene amine, ring ketone amine, and amorolfine, the overall effect is not stable [2]. Therefore, discovery of novel effective medicine against rabbit dermatosis is urgent. The interest in natural medicine has remarkably increased in recent years due to the side effects of conventional drugs, as well as the emergence of antimicrobial resistance to available drugs. Plants and their derivatives are therefore an important alternative in the development of novel drugs. The bark from the Phellodendron tree has been used as traditional Chinese medicine for thousands of years. The species, Phellodendron amurense, is widely used to treat gastroenteritis, abdominal pain or diarrhea, and various inflammatory diseases including arthritis and dermatophytosis [3]. Recently the ethanol extract from Phellodendron amurense has been found to inhibit the growth of certain kinds of fungi such as Trichophyton mentagrophytes [4] and Staphylococcus [5], and in our previously experiment, the MIC value of berberine against M. Canis was 2 mg/mL. Berberine, a major component of M. Canis, is a natural isoquinoline alkaloid found in medicinal herbs. Previous studies have implied a number of biological activities of berberine, including anti-secretory, anti-inflammatory, anti-bacterial, anti-malarial, anti-mycobacterial [6], anti-tumor and anti-cholesterol activities. Therefore, berberine has been widely used in the treatment of bacterial diarrhea and intestinal parasite infections in many countries.
The digital gene expression (DGE) system is a tag-based transcriptome sequencing approach in the production of technical innovation of RNA deep sequencing technologies, which is mainly used to study the comparative gene expression. When DGE is applied, clean tags strained through raw tags with certain end nuclease should be mapped to the reference database. The degrees of expression of the actual total genes in the sample are defined by the quantity of each mRNA yield from each gene. Therefore, in this study, by using DGE, we performed a large-scale analysis of gene expression changes when M. canis was exposed to berberine, to better understand how berberine inhibits the growth of M. canis.

Fungi growth and berberine hydrochloride treatment
The M. canis strain (CMCC(F)M3h) used in this study was provided by the Institute of Internal Medicine, Chinese Academy of Medical Sciences (Nanjing, China) National Center. The M. canis was cultured in Tryptic Soy Broth (TSB) at 28°C for 96 h. The densities of the suspension were adjusted to primary inoculum of 1.0 × 10 7 CFU/mL in the normal saline solution with total volume of 100 mL in conical flasks. Berberine hydrochloride (BC) (HPLC>98%) was purchased from Yuan Ye Biological technology co., Ltd, Shanghai, China. The BC treated group (three samples) was subject to medicinal stress dose which means 1/2 MIC (minimum inhibitory concentration) dose (1 mg/mL berberine in DMSO), while the control group was fed with DMSO (Purity>99%) which was purchased from Sangon Biotech co., Ltd Shanghai, China only. Each group has three repeats. After 6 h incubation (28°C,150 rpm), the seeds were filtered with gauze and frozen after treatment for RNA isolation and further analysis.

RNA isolation and library preparation for DGE
Total RNA was extracted using Trizol reagent (Invitrogen, CA, USA) following the manufacturer's procedure. Total RNA quantity and purity were analyzed by Bioanalyzer 2100 and RNA 6000 Nano LabChip Kit (Agilent, CA, USA) with RIN number >8.0. Approximately 10 μg of total RNA from each sample was subjected to isolate Poly (A) mRNA by using poly-T oligo attached magnetic beads (Invitrogen). Following purification, the mRNA is fragmented into small pieces using divalent cations under elevated temperature. Then the cleaved RNA fragments were reverse-transcribed to create the final cDNA library in accordance with the protocol for the mRNA-Seq sample preparation kit (Illumina, San Diego, USA), the average insert size for the single-end libraries was 200bp-350bp. Then we performed the single-end sequencing (36 bp) on an Illumina Hiseq2500 at the LC-BIO (Hangzhou, China) following the vendor's recommended protocol. Data has been deposited to the server of National Institute of Health (NIH) (http://trace.ncbi.nlm.nih.gov/Traces/sra_sub/sub.cgi?subid=298920&from = list&action = show:submission), and included the accession number (SRP044202).

DGE tag profiling
10μg purified mRNA of each sample was used for synthesizing the first-and second-strand cDNA. A library containing all the possible fragments of CATG+17 bases in length among the reference gene sequences was generated. All clean tags were mapped to the reference sequences allowing for only 1 bp mismatch. Clean tags that mapped to reference sequences from multiple genes were filtered. The remaining clean tags were designed as unambiguous clean tags, the number of which was calculated for each gene and normalized to the number of transcripts per million clean tags (TPM).

Analysis of differential gene expression
The detection of differential gene expressions across samples was performed using a rigorous algorithm method [7] and the differences of gene expression profile between berberinetreatment group and control group were calculated and compared by DESeq software which is good for differential expression analysis in biological samples. False discovery rate (FDR) was used to determine the P value threshold in multiple tests and analyses. For our experiment, we had three repeats for each sample group, and we only considered the significant different gene expression with a P value equal or less than 0.05 and the value of |log 2 ratio | 1, which means that the fold change of gene expression between berberine-treatment group and control group is equal to or greater than 2.

Gene ontology and Kyoto Encyclopedia of Genes and Genomes
The gene ontology (GO) classification system was used to determine the possible functions of all differentially expressed genes. P-value was calculated for each changed gene. A corrected P-value 0.05 was selected as a threshold for significant enrichment of the gene sets. Pathway enrichment analysis can further identify significantly enriched metabolic pathways or signaling transduction pathways using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database. Pathways with a P-value 0.05 are considered as significantly enriched pathways in DGE. GO enrichment analysis of functional significance applies a hypergeometric test to map all DEGs to the GO database, identifying significantly enriched GO terms compared to the genomic background. Pathway analysis was mainly based on the KEGG database. Two-sided Fisher's exact tests with multiple testing and the X 2 test were used to classify the pathway categories. The FDR was used to correct P values, and only pathway categories with P 0.05 were chosen.

Real-time PCR analysis
One μg of total RNA from each sample was used to synthesize the first strand cDNA using the kit from Promega (Madison, WI, US) according to the protocol of the manufacturer.
Quantitative RT-PCR was carried out in an ABI StepOnePlus (Applied Biosystems) using SYBR Green Supermix (TaKaRa) according to manufacturer instruction. The thermal cycle conditions were 94°C for 10 min for denaturation, followed by 40 cycles of 94°C for 15 s, and 60°C for 31 s for annealing, and then an extension. The primers of examined genes are listed in Table 1. Each gene was analyzed in triplicate, and the average threshold cycle (CT) was calculated by 2−ΔΔCt method [8], the comparative Ct method we used in the manuscript is known as the 2-ΔΔCt method, where ΔΔCt = ΔCt, sample -ΔCt, reference. Here ΔCt, sample is the Ct value for any sample (berberine-treated group) normalized to the endogenous housekeeping gene and ΔCt, reference is the Ct value for the calibrator (control group) also normalized to the endogenous housekeeping gene. The expression of 18S was as an endogenous housekeeping control. The value stood for the fold difference relative to the calibrator. All data were given in terms of relative mRNA expression as mean ± SD.

Statistic analysis
All data were analyzed using a Duncan's test with a minimum n = 3. P-values equal to or less than 0.05 were considered statistically significant.

Analysis of DGE libraries
We obtained 8476945 clean reads from control group 1 (C1), 14256722 clean reads from control group 2 (C2), 7708575 clean reads from control group 3 (C3), 5669955 clean reads from

GO analysis and KEGG pathways
A total of 1890 differentially expressed genes were detected, of which, 1030 were up-regulated and 860 were down-regulated in berberine-treated samples (Fig 1A and 1B). To understand their  functions, all of the 1890 differentially expressed genes were mapped to terms in GO database and compared with the whole genome background. GO has three ontologies, namely, molecular function, cellular component, and biological process. 78 differentially expressed genes including 35 up-regulated genes and 43 down-regulated genes have a GO ID and can be categorized into a total of 23 functional groups using the WEGO software [9] (Fig 2A and 2B) (S1 Dataset). Out of the three main categories of the GO classification, ''biological process" is dominant. Interestingly, some genes, related to the rRNA processing, translation, and nucleolus, were all down-regulated in berberine-treated samples (Fig 2A). KEGG (http://www.genome.jp/kegg) ontology assignments were used to classify the functional annotations of the identified genes to further understand their biological functions (Fig 2B). Among the differentially expressed genes, 128 were mapped to 76 pathways in the KEGG database in Microsporum canis DGE class, and six terms were significantly enriched (P0.05) (S2 Dataset).

Identification and verification of differentially expressed genes
The differentially expressed tags between the two samples were identified by P 0.05. A total of 1860 significantly changed genes were found with 1030 up-regulated and 860 downregulated genes in DGE libraries (Fig 1). A total of 11 differentially expressed genes including 4 up-regulated and 7 down-regulated genes from DGE libraries were selected for real-time PCR analysis to validate the DGE data. The results showed that 11 genes were demonstrated to have a consistent change for both DGE and real-time PCR while cytochrome P450, lipase 2, zinc transporter zupT and symbiotic chitinase genes had no significant difference in real-time PCR (P >0.05) (Fig 3A and 3B).

Discussion
Efficiency of medicinal plants against fungi has been reported [2,4], but studies on their underlying mechanisms are very few. Berberine chloride, a member of alkaloids and the major bioactive component of Cortex phellodendri, has been reported to display various biological and pharmacological activities [10]. Park etc. has reported that berberine could inhibit Candida krusei with MIC <1 μg/ml [11]. To determine the anti-fungal mechanism of Phellodendron amurense against M. canis, we performed in vitro experiments with berberine chloride. To comprehensively evaluate gene profile, we performed a large-scale analysis of gene expressions in M. canis upon berberine chloride treatment using DGE approach. KEGG analysis suggested that several signaling pathways are affected in M. canis upon berberine chloride treatment, including metabolism, genetic information processing, environmental information processing, cellular processes, and organismal system. The results of qRT-PCR are consistent with DGE data, confirming that changes of gene expression upon berberine chloride treatment.
Most genes (32 genes) identified were connected to cellular metabolism, biosynthesis, and stress responses. Among these genes, glycerol-3-phosphate phosphatase (XM_002848236.1) (GPPs) and Na + /H + anti-porter (XM_002845425.1) were down-regulated and up-regulated respectively upon berberine treatment. The hydrolysis of glycerol-3-phosphate to glycerol by GPPs is important for hyperosmotic stress responses, as described in the eukaryotes Saccharomyces cerevisiae and Candida albicans [12,13]. Our data suggested that increased expression of GPPs in M. canis by berberine treatment might be one of the anti-fungal mechanisms. Berberine may modulate the expression of GPPs to weaken the responses to hyperosmotic stresses, which leads to dramatic changes in circumstances around M. canis and the following apoptosis.
The phenomenon of up-regulation of Na + /H + anti-porter was consistent with the results of Wang etc. [14] Na + /H + anti-porters play a key role in maintenance of the cytoplasmic K + /Na + ratio, through pumping Na + either into or out of organelles, mainly vacuole, to exchange for H + . The over expression of Na + /H + anti-porter genes could increase tolerance under saline conditions in buckwheat [15] and Petunia hybrida [16]. Up-regulation of Na + /H + anti-porter from our experiments suggested that M. canis might modulate the expression of Na + /H + antiporter genes to tolerate the berberine exposure.
We also identified that heavy metal tolerance protein (XM_002850584.1), a member of ATP-binding cassette (ABC) transporter pathway, is down-regulated upon berberine treatment. ABC transporters constitute one of the largest protein families in nearly all the living organisms, and are involved in some detoxification processes, which help cell recovery from various stress conditions. Our results suggested that berberine treatment limited ABC transporter-related gene expressions to increase berberine accumulation in M. canis, leading to the death of M. canis.
Another important identified KEGG pathway is the steroid biosynthesis pathway. Among the genes involved in steroid biosynthesis pathway (Fig 4), cytochrome P450 enzymes (P450s), members of CYP family, are the most important. P450s are hemethiolate proteins found in all life forms from prokaryotes to eukaryotes. P450s function in the metabolism of aliphatic, alicyclic, and aromatic molecules to regulate hydroxylation, epoxidation, dealkylation, sulfoxydation, deamination, desulphuration, dehalogenation, and noxide reduction [17]. Since CYP65 family includes enzymes can transform terpenoids, members of CYP subfamily in G. clavigera may play a role in terpene modification. For example, CYP65A hydroxylates an intermediate in the biosynthetic pathway of the sesquiterpenoid mycotoxin trichothecene in Fusarium species [18]. CYPome is involved in tolerance of host defensive chemicals as well as stress responses in G. clavigera [19]. Consistently, our results showed that cytochrome P450 was upregulated upon berberine treatment, but with no significance, the results of cytochrome P450 from DGE suggested that it might be involved in tolerance of toxic chemicals (Fig 3A and 3B).
In GO term, fungal-type vacuole membrane is important for the responses of fungi exposed to berberine chloride (Fig 5). Being the largest organelles in cells of eukaryotic microorganisms, vacuoles are known to function in intracellular ion homeostasis, storage of metabolites, intracellular protein transport, and degradation of proteins [20,21]. Moreover, vacuoles could store organic and inorganic nutrients, which could detoxify the inside toxic substances. In our results, zupT (zinc transporter), the major gene related to fungal type vacuole membrane, is upregulated upon berberine chloride treatment (Fig 5). Zinc transporters reported a role for zinc in maintaining the cellular balance between cell growth and cell death [21]. Some studies reported that a role for zinc transporter gene expression in zinc-mediated apoptosis and cytoprotection, which are vital following airway inflammation in mouse-lung or cultured-sheep pulmonary artery endothelial cells [22,23]. In addition, some reports suggested that Zn is required for Akt-dependent cell survival in many types of apoptosis [24,25]. Furthermore, Zhang's investigation showed that the induced expressions of ZnT5 and ZnT7 in rat peritoneal mesothelial cells (RPMCs) are essential for the mobilization of zinc into the Golgi under highglucose conditions, thereby protecting cells from apoptosis [26]. Our results suggested that upregulation of zinc transporter zupT might be related to the compensatory ability of M. canis upon berberine treatment.
Post-transcriptional RNA regulons consist of RNP complexes that coordinately regulate the production of functionally related proteins [27,28]. RNA regulons efficiently and economically coordinate the production of components of biochemical pathways or macromolecular protein complexes in respond to endogenous and exogenous signals. Recent evidence indicates that post-transcriptional partitioning of messenger RNA subsets by RNA-binding proteins help physically localize, temporally coordinate, and efficiently translate cell cycle-related proteins. This dynamic organization of mRNAs encoding cell cycle components contributes to the overall economy of the cell cycle and consistent with the post-transcriptional RNA regulon model of gene expression. From our results, the expression of RNP was down-regulated upon berberine treatment, which means that berberine might inhibit the growth of M. canis through inhibiting the expression of RNP, the key player in cellular post-transcriptional activity.

Conclusion
DGE can be used for analyzing variation in transcriptome of the M. canis cellular responses upon berberine treatment. Our results comprehensively showed the whole gene expression profile of M. canis upon berberine treatment for the first time. Most of the expression-changed genes were identified to cellular metabolism, biosynthesis, and stress responses. The expressions of glycerol-3-phosphate phosphatase, Na + /H + anti-porter, heavy metal tolerance protein, cytochrome P450 enzymes, zinc transporter zupT, and RNP were dramatically changed upon berberine treatment and might play key roles in anti-fungal activity. Our results provided the first step toward a global understanding of the mechanisms that are involved in the exposure of filamentous fungi to azoles. Furthermore, this work demonstrates the potential utility of gene expression profiling in anti-fungal studies.