RNA-Seq data for C. albicans grown in a number of conditions, stresses and environments were downloaded from the Gene Expression Omnibus (GEO) database. These data include C. albicans grown in YPD and YPS (GSE41749, Grumaz et al., 2013), in planktonic and biofilm states (GSE45141, Desai et al., 2013), with a range of weak organic acids (GSE49310, Cottier et al., 2015), and in contact with host cells (GSE56091, Liu et al., 2015). In all cases all replicates, treatments and control data for strain SC5314 were used. In total, 130 RNA-Seq samples were used, as maximizing the number of samples used allows us to create the most robust NeXO possible. We note that the majority of samples were taken from an infection-related environment composed of C. albicans cells exposed to endothelial and epithelial cells with RNA taken at 1.5, 5, and 8 h during infection (Liu et al., 2015). Therefore, the NeXO built in this study will predominantly represent C. albicans infection processes related to these cell types and timepoints and is suitable for the application in this study. The addition of more expression data, from varied environments or infection stages, would allow the NeXO to represent a broader range of functions and processes.

Raw data for all publicly available data sets was re-analyzed as experimental setups and conditions may vary between sources. The C. albicans (strain SC5314, assembly 22, haplotype A) reference genome was downloaded from the Candida Genome Database (CGD, Inglis et al., 2012). RNA-Seq reads were aligned to the reference genome using bowtie (Langmead et al., 2009) and bowtie2 (Langmead and Salzberg, 2012) for single-end and paired-end reads, respectively. Alignment was performed using default parameters on a single processor and alignment files were sorted and converted to binary format using samtools (Li et al., 2009). Reads aligned to known genes were counted using the union model with HT-Seq (Anders et al., 2014) and these counts were used to estimate gene expression as reads per kilobase per million mapped reads (RPKM) with edgeR (Robinson et al., 2010). By generating estimates of expression from raw data we ensure that the data is as comparable as possible. Finally, correlations were calculated for all pairs of genes across all samples using the Pearson's correlation coefficient (PCC).

CliXO was used to create a C. albicans NeXO from correlations in gene expression (Kramer et al., 2014). CliXO requires a user-defined noise parameter α and a parameter to infer missing edges β. CliXO has been benchmarked on several 'omics data sets for the yeast Saccharomyces cerevisiae, including gene expression profile correlations, and it has been shown that α = 0.01 and β = 0.5 produces ontologies that are very similar to that of the Gene Ontology (GO) (Kramer et al., 2014). Therefore, in this study α and β were set to 0.01 and 0.5, respectively. Varying these parameters produces similar NeXO structures. Altering the α parameter to α = 0.02 produces a NeXO structure where entities have an average similarity to the reference C. albicans NeXO of 0.66 ± 0.28 (max 1.0) using the alignment and similarity metrics of Kramer et al. (2014) with a the strict hierarchy model. Varying the β parameter to β = 0.4 and β = 0.6 also produced a similar NeXO structure with average entity similarity of 0.83 ± 0.22 and 0.75 ± 0.23, respectively.

The Gene Ontology (GO) structure (OBO file) and annotations for C. albicans were downloaded from geneontology.org. To prepare the GO for alignment to the C. albicans NeXO, the structure of the GO was extracted using “is_a” GO relations, terms with no annotations to known C. albicans genes were removed, and redundant terms were removed. Redundant terms were defined as terms that share the same gene content as their direct descendants. When redundant terms were identified the more general term, i.e., the parent, was removed. After processing of the GO there were 776, 1,919, and 3,804 terms in the Cellular Component, Molecular Function, and Biological Process ontologies, respectively. The method of Kramer et al. (2014) was used to align the C. albicans NeXO to the GO ontologies using the strict hierarchy model and to generate an alignment for all nodes.

