Persistence of an epidemic cluster of Rhodotorula mucilaginosa in multiple geographic regions in China and the emergence of a 5-flucytosine resistant clone

ABSTRACT Rhodotorula mucilaginosa, an environmental yeast widely used in industry and agriculture, is also an opportunistic pathogen resistant to multi-antifungals. During the national surveillance in China, R. mucilaginosa has been documented from various hospitals and regions. At present, the molecular epidemiology of invasive infections caused by R. mucilaginosa and their resistance profiles to antifungals were unknown. Here we collected 49 strains from four hospitals located in different geographic regions from 2009 to 2019 in China, determined their genotypes using different molecular markers and quantified susceptibilities to various antifungals. Sequencing of ITS and D1/D2 regions in rDNA indicated that 73.5% (36/49) of clinical strains belong to same sequence type (rDNA type 2). Microsatellite (MT) genotyping with 15 (recently developed) tandem repeat loci identified 5 epidemic MT types, which accounted for 44.9% (22/49) of clinical strains, as well as 27 sporadic MT types. Microsatellite data indicated that the presence of an epidemic cluster including 35 strains (71.4%) repeatedly isolated in four hospitals for eight years. Single nucleotide variants (SNVs) from the whole genome sequence data also supported the clustering of these epidemic strains due to low pairwise distance. In addition, phylogenetic analysis of SNVs from these clinical strains, together with environmental and animal strains showed that the closely related epidemic cluster strains may be opportunistic, zoonotic pathogens. Also, molecular data indicated a possible clonal transmission of pan echinocandins-azoles-5-flucytosine resistant R. mucilaginosa strains in hospital H01. Our study demonstrated that R. mucilaginosa is a multi-drug resistant pathogen with the ability to cause nosocomial infection.


Introduction
Invasive fungal diseases (IFD) are associated with high mortality and morbidity [1]. In the last two decades, reports of invasive infections caused by Rhodotorula mucilaginosa have increased [2]. Although R. mucilaginosa is an environmental yeast occurring in the soil, lakes and deep-sea [3], it has emerged as an opportunistic pathogen in cases of fungemia, central nervous system infections, ocular infections, peritonitis and endocarditis with 9.1%-15% mortality [4,5]. In terms of antifungal susceptibility profiles, echinocandins and azoles including fluconazole, the most widely used antifungal agent, are not recommended for treatment because of high minimum inhibitory concentrations (MICs) against R. mucilaginosa [6]. Therefore, only amphotericin B (AMB) and 5-flucytosine (5-FC) can be used in treatment [7].
In the ARTEMIS surveillance project, Asia surveillance and National China Hospital Invasive Fungal Surveillance Net (CHIF-NET) programs, Rhodotorula species accounted for 4.2-6.0% among invasive noncandidal yeast infections [5,8,9]. While cases of Rhodotorula-related infections from China were higher than cases from other 23 countries, such as Brazil, Spain and India [7]. Pulsed-field gel electrophoresis (PFGE) and random amplification of DNA (RAPD) have been utilized in the genotyping of R. mucilaginosa [10,11]. However, no universal and high-resolution molecular typing method for R. mucilaginosa has been established worldwide, and molecular epidemiological investigations have been rarely performed.
In the CHIF-NET program (2009-2019), we discovered that the number of R. mucilaginosa strains from one hospital (hospital H01) was higher than strains from other hospitals located in seven geographic regions in China. To determine if these strains in hospital H01 belonged to a clonal lineage and were part of a hospital-associated outbreak, we investigated the molecular epidemiology and genetic variability of the R. mucilaginosa strains, including those collected from other three hospitals using microsatellite markers and whole genome sequencing (WGS) approaches. We included four environmental strains and one animal strain to infer the genetic relatedness among these clinical and environmental strains in China. Presently, there was no available breakpoint of antifungal agents for R. mucilaginosa, we investigated the antifungal susceptibility profiles on these Chinese isolates based on M60 document approved in 2020 for interpretive breakpoints among Candida species. Understanding the genetic diversity and antifungal susceptibilities of R. mucilaginosa could help clinicians to implement appropriate diagnostic, therapeutic and preventive strategies.

