Gene expression in human fungal pathogen Coccidioides immitis changes as arthroconidia differentiate into spherules and mature

Growth of mycelia and spherules

C. immitis RS strain directly isolated from infected mice was grown on Mycosel agar (3.6% Mycosel agar, 0.5% yeast extract, and 50 μg/ml gentamicin). The animal protocol for infecting mice was approved by the Subcommittee on Animal Studies #07-017. The plates were incubated at 30°C until the mycelia covered the surface of the agar. Arthroconidia were harvested from the plate after 6 weeks of incubation at 25°C by adding 25 ml of saline. The plate was gently scraped using cell scraper and the fluid transferred to a 50 ml tube that was then vigorously mixed for 10 sec and centrifuged at 3000 rpm for 10 min at 4°C. The supernatant containing floating mycelia was discarded. The pellet containing arthroconidia was re-suspended in 10 ml saline and the suspension was passed through a 40 μM nylon cell strainer (BD Falcon). The strained suspension was centrifuged again and the pellet used to produce mycelia and spherules.

To grow mycelia, arthroconidia were washed 2 times with glucose-yeast extract (GYE) media and 2×10⁶ spores/ml were incubated in 250 ml flat-bottom Erlenmeyer flasks (Corning) in 50 ml GYE media. Four flasks were cultured in a 30°C incubator without shaking for 5 days. To grow spherules, arthroconidia were washed 2 times in modified Converse media [12]. The spores were inoculated at 4×10⁶ arthroconidia/ml into a 250 ml baffled Erlenmeyer flask containing 50 ml of modified Converse media. Eight identical flasks were set up and grown on a shaker at 160 rpm, in 14% CO₂ at 42°C. Four flasks were harvested 2 days after inoculation and the remaining four flasks after 8 days. The spherules did not rupture and release endospores within that time in this culture system.

Inhibition of growth with nitisinone

Nitisinone, 2-(2-nitro-4-trifluoromethylbenzoyl)-cyclohexane-1, 3 dione, a potent specific inhibitor of 4-HPPD was purchased from Swedish Orphan Biovitrum, Sweden. A stock solution of 30 mg/ml was made in 0.2 M NaOH. Nitisinone was added at several concentrations to glucose yeast extract media (GYE) or modified converse media in the presence of 2×10⁶ spores/ml in a 15 ml round-bottom tissue culture tubes (BD Falcon). The culture was grown as described above for mycelial and spherule growth. The control tubes contained equal amounts of 0.2 M NaOH without Nitisinone. For microscopy, 1% formaldehyde was added to the culture overnight and the tubes were centrifuged 10,000 rpm for 10 min. The pellet was re-suspended in Lactophenol Aniline blue stain (Remel) and examined microscopically.

RNA isolation

C. immitis mycelia were harvested by straining the media from four cultures through a 40 μM nylon cell strainer (BD Falcon). The mycelia were picked up from the cell strainer using a sterile disposable loop (BD Falcon) and dropped in a 2 ml ZR BashingBead lysis tube with 0.5 mm beads (Zymoresearch) and 0.5 ml Qiazol reagent (Qiagen). The tubes were arranged in a pre-cooled Tissuelyzer II adapter (Qiagen) and mycelia was disrupted by shaking at 50 Hz for 25 min. Spherules in Converse media were harvested from four 2 day cultures and four 8 day cultures. The cell concentration was determined by counting the spherules in Lactophenol Aniline blue stain. The media was centrifuged at 10,000 rpm for 10 min at 4°C. Qiazol (Qiagen) was added to the cell pellet at 4×10⁶ spherules/ml and 0.5 ml of the mixture added to a 2 ml ZR BashingBead lysis tube with 0.5 mm beads (Zymoresearch). Total RNA was purified from mycelia and spherule samples (4 replicates/condition) using the RNeasy Microarray tissue mini-kit (Qiagen) in a Qiacube machine (Qiagen). If necessary RNA was concentrated or re-purified using RNeasy Minelute Cleanup kit (Qiagen) according to the manufacturer’s protocol. Total RNA concentrations were initially determined using a NanoDrop spectrophotometer (NanoDrop) and RNA quality was assessed using a 2100 Bioanalyzer (Agilent). All samples had a RNA integrity number greater than 7.

