A New Variant of Mutational and Polymorphic Signatures in the ERG11 Gene of Fluconazole-Resistant Candida albicans

Plain Language Summary

Most vaginal infections are very uncomfortable, with symptoms that include itching, pain during sexual intercourse, and soreness in the vagina. Vulvovaginal candidiasis is one such infection in which a serious condition could lead to vaginal cracks. This infection is caused by different species of the microorganism known as Candida, of which the most common is Candida albicans. Several drugs, including fluconazole, have been developed to effectively treat this infection. However, these drugs are not effective in treating some cases because some isolates of this organism have developed a mechanism called mutations that enables them to become resistant to the drugs.

Understanding the cause of this resistance will help health professionals combat this challenge. Hence, many researchers have studied many causes of this resistance but have not sufficiently understood all possibilities. In this study, we used a variety of proven methods to discover a new cause of fluconazole resistance by Candida found in women and dogs. Our results showed that the resistant strain of the organism altered its DNA in a way that affected two vital anchors of the docking area for the drug. Consequently, the drug was rendered ineffective. This knowledge will arm pharmaceutical and medicinal chemists with information to create new designs that fit into such an abridged version of the target or redirect to other more stable spots that are potential drug targets within the organism.

Introduction

The incidence and prevalence of fungal infections have increased since the 1980s, especially in the population of immunocompromised patients.^1–4 One of the most common Candida infections is vulvovaginal candidiasis, with symptoms that include vulvovaginal itching, friction and pain during sexual intercourse, inflammation, soreness in the vagina, burning discomfort when urinating, vaginal discharge, and vaginal cracks.^5–7 According to a survey, vulvovaginal candidiasis affects 70–75% of women of childbearing age at least once in their lifetime,⁸ and 40–50% will experience reoccurrence.⁵ More than 17 different Candida species are known to be etiological agents of human fungal infection, and about 90% of invasive Candida infections are caused by C. albicans, C. glabrata, C. Parapsilosis, C. tropicalis, and C. krusei.^4,9,10 The predominant cause of invasive fungal infection is C. albicans, of which the clinical indicators vary from superficial mucosal infection to deep invasive organ or systemic infection.¹¹

Despite the progress made in the treatment of Candida infections, recent epidemiological studies indicate that the increasing incidence of antifungal resistant isolates of Candida in clinical settings poses a great challenge worldwide, which necessitates the search for new therapeutics.^12–15 Within the limited range of antifungal agents available, the azoles, often referred to as the “workhorse” of antifungal treatment, have been the most frequently used option for Candida infections for decades.^16–20 This is due to their effectiveness, relatively low cost, limited toxicity, and oral availability for systemic antifungal treatment. The azoles are fungistatic^11,16 to Candida, hence lengthy, repeated, and inconsistencies with prophylaxis and treatment courses are the most probable causes of azole resistance in clinical isolates. Azoles work by inhibiting the biosynthesis of ergosterol,²¹ an indispensable component for maintaining the fluidity in the membranes of eukaryotic cells, which leads to the toxic accumulation of its precursor, lanosterol. One of the most commonly prescribed azole antifungals used for the treatment of C. albicans infections is fluconazole.^22–24 Fluconazole inhibits the enzyme sterol 14–α demethylase (CYP51), which is essential for the biosynthesis of the fungal–specific membrane ergosterol.^25,26 C. albicans CYP51 (CaCYP51) whose normal substrate is lanosterol, catalyzes the demethylation reaction of lanosterol in a three-step process towards producing ergosterol. Azole antifungal drugs inhibit this enzyme by binding their nucleophilic N–4 atom to the enzyme’s heme Fe (iron) at the sixth coordinate position, hence occupying the binding pocket competitively.^27–29 CaCYP51 is encoded by the ERG11 gene,³⁰ and altering the nature of the enzyme at the sequence level is a common strategy for evading the inhibitory action of the azoles, leading to resistance.^27,31 Earlier reports have linked ERG11 mutations to three main amino acid sequence hot spots, namely 105–165, 266–287, and 405–488.^32,33 In accord, Kumar et al reported an isolate with triple alterations (Y132C, F145L, A114V), of which Y132 and F145 sites were previously linked to azole resistance.³⁴ They also discovered that one of the isolates had a novel mutation, G129R. In this study, we combined efficacy studies with advances in sequence analysis to characterize the mutational and polymorphic signatures in the ERG11 gene of clinical isolates of C. albicans strains with varying susceptibilities to fluconazole, isolated from women and dogs of reproductive age.

Materials and Methods

Sources of Candida Isolates

Candida isolates were obtained from freshly cultured high vaginal swabs (HVS) samples banked at the repository of the University of Nigeria Medical Center, Nsukka (UNNMC). These samples were collected from anonymous women of reproductive age, between the ages of 20 and 35 years, with symptoms of vulvovaginal candidiasis. Similarly, HVS samples from dogs (1–2 years) were obtained from the repository of the University of Nigeria Veterinary Teaching Hospital, Nsukka (UNVTH). A total of 57 human HVS (h–HVS) and 7 dog HVS (d-HVS) samples were collected, respectively.

Culture Condition and Isolation of Candida Strains

The HVS samples obtained were screened for yeast strains by streaking on Sabouraud Dextrose Agar (SDA),^35,36 Thermo Scientific™ Oxoid™(40 g/L dextrose, 10 g/L peptone, 20 g/L agar, pH 5.6). The SDA was prepared according to the manufacturer’s instructions. To process the specimen, direct streaking with an inoculating loop was used to inoculate samples on SDA. Plates were incubated at 25–30°C in an inverted position for 24–48 h and monitored for growth.

Microscopic Identification and Characterization of Isolates

Isolates were phenotypically identified based on their microscopic and macroscopic characteristics as described by classical methodologies, including colony formation, and the production of germ tubes.^37–39 The clinical samples were incubated in normal saline for 2 to 3 h at 30°C to 37°C. Samples were observed on slides under a light microscope (Motic^TM, USA) and viewed for the short, slender, tube-like structures (germ tube). The presence of germ tubes is indicative of fungal species.

