A Comparative Transcriptome Between Anti-drug Sensitive and Resistant Candida auris in China

Candida auris emerged as a pathogenic species of fungus that causes severe and invasive outbreaks worldwide. The fungus exhibits high intrinsic resistance rates to various first-line antifungals, and the underlying molecular mechanism responsible for its multidrug resistance is still unclear. In this study, a transcriptomic analysis was performed between two C. auris isolates that exhibited different anti-drug patterns by RNA-sequencing, namely, CX1 (anti-drug sensitive) and CX2 (resistant). Transcriptomic analysis results revealed 541 upregulated and 453 downregulated genes in the resistant C. auris strain compared with the susceptible strain. In addition, our findings highlight the presence of potential differentially expressed genes (DEGs), which may play a role in drug resistance, including genes involved in ergosterol and efflux pump biosynthesis such as SNQ2, CDR4, ARB1, MDR1, MRR1, and ERG genes. We also found that Hsp related genes were upregulated for expression in the anti-drug-resistant strain. Biofilm formation and growth conditions were also compared between the two isolates. Our study provides novel clues for future studies in terms of understanding multidrug resistance mechanisms of C. auris strains.


INTRODUCTION
Candida auris is a species of fungus that was firstly isolated from a patient's external ear canal in Tokyo, Japan and represented a novel species within this genus in 2009 (Satoh et al., 2009). However, a study performed in 2011 revealed that C. auris already existed as a nosocomial infective agent in South Korea since 1996, in a case that the infection was treated as unidentified yeasts and which eventually caused nosocomial fungemia in a 1-year-old girl (Lee et al., 2011). Since its first description in 2009, the presence of C. auris was soon reported in all continents, except Antarctica, and which caused infections in more than 40 countries Rhodes and Fisher, 2019;Vila et al., 2020). Genomic epidemiology revealed that C. auris exhibits five genetically diverse clades by genome sequencing. Consequently, this species was divided according to the following geographical clusters: Clade I (South Asian), Clade II (East Asian), Clade III (African), Clade IV (South American), and Clade V (Iran). Clade V was isolated from a patient in Iran and is the most recently identified clade Chow et al., 2019). Genetic sequences of isolates from patients or environments are very important to trace the transmission of C. auris among different countries, hospitals, and perhaps different patients' parts. What makes this species so important is its ability to spread and cause nosocomial outbreaks in hospitals and healthcare facilities in a similar manner to bacteria (Schelenz et al., 2016;Adams et al., 2018;Eyre et al., 2018;Ruiz-Gaitan et al., 2018Armstrong et al., 2019;O'Connor et al., 2019). In fact, a recent review showed that the majority of the recovery sites in C. auris outbreaks involved candidemia (blood) and skin, which in turn underlines the ability of C. auris to induce invasive infections and colonization (Vila et al., 2020). Misidentification of C. auris, which unavoidably induces inappropriate and delayed treatment, is another important reason for these outbreaks. Unlike other yeasts, C. auris has often been confused with other pathogens such as C. parapsilosis, C. haemulonii, and C. sake by commercial systems (Lee et al., 2011;Kathuria et al., 2015;Mizusawa et al., 2017;Adams et al., 2018;Snayd et al., 2018;ElBaradei, 2020). Reliable identification methods involve systems with updated databases and molecular methods, which target the specific sequences of C. auris (Kathuria et al., 2015;Kordalewska et al., 2017;Sexton et al., 2018;Forsberg et al., 2019;Lima et al., 2019;Lone and Ahmad, 2019). However, the most worrying characteristic of C. auris pertains to its intrinsic or required resistance to a variety of first-line antifungals commonly used in clinical settings, including the three main categories of antifungals, i.e., azoles, echinocandins, and polyenes, thus limiting the underlying treatment options (Adams et al., 2018;Chaabane et al., 2019;Wickes, 2020). As a result, the misidentification and multidrug resistance of C. auris have both been prominent factors that enhance the ability of this fungus to cause infections with significant patient mortality, especially in patients that already suffer from concurrent diseases or receive invasive treatments Armstrong et al., 2019;de Jong and Hagen, 2019;Taori et al., 2019).
Despite the fact that the multidrug resistance of C. auris has been a prevalent concern, there have only been a limited number of research studies performed to expound its mechanisms of antifungal resistance, and thus, these molecular mechanisms remain unknown (Chaabane et al., 2019;Kean and Ramage, 2019). For instance, several studies have focused on its phenotype. However, the majority of studies on antifungal resistance mechanisms are based on drug resistance-related genes in non-auris Candida species that have been previously reported (Arastehfar et al., 2020;Lee et al., 2020). Pertaining to azoles resistance, mutations in ERG11 and ERG11 overexpression play a distinct role in C. auris that facilitate drug target alteration and overexpression, respectively (Pristov and Ghannoum, 2019). Another common mechanism for azole resistance in C. auris involves efflux pumps enhanced overexpression including Major Facilitator Superfamily (MFS) and ATP Binding Cassette (ABC) that accelerate drug efflux (Kean et al., 2018;Lee et al., 2020). Echinocandin resistance in C. auris is the result of mutations in FKS1 due to drug target alterations (Pristov and Ghannoum, 2019). However, the resistance mechanism to amphotericin B is still not confirmed, and it cannot be explained by any of the proposed mechanisms thus far. Apart from the pathways that the above genes participate in, certain genes play a role in triggering stress responses such as HSP90, which also confer resistance (Lee et al., 2020).
In this study, we performed transcriptome analysis on two C. auris strains using RNA-sequencing (RNA-seq). One strain showed elevated minimum inhibitory concentration (MIC) in fluconazole and micafungin, whereas the other strain was a susceptible strain. Until now, a limited number of research studies have been performed with RNA-seq between susceptible and resistant C. auris strains without imposing any conditions in China. Consequently, through the analysis of gene expression differences, we aimed to explore the resistance mechanism and identify genes that may play an important part of this process.

Bacterial Strains and Identification
The first C. auris strain used in this study was isolated from the environment in China and was named CX1. The other C. auris strain was acquired from the NCCLs (National Center for Clinical Laboratories), was isolated from a patient's ascitic fluid, and was named CX2. Two isolates were identified as C. auris by sequencing ribosomal DNA internal transcribed spacer (ITS) combined with using matrix-assisted laser desorption/ionization time of flight mass spectrometry (Bruker Daltonics, Germany) (Schoch et al., 2012;Fraser et al., 2016). The Ethics Committee of The First Affiliated Hospital of Nanchang University (approval no. 2016026) approved the present study. All participants provided a written informed consent to participate in this study.

Antifungal Susceptibility Testing
In vitro antifungal susceptibility testing was performed with YeastOne plate (Thermo Fisher, United States) on three replicates of two strains according to the Clinical and Laboratory Standards Institute (CLSI) broth microdilution method M27-A3 (Clinical and Laboratory Standards Institute [CLSI], 2008). Several fully isolated strains were also selected from a yeast isolate of 24 h pure yeast culture species, emulsified in sterile water, and mixed with vortex. Then, we adjusted yeast suspension to 0.5 McFarland and Pipetted 20 µl to 11 ml fungus inoculation broth. Broth suspension (100 µl) was pipetted onto the drug sensitive plate and incubated in 35 • C for 24 h. Finally, SensititreVizion system (Thermo Fisher, United States) was used to read the test results of the drug sensitive plate. The activities of nine antifungals against two C. auris isolates were tested, including fluconazole, itraconazole, voriconazole, posaconazole, micafungin, anidulafungin, caspofungin, amphotericin B, and 5flucytosine. The MIC endpoints were interpreted in tentative breakpoints proposed by CDC 1 , as follows: ≥32 for fluconazole, ≥2 for amphotericin B, ≥4 for anidulafungin and micafungin, and ≥2 for caspofungin. 30S • C for 18 h. The fungus was collected at an approximate OD 600 = 0.7, transferred to the EP tube, and resuspended in 50 µl of sterile water preheated at 30 • C. Cells were cooled quickly in the liquid nitrogen and grinded to powder in the pre-cooled grinding tool. Then, we added 1 ml of Trizol solution (Invitrogen, United Kingdom), grinded, sealed the mortar with tin foil, and let it stand. When the Trizol-bacteria mixture became liquid, we gently grinded this mixture again. Ribozyme-free pipette tips were used to suck the Trizo-bacterial mixture into new ribozymefree 1.5 ml EP tubes, which were subsequently centrifuged at 12,000 rpm for 10 min at 4 • C. We pipetted the water phase to new 1.5 ml EP tubes and added an equal volume of 25:24:1 phenol/chloroform/isoamyl alcohol, which was vigorously shook for 10 s and centrifuged at 4 • C, 12,000 rpm for 5 min. Water phase was transferred to new 1.5 ml EP tubes, and an equal volume of isopropanol was added into an ice-bath for 10 min. We then centrifuged the tubes at 12,000 rpm at 4 • C for 5 min, and we discarded the supernatant. The total RNA was then washed with 75% ethanol and centrifuged at 12,000 rpm for 5 min, and the supernatant was discarded. Total RNA was dried for 5-10 min, resuspended in 25 µl DEPC water, and temporarily stored at −20 • C. Following this protocol for obtaining the total RNA, we used electrophoresis to observe the integrity of the RNA using 1.0% agarose gel. Quality and concentration of the isolated RNA were assessed by NanoDrop 2000c (Thermo Scientific, United Kingdom). CX1 and CX2 isolates were cultured in triplicates named CX1-1, CX1-2, CX1-3 and CX2-1, CX2-2, CX2-3, respectively.

Transcriptome Analysis
In total, six samples and three biological replicates for two strains were sent for cDNA library construction, transcriptome sequencing, and analysis conducted by OE biotech Co., Ltd. (Shanghai, China). Furthermore, TruSeq Stranded mRNA LT Sample Prep Kit (Illumina, San Diego, CA, United States) was used to construct cDNA library in accordance with the manufacturer's instructions. After the constructed library was qualified with Agilent 2100 Bioanalyzer, it was sequenced using Illumina HiSeq X Tento to generate 150 bp pairedend reads. Raw reads of fast format were preprocessed using Trimmomatic (Bolger et al., 2014), and the number of reads in the whole quality control was statistically summarized. Each sample's clean reads remained after quality pretreatment steps, including the removal of reads containing adaptor, low quality reads, and bases arranged in different ways from the 3 end and 5 end. HISAT2 (Kim et al., 2015) was used to compare the clean reads with the specific reference genome 2 to obtain the position information on the reference genome and the unique sequence feature information of the sequenced samples. The STRING (Search Tool for the Retrieval of Interacting Genes) database (Szklarczyk et al., 2019) was used to construct a protein-protein association network and to visualize the interactome network. 2 http://www.candidagenome.org/download/chromosomal_feature_files/C_auris_ B11221/

DEG Analysis
Each gene's FPKM (Roberts et al., 2011) value was calculated using Cufflinks (Trapnell et al., 2010), and its read counts were obtained by HTSeq-count (Anders et al., 2015). Cluster analysis was performed with the "pheatmap" package to explore gene expression pattern. Correlation test between samples was carried out using R language to calculate the Pearson correlation coefficient. DEseq (2012) R package was used to perform differential expression analysis. The threshold of significantly differential expression between samples was a false discovery rate (FDR) < 0.05 and fold change > 1.5.

Functional Annotation
Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Gene and Genomes (KEGG) (Kanehisa et al., 2008) pathway enrichment analysis of differentially expressed genes (DEGs) were performed using R based on the hypergeometric distribution. Biological process, cellular component, and molecular function enrichment analysis in GO level 2 were performed on DEGs using the fisher algorithm.

Virulence Factor Prediction
We initially searched for genes associated with phospholipase, proteinase, hemolysin, adhesin, and biofilm formation in our transcriptome results. Subsequently, we searched for phospholipase, proteinase, hemolysin, and adhesin in Candida genome database (CGD). After retrieving related genes in other Candida species, we then searched for orthologous genes or best hits in Candida auris, shown in Supplementary Table 2. For biofilm formation, we downloaded the biofilm formation phenotype from CGD and performed a blastp with DEGs in the transcriptome. The criteria for filtering the blastp results were over 50% coverage and identity.

Biofilm Formation Experiment and Growth Experiment
Prior to developing the biofilm, we treated the 96-well plate (Corning, United Kingdom) with fetal bovine serum, blocked at 4 • C for 72 h, and then washed with sterile phosphate buffer saline. Candida was inoculated in YPD medium broth, incubated with shaking at 220 rpm at 30 • C overnight, and then resuspended with Spider medium to make OD 600 = 0.5. We added 200 µl bacterial solutions to each well in the experimental group and 200 µl Spider medium in the control group. The plate was incubated at 37 • C under shaking at 200 rpm for 90 min and then washed with PBS. Each hole was added with 200 µl Spider medium and sealed with sealing film. After being incubated at 37 • C under shaking at 200 rpm for 48 h, the culture medium was discarded, and each hole was washed with PBS. We fixed each well with 200 µl methanol for 30 min and then used 200 µl 11% crystal violet to stain after discarding fixative. Then, we absorbed the crystal violet and washed each well with PBS after washing with slow water flow and then decolorized each well with glacial acetic acid for 30 min. Finally, biofilm formation was measured according to spectrophotometric methods using microplate reader (Biotek, United States). A spot dilution assay was performed to compare the growth status of the two C. auris isolates and C. albicans. Candida was incubated in YPD liquid medium with constant shaking at 220 rpm at 30 • C overnight. Then, we collected the bacteria by centrifugation at 3,000 rpm, washed with sterile PBS, and resuspended with PBS. The yeast suspension was adjusted to an optical density (OD 600 ) of 0.1 and was diluted by 10-fold serial in sterile PBS to a final OD 600 of 10 −2 , 10 −3 , 10 −4 , 10 −5 , 10 −6 , and 10 −7 . A total of 1 µl suspension of each dilution was spotted on the YPD agar plate, and the plate was cultured at 30 • C for 4 days before observing the growth differences of the two strains' colonies.

Biofilm Formation and Growth Condition
As opposed to the two C. auris isolates that did not form any biofilm and were thus the same as the blank control, C. albicans (SC5314) formed biofilm in the 96-well plate (Supplementary Figure 1A). The spot assay was used to assess the growth status of the two C. auris isolates and C. albicans (SC5314) (Supplementary Figure 1B). No growth differences were found between the two C. auris strains and between C. auris and C. albicans strains.

Transcriptome Analysis
In this study, we performed a comparative RNA-seq analysis on two clinical types (CX1 and CX2) including six samples (CX1-1, CX1-2, CX1-3, CX2-1, CX2-2, and CX2-3) to identify DEGs that may potentially facilitate drug-related resistance. The Illumina sequencing of the transcriptome of these six samples produced raw RNAseq reads. To minimize the effect of data error in our results, we used the Trimmomatic software to pre-process the original data, and the number of reads in the whole quality control process was statistically summarized (Supplementary Table 1).
Hisat2 was used to sequence the aligned Clean Reads with the specified reference genome to obtain the location information on the reference genome or gene, as well as the sequence characteristic information that was unique to the sequencing sample (Supplementary Table 2). The percentage of raw reads of samples mapped to the C. auris genome was high (over 98% in total mapped reads), thus indicating that these samples were consistent with C. auris.
Principal-component analysis and cluster analysis were used to demonstrate a visual representation of the correlation in transcriptome between the different replicates ( Figures 1A,B). Samples from different isolates were clustered separately, whereas samples from the same isolates were clustered together. Furthermore, the protein coding gene expression level was used to generate a heatmap of correlation coefficient between samples ( Figure 1C). The results obtained reflected a high-level of relevance among samples from the same strain and slight differences between CX1 and CX2.
RNA-seq data analysis clearly displayed a wide range of differences between CX1 and CX2 gene expression in the transcriptome. More specifically, genes with a 1.5-fold change (up or down) in the level of expression were considered to reflect differences in gene expression, and they were thus regarded as significantly regulated genes. A negative binomial distribution test revealed a P-value < 0.05. Nine hundred ninety-four statistically significant variation genes were found in expression between our samples, of which 541 were upregulated and 453 were downregulated (Supplementary Table 3). The volcano plot was used to demonstrate the number of significantly DEGs between CX1 Every dot represents a differentially expressed gene. Gray dots were non-significantly different genes, and the green and red dots were significantly upregulated different genes in CX1 and CX2, respectively. and CX2 (Figure 1D). Clustering analysis of the differential expression patterns showed that the DEGs were consistent across replicates, with a significant variation between CX1 and CX2 (Figure 2).

Enrichment Analysis of DEGs
After obtaining the statistically DEGs, the GO enrichment analysis was carried out on the DEGs to describe their respective functions (combined with the GO annotation results). The DEGs could be divided into three main GO categories, i.e., biological process, cellular component, and molecular function (Figure 3). The 669 downregulated differential expression genes were assigned to 43 GO terms, including 21 biological processes, 11 cellular components, and 11 molecular functions. More specifically, these genes were distributed as follows: biological processes included cellular processes (70.3%), single-organism processes (62.9%), and metabolic processes (53.7%); cellular components associated with cells (85.2%), cell parts (84.6%), and organelles (60.5%); and molecular functions containing catalytic activity (49.3%) and binding (49.0%). In contrast, the 885 upregulated differential expression genes were assigned to 41 GO terms, including 20 biological processes, 10 cellular components, and 11 molecular functions. More specifically, these genes were distributed as follows: biological processes, including cellular processes (84.0%) and metabolic processes (74.4%); cellular components associated with cells (93.6%), cell parts (93.6%), and organelles (76.0%); molecular functions contributing to binding (58.8%) and catalytic activity (42.7%).
Moreover, the pathway-analysis of differentially expressed protein coding genes using the KEGG database (combined with KEGG annotation results) could reveal the relationship between drug resistance and cellular pathway changes. A total of 427 downregulated genes and 558 upregulated genes were categorized into known KEGG pathways. Among the 427 downregulated genes, 82 DEGs were distributed in cellular processes that were mainly sub-categorized under transport and catabolism (17.1%), and cell growth and death (10.2%). In addition, 44 DEGs were distributed in environmental information processing that were mainly sub-categorized under signal transduction (17.1%), whereas 41 DEGs were distributed in genetic information processing that were mainly sub-categorized under folding, sorting, and degradation (9.8%). Finally, 260 DEGs were   distributed in metabolism, and they were mainly sub-categorized under carbohydrate metabolism (20.7%), amino acid metabolism (19.9%), and global and overview maps (14.6%). In contrast, among the 558 upregulated genes, 29 DEGs were distributed in cellular processes that were mainly sub-categorized under cell growth and death (3.4%), and 27 DEGs were distributed in environmental information processing that were mainly sub-categorized under signal transduction (6.1%). Furthermore, 256 DEGs were distributed in genetic information processing that were especially sub-categorized under translation (48.7%), and 246 DEGs were distributed in metabolism that were mainly sub-categorized under nucleotide metabolism (9.5%), amino acid metabolism (10.8%), and global and overview maps (10.3%) (Figure 4).

Virulence Factors
A previous study performed showed that C. auris exhibit phospholipase, proteinase, and hemolysin activity in vitro . In another study, C. auris were found to produce phospholipase and proteinase that varied by strain (Larkin et al., 2017). The draft genome of C. auris revealed that virulence may be caused as a result of its diverse transporters, secreted aspartyl proteinases, secreted lipases, phosphatases, phospholipase, mannosyl transferases, integrins, adhesins, and transcription factors (Chatterjee et al., 2015). In this study, we investigated phospholipase, proteinase, hemolysin, adhesin, and biofilm formation-related genes in CGD. Initially, we did not identify hemolysin-related genes in DEGs between sensitive and resistant strains. However, a certain number of DEGs exist, which were related with phospholipase, proteinase, adhesion, and biofilm formation ( Table 2). In fact, Table 2 shows that a greater number of downregulated than upregulated genes were related with phospholipase, proteinase, and adhesin. These findings are consistent with the common knowledge that drug resistance strains possess fewer virulence factors and thus demonstrate wear virulence. Meanwhile, our findings regarding DEGs related with biofilm formation were inconsistent with previous results that suggested a slightly greater number of upregulated genes compared to downregulated genes.

Protein Interaction Network
In order to acquire a more accurate visualization of the molecular mechanism of multidrug resistance involved in C. auris, we used the STRING (Search Tool for the Retrieval of interacting Genes) database to generate a predicted protein interaction network that contained the top 20 downregulated DEGs and top 20 upregulated DEGs (Figure 5). The red and green nodes with gene names represented the upregulated and downregulated DEGs, respectively. The obtained picture shows that the proteins were divided into two clusters. Consequently, the core genes that more predicted associations with other genes and other remaining genes within the network were then investigated.

DISCUSSION
Candida auris is a pathogen that has been known for more than 10 years, yet it continuously causes outbreaks and exhibits alarmingly high drug resistance rates and even pan-drug resistance, which has not been efficiently studied (Wickes, 2020).
Despite the high mortality rates caused by Candida auris infections, due to the rapid spread and frequent worldwide outbreaks caused by this fungus, only a small number of drugs can appropriately treat fungal infections. To understand the mechanism that facilitates drug resistance, we performed transcriptome of two Candida auris isolates, one of which was found to be resistant to fluconazole (MIC = 64) and micafungin (MIC = 8), and the other was found to be susceptible to antifungal drugs based on CDC reports (see text footnote 1). Therefore, we focused on the underlying mechanism of resistance to azoles and echinocandins in Candida auris. Due to the fact that the two Candida auris strains were unable to form biofilms, we searched the agglutinin-like sequence (ALS) genes in DEGs that play an important role in biofilm formation after attachment to abiotic surfaces (Hoyer and Cota, 2016). As expected, we only found ALS4, which is the most frequently expressed gene in the ALS gene family and which is differentially expressed between these two strains (but with low expression levels) (Monroy-Perez et al., 2012). The average value of ALS4 gene expression in three biological replicates of CX1 was 0.012, whereas the average value of ALS4 gene expression in three biological replicates of CX2 was 4.069. The significantly low expression level of ALS genes may be the reason for their inability to form biofilms. Enrichment analysis of DEGs revealed abundant differential pathways and gene functions between strains. The comparison of the GO classification enriched by upregulated and downregulated genes revealed that more upregulated genes were found than downregulated genes involved in electron carrier activity, structure molecule activity, and translation regulator activity, which may be related with transmembrane transporters such as efflux pumps. In addition, more downregulated genes were involved in immune system processes, extra cellular regions, extracellular region parts, and metallochaperone activity; however, further research is required to validate these findings. Nonetheless, pertaining to the KEGG pathway classification, more upregulated genes were found in resistant Candida auris with respect to translation, transduction, and nucleotide metabolism. However, more downregulated genes were involved in transport and catabolism, signal transduction, folding, sorting and degradation, and most metabolism pathways, a finding indicating that the metabolism of resistant strain is less active compared to susceptible strains or that a lower number of genes can in fact play a more significant role. The azoles inhibit the activity of the lanosterol 14-αdemethylase encoded by ERG11, thus blocking ergosterol synthesis. Known mechanisms of resistance to azole include point mutations in ERG11 that lead to a decrease in the affinity of drugs and enzymes. Furthermore, overexpression of ERG11 and efflux pump genes can also cause a decrease in azole susceptibility (Lee et al., 2020). In our results, we investigated genes that are involved in steroid biosynthesis and efflux pumps pathways, which may in turn facilitate resistance to fluconazole. We identified genes that are involved in ergosterol and efflux pump biosynthesis that were upregulated expressed in resistant Candida auris strains including SNQ2, CDR4, ARB1, MDR1, MRR1, and 9 (ERG1, ERG7, ERG11, ERG24, ERG25, ERG6, ERG2, ERG3, and ERG5) of 13 of the ergosterol synthesis genes, which were consistent with previous studies performed Kean et al., 2018;Rybak et al., 2019). SNQ2, CDR4, MDR1, and ARB1 were found to be present and upregulated in drug-resistant isolates of C. auris (Munoz et al., 2018;Wasi et al., 2019). The ATP binding cassette (ABC) superfamily and major facilitator (MF) superfamily are the two main classes of efflux pumps. RNA-seq analysis indicated that the downregulated efflux pump genes such as an ABC transporter gene (ABCF3) and two MFS transporter genes (hxnP and ecdD) may play an atypical role. However, these genes have not been previously reported, and thus, further investigating their properties is imperative. In addition to ERG11 and efflux pump overexpression, we also identified point mutations in ERG11. The amino acid substitutions at Y132F, K143R, and F126L were considered the cause of the azole resistance in C. auris, with Y132F being the most widespread mutation associated with azole resistance (Lee et al., 2021). We found point mutations at T395A which led to the amino acid substitution at Y132F in all three biological replicates of drugresistant strains; this was surprisingly consistent with previous studies. However, we did not identify any point mutations in ERG11 in the susceptible strains.
Echinocandins exert their antifungal effect by inhibiting beta (1,3) D-glucan synthase encoded by FKS, thus causing a defective cell wall. Unlike azole resistance, several studies exist, which demonstrate that amino acid substitutions at S639F, S639P, and S639Y in FKS1 simply lead to elevated MIC levels of echinocandins, as opposed to the role played by efflux pumps (Berkow and Lockhart, 2018;Chowdhary et al., 2018;Kordalewska et al., 2018;Lee et al., 2021). However, our study detected one amino acid substitution at R464P in the sensitive strain, in addition to three point mutations, T588C, A897G, and G2458T (same-sense mutations), in the resistant strain. We tried to explain the high MIC of micafungin (MIC = 8 µg/ml) in CX2 with HSP-associated genes. Heat shock proteins (HSPs) contain a large number of proteins that are distributed widely and are involved in many cellular pathways to modulate stress responses (Gong et al., 2017;Lee et al., 2020). In transcriptome analysis, we found that all HSP related genes were upregulated in resistant Candida auris strains such as HSP90, a strain that has been studied before. However, the exact role of HSPs in the resistance mechanism needs to be further studied.
In addition, this study assessed the potential role of protein modification genes in the underlying resistance mechanism.
Epigenetics has been continuously studied in mammals; however, few studies have been performed that focus on fungal resistance mechanisms and epigenetics (Madhani, 2021). Therefore, we searched for genes that are involved in methylation, SUMOylation, ubiquitination, acetylation, glycosylation, and phosphorylation processes. Our results indicate that the number of DEGs related with methylation was the largest among the protein modification genes, whereas SUMOylation and ubiquitination did not reveal any significant differences. Although multi-omics studies of fungi can quickly obtain a lot of information, more experimental studies are still needed to evaluate the biological effects of these genes.
After searching for genes with relatively large differential expression in the network, we identified that the functions of most genes have not been individually studied in C. auris before. Consequently, we searched their orthologous genes in other Candida species and yeasts in Pubmed, CGD (Candida Genome Database), and SGD (Saccharomyces Genome Database). FCY2 participated in the process of transmembrane transporting and converting 5-fluorocytosine (5FC) into toxic 5-fluorouracil (5FU). Furthermore, FCY2 mutations can mediate resistance to 5FC in both Candida species and Cryptococcus (Billmyre et al., 2020). In this study, although we did not know whether genetic mutations are present, we found that FCY2 (B9J08_002435) was upregulated in the resistant C. auris strain, a finding that is consistent with the drug sensitivity results obtained for 5FC (MIC≤0.06 µg/ml). EBP1 was also upregulated in the resistant C. auris strain, and previous studies identified that it plays an essential role in the yeast-to-hypha transition (Kurakado et al., 2017). Furthermore, pertaining to genes involved in CO 2 signaling in fungi, NCE103 and UME6 were upregulated in the resistant C. auris strain. NCE103 and UME6 not only are essential for ensuring the carbon supply required for cell metabolism but also play an important role in signal transduction process such as fungi's morphology and communication (Martin et al., 2017;Lu et al., 2019). Therefore, we believe that these two genes are perhaps involved in the drug resistance mechanism, yet their specific roles need to be further determined. Among the downregulated genes, we found certain genes or orthologous genes with transmembrane transporter activity, including B9J08_005345, B9J08_004188, B9J08_005571, B9J08_00002660, and B9J08_004099. These transmembrane transporter genes may reduce drug uptake and lead to drug resistance through a different approach than drug efflux. In addition, ADY2, which is assigned to the acetate uptake transporter family, was found to be downregulated in the resistant C. auris strain. So far, no pleiotropic drug resistance (PDR) transporters belonging to the ABC superfamily exist, which are found to be related with the export of carboxylates under acid stress conditions (Alves et al., 2020). Our research was also in line with a previous study that verified that PKC1 downregulation leads to defects in the formation of biofilms (Heinisch and Rodicio, 2018). Moreover, B9J08_004787 (orthologous gene: YPR013C) and B9J08_004804 (orthologous gene: At4g20930) downregulation may lead to reduced virulence, which is consistent with the weakened virulence of the resistant strain (Mao et al., 2008;Otzen et al., 2014).
Through transcriptome analysis between resistant and susceptible Candida auris strains, our study may provide some novel ideas for future studies with respect to understanding the mechanisms of drug resistance in Candida auris. The limitation of this study is that only one drug-resistant C. auris strain was compared with one susceptible C. auris strain. Future research should be performed to confirm the exact function of a certain DEG involved in molecular mechanism of multidrug resistance.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm. nih.gov/, PRJNA735406.

ETHICS STATEMENT
The Ethics Committee of The First Affiliated Hospital of Nanchang University (approval no. 2016026) approved of the present study. The participants provided written informed consent to participate in this study.