Draft genome of a commonly misdiagnosed multidrug resistant pathogen Candida auris
- Sharanya Chatterjee†1,
- Shuba Varshini Alampalli†1,
- Rishi Kumar Nageshan1,
- Sivarajan T. Chettiar1,
- Sangeeta Joshi2 and
- Utpal S. Tatu1Email author
© Chatterjee et al. 2015
Received: 3 February 2015
Accepted: 18 August 2015
Published: 7 September 2015
Candida auris is a multidrug resistant, emerging agent of fungemia in humans. Its actual global distribution remains obscure as the current commercial methods of clinical diagnosis misidentify it as C. haemulonii. Here we report the first draft genome of C. auris to explore the genomic basis of virulence and unique differences that could be employed for differential diagnosis.
More than 99.5 % of the C. auris genomic reads did not align to the current whole (or draft) genome sequences of Candida albicans, Candida lusitaniae, Candida glabrata and Saccharomyces cerevisiae; thereby indicating its divergence from the active Candida clade. The genome spans around 12.49 Mb with 8527 predicted genes. Functional annotation revealed that among the sequenced Candida species, it is closest to the hemiascomycete species Clavispora lusitaniae. Comparison with the well-studied species Candida albicans showed that it shares significant virulence attributes with other pathogenic Candida species such as oligopeptide transporters, mannosyl transfersases, secreted proteases and genes involved in biofilm formation. We also identified a plethora of transporters belonging to the ABC and major facilitator superfamily along with known MDR transcription factors which explained its high tolerance to antifungal drugs.
Our study emphasizes an urgent need for accurate fungal screening methods such as PCR and electrophoretic karyotyping to ensure proper management of fungemia. Our work highlights the potential genetic mechanisms involved in virulence and pathogenicity of an important emerging human pathogen namely C. auris. Owing to its diversity at the genomic scale; we expect the genome sequence to be a useful resource to map species specific differences that will help develop accurate diagnostic markers and better drug targets.
KeywordsNosocomial infections Drug resistance Virulence Fungemia Misdiagnosis Candida haemulonii Next Generation Sequencing (NGS)
Hospital acquired infection (HAI) also known as nosocomial infections are gaining momentum and Centre for Disease control USA, estimates about 99,000 deaths a year due to infections acquired from hospital . HAI are caused by organisms that include bacteria and fungi entering via surgical sites, urinary and other catheters. Nosocomial associated invasive fungal diseases are of public health concern and candidemia is becoming very prevalent in European countries . The fourth most common cause of bloodstream infection is Candida, accounting for more than 85 % of all fungemias in USA and Europe [3, 4]. However the rising number of immunocompromised people, unwarranted use of multiple broad spectrum antibiotics and the advent of implanted medical devices  has paved the way for rare non albicans Candida species  as agents of invasive mycoses and nosocomial bloodstream infections . Recently invasive non albicans candidiasis cases have been reported from many parts of the world [8–11]. Among the non albicans Candida species, C. tropicalis and C. glabrata have emerged as important opportunistic pathogens .
Most recently, species belonging to the Candida haemulonii complex  has been described as important agent of fungemia with a significant global distribution [14, 15] and Lehman et al.  categorised these species belonging to C. haemulonii into two genetically distinct groups. Furthermore infections caused by two phenotypically related species – C. pseudohaemulonii  and C. auris are on the rise . First described in 2009 by Satoh et al.  in a Japanese patient, it is striking to see the aggressive pace at which C. auris has expanded its clinical spectrum worldwide from minor cases of superficial infections such as ear canal infections to highly invasive cases of bloodstream infections . Previous studies  as well as our study demonstrate that all these clinical isolates have a precociously high tolerance to AmphotericinB (AmB)  and Fluconozole (Fcz) [14, 15], the first line treatment antifungals. Even more concerning is the rapid emergence of resistance to echinocandins , the newest class of antifungals which may leave no treatment option available leading to clinical failure.
Many pathogenic species within the Candida clade such as Candida albicans and Candida glabrata have been extensively studied at the genome level, while emerging fungal pathogens Candida auris and Candida haemulonii remains unexplored. The basic characteristics of the genome of C. auris was recently made available . However detailed information regarding the genome architecture, virulence and mechanisms of multidrug resistance of these emerging novel complexes of pathogenic yeasts are lacking. Furthermore, the commercial automated systems routinely fail to identify C. auris correctly; thereby its actual occurrence is underreported. Even more alarming is the fact that misdiagnosis may lead to incorrect treatment or delay of proper treatment, thereby increasing the chances of fatalities. As expected, C. auris fungemia is associated with a high mortality rate (66 %) and therapeutic failure . It also does not exhibit the known attributes responsible for virulence in Candida species such as hyphae formation and the cells are much smaller in size than that of C. albicans (Additional file 1: Figure S2). Towards understanding the basic biology of the multidrug resistant pathogen, we have carried out whole genome sequencing of a multidrug resistant clinical isolate of C. auris using Illumina sequencing technology and report that C. auris has a highly divergent genome. Analysis using C. albicans as a reference genome revealed a set of orthologs such as drug transporters, oligopeptide transporters, secreted proteinases and mannosyl transferases which may play a role in virulence and drug resistance. However most of the genome is uncharacterized and we speculate that some of these hypothetical proteins may be involved in species specific characteristics which promote its aggressiveness as a pathogen.
Results and discussion
Clinical isolates of Candida show multi drug resistance
In vitro antifungal susceptibility pattern of pathogenic Candida clinical isolates (from bloodstream) to the most commonly used drugs belonging to the four different classes of antifungals
MIC50 range (μg/ml)
Clinical isolates_ Cha (34)
C. albicans (20)
C. tropicalis (34)
C. glabrata (9)
C. lusitaniae (2)
C. parapsilosis (23)
Complete genome sequence of the clinical isolate Ci 6684
We sequenced the genome of Ci 6684 using Illumina sequencing technology. A high-quality reference genome using Illumina reads was assembled de novo as described in Methods (Additional file 3: Figure S1). The assembled draft genome of Ci 6684 comprises 99 scaffolds with an estimated genome size of 12,498,766 bp, 44.53 % GC and 1.327 % Ns. The average base is found in the scaffold with a scaffold N50 of 279 Kb. A total of 8358 protein coding genes, 7 rRNAs and 189 tRNAs were predicting using different tools (Description in Methods and Additional file 3: Figure S1).
In vitro antifungal susceptibility profile of Candida clinical isolate Ci 6684
General features of Candida species and clinical isolate Ci 6684 genome
Number of chromosomes (or scaffolds)
GC content (%)
Number of genes
Average CDS size (bp)
Average intergenic distance (bp)
C. albicans SC5314
C. dubliniensis CD36
C. orthopsilosis Co 90–125
C. tropicalis MYA-3404
C. guilliermondii ATCC 6260
C. lusitaniae ATCC 42720
C. glabrata CBS 138
S. cerevisiae S288C
Phylogenetic analysis reveals Ci 6684 is closely related to Candida auris
Candida auris has a highly divergent genome
Genomic relatedness calculated using Genome-to-Genome Distance Calculator
Formula 1 (HSP length/total length)
Formula 2 (identities/HSP length)
Formula 3 (identities/total length)
Prob. DDH ≥ 70 %
Prob. DDH ≥ 70 %
Prob. DDH ≥ 70 %
C. auris 6684
C. glabrata CBS 138
C. auris 6684
C. albicans SC5314
C. auris 6684
C. albicans WO-1
C. auris 6684
C. lusitaniae ATCC 42720
C. auris 6684
S. cerevisiae S288c
Functional annotation of the C. auris 6684 genome
We also performed enzyme classification analysis based on Enzyme Commission (EC) numbers predictions for each sequence. We found that hydrolases are the largest group of C. auris 6684 enzymes (42 %), followed by transferases (25 %) and oxidoreductases (19 %). Blast2GO identified 466 enzyme (Fig. 6d) out of which 329 enzymes got mapped to KEGG pathways. BlastKOALA was used to reconstruct KEGG pathways for C. auris 6684. 2775 proteins (out of 8358 predicted proteins) got annotated into various pathways. This analysis revealed that the central pathways pertaining to carbohydrate, lipid and amino acid metabolisms are conserved.
Core circuitry related to virulence is conserved in C. auris 6684
Considering the high genomic variability of C. auris 6684, we asked the question that whether gene families that are known to have a role in pathogenicity of Candida species  are also conserved in C. auris 6684? We used the genome of C. albicans SC5314 as the template gene model to predict orthologs in our isolate as it is well annotated at the experimental level. This approach yielded 1988 orthologous proteins with functional annotations. Our analysis predicted an arsenal of transporters orthologous to that of C. albicans, belonging primarily to the major facilitator superfamily and ABC (ATP binding cassette) superfamily  (Fig. 6d). The up regulation of these multidrug efflux pumps may explain the intrinsically low susceptibility of C. auris 6684 to antifungal drugs. Apart from the general transcription factors, 193 proteins were predicted to have DNA binding/sequence specific DNA binding/transcription factor activity. We also predicted a multitude of zinc finger transcription factors orthologous to those present in Saccharomyces cerevisiae, Candida albicans and Scheffersomyces stipites. Notably the Zn (II) 2 Cys 6 transcription factor family is enriched in our isolate (26 in number). Four of these are known to be key regulators of MDR1 transcription in C. albicans; gain-of-function mutations of which leads to up regulation of multidrug efflux pump MDR1, thereby leading to multidrug resistance [32, 33].
The genome was found to contain transcription factors like STE-related and MADS-box proteins which have been previously shown to be involved in the virulence of human fungal pathogens [34, 35] and plant fungal pathogens [36, 37] respectively. Ste12p is conserved in many fungi, regulating processes involved in mating, filamentation, substrate invasion, cell wall integrity and virulence , while MADS-box proteins bind to DNA and have dimerization activity . Our analysis also indicated conservation of the Rim101 transcriptional pathway that is known to respond to alkaline pH in Saccharomyces cerevisiae. 122 proteins were predicted to have kinase/phosphorylation activity. Out of this, 93 proteins have the serine/threonine kinase domain and the rest were predicted to be involved in protein phosphorylation due to the presence of putative kinase domain/ATP binding domain. C. auris 6684 draft genome encodes for kinases like Hog1, Protein Kinase A (PKA) and two-component histidine kinase. Activation of stress signaling pathways regulated by these protein kinases have been implicated to enhance tolerance of pathogenic fungi to chemical fungicides and antifungal peptides . HOG1 protein is a fungal mitogen-activator protein (MAP) kinase which has been implicated in responses to oxidative and hyperosmotic stresses in a few human pathogens including C. albicans . PKA is shown to be activated in response to extracellular nutrients and subsequently regulates metabolism and growth, while two-component histidine kinase is shown to be critical to morphogenesis and virulence [31, 40, 41].
We also identified eight OPT genes encoding putative oligopeptide transporters which have been implicated in the acquisition of nutrient versatility thereby helping the pathogen to adapt to various host niches . Interestingly it has been reported that in C. albicans, these genes are also induced upon phagocytosis by macrophages . We also found orthologs of genes predicted to be hexose transporters, maltose transporters and permeases (amino acid permeases, sulfur permeases, allantoate permeases, glycerol permeases and iron permeases) which further expands its nutrient assimilation machinery, thereby helping it to acclimatize to diverse host niches.
Our next step was to hunt down the attributes that may explain the aggressive behavior of the pathogen. Our analysis indeed predicted many known virulence associated genes (Fig. 6e). Since the cell wall serves as the interface between the pathogen and the host immune defense, components of the cell wall serve as pathogen associated molecular patterns and virulence factors. Our analysis indicated that the family of mannosyl transferases is conserved in C. auris 6684 with many predicted orthologs. Apart from maintaining cell wall architecture by coordinating glycan synthesis, these enzymes play a very important role in immune recognition, host cell adherence and virulence in C. albicans . Integrins and adhesins are the other two gene families which have a crucial role in adherence and virulence of C. albicans [45, 46]. However our annotation predicted only two proteins, one structurally similar to alpha-subunit of human leukocyte integrin; predicted to play a role in morphogenesis, adhesion, mouse coecal colonization and virulence; and another secreted protein similar to alpha agglutinin anchor subunit which has been previously shown to be induced upon exposure to fluconazole. This clearly suggests that C. auris employs distinct mechanisms for host cell adhesion.
We also found four orthologs of secreted aspartyl proteases (SAP) two of which were predicted to have greater expression upon deep epidermal invasion; greater expression in vaginal than oral infection  and prominent role in biofilm formation. We also found two genes annotated as vacuolar aspartic proteinases. The secreted aspartic proteinases help the fungus to digest host proteins and the resulting peptides are taken up into the cell by specific transporters like the oligopeptide transporters family mentioned above . Our results also annotated eight genes orthologous to secreted lipases. In all, our analysis revealed that enzyme families implicated in invasiveness like mannosyl transferases, secreted aspartyl proteases and lipases are enriched in our clinical isolate. However the adhesion and integrin gene families are ill represented. This information has been categorized in Additional file 2: Table S3. Our analysis also revealed 686 proteins predicted to be induced or repressed upon rat catheter or biofilm formation. This includes a multitude of enzymes, transcription factors, ribosomal proteins and transporters. This clearly indicates that C. auris 6684 has significant ability to form biofilms since the core genes involved in biofilm formation are conserved. However experiments need to be done to validate the same.
Structure of mating loci in C. auris and PCR based diagnostic test to differentiate between C. auris and C. haemulonii
Another peculiarity seen in Candida species is the highly diverse nature of sexuality. Diploids like C. tropicalis and C. parapsilosis are unable to mate while C. albicans shows a parasexual cycle. Haploids like C. lusitaniae and C. gulliermondii are heterothaliic in nature . It is interesting to note that virulence and mode of reproduction are being analysed as linked phenomenon in recent years. C. lusitaniae is a heterothallic species known to be involved in sexual reproduction. On the other hand certain Candida species are either parasexual or asexual. Considering the high similarity shared by C. auris 6684 and C. lusitaniae, we speculated that C. auris 6684 might have a sexual stage similar to the latter. Sexual mating is controlled by a single genetic locus called the MAT locus consisting of two alleles-MATa and MATα.
Opportunistic infections caused by Candida are on the rise globally and newer pathogenic species are emanating at an unprecedented rate. What are not evolving at the same pace are the current methods of diagnosis and treatment options leading to misdiagnosis and clinical failure. The last decade has witnessed the emergence of newer species called C. haemulonii from being the causative agent of minor infections to one of the leading causes of invasive infections. It is currently increasing in prevalence, with several ongoing outbreaks in developing and underdeveloped countries. The actual incidence rate is however misleading because of the inability of the current automated systems used for screening of fungal species to identify novel emerging fungal pathogens such as C. auris, C. pseudohaemulonii and other related species due to striking similarities in biochemical characters and the unavailability of molecular markers for accurate identification. We have generated the first draft genome sequence of a commonly misdiagnosed, emerging pathogen C. auris. The isolate was identified as C. haemulonii by Vitek2. However PFGE analysis revealed 5 bands similar to that of C. auris and accurate species identification was done by phylogenetic analysis based on the partial sequence of 18S rRNA, ITS1, 5.8S rRNA complete sequence, ITS2 and 28S rRNA partial sequence. Genome sequencing will highlight important differences which may act as accurate identification markers for this group of emerging pathogens at the species level. Towards this we have developed a PCR based diagnostic test to distinguish between these two pathogens.
The genome of C. auris spans about ~12.5 Mb with 8358 predicted protein coding genes. Strikingly, at the genomic level, C. auris shows a highly divergent relationship with other pathogenic Candida species as indicated by a meagre 0.5 % alignment of the sequencing reads to other Candida genomes and supported by lack of linear synteny of genomic dot plots. C. auris is phylogenetically closest to C. haemulonii whose genome sequence is unavailable. Among the sequenced yeast species, it is closest to C. lusitaniae; however, its genome is also not well annotated functionally. Therefore majority of the protein coding genes were predicted to be hypothetical/functionally uncharacterized. The role of each of these unique candidate proteins demands for urgent functional studies. Hence accurate identification and de novo assembly and annotation still remains a challenge for divergent sequences among emerging pathogenic species. 37.71 % of the protein coding genes showed no sequence similarity to genes available in public database, thus indicating that speciation genes are embedded within the genome which may be involved in grooming it as an aggressive pathogen. With the limited data available, it is difficult to comment about the genomic architecture of speciation and how it facilitates or impedes further divergence. To further probe into the difference at the functional level we resorted to synonymous codon usage plots which distinguish ways by which translational selection of protein coding genes occurs among related species. The above observation is supported by GGDC that calculated the in silico relatedness of C. auris and sequenced Candida pathogens, surprisingly the logistic regression quantifies no relatedness among the species. The ecological niche of most of these Candida species is known, that may throw light on the evolutionary forces grooming these organisms at the species level. However till date there are no reports of naturally occurring C. auris species. C. auris can grow at elevated temperatures of 42 °C whereas C. haemulonii cannot. This gives us a hint that C. auris has the potential to infect the avian fauna whose body temperature is in the range of 40 °C to 42 °C. However, additional experiments need to be done in order to validate this phenomenon.
The foremost criterion to be a successful Candida pathogen is the ability to colonize diverse anatomical niches within the host such as skin, oral cavity, gastrointestinal tract, vagina and the vasculature. Each Candida pathogen has its own machinery dedicated to host cell adhesion, recognition, invasion and colonization. We compared C. auris genome with that of C. albicans since it is well annotated and well-studied as well as distantly related to C. auris. While the spectrum of virulence traits like hyphae formation, white opaque switching is quite different between these two species, we found that C. auris still shares some common virulence traits with C. albicans. Our analysis highlights that a significant portion of C. auris genome encodes for transporters belonging to the ABC transporter family and major facilitator superfamily. This may partly explain its increased tolerance to antifungal drugs. The multidrug resistant nature of the pathogen and the limited arsenal of antifungal agents indicate that there is a critical need for finding new drug targets and genome sequence of C. auris therefore may prove useful in finding alternative targets that can augment the existing antifungal therapy. Our analysis also provides a snapshot of the potential genetic attributes that may explain its virulent nature. The genome of the pathogen harbours gene families such as lipases, oligopeptide transporters, mannosyl transferases and transcription factors which play a multitude of roles in colonization, invasion and iron acquisition. Also majority of genes known to be involved in formation of biofilm appears to be conserved. In all, we see that C. auris shares many genes with C. albicans and C. lusitaniae indicating a common ancestry; however it may have acquired novel genetic traits that have groomed it as a specialist pathogen. It is possible that the indiscriminate use of antibiotics shaped its genome to expand not only its clinical spectrum of infection but also to emerge as a successful multidrug resistant pathogen.
In all, our study provides the first whole genomic overview of C. auris, the first member of the Candida haemulonii and related pathogenic fungi complex to be sequenced. This report is a major step toward the initiation of genomic studies of this complex group of fungi which are fast turning drug resistant and may be a menace with limited treatment options available in the future.
Strain and growth conditions
All clinical isolates were obtained from Manipal Hospital, Bengaluru and the ethical approval was obtained from Ethics Committee of Manipal Hospitals, Bengaluru and informed consent was taken as required during the study. Ci 6684 was isolated from a patient who had sepsis with multiorgan dysfunction. C. haemulonii 8176 was obtained from MTCC, IMTECH Chandigarh, India. Strains were routinely grown in Yeast Peptone Dextrose (YPD) medium at 37 °C.
Minimum inhibitory concentration and growth assays
To determine the in vitro susceptibility to antifungal drugs, broth microdilution protocol  was used. Overnight cultures were grown at 37 °C in YPD. Approximately 103 cells per well in YPD media at 37 °C. Minimum inhibitory concentration (MIC) tests were set up in a total volume of 0.2 ml/well with 2-fold dilutions of drugs. Fluconazole gradients where in the following concentration steps in μg/ml: 64, 32, 16, 8, 4, 2, 1, 0.5, 0.25, 0.125, 0.0625 and 0.03125. For Amphotericin B, gradients where in the following the concentration steps in μg/ml were: 16, 8, 4, 2, 1, 0.5, 0.25, 0.125, 0.0625, 0.03125 and 0.015625. 24 or 48 h post incubation, growth was measured by reading the optical density at 600 nm after agitation using a spectrophotometer (Tecan). MIC50 was defined as the concentration of drug reducing growth by 50 % relative to the wells containing no drug. Sterile water was the vehicle for Fcz and AmB.
Short reads and long reads library preparation was performed at Genotypic Technology’s Genomics facility following NEXTFlex DNA library protocol outlined in “NEXTFlex DNA sample preparation guide (Cat # 5140–02). ~3 μg of genomic DNA was sonicated using Bioruptorto and 300 to 600 bp sized fragments were obtained. The size distribution was checked by running an aliquot of the sample on Agilent HS DNA Chip. The resulting fragmented DNA was cleaned up using Agencourt AMPure XP SPRI beads (Beckman Coulter). Fragmented DNA was subjected to a series of enzymatic reactions that repair frayed ends, phosphorylate the fragments, and add a single nucleotide A overhang and ligate adaptors (NEXTFlex DNA Sequencing kit). Sample cleanup was done using AMPure SPRI beads. After ligation-cleanup, ~300–600 bp fragments was size selected on 2 % low melting agarose gel and cleaned using MinElute column (QIAGEN). PCR (10 cycles) amplification of adaptor ligated fragments was done and cleaned up using AMPure SPRI beads. The prepared libraries were quantified using Qubit flourometer and validated for quality by running an aliquot on High Sensitivity Bioanalyzer Chip (Agilent). The short read inserts were sequenced in Illumina MiSeq and long read inserts were sequenced in Illumina NextSeq 500.
Mate-pair reads library preparation was performed at Genotypic Technology’s Genomics facility following Nextera Mate Pair Gel Plus protocol outlined in “Illumina Nextera Mate Pair library preparation guide (Cat# FC-132-9001DOC, Part#15035209 Rev D.)”. ~4 μg of Qubit quantified DNA was taken for Tagmentation. The tagmented sample was cleaned up using AMPure beads and subjected to strand displacement. 3–5 kb range of the strand displaced sample was size selected on 0.6 % agarose gel. Size selected sample was taken for circularization overnight, followed by linear DNA digestion with DNA Exonuclease. The circularized DNA molecules were sheared using Covaris to obtain fragments in the size range of 300 to 1000 bp. Sheared DNA was subjected to bead binding with M280 Streptavidin beads to isolate biotinylated molecules. End repair, A-Tailing and adapter ligations were performed on the bead-DNA complex. Adaptor ligated sample was amplified for 15 cycles of PCR followed by AMPure XP bead clean up. The prepared library was quantified using Qubit and validated for quality by running an aliquot on High Sensitivity Bioanalyzer Chip (Agilent). The mate-pair reads were sequenced using Illumina NextSeq 500.
Assembly, annotation and analysis
The qualities of the reads were checked using Genotypic proprietary tool SeqQC v2.21. The average sequencing depth (coverage) for short paired-end reads is 158.19x, long paired-end reads is 175.51x and mate-pair reads is 205.78x. Processed short paired-end reads (3.27 million) were used to generate (250–400) long fragments using ARF-PE v0.2. 467178 long fragments were generated using 467178*2 paired end reads (ie, 14.29 % reads were used in long read generation). 467178 long fragments and 3269025*2 paired end reads used for Newbler Genome assembly. Newbler version 2.8’s default assembly parameters were used for the assembly and 721 scaffolds were generated. The paired-end long insert reads and mate-pair reads were used to gap fill using SSPACE-STANDARD v3.0  and the contigs were reduced to 65 scaffolds. Using Reapr v1.0.17 , the 65 scaffolds were corrected, removing the erroneous bases and the final number of scaffolds was 97. These 97 scaffolds were used as input for GeneMarkS  to predict protein-coding genes with –eukaryotic as the main option. The resulting 8388 proteins were subjected to local blastp, resulting in 5175 proteins being annotated to RefSeq fungal protein database. Proteins having query coverage of greater than 40 % were only considered from this blast results. An InterproScan  was carried out using the tool Blast2GO  v3.0 to group the predicted proteins according to the presence of domain/motif in their sequences. GO terms were assigned through Blast2GO tool based on NR Database orthologs (blastp with Evalu > e−10). Proteins involved in various KEGG pathways were assigned using BlastKOALA . Transfer RNAs were identified using the tRNAScan-SE program . Ribozomal RNAs were predicted by RNAmmer . The sequenced reads were mapped to various pathogenic Candida genome using Bowtie2 v2.2.3  with default parameter. The generated SAM files were used to calculate the percent of reads aligned using R.
Modified PFGE, Counter-clamped homogeneous electrical field (CHEF) (BIO-RAD) was used for electrophoretic karyotyping of C. auris 6684 and C. albicans. The protocol was adapted from Iadonato et al. 1996 . Briefly 5 ml yeast cultures were grown in YPD medium at 30 °C. The cells were the harvested and washed with 50 mM EDTA. Approximately 2× 109 cells/ml were added to equal volumes of 1 % (w/v) low melt Pulse Field certified Agarose (BIO-RAD), prewarmed at 45 °C. The mixture was then transferred in to disposable plug moulds to harden. Plugs were then extruded and suspended in freshly prepared spheroplasting solution containing Zymolase, and incubated at 37 °C for 4 h. After this the plugs were washed with 1 % Lithium dodecyl sulfate (LDS) (2X 30 min) buffer followed by cell lysis with 1 % N-lauryl sarcosine (NDS) (3X 30 min) buffer. Finally the plugs were rinsed (6x 30 min) with TE buffer pH 8. Agarose plugs containing yeast DNA was then loaded into 0.8 % low melt Pulse Field certified Agarose (BIO-RAD) prepared with 0.5X TBE buffer. The DNA samples were resolved by running the gel in CHEF-DR® III system with 5 V/cm2 with pulse time of 120 s and total run time of 36 h at 12 °C. Gel was then stained with ethidium bromide (1ug/ml) for 30 min and visualized at ImageQuant LAS 4000 transilluminator (GE).
Phylogenetic tree and evolutionary analysis
The partial sequence of 18 s rRNA, ITS1, 5.8 s rRNA complete sequence, ITS2 and 28 s rRNA partial sequence retrieved from NCBI (Additional file 2: Table S5) were used to categorise Clinical isolate 6684 with Candida auris clade. The evolutionary tree was inferred using the Maximum Likelihood method based on the Tamura-Nei model . The tree with the highest log likelihood (−307.3435) is shown. The percentage of trees in which the associated taxa clustered together is shown next to the branches. Initial tree(s) for the heuristic search were obtained automatically by applying Neighbor-Join and BioNJ algorithms to a matrix of pairwise distances estimated using the Maximum Composite Likelihood (MCL) approach, and then selecting the topology with superior log likelihood value. The tree is drawn to scale, with branch lengths measured in the number of substitutions per site. The analysis involved 48 nucleotide sequences. All positions containing gaps and missing data were eliminated. There were a total of 167 positions in the final dataset.
95 conserved proteins (Additional file 2: Table S2) from Saccharomyces cerevisiae S288c were retrieved using YGD, CGD and BLASTn for the following organisms: Saccharomyces cerevisiae S288c, Candida albicans SC-5314, Candida dubliniensis CD-36, Candida glabrata CBS 138,,Candida isolate 6684, Candida tropicalis MYA-3404, Candida lusitaniae ATCC 42720, Candida gulliermondii ATCC 6260, Candida orthopsilosis Co-90–125, Ashbya_gossypii and Histoplasma capsulatum. The phylogenetic tree was constructed using the Neighbor-Joining method. The optimal tree with the sum of branch length = 1.22757517 is shown. The percentage of replicate trees in which the associated taxa clustered together in the bootstrap test (2000 replicates) is shown next to the branches. The tree is drawn to scale, with branch lengths in the same units as those of the evolutionary distances used to infer the phylogenetic tree. The evolutionary distances were computed using the p-distance method and are in the units of the number of amino acid differences per site. The analysis involved 11 amino acid sequences. All positions with less than 95 % site coverage were eliminated. There were a total of 51712 positions in the final dataset.
Tajima’s neutrality analysis involved concatenated amino acid sequences from the 11 species. All positions with less than 95 % site coverage were eliminated. There were a total of 51712 positions in the final dataset. The equality of evolutionary rate between Candida lusitaniae, Clinical isolate 6684 with Candida albicans as an out-group was determined by Tajima’s relative rate test [64, 65]. All positions containing gaps and missing data were eliminated. There were a total of 56989 positions in the final dataset. All the phylogenetic trees and evolutionary analyses were conducted in MEGA6  .
For genome comparison the current genome sequences (whole or draft) were downloaded from Broad Institute (https://www.broadinstitute.org/scientific-community/science/projects/fungal-genome-initiative/fungal-genomics) and CGD (www.candidagenome.org/). The analysis was carried out using GFFex v2.3 and Biostrings package of Bioconductor in R v3.1. The DNA-DNA hybridizations (DDH) distances were calculated using the online tool Genome-to-Genome Distance Calculator (GGDC 2.0) (http://ggdc.dsmz.de/). Dot plot were done in an online tool called YASS  by setting the e-value to e-10 and the synonymous codon usage plots were done in R (v3.1) using ape4 and seqinr packages  of Bioconductor.
Polymerase chain reaction
Genomic DNA was isolated as described previously. Based on the MFα region sequence from C. auris, a specific PCR-based method was developed for the direct detection of C. auris DNA by using a C. auris -specific primer (CaMF [5′- GAGAAAAGAGACGCTGAAGCTGAG-3′]) designed using the gene sequence which codes for the unique pheromone together with reverse primer (CaMR [5′- TCAACCTTCGAGGTCAGCTTCA-3′]).
Ploidy analysis by FACS
Cultures were grown in YPD till A600 of 1.0. The cells were washed in 1X PBS (137 mM NaCl, 2.7 mM KCl, 10 mM sodium phosphate dibasic (NaH2PO4), 2 mM potassium phosphate monobasic (K2HPO4), pH of 7.4) and fixed in 70 % ethanol for 1 h at room temperature or kept at 4 °C overnight. The cells were suspended in 1X PBS and incubated with RNase A (1 mg/ml) at 37 °C for 4 h in the same buffer. Cells were subsequently washed with PBS, and finally stained with propidium iodide (PI, 16 μg/ml) for flow cytometric analysis in BD FACS Canto.
Availability of supporting data
The whole genome sequencing data can be accessed through BioProject accession number PRJNA267757. The respective BioSample accession numbers is SAMN03200169. The SRA reference numbers of the whole genome sequencing are SRX766223 (Illumina MiSeq short paired-end reads), SRX766234 (Illumina NextSeq 500 mate-pair reads) and SRX766231 (Illumina HiSeq2500 long paired-end reads). This Whole Genome Shotgun project has been deposited at DDBJ/EMBL/GenBank under the accession LGST00000000. The version described in this paper is version LGST01000000.
Internal transcribed spacer
Hospital acquired infections
Minimum inhibitory concentration
Pulse field gel electrophoresis
Genomic to genomic distance calculator
Kyoto encyclopedia of genes and genomes
ATP binding cassette
Secreted aspartyl proteinases
The authors would like to acknowledge Genotypic Technology, Bangalore, India. We acknowledge funding from the DBT-IISc partnership program and Grant Challenge Canada (Sub-grant fund: 494417). Research fellowship from DST INSPIRE for Sharanya Chatterjee is acknowledged.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- Klevens RM, Edwards JR, Richards Jr CL, Horan TC, Gaynes RP, Pollock DA, et al. Estimating health care-associated infections and deaths in U.S. hospitals, 2002. Public Health Rep. 2007;122:160–6.PubMed CentralPubMedGoogle Scholar
- Lass-Florl C. The changing face of epidemiology of invasive fungal disease in Europe. Mycoses. 2009;52:197–205.View ArticlePubMedGoogle Scholar
- Quindos G. Nosocomial candidemias and invasive candidiasis. Med Clin (Barc). 2010;134:17–9.View ArticleGoogle Scholar
- Tortorano AM, Kibbler C, Peman J, Bernhardt H, Klingspor L, Grillot R, et al. Candidaemia in Europe: epidemiology and resistance. Int J Antimicrob Agents. 2006;27:359–66.View ArticlePubMedGoogle Scholar
- Adhikary R, Joshi S. Species distribution and anti-fungal susceptibility of Candidaemia at a multi super-specialty center in Southern India. Indian J Med Microbiol. 2011;29:309–11.View ArticlePubMedGoogle Scholar
- Pfaller MA, Andes DR, Diekema DJ, Horn DL, Reboli AC, Rotstein C, et al. Epidemiology and outcomes of invasive candidiasis due to non-albicans species of Candida in 2,496 patients: data from the Prospective Antifungal Therapy (PATH) registry 2004–2008. PLoS One. 2014;9, e101510.PubMed CentralView ArticlePubMedGoogle Scholar
- Pfaller MA, Diekema DJ, Procop GW, Rinaldi MG. Multicenter comparison of the VITEK 2 antifungal susceptibility test with the CLSI broth microdilution reference method for testing amphotericin B, flucytosine, and voriconazole against Candida spp. J Clin Microbiol. 2007;45:3522–8.PubMed CentralView ArticlePubMedGoogle Scholar
- Papon N, Courdavault V, Clastre M, Bennett RJ. Emerging and emerged pathogenic Candida species: beyond the Candida albicans paradigm. PLoS Pathog. 2013;9, e1003550.PubMed CentralView ArticlePubMedGoogle Scholar
- Colombo AL, Nucci M, Park BJ, Nouer SA, Arthington-Skaggs B, da Matta DA, et al. Epidemiology of candidemia in Brazil: a nationwide sentinel surveillance of candidemia in eleven medical centers. J Clin Microbiol. 2006;44:2816–23.PubMed CentralView ArticlePubMedGoogle Scholar
- Colombo AL, Garnica M, Aranha Camargo LF, Da Cunha CA, Bandeira AC, Borghi D, et al. Candida glabrata: an emerging pathogen in Brazilian tertiary care hospitals. Med Mycol. 2013;51:38–44.View ArticlePubMedGoogle Scholar
- Hachem R, Hanna H, Kontoyiannis D, Jiang Y, Raad I. The changing epidemiology of invasive candidiasis: Candida glabrata and Candida krusei as the leading causes of candidemia in hematologic malignancy. Cancer. 2008;112:2493–9.View ArticlePubMedGoogle Scholar
- Turner SA, Butler G. The Candida pathogenic species complex. Cold Spring Harb Perspect Med. 2014;4:a019778.View ArticlePubMedGoogle Scholar
- Cendejas-Bueno E, Kolecka A, Alastruey-Izquierdo A, Theelen B, Groenewald M, Kostrzewa M, et al. Reclassification of the Candida haemulonii complex as Candida haemulonii (C. haemulonii group I), C. duobushaemulonii sp. nov. (C. haemulonii group II), and C. haemulonii var. vulnera var. nov. three multiresistant human pathogenic yeasts. J Clin Microbiol. 2012;50:3641–51.PubMed CentralView ArticlePubMedGoogle Scholar
- Khan ZU, Al-Sweih NA, Ahmad S, Al-Kazemi N, Khan S, Joseph L, et al. Outbreak of fungemia among neonates caused by Candida haemulonii resistant to amphotericin B, itraconazole, and fluconazole. J Clin Microbiol. 2007;45:2025–7.PubMed CentralView ArticlePubMedGoogle Scholar
- Kim MN, Shin JH, Sung H, Lee K, Kim EC, Ryoo N, et al. Candida haemulonii and closely related species at 5 university hospitals in Korea: identification, antifungal susceptibility, and clinical features. Clin Infect Dis. 2009;48:e57–61.View ArticlePubMedGoogle Scholar
- Lehmann PF, Wu LC, Pruitt WR, Meyer SA, Ahearn DG. Unrelatedness of groups of yeasts within the Candida haemulonii complex. J Clin Microbiol. 1993;31:1683–7.PubMed CentralPubMedGoogle Scholar
- Sugita T, Takashima M, Poonwan N, Mekha N. Candida pseudohaemulonii Sp. Nov. an amphotericin B-and azole-resistant yeast species, isolated from the blood of a patient from Thailand. Microbiol Immunol. 2006;50:469–73.View ArticlePubMedGoogle Scholar
- Lee WG, Shin JH, Uh Y, Kang MG, Kim SH, Park KH, et al. First three reported cases of nosocomial fungemia caused by Candida auris. J Clin Microbiol. 2011;49:3139–42.PubMed CentralView ArticlePubMedGoogle Scholar
- Satoh K, Makimura K, Hasumi Y, Nishiyama Y, Uchida K, Yamaguchi H, et al. Candida auris sp. nov. a novel ascomycetous yeast isolated from the external ear canal of an inpatient in a Japanese hospital. Microbiol Immunol. 2009;53:41–4.View ArticlePubMedGoogle Scholar
- Chowdhary A, Sharma C, Duggal S, Agarwal K, Prakash A, Singh PK, et al. New clonal strain of Candida auris, Delhi, India. Emerg Infect Dis. 2013;19:1670–3.PubMed CentralView ArticlePubMedGoogle Scholar
- Sarma S, Kumar N, Sharma S, Govil D, Ali T, Mehta Y, et al. Candidemia caused by amphotericin B and fluconazole resistant Candida auris. Indian J Med Microbiol. 2013;31:90–1.View ArticlePubMedGoogle Scholar
- Rodero L, Cuenca-Estrella M, Cordoba S, Cahn P, Davel G, Kaufman S, et al. Transient fungemia caused by an amphotericin B-resistant isolate of Candida haemulonii. J Clin Microbiol. 2002;40:2266–9.PubMed CentralView ArticlePubMedGoogle Scholar
- Muro MD, Motta Fde A, Burger M, Melo AS. Dalla-Costa LM Echinocandin resistance in two Candida haemulonii isolates from pediatric patients. J Clin Microbiol. 2012;50:3783–5.View ArticlePubMedGoogle Scholar
- Sharma C, Kumar N, Meis JF, Pandey R, Chowdhary A. Draft genome sequence of a fluconazole-resistant Candida auris strain from a candidemia patient in India. Genome Announc. 2015;3:e00722–15.PubMed CentralView ArticlePubMedGoogle Scholar
- Chowdhary A, Anil Kumar V, Sharma C, Prakash A, Agarwal K, Babu R, et al. Multidrug-resistant endemic clonal strain of Candida auris in India. Eur J Clin Microbiol Infect Dis. 2014;33:919–26.View ArticlePubMedGoogle Scholar
- Oh BJ, Shin JH, Kim MN, Sung H, Lee K, Joo MY, et al. Biofilm formation and genotyping of Candida haemulonii, Candida pseudohaemulonii, and a proposed new species (Candida auris) isolates from Korea. Med Mycol. 2011;49:98–102.View ArticlePubMedGoogle Scholar
- Kim HY, Huh HJ, Choi R, Ki CS, Lee NY. Three cases of candidiasis misidentified as Candida famata by the Vitek 2 system. Ann Lab Med. 2015;35:175–7.PubMed CentralView ArticlePubMedGoogle Scholar
- Ochiuzzi ME, Cataldi S, Guelfand L, Maldonado I, Arechavala A. Evaluation of Vitek 2 for the identification of Candida yeasts. Rev Argent Microbiol. 2014;46:107–10.PubMedGoogle Scholar
- Kathuria S, Singh PK, Sharma C, Prakash A, Masih A, Kumar A, et al. Multidrug-resistant Candida auris misidentified as Candida haemulonii: characterization by matrix-assisted laser desorption ionization-time of flight mass spectrometry and DNA sequencing and its antifungal susceptibility profile variability by vitek 2, CLSI broth microdilution, and etest method. J Clin Microbiol. 2015;53:1823–30.View ArticlePubMedGoogle Scholar
- Kent WJ. BLAT--the BLAST-like alignment tool. Genome Res. 2002;12:656–64.PubMed CentralView ArticlePubMedGoogle Scholar
- Calderone RA, Fonzi WA. Virulence factors of Candida albicans. Trends Microbiol. 2001;9:327–35.View ArticlePubMedGoogle Scholar
- Ramage G, Bachmann S, Patterson TF, Wickes BL, Lopez-Ribot JL. Investigation of multidrug efflux pumps in relation to fluconazole resistance in Candida albicans biofilms. J Antimicrob Chemother. 2002;49:973–80.View ArticlePubMedGoogle Scholar
- Sanglard D, Kuchler K, Ischer F, Pagani JL, Monod M, Bille J, et al. Mechanisms of resistance to azole antifungal agents in Candida albicans isolates from AIDS patients involve specific multidrug transporters. Antimicrob Agents Chemother. 1995;39:2378–86.PubMed CentralView ArticlePubMedGoogle Scholar
- Qu X, Yu B, Liu J, Zhang X, Li G, Zhang D, et al. MADS-box transcription factor SsMADS is involved in regulating growth and virulence in Sclerotinia sclerotiorum. Int J Mol Sci. 2014;15:8049–62.PubMed CentralView ArticlePubMedGoogle Scholar
- Calcagno AM, Bignell E, Warn P, Jones MD, Denning DW, Mühlschlegel FA, et al. Candida glabrata STE12 is required for wild-type levels of virulence and nitrogen starvation induced filamentation. Mol Microbiol. 2003;50:1309–18.View ArticlePubMedGoogle Scholar
- Ortiz CS, Shim WB. The role of MADS-box transcription factors in secondary metabolism and sexual development in the maize pathogen Fusarium verticillioides. Microbiology. 2013;159:2259–68.View ArticlePubMedGoogle Scholar
- Mehrabi R, Ding S, Xu JR. MADS-box transcription factor mig1 is required for infectious growth in Magnaporthe grisea. Eukaryot Cell. 2008;7:791–9.PubMed CentralView ArticlePubMedGoogle Scholar
- Hayes BM, Anderson MA, Traven A, van der Weerden NL, Bleackley MR. Activation of stress signalling pathways enhances tolerance of fungi to chemical fungicides and antifungal proteins. Cell Mol Life Sci. 2014;71:2651–66.View ArticlePubMedGoogle Scholar
- Alonso-Monge R, Navarro-Garcia F, Molero G, Diez-Orejas R, Gustin M, Pla J, et al. Role of the mitogen-activated protein kinase Hog1p in morphogenesis and virulence of Candida albicans. J Bacteriol. 1999;181:3058–68.PubMed CentralPubMedGoogle Scholar
- Calera JA, Choi GH, Calderone RA. Identification of a putative histidine kinase two-component phosphorelay gene (CaHK1) in Candida albicans. Yeast. 1998;14:665–74.View ArticlePubMedGoogle Scholar
- Yamada-Okabe T, Mio T, Ono N, Kashima Y, Matsui M, Arisawa M, et al. Roles of three histidine kinase genes in hyphal development and virulence of the pathogenic fungus Candida albicans. J Bacteriol. 1999;181:7243–7.PubMed CentralPubMedGoogle Scholar
- Reuss O, Morschhauser J. A family of oligopeptide transporters is required for growth of Candida albicans on proteins. Mol Microbiol. 2006;60:795–812.View ArticlePubMedGoogle Scholar
- Lorenz MC, Bender JA, Fink GR. Transcriptional response of Candida albicans upon internalization by macrophages. Eukaryot Cell. 2004;3:1076–87.PubMed CentralView ArticlePubMedGoogle Scholar
- Hall RA, Bates S, Lenardon MD, Maccallum DM, Wagener J, Lowman DW, et al. The Mnn2 mannosyltransferase family modulates mannoprotein fibril length, immune recognition and virulence of Candida albicans. PLoS Pathog. 2013;9, e1003276.PubMed CentralView ArticlePubMedGoogle Scholar
- Hostetter MK. Adhesins and ligands involved in the interaction of Candida spp. with epithelial and endothelial surfaces. Clin Microbiol Rev. 1994;7:29–42.PubMed CentralPubMedGoogle Scholar
- Kinneberg KM, Bendel CM, Jechorek RP, Cebelinski EA, Gale CA, Berman JG, et al. Effect of INT1 gene on Candida albicans murine intestinal colonization. J Surg Res. 1999;87:245–51.View ArticlePubMedGoogle Scholar
- Naglik JR, Rodgers CA, Shirlaw PJ, Dobbie JL, Fernandes-Naglik LL, Greenspan D, et al. Differential expression of Candida albicans secreted aspartyl proteinase and phospholipase B genes in humans correlates with active oral and vaginal infections. J Infect Dis. 2003;188:469–79.View ArticlePubMedGoogle Scholar
- Naglik JR, Challacombe SJ, Hube B. Candida albicans secreted aspartyl proteinases in virulence and pathogenesis. Microbiol Mol Biol Rev. 2003;67:400–28. table of contents.PubMed CentralView ArticlePubMedGoogle Scholar
- Reedy JL, Floyd AM, Heitman J. Mechanistic plasticity of sexual reproduction and meiosis in the Candida pathogenic species complex. Curr Biol. 2009;19:891–9.PubMed CentralView ArticlePubMedGoogle Scholar
- Srikantha T, Daniels KJ, Pujol C, Sahni N, Yi S, Soll DR, et al. Nonsex genes in the mating type locus of Candida albicans play roles in a/alpha biofilm formation, including impermeability and fluconazole resistance. PLoS Pathog. 2012;8, e1002476.PubMed CentralView ArticlePubMedGoogle Scholar
- Rustad TR, Stevens DA, Pfaller MA, White TC. Homozygosity at the Candida albicans MTL locus associated with azole resistance. Microbiology. 2002;148:1061–72.View ArticlePubMedGoogle Scholar
- Cowen LE, Lindquist S. Hsp90 potentiates the rapid evolution of new traits: drug resistance in diverse fungi. Science. 2005;309:2185–9.View ArticlePubMedGoogle Scholar
- Boetzer M, Pirovano W. SSPACE-LongRead: scaffolding bacterial draft genomes using long read sequence information. BMC Bioinformatics. 2014;15:211.PubMed CentralView ArticlePubMedGoogle Scholar
- Hunt M, Kikuchi T, Sanders M, Newbold C, Berriman M, Otto TD, et al. REAPR: a universal tool for genome assembly evaluation. Genome Biol. 2013;14:R47.PubMed CentralView ArticlePubMedGoogle Scholar
- Borodovsky M, Lomsadze A. Gene identification in prokaryotic genomes, phages, metagenomes, and EST sequences with GeneMarkS suite. Curr Protoc Microbiol. 2011;32(Unit 1E):7.Google Scholar
- Mulder N, Apweiler R. InterPro and InterProScan: tools for protein sequence classification and comparison. Methods Mol Biol. 2007;396:59–70.View ArticlePubMedGoogle Scholar
- Gotz S, Garcia-Gomez JM, Terol J, Williams TD, Nagaraj SH, Nueda MJ, et al. High-throughput functional annotation and data mining with the Blast2GO suite. Nucleic Acids Res. 2008;36:3420–35.PubMed CentralView ArticlePubMedGoogle Scholar
- Kanehisa M, Goto S, Sato Y, Kawashima M, Furumichi M, Tanabe M, et al. Data, information, knowledge and principle: back to metabolism in KEGG. Nucleic Acids Res. 2014;42:D199–205.PubMed CentralView ArticlePubMedGoogle Scholar
- Lowe TM, Eddy SR. tRNAscan-SE: a program for improved detection of transfer RNA genes in genomic sequence. Nucleic Acids Res. 1997;25:955–64.PubMed CentralView ArticlePubMedGoogle Scholar
- Lagesen K, Hallin P, Rodland EA, Staerfeldt HH, Rognes T, Ussery DW, et al. RNAmmer: consistent and rapid annotation of ribosomal RNA genes. Nucleic Acids Res. 2007;35:3100–8.PubMed CentralView ArticlePubMedGoogle Scholar
- Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012;9:357–9.PubMed CentralView ArticlePubMedGoogle Scholar
- Iadonato SP, Gnirke A. RARE-cleavage analysis of YACs. Methods Mol Biol. 1996;54:75–85.PubMedGoogle Scholar
- Tamura K, Nei M. Estimation of the number of nucleotide substitutions in the control region of mitochondrial DNA in humans and chimpanzees. Mol Biol Evol. 1993;10:512–26.PubMedGoogle Scholar
- Tajima F. Simple methods for testing the molecular evolutionary clock hypothesis. Genetics. 1993;135:599–607.PubMed CentralPubMedGoogle Scholar
- Tamura K, Battistuzzi FU, Billing-Ross P, Murillo O, Filipski A, Kumar S, et al. Estimating divergence times in large molecular phylogenies. Proc Natl Acad Sci U S A. 2012;109:19333–8.PubMed CentralView ArticlePubMedGoogle Scholar
- Tamura K, Stecher G, Peterson D, Filipski A, Kumar S. MEGA6: Molecular Evolutionary Genetics Analysis version 6.0. Mol Biol Evol. 2013;30:2725–9.PubMed CentralView ArticlePubMedGoogle Scholar
- Noe L, Kucherov G. YASS: enhancing the sensitivity of DNA similarity search. Nucleic Acids Res. 2005;33:W540–3.PubMed CentralView ArticlePubMedGoogle Scholar
- Charif D, Thioulouse J, Lobry JR, Perriere G. Online synonymous codon usage analyses with the ade4 and seqinR packages. Bioinformatics. 2005;21:545–7.View ArticlePubMedGoogle Scholar