The C. albicans NeXO is split into a number of entities with hierarchical relationships between them. As these entities will represent biological pathways and processes it was important to estimate the robustness of the entities with regard to the underlying experimental data used to construct the NeXO. To score the robustness of each term a measure of entity quality that considers the network support of a term and its robustness to perturbations of input data was taken from Dutkowski et al. (2013). In this measure, network support NS(e) for an entity e is defined as the enrichment for co-expression between genes within an entity where co-expression is defined as PCC > 0.2. Enrichment [-log(P-value)] is estimated using the hypergeometric distribution where the sample size is equal to all possible co-expressed genes in an entity, sample successes are the observed co-expressed genes in an entity, the population size is all possible co-expressed genes in the NeXO, and population successes are the observed co-expressed genes in the NeXO.

To quantify an entity's robustness to changes in the input data, a bootstrapping approach was employed. Samples were randomly selected with replacement to build an expression set of 130 samples. Correlations for all pairs of genes were calculated and a NeXO was inferred as described above. This process was repeated to produce 100 bootstrapped ontologies. Each bootstrapped NeXO was then aligned to the original NeXO using the method described in Kramer et al. (2014). Alignment was performed using the strict hierarchy model and all nodes were aligned. The bootstrap score B(e) for an entity (e) was then calculated as:

where n is the number of bootstrapped NeXOs and S_i is the alignment score for entity e when the original NeXO is aligned to the i-th bootstrapped NeXO. Finally, the robustness score for each entity was calculated as a geometric mean of the entity's network support and bootstrap score:

Entities, clusters of genes identified from their correlated expression profiles, potentially important for disease, were identified as those containing more than two members and that show enrichment for genes implicated in disease. Two sources were used to annotate disease genes; (i) 167 disease genes identified from the pathogen-host interaction database (Winnenburg et al., 2006) and (ii) 255 genes that produce secreted and cell wall proteins thought to be associated with disease were taken from Butler et al. (2009). Fisher's exact test was used to identify enrichment by comparing the proportion of genes within each entity associated with disease to the proportion of genes in the whole NeXO associated with disease. The method of Benjamini and Hochberg (1995) was used to control for a false discovery rate of 0.05.

Entities within the NeXO were annotated by identifying enriched GO terms (Ashburner et al., 2000). The GO slim ontology and annotations were downloaded from the CGD and Fisher's exact test was used to identify enriched GO terms in each entity. Briefly, for each entity each annotated GO term was tested for enrichment where the population size was the number of genes in the C. albicans NeXO and the sample size was the size of the entity. The population successes are the number of genes in the C. albicans NeXO (i.e., whole network) annotated with the specific GO term and the sample successes are the number of genes within the entity (i.e., a single cluster) annotated with that term. The method of Benjamini and Hochberg (1995) was used to control for a false discovery rate of 0.05.

Homozygous deletion strains for genes identified in entities associated with disease were assembled from Rauceo et al. (2008) and Noble et al. (2010) in VB GKO plate 1 obtained from the Fungal Genetics Stock Center (McCluskey et al., 2010). A full list of strains used in this study can be found in Supplementary File 1. Strains were stored at −80°C in 96-well plates, replicated on agar plates and subjected to a variety of phenotypic assays including; (i) thermotolerance on YPD at 25, 30, 37, and 42°C, (ii) osmotic stress with YPD+NaCl and YPD+KCl at 30°C at both 0.5 and 1 M; (iii) YPD+H₂O₂ at 30°C with concentrations of 2.5, 5, 7.5, and 10 mM; (iv) acidic environments from pH 4.0 to pH 7.0 at both 30 and 37°C; and (v) YPD+Fluconazole and YPD+Caspofungin at 30°C with concentrations of 1, 2 μg/ml and 0.064 and 0.128 μg/ml, respectively. Strains were incubated under these conditions for 48 h.

Hyphal growth was assessed by first examining colony morphology on agar plates containing YPD+20% serum and in Spider medium incubated at 30 and 37°C for 48 h (Liu et al., 1994). Second, overnight cultures were diluted to OD₆₀₀, incubated in a shaking incubator in liquid YPD containing 20% serum at 37°C and cells were fixed in 10% formalin after 90 m so that hyphal growth could be observed by microscopy.