Differential Identification and Characterization of Candida spp. Using Chromogenic Media

The fungal–positive isolates were characterized using differential (Chromogenic Agar) media for Candida;^40–43 Chromogenic Candida Agar (CHROMagar^®; Paris, France).^44–47 The agar plates were brought to room temperature before use. Inoculation was done by streaking the sample onto the plates. The plates were incubated in aerobic conditions at 30–37°C for 48 h and the different species of Candida identified by their colors.^44–47

Antifungal Drug Susceptibility Testing

The susceptibility of the Candida isolates to the selected antifungal agent (fluconazole) was evaluated using disk diffusion test (DDT)^48,49 and minimum inhibitory concentration (MIC),⁵⁰ following guidelines of the Clinical and Laboratory Standards Institute (CLSI) document M27–A316.^51,52 Suspension of inoculums were prepared in 5 mL of sterile saline (0.85%) and the turbidity was adjusted to 0.5 McFarland standards.^53–58 Within 15 minutes of adjusting the turbidity, each isolate was inoculated onto a dried surface of prepared Mueller–Hinton agar plates using sterile cotton swabs. Antimicrobial disks containing 25 μg of fluconazole (Thermo Scientific^TM; Oxoid^TM, USA) were dispensed and pressed firmly onto the inoculated agar plate’s surface to ensure complete contact with the agar surface. Plates were incubated at 37°C for 24 and 48 h and afterward examined. The inhibition zones were measured in millimeters and the results were interpreted using validated CLSI interpretive breakpoints for fluconazole in vitro susceptibility testing.^55–58 Candida species have been classified as susceptible S (zone diameter ≥19 mm), susceptible dose-dependent SDD (15 to 18 mm), and resistant R (≤14 mm). The susceptibility data were acquired in triplicates.

DNA Isolation

ZYMO RESEARCH Fungal/Bacterial DNA Miniprep Kit was used to isolate the Candida DNA following the manufacturer’s instructions. Beta–mercaptoethanol was added to the Genomic Lysis Buffer to a final dilution of 500 μL per 100 mL. In a ZR Bashing Bead^TM Lysis tube, 50–100 mg (wet weight) of Candida cells were weighed into 200 μL of water, followed by 750 μL of lysis solution, and centrifuged at 10,000 x g for 6 minutes. The supernatant was transferred to a Zymo–Spin™ IV Spin Filter in a collection tube and centrifuged further at 7000 x g for 1 minute. Genomic Lysis Buffer (1200 μL) was added to the filtrate in the collection tube. About 800 μL of the mixture was transferred to a Zymo–Spin™IIC Column in a collection tube and centrifuged at 10,000 x g for 1 minute. The Zymo–Spin^TM IIC Column was loaded with 200 μL of DNA Pre–Wash Buffer and centrifuged at 10,000 x g for 1 minute. The Zymo–Spin^TM IIC Column was loaded with 500 μL of g–DNA Wash Buffer and centrifuged for 1 minute at 10,000 x g. The Zymo–SpinTM IIC Column was transferred to a clean 1.5 mL microcentrifuge tube, 100 μL of DNA Elution Buffer was added directly to the column matrix, and the DNA was eluted by centrifugation at 10,000 x g for 30 seconds. The pure DNA was stored at –20°C for further use.

Molecular Identification of Candida Species by PCR–RFLP

Molecular identification of Candida species by PCR–RFLP based on the ribosomal region consisting of internal transcribed spacers ITS1 and ITS2 as well as the 5.8S gene is often used as a confirmatory test for the pure isolate of the organism (Supplementary Table 1).⁵⁹ Amplification was performed in 50 μL volume of 25 μL of 2X PCR Master Mix, 1 μL forward and reverse primers (10 μM each), 1 μL of DNA template (40 ng/μL) and 28 μL water. The universal primer sequences were: ITS1– 5’–TCCGTAGGTGAACCTGCGG–3’, and ITS4– 5–TCCTCCGCTTATTGATATGC–3’. The samples were run on a 96–well Thermal Cycler (Thermo Fisher Scientific) as follows: 95°C (5 min), 95°C (30s), 56°C (30s), 72°C (30s), 37 cycles (step-2 to step-4), and a final extension at 72°C (5 min). Amplicons from the reaction were resolved on 1.5% agarose gel stained with ethidium bromide and analyzed using electrophoresis. A volume of 10 μL amplificons of the ribosomal region were digested with MspI restriction enzyme at 37°C for 2 h. The fragments from the digest were resolved on 2% agarose gel stained with ethidium bromide and analyzed using electrophoresis. The gel was read on BioDoc–It Imaging System (CA, USA). Purity of DNA was verified by amplifying the human beta–globin gene⁶⁰ using the primers PCO3: 5’–ACACAACTGTGTTTCACTAGC–3’ and PCO5: 5’–GAAACCCAAGAGTCTTCTCT–3’.

PCR Amplification of the ERG11

The PCR reaction mixture (25 µL total volume) contained 12.5 µL PCR master mix, 1 μL forward and reverse primers (10 μM each), 1 μL of DNA template (40 ng/μL), and 10.5 μL deionized water. PCR was carried out using two pairs of primers that span the ERG11open reading frame: F1 (5′–AT GGCTATTGTTGAAACTGTCATT–3′), R1 (5-′GGATCAATATCACCACGTTCTC–3′); F2 (5′–ATTGGAGACGTGATGCTGCTCAA–3′), R2 (5’–TTAAAACATACAAGTTTCTCTT–3’). The PCR was performed in a 25–well thermocycler (Eppendorf, Germany). The amplification program for all reactions was as follows: initial denaturation at 94°C (5 min), denaturation at 94°C (30 s), annealing at 55°C (40 s), extension at 72°C (50 s), 30 cycles (step-2 to step-4) followed by final extension at 72°C (10 min). The PCR products were resolved on 2% agarose gel electrophoresis to assess their quality and integrity.

