We manually collected information on 150 virulence related factors from peer-reviewed publications on PubMed Literature which have been confirmed in a mammalian model of infection (Figure 1). For each VAF, we found and recorded PMID(s)¹ corresponding to a peer-reviewed article(s) that provided evidence for asserting a Cryptococcus protein as a VAF (Supplementary Figure S1). For each VAF, gene information was identified in the Broad Institute C. neoformans databases² and/or the NCBI Gene website³, and its protein sequence recorded for further analyses.

To predict interaction networks between VAFs, protein sequences of VAFs were put into the STRING database⁴. In the STRING database, the Cryptococcus organism was selected. In the current study, those VAFs related to brain and lung infections were used to generate interaction networks for C. neoformans B-3501. For generating the figures a confidence cutoff of 0.4 was used. The resulting network view, under evidence view, was downloaded for the figure. The database view contains other views of predicted functional partners, for example, neighborhood and fusion of these VAFs. To identify different VAFs involved in MAPK pathways of the pathogenic Cryptococcus species, the enrichment function was used, which contains various pathways found on the Kyoto Encyclopedia of Genes and Genomes.

To determine transmembrane probability for each VAF, its protein sequence was pasted to TMHMM Server v. 2.0⁵) Library in the FASTA format and analyzed using standard settings. To further determine if there would be a difference in predicted membrane location among factors associated with lung vs. brain invasiveness, a separate TMHMM analysis was performed for four groups of VAFs: lung alone, brain alone, factors involved in both brain and lung, and other factors not yet linked to lung or brain invasion. Furthermore, the annotated genome was downloaded from the Broad Institute’s Cryptococcus database and input into STRING-db. The output file was then analyzed using a Python script to log total number of transmembrane proteins and the number of transmembrane helices in each protein.

To identify common motifs from the VAFs, protein sequences were queried against http://www.genome.jp/tools/motif/MOTIF.html Library using Pfam (determining positions and description for each VAF) as shown in Figure 1. Default genome.jp settings were used for this, with an E-value cutoff score of 1.0.

We used C. neoformans var. grubii wild-type (H99 strain) and the sho1Δ mutant strains (YSB1719) (Kim et al., 2015). The sho1Δ::SHO1 complemented strain was constructed as follows. The terminator region of Sho1 was amplified via PCR by using B4796 and B4797 containing NotI, SacII restriction sites, respectively, and cloned into a plasmid pTOP-V2 (Enzynomics). The promoter and ORF region of Sho1 were amplified via PCR by using B4794 and B4795 including SalI, NotI restriction sites, respectively, and cloned into a plasmid pTOP-V2 (Enzynomics). The B4795, containing FLAG sequence (5′-CTTATCGTCGTCATCCTTGTAATC-3′). After sequencing, first the NotI and SacII digested SHO1T was subcloned into a plasmid pJAF12 containing NEO (nourseothricin) resistance marker. Second, the SalI and NotI digested insert was subcloned into the SalI and NotI digested pJAF12-SHO1T to produce a plasmid pJAF12-Sho1-FLAG. The HindIII-digested linearized pJAF12-SHO1-FLAG was transformed into the sho1Δ mutant (YSB1719) and targeted re-integration of the SHO1-FLAG gene into its native locus was confirmed by diagnostic PCR. Functionality of the SHO1-FLAG fusion gene in the new strain (YSB2141) was confirmed by phenotypic analysis.

C57Bl/6 (Jackson Labs, Sacramento, CA, USA) were housed at the Veterans Affairs Ann Arbor animal care facility, and were used at around 7–9 weeks of age. At the time of data collection, mice were humanely euthanized by CO₂ inhalation. All experiments were approved by the University Committee on the Use and Care of Animals and the Veterans Administration Institutional Animal Care and Use Committee.

Mice were anesthetized via intraperitoneal injection of ketamine/xylazine (100/6.8 mg/kg body weight). The trachea was exposed by making a small incision on the throat and moving aside underlying salivary glands and cutting away muscle around the trachea. 10⁴ cells of C. neoformans were inoculated into the lungs by injecting 30 μL of infection mix (see section above) via 30 gauge needle. After inoculation the salivary gland was moved back into place and the incision was closed using a cyanoacrylate adhesive.

