Pathogenic genes implicated for oral candidiasis were taken from three research articles (Supplementary Material) (Park et al., 2009; Liu et al., 2015; Lemberg et al., 2022) in which C. albicans was cultivated in three oral cell lines. We retrieved 59, 34, and 873 differentially expressed genes (DEGs) of C. albicans from the FaDu cell line GSE5340 (Park et al., 2009), OKF6/TRET cell line GS56093 (Liu et al., 2015), and oral keratinocyte cell line (Lemberg et al., 2022), respectively. We prepared an exhaustive list of 948 pathogenic genes (Supplementary Table S1). To obtain host DEGs, we sourced data from two oral squamous cell carcinoma (OSCC) cell lines, namely, HO-N1-1 and HSC, which had been infected with C. albicans. These data were extracted from Supplementary Material of GSE169728 (Vadovics et al., 2022). We made an exhaustive list of 6,806 significant human DEGs (Supplementary Table S2).

We used STRING v11.5 (Szklarczyk et al., 2019) to build protein–protein interaction networks with the 948 DEGs of C. albicans. These networks contain only physical interactions with a cut-off of interaction confidence score ≥ 0.4. The pPPIN was observed to consist of 559 nodes and 4,167 edges (Supplementary Figure S1). We used Cytoscape 3.9.1 (Shannon et al., 2003) to visualize the networks. Furthermore, we also utilized PathoYeastract (Monteiro et al., 2020) databases for predicting TF–target gene (TG) interactions of 559 nodes of the pPPIN under biofilm-forming environmental conditions. The PathoYeastract database contains experimentally validated (RNA-seq—WT vs. TF mutant, microarray WT vs. TF mutant; chip on chip, ChIP) TF–TG association data (Monteiro et al., 2020). For the extraction of TF–TF interaction, we also used the PathoYeastract (Monteiro et al., 2020) database under biofilm-forming conditions. For TG–TG interaction, we used STRING v11.5 (Szklarczyk et al., 2020) with a cut-off of interaction confidence score ≥ 0.4. We obtained a gene regulatory network (GRN) with 9 TFs and 287 TGs (Supplementary Figure S2).

For constructing a cross-species interaction network, we extracted the whole-genome sequence of C. albicans from Ensembl Fungi (Yates et al., 2021) by choosing the Ensembl Fungi Genes 55 database. Since there is a paucity of confirmed experimental data on protein–protein interactions between C. albicans and host proteins in the context of oral candidiasis owing to the time-consuming and resource-intensive experimental techniques, we were constrained to utilize a prediction-based technique. Although many studies employed protein structure-based predictions to infer host–pathogen protein interactions (Mariano and Wuchty, 2017), this approach cannot assure the discovery of a significant similarity between proteins from pathogens and those with known structures. Thus, to predict through the best possible way, we chose HPIDB 3.0 (Ammari et al., 2016) that aids the annotation, prediction, and visualization of host–pathogen interactions (HPIs) by retrieving all experimentally established pathogen–host PPI data (Ammari et al., 2016). Moreover, this technique was widely used in various research articles to envision the host–pathogen interaction, e.g., Leptospira interrogans and Homo sapiens (Kumar et al., 2019), Burkholderia pseudomallei and Homo sapiens (Loaiza et al., 2020), and Klebsiella pneumoniae and Homo sapiens (Saha and Kundu, 2021). We initially entered the sequences of C. albicans in FASTA format in the HPIDB 3.0. Only those PPIs inferred from experiments such as Y2H, co-immuno-precipitation, and other experimentally robust protocols were considered the template for predicting human–C. albicans PPIs. Since false positive interaction is routinely produced using interolog PPI prediction algorithms, to decrease the retention of false positives in our dataset, we applied a two-step procedure. First, we refine the BLAST alignments and generated random sets of sequence identity, e-value, and sequence coverage combinations, i.e., sequence identity (30%, 40%, 50%, and 60%), e-value (1 × 10⁻⁴, 1 × 10⁻⁵, 1 × 10⁻¹⁰, 1 × 10⁻²⁰, and 1 × 10⁻²⁵), and sequence coverage (40%, 50%, 60%, 80%, and 90%). Next, we chose the combination with sequence identity ≥ 50%, e-value ≤ 1 × 10⁻¹⁰, and coverage ≥ 90% for interolog detection as it conferred an optimal number of interactions. Other combinations provided either profuse or trivial interactions which were not suitable for our downstream analysis. Furthermore, we noticed that these screening criteria had also been utilized as the best reciprocal blast hit criteria for deducing human–Klebsiella pneumoniae interaction (Saha and Kundu, 2021). Subsequently, considering these stringent criteria, we retrieved 822 interactions, consisting of 412 C. albicans and 353 human nodes. Second, we used GO functional similarity enrichment analysis to remove any remaining false positive interolog pairs between the host and pathogen (Saha and Kundu, 2021). We extracted the GO functional annotation of Homo sapiens (host) and C. albicans from the Ensembl database(Yates et al., 2021) and Ensembl Fungi database (Cunningham et al., 2021), respectively. We classified a pair of homologs to be functionally comparable if they shared at least one GO term. We performed the GO functional similarity enrichment analysis using an in-house perl script and received 119 host proteins interacting with 317 C. albicans proteins, constituting 432 interactions (Supplementary Table S3).