Sequencing the ERG11 Gene

PCR products were cleaned using ExoSAP according to the manufacturer’s protocol. Briefly, the Exo/SAP master mix was prepared by adding 50 µL Exonuclease I (NEB M0293) 20U/ul and 200 µL Shrimp Alkaline Phosphatase (NEB M0371) 1U/ul to a 0.6mL microcentrifuge tube. This was followed by mixing 10µL of PCR Mixture and 2.5µL Exo/SAP Mix. The mixture was properly mixed and incubated at 37°C for 30 min. The reaction was then stopped by heating the mixture at 95°C for 5 min. Gene sequencing was done with the ABI V3.1 Big dye kit following the instructions stipulated in the manufacturer’s manual. The labelled products were cleaned with the Zymo Seq clean–up kit following the manufacturer’s instruction. The cleaned products were injected on the ABI3500XL analyzers with a 50 cm array, using POP7.

Sequence–Level Investigation of the Resistance–Inducing Genetic Variants

In order to understand the observed variability in the pattern of resistance and susceptibility to fluconazole exhibited by the clinical isolates of interest, we employed several structural analytic approaches, which include multiple sequence and structural alignment, in silico site-directed mutagenesis and protein–ligand interaction profiling. The reference sequences and 3-dimensional structures that were used for the purpose of this study were obtained from the National Center for Biotechnology Information (NCBI)⁶¹ and Protein Data Bank (PDB),⁶² respectively. XM_711668.2 represents the accession number of the wild–type C. albicans lanosterol 14–α demethylase (CYP51) nucleotide sequence, while 5FSA and 4WMZ represent the PDB codes for both C. albicans and S. cerevisiae CYP51 proteins, in complex with Posaconazole and fluconazole, respectively.^63,64 To study the variations in the nucleotide composition of the C. albicans drug-resistant and susceptible strains, different multiple sequence alignments (MSA) were conducted using the Clustal Omega software.⁶⁵ The same platform was used for the comparative amino acid sequence of the C. albicans and S. cerevisiae CYP51 in an attempt to evaluate the degree of conservation of the selected residues of interest, while the 3D structural alignment was carried out using the alignment/superposition function of the pymol plugin manager.⁶⁶ In addition, the C. albicans CYP51 mutant was computationally designed using the Chimera–curated Dynameomics rotamer library.⁶⁵ All interactions were predicted using the protein–ligand interaction profiler,⁶⁷ while saved outputs were visualized using the Pymol molecular visualizer.⁶⁶

Molecular Dynamics Simulation

Using the Charmm36 force field and the 2019 version of GROMACS,^68,69 we conducted molecular dynamics simulations on the wild-type and mutant C. albicans CYP51 in complex with their respective prosthetic heme. Solvation was carried out on both systems using the SPC (simple point charge) water molecules, and system neutralization were achieved with the inclusion of adequate chloride and sodium ions. A 5000 step steepest descent was used for the energy minimization of the systems and equilibration was achieved at 1 bar NPT and 300K NVT. With a production run of 100 ns, the comparative analysis of the conformational stability of both systems was estimated in reference to the output of equilibrium properties such as the RMSD (root mean square deviation), RMSF (root mean square fluctuation), Rg (radius of gyration), SASA (solvent accessible surface area), and the intramolecular hydrogen bonding.^68,69

Sequence Retrieval and Phylogenetic Analysis

Additional nucleotide sequences of CandidaERG11gene were retrieved from NCBI–GeneBank⁶¹ for phylogenetic analysis. The criteria for selecting the additional 16 sequences are that they are closest to the four isolates we studied based on NCBI BLAST search, with a cut–off similarity score of 95% and 100% query coverage, they have 1587 bp coding sequence length, they only differ from the wild-type by at least a point mutation, and have been implicated in azoles susceptibility studies. The construction of a phylogenetic tree to depict the evolutionary relationship between the Candida strains based on the ERG11 sequence was carried out using Geneious Prime software⁷⁰ and Molecular Evolutionary Genetic Analysis Software 7.0 (Mega7).⁷¹ The analysis was done using the Maximum Likelihood method (–10997.03) and the Tamura–Nei model with 1000 bootstrap replications.

Statistical Analysis

Other statistical analysis and building of graphs and charts were done using GraphPad Prism version 8.1.0 (221).

Results

Identification and Characterization of Candida Species

Approximately half of the h-HVS samples grew on SDA with a germ tube, indicating fungal presence (Figure 1A and C; Supplementary Table 2), and a similar result was obtained for d-HVS samples (Figure 1B and C; Supplementary Table 2). After considering the phenotypic characteristics of the dog samples, such as growth rate, color, and colony formation (observation not shown), the d-HVS–001 sample showed the right characteristics that necessitated its pick for further studies. We further characterized the Candida isolates at the species level, which was identified by varying colors on chromogenic (differential) media (Figure 1D), and the result showed different species of Candida, including C. albicans, C. tropicalis, C. krusei, C. glabrata, and other unidentified species.^44–47 Some displayed mixed cultures, while others were pure cultures. The result of the study showed mostly pure cultures and a relatively lower percentage of mixed cultures (Figure 1E). Species-level screening of the 38 pure isolates obtained from the h-HVS showed a predominantly C. albicans species, followed by C. glabrata, C. krusei, and C. tropicalis (Figure 1F). The result further showed that the percentage of other unidentified species was as high as 37.84%. We identified the isolate from the dog (Can Iso–dog/029) as a C. albicans strain.

