key: cord-0916388-h6q7ktn7
authors: Robinson, Sarah; Peterson, Christine B.; Sahasrabhojane, Pranoti; Ajami, Nadim J.; Shelburne, Samuel A.; Kontoyiannis, Dimitrios P.; Galloway-Peña, Jessica R.
title: Observational Cohort Study of Oral Mycobiome and Interkingdom Interactions over the Course of Induction Therapy for Leukemia
date: 2020-04-15
journal: mSphere
DOI: 10.1128/msphere.00048-20
sha: dc2e9ca36ae9e4514dddde0980925658fe2accd6
doc_id: 916388
cord_uid: h6q7ktn7

Although the term “microbiome” refers to all microorganisms, the majority of microbiome studies focus on the bacteriome. Here, we characterize the oral mycobiome, including mycobiome-bacteriome interactions, in the setting of remission-induction chemotherapy (RIC) for acute myeloid leukemia (AML). Oral samples (n = 299) were prospectively collected twice weekly from 39 AML patients during RIC until neutrophil recovery. Illumina MiSeq 16S rRNA gene (V4) and internal transcribed spacer 2 (ITS2) sequencing were used to determine bacterial and fungal diversity and community composition. Intrakingdom and interkingdom network connectivity at baseline (T1) and at midpoint (T3) and a later time point (T6) were assessed via SPIEC-EASI (sparse inverse covariance estimation for ecological association inference). In this exploratory study, mycobiome α-diversity was not significantly associated with antibiotic or antifungal receipt. However, postchemotherapy mycobiome α-diversity was lower in subjects receiving high-intensity chemotherapy. Additionally, greater decreases in Malassezia levels were seen over time among patients on high-intensity RIC compared to low-intensity RIC (P = 0.003). A significantly higher relative abundance of Candida was found among patients who had infection (P = 0.008), while a significantly higher relative abundance of Fusarium was found among patients who did not get an infection (P = 0.03). Analyses of intrakingdom and interkingdom relationships at T1, T3, and T6 indicated that interkingdom connectivity increased over the course of IC as bacterial α-diversity diminished. In (to our knowledge) the first longitudinal mycobiome study performed during AML RIC, we found that mycobiome-bacteriome interactions are highly dynamic. Our study data suggest that inclusion of mycobiome analysis in the design of microbiome studies may be necessary to optimally understand the ecological and functional role of microbial communities in clinical outcomes. IMPORTANCE This report highlights the importance of longitudinal, parallel characterization of oral fungi and bacteria in order to better elucidate the dynamic changes in microbial community structure and interkingdom functional interactions during the injury of chemotherapy and antibiotic exposure as well as the clinical consequences of these interrelated alterations.

