Analysis of antifungal resistance genes in Candida albicans and Candida glabrata using next generation sequencing

Introduction/Objectives An increase in antifungal resistant Candida strains has been reported in recent years. The aim of this study was to detect mutations in resistance genes of azole-resistant, echinocandin-resistant or multi-resistant strains using next generation sequencing technology, which allows the analysis of multiple resistance mechanisms in a high throughput setting. Methods Forty clinical Candida isolates (16 C. albicans and 24 C. glabrata strains) with MICs for azoles and echinocandins above the clinical EUCAST breakpoint were examined. The genes ERG11, ERG3, TAC1 and GSC1 (FKS1) in C. albicans, as well as ERG11, CgPDR1, FKS1 and FKS2 in C. glabrata were sequenced. Results Fifty-four different missense mutations were identified, 13 of which have not been reported before. All nine echinocandin-resistant Candida isolates showed mutations in the hot spot (HS) regions of FKS1, FKS2 or GSC1. In ERG3 two homozygous premature stop codons were identified in two highly azole-resistant and moderately echinocandin-resistant C. albicans strains. Seven point mutations in ERG11 were determined in azole-resistant C. albicans whereas in azole-resistant C. glabrata, no ERG11 mutations were detected. In 10 out of 13 azole-resistant C. glabrata, 12 different potential gain-of-function mutations in the transcription factor CgPDR1 were verified, which are associated with an overexpression of the efflux pumps CDR1/2. Conclusion This study showed that next generation sequencing allows the thorough investigation of a large number of isolates more cost efficient and faster than conventional Sanger sequencing. Targeting different resistance genes and a large sample size of highly resistant strains allows a better determination of the relevance of the different mutations, and to differentiate between causal mutations and polymorphisms.

Introduction study a set of genes that are known to be involved in antifungal-resistance in a comprehensive strain set. Therefore, NGS based on targeted resequencing design was used to investigate the underlying resistance mechanisms in our clinical isolates.
We sequenced the genes involved in echinocandin resistance (GSC1 in C. albicans and FKS1 and FKS2 in C. glabrata) and the genes involved in azole resistance (ERG11, TAC1 and ERG3 in C. albicans and ERG11 and CgPDR1 in C. glabrata). The aim of this study was not only to detect established but also to identify novel point mutations associated with echinocandin or azole-resistance in the described resistance genes.

Sampling and antifungal susceptibility testing
Forty isolates obtained from specimens such as swabs (4), sterile fluids (8), blood cultures (12), a central venous catheter (1), as well as urine (3), feces (4), sputum (1) as well as not specified (7) from various centres in Austria and Germany were investigated. In addition, the susceptible strains ATCC 90030 and ATCC Y33.90 for C. glabrata, and ATCC 90028, ATCC 10231, as well as the azole-resistant ATCC 64124 for C. albicans were included as controls for the validation of the sequencing process. Antifungal susceptibility testing using the broth microdilution method was performed for all strains and C. parapsilosis ATCC 22019 and C. krusei ATCC 6258 as control strains as described by the European Committee of Antimicrobial Susceptibility Testing (EUCAST E.DEF 7.3 December 2015) [15]. Minimal inhibitory concentrations (MIC) were determined for anidulafungin, micafungin and caspofungin, as well as for fluconazole, posaconazole, voriconazole, itraconazole and isavuconazole. Strains with a MIC one to twofold dilutions above the clinical breakpoint for echinocandins were classified as borderline echinocandin resistant. Multi-resistant isolates were defined as resistant to all tested echinocandins and azoles.

