Functions of thga1 Gene in Trichoderma harzianum Based on Transcriptome Analysis

Trichoderma spp. are important biocontrol filamentous fungi, which are widely used for their adaptability, broad antimicrobial spectrum, and various antagonistic mechanisms. In our previous studies, we cloned thga1 gene encoding GαI protein from Trichoderma harzianum Th-33. Its knockout mutant showed that the growth rate, conidial yield, cAMP level, antagonistic action, and hydrophobicity decreased. Therefore, Illumina RNA-seq technology (RNA-seq) was used to determine transcriptomic differences between the wild-type strain and thga1 mutant. A total of 888 genes were identified as differentially expressed genes (DEGs), including 427 upregulated and 461 downregulated genes. All DEGs were assigned to KEGG pathway databases, and 318 genes were annotated in 184 individual pathways. KEGG analysis revealed that these unigenes were significantly enriched in metabolism and degradation pathways. GO analysis suggested that the majority of DEGs were associated with catalytic activities and metabolism processes that encode carbohydrate-active enzymes, secondary metabolites, secreted proteins, or transcription factors. According to the functional annotation of these DEGs by KOG, the most abundant group was “secondary metabolite biosynthesis, transport, and catabolism.” Further studies for functional characterization of candidate genes and pathways reported in this paper are necessary to further define the G protein signaling system in T. harzianum.


Introduction
Trichoderma is an important fungal genus, whose species exhibit favorable properties, such as diverse mechanisms of antagonistic action, a broad spectrum of activity in plant disease prevention and control, survival under unfavorable conditions, and environmental friendliness. Trichoderma harzianum is one of the most commercially used biofungicides, particularly T. harzianum T22 and T39 [1]. Growth, conidiation, secondary metabolism, and mycoparasitism are all important processes that contribute to biofungicidal property [2]. An enhanced understanding of the signal regulatory mechanism in T. harzianum is necessary to further explore the fungi's extraordinary biocontrol potential.
G protein-mediated signal transduction system is an important transmembrane signaling system in eukaryotic cells. It plays a key role in the regulation of cellular reactions and the transmission of extracellular signals to cells. The role of G protein in the transmission of external stimuli has been studied in detail in genetic fungi models, such as Neurospora, Penicillium [3], and Aspergillus [4]. Heterotrimeric G protein is composed of , , and subunits; each subunit is encoded by independent genes. Among the G protein subunits, the G subunits were found more frequently and reported to regulate vegetative growth, conidiation, and the mycoparasitic responses in fungi [5]. However, current observations showed that a particular G subunit may have different functions in different fungal species [6]. G subunits are classified into three groups according to conserved motif sequences [7], and several subgroup G I subunits from T. atroviride and T. virens have been studied [4,7]. G from T. atroviride, tga1, was reported to negatively regulate conidiation but had no effect on hyphal growth [8]. TgaA is a homology gene of tga1 from T. atroviride, which does not influence growth and conidiation in T. virens. The ΔtgaA mutants grow normally and sporulate like the wild type but possess a reduced ability to colonize Sclerotinia sclerotiorum, whereas they are fully pathogenic against Rhizoctonia solani [4].

BioMed Research International
In our previous study, a subgroup I G gene (thga1) was cloned from the T. harzianum Th-33 genome. Results showed that THGA1 has the same amino acid sequence as that of TGA from T. atroviride but with different functions. Compared with the wild-type Th-33, the Δthga1 mutants changed significantly in biological characteristics and physicochemical properties. In particular, the hyphal growth rate dropped by about 40%, conidiophore branches became sparse, secondary branch and phialide numbers reduced, conidiation was delayed for about 20 h, and conidia yield declined by about 300-fold. In addition, the hydrophobicity of the mutants weakened, intracellular cAMP levels decreased by about 50%, and the inhibition of the mutant against plant pathogen of R. solani reduced significantly. The results showed that the thga1 gene positively affected the growth, conidial production, hydrophobicity, and antagonism of T. harzianum to R. solani.
The present study characterized the genes under-or overexpression in the mutant compared to WT associated with the G subunit, thga1, using next-generation RNA sequencing (RNA-seq) technology, to gain insight into the regulatory mechanism of G subunits thga1 in Th-33. This study is the first initiative to use RNA-seq for identifying differentially expressed genes (DEGs) to clarify the function of G in Trichoderma by genome-wide transcriptional analysis of hyphal cells of the wild-type Th-33 and its Δthga1 mutants. The results may reveal that G regulated a series of biological processes as well as new targets of thga1 function.