the core set of oral mycobiome taxa previously identified among healthy individuals, which includes Candida, Cladosporium, Saccharomycetales, Aspergillus, Fusarium, Cryptococcus, and Malassezia (22, 23) . Although Aspergillus was seen in at least one oral sample in over 70% of the patients, it had a detectable abundance in only 13% of all longitudinal samples. Additionally, among the samples which had a detectable abundance of Aspergillus, the median relative abundance was only 1.8% (interquartile range [IQR] , 0.03% to 8.4%). This observation was interesting given that Aspergillus is the primary cause of invasive mold infection among patients with hematologic malignancy (24) . The abundances of the top ten most abundant genera showed significant intrapatient and interpatient variability across time (Fig. 1) . However, the ␤-diversity appeared to be more homogeneous within each patient across time, with the mean intrapatient Bray-Curtis dissimilarity between samples within the same patient being 0.649 among the cohort members versus an average interpatient ␤-diversity value of 0.745 for all samples across time points. This is consistent with another study that analyzed the oral mycobiome stability over 30 days, which revealed high interindividual diversity and yet indicated intraindividual stability (25) . Using the results of principal-coordinate analysis (PCoA) shown in Fig. 2 , we found no distinct grouping by patient or time point based on Bray-Curtis dissimilarity. We next sought to describe the oral mycobiome changes over the course of RIC. Through mixed-effects models, we found that mycobiome ␣-diversity did not significantly change over time (Chao1, P ϭ 0.4639; Shannon, P ϭ 0.8049; Simpson, P ϭ 0.6637) (Fig. 3) , unlike what has been previously been noted for the bacteriome, where diversity decreased from baseline until neutrophil recovery among leukemia patients receiving RIC (2, 26) . Mycobiome characteristics correlate with chemotherapy subtype and antifungal administration. We next sought to determine if mycobiome profiles were associated with clinical variables among the patients in the cohort, to include chemotherapeutic or antimicrobial administration. In this sampling of patients, 72% were on high-intensity chemotherapy, while 28% received low-intensity chemotherapy (Table 1). The majority of patients (44%) received fludarabine-based regimens, with the majority of the regimens being FIA (idarubicin, cytarabine, and fludarabine) or FLAG-Ida (fludarabine, idarubicin, and granulocyte-colony-stimulating factor). Twenty-eight percent of patients received nonfludarabine high-intensity regimens, which typically included CIA (idarubicin, cytarabine, and clofarabine)-based, IA (idarubicin and cytarabine)-based, or CLIA (cladribine, idarubicin, and cytarabine)-based regimens. Patients also receive hypomethylator-based treatments (such as decitabine or azacitadine) (13%) or other low-dose therapies to include LDAC (low-dose antimetabolite cytarabine)-based regimens (14%). Given the uneven sampling time (time to neutrophil recovery) among patients, we evaluated the mycobiome structure and composition at the baseline (T1), a midpoint (T3), and a later time point (T6), common among most patients. At T6, patients who received low-intensity RIC had significantly higher mycobiome ␣-diversity than patients who received high-intensity RIC (observed operational taxonomic units [OTUs], P ϭ 0.02; Shannon, P ϭ 0.054) (see Fig. S1 in the supplemental material). Although Kruskal-Wallis tests indicated that the relative abundances of the five most abundant fungal genera did not significantly differ between the four different chemotherapy subtypes at any single time point (see Table S1 in the supplemental material), we did determine chemotherapy subtypes to be significantly associated with the relative abundance of specific taxa over time. Specifically, patients who received high-intensity chemotherapies experienced a decrease in Malassezia abundance over time (P ϭ 0.003, Fig. 4 ). This effect was not seen in other taxa analyzed (data not shown). To our knowledge, the changes in the oral mycobiome occurring as a result of chemotherapy have been studied only in regard to their role in chemotherapy-induced mucositis in the setting of 5-fluorouracil (5-FU) or doxorubicin-based chemotherapy in mice and humans (13, 14) . In patients, although the authors expected chemotherapy to affect the oral mycobiome, no major changes in the composition of fungal communities were observed (13) . Thus, our study would be the first to show specific changes in the oral mycobiome associated with chemotherapeutic intensity and type.

We next sought to determine if antimicrobial administration resulted in significant mycobiome changes during RIC for AML. We examined the three primary empirical broad-spectrum antibiotics received in this cohort, which were carbapenems, cephalosporins, and piperacillin-tazobactam, and found that 87% of the patients were administered at least one of these three broad-spectrum antimicrobials. Some patients received more than one broad-spectrum antibiotic prior to neutrophil recovery, with 67% receiving a carbapenem, 44% a cephalosporin, and 28% piperacillin-tazobactam. 

Changes in Malassezia over time stratified by chemotherapy intensity and chemotherapy type. Centered log ratio (clr) data represent transformed relative abundances of Malassezia with fitted lines from the mixed-effects model over time according to each variable of interest. P value data were derived from ANOVA F-test results. Panel A has fitted lines representing the intensities of the chemotherapy regimens, with high-intensity regimens indicated in coral and low-intensity regimens in aqua. Panel B further dissects these by chemotherapy subtype, with high-intensity fludarabine-containing regimens indicated in coral, hypomethylator-based therapies in green, non-fludarabine-containing high-intensity regimens in aqua, and other low-intensity/low-dose cytarabine-based therapies in purple.

We also evaluated the three primary treatment antifungals, consisting of echinocandins, amphotericin B, and triazoles. All patients were administered at least one of these three antifungals, with 79% receiving an azole, 77% an echinocandin, and 15% amphotericin B. We analyzed the later time point during RIC (T6) given that all antimicrobials analyzed were administered by T6 and that it was the latest time point common among most patients. Surprisingly, there was no significant association between the relative abundance of the top 5 most abundant fungal taxa at time point 6 and receipt of specific antimicrobials (Tables S2 toS7) . However, amphotericin B seemed to have a significant effect on the composition of the oral fungal community over time. Patients who received amphotericin B experienced an increase in the relative abundance of Fusarium over time (P ϭ 0.034) (Fig. 5 ). There is biological plausibility to this observation inasmuch as Fusarium species are associated with high amphotericin B MICs (27) . Amphotericin B had no significant effect on the relative abundance of other taxa over time.