For the elucidation of the target network of C. albicans, we extracted 3,888 genes expressed in the adult tongue from Expression Atlas (Mabbott et al., 2013) with TPM value ≥ 11. Then, we reconstructed a tPPIN (Supplementary Figure S3) in stringApp in Cytoscape plugin v3.9.1 (Shannon et al., 2003) by considering physical interaction and a cut-off score of interaction confidence score ≥ 0.4.

We executed a module analysis by using ModuLand (Palotai et al., 2012), Cytoscape plugin v3.9.1 (Assenov et al., 2008). ModuLand helps in determining essential network positions (such as module cores and bridges), physiologically significant groupings, module cores that aid in the identification of biological functions, and inter-modular nodes that play an important role in a range of biological networks (Rodriguez et al., 2020).

We performed GO biological function and KEGG pathway enrichment analyses of the DEGs pPPIN and GRN of C. albicans using YeastEnrichr (Kuleshov et al., 2016) (Supplementary Table S4). GO biological function and Reactome pathway enrichment analyses for human proteins were performed by using Enrichr (Xie et al., 2021) (Supplementary Table S5). YeastEnrichr and Enrichr are web-based enrichment analysis applications that provide a complete range of functional and pathway annotation tools to assist researchers in comprehending the biological importance of lengthy gene lists (Kuleshov et al., 2016; Xie et al., 2021). We considered those GO biological and pathway enrichment terms having adj. p value < 0.05.

The graphical plots were prepared using Cytoscape 3.9.1 (Shannon et al., 2003; Assenov et al., 2008) and R v4.2.2.