Bone marrow was harvested from three eight-week old male C57BL/6 mice, from which bone marrow derived macrophages (BMDMs) were differentiated, as described previously (Davies and Gordon, 2005). BMDMs were transferred to a 96-well plate at 2.4 ×10⁴ cells per well. C. albicans cells were grown overnight in YPD and then washed three times in phosphate-buffered saline (PBS) and resuspended in 1 ml PBS. Following counting, 7.2 ×10⁴ of C. albicans cells were added to each well, ensuring a 3:1 ratio of yeast to macrophage. YOYO-1 dye (150 μl) was then added to each well. Fluorescence arising from lysed macrophages was measured at 30 m intervals for 12 h in a FLUOstar Optima plate reader (BMG Labtech) using excitation 485 nm/emission 520 nm at 37°C (Bain et al., 2014).

BMDMs were prepared from three 8-week old male C57BL/6 mice, which were selected randomly, bred in-house, and housed in stock cages under specific pathogen-free conditions. The mice did not undergo surgical procedures prior to culling by cervical dislocation. All animal experimentation was approved by the UK Home Office and by the University of Aberdeen Animal Welfare and Ethical Review Body.

Expression of C. albicans genes from different growth states, environments, stresses, and infection models representing 130 individual samples were used to characterize the C. albicans system. Correlated gene expression profiles for 5,698 C. albicans genes were used to infer a network-extracted ontology (NeXO) with the CliXO method (Kramer et al., 2014). The resulting NeXO consists of 1,366 entities connected by 1,801 relationships where larger entities split into smaller subsets representing functions on multiple scales. The C. albicans NeXO structure can be found in Supplementary File 2.

Aligning this C. albicans NeXO to the GO shows that 1,033/1,366 (75.6%) entities can be aligned to at least one ontology within the GO, with many terms mapping to more than one (Figure 1). Of the 496 entities that contain more than 4 genes, 336 can be aligned to GO terms. The C. albicans NeXO captures 70.6% of terms in the Cellular Component ontology but only captures 18.1 and 20.9% of terms in the Biological Process and Molecular Function ontologies, respectively (Supplementary File 3). This suggests that the C. albicans NeXO largely describes subcellular structures and macromolecular complexes.

Aligning the C. albicans NeXO to the GO provides useful information about the similarity between the structure of the two ontologies. To gain more information about the function, or multiple functions, of entities within the C. albicans NeXO, enriched GO terms for each entity were identified using Fisher's Exact Test (adjusted P <0.05) and used as a functional annotation. There are 333/1,366 (24.3%) entities with significant GO term enrichment and 159/496 (32.0%) entities with significant enrichment when examining entities with greater than four genes. GO term enrichment covers terms from all three GO ontologies with 35, 20, and 23 unique GOslim terms enriched from the biological process, molecular function, and cellular component ontologies, respectively (Supplementary File 4). This suggests that our C. albicans NeXO focusses on a subset of functions across all GO ontologies with many entities showing no functional enrichment. These may represent noise in NeXO generation or the identification of currently unknown or unannotated pathways and functions.

To validate the C. albicans NeXO, robustness scores for each entity were calculated based on network support from co-expression and bootstrapping of the input expression data. We find a strong correlation between term size and robustness with larger terms being more robust (R² = 0.87, P = <0.0001 on log transformed data, Figure S1). Therefore, the larger entities in the NeXO are more robust. Excluding entities with fewer than 10 members yields a mean robustness score of 4.47 compared to 0.82 when all entities are analyzed. There were two clear outliers: these were two entities immediately under the root, which both contain >5,400 genes and are far larger than any other entities in the NeXO (Figure S1), suggesting that these may be spurious entities. Overall, entities within the C. albicans NeXO, particularly larger entities, are robust to the underlying data suggesting that the addition or removal of data wouldn't change their composition and so they are likely to represent biological functions.