Mycobiome characteristics are correlated with bacterial infection outcomes. In light of the complex interplay between fungal populations and opportunistic bacterial pathogens (28), we sought to determine the relationships between the mycobiome and infections in AML patients during RIC. A total of 23% of patients had a microbiologically defined infection with a positive bacterial culture, which included mostly bacteremia or urinary tract infections (Table 1) . We also considered the 31% of patients who had a clinically defined infection, such a pneumonia or cellulitis, in the absence of positive culture in our analysis. In analyzing the three major time points (T1, T3, and T6), we found that the relative abundance of Candida at T3 was significantly higher in patients that contracted microbiologically or clinically defined infections prior to neutrophil recovery (P ϭ 0.008), while there was a significantly higher relative abundance of Fusarium (P ϭ 0.03) at T6 in those that did not develop an infection prior to neutrophil recovery (Fig. 6) . These results were corroborated using the LefSe (linear discriminant analysis of effect size) algorithm, where Candida was more differentially abundant at T3 among patients who developed an infection prior to neutrophil recovery (Fig. S2) and Fusarium more abundant at T6 among patients who did not have an infection (Fig. S3) . No other statistical associations were observed for other genera at the three time points or at T1 (data not shown). Interkingdom interactions of the oral microbiome over time in AML patients receiving RIC. Cross-domain association networks were constructed at T1, T3, and T6 using matching bacterial and fungal microbiome samples (n ϭ 272). The three SPIEC-EASI (sparse inverse covariance estimation for ecological association inference) networks show vast differences in intra-and interkingdom connectivity between time points, highlighting the dynamic nature of bacterial and fungal interactions induced by changes in the community during RIC (Fig. 7) . The average node degrees were 0.43 (standard deviation [SD], Ϯ0.48), 0.84 (SD, Ϯ0.92), and 0.5 (SD, Ϯ0.56) for T1, T3, and T6, respectively.

Relative abundance differences of taxa at individual time points among AML patients who did and did not experience bacterial infection during RIC. Box plots depict the taxon abundance differences of Candida at T3 and Fusarium at T6 between those patients that did (aqua) and did not (coral) contract bacterial infections prior to neutrophil recovery following RIC. P values are derived from Mann-Whitney test. Some of the intrakingdom and interkingdom interactions have been previously observed. Specifically, Malassezia and Candida are often coisolated from infectious samples from the skin infections (29) and onychomycosis (30) , but to our knowledge cooccurrence has not been shown in the oral cavity. Interestingly, we found a negative correlation between Malassezia and Candida at T1 and T3 in the oral cavity. In a study of the oral microbiome of HIV patients versus uninfected individuals (31), it was found that Cladosporidium and Rothia had a positive correlation in uninfected individuals. We observed the same cooccurrence at T6 in treated leukemia patients. Also, in individuals uninfected with HIV, there was a negative correlation between Aspergillus and Prevotella (31) . We observed the same cross-domain connection, but as a cooccurrence. In the same study, the authors observed cocolonization of Cladosporidium and Pichia in uninfected individuals (31) , which is a fungal relationship that we detected at T3. Additionally, we saw a bacterial cooccurrence relationship between Prevotellaceae and Leptotrichia at T6, which has also been observed in subjects with apical periodontitis (32) and in subjects with halitosis, as well as in healthy individuals (33) .

Only one connection overlapped any two networks. Interestingly however, time point 3 showed more-dependent connectivity (multiple nodes linked to each other), while time points 3 and 6 had increased connectivity overall. As we had previously reported loss of bacterial diversity over time among the patients included in these analyses (2) , this suggests that either loss of bacterial diversity or changes in the oral environment under conditions of chemotherapeutic and antimicrobial pressure or both affect the bacterium-fungus interplay. This may potentially result in clinical consequences, as organisms that do not typically interact would inhabit the niche. This could potentially have acute consequences for functional metabolites, genetic exchange, and alterations in the virulence of the pathogen or in the immune status of the host (28, 34) .