In the oral cavity of healthy humans, C. albicans is commonly established as a normal commensal. The host immune system continuously tracks the growth rate of opportunistic pathogens. Under stress conditions, they acclimatize themselves by avoiding host immunological challenges by forming biofilm layers, changing their morphology, i.e., yeast to hyphal, which is considered a pathogenic form as well as rapidly growing filamentous form, assisting the progression of pathogenesis. To explore the pathogenesis of C. albicans during oral candidiasis, we extracted DEGs of C. albicans from three literature supplementary datasets in which C. albicans has infected the FaDu cell line (GSE5340) (Park et al., 2009), OKF6/TRET cell line (GSE56093) (Liu et al., 2015), and oral keratinocyte cell line TR146 (Lemberg et al., 2022). We prepared an exhaustive list of 948 DEGs by taking adj. p value < 0.05 as a cut-off (Supplementary Table S1). Furthermore, to deduce the relevance of the DEGs extracted from oral cell lines in our study, we employed gene ontology enrichment analysis by YeastEnrichr and found that “cellular response to oxidative stress” (GO:0034599), “cellular response to glucose starvation” (GO:0042149), and “glutathione metabolism”-like functions were enriched (Supplementary Table S4). It was demonstrated that glucose starvation induces the expression of pathogenic genes responsible for oral candidiasis to counter the stress (Solis et al., 2023), and in response to oxidative stress induced in the host during oral candidiasis infection, C. albicans utilized the glutathione metabolism pathway (Miramón et al., 2014; Naglik et al., 2017). Moreover, among these DEGs, SAP4, SAP5, and SAP6 are proved to be expressed in the oral epithelia during C. albicans infection in immunocompromised individuals, such as HIV-positive patients, and absent in commensal vicinity (Naglik et al., 2008). Moreover, SAP5 was also proven to help in colonization, penetration, infection, target E-cadherin (a major protein in epithelial cell junction), and evasion in host tissues (Naglik et al., 2008; Meylani et al., 2021). These lines of evidence dictate the robustness of our dataset for pursuing further downstream analysis.

Next, we reconstructed a pathogenic pPPIN, consisting of 559 nodes and 4,167 edges using the STRING database v11.5 (Supplementary Figure S1) (Szklarczyk et al., 2020). Consequently, using YeastEnrichr (Kuleshov et al., 2016), we performed GO biological function enrichment analysis and observed that “cellular response to oxidative stress” (GO:0034599), “regulation of response to endoplasmic reticulum stress” (GO:1905897), and “response to unfolded protein” (GO:0006986) were enriched in the pPPIN (Supplementary Table S4). In response to oral candidiasis infection, antimicrobial peptides (AMPs) such as LL-37 have been released by the host in the oral cavity to challenge pathogen infection, which creates stress in the endoplasmic reticulum, unfolding of proteins, and affects cell adhesion of C. albicans (Tokajuk et al., 2022). In response to endoplasmic reticulum stress, C. albicans regulates the cell wall integrity and endoplasmic reticulum (ER) homeostasis (Hsu et al., 2021) through the Sfp1p protein and aids in oral candidiasis. Moreover, we also performed KEGG pathway enrichment analysis and noticed that “autophagy,” “endocytosis,” “arginine biosynthesis,” “pyruvate metabolism,” and “peroxisome” were enriched in the pPPIN interacting proteins (Supplementary Table S4). Previous reports have indicated that pyruvate metabolism serves as a virulence factor in C. albicans, contributing to localized tissue ketosis in the host as well as compromising the typical defensive action of the host raised via neutrophil activation during OPC (Saeed, 2000; Huppler et al., 2015; Altmeier et al., 2016), while arginine biosynthesis is implemented by the pathogen for hyphal formation to escape host phagocytosis aiding in eliminating C. albicans in the oropharyngeal cavity (Ghosh et al., 2009; Jiménez-López et al., 2013; Černáková and Rodrigues, 2020). Thus, it could be inferred from the functional analyses of the pathogenic protein that they are solely implicated for infection and withstanding the hostile environment in the host.