Microarray design and hybridization

Known and predicted ORFs from the C. immitis genome (RS strain) were previously identified using sequence data available at the Broad Institute [14]. This information was supplied to Roche Nimblegen in order to manufacture a custom oligonucleotide array consisting of 68,927 probes (Nimblegen custom array OID30589). Probes were 60 nucleotides in length and the expression of the majority of genes was assayed using 7 different probes printed in duplicate. The expression of small genes was assayed with fewer probes. Twelve custom microarrays fit on a single slide such that all the samples in this study (4 × mycelia, 4 × day 2 spherule, and 4 × day 8 spherule) could be assayed for gene expression in a single experiment to eliminate technical batch effects. Ten μg of total RNA at a concentration greater than 1 μg/ml from each sample was used for microarray hybridization. Total RNA was converted to cDNA, labeled with dye, and hybridized to the microarray by the VA San Diego Gene Chip Microarray Core according to the Nimblegen protocol. All C. immitis genes are referred to by their locus tag and further information about these genes can be found at the Coccidioides group database at the Broad Institute http://www.broadinstitute.org/annotation/genome/coccidioides_group/MultiHome.html. FungiDB (http://fungidb.org/fungidb/) was also used for annotation because it has more informative gene names for many genes.

Microarray data analysis

Quality control analysis and normalization of microarray gene expression data were performed as previously described [15]. Briefly, several quality control assessments (e.g., boxplots and volcano plots) were applied to assess microarray data quality. Unsupervised clustering was also performed using the web-based tool ANAIS [16] to determine if samples clustered as expected based on the expression of genes in each sample. All arrays passed quality control filters and no outliers were found. Differentially expressed probes were identified between mycelium, day 2 spherule and day 8 spherule conditions using a one-way ANOVA and the Tukey post hoc test implemented in GeneSpring GX version 11.5 (Agilent Technologies Inc.). The false discovery rate (FDR) associated with multiple tests was corrected for using the Benjamini-Hochberg method [17]. In a conservative approach, a gene was only identified as differentially expressed if all probes for that gene had a fold change greater than 2 or less than −2 and an ANOVA p-value (Tukey and FDR corrected) less than 0.05. Fold changes were calculated for each gene that passed this filter by averaging across the seven probes. The overlap in up and downregulated genes between day 2 and day 8 spherules compared to mycelia was visualized in Venn diagrams using BxToolBox (http://bioinforx.com/). The ANOVA analysis was also used to identify genes that were differentially expressed between day 2 and day 8 spherules where positive fold changes are indicative of greater expression at day 8 compared to day 2 and negative fold changes suggest decreased expression. Gene expression data are available at the Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/) under accession number GSE44225.

PFAM and GO analysis

PFAM enrichment was determined using a tool at the Broad Institute http://www.broadinstitute.org/annotation/genome/coccidioides_group/BatchSelect.html?target=GeneEnrichment.html. This tool looks for over-representation of PFAMs in up- or downregulated genes using a hypergeometric test, and only PFAMs with an FDR-corrected p-value <0.05 were considered significant.

GO terms were assigned to C. immitis genes by reciprocal homology searches at the protein level against the Saccharomyces cerevisiae proteome using BLAST (Additional file 1: Table S1). UniProt IDs were obtained using the C. posadasii homologs of C. immitis genes because many more C. posadasii genes have UniProt IDs. The Biological Networks Gene Ontology (BiNGO) plugin (version 2.441) [18] for Cytoscape (version 2.8.3) was used to identify those GO terms related to biological processes that were over-represented for differentially expressed genes identified between each of the three comparison groups (mycelia vs. day 2 spherule, mycelia vs. day 8 spherule, day 8 vs. day 2 spherule). BiNGO preserves the hierarchical relationship between GO terms. Significance was assessed with a hypergeometric test and only GO terms with an FDR-corrected p-value <0.05 were considered significant.

Gene annotation

C. immitis protein kinases were identified and classified by orthology with the curated Trichophyton rubrum kinome [19]. Non-orthologous kinases were identified and classified by searching the proteome with a protein kinase HMM built from an alignment of Dictyostelium kinases [20] followed by a BLAST against the curated kinase database (http://kinase.com/) [21]. Kinase abbreviations are provided in Additional file 2: Table S3. Signal peptides in the proteins coded for in the C. immitis genome were identified using artificial neural networks implemented in SignalP version 4.0 [22].

RT-qPCR confirmation of gene expression

Microarray gene expression was confirmed by RT-qPCR for 24 genes. Three highly expressed genes with low standard deviation across the 12 samples were selected as normalizers (CIMG_01599, CIMG_10083 and CIMG_12902). SYBR® Green primers were designed using Primer Express version 3.0 (Applied Biosystems Inc.) and obtained from Integrated DNA Technologies, Inc. (Coralville, IA). Reverse primers were designed to span a splice site in the same region of the gene probed by the microarray. Duplicate RT-qPCR reactions were performed for each sample using 50 ng of reverse transcribed RNA per reaction. Fold changes were calculated as described previously using the 2^-∆∆CT method [23] implemented in the DataAssist software version 3.0 (ABI), and significance was determined using one-way ANOVA in the R statistical package (version 2.13.2).

Genes differentially expressed in mycelia and spherules

Gene expression was assessed in a total of 12 samples derived from 4 replicate samples isolated from the following three growth phases: mycelia, day 2 spherules, and day 8 spherules. A photograph of mycelia and day 2 and day 8 spherules grown in Converse medium is shown in Figure 1. The image shows the difference in shape and size between spherules and mycelia and the increase in spherule size between 2 and 8 days of culture. A custom oligonucleotide microarray (Nimblegen), which contained probes for all predicted ORFs of the RS strain of C. immitis was used to assess gene expression. 91% of the predicted ORFs were expressed in either mycelia or spherules, suggesting that the annotation and the detection of hybridization were robust. Unsupervised clustering using the expression of all genes on the microarray revealed that mycelia samples clustered distinctly from spherule samples. Furthermore, spherule samples formed two sub-clusters based on the number of days in culture. A dendrogram showing that the four replicate samples cluster together is shown in Additional file 3: Figure S1. Fungal morphologic stage was the dominant determinant of the pattern of gene expression.

Genes that were significantly differentially expressed (p < 0.05) between the three conditions (mycelia, spherules on day 2 and 8) were identified in a supervised approach using a one-way ANOVA with appropriate corrections for multiple testing (see Methods). All the up- and downregulated genes differentially expressed between each of the three conditions are detailed in Additional file 4: Table S2. A heatmap depicting expression levels in each sample for the top 100 differentially expressed genes is presented in Figure 2. The heatmap indicates there was limited variation in gene expression across the four replicates within each of the three conditions, suggesting that the data was highly reproducible. Multiple patterns of gene expression are evident comparing the three different conditions we studied. One cluster of genes was expressed to a lesser extent in the mycelia condition and a greater extent in both spherule conditions and another cluster of genes were expressed at a higher level in mycelia than in spherules. The expression of four genes (CIMG_08103, CIMG_09765, CIMG_10037, CIMG_10264) exhibiting the upregulated pattern was confirmed by RT-qPCR (see Figure 3 below).

There were a total of 2208 genes (22% of the genome) that were differentially expressed between spherules at one or both the time points we studied and mycelia. Figure 4 shows Venn diagrams depicting up- and downregulated genes in day 2 and day 8 spherules compared to mycelia. About a third of the differentially expressed genes were up- or downregulated in both day 2 and day 8 spherules compared to mycelia. However, similar numbers of genes were exclusively upregulated in either day 2 (N = 443) or day 8 (N = 319) spherules, or exclusively downregulated at either day 2 (N = 565) or day 8 (N = 233) spherules. The difference in gene expression between day 2 and day 8 spherules was apparent when we compared day 2 and 8 spherules directly to each other; 1,197 differentially expressed genes (12% of the total genome) were identified (Additional file 4: Table S2). Therefore, although gene expression by environmental form of the fungus and the parasitic form were quite distinct as might be expected, gene expression by young and mature spherules was also quite different from each other. Not only were there differences in which genes were expressed at each stage, but also the degree of modulation was large. For example, the maximum difference in expression of a gene (CIMG_10264) between day 2 spherules and mycelia was 48.6 fold and the median modulation between mycelia and day 2 spherules was 3.26.

A recent paper by Whiston et al. assessed transcription in C. immitis and C. posadasii mycelia and day 4 spherules by RNA-seq [13]. We have compared our results to theirs. The two studies used different methods for assessing changes in gene expression. We used microarray technology to estimate transcript abundance while Whiston et al. used RNA-seq to estimate transcript abundance [13]. The literature suggests that these methods should yield comparable results [24]. Despite this difference in methodology, we confirmed the upregulation of 25% of the genes that Whiston found to be upregulated in spherules. Conversely, 43% of genes that we have found to be upregulated in day 2 and day 8 spherules were also upregulated in day 4 spherules in the Whiston study (Additional file 5: Figure S2). Despite the differences in the two studies many of our conclusions are similar (see below).

We know from previous experiments that some genes are overexpressed in spherules compared to mycelia. Some of these genes, such as the spherule outer wall glycoprotein (CIMG_04613) [25] and the parasitic-phase specific protein PSP-1 (CIMG_05758) [26] were up regulated more than four fold in spherules in this experiment (Additional file 4: Table S2). Other genes, such as the metalloproteinase Mep1 (CIMG_06703), which has been found to be expressed at high levels in endosporulating spherules in C. posadasii was not found to be over-expressed in this experiment [27]. We also examined the expression level of the Mep1 gene by RT-qPCR and found that its expression was slightly downregulated in spherules compared to mycelia, rather than upregulated as previously reported (see below). Whiston et al. also examined the expression of this gene and found that it was upregulated in C. posadasii spherules but not C. immitis spherules [13].

Confirmation of differential expression by RT-qPCR

Twenty-four differentially expressed genes as detected by microarray analysis were selected for confirmation by RT-qPCR (Figure 3). Genes were selected for RT-qPCR confirmation of gene expression based on the magnitude of fold change (up- or downregulation) between mycelia and day 2 spherules, mycelia and day 8 spherules, and day 2 and day 8 spherules, and their identification in the PFAM or GO analysis. The significant differential expression (p < 0.05, t-test) of each of these 24 genes was confirmed for at least one of the three comparison groups. In the majority of cases, RT-qPCR analysis detected larger differences in expression than microarray analysis. This is probably because RT-qPCR has a greater dynamic range than microarray does [28].

PFAM analysis

Gene ontology analysis

Upregulated or downregulated genes in day 2, day 4 and day 8 spherules

We identified 153 genes that were upregulated more than two fold in day 2 spherules, day 8 spherules and day 4 spherules in the Whiston study [13]. 140 genes were downregulated more than two fold in all three spherule samples. 57% of the upregulated genes and 42% of the downregulated genes had no function annotation (Additional file 7: Table S4). Many of these unannotated genes were highly differentially expressed suggesting that they may be important for spherule development. One upregulated gene that has not been discussed is opsin-1 (CIMG_08103, maximal upregulation 25.72). This gene is closely related to the bacterial rhodopsin gene coding for G protein coupled receptors; its function in fungi has not been determined [59]. Another gene that was upregulated at all three spherule time points was the sulfite transporter Ssu1 (CIMG_05899, maximum upregulation 6.37). Downregulated genes of interest include several glucosidases, glucanases and a chitosanase. Two septin genes are downregulated in spherules. Septin genes are important regulators of the cytoskeleton and play a role in determining cell shape [60]. Why these genes are downregulated is unclear since the spherule is undergoing extensive cellular remodeling. Perhaps septins are required for polar growth and other regulators are needed for isotropic spherule growth. Further analysis of the relatively small group of genes that are consistently up- or downregulated throughout spherule development may be useful for understanding the pathogenic phase of this organism.

Comparison of results to those obtained in other pathogenic fungi

The dimorphic pathogenic fungi are phylogenetically closely related [61] so it is reasonable to ask if genes important for conidium to yeast transformation in those pathogens are also important for arthroconidia to spherule transformation in Coccidioides. One H. capsulatum gene that is required for mycelium to yeast transformation is the alpha amylase (AMY1) gene. This gene affects the amount of α (1,3) glucan in the yeast cell wall, and that carbohydrate appears to be important for the ability of the yeast to survive in macrophages and in mice [62]. This gene has a nearly identical homolog in C. immitis, CIMG_03142, that was upregulated 3.6 fold in day 2 spherules and 3.39 fold in day 8 spherules. Whiston et al. also found it to be upregulated in spherules [13]. Another H. capsulatum gene that is required for yeast formation is α glucan synthase (AGS1) gene [62]. This enzyme catalyzes the production of α (1,3) glucan in the cell wall that obscures the β (1,3) glucan and prevents activation of innate immunity via the dectin-1 receptor [62]. C. immitis has an AGS1 gene (CIMG_13256) that was upregulated in the day 8 spherule (2.48 fold) but not day 2 spherules. Whiston et al. found this gene to be upregulated 1.93 fold in spherules compared to mycelia [13]. There is no literature describing the relative amounts of α (1,3) glucan and β (1,3) glucan in C. immitis mycelia or spherules. We know, however, that there is enough exposed β (1,3) glucan in Coccidioides spherules to stimulate macrophages to produce cytokines via dectin-1 [63].

Two genes coding for transcription factors, Ryp2 and Ryp3, have been found to be essential for conversion from filaments to yeast in H. capsulatum[64]. These genes are overexpressed in the yeast phase of H. capsulatum[64]. C. immitis has nearly identical homologs of these genes but they were not overexpressed in either day 2 or day 8 spherules, suggesting that they may not be required for the transformation from mycelium to spherule.

Gene disruption experiments in B. dermatitidis have shown that a histidine kinase, DRK1, is required for the transformation from filaments to yeast [65]. It is not clear from the literature whether or not this gene is overexpressed in the B. dermatitidis yeast phase. C. immitis has a very closely related homolog of this gene (CIMG_04512) but it was not up or down regulated in day 2 or day 8 spherules. In another dimorphic pathogenic fungus, S. schenckii, the calcium/calmodulin kinase I gene (SSMK1) was found to be required for formation of yeast [53]. There are two genes in C. immitis that are highly homologous to the S. schenckii SSMK1 gene; neither one of these was up- or downregulated in day 2 or day 8 spherules.

A number of studies have been done studying the transcriptome of P. brasiliensis[66, 67]. One study identified the 4-HPPD gene to be required for P. brasiliensis conidia to convert to yeast [66]. They found that the 4-HPPD gene expression was upregulated in the yeast form and that a biochemical inhibitor of this enzyme, nitisinone, inhibited mycelium conversion to yeast. 4-HPPD (E.C. 1.13.1127) is an enzyme that converts 4-hydroxyphenylpyruvate to homogentisate that is involved in the synthesis of tyrosine, phenylalanine, and ubiqinone (KEGG, whttp://www.genome.jp/keg).

There are two homologs of the 4-HPPD in the C. immitis genome, which have significantly different sequences. One of these, CIMG_01466 was upregulated 5.27 fold in day 2 spherules and 3.80 fold in day 8 spherules. This gene was also found to be upregulated in spherules by Whiston et al. [13]. The other homolog, CIMG_01310, was downregulated −23.67 fold in day 2 spherules and −6.09 fold in day 8 spherules. The biggest difference in sequence is that CIMG_01466 has two substantial deletions compared to CIMG_01310. These deletions flank the highly conserved site that is predicted to contact the active site metal ion [68]. Furthermore, CIMG_01466 had substitutions in the predicted metal ion contact site, suggesting that it may not be an active enzyme. Nevertheless, we tested the effect of nitisinone on mycelial growth and mycelium to spherule conversion. We found that nitisinone inhibits mycelial growth at concentrations as low as 1 μg/ml (Figure 5). Surprisingly, there was no effect on mycelium to spherule conversion (data not shown). This is distinctly different from the results seem in P. brasiliensis. Our data suggests that 4-HPPD enzyme activity is not required for mycelium to spherule conversion or the growth of spherules but it is important for mycelial growth.