Limitations. Some fungus-bacterium interactions in the oral cavity of diseased and healthy individuals (namely, Candida-bacterium interactions) have been intensely studied; however, interactions of bacteria with other known fungal inhabitants of the oral cavity have yet to be described. Thus, the impact of fungus-bacterium interactions in the oral cavity and on overall health is not completely understood. Although this exploratory study had a relatively small sample size and assessed only oral samples, the data shown here still highlight the dynamic nature of mycobiome and bacteriome communities and their interactions over the course of RIC. Unfortunately, the study design did not include information on diet, the metabolic milieu (e.g., other relevant comorbidities such as diabetes), oral health information (including oropharyngeal candidiasis), or occurrence of mucositis. These clinical factors/outcomes could potentially account for the heterogeneity seen within the data or might also be associated with specific mycobiome or bacteriome compositions/relationships. Despite not includ-ing a comprehensive record of clinical covariates, we were able to determine that specific changes in the mycobiome community structure appear to have clinical consequences in terms of bacterial infections. Moreover, interkingdom dynamics are likely driven by changes in the mycobiome due to chemotherapy intensity and subtype, as well as bacteriome changes driven by antibiotic exposures. Although this report highlights the underappreciated complex and dynamic nature of the cross-kingdom interactions that occur over time and result from the stress of chemotherapy and antibiotic therapy, these data are limited by our having used only marker gene analysis (16S rRNA gene and internal transcribed spacer 2 [ITS2] sequencing) to identify bacterium-fungus dynamics. Target gene analysis is limited in its ability to differentiate individual taxa, specifically to lower taxonomic levels such as the species level. Although the ability to assemble and differentiate fungal reads at the whole-genome level is currently restricted by the number of fungal genomes available in public databases, this information could be vastly improved in the future via shotgun metagenomics sequencing, particularly in regard to functional pathway analysis. An alternative approach would be to layer on metabolomics or metaproteomics to understand the functional consequences of bacterial and fungal changes as well as the interkingdom metabolic pathway interactions. Again, however, methods for data integration of shotgun metagenomic and metabolomics data, particularly in regard to cross-domain networks, are still in their infancy.

Conclusions. Here, we provide the first longitudinal analysis of the oral mycobiome in patients with acute leukemia combined with parallel oral bacteriome data. Incorporation of fungal analyses into microbiome studies may help improve mechanistic understanding of how commensal flora impact clinical outcomes. These data suggest that simultaneous mycobiome and bacteriome longitudinal analyses may add additional insight in order to fully understand how the microbiome impacts infectious risk in immunocompromised patients. We anticipate these data will inspire future studies aimed to further elucidate these relationships in the cancer treatment setting.

Patient sample collection, DNA extraction, sequencing, and bioinformatics. In this pilot exploratory study, we collected oral (buccal swab) samples (n ϭ 299) from 39 newly diagnosed adult AML patients undergoing IC at the MD Anderson Cancer Center (MDACC) in Houston, TX, from September 2013 to August 2015. These oral samples and accompanying 16S rRNA gene sequences were derived from a previously reported cohort (2) . The study protocol was approved by the MDACC Institutional Review Board (PA13-0339), and the study was conducted in compliance with the Declaration of Helsinki. Written informed consent was obtained from all participants before enrollment. Buccal swabs were collected by swabbing the inside the cheek 3 times on each side using a Catch-All sample collection swab (Epicentre), and the swab tips were then placed in sterile 2-ml cryovials and maintained at Ϫ80°C until DNA extraction. Collection of samples from each patient was performed at baseline, continued approximately every 96 h, and stopped upon neutrophil recovery. Clinical data, including antibiotic and antifungal administration, chemotherapy subtypes, and infection status, were extracted from the electronic medical records. Infections prior to neutrophil recovery were classified as microbiologically defined infections (MDI) or clinically defined infections (CDI) using established guidelines (35) . Patients were categorized as receiving (i) non-fludarabine-based high-intensity regimens, (ii) high-intensity fludarabinecontaining regimens, (iii) low-intensity hypomethylating-agent-based therapy, or (iv) other low-intensity/ low-dose cytarabine-based regimens.