Now, in response to environmental cues, the transcriptional circuit controls gene expression and its regulation. Therefore, the TFs that regulate the 559 TGs in the pPPIN might play a crucial role in their TG expression under stress conditions. Thus, we constructed the GRN (Supplementary Figure S2), which consists of three types of interactions: TF–TG, TF–TF, and TG–TG. Since it was evidenced that biofilm development plays a crucial role in oral candidiasis (Cho et al., 2021), we searched those TFs which control gene expression during biofilm-forming conditions using the PathoYeastract database. Finally, we constructed a GRN with 300 nodes and 1,232 edges in which 9 TFs control the expression of 287 TGs. Afterward, using YeastEnrichr, we performed GO biological function enrichment analysis of the GRN. We noticed that several stress-responsive functions such as “cellular response to alkaline pH” (GO: 0071469), “cellular response to glucose starvation” (GO: 0042149), “cellular response to starvation” (GO: 0009267), and “positive regulation of transcription from RNA polymerase II promoter in response to stress” (GO: 0036003) were enriched in the GRN (Supplementary Table S4). The functions revealed in our analysis were also found to be relevant in establishing oral candidiasis as it was found that a low oxygen level, decrease in the flow rate of saliva, and pH dysbiosis in the oral cavity help in the expression of the SAP gene in C. albicans to thrive in the oral niche and also help in adherence (Schaller et al., 2005; Vila et al., 2020). Moreover, acidic and alkaline pH decrease the expression of cell wall polysaccharides and interaction with sIgA, which aid C. albicans to proliferate and evade host response (Bikandi et al., 2000). Thus, the pH-sensing pathway Rim101 of C. albicans plays a pivotal role during oropharyngeal candidiasis, and the rim101Δ/Δ mutant significantly reduces the virulence activity of C. albicans in a murine model of oropharyngeal candidiasis (Nobile et al., 2008; Mayer et al., 2013).

Concurrently, we performed KEGG pathway enrichment analysis of the GRN and observed that “Glycolysis/Gluconeogenesis,” “Endocytosis,” “Nitrogen metabolism,” and “Autophagy” pathways were enriched (Supplementary Table S4). It was observed in a previous report (Villar and Dongari-Bagtzoglou, 2008) that during oropharyngeal candidiasis, the macrophage helps identify the mannose ligand of C. albicans through dectin-1. Endocytosis of the pathogen through the macrophage creates a nutrient starvation condition, which makes them shift from glycolysis to gluconeogenesis, and this nutrient flexibility helps them form hypha and aid in escaping from the macrophage (Mayer et al., 2013). Thus, these function and pathway enrichment analyses clearly depict that the proteins involved in the gene regulatory network of C. albicans play a key role in stress-induced pathogenesis during oral candidiasis.

To detect the C. albicans–human cross-species PPI, we created an in silico methodology. First, we extracted the whole-genome sequence of C. albicans from EnsemblFungi (Yates et al., 2021), and afterward, using an interolog approach, we identified the C. albicans–host protein–protein interaction network (CHPPIN) by a rigorous homology search against the host proteome by using the HPIDB 3.0 database (Ammari et al., 2016). We found a putative PPI between 412 C. albicans proteins and 353 human proteins constituting 822 PPIs, using a strict reciprocal blast search. Next, we attempted to demonstrate that the C. albicans proteins involved in interactions with the human host proteins in the interolog PPI are true homologs of pathogen proteins to mediate potential interactions with the host proteins. For this, we used an in silico GO function similarity analysis of the homologous proteins between humans and C. albicans, in which we retrieved the functional annotations from the Ensembl genome browser (Cunningham et al., 2021) and Ensembl Fungi (Yates et al., 2021) for humans and C. albicans, respectively. We observed that 317 of 412 proteins (∼77%) from C. albicans share at least one GO function with the relevant homologs to the human protein. Therefore, it might be concluded that the 317 C. albicans homologs that make up the CHPPIN were likely to be enriched in real functional homologs. We mapped 317 C. albicans proteins with a putative PPI to identify the human counterpart proteins. As a result, we obtained 317 C. albicans homologs interacting with 119 human proteins, constituting 432 PPIs (Supplementary Table S3).

In the next step of the in silico workflow, we attempted to extract those C. albicans proteins from the interologs which have the potential to mediate oral candidiasis as previously mentioned. To pursue this, we mapped 317 C. albicans interologs with 300 proteins of the GRN implicated in oral candidiasis. We obtained 44 proteins intersecting between these two pools of proteins. Thus, it could be stated that these 44 proteins help C. albicans establish a pathogenic repertoire in the oral cavity of the host. We then mapped these C. albicans proteins to the interologs to obtain the corresponding human proteins. Meanwhile, we received 24 human proteins which directly interact with 44 C. albicans proteins. Thus, it is fascinating to know and explore the molecular crosstalk between 44 C. albicans proteins interacting with 24 human proteins, constituting 62 PPIs in the oral mucosa during oral candidiasis (Figure 1).