Dissected lungs were homogenized in 2 ml of sterile water with protease inhibitor and homogenates spun at 10000 RPM for 10 min to collect supernatants that were stored frozen at -80°C until cytokines were analyzed. Quantification of levels was performed using LEGENDplex cytometric bead array (CBA) and analyzed as described previously (Xu et al., 2016). Lung associated lymph nodes were dissected and processed, as described by Xu et al. (2016). Total RNA from node cells was isolated, cDNA converted and qPCR analysis performed with data normalized to 18S mRNA, as described by Xu et al. (2016).

Lung homogenates were serially diluted in sterile water. Series of 10-fold dilutions of the samples were plated on Sabouraud dextrose agar plates in duplicates in 10-μl aliquots and incubated at room temperature. C. neoformans colonies were counted 2 days later, and the number of CFUs was calculated on a per whole lung basis.

Through manual literature annotation, we identified reports on 150 C. neoformans genes which have an effect as established through in vivo mammalian infection models. This group is defined by experimental evidence of attenuation of virulence in a mouse model. For each gene on our list, the following information was annotated: (1) paper PMID; (2) protein sequence; (3) gene symbol; (4) protein sequence; and (5) for confirmed VAF’s effect on infected organ burden in various articles.

The 150 C. neoformans VAFs identified in this study were all found in STRING db10 from the total identified genes of C. neoformans. The network is composed of the interactions between 150 VAFs that had an impact on the survival of infected mice. Of these VAFs, 8 had documented effects in promoting brain infection, 28 were related to lung infection, 31 genes had dual effects for fungal invasiveness in both organs, and 83 affected survival with target organ partially or completely unknown (Figure 2). In some instances, an assay for brain-virulence was conducted without attention to cryptococcal level or pathology in the lungs, providing evidence for brain virulence but without elucidating the role of the gene in lung infection.

While CryptoNet has already been used to predict protein–protein interactions in C. neoformans virulence genes, STRING offers several additional features. STRING is more liberal with assigning interactions, when compared to CryptoNet, because it uses data from homologous protein interactions in other organisms. A homology comparison feature is incorporated into STRING, which makes it easier to determine the function of un-identified C. neoformans genes. In addition to these features which we used in this study, STRING also has other useful features, including a feature that allows for homology comparisons in a phylogenetic context, and protein family analysis.

Our STRING analysis of C. neoformans B-3501 strain generated networks for 150 VAFs associated with lung and brain virulence. These networks include 410 interactions while 81 interactions could be expected if the collection of genes were randomly chosen among C. neoformans genes (Figure 3). This indicates that VAFs act in orchestration.

The nodes on Figures 3 can be loosely classified into two groups – one group defined by green text-mining interactions and another group of nodes connected by interactions of all colors. The area that is defined by text-mining interactions is primarily composed of proteins involved in metabolism. The area which is densely filled with multi-colored interactions is composed primarily of transcription factors and kinases. Additionally there are many disconnected nodes, which may be loosely regulated constitutively expressed genes or genes whose interacting partners have not been discovered yet in C. neoformans, or proteins whose homologs have not been detected in other species by STRING.

Our STRING analysis of C. neoformans VAFs suggests that there are still many virulence factors to be discovered, and that not all protein interactions have been elucidated. One of the disconnected nodes in the STRING analysis, Cuf1, is a copper dependent transcription factor (Waterman et al., 2007). It is expected that functioning as a transcription factor, Cuf1 regulates transcription of genes including some genes that are involved in mediating pathogenicity. Therefore analysis of genes interacting with CUF1 may help to identify factors that have not been yet recognized as VAFs.

Alternatively, it could mean that STRING has not yet effectively incorporated all possible protein–protein interactions from the literature. For instance, the CAP genes, which are important for synthesis of the Glucuronoxylomannan (GXM) capsule, do not yet have any functional interacting partners in STRING. It has been shown, however, that capsule expression depends on high CO2 and low-iron conditions (O’Meara and Alspaugh, 2012). This would indicate that the CAP genes are interacting partners with CO2 response and iron regulation genes cac1 and cir1, or some intermediate gene. It appears that these interactions were missed by the STRING algorithm when it was mining its sources.

In addition, using the analysis functions of STRING allows us to highlight the importance of kinases and the galactoxylomannan (GalXM) capsule components. The GXM/GalXM capsule is a defining feature of C. neoformans infection and its prominence is highlighted using STRING. The expression of the capsule is regulated by several pathways, discussed above and below, which incorporate numerous kinases. The CAP genes are identified as mannosyl-transferases, crucial for creating cross-links of GXM. The GDP-mannose transporters are primarily important for export and import of the polysaccharide capsule (Cottrell et al., 2007). Finally, four (Man1, Uxs1, Ugd1, and Uge1) of the highlighted sugar metabolism proteins are involved in synthesizing the saccharide components of the capsule.