DNA extraction
Due to better quantitative and qualitative results in comparison to commercial DNA extraction kits, a modified SDS CTAB chlorophorm based method was used [16]. DNA extraction was performed from a 24h Candida culture on Sabouraud-dextrose agar (SAB). Mechanical lysis was carried out using 1mm silica spheres under addition of the detergents SDS (sodium dodecyl sulfate) and CTAB (cetyltrimethylammonium bromide) as well as proteinase K (Qiagen, Venlo, Netherlands). After adding chloroform-isoamylalcohol 24:1, the water-soluble polar layer was transferred to a new tube followed by precipitation with ammonium acetate and isopropanol and was subsequently washed using ethanol. The air-dried DNA was then resuspended in 10 mM Tris-EDTA buffer. After extraction, the amount of DNA was determined with Qubit 2.0 via the dsDNA HS kit (Life technologies, Carlsbad, California) and NanoDrop 2000c spectrophotometer (Thermo Scientific, Waltham, Massachusetts). Additionally, the ratios A260/280 and A260/230 were used to estimate the purity of the DNA. additional overhang sequence that is mandatory for sequencing with Illumina technology. Four of the 26 primer pairs published by Garnaud et al. were newly designed by Primer3 Tool because of the formation of hairpins and primer dimers (S1 and S2 Tables) [18,19]. The library amplification was performed using the KAPA HiFi Hot Start Ready Mix Kit (Kapa Biosystems, Wilmington, Massachusetts), a high fidelity polymerase with proofreading activity which is well suited for the production of NGS-libraries [20]. 12.5 ng genomic DNA was added to the PCR-mix. Afterwards, the amplified PCR products of each isolate were pooled. The washing steps were based on Ampure Beads (Beckman Coulter, Brea, California). Subsequently, an index PCR was performed to tag the amplicons for identification of the different isolates after pooling. The DNA quantification of the PCR products was carried out using Qubit 2.0 via the dsDNA HS kit (Life technologies, Carlsbad, California). DNA was diluted to a concentration of 8pM for sequencing on the V2-Flowcell 2x250bp (Illumina, San Diego, California) and all isolates were pooled. The DNA library was denatured according to the protocol [17]. For quality control, the library was spiked with 5% PhiX DNA.

Bioinformatic analysis
The quality of the NGS run was verified using the software FASTQC 0.11.4 [21]. The removal of low-quality bases was carried out with the Trimmomatic-0.35 software [22]. This tool also removed all reads under a minimum length of 90bp. In addition, the first 24bp were removed to exclude the primer sequences. The reads were assembled with Bowtie2-2.2.7 [23]. Subsequently, alignment to the reference sequence was carried out. The strains SC5314 for C. albicans and CBS138 for C. glabrata were used as reference sequences. The gene sequences were downloaded from www.candidagenome.org [24]. To determine variants from the reference sequence Samtools 0.1.19 and VarScan.v2.3.9 [25] were used. After this, SnpEff 4.270 was used to detect alterations causing amino acid substitutions. Finally, a visual validation of the mutations in the assembly files was performed to exclude bias variants. Table 2 shows the sources of the reference sequences. The sequences of the isolates have been deposited in the BioProject database under accession number PRJNA510782.

EUCAST microdilution
Among the 19 C. albicans isolates, two were susceptible control strains, seven were resistant to azoles, six were echinocandin-resistant, two borderline echinocandin-resistant, and two were

Validation of the sequencing process
For the validation of the sequencing process, different control strains were sequenced and no discrepancies to the published sequences were found. The Q30 quality score was 87% and the total number of reads was 11.7 million. The number of reads of each isolate was between 250,000 and 400,000. The minimal and maximal mean coverages of the genes were 2,250x and 32,000x respectively. NGS data revealed the presence of heterozygous mutations in C. albicans with a variant (frequency rate) in approximately 50% of the reads (48.2%-51.7%). For homozygous mutations close to 100% (98.6%-100%) of the reads showed the presence of the variant. As C. glabrata is haploid, variants in the reads of this species had a frequency of approximately 100% (97.9%-100%). In isolate Cg28, we found two mutations with a frequency of 26% and 73%, which could be caused by a subpopulation. This observation was verified with Sanger sequencing. The frequencies of the mutations of all other isolates indicate that they are clonal and without subpopulations.

Non causal polymorphisms
Mutations present in both susceptible and resistant isolates were defined as polymorphisms non-causal for antifungal resistance development. Of the detected 54 missense mutations, seven were already defined as polymorphisms and four mutations were displayed by susceptible isolates and thus were classified as polymorphisms. The missense mutations S935L [18] and S941P [18] in TAC1 have already been described as polymorphisms. The heterozygous mutation S937L was also present in two azole-susceptible isolates. This mutation has not been described, and causality in homozygous cases cannot be ruled out. In heterozygous cases, however, this mutation does not cause azole resistance. In ERG11, the mutations D116E [18], K128T [26] and E266D [26] were observed homozygous as well as heterozygous in susceptible strains and have already been described as non-causal mutations. The mutation V488I in ERG11 was described as causal by Manastir et al., and as non-causal by Wang et al. [26,27]. We found this mutation homozygous in an azole susceptible strain, which supports the study published by Wang et al. In ERG3 we detected the hitherto unknown homozygous mutations A138T, P181A and the heterozygous P267S, which could also be found in susceptible isolates in our study.