Figure 1 Prevalence and species characteristics of Candida isolates in h-HVS and d-HVS samples. (A and B) represent the prevalence of yeast growth in h-HVS and d-HVS samples from women and dogs, respectively showing symptoms of vulvovaginal candidiasis. (C) A representative image of germ tube indicative of fungal species. Arrows show mother cell (blue), daughter cell (yellow), and germ tube (white). DIC image was taken using a 40X objective at a scale of 10 μm using an upright light microscope (MoticTM, USA). (D) Representative images of the chromogenic display of different strains of Candida on the differential media. Arrows indicate different color variations spanning C. albicans → green; C. tropicalis → metallic blue; C. krusei → pink, fuzzy; C. glabrata → mauve–brown; other species → white to mauve. (E) Some samples have mixed cultures of isolates while some were pure cultures. (F) The species distribution of the isolates is shown, with C. albicans exhibiting the highest prevalence.

Isolates Exhibited Varying Susceptibility to Fluconazole

Antifungal susceptibility testing on the h-HVS samples revealed that the majority of the Candida isolates were susceptible (≥19 mm) to 25 g fluconazole, with a lower proportion being resistant (≤14 mm) and SDD (15–18 mm) strains (Table 1; Figure 2A–G). As shown in the result, the isolate from dog was susceptible. Can Iso–001, a C. albicans strain, was the most resistant (0 mm), and Can Iso–028, also a C. albicans strain, was the most susceptible. We extended the drug exposure time to 48 hours (Table 1), and compared it with the outcome at 24 hours. However, the result showed that there was no significant difference (p > 0.05) between the zone of inhibition and MIC values exhibited by the Candida isolates on 25 μg fluconazole at 24 h and 48 h (Figure 2H and I, respectively).

Table 1 Susceptibility of Candida Isolates to 25 μg Fluconazole

Figure 2 Representative images of C. albicans isolates displaying a varying degree of susceptibility to 25 μg fluconazole including (A) SC5314, (B) the most resistant isolate, (C) SSD isolate, (D) the most susceptible isolate, and (E) the isolate from a dog. (F) and (G) shows the fluconazole susceptibility distribution between the Candida species. (H and I) are heatmap images of the comparative analysis between the Candida strains susceptibility to 25 μg fluconazole at 24 h and 48 h using zone of inhibition and MIC, respectively. Student’s t-test analysis shows that there was no significant difference (p > 0.05) between the susceptibility at 24 h and 48 h.

Disparities in the ERG11 Gene Sequence Elucidates the Varying Susceptibility to Fluconazole

The ERG11 gene codes for the Candida protein lanosterol 14–α demethylase (CYP51).⁷² Four (4) isolates from the experimental isolates, in addition to SC5314 were the strains of interest for ERG11 gene sequencing, including the most resistant (Can Iso–001), SSD (Can Iso–17), the most susceptible (Can Iso–028) and the dog isolate (Can Iso–029/dog). The nucleotide sequences of the ERG11 gene of the various samples showed that the primary sequence structure varied. It has been previously reported that the inhibition of fungal CYP51 by azoles occurs through the formation of an axial coordination bond between the azole drugs and the prosthetic heme iron of the protein, thereby affecting the reduction potential of the iron moiety.⁶³ To further demonstrate this as related to this study, we analyzed the binding modes and interactions of two azole drugs (posaconazole and fluconazole) that have been complexed with the binding pocket of CYP51 in C. albicans and S. cerevisiae, respectively. A keen study of the output by the protein–ligand interaction profiler revealed various forms of interaction between the bound azole drugs and different residues of the CYP51 binding pocket and, most importantly, with the prosthetic heme iron (Figure 3A and B; Supplementary Tables 3 and 4).⁷²

Figure 3 A depiction of the binding poses of two different azole drugs in the CYP51 drug-binding pocket. (A) shows the interaction of posaconazole (yellow) with different active site residues (blue) and prosthetic heme group (brown) of the C. albicans CYP51. (B) shows the interaction of fluconazole (green) with different active site residues (blue) and prosthetic heme group (brown) of the S. cerevisiae CYP51. Notes: A similar study by Zhang et al serves as additional evidence/support for our original analysis.72

Several studies have reported that the prosthetic heme group of the fungal CYP51 protein is anchored by different active site residues that when mutated have the potential to disrupt the position of the heme, hence conferring drug resistance on the mutated fungal species.^72–79 Based on this knowledge, we conducted a MSA between the nucleotide sequence (ERG11) of the wild-type C. albicans CYP51 protein and the nucleotide sequences of the mutated clinical isolates. Different levels of nucleotide and protein sequence similarity studies were conducted in this regard, in order to separate the resistant strains from the susceptible ones. The first alignment was directed towards the detection of specific mutations in the nucleotide sequences that code for the heme-binding amino acid residues in the C. albicans CYP51 active site. This was conducted on the assumption that mutations in the drug-binding site residues of CYP51 that reduce the stability of the prosthetic heme group affect the binding affinity of antifungal azole drugs, thus conferring resistance to various degrees. We aligned the nucleotide sequence of the wild-type CYP51 (XM_711668.2) with the sequences of the most resistant and SDD clinical isolates (Can–iso–001 and Can–iso–017, respectively). The results of this alignment study revealed two mutated codons (LEU–370 and PRO–375) that code for heme–binding residues in the active site of C. albicans CYP51 (Figure 4A and B; Supplementary Figures 1 and 2). We, therefore, speculated that the observed differences in the degree of resistance by both strains are linked to the mutation of these two residues to SER–370 (in both strains) and HIS–375 (only in Can–iso–001), as revealed in the amino acid sequence alignment.