Additional information about each protein from our database cannot be described in detail here, but can be gained by interactive viewing of our STRING network at⁶ when our database (Supplementary Table S1) is imported.

Using STRING analysis we are able to illustrate a key kinase pathway that is central to the C. neoformans virulence network (Figure 4). Hog1 (located in the center of the interaction chart) is a MAP kinase involved in a signal transduction pathway that is activated by changes in the osmolarity in the extracellular environment (Bahn et al., 2005). Hog1 controls osmotic regulation of transcription of target genes. Hog1 is also involved in response to UV radiation and mediates sensitivity to fludioxonil. Hog1 has several strong predicted functional interactions with: Ste7 (a MAP kinase modulating capsule formation; Feretzaki et al., 2014), Sch9 (a protein kinase Sch9 modulating cell wall formation; Wang et al., 2004), Gpa1 (a G-protein alpha subunit homolog required for melanin synthesis, capsule formation and other functions) (Alspaugh et al., 1997; O’Meara et al., 2010a), and Ras1 (Ras superfamily GTPase that controls thermotolerance) (Waugh et al., 2002; Ballou et al., 2013). The most commonly reported biological functions associated with these proteins are linked to regulation of C. neoformans pathogenicity, although they appear to regulate many other aspects of fungal life cycle.

Analysis of total VAFs using Kyoto Encyclopedia of Genes and Genomics (KEGG) tools in STRING-db indicates that the MAPK pathway is the most prominent among virulence factors. Eight VAFs constitute the MAPK signaling pathway with (p = 8.6 e^-4). The proteins are, Hog1, Ste12alphap, MPK1, Pbs2, Pak1, Pck1, CDC42, and CDC24. The KEGG is a collection of databases dealing with genomes, biological pathways, diseases, drugs, and chemical substances, and has been incorporated into the analysis function of STRING (Kanehisa and Goto, 2000). KEGG has been widely utilized for bioinformatics research and education, including data analysis in genomics, metagenomics, metabolomics, and other -omics studies, modeling and simulation in systems biology, and translational research in drug development (Kanehisa et al., 2010).

In order to more closely investigate this interaction within MAPK gene network the eight regulators of the MAPK pathway were inserted into STRING-database (Figure 5A). The top 20 predicted interacting partners include several genes that are confirmed virulence factors, and many genes which have not yet been investigated in the context of virulence. Among the four functional partners in this group that have already been tested, three are confirmed to be VAFs in vivo at the time of publishing. Most have not yet been named or studied in C. neoformans, however. Of the predicted VAFs that have not been tested in vivo for virulence (Figure 5A) only three have been knocked out in: RHO. STE20, and SHO1. However, RHO has three close homologs, and thus determination of its role in pathogenesis will require extensive studies with many combinations of knockouts. On the other hand, SHO1 is a gene already known to be a protein upstream of the MAPK pathway in another pathogenic fungus Candida albicans, where it serves as cell membrane-bound receptor protein sensing external osmolality (Biswas et al., 2007). STE20 is also a good candidate for study, however, it is down-stream from SHO1 in the HOG pathway (Figure 5B). Considering all these factors, SHO1 was selected as the best gene target for verification of virulence in C. neoformans infection model.

Initial studies lead to the conclusion that Sho1 was absent in C. neoformans. Analysis using STRING and using the KEGG indicated that Sho1 is present in serotype A strain H99 and the serotype D strain B-3501 (Figure 5B), but not present in the serotype D strain JEC21. We suspect that this was the reason why a Sho1 homolog has not been detected in C. neoformans until 2015, when it was characterized in H99 as potential factor in virulence pathways by another bioinformatics study. Kim et al. (2015) conducted an informatics-based genome-wide study to identify novel genes for virulence and adaptation to chemical stresses, pinpointing Sho1 among other VAF candidates. Their study, using A/J mice, has not shown significant effect of Sho1 on the survival rate of infected mice. Due to the prominence of this gene in, and in relation to other genes upstream of, the important HOG pathway, we decided to conduct a study of fungal burdens in an experimental mouse model of cryptococcosis to possibly determine Sho1 role as a factor associated with fungal growth in the lungs.