Microbial DNA was extracted from oral swabs using the MO BIO PowerSoil DNA isolation kit protocol as described previously (36) . We performed ITS sequencing on a subset of patients (the first sequential 39 patients enrolled on study with available residual sample) in which we had previously characterized the microbiome using 16S rRNA gene sequencing (2, (37) (38) (39) . The internal transcribed spacer 2 (ITS2) region was amplified from oral swab DNA using primers ITS3 and ITS4 as described previously (36) . Each primer included an Illumina adapter and linker sequence, and the reverse primer (ITS4) also contained a unique 12-bp Golay barcode (40) . Amplicons were sequenced by the use of an Illumina MiSeq platform with the 2 ϫ 300-bp paired-end protocol. The ITS2 read pairs were demultiplexed based on the unique barcodes. Reads were merged and filtered using default settings in USEARCH v7.0.1090 (41) . ITS2 sequences were clustered into operational taxonomic units (OTUs) at a similarity cutoff of 97% using the UPARSE pipeline (42) . Chimeras were removed using USEARCH v8.0.1517 and UCHIME (43) . OTUs were aligned against sequences from the NCBI GenBank Plant and UNITE databases (44) . Abundance data were recovered by mapping the demultiplexed reads to the UPARSE OTUs. A custom script developed at the Center for Metagenomics and Microbiome Research at Baylor College of Medicine was used to construct an OTU table. Unmapped (Ͻ80% identity or Ͻ95% coverage) OTUs were manually analyzed by BLASTN (37) . After we performed filtering that included only the samples which had corresponding clinical and 16S data to match the ITS data, 290 samples were obtained for further analyses.

Statistical methods. For ␤-diversity, principal-coordinate analysis (PCoA) plots were produced using Bray-Curtis distances. A mixed-effects model was applied to determine if time was a significant predictor of ␣-diversity. A random intercept for patients was included. A Kruskal-Wallis test was performed to compare ␣-diversity and taxon abundance data between chemotherapy subtype groups. In association testing, we focused on the top five most abundant fungal genera based on average abundance: Candida, Cladosporium, Fusarium, Malassezia, and Saccharomyces. Similarly, a Mann-Whitney test was performed to compare both alpha-diversity data and taxon abundance data corresponding to chemotherapy intensity (high versus low), use of amphotericin B (yes versus no), and infection status (yes versus no). The linear discriminant analysis effect size (LEfSe) algorithm was used to test whether taxon abundance predicted infection status. The LEfSe algorithm was run using the Galaxy module with default settings. Mixed-effects models were used to determine whether taxon abundance significantly changed over time. Analyses based on these models were performed on the centered log ratio (clr) transformed relative abundance data to address the compositional characteristics of the data. Each model contained the following variables: time point, the variable of interest, and an interaction between the two. To account for repeated measures, a random intercept for patient was included in the model. An analysis of variance (ANOVA) F-test was applied to determine whether time and the variable of interest were significantly associated. The statistical analyses performed as a part of this study were exploratory in nature. In general, P values of less than 0.05 were considered significant; for analyses across the mycobiome, we computed false-discovery-rate (FDR)-adjusted P values, which are provided in the supplemental tables. Data analysis and visualization were performed using R, mainly using functions from the R packages ggplot2 and phyloseq (26, 45) .

Cross-domain association networks. SPIEC-EASI was used to construct the cross-domain microbial networks at the genus level for both bacteria and fungi. For the tables listing bacterial and fungal genera, the data corresponding to time points 1 (T1), 3 (T3), and 6 (T6) were separated (46) . As statistically sound methods for longitudinal analyses/trends are underdeveloped for microbiome network analyses, we aimed at analyzing samples at baseline prior to chemotherapy and antimicrobial administration (T1), at a time point after chemotherapy administration (T3), and at a time point after antimicrobial administration (T6). As we had significant sample dropout after T6, and given our already small sample size, T6 was the last time point which retained most samples. In addition, all patients who received broadspectrum antibiotics had received at least one type of broad-spectrum antibiotic by T6. For each OTU table, taxa are included in the network analysis if the mean abundance was greater than 0.01. T1 had 35 samples with 16 bacterial taxa and 13 fungal taxa. T3 was represented by 36 samples, with 22 bacterial taxa and 16 fungal taxa. T6 contained 31 samples with 21 bacterial taxa and 15 fungal taxa. The sparsity parameter was set to 0.3, and the value corresponding to was set to 50. We verified the stability of our results for a single value of . We also examined sensitivity across multiple values and found that the networks were not sensitive to the choice of .

Data availability. All 16S rRNA gene sequences were deposited in the NCBI Sequence Read Archive (https://www.ncbi.nlm.nih.gov/sra) under BioProject identifiers (IDs) PRJNA352060 and PRJNA526551, whereas ITS2 sequences were submitted under the BioProject ID PRJNA599151 (the authors could not make this BioProject record available at the time of this paper's publication due to circumstances related to the COVID-19 pandemic, but it will be made accessible as soon as possible after publication).

Supplemental material is available online only. Any opinions, findings, and conclusions or recommendations expressed in this material are ours and do not necessarily reflect the views of the National Science Foundation.

We have no relevant disclosures with respect to conflicts of interest.

Gut microbiota injury in allogeneic haematopoietic stem cell transplantation

Gut microbiome signatures are predictive of infectious risk following induction therapy for acute myeloid leukemia

Gut microbiome composition predicts infection risk during chemotherapy in children with acute lymphoblastic leukemia

Intestinal microbiota and relapse after hematopoietic-cell transplantation

The effects of intestinal tract bacterial diversity on mortality following allogeneic hematopoietic stem cell transplantation

The mycobiome: impact on health and disease states

Role of the mycobiome in human acute graftversus-host disease

Implementation of a pan-genomic approach to investigate holobiont-infecting microbe interaction: a case report of a leukemic patient with invasive mucormycosis

Integrated analysis of clinical and microbiome risk factors associated with the development of oral candidiasis during cancer chemotherapy

Gastrointestinal microbiome and mycobiome changes during autologous transplantation for multiple myeloma: results of a prospective pilot study

Chemotherapy-driven dysbiosis in the intestinal microbiome

Direct effects of nonantifungal agents used in cancer chemotherapy and organ transplantation on the development and virulence of Candida and Aspergillus species

Chemotherapy-induced oral mucositis is associated with detrimental bacterial dysbiosis

Candida albicans induces mucosal bacterial dysbiosis that promotes invasive infection

Antibiotic exposure as a risk factor for fluconazole-resistant Candida bloodstream infection

The effect of antibacterial and non-antibacterial compounds alone or associated with antifungals upon fungi

Candida spp. colonization significance in critically ill medical patients: a prospective study

Effects of carbapenems and their combination with amikacin on murine gut colonisation by Candida albicans

The human mycobiome in health and disease

Biodiversity of the human oral mycobiome in health and disease

Pneumonia during remission induction chemotherapy in patients with acute leukemia

Redefining the human oral mycobiome with improved practices in amplicon-based taxonomy: discovery of Malassezia as a prominent commensal

Characterization of the oral fungal microbiome (mycobiome) in healthy individuals

The epidemiology of fungal infections in patients with hematologic malignancies: the SEIFEM-2004 study

Interindividual variability and intraindividual stability of oral fungal microbiota over time

Bioconductor workflow for microbiome data analysis: from raw reads to community analyses

In vitro resistance of clinical Fusarium species to amphotericin B and voriconazole using the EUCAST antifungal susceptibility method

Fungal-bacterial interactions and their relevance in health

Identification of Malassezia species from patient skin scales by PCR-RFLP

The prevalence and species composition of Malassezia yeasts in patients with clinically suspected onychomycosis

Oral mycobiome analysis of HIV-infected patients: identification of Pichia as an antagonist of opportunistic fungi

Apical periodontitis: preliminary assessment of microbiota by 16S rRNA high throughput amplicon target sequencing

Supragingival plaque microbial community analysis of children with halitosis

Fungal-bacterial interactions in health and disease

Clinical practice guideline for the use of antimicrobial agents in neutropenic patients with cancer: 2010 update by the Infectious Diseases Society of America

The gut mycobiome of the Human Microbiome Project healthy cohort

Basic local alignment search tool

The role of the gastrointestinal microbiome in infectious complications during induction chemotherapy for acute myeloid leukemia

Characterization of oral and gut microbiome temporal variability in hospitalized cancer patients

Ultra-high-throughput microbial community analysis on the Illumina HiSeq and MiSeq platforms

Search and clustering orders of magnitude faster than BLAST

UPARSE: highly accurate OTU sequences from microbial amplicon reads

UCHIME improves sensitivity and speed of chimera detection

UNITE: a database providing Web-based methods for the molecular identification of ectomycorrhizal fungi

phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data

Sparse and compositionally robust inference of microbial ecological networks