Potentially causal resistance mutations
In all strains, 87 different silent mutations compared to the reference sequences of the strains CBS138 and SC5314 were found. Fifty different missense mutations as well as an in frame deletion and two premature stop codons were detected. All mutations including silent mutations are listed in the supplementary data (S3 Table). We identified 43 missense mutations as potentially causal, as they were only present in resistant isolates [28]. Of these, 30 have already been described as causative for resistance acquisition. In our study, 13 potentially causal mutations are reported for the first time. Tables 3, 4, 5 and 6 show the potentially causal mutations of the respective genes. Only missense mutations, which are classified as potentially causal, are listed. In C. albicans, all seven azole-resistant strains showed a mutation in the target gene ERG11. In ERG3 and TAC1 four and five potentially causal mutations could be found respectively. In comparison, no mutations in ERG11 could be detected in C. glabrata, but 10 of 13 azole-resistant strains showed a mutation in CgPDR1. GSC1 or FKS mutations were found in all echinocandin-resistant isolates. No FKS mutations could be found in the four borderline echinocandin-resistant C. albicans and C. glabrata isolates. All eight multi-resistant isolates had at least one mutation possibly leading to azole or echinocandin resistance (Table 5). Analysis of antifungal resistance genes Table 5. Potentially causal missense mutations in C. albicans isolates, homozygous mutations in bold.
In contrast to C. albicans, none of the azole-resistant C. glabrata isolates showed any mutations in ERG11. In CgPDR1, 10 different mutations were observed. H288D has not yet been described. The mutations L291P, G346D, G346A, G348S, G348D and Y372H have been described in these positions but with different amino acid substitutions [33,34]. The mutations D876Y, G1079R and G1079V have been described as causal by Ferrari et al. [33]. No mutations in the sequenced regions of CgPDR1 have been detected in the isolates Cg27, Cg32 and Cg39 (Table 6).

Mutations in echinocandin-resistant C. albicans and C. glabrata
In all six echinocandin-resistant C. albicans isolates, a mutation was found in the target gene GSC1. All mutations were detected in GSC1 HS 1 or its immediate vicinity. The mutations observed were S645P [18], which was detected three times, F641C, F641S [7,35], P649H [36][37][38] and M696V. For the isolates Ca13 and Ca16 the mutations F641C and S645P, respectively, are shown as heterozygous. Isolate Ca14 shows the mutation F641S, which was homozygous in this strain and was associated with higher MIC than for isolate Ca13. Isolate Ca18 showed the two homozygous mutations P649H and M696V, which are located in HS 1 and 90 nucleotides downstream of HS 1, respectively. The borderline echinocandin-resistant isolates Ca19 and Ca21, which exhibit a MIC of micafungin of 0.064 mg/l and 0.032 mg/l respectively, did not have any mutations in GSC1. Our strains showed no missense mutations in the HS regions of GSC1 genes whenever the MIC was lower than 0.125 mg/l in anidulafungin and caspofungin, and lower than 0.064 mg/l in micafungin (Table 5).

Mutations in multi-resistant C. albicans and C. glabrata
The multi-resistant C. albicans isolates Ca12 and Ca22 showed the homozygous premature stop codons Y325 � and Y190 � in ERG3, respectively. In both GSC1 HS regions, no mutations were detected despite elevated echinocandin MIC values (Table 5).
In the six multi-resistant C. glabrata isolates, five mutations were verified in FKS genes. In FKS1, the mutations K1323E and F625S [18] could be detected in isolates Cg29 and Cg46, respectively. In FKS2, the mutations F659S and S663P could be detected in isolates Cg50 and Cg47 respectively. For three of the six multi-resistant isolates, potentially causal mutations were found in CgPDR1: W297R [33,34] for isolate Cg46, V329F [44] for isolate Cg49 and G1088E for isolate Cg50. Only isolates Cg46 and Cg50 showed a mutation in FKS genes and CgPDR1. Isolate Cg49, which displayed moderately elevated MIC values for echinocandins showed no mutations in FKS1 and FKS2 (Table 6).