A common mechanism by which cryptococcal virulence factors confer better survival of the pathogen in the host is through interference with either innate or adaptive host immune response or both (Olszewski et al., 2010). We analyzed fungal burdens in the infected lungs of C57BL/6 mice at time points consistent with the immediate adaptation to the host tissue conditions (day 3), interference with the innate (day 7), and with the adaptive phase of the immune response (day 14). Lungs were collected and homogenized to assay for fungal burdens. Cytokine response was also assessed for possible effects of Sho1 on Th1/ Th2 bias, which is an important index of protective vs. non-protective response used in C. neoformans virulence studies (Osterholzer et al., 2009; Zhang et al., 2010; He et al., 2012; Qiu J. et al., 2012).

This fungal burden assay indicated no significant difference between wild type strain and the sho1Δ and sho1Δ::SHO1 mutants on day 3, but significantly lower fungal burdens in the sho1Δ-infected lungs were detected on days 7 and 14 compared to wild type and sho1Δ::SHO1 infections (Figure 6A). These results indicate that SHO1 expression is not required for the immediate adaptation of C. neoformans to lung environment but it supports optimal growth rate of the pathogen at time points consistent with the innate fungal control (Eastman et al., 2015) and the development of the polarized adaptive immunity (Osterholzer et al., 2009; Eastman et al., 2015), suggesting that SHO1 may have a role in promoting immunomodulatory effects of C. neoformans.

In support of this, we observed greatly enhanced levels of non-protective Th2 cytokines (IL-4, IL-5, and IL-13) in lungs infected with SHO1-sufficient strains compared to lungs infected with sho1Δ at day 14 (Figure 6B). Since the much higher fungal burdens in H99 and sho1Δ::SHO1 at day 14 likely contributed to overall higher cytokine levels in these groups, we analyzed Th1 (IFNγ and TNFα) vs. Th2 (IL-4, IL-5, and IL-13) cytokine ratios generating combined Th1/Th2 cytokine index ratio (Figure 6C). This index reflects overall balance of these cytokines in each infected animal, rather than individual cytokine levels and thus is much less affected by the fungal load. This analysis revealed that the Th1/Th2 cytokine balance was consistently and significantly more skewed toward Th1 in sho1Δ-infected lungs and more toward Th2 in the lungs infected with SHO1-sufficient strains (Figure 6C), suggesting that cryptococcal SHO1-expression contributes to the development of non-protective Th2 bias.

To further support this we assessed gene expression in lung associated nodes on day 7, a time point consistent with early priming of the adaptive immune response in the nodes (Xu et al., 2016). Accordingly, we have found that SHO1-deletion enhanced early expression of protective cytokines IFNγ and IL-17 (Supplementary Figure S1) and resulted in strong trends in upregulation of transcription factors Tbet and RorγT, responsible for the development of protective Th1 and Th17 cells. Thus, cryptococcal SHO1 expression clearly promoted expansion of C. neoformans in the infected lungs and contributed to immunomodulatory potential of C. neoformans.

Both soluble and transmembrane proteins should be important for the pathogenicity of C. neoformans; however, these proteins have not been systematically studied as a group. Multiple in vivo studies confirm that intracellular proteins involved in metabolism and catabolism are vital to C. neoformans proliferation in the lungs and dissemination to the brain (Olszewski et al., 2004; Qiu Y. et al., 2012). On the other hand, the outer capsule is considered one of the defining pathogenicity factors of C. neoformans. Acapsular strains of C. neoformans are generally categorized as non-pathogenic (Okabayashi et al., 2005; O’Meara and Alspaugh, 2012). Also, many studies indicate membrane composition in C. neoformans determines various aspects of pathogenicity, for instance ability to cross the blood brain barrier, and organization of VAFs on the fungal membrane (Siafakas et al., 2006; Singh et al., 2012; Stie and Fox, 2012). Here we classify our collection of VAFs into soluble and transmembrane VAFs using TMHMM.

TMHMM v 2.0 utilizes a hidden Markov model to determine transmembrane domains on proteins, with the ability to distinguish cytoplasmic and outer domains (Krogh et al., 2001), and is currently one of the most accurate membrane protein topology prediction methods. To validate the accuracy of TMHMM, we tested the transporters and transcription factor VAFs of C. neoformans. Because transporters are inherently membrane-associated and transcription factors localize in the nucleus, these are positive and negative controls, respectively (Figure 7A).