Strains and cultures
We included 49 strains from 4 hospitals, from which more than 5 Rhodotorula mucilaginosa strains had been isolated between August 2009 and August 2019, as a part of the CHIF-NET program. Basic patients and sampling information were collected, while detailed clinical data were not included in this retrospective investigation. Four hospitals (named H01, H02, H03 and H04) are in the north, east, northeast and southwest parts of China. The inclusion/ exclusion criteria of strains in CHIF-NET program were previously described [12]. In addition, we included four environmental strains (one of them was the type) and one animal strain of R. mucilaginosa to study the genetic relationships among clinical, animal and environmental strains. The type strain (R. mucilaginosa CBS 316) was purchased from Westerdijk Fungal Biodiversity Institute (Utrecht, The Netherlands). Three other environmental strains isolated from marine plant, leak and sea in Shandong, Inner Mongolia and Hebei provinces, respectively, were purchased from China General Microbiological Culture Collection Center (CGMCC) in Beijing, China. The animal strain was isolated from the nasal secretions of a dog in a pet hospital in Beijing. Yeast isolates were inoculated onto Sabouraud dextrose agar (SDA) at 35°C for 2 days. The study was approved by the Human Research Ethics Committee of Peking Union Medical College Hospital (no. S-K690).
rDNA typing based on SNPs of ITS and D1/D2 regions Yeast cells of R. mucilaginosa were grown on Yeast Extract Peptone Dextrose (YPD) medium at 30°C overnight prior to DNA extraction. Genomic DNA was extracted using the Fungi Genomic DNA Extraction Kit (Solarbio Science & Technology, Beijing, China) according to the company's recommended protocols. The internal transcribed spacers (ITS), 5.8S region of the ribosomal DNA (rDNA) and the D1/D2 domains of the largest subunit rDNA were amplified using primer pairs and PCR conditions as previously described [13]. The isolates were identified and confirmed based on sequence analysis of the ITS and D1/D2 regions. Corresponding sequences (NR_073296.1 and NG_055716.1) from the type strain (R. mucilaginosa CBS 316) were retrieved from NCBI database. Sequence alignment was carried out using CLC Sequence Viewer software v8.0 (Qiagen, Denmark).