Figure 4 CYP51 nucleotide and amino acid MSA for the wild-type (XM_711668.2) and clinical isolates of C. albicans. (A) Nucleotide MSA between the wild-type (XM_711668.2) and resistant strains (Can–iso–001 and Can–iso–017). (B) Amino acid MSA between the wild-type (XM_711668.2) and the resistant strains. (C) Nucleotide MSA between the wild-type (XM_711668.2) and susceptible strains (Can–iso–028 and Can–iso–dog). (D) Amino acid MSA between the CYP51 of the wild-type and the susceptible strains. (E) Alignment of the amino acid sequences of C. albicans (5FSA) and S. cerevisiae (4WMZ) CYP51. The highlighted residues (LEU–370 and PRO–375 in the C. albicans CYP51) are the amino acid residues of interest and both are equivalent to LEU–374 and PRO–379 respectively in the S. cerevisiae CYP51. (F) Structural alignment of the C. albicans (grey cartoon) and S. cerevisiae (turquoise cartoon) CYP51 3D structures as visualized using the Pymol molecular visualizer. The prosthetic heme molecule of the C. albicans CYP51 is colored in purple while that of the S. cerevisiae is colored in blue. The highlighted regions depict the codons and amino acid residues of interest while the asterisks denote the degree of conservation for within residues.

We then asked if the nucleotide sequences from the susceptible strains exhibited any of the mutations that were observed in the resistant strains. For this reason, a separate MSA was conducted to include only nucleotide sequences of the wild-type and the susceptible strains of the C. albicans CYP51 (Can–iso–028 and Can–iso–dog). Our observation, which supports the earlier speculation, shows that the mutations of the CTC and CCA codons are directly linked to the drug-resistant attributes of the Can–iso–001 and Can–iso–017 isolates, as neither of the two susceptible strains (Can–iso–028 and Can–iso–dog) exhibited these mutations (Figure 4C and D; Supplementary Figures 3 and 4).

Additionally, we performed a pairwise alignment between the amino acid sequences of both C. albicans and S. cerevisiae CYP51 proteins and confirmed the conservation of these heme-binding residues of interest (Figure 4E; Supplementary Figure 5). To better understand the dynamics of the mutations linked to drug resistance, we carried out a structural alignment study on the 3D structures of C. albicans and S. cerevisiae CYP51 accessed from the Protein Database (PDB). The structural alignment showed a high degree of similarity in the binding orientation of the prosthetic heme group in both proteins (Figure 4F, Supplementary Figure 6A and 6B).

Using the Dynameomics rotamer library,⁸⁰ we created an in-silico mutant of C. albicans CYP51 by replacing serine and histidine with LEU–370 and PRO–375, respectively. The result reveals a difference in the orientation of these heme-interacting residues in the wild-type (Figure 5A) and the mutant (Figure 5B). Additionally, this site-directed mutagenesis protocol was conducted in an attempt to also evaluate the effect of both mutations on the observed hydrophobic interaction with the prosthetic heme group. Upon the mutation of both residues to serine and histidine, respectively, the hydrophobic interaction between the wild-type residues (LEU–370 and PRO–375) and the prosthetic heme moiety was lost (Figure 6A and B). This interaction loss speculatively leads to the destabilization of the prosthetic heme, resulting in an eventual increase in resistance to the heme–binding antifungal azole drug (fluconazole).

Figure 5 Wild–type and mutant C. albicans CYP51 drug–binding site. (A) depicts the binding orientation of the prosthetic heme (purple) and the two wild-type residues of interest (LEU–370 and PRO–375) in yellow sticks near the drug-binding pocket (Orange surface representation). (B) shows the binding orientation of the prosthetic heme (purple) and the two mutant residues of interest (SER–370 and SER–375) in red sticks near the drug-binding pocket (Orange surface representation).

Figure 6 A view of the C. albicans CYP51 residues interacting with the prosthetic heme in the drug–binding pocket of the protein. (A) The wild-type CYP51 shows residues LEU–370 and PRO–375 as part of the heme interacting residues. (B) Mutation in these residues to SER–370 and HIS–375, respectively led to the loss of the hydrophobic interaction with these residues hence, these residues are lost in the interaction profile.

The Resistance-Inducing Mutation Stabilized the CYP51 Protein at the Expense of the Prosthetic Heme

We used a comparative molecular dynamics simulation to study the movement of molecules and atoms within the wild-type and mutant C. albicans CYP51 proteins over a period of 100 ns after observing the interaction disparities caused by mutations affecting the prosthetic heme-binding pocket residues. This is to better understand their conformational stability changes. In molecular dynamics simulation, the molecular and atomic trajectories are determined through the numerical simplification of Newton’s equation of motion for systems of interacting particles.^68,69 Forces between potential energies and their corresponding particles are usually estimated using molecular mechanics force fields or interatomic potentials.^68,69 In this method, the RMSD measures the differences in the protein backbones from the initial to the final structural conformation.^81,82 The protein conformation relative to its stability can be estimated through the observed deviations during the simulation. The larger the deviations, the less stable the protein structure.^68,69 The RMSD plots showed that the wild-type protein exhibits a gradual increase in deviation up until 75 ns, after which the system stabilizes until the end of the simulation (Figure 7A). A different trajectory was observed in the mutant protein, which experienced a gradual increase in deviation till 30 ns, after which the system equilibrated and stabilized at an RMSD of 3 nm till the end of the simulation time. The regions of both proteins (wild-type and mutant C. albicans CYP51) that display a high degree of flexibility were also estimated using the residue-based RMSF. A comparative study of the RMSF plots showed that the mutant protein displayed fewer flexible residues compared to the wild-type (Figure 7B), indicating increased stability upon mutation.

Figure 7 Analysis of the molecular dynamics simulations and trajectory of the wild-type and mutant proteins. (A) plots of the root mean square deviation, (B) root mean square fluctuation, and (C) the radius of gyration for the wild-type and mutant proteins, respectively.

The packing of amino acid residues is also known to have a great impact on the stability of a protein,⁸³ and the measure of protein compactness has often been estimated using the Rg.⁸⁴ We observed a lower value for the mutant protein’s radius of gyration (approximately 2.2 nm) compared to the wild-type (approximately 2.4 nm) (Figure 7C). This result is consistent with the outcome of the RMSD and RMSF, which indicates its increased stability over the wild-type. The low fluctuation as observed in the Rg plot was also maintained throughout the simulation. The compactness of a protein has also previously been defined as the ratio of its accessible surface area to the surface area of an ideal sphere of the same volume.⁸⁵ Using the SASA plot, we estimated the solvent-exposed surface area of both the wild-type and mutant proteins.^86,87 A comprehensive study of both plots reveals a lower degree of solvent accessibility upon mutation (Supplementary Figure 7), hence a higher degree of compactness and structural stability in comparison to the wild-type protein.