Strains and Culture
Media. The wild-type T. harzianum Th-33 was isolated from soil samples in the Beijing region as described previously. The Δthga1 knockout mutant was created with hygromycin B resistance by homologous recombination and then purified by isolation of single conidia [9]. Mutant colony was identified by southern hybridization. The fungi were grown on potato dextrose agar (PDA) at 28 ∘ C for conidia production. The PDA medium was covered with cellophane for mycelia collection. Conidial suspensions were prepared by adding sterilized distilled water to the PDA plates. The conidial suspensions were inoculated on the PDA medium covered with cellophane and cultured for 32 h prior to conidia formation at the tip of the mycelia. The mycelia were then scraped off the cellophane, washed with cold distilled water, frozen in liquid nitrogen, and ground to a fine powder. Equal mixture of three biological samples was sequenced.

Preparation of the cDNA Library and Illumina
Sequencing for Transcriptome Analysis. Total RNA was extracted using Trizol reagent (Invitrogen, CA, USA) according to the manufacturer's instruction. The RNA quality and quantity were determined using an Agilent 2100 Bioanalyzer. The Qubit RNA Assay Kit was used for accurate quantification of the initial total RNA. After total RNA extraction and DNase I treatment, magnetic beads with oligo (dT) were used to isolate mRNA. The mRNA was mixed with the fragmentation buffer to promote fragmentation into short segments. Subsequently, cDNA was synthesized using SuperScript II reverse transcriptase following the manufacturer's protocol. The short fragments were connected with adapters. The suitable fragments were selected as templates for PCR amplification. During the QC steps, a 2100 Bioanalyzer High Sensitivity DNA chip was used for quantification and qualification of the sample library. RNA libraries were sequenced by paired-end mode using an Illumina HiSeq2000 system.

Sequence Alignment.
Clean reads were obtained after removal of low quality reads, reads with adaptors, and reads with unknown nucleotides larger than 5%. All reads were deposited in the National Center for Biotechnology Information (NCBI) database and can be found under accession number SRS823675. The reads were mapped to the reference genome [10] of T. harzianum Th-33 using TopHat (Version: 2.0.11) program [11].

Expression Analysis.
Expression values were obtained by calculating the fragments per kilobase of transcript per million mapped reads, and differential gene expression was analyzed using the Cufflinks program [12]. Genes differentially expressed with more than twofold changes and at values less than 0.05 were identified as DEGs. The threshold of the value in multiple tests was determined by the value for the false discovery rate (FDR) [13]. We used "FDR ≤ 0.001 and the absolute value of log 2 fold change (log 2 FC) ≥ 2" as the threshold to assess the significance of gene expression differences.