Establishment of microsatellite typing system
Tandem repeat loci were screened from the reference genome (GenBank assembly accession: GCA_000931965.1) of R. mucilaginosa strain C2.5t1 using Tandem Repeats Finder v 4.09 (https:// tandem.bu.edu/trf/trf.basic.submit.html). A total of 161 potential loci (RM001 to RM161) were discovered, and non-labelled forward and reverse primers targeting these loci were designed using Primer3Plus (http://www.primer3plus.com/). Gradient PCR was performed to determine the optimal annealing temperature for all primer pairs. A panel of 15 R. mucilaginosa strains isolated in different geographic regions and years (listed in Supplementary Table 1) were used in preliminary screening to verify the universality of these 161 loci. Each PCR mix (25 µL) contained 5 µL of 5 × PrimeSTAR Buffer (Mg 2+ Plus) (TAKARA Biomedical Technology, Beijing, China), 2 µL of dNTP Mixture, 0.25 µL of PrimeSTAR HS DNA Polymerase, 0.5 µM of forward and reverse primers, 0.5 µL of DNA template, and molecular biology grade water. The reactions were performed as follow: initial denaturation at 98°C for 2 min, followed by 30 cycles of denaturation at 98°Cfor 10 sec, annealing at 60°C for 30 sec, extension at 72°C for 30 sec, with a final extension at 72°C for 5 min.
The forward primers of 136 loci were then labelled with either 6-carboxyfluores-cein (FAM) or hexachlorofluorescein (HEX). Following PCR, amplicons were sized on an ABI 3730XL DNA Analyzer (Applied Biosystems, Foster City, CA, USA) coupled with Gen-eMarker v2.2 software (SoftGenetics LLC, State College, PA, USA). Allele sizes were scored with respect to GeneScan™ 500 LIZ® Size Standard (Applied Biosystems) in the 35-500 bp range. Polymorphic alleles of each potential microsatellite locus were recorded and interpreted manually. Finally, 57 loci labelled four kinds of fluorescent (FAM, HEX, 6-carboxytetramethylrhodamine [TAMRA] or carboxy-X-rhodamine [ROX]) were used to evaluate the polymorphism of all R. mucilaginosa isolates. Among them, 15 loci were selected and combined to achieve the highest resolution. The discriminatory power (DP) for 15 loci was calculated using the Simpson index as previously reported [14].

Genetic relationship analysis based on microsatellite
The genetic relationship of Rhodotorula strains was inferred by constructing a minimum spanning tree using the BioNumerics software v7.6 (Applied Maths, Sint-Martens-Latem, Belgium). Samples with genotypes showing the same alleles for all loci were considered identical clones. Epidemic genotypes were defined as genotypes infecting ≥3 different patients in one hospital [15]. DNA concentration was quantified, and sequencing libraries were constructed using NEBNext® Ultra™ DNA Library Prep Kit for Illumina (NEB, MA, USA) following manufacturer's recommendations, and different barcodes were added to each sample. DNA libraries were sequenced on the Illumina Nova-Seq platform in Novogene Co. (Beijing, China) with 2 × 150-bp (paired-end) reads.
Using fastx_toolkit tools, quality control was conducted for the sequencing data [16]. A threshold of 0.01 (Phred score of 20) was used for the trimming of the raw Illumina sequencing reads. Paired-end reads were mapped to the R. mucilaginosa JGTA-S1 reference genome (GCA_003055205.1) using BWA 0.7.17 [17]. Variant calling was performed for all 32 R. mucilaginosa. Single Nucleotide Polymorphisms (SNPs) were analyzed using SAMtools 0.1.18 and bcftools (part of the SAMtools package) [18]. We filtered out low-quality SNPs with a sequencing depth less than 30 and a ratio of reads supporting the mutation less than 0.8. A phylogenetic tree of genomic SNPs was constructed with software MEGA X [19] using the Maximum Likelihood method. The confidence of the topology was evaluated by bootstrapping with 1000 randomizations. Heatmap of comparison between microsatellite genotypes and genomic SNPs were generated by "pheatmap" of R program based on "complete" clustering method [20]. Gradient colour represented the numbers of pairwise genomic SNPs (bp). Illumina reads from this study have been deposited in National Centre for Biotechnology Information (NCBI) under BioProject accession number PRJNA794005.

In vitro susceptibility testing
The in vitro susceptibility to nine antifungal drugs, including three echinocandins (caspofungin, micafungin, and anidulafungin), four azoles (fluconazole, voriconazole, itraconazole, and posaconazole), amphotericin B, and 5-flucytosine, was determined using broth microdilution method with Sensititre YeastOne™ YO10 methodology (Thermo Scientific, Cleveland, OH, USA) following the manufacturer's instructions. Candida parapsilosis ATCC 22019 and Candida krusei ATCC 6258 were used as quality control. Results were interpreted in accordance with cutoff values of MICs for Candida spp. to initially determine antifungal susceptibility following the manufacturer's instructions.

Statistical analyses
Statistical analyses were performed using one-way ANOVA to compare differences between strains of different genotypes by SPSS 22.0 (IBM Corp. Released 2013. IBM SPSS Statistics for Windows, Version 22.0. Armonk, NY: IBM Corp.). The differences were considered statistically significant if p was <0.05.

Clinical characteristics of patients
In the 49 clinical cases, 63.3% were males. The patient age ranged from 23 to 89 years, and the average was 55 years, and the quartiles were 39, 51 and 66 years, respectively. 38.8% and 42.9% were from patients admitted to the ICUs and surgical departments, respectively. The most common specimen was venous blood (85.7%), followed by ascitic fluid (4.1%), pus (4.1%), pleural fluid (2.0%), catheter (2.0%) and tissue (2.0%). We examined the medical records of 43 invasive infection cases, which revealed that 27.9%, 27.9% and 20.9% of the patients had hypoproteinemia, tumour and diabetes, respectively (Table 1). Most patients (88.4%) were prescribed prophylaxis with antimicrobial agents, but only 46.5% of the patients received systematic antifungal treatment (azoles and/ or echinocandins), and the overall mortality was 11.6% (5/43).

Distribution of rDNA types of strains
All strains from 49 patients were identified as R. mucilaginosa by sequencing ITS and D1/D2 regions of rDNA. Their polymorphic sequences can be divided into four and three molecular types of ITS and D1/D2 sequences, respectively. Combining sequences of the two regions of rDNA, six rDNA genotypes were identified among the clinical strains ( Table 2). Strains of the rDNA type 2 were prevalent, including 36 isolates (73.5%) from four hospitals. In hospital H01, 92.3% (24/26) of strains belonged to rDNA type 2. To further characterize the genetic relationship of the clinical strains, especially strains from hospital H01, a highresolution molecular typing was performed with microsatellites.

Phylogenetic analysis based on microsatellite genotypes
The microsatellite genotyping method with 15 tandem repeat loci has been established after studying the polymorphism of 57 markers with 15 representative strains ( Table 3, Table S1). In total, 32 different microsatellite (MT) types were identified for 49 clinical strains with 0.97 discriminatory power (Figure 2(A)). 26 strains of hospital H01 belonged to 12 MT types and 23 strains from other three hospitals had 20 MT types (Figure 2(B)). Among 5 clinical epidemic MT types (one MT type strains infected ≥ 3 patients), 4 epidemic MT types (MT03, MT04, MT05 and MT06) including 18 strains were from hospital H01, and one MT type (MT01) including 4 strains were from hospital H02 (Table S2). The remaining 27 MT types were distributed sporadically over different hospitals and time (Figure 2(B)).
Overlaid with the rDNA types on the minimum spanning tree of the microsatellite, we found that 35 clinical rDNA type 2 strains clustered together which included 4 epidemic MT types (MT03 to MT06) from hospital H01 and 17 sporadic MT types from 4 hospitals (epidemic cluster, Figure 2(C)). Twenty-four rDNA type 2 strains from hospital H01 all belonged to the epidemic cluster. Only one rDNA type 2 strain from hospital H03 did not belong to the epidemic cluster. However, two rDNA type 4 strains and two rDNA type 5 strains from hospital H02 were identified as the same MT type -MT01 (Figure 2(B,C)). Comparing genomic SNPs of strains with different MT types, it was found that the pairwise SNPs difference of representative strains was consistent with genetic relationship of strains in minimum spanning tree of the microsatellites (Figure 4(A)). Strains of different MT genotypes, belonging to epidemic cluster in MT tree, clustered together with small pairwise SNPs difference (within 1,200 bp, Figure 4(B)). Particularly, in strains of epidemic MT types (MT03, MT04, MT05 and MT06) from hospital H01, pairwise genomic SNPs difference of the earliest and the latest strains with same MT type was less than 200 bp (36 bp-188 bp, Figure 4(B)). In addition, SNPs difference of the two strains of MT01 with different rDNA types (rDNA type 4 and 5) from hospital H02 was 4441 bp.  (Table S3). However, in hospital H02, case 5 (MT30) and case 6 (MT28) were also collected in a short period, and their pairwise SNPs was 668 bp. In addition, although four patients (case 3, 4, 7 and 8) from hospital H02 all underwent being transferred between ICU and surgical wards in six months ( Figure  1), their isolates did not cluster together in MT tree.

Relationship among strains from patients, environment and pet
In order to find out the genetic relationship of strains from different sources, we included four environmental strains and one animal strain of R. mucilaginosa in the phylogenetic analysis. The environmental strain CGMCC 2.5541 isolated from a lake in Inner Mongolia and the animal strain A343 isolated from nasal secretions of a pet dog nested within the epidemic cluster. Together with the clinical strain (H01-06) from hospital H01, they all belonged to the same MT type (MT02) (Figure 4(B), Table S3). In addition, an endophytic strain JGTA-S1 of Typha angustifolia grown in India (from GenBank) also clustered with them in the tree (Figure 3). The other three environmental strains CGMCC 2.2506 (MT17), CGMCC 2.5690 (MT11) and CBS 316 (MT07) were divergent and did not belong to the epidemic cluster, having pairwise SNPs difference greater than 10,000 bp ( Figure 3).

Antifungal susceptibility profiles
Antifungal susceptibility profiles were interpreted after 48 h incubation. All clinical strains had identical susceptibility to echinocandins (caspofungin, Figure 1. Hospital course, history and distribution of microsatellite types included in the study. Different colours represent different clinical departments. Number in colourful square means the serial number of cases in each hospital. The red dot on the right side of the number indicates that the strain of the case belongs to the epidemic cluster. The red rectangle shows that patients underwent a transfer between two wards. *, the area where the colour block is located shows the part time of the patient in hospital because of unavailable clinical records. **, the case was diagnosed in an outpatient clinic. micafungin and anidulafungin), fluconazole and amphotericin B, with MICs of ≥8, ≥64 and ≤1 µg/ mL, respectively (Table 4). According to CLSI M27 guidelines, ≤4, 8-16, and ≥32 µg/mL of MICs were interpreted as susceptible, intermediate and resistant to 5-FC against Candida spp., respectively. In this study, MICs of 81.6% (40/49) of clinical strains were ≤0.12 µg/mL. While MICs of the remaining 9 strains of 5-FC were all ≥8 µg/mL. The correlation analysis between antifungal susceptibility phenotypes and MT types found that all 9 strains with MIC of 5-FC at ≥8 µg/mL belonged to MT05 and MT06 types in the epidemic cluster (Figure 2(D) and Figure 3), also had the same colony morphology ( Figure S1). Strains from hospital H01 were more susceptible to voriconazole than strains from the other three hospitals in the epidemic cluster (p = 0.001).

Discussion
Rhodotorula mucilaginosa is an opportunistic pathogenic yeast [5,21,22]. It is oligotrophic and can grow over a broad temperature range (from 0.5°C to 37°C ) and even survive in heavy metal-rich sludge [23][24][25]. As the main carotenoid-forming yeast, R. mucilaginosa was widely used to biosynthesize carotenoids with low-cost substrates for industry. Good N 2 -fixing ability and biocontrol to postharvest decay of fruits make it popular in agriculture [26][27][28]. The extensive use of industry and agriculture makes it highly adaptable to the environment. Although R. mucilaginosa is a rare, emerging pathogen, reports of invasive infections have been increasing, especially in Asia [7].
Analyses of microsatellites and SNPs indicated the presence of the epidemic cluster of R. mucilaginosa in China. This opportunistic pathogen is widely distributed in China because highly genetically similar strains can be collected from four hospitals in various geographic regions from 2012 to 2019 (Figure 1 and Figure 2(B)). We hypothesize that the pathogenic strain may be originated from the natural environment, and the patients acquired the colonization/ infection during their interaction with soil/plants.  According to the phylogenetic analysis, our clinical strains clustered strongly with the environmental strain from an endophytic strain of Typha angustifolia grown in India, and the lake in Inner Mongolia, as well as an animal strain from a dog (see Figure 4), suggesting the fungus can be dispersed from the environment to humans/animals. Indeed, Rhodotorula species have been isolated from fruit salads and juice purchased from the supermarkets, and they were reported in the microbial communities of Chinese milk fan and smear surface-ripened cheeses as early as the beginning of twenty-first century [29][30][31].
Once the pathogen persists in the hospitals, it will pose a public health threat. Many fungi can cause zoonotic infections, such as Sporothrix schenckii, Coccidioides posadasii, Enterocytozoon bieneusi, and most of them can cause infections to dogs and cattle through soil-transmission [32][33][34]. R. mucilaginosa has been isolated from cattle and chickens, and it has caused respiratory infections in dogs [35,36]. In this study, the strain was isolated from the nasal secretions of a pet dog and clustered with many clinical strains (Figure 2(A) and Figure 3).
Although the pet dog did not belong to any patients included in the investigation, the same fungus has the capacity to cause infection in both humans and animals.
Since multiple strains in the epidemic cluster, especially strains of the same MT types were isolated in the hospital H01 in different years (see Figure 1 and Figure 2(B)), we speculate that R. mucilaginosa may be able to persist/survive in the hospital environment. The pairwise SNPs differences were only 36 and 104 bp within strains of MT05 and MT03 respectively, with MIC towards 5-FC at ≥8 µg/mL and ≤4 µg/mL, respectively. Given the low genetic polymorphisms and consistent phenotype of 5-FC resistance, there may be clonal expansion/transmissions in hospital H01 since 2012.
We speculate that strains in the epidemic cluster had persisted in the hospital (H01) environment and entered the bloodstream via invasive operations to cause fungemia. Although patients of case 2 to case 4 (MT05) were admitted to different wards (for instance, cardiac surgery, general surgery and general ICU), they all went through invasive operations, placed with central venous catheters, so did cases with strains of MT03 type (case 8 and 9). In fact, catheterrelated bloodstream infections were the most common form of infection caused by R. mucilaginosa [37]. Although only the first isolate of each patient was included in the investigation, we found that five patients (case 2 in 2012, case 7 in 2013 and case 18, 20, 21 in 2016) had multiple positive blood cultures of R. mucilaginosa in hospital H01 after reviewing their clinical records.
Even though WGS data has high-resolution power, it is rather expensive to utilize the technique in a routine operation. The microsatellite genotyping scheme (developed in this investigation) could be a useful alternative to differentiate closely related strains of R. mucilaginosa. In this investigation, genetic relationships could be resolved for epidemic strains (MT02 to MT06) and those beyond the epidemic cluster (MT01) when their rDNA had low resolution. For molecular epidemiological investigation in a clinical microbiology laboratory, rDNA typing method can be used for preliminary screening of epidemic cluster followed by microsatellite genotyping to differentiate their relatedness.  Antifungal susceptibilities of R. mucilaginosa to echinocandins, azoles and AMB in this study were consistent with other report [6]. However, the strains were also resistant to 5-FC, there was limited choice in antifungal agents for treating the infections (Table 4, Figure 2(D)). Resistance to 5-FC may be associated with the pathway critical for the formation of capsule in Cryptococcus [38]. Since R. mucilaginosa has a thin capsule, further study would be required to understand the mechanism of 5-FC resistance and its relationship to capsule production.
Global guideline recommended AMB with or without 5-FC (moderately recommended) as first-line therapy for the treatment of Rhodotorula spp. infections [7]. Reviewing the medical records indicated that all patients (n = 20) receiving medical treatment were not treated by AMB, which have good in vitro activity against R. mucilaginosa, with overall mortality rate of 15.0% (3/20). If these patients had received treatment of AMB, the outcomes may be improved.
This study has a few limitations. First, we cannot clarify the route of transmissions and the specific site/environment where these strains attach to from the original patients. In future, we should collect more samples from the hospital environment to determine the reservoir of these persistent strains. Second, we had only four environmental strains and an animal strain. Additional samples of environmental and zoonotic sources would be beneficial to determine their relationship with the clinical strains. Third, we used cut-off values of MICs for Candida spp. to initially determine antifungal susceptibility due to the absence of clinical breakpoint for R. mucilaginosa.

Conclusion
Our study characterized the molecular epidemiology of an emerging pathogen R. mucilaginosa in China. Both microsatellites and genomics SNPs showed the presence of an epidemic cluster of strains collected from multiple hospitals in different geographic regions and the occurrence of a 5-FC resistant clone. Inferring from the genetic relatedness among clinical, animal and environmental data, the pathogen could spread from the environment to humans and animals, and pose public health risk to some degree because it causes invasive infections and may be able to persist in the hospital environment. The funding bodies had no role in study design, data collection and analysis, interpretation of data, decision to publish, or preparation of the manuscript.

Data availability
Genome raw reads are available at National Centre for Biotechnology Information (NCBI) under BioProject accession number PRJNA794005.
Microsatellite genotypes with fragment sizes of 15 tandem repeat loci were listed in Table S4. Minimum inhibitory concentration of 9 antifungal agents against R. mucilaginosa clinical strains was listed in Table S5. Phenotypic characteristics of colonies of R. mucilaginosa after 7 days incubation were showed in Figure S1.