We also compared the RMSD trajectories of the bound prosthetic heme of the wild-type and mutant proteins. The wild-type heme equilibrated at about 60 ns with an approximate RMSD of 3 nm, unlike the mutant heme that maintained constant fluctuations all through the simulation period (Figure 8A). Protein-ligand hydrogen bonding provides specificity and directionality of interaction, which is a crucial molecular recognition aspect.⁸⁸ While the wild-type heme ligand exhibited a drastic increase in hydrogen bonds throughout the simulation, the mutant heme experienced a drastic reduction in hydrogen bonds from 60 ns till the end of the simulation period (Figure 8B). This is in contrast to the wild-type CYP51 protein, which displays less hydrogen bond stability when compared to that of the mutant (Supplementary Figure 7). This result, supported by the ligand RMSD output, both suggest a reduction in the stability of the bound prosthetic heme upon the mutation of the active site residues (LEU–370 and PRO–375).

Figure 8 Analysis of the molecular dynamic simulations and trajectory of the wild-type and mutant bound prosthetic heme. (A) Plots of the root mean square deviation, and (B) the intramolecular hydrogen bonds for the wild-type and mutant bound prosthetic heme, respectively.

Evolutionary Relatedness of the Candida Strains Based on the ERG11 Gene Sequence

We analyzed the ancestrally unrooted phylogenetic relationship between 22 different C. albicans isolates based on the ERG11 gene sequence, and an additional ERG16 sequence that came up as a significant hit to all the ERG11 sequences of our isolates in this study (Figure 9). These strains include XM_711668 (SC5314), Can–iso–001, Can–iso–017, Can–iso–028, Can–iso–029, and 16 other pre-selected isolates from previous studies deposited in the GenBank, as well as the C. albicans ERG16 sequence (Supplementary Figure 8; Supplementary Table 5). The study showed that the new isolates studied here are categorically distinct from the other isolates analyzed. Our findings also showed that the most resistant isolate, Can–iso–001, is most closely related to the susceptible isolate from a dog (Can–iso–029). Our results further showed that the most related azole-resistant isolate to our isolates is HM194202.1 (SZ169) isolated in China, with the following nucleotide changes; T693A, C341T, and G1400R (a non–specific nucleotide encoded by the R). The least related isolate is AF153849.1 (J913004/1) isolated in Belgium, with the following sequence changes; A348T, C383A, C1355T, and A1390G. All these spots are different from the drug resistance-related position discovered in this work. We also conducted an XM_711668 (SC5314) ancestrally rooted phylogenetic analysis of the same set of sequences. We found that Can–iso–001, Can–iso–017, Can–iso–028, and Can–iso–029, are distant from SC5314 by 9 strains but are closer to the SC5314 than seven 7 strains of the 16 additional strains (Supplementary Figure 9).

Figure 9 Unrooted phylogenetic relationship between the ERG11 genes of the 4 Candida isolates sequenced in this study (in purple), and 17 pre-selected nucleotide sequences from GenBank repository. Wild–type is in brown, and ERG16 in pink.

Discussion

This study investigated the species distribution, azole susceptibility, ERG11 gene profile, implications for the CYP51 protein properties, and the molecular phylogenetic relationship (based on the ERG11 sequence) of Candida strains isolated from the vagina of women and dogs of reproductive age. From the h-HVS samples, we found out that about half of the patients showing symptoms similar to those of vaginal candidiasis do not have Candida infection as the swab from such showed no growth or germ tube typical of Candida (Figure 1A and C). A similar observation was obtained for the d-HVS samples (Figure 1B and C). This implies that other infections may display the same symptoms as those exhibited in vulvovaginal candidiasis. Some of the samples showed mixed cultures, while others showed pure cultures. Different species of Candida, including C. albicans, C. tropicalis, C. krusei, and C. glabrata, were identified in this study (Figure 1D), and many other isolates could not be identified by the differential method employed. We showed that most of the samples from women that had Candida infection had pure cultures and a lower number of mixed cultures, not exceeding three species (Figure 1E). A similar study identified mixed cultures in h-HVS samples,⁸⁹ however, the study did not quantify how many isolates are in mixed cultures. Mixed cultures can further complicate drug efficacy, as slight molecular variations may exist between the different species at the drug target site.

Species characterization of our isolates by differential media showed C. albicans to be the predominant species, followed by C. glabrata, C. krusei, and C. tropicalis (Figure 1F). The proportion of the isolates that were unidentified NAC species was equated with the proportion of C. albicans. The relationship between the C. albicans and the other unidentified NAC species may be difficult to explain as the method of identification used in this study is limited. The only successful isolate from the dog was identified as a C. albicans strain. Other studies have reported a similar distribution as seen in our study. For example, Khadka et al reported 56% C. albicans (as the predominant species), 20% C. tropicalis, 14% C. glabrata, and 10% C. krusei in clinical isolates from various specimens including HVS, urine, blood, sputum, tracheal tubes, and catheter tips in a hospital in Nepal.⁹⁰ Likewise, related studies also reported a comparable prevalence to our findings.^91,92 In contrast to our study, which showed that C. glabrata is the most predominant non–albicans (NAC) species, a study reported that C. tropicalis was the most prevalent NAC species in clinical isolates obtained from a hospital in Hyderabad.⁹³ Combining the prevalence of the other unidentified NAC species and the NAC (identified) species, we show that NAC has a greater prevalence than C. albicans, and other studies have corroborated this observation,^94,95 making the emergence of NAC species of important concern. The most predominant species in mixed culture with C. albicans were other unidentified NAC species, explaining the substantial thriving of C. albicans along with other species (Table 1). This indicates that there are diverse Candida species that are harbored in the vagina of women apart from the common species such as C. albicans, C. krusei, C. glabrata, C. tropicalis, and C. parapsilosis. The significance of this finding is that other unidentified Candida species, such as the recently trending C. dubliniensis,⁹² for which resistance to fluconazole has also been identified^96–98 may substitute for C. albicans under favorable conditions.