Discussion
NGS has previously been shown to successfully detect antifungal resistance mutations in clinically important Candida species [18]. With targeted resequencing resistance genes of a high number of isolates can be studied simultaneously. Compared to whole genome sequencing (WGS), this approach reduces sequencing costs, generates more manageable raw data and reduces the burden of data analysis. Investigating 50 isolates the costs would reach 16.000 € using WGS. In this study, the costs were only 3000 €. Thus, the sequencing costs could be reduced to at least 5x due to the targeted resequencing design. Additionally, targeted resequencing results in coverage levels much higher than those achieved with WGS. Thus, this method is very reliable and allows the detection of low frequent variants, e.g. resistant subpopulations. Furthermore, sequencing runs are much faster than with the conventional Sanger method. The sequencing process in our project took about 22 hours for 51 strains with 13 amplicons each. In comparison, using Sanger sequencing on a 4-capillary sequencer, the run time would have been about four weeks. Sanger sequencing is more expensive, and the analysis of this extensive sequencing data would be time-consuming.
In this study, a high number of phenotypically resistant isolates obtained from various centres in Austria and Germany were investigated in order to find mutations causing resistance, and to examine if strains categorized as resistant using the clinical breakpoints possess already described or hitherto unknown resistance mutations. Due to high throughput sequencing, we were able to detect established as well as novel resistance mutations in a large sample of antifungal-resistant strains.

Azole resistance
Target mutations in ERG11 were detected only in C. albicans. In every azole-resistant isolate, potentially causal mutations were detected in this gene. Five of these mutations have already been described as inducing resistance, based on the putative mechanism that the change in the protein sequence leads to a reduced binding affinity of azoles [27,29,30]. In contrast to C. albicans, there were no mutations in ERG11 in any of the 16 azole-resistant C. glabrata strains. The absence of ERG11 mutations in azole-resistant C. glabrata has already been described [45,46]. This suggests that our results are in concordance with other investigations and ERG11 mutations have no impact on azole resistance in our C. glabrata isolates.
CgPDR1 appears to be a more important cause for azole resistance in C. glabrata. CgPDR1 is a transcription factor, which induces the gene expression of the efflux pumps CgCDR1/2p and CgSNQ2p. In CgPDR1 67 gain-of-function mutations have been described hitherto [33]. These mutations were associated with intrinsically high expression of the efflux pumps and specifically were related to azole resistance [12]. In our study, 13 potentially causal mutations were found. In 10 out of 13 azole-resistant strains, at least one mutation was found in CgPDR1. Eleven of these mutations are known, the remaining mutations H288D and G1088E have been detected for the first time. Tsai et al. described the domains of CgPDR1 based on the homology between S. cerevisiae PDR1 and C. glabrata CgPDR1 [34]. The DNA-binding domain is positioned at residues 26-59, the regulatory domain at 322-465 and the activation domain at 903-1107. They found four mutations at residues 280-391. Nine of 13 mutations in CgPDR1 found in our study were located at residues 288-372, near the putative regulatory domain. Accordingly, these mutations could be associated with altered expression of efflux pumps. Three mutations (G1079R, G1079V and G1088E) were identified in the putative activation domain.
Thus, mutations in these regions could be associated with overexpression of drug efflux pumps. To confirm the impact of these mutations, the analysis of the expression levels of efflux pumps could be performed in future studies.
Gain-of-function mutations in the transcription factor TAC1 lead to inherently increased expression of efflux pumps CDR1p and CDR2p and were associated with azole resistance [47]. In our study, five potentially causal mutations in the sequenced areas of TAC1 could be identified. Of these, G649D and F974L have not been described so far. However, since none of the sequenced strains showed mutations in TAC1 only, the relevance of this mechanism could not be clarified in this study.
We could also detect loss-of-function mutations in ERG3 which presumably leads to azole resistance in combination with moderate echinocandin resistance and is described in the multi resistance section.