Functional Annotation of T. harzianum Transcriptome
and Classification of DEGs. The T. harzianum transcriptomes were compared against amino acid sequences available at the UniProt database using BLASTx algorithm. For each query sequence, the associated hits were searched for their respective gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) results. The highest bit score was selected, with -value threshold of 1 201. Protein families were classified by searching the assembled transcripts against Pfam and InterProScan. Potential transmembrane domains of the G protein-coupled receptors (GPCRs) selected from InterProScan and Pfam analysis were predicted using transmembrane domain prediction tool TMHMM v2.0. GPCRs predicted to contain the seven-transmembrane domain were used for further analysis. We used the CAZy database to annotate the carbohydrate-active enzymes (http://www.cazy .org/). We annotated transcription factors (TFs) of DEGs based on protein sequence homology using the Fungal Transcription Factor Database (http://ftfd.snu.ac.kr/). The amino acid sequences of DEGs were further analyzed to predict secreted proteins. Sequences smaller than 70 amino acids were not considered for further analysis. The remaining sequences with positive SignalP prediction for signal peptide cleavage site at the N-terminal region between 10 and 40 amino acids, without any transmembrane region, were selected as the candidate secreted proteins. Secreted proteins with lengths less than 200 amino acids and cysteine content more than 4% of the protein were identified as small molecule cysteine-rich secreted proteins (SSCPs).

Validation of RNA-Seq Results by Quantitative Real-
To validate the expression profiles of the assembled genes obtained through sequencing data analysis, qRT-PCR was performed for selected genes. Genes were randomly selected, and four reference genes were used, namely, beta-tubulin 1, actin, transcription elongation factor, and ubiquitin-conjugating enzyme (UCE). Primers used in qRT-PCR were designed using Primer Express Software v2.0 (see S2 Table in    Among the 20 significantly upregulated genes (S1 Table), 13 genes exhibited defined functions, including six metabolismrelated genes and two genes predicted to encode enzymes. Among the significantly downregulated genes (S1 Table), seven were metabolism-related genes and four were catalytic activity genes. The changes in expression were evidently connected to metabolism. Hence, thga1 might regulate a series of biological processes through metabolism.

Real-Time qRT-PCR Analyses.
To verify the quality of the assembly, cDNA fragments of 15 randomly selected unigenes were amplified using unigene-specific primers (S2 Table) and then sequenced. Four reference genes (betatubulin 1, actin, transcription elongation factor, and UCE) were used as endogenous control. Among the reference genes, UCE remained constant in all treatments and showed the best performance in qRT-PCR analysis using Best Keeper program. Expression patterns of the tested genes are shown in Figure 1. Three genes with infinite fold changes in the RNAseq results were not shown in the figure, as they cannot be shown in the bar chart. Expression patterns determined by real-time qRT-PCR were consistent with those obtained by RNA-seq, thereby confirming the accuracy of the RNA-seq results reported in this study.

GO Enrichment Analysis.
Blast2GO [15] program was used for GO analysis. The program extracted the GO terms associated with the homologies identified by BLAST and returned a list of GO annotations, which were presented as hierarchical categories of increasing specificity. GO enrichment analyses were performed using Fisher's exact test with multiple testing corrections and an FDR of 0.05. A total of 517 DEGs were categorized into 707 functional groups in three main categories, namely, "cellular component," "molecular function," and "biological process" (S4 Table and Figure 3). Some unigenes were assigned to multiple categories of GO terms, whereas others could not be assigned to a given GO term. In the biological process category, "metabolic process" (381, 73.69%), "cellular process" (201, 38.87%), "singleorganism process" (115, 22.24%), "localization" (100, 19.34%), and "establishment of localization" (99, 19.15%) were the most abundant terms. In the molecular function category, genes associated with "catalytic activity" (351, 67.89%), "binding" (191, 36.94%), and "transporter activity" (54, 10.44%) were the most abundant. In the cell component category, "membrane" (124, 23.98%) and "membrane part" (100, 19.34%) were the most abundant terms. These findings indicated that the main changes in expression between the wild-type Th-33 and the thga1 mutant were those related to membrane parts, metabolic processes, and catalytic activities.

Genes Related to G Protein Signal
Pathway. G protein transduced signals received by the heptahelical GPCRs from outside the cell influence numerous regulatory pathways via their respective effectors, which in turn affect the activities of secondary messengers [16,17]. Three G -coding genes, one beta, and one gamma subunit genes were found in the transcriptomes and the genome of Th-33, which corresponded well with the data of T. virens [18], Neurospora crassa [19], and many other filamentous fungi [6]. Thga1, the subgroup I G gene, was not expressed in the Δthga1 mutant. Another G gene, a subgroup II tgaB homologous gene in T. virens [4], exhibited downregulated expression (twofold), whereas the third G gene and G gene were not differentially expressed. These results showed that knockout of thga1 caused the downregulated expression of another G gene.
The signals transmitted through a heterotrimeric G protein signaling cascade originate from the activation of plasma membrane-localized GPCRs. The identification and characterization of GPCRs will provide insights into how Trichoderma communicates with G protein and downstream genes. Thirty-nine GPCRs or putative GPCR genes were found in the transcriptome of Th-33. Six of them were differentially expressed in the mutant (five of them belong to the PTH11-type; Table 2). PTH11-type GPCRs were reported to influence light responsiveness, glycoside hydrolase gene transcription, sexual development [20,21], surface recognition, and pathogenicity [22]. PTH11 GPCR genes were reported to be upregulated in the mycoparasite Coniothyrium minitans during colonization of S. sclerotiorum [23]. In the current study, five PTH11-type GPCR-encoding genes were downregulated in the mutant. This observation was consistent with the phenomenon that the inhibitory effect of the mutant against R. solani decreased significantly compared with that of the wild-type Th-33. cAMP acts as a secondary messenger for morphogenic signals in both prokaryotes and eukaryotes, and it is generated by adenylate cyclase and degraded by phosphodiesterase [24]. In T. virens, the intracellular cAMP levels can be regulated by an adenylate cyclase gene, tac1 [25], and it has been reported to influence conidiation [26]. In Aspergillus nidulans, a G -subunit GanB mediates a rapid and transient activation of cAMP synthesis in response to glucose during the early period of germination [27]. In our previous study, intracellular cAMP levels in the Δthga1 mutant decreased by about 50% of that in the wild-type Th-33. According to the transcriptomes of the Δthga1 mutant and the wild-type Th-33, the adenylate cyclase-encoding gene was not expressed differentially, whereas the cAMP phosphodiesterase gene was found to be upregulated in the mutant. These findings indicated that the decrease in intracellular cAMP levels was caused by the increase in cAMP phosphodiesterase activity and not by the decrease in adenylate cyclase activity. This study showed that the G protein subunit in Trichoderma may be involved in the mechanisms regulating intracellular cAMP levels.

Carbohydrate Activity Enzymes (CAZymes).
CAZymes in fungi are correlated with their nutritional modes and infection mechanisms [28]. Various CAZyme genes exist in the Th-33 genome, of which 66 differentially expressed CAZyme genes were identified in the Th-33 transcriptome, including 30 genes from the glycoside hydrolase family (GH), 14 from the auxiliary activity family, 13 from the carbohydrate esterase family, 7 from the glycosyltransferase family, and 2 from the carbohydrate-binding module family (S6 Table). In addition, the number of downregulated CAZyme genes (43 genes) in the mutant was greater than that of upregulated genes (23 genes) ( Figure 5). The GH family has been described to play a central role in mycoparasitism in Trichoderma. For example, GH18 family is highly expanded in T. atroviride and T. virens during mycoparasitism [29]. GH16 was reported to be upregulated during the mycoparasitic reaction of T. atroviride against R. solani [30]. In the present study, the inhibition of the Δthga1 mutant T. harzianum against R. solani and Phytophthora capsici was reduced significantly, and the mycelia of Δthga1 could not cover and parasitize the pathogenic fungi's mycelia in confrontation tests. The heterotrimeric G protein pathway is crucial in the interconnections of nutrient signaling in Trichoderma [21]. Results implied that the knockout of thga1 caused the decrease in the activities of the CAZymes, which consequently influenced nutrient condition, further retarded the growth rate, and decreased the mycoparasitic ability of Trichoderma.
3.9. TFs. Twenty-nine TFs were identified in the Δthga1 mutant and compared with that of the wild-type Th-33. Of these TFs, 15 were upregulated and 14 were downregulated TFs (S7 Table). Among those differentially expressed TFs, the Zn 2 Cys 6 family (53%) was the most dominant TF family, followed by the C 2 H 2 zinc finger (31%) TF family ( Figure 6). The functions of most TFs in fungi appear to be highly diverse and remain rather unclear. Of the known TF families, the Zn 2 Cys 6 TF family comprises TFs that are fungal-specific and dominant [31]. These TFs are involved in primary and secondary metabolism, drug resistance, and meiotic development (http://ftfd.snu.ac.kr/). BglR, a Zn 2 Cys 6 -type TF in T. reesei, can increase the expression of the -glucosidase gene [32]. C 2 H 2 zinc finger in fungi is involved in pathogenicity, carbon catabolite repression, acetamide regulation, and differentiation of ascogenous or fruiting body. GipA in A. fumigatus was reported to be involved in secondary  metabolic processes [33]. In A. nidulans, brlA defined the central regulatory pathway that controls conidiation-specific gene expression, but no orthologue of brlA in Trichoderma conidiophores currently exists [34]. TFs can control gene transcription rates and are key mediators of cellular function [35]. In the present report, differentially expressed TFs might play key roles in cell processes in Trichoderma under the control of the G protein signal pathway. Understanding the regulatory mechanisms of these TFs and their functions could provide valuable insight into the signal control pathways of cell processes in Trichoderma.
3.10. Secretome of T. harzianum. A total of 137 genes were annotated as secreted proteins, including 51 upregulated and 86 downregulated genes among the DEGs. Secreted proteins play an important role in the process of cell signaling, cell proliferation, differentiation, apoptosis regulation, development, and other important biological processes. In Botrytis cinerea, bcg1 encodes subunits of heterotrimeric GTPbinding proteins; the bcg1 null mutants differ in colony morphology from the wild-type strain, which do not secrete extracellular proteases [36]. SSCPs comprise one of the largest groups of proteins secreted by Trichoderma [29]. Hydrophobins, which are probably the most widely known SSCPs, are found on the outer surfaces of cell walls of hyphae and conidia. Trichoderma genomes encode an unusually large number of hydrophobins, possibly for the purpose of providing flexibility in surface properties needed for conidiation, mycoparasitism, and interaction with plant roots [2,37]. In the present report, two SSCPs were found among the DEGs, and these genes were downregulated in the mutant strain, which might account for the phenomenon of reduced hydrophobicity in the mutant. This was in accordance with the report that G alpha subunit can mediate fungal conidiation and hydrophobin synthesis in Cryphonectria parasitica [38]. Furthermore, SSCPs were shown to be highly expressed in T. atroviride during colonization of R. solani, which suggested a potential role in mycoparasitism [30]. These genes may also explain the significant reduction in the inhibition of the mutant against the plant pathogen R. solani. 3.11. Secondary Metabolism. KOG functional classification of DEGs revealed that the number of genes involved in secondary metabolite biosynthesis, transport, and catabolism was 71, in which cytochrome P450s (CYPs) accounted for nearly 24% (S8 Table). CYPs play an important role in the physiology of fungi and are involved in the biosynthesis of secondary metabolites (SMs) and in detoxification [39]. Fungi produce a wide range of SMs that are not directly essential for growth and yet have important roles in signaling, asexual conidiation [40], and interactions with other organisms [41][42][43]. Trichoderma spp. are a rich source of SMs, such as nonribosomal peptides, polyketides, terpenoids, and pyrones. Although a clear correlation exists between conidiation of the fungus and the secretion of antifungal metabolites in the T. atroviride parental strain, its type G III gene Δtga3 mutants were fully impaired in the production of peptaibols despite exhibiting a hypersporulating phenotype [5]. The results in the current report implied that thga1 regulated certain secondary metabolism processes, and analysis of SMs should be carried out for further identification.

Conclusion
Overall, this study identified transcripts that were possibly involved in several important molecular and cellular functions of T. harzianum. However, additional studies aimed at the functional characterization of the genes reported herein will aid in further defining the pathways mediated by G protein in T. harzianum. An enhanced understanding of the expression profiles of these genes could improve T. harzianum performance, either by predicting the regulation of the genes involved in sporulation or by improving their use in biotechnology processes.