Antifungal susceptibility testing of Candida is largely recognized as a useful aid in optimizing the treatment of Candida infections.^51,52,57 This is essential due to the emergence of resistant strains that continuously threaten azole therapy as recently reported in a surveillance study on fluconazole resistance in both albicans and NAC species.⁹⁹ Fluconazole is often the first choice for treating vaginitis in many parts of the world.^100–104 In our study, the antifungal susceptibility testing on human isolates revealed that most of the Candida isolates were susceptible to 25 μg fluconazole (Table 1 and Figure 2). This observation shows that fluconazole remains a very useful antifungal drug for Candida infection. Overall, only about 5% of the isolates from humans investigated were resistant to fluconazole, and about 5% were SSD. This low percentage suggests that the need for new antifungals for Candida infection, though very important, may not be pressing currently. The SDD strains could also become susceptible by increasing the dosage or usage frequency of fluconazole.¹⁰⁵ Our results showed that the most resistant isolate, a C. albicans strain, was completely resistant to fluconazole. We speculated that this strain may have undergone mutations as well as evolutionary adaptation in the ERG11 gene such that it has countered all possible mechanisms of fluconazole action.^{24,26,27,31,75} However, other resistance mechanisms (not investigated in this present study) apart from mutations in ERG11 may also be responsible.^106–108 We extended the duration of the susceptibility testing from 24 h to 48 h to check for any significant improvement in efficacy.^109,110 Our findings showed that there was no significant difference (p > 0.05) between the zone of inhibition at 24 h and 48 h (Figure 2H and I, respectively).

Four (4) isolates, in addition to the SC5314 strain, were the strains of interest selected for ERG11 gene sequencing and profiling, comprising of the most resistant isolate (Can Iso–001), an SDD isolate (Can Iso–17), the most susceptible isolate (Can Iso–028), and the dog isolate (Can Iso–029). The sequence analysis of the ERG11 gene for the various isolates showed that the nucleotide compositions of each strain varied. Differences observed between sequences of strains from the same species could be attributed to mutations or single nucleotide polymorphism (SNPs), and may also be due to changes resulting from counter mechanisms against antifungal agents, as the case may be.^111–113 Variations in the ERG11 gene sequence have very diverse effects on the corresponding protein properties, including stability, electrostatic effects, protein packing, local and global structure, interactions, abundance, catalytic function, protein recognition, transport function of protein and other properties.^111–113 The effects on protein tertiary structure further includes effect on protein dynamics, which covers those affecting allosteric sites, structural disorder, induced fit, and structural flexibility of the enzyme.^114–116 Many deleterious variants, such as large deletions, protein truncations, and amphigoric amino acid insertions, are straightforward to explain. Usually, the most difficult ones are minor sequence alterations, most often resulting in amino acid substitutions in the protein. Therefore, variations in primary, secondary, and tertiary structures of CYP51 proteins (encoded by ERG11) are more difficult to explain than those of the ERG11 DNA sequence itself.

Previous studies have reported several mutations in the ERG11 gene of C. albicans that are linked to resistance to fluconazole and other azoles–based antifungal drugs in clinical isolates.^{27,33,117,118} The exploration of these mutations and their impacts could therefore provide major insights into the mechanisms that underlie their resistance to the azole drugs. In–silico techniques in the modern era play a crucial role in our understanding of biological systems.¹¹⁹ As part of these techniques, mutations are often generated in the proteins in order to elucidate the effects on the stability and conformation of the tertiary or quaternary structures.⁸² These techniques are currently used in the prediction of different biological mechanisms, such as protein–protein interactions (PPI), DNA-protein interactions, and drug resistance and sensitivity.⁸² The substitution of different amino acids in the C. albicans CYP51 protein is one of the most critical mechanisms that contribute to azole drug resistance in the organism.^{27,33,117,118} Studying the effect of residue substitution in this protein is therefore crucial to the identification of important residues that may be involved in the resistance of C. albicans to drugs.

Fluconazole has been used for more than 30 years as the major therapeutic agent for the treatment of C. albicans–linked infections, and there is substantial data available on clinical isolates of C. albicans that have displayed varying degrees of resistance to the drug.^{12,13,74,75,106,108,118} Various mechanisms of acquired resistance towards antifungal azoles have been proposed, and mutations in CYP51 are listed among such mechanisms, although their contribution to the resistance phenomenon remains inadequately understood.¹²⁰ Considering the close structural evolutionary similarity of the CYP51 family across fungal species, and the property of the protein towards the preservation of their conserved biological function,^121–124 it is unlikely to experience a high mutational frequency in the protein. Earlier studies have summarized that most of the reported mutations in the CYP51 from strains of C. albicans that displayed resistance against fluconazole have also been reported in susceptible strains, possibly making them irrelevant to the mechanism of resistance.¹²⁰ Morio et al reported 10 mutations in the C. albicans CYP51 structure, which has so far not been detected in sequences of susceptible strains and 5 of these mutations, which include Y132H, Y132F, K143R, G307S, and S405F were reported as part of the exposed residues in the drug-binding pocket of the protein, suggesting that they might potentially affect the binding affinity of fluconazole. In agreement with our protein–ligand interaction prediction report (Figure 3A and B; Supplementary Tables 3 and 4), Morio et al reported that 2 of the 5 listed active site residues of the C. albicans CYP51 structure (TYR–132 and LYS–143) form hydrogen bond interactions with the bound prosthetic heme group. Abolishing the hydrogen bond interaction or changes in the hydrogen bond strength as a result of the mutation of both residues to histidine and arginine, respectively, are suggested to have a possible effect on the heme iron redox potential;¹²⁰ thereby affecting the ability of the iron to coordinate the basic nitrogen of the azole ring.¹²⁰ Several mutations (spontaneous or SNPs) were observed in the nucleotide sequence of the four clinical isolates that were used in this study, but as a result of the reported relevance of the heme-binding active site residues to the drug resistance attribute of the C. albicans CYP51,^{29,74,77–79,116} our study was targeted at these set of residues in order to investigate the impact of mutations on the observed variability in resistance and susceptibility to fluconazole.

Through MSA protocols, we detected two codons in the nucleotide sequence of the wild-type (CTC and CCA) that codes for LEU–370 and PRO–375. An alignment of the wild-type sequence (XM_711668.2) with sequences of the resistant clinical isolates (Can–iso–001 and Can–iso–017) revealed that these residues were mutated to serine and histidine, respectively (Figure 4A and B). In order to validate the relevance of the observed mutations to drug resistance as previously reported, we repeated the alignment using the wild-type nucleotide sequence and the nucleotide sequences of the drug-sensitive isolates (Can–iso–028 and Can–iso–dog). The output of this alignment (Figure 4C and D) showed that the codons for the heme-binding residues were not mutated in the susceptible strains, which further supports our observation. So far, the L370S and P375H resistance–linked mutations have not been previously reported in any literature. The alignment of their amino acid sequences and 3D structural alignment between the C. albicans and S. cerevisiaeCYP51 showed a high degree of conservation both generally and for the two residues of interest (LEU–370 and PRO–375) as well as a conserved binding mode of the prosthetic heme group (Figure 4E and F). To further elucidate the relevance of these residues to the observed resistance to fluconazole by the clinical isolates, the output from our site-directed mutagenesis study also revealed a loss of hydrophobic interaction between the residues of interest and the prosthetic heme group (Figure 5). In the wild-type caCYP51, LEU–370, and PRO–375 form part of the heme interacting residues however, the mutation in these residues to SER–370 and HIS–375, respectively, led to the loss of these residues in the interaction profile (Figure 6). This possibly leads to the heme moiety not being effectively held in place.¹²⁰

More insight into the mechanism of resistance of the mutant was further revealed through the molecular dynamics simulation, involving the comprehensive analysis of the RMSD, RMSF, Rg, SASA, and the hydrogen bond plots of the wild-type and mutant protein–ligand interaction. The directional interaction that supports molecular recognition, protein folding, and protein structure is provided by hydrogen bonds.⁸⁸ Most protein structural cores are made up of secondary structures, such as the alpha-helix and beta-sheet. This satisfies the potential of hydrogen bonding between the carbonyl oxygen of the main chain and the buried nitrogen in the protein hydrophobic core.⁸⁸ The results of the molecular dynamics simulation collectively suggested that the novel mutant strain of the pathogen in order to resist the antifungal activity of fluconazole, mutated the LEU–370 and PRO–375 to ultimately alter the hydrogen bonding pattern, and stabilize CYP51 at the expense of the bound prosthetic heme (Figures 7 and 8; Supplementary Figure 7).

Environment could play a major role in the biological evolution of organisms in order to adapt and survive. This has led to highly divergent variations in genetic sequences within the same gene of the same species of organisms. Some of these genetic variations are associated with drug resistance. ERG11 gene, being a gene encoding a major protein (caCYP51) essential for the survival of Candida, we were interested in analyzing its evolutionary divergence in isolates that have been implicated in resistance. It will be important to understand the relatedness of the isolates from our environment to other isolates obtained in other regions. We included the criteria for the selection in the materials and methods section. The result of our unrooted phylogenetic analysis between our isolates of interest and other pre-selected isolates based on the ERG11 gene sequence revealed that the most resistant isolate Can–iso–001, is most closely related to the susceptible isolate from dogs (Can–iso–029). This observation, though yet to be empirically proven, suggests that Candida strains could transfer between humans and animals via surfaces or other means of contact and assume different behaviors in different environments. We further showed that the most related azole-resistant isolate to our isolates is HM194202.1 (SZ169) isolated in China, and the least related isolate is AF153849.1 (J913004/1) isolated in Belgium with different mutation spots in their ERG11 sequence. This strongly supports the knowledge that other parts of the caCYP51 protein can determine the susceptibility to azoles, and that our variants of the resistant isolates are new. These findings add to the knowledge that multiple regions of the CYP51 protein either singly, or in combination can determine the susceptibility of Candida strains to azoles.

Conclusion

This study attempted to characterize the species distribution, azole susceptibility, ERG11 gene, and caCYP51 protein profiles of clinical isolates of Candida from the vagina of women and dogs of reproductive age. We focused more on exploring the patterns in sequence variations in the ERG11 gene, and the consequences on the caCYP51 protein, that may determine the susceptibility and resistance of Candida strain to azoles. We achieved this by applying proven classical susceptibility testing methods and molecular sequence analysis techniques to the clinical isolates studied. The result further established that the novel mutant strain mutated the LEU–370 and PRO–375 to stabilize caCYP51 at the expense of the bound prosthetic heme in order to resist the antibiotic activity of fluconazole. Taken together, our result showed new mutations in the heme-binding pocket of caCYP51 that explains the resistance to fluconazole, exhibited by the Candida isolates.

Limitations of the Study

This study provided robust data to support our findings however, a few limitations are noted. First, the differential method for identifying the specific species of Candida is narrow, although highly efficient. Hence, the unidentified NAC species could not be sufficiently characterized. Secondly, other metabolic mechanisms, such as drug efflux not investigated here may also be contributory to the observed resistance to fluconazole however, this study focused on the mechanism of action of fluconazole at the drug target site. Additionally, reverse genetics targeting these points of interest is yet to be applied in order to rule out all other metabolic and genetic possibilities that could be contributory to the resistance.