Although many entities in the NeXO were annotated by identifying enriched GO terms, many remain unannotated (Figure 1). The unannotated entities may represent understudied, functionally-related genes with roles in pathogenicity that are lacking GO annotation. To identify entities related to pathogenesis irrespective of GO annotation, we identified entities enriched for known disease genes in the pathogen-host interaction database (Winnenburg et al., 2006) (Figure 2). This is a major advantage of using a data-driven NeXO that can complement existing GO annotation but can also discover new entities with roles in specific biological processes.

Eight entities within the C. albicans NeXO are enriched for known disease genes. Six of these eight entities are also enriched for the GO term “pathogenesis” (GO:0009405) indicating that the genes in these entities are related to disease. Additionally, several other terms related to C. albicans disease processes are also enriched in these entities. GO terms “interspecies interaction between organisms” (GO:0044419), “filamentous growth” (GO:0030447), and “growth of unicellular organism as a thread of attached cells” (GO:0070783) are all regularly enriched in putative disease-related entities. However, entities CLX:6939 and CLX:6889 show no GO enrichment and may represent areas of unstudied function associated with disease that have only been detected using this approach.

In those entities that are enriched for known disease genes, many other genes are uncharacterized and could potentially be involved in disease processes. In this way we can predict functional annotations for genes not present in literature-curated sources such as the GO. The entity CLX:6916 (Figure 3), contains 24 genes with 5 genes annotated as being part of the disease process in the pathogen-host interaction database (Winnenburg et al., 2006). Of the other 19 genes in this entity 6 genes have experimental annotation in the CGD and the remaining 11 genes are uncharacterized. CLX:6916 is a highly connected entity with 189 edges with a Pearson correlation coefficient >0.2, suggesting that the genes in this entity may be functionally related. The entity CLX:6916 contains genes involved in the disease process. These include RCN1 (C6_01160W), which encodes a calcineurin-dependent signaling protein that has a role in controlling the stress response and virulence (Reedy et al., 2010). The genes VPS24 (C2_00930C) and C6_02690C have been linked to the processes of adherence Marchais et al. (2005) and hyphal growth (Bensen et al., 2004), respectively. Also, the entity CLX:6857 (Figure S2) contains the gene TRY5 (C6_01500C), which acts as a regulator of adhesion genes (Finkel et al., 2012). These entities display tight clustering suggesting that the genes they contain are functionally related. Since these entities are enriched for known disease genes, we reasoned that the uncharacterized genes in these entities may also play roles in disease and C. albicans pathogenicity.

Entity CLX:6889 (Figure 4) shows enrichment for known disease genes but no significant enrichment for any GO terms. CLX:6889 may represent an understudied process, detectable directly from transcriptomic data using the NeXO approach, where incorrect or missing GO term annotation means we cannot detect any significant GO term enrichment or existing functional annotation. Interestingly, CLX:6939 (Figure S3) shows no GO term enrichment. However, its descendent entities, CLX:6819 and CLX:6036, are enriched for multiple GO terms. Possible interpretations again include the possibility of missing GO term annotation, or that CLX:6939 represents multiple disease-related processes with no coherent GO description while the descendent entities CLX:6819 and CLX:6036 represent singular functions with adequate GO descriptions. In these entities we find that the majority of genes annotated with GO terms have been annotated on the basis of sequence similarity or inferred from electronic annotation, with few terms annotated by mutant phenotypes or from direct assays. This suggests that the GO annotations for these genes may be inaccurate or incomplete. Indeed, ~70% of C. albicans genes are uncharacterized (Candida Genome Database, Nov 2019) and have no experimental evidence for functional annotation (Skrzypek et al., 2016). Together, the enrichment of known disease genes in this entity and the lack of experimental annotation of GO terms, may suggest that these entities, and their genes, represent coherent, but under-characterized, functional units related to disease.

The entity CLX:6889 (Figure 4) contains characterized genes involved in biofilm production (MRPL33—C3_01080W) (Nett et al., 2009) and response to mouse macrophages (RPS14B—C1_06450C) (Lorenz et al., 2004). Both of these genes have not been shown to contribute to the disease process, but have functions associated with disease. CLX:6889 also contains uncharacterized genes that are induced in response to stress (C5_04300C), involved in protein secretion (C5_01900C) and involved in the response to macrophages (THS1) (Nobile et al., 2012). These have inferred functions potentially related to disease. These entities suggest that some genes, that may have been partially characterized or annotated with GO terms, might yet have undiscovered functions involved in the disease process, which might be niche-specific, for example.

To test whether some partially characterized genes might have undiscovered functions involved in the disease process, we selected available, homozygous null mutants corresponding to genes associated with disease-related entities in the NeXO (AHR1, WOR2, CST20, GRR1, HGC1, HOG1, HST7, HWP1, INT1, PBS2, PEP8, RIM101, RSR1, SEF1, TRY5 HIT1, NVL6) and irrespective of GO annotation. The mutants were subjected to a range of phenotypic assays alongside their respective control strains (CAI4+Clp10, RM1000+Clp20, SN250). These virulence-related assays included yeast-hypha morphogenesis, macrophage phagocytosis, thermotolerance, resistance to osmotic, oxidative and acid stresses, and sensitivity to the antifungal drugs, fluconazole, and caspofungin (Figure S4).

There were detectable phenotypes for mutant strains in a range of conditions. Several mutants (HIT1, NVL6, HOG1, PBS2) showed differing levels of sensitivity to YPD+H₂O₂, particularly at concentrations of 7.5 and 10 mM. We also observed potential reduced growth of hit1Δ when exposed to fluconazole and rsr1Δ when treated with caspofungin. Finally, pep8Δ, showed an aberrant colony morphology on Spider medium (Figure 5). Whilst the control strain, SN250, displayed the distinctive wrinkly appearance associated with hyphal development, pep8Δ colonies were smooth, similar to the appearance of hwp1Δ colonies with known hyphal abnormalities, suggesting a potential defect in hyphal development (Figure S4).

To explore the potential defect of pep8Δ cells in hyphal development, we compared their morphology to control cells (SN250) after growth in YPD containing 20% serum at 37°C (Figure 5B). The pep8Δ cells formed stunted hyphae compared to the control. This confirmed that the loss of Pep8 either delays or stunts hyphal development in C. albicans, suggesting that Pep8 plays a role in hyphal development.

The formation of hyphae after phagocytosis is known to promote the escape of C. albicans cells from innate immune cells (Lo et al., 1997; McKenzie et al., 2010). Therefore, we reasoned that, given their attenuated hypha formation, pep8Δ cells might be less able to lyse and escape macrophages following phagocytosis. We tested this by comparing the ability of C. albicans strains to promote the release of the fluorescent dye, YOYO-1, from macrophages (Figure 6). Macrophages infected with wild type C. albicans control strains (SN250 or RM1000+CIp20) released significant amounts of the dye over the time course, whereas uninfected macrophages did not. Also, macrophages infected with the hypha-defective C. albicans cph1Δ efg1Δ mutant did not display significant levels of lysis (Figure 6), which was consistent with the role that hypha formation plays in fungal escape from macrophages (Lo et al., 1997; McKenzie et al., 2010). Significantly, the pep8Δ cells displayed a reduced ability to lyse the macrophages compared to the wild-type controls RM1000 (Mann–Whitney U, P = 5.875e-11) and SN250 (Mann–Whitney U, P = 0.002). However, the pep8Δ cells were able to escape these innate immune cells to some degree (Figure 6). This was entirely consistent with our observation that loss of Pep8 partially attenuates, but does not completely block, hypha formation (Figure 6). We conclude that Pep8 plays a role in hyphal development and defects in this gene's function affect the ability of C. albicans to evade the host immune system.