Echinocandin resistance
In our study, mutations in GSC1, FKS1 and FKS2 were detected in all isolates showing an increase of MIC values (>2 dilution folds above CB), which was associated with complete cross-resistance with the exception of one C. glabrata isolate. Several missense mutations were found in these isolates. For the isolates with moderately elevated MICs, i.e. one to twofold dilutions above the clinical breakpoint and classified as borderline resistant, neither cross-resistance within the echinocandins nor FKS mutations were found.
In our study, a single mutation in the HS regions of the target genes led to complete crossresistance in the class of echinocandins and was seen for both heterozygous and homozygous mutations. In contrast to this observation, the C. glabrata isolate Cg41 displayed an isolated susceptibility to micafungin (MIC 0.016 mg/l) whereas anidulafungin showed an elevated MIC of 0.25 mg/l and caspofungin 2 mg/l. In this isolate the already known mutation P667T in FKS2 HS 1 was detected. The absence of cross-resistance has already been reported [3,48]. This is of particular interest as anidulafungin has been discussed to serve as a surrogate marker for all echinocandins [49,50]. In our case, however, micafungin could have been an important therapeutic alternative. These results indicate that changes in conformation of 1,3-β-D-glucan synthase may lead to an incomplete cross-resistance in the class of echinocandins depending on the specific structure of the respective agent. Therefore, susceptibility testing of every echinocandin seems to be preferable to using anidulafungin as an indicator for echinocandin resistance.
The external localization of some regions of the putative transmembrane protein 1,3-β-Dglucan synthase seems to reflect the HS regions [51]. All of our six echinocandin-resistant C. albicans isolates displayed a mutation in GSC1 HS 1. No mutations were detected in HS 2. Therefore, HS 1 appears to play a more important role in the development of echinocandin resistance in our C. albicans strains. The C. albicans isolate Ca18 showed two homozygous mutations in GSC1. P649H is located at the downstream end of HS 1 and M696V is located 90 nucleotides downstream HS 1 and 30 nucleotides upstream HS 3, which is described by Johnson et al. [51]. Thus, it could be assumed that in rare cases other HS regions may be involved in acquisition of echinocandin resistance.
Out of six echinocandin-resistant C. glabrata strains, four showed a mutation in FKS2 HS 1. The multi-resistant isolates Cg29 and Cg48 showed the mutation K1323E in FKS1, which is five amino acids upstream HS 2. These isolates displayed a minimal rise of MIC values being only onefold dilution above the clinical breakpoint. In these cases, the mutation outside the HS was associated with minimally elevated MIC values.
In isolate Cg51 the in frame deletion F659del in FKS2 HS 1 was detected. This strain showed MIC values in the resistant range for anidulafungin (2 mg/l), caspofungin (> 16 mg/l), and micafungin (4 mg/l). This mutation has already been described associated with very high as well as low MIC values, which might be caused by different expression rates [41][42][43].

Multi-resistance
Out of 40 clinical strains, two C. albicans and six C. glabrata were multi-resistant. In C. glabrata strains, only two isolates showed mutations in both FKS1 and CgpDR1, which explains azole and echinocandin resistance. In the other four isolates, only a mutation in either the FKS genes or in CgPDR1 could be detected.
The multi-resistant C. albicans strains Ca12 and Ca22 only showed a loss-of-function mutation in ERG3 due to the premature stop codons Y325 � and Y190 � , which presumably leads to azole resistance and moderate resistance to echinocandins without displaying FKS mutations. Although the molecular mechanism for the moderate echinocandin resistance is not clear, this observation is in concordance with Rybak et al. [52]. Thus, the importance of ERG3 loss-offunction mutations for resistance development should not be overlooked.
The fact that a potentially causal mutation for azole and echinocandin resistance is not present in every resistant isolate suggests that different-hitherto unknown-resistance mechanisms are involved.

Conclusion
NGS proved to be a suitable method to detect resistance mutations. This technique allows a thoroughly, more cost-efficient and much faster sequencing method than conventional Sanger sequencing. Furthermore, the targeted resequencing design enables the investigation of a larger sample size than WGS. In combination with the phenotypic analysis of resistance patterns, conclusions can be drawn about the underlying molecular mechanisms.
We investigated 40 resistant clinical C. albicans and C. glabrata isolates and found 30 described and 13 novel mutations in six resistance genes. In addition, a high rate of polymorphisms was found in the coding sequences. This observation underlines the importance of the differentiation between polymorphisms and causal mutations. This applies especially for azole resistance, where several mechanisms can lead to resistance and interact with each other. As a consequence, a SNP database, which includes each variant as well as the phenotype, would be helpful to distinguish between polymorphisms and relevant mutations.
In conclusion, an association between mutations in FKS genes and echinocandin resistance can be confirmed. However, the acquisition of azole resistance seems to have multifactorial causes. The mutations in ERG11 appear to play a role only in C. albicans. In C. glabrata, overexpression of efflux pumps is conceivable instead. In C. albicans, homozygous ERG3 nonsense mutations seem to be associated with azole resistance and moderately elevated echinocandin MICs. Four of the multi-resistant C. glabrata isolates showed no underlying mutations for both echinocandin and azole resistance. Hence, there are still other cellular mechanisms, which require further investigations.
Supporting information S1