Overall, the majority of VAFs showed characteristics of soluble proteins regardless of their association with target organs. Using the TMHMM method, we found that of the 150 known organ specific VAFs of C. neoformans B-3501 77.3% (116) are soluble VAFs and 22.7% (34) are transmembrane (Figure 7A). Among CNS and lung-virulence-related VAFs (31): 81% (26) were soluble and 19% (9) were membrane associated. Among VAFs related to pulmonary infection (28): 68% (19) were soluble proteins, 32% (9) are transmembrane proteins. The VAFs which promote infection in brain (8): consisted of 63% (5) soluble proteins and 37% (3) transmembrane proteins. Finally, among VAFs which have incompletely known target organ (83): 76% (63) were non-membrane proteins and 24% (20) were transmembrane proteins.

Additionally, we conducted a TMHMM study of all annotated C. neoformans genes to elucidate whether VAFs were enriched in transmembrane domains. Of a total 6962 genes, 1354 (19.4%) contain transmembrane domains, compared to 22.6% in our VAF database (Figure 7B). This suggests that both transmembrane and soluble proteins are equally important in the pathogenicity of C. neoformans.

In order to define the functionality of transmembrane VAFs we recorded their functions (Table 1). The largest group is ion sensors and transporters, which play a key role in regulating ion homeostasis and expression of virulence genes within the host (Silva et al., 2011; Jung et al., 2012). Proteins related to mannose and xylose transport and crosslinking are also numerous and these are the proteins required for construction of the immune-modulating capsule of C. neoformans. Another important group of transmembrane VAFs are involved in lipid metabolism, particularly sterols and sphingolipids, which are an important multi-functional membrane component of this fungal pathogen (Siafakas et al., 2006).

To further analyze transmembrane VAFs as a group the number of transmembrane helices present in each protein was determined, in the genome and on our list of VAFs. Although our study indicates that proteins with six transmembrane domains are more frequent among VAFs than the general gene pool (Figure 7B), the proteins in this group function independently of one another and do not possess functional interacting partners in STRING, so we are assuming this is only a coincidence. Although not significant at the statistical level, no VAFs with seven or eight transmembrane domains have been found. In conclusion, the proportion of transmembrane helices in VAFs is similar to that in total genomic proteins.

To further expand our analysis and sort the VAFs based on their expression of functional domains, all 150 organ-specific VAFs were queried against the Motif Search Library (Kanehisa, 1997). Table 2 lists the top common motifs in VAFs of Cryptococcus resulting from this analysis.

The protein kinase (PK) family motif (PF00069) is present in 14 VAFs and is the most common motif among all VAFs. Since C. neoformans is an opportunistic pathogen, it requires efficient upregulation of pathogenicity genes once it enters the stressful host microenvironment, a role partially filled by kinase phosphorylation. Many VAFs from this pathway have been identified, for example, Tyrosine PK, Kinase-like proteins, APH (Phosphotransferase enzyme family) present in MAP kinase pathway, and cAMP-PKA pathway.

Concurrent with the high frequency of kinases, DNA binding proteins are also very frequent among VAFs. Two distinct groups of zinc fingers are found in the VAFs of C. neoformans: GATA and the C2H2, H2C2_2, H2C2_4 proteins. Most of these DNA binding proteins are transcription factors, which is consistent with importance of kinases as regulators in virulence. Efficient upregulation of pathogenicity factors and stress response genes is key to C. neoformans success in evading host immune responses. This might suggest that disrupting the fine-tuned regulatory process in C. neoformans may be a viable alternative to antibiotic treatments.

Finally, we analyzed a less frequent but equally interesting set of genes linked by one motif: the FoxP coiled coil domain, found in three proteins of C. neoformans (Znf2, Zap103, and Rim101), all of which decrease virulence in a mouse model of infection. Although the FoxP domain in mammalian leukocytes facilitates dimerization of FoxP transcription factors and commitment to a T regulatory lineage, Znf2 and Rim101 have very disparate roles and different impacts on cell morphology (Lin et al., 2010; O’Meara et al., 2010b), so it seems unlikely that they dimerize. Zap103 was only recently shown to have impact on virulence in a mouse survival screen (Kim et al., 2015); knockout of Zap103 also leads to increased virulence. The immediate question is whether the FoxP-cc domain is the reason these three proteins are all negative regulators of virulence, or whether it is just an unexpected coincidence. Additionally, as an input for STRING-db, these genes are very effective predictors of other VAFs. Of the top 20 interactors, 9 have a confirmed direct impact on organ fungal burden or survival, while 11 have not yet been tested.