It is well known that oral candidiasis involves C. albicans infections in the tongue and other oral mucosal regions and is characterized by fungal overgrowth and penetration of surface tissues (Singh et al., 2014; Millsop and Fazel, 2016). Thus, all the 24 human proteins previously mentioned might not be involved in oral candidiasis if they are not expressed in the oral tissues. Thus, we built a tongue-specific PPI network to obtain candidate human proteins in oral candidiasis. First, we extracted 3,888 genes of the tongue with TPM values ≥ 11 (medium expression level) from the Expression Atlas (Mabbott et al., 2013) and then reconstructed the tPPIN comprising 2,775 nodes and 24,889 edges using stringApp in Cytoscape plugin v3.9.1 by considering physical interactors and cut-off score ≥ 0.4 (Assenov et al., 2008) (Supplementary Figure S4). Next, to identify or locate the target network of C. albicans, we mapped 2,775 tPPIN nodes with 24 proteins of the host which directly interact with 44 C. albicans proteins (Figure 2). We acquired a C. albicans-targeted PPIN (CTPPIN), which we can contemplate as a sub-PPI network in the tPPIN, in which we perceived that either of the interacting nodes interacts with C. albicans proteins. Since the CTPPIN comprising 99 nodes and 119 edges is tongue specific, it could be inferred that it might be implicated in oral candidiasis (Figure 2). To substantiate this, we checked the underlying biological processes and metabolic pathways of the proteins of the CTPPIN. We performed GO biological function enrichment analysis using Enrichr (Xie et al., 2021). Figure 3 depicts that the sub-tPPIN targeted by Candida is enriched with several crucial immunological functions and pathways which confer the potential to the host to resist the pathogen infection. Interestingly, apart from the immunological functions, host cell repair functions such as “nucleotide-excision repair DNA damage recognition” (GO:0000715), “apoptotic process” (GO:0006915), “cellular response to hypoxia” (GO:0071456), “modulation by the host of symbiont process” (GO:0051851), “positive regulation of type I interferon production” (GO:0032481), “DNA Double-Strand Break Response” (R-has-5693606), “PIP3 Activates AKT Signaling” (HASHSA-1257604), “Clathrin-mediated Endocytosis” (HASHSA-8856828), “Class I MHC-Mediated Antigen Processing And Presentation” (HASHSA-983169), and “Toll-Like Receptor 3 (TLR3) Cascade” (HASHSA-168164) were enriched in Candida-targeted proteins (Supplementary Table S5). Villar and Zhao (2010) demonstrated that during C. albicans infection, apoptosis occurred for the removal of apoptotic bodies containing pathogens by secondary phagocytes (Doran et al., 2020). The occurrence of epithelial damage and the manipulation of macrophage apoptosis create conditions conducive to fungal colonization and infection (Camilli et al., 2021), which corroborates our observation on the involvement of proteins in the candida-targeted sub-tPPIN. Thus, it could be stated that the proteins engaged in the CTPPIN accomplish the functions required to combat or provide early symptoms to detect the pathogenesis stage of C. albicans in humans.

Deciphering the molecular link between C. albicans pathogenic genes and human genes would provide a holistic view to know the exact scenario during oral candidiasis. To identify the key molecular players in the CTPPIN, we performed module analysis using ModuLand (Kovács et al., 2010) in Cytoscape plugin v3.9.1 (Assenov et al., 2008) and obtained three core modules—DNMT1 (consists of 72 proteins), SGTA (consists of 17 proteins), and TOR1A (consists of 10 proteins) (Figure 2). Next, using Enrichr (Xie et al., 2021), we shortlisted among these module core genes which are well-characterized immunological players. We found that the DNMT1 module significantly enriched (p. adjust <0.05) with 24 immunological function and 85 immunological pathway terms (Figure 4; Supplementary Table S5) consists of 46 proteins. Similarly, the SGTA core module was enriched with 34 immune-related function and 4 pathway terms (Figure 4; Supplementary Table S5), encoded by 13 core module proteins, and the TOR1A module was enriched with seven immunological function and four pathway terms (Figure 4; Supplementary Table S5), composed of eight proteins. From the function and pathway enrichment analyses, we could confer that a considerable fraction of proteins, i.e., 63.8 % in the DNMT1 module, 76.4% in the SGTA module, and 80% in the TOR1A module, are implicated for first-line immunological response against the pathogen (Table 1). Thus, dysregulation in the expression of three core module nodes could hamper the host response against C. albicans during oral candidiasis. To know the expression of core module nodes during a diseased state, we mapped it with 6,806 DEGs extracted from the oral carcinoma cell lines (HO-1-N1 and HSC) GSE169278 (Vadovics et al., 2022) infected with C. albicans (Supplementary Table S2). Subsequently, functional analyses on the DEGs showed the enrichment of “Regulation of Epithelial to Mesenchymal Transition (EMT)” (GO:0010717), “Regulation of Mesenchymal Cell Proliferation” (GO:0010464), “PI3K/AKT Signaling in Cancer” (HASHSA-2219528), “Signaling By WNT in Cancer” (HASHSA-4791275), “Signaling By TGF-beta Receptor Complex In Cancer” (HASHSA-3304351), and “Constitutive Signaling By AKT1 E17K In Cancer” (HASHSA-5674400) (Supplementary Table S5). It was previously reported that infection of C. albicans in the HSC-2 cell line upregulates a set of genes associated with EMT, which in turn help in the progression of OSCC (Krisanaprakornkit and Iamaroon, 2012; Jayanthi et al., 2020; Vadovics et al., 2022). Moreover, the PI3K/AKT pathway, Wnt pathway, and TGF-β pathway play a crucial role in the EMT process and could be considered indicators of metastasis (Chen et al., 2008; Krisanaprakornkit and Iamaroon, 2012; Lakshminarayana et al., 2018; Sasahira and Kirita, 2018). From the aforementioned ontology analysis of diseased condition DEGs, we could decipher that they help in the progression of oral cancer when infected with C. albicans.

Subsequently, after mapping of three core module nodes with diseased condition DEGs, we observed that 54.3% of nodes (n = 46) of the DNMT1 module, 69.23% of nodes (n = 13) of the SGTA core module, and 37.5% of nodes (n = 8) of the TOR1A core module were up- or downregulated in diseased conditions (Table 1). Thus, it could be demonstrated that the deregulation of these crucial network proteins during oral candidiasis facilitates the pathogen to subvert host defense to establish pathogenesis in oral tissues in diseased condition. Now, it would be imperative to examine which of these deregulated nodes of the modules are directly targeted by the pathogen. From the list of potential interologs, we obtained three core module proteins, i.e., ING4, SGTA, and TOR1A, which interact with the HHT21, CYP5, and KAR2 protein of C. albicans, respectively. Intriguingly, it is also evident from our DEG analysis (Supplementary Table S1) that these three genes HHT21 (log2FC = 1.80), KAR2 (log2FC = 2.46), and CYP5 (log2FC = 1.37) of C. albicans are upregulated during oral candidiasis. Searching for the TFs regulating the expression of these three pathogenic genes from the GRN (Supplementary Figure S2) revealed that HHT21 is tuned by Tec1 and Wor1p, CYP5 is regulated by Rob1p, and KAR2 is controlled by Tye7 and Ume6p. Moreover, we perceived that some of the TFs of these three C. albicans genes are up/downregulated during infection (Table 2). Thus, it could be worthwhile to confer that these module-targeting pathogenic proteins, along with their regulators, play significant roles in establishing the infection in the human oral mucosa.