|Home | About | Journals | Submit | Contact Us | Français|
Bone and joint infections (BJI) are severe infections that require a tailored and protracted antibiotic treatment. Yet, the diagnostic based on culturing samples lacks sensitivity, especially for hardly culturable bacteria. Metagenomic sequencing could potentially address those limitations. Here, we assessed the performances of metagenomic sequencing on 24 BJI samples for the identification of pathogens and the prediction of their antibiotic susceptibility. For monomicrobial samples in culture (n=8), the presence of the pathogen was confirmed by metagenomics in all cases. For polymicrobial samples (n=16), 32/55 bacteria (58.2%) were found at the species level (and 41/55 [74.5%] at the genus level). Conversely, 273 bacteria not found in culture were identified, 182 being possible pathogens and 91 contaminants. A correct antibiotic susceptibility could be inferred in 94.1% and 76.5% cases for monomicrobial and polymicrobial samples, respectively. Altogether, we found that clinical metagenomics applied to BJI samples is a potential tool to support conventional culture.
Bone and joint infections (BJI) are severe infections that affect a growing number of patients1. Along with the surgical intervention, the microbiological diagnosis is a keystone of the management of BJI in (i) identifying the bacteria causing the infection and (ii) assessing their susceptibility to antibiotics. Currently, this is achieved by culturing surgical samples on various media and conditions, together with a long time of incubation to recover fastidiously-growing bacteria that can be involved in BJI. Still, some bacteria will not grow under these conditions because of extreme oxygen sensitivity, a prior antibiotic intake or metabolic issues (e.g. quiescent bacteria in chronic infections). Consequently, the antibiotic treatment may not span all the bacteria involved in the infection, which can favour the relapse and the need for a new surgery.
Clinical metagenomics refers to the concept of sequencing all the DNA (i.e. all the genomes) present in a clinical sample with the purpose of identifying pathogens and inferring their antibiotic susceptibility pattern2. This new, culture-independent method takes advantages of the thrilling development of next-generation sequencing (NGS) technologies since the mid-2000s. The NGS platforms typically yield thousands to millions of reads (sequences of size ranging from 100bp to a few kbp), which virtually enables to recover the sequences of all the genes present in the sample, yet in a disorganised fashion. Substantial bio-informatics efforts are thereby needed to re-construct and re-order the original sequences in genomes, and are referred to as the assembly process. Hence, various information such as the taxonomic identification of the present species, antibiotic resistance determinants (ARDs), mutations (as compared to a reference genome or sequence), single nucleotide variants (SNVs, for clonality assessment) and virulence genes can be determined.
Clinical metagenomics is an emerging field in medicine. So far, a few attempts to use metagenomics on clinical samples have been performed (on urines3, 4, cerebrospinal fluid or brain biopsy5, 6, blood7 and skin granuloma8) likely because of the high price of metagenomics and the complexity of the management of sequence data for clinical microbiologists. To the best of our knowledge, metagenomics has never been applied to BJI samples.
As for the inference of antibiotic susceptibility testing from the genomic information, a few studies focusing on Escherichia coli, Klebsiella pneumoniae, Mycobacterium tuberculosis and Staphylococcus aureus have constantly shown excellent correlations between the analysis of the genomic content of antibiotic resistance determinants (ARDs) and the phenotype9–15 while performances were not as good for Pseudomonas aeruginosa 16. Furthermore in metagenomic data, the possible presence of multiple pathogens raises the issue of linking ARDs to their original host in order to infer its antibiotic susceptibility pattern3. So far no method has been proposed to address this question.
Applying metagenomics in the context of BJI is thus seducing in that (i) there is no limit in the number of species and ARDs that could be detected (as opposed to PCR-based methods which detect targeted ARDs), (ii) unculturable bacteria, fastidious growers (such as Propionibacterium sp.) or bacteria altered by prior antibiotic use would be recovered, and (iii) the antibiotic susceptibility inference would benefit from both the detection of ARDs (such as mecA, qnr, dfr, erm, etc.) and the identification of mutations leading to resistance to key antibiotics used in BJI. Here, our main objective is to assess the performances of clinical metagenomics in BJI in terms of pathogen identification and inference of antibiotic susceptibility, as compared with conventional microbiology (gold standard).
We first extracted the samples for which the quantity of material exceeded or was equal to 1mL (n=77). We recovered more than 1pg/µL bacterial DNA mostly for samples that had grown at>100 CFUs (Extended Data Fig. 1, panel A), while the concentration of human DNA did not correlate with the bacterial load (Extended Data Fig. 1, panel B). The remaining samples (<1mL material), that had grown at least 100 CFUs (n=25, see Extended Data), were submitted to extraction. In total, the DNA of 104 samples was extracted, from which 24 met our internal requirements for being sequenced (i.e. contained at least 1pg/µL bacterial DNA, and less than 99% human DNA, Supplementary Table 1). Other samples could not be sequenced because of the low bacterial DNA concentration and/or a high human DNA concentration (Supplementary Table 1). These 24 samples were obtained from 14 patients (Table 1). All throughout the manuscript, we will refer as monomicrobial (n=8) and polymicrobial (n=16) samples which respectively yielded one and more than one bacterial species in culture.
The bioinformatic pipeline is depicted in the Fig. 1. After trimming, we obtained a mean number of 20,092,168 paired reads per sample (range 8,256,850–29,099,374, Supplementary Table 2). With the Kraken classifier, the mean rate of classified reads (as bacteria, archea or virus) was 27.9% (range 1.8–85.7, Supplementary Table 2). Of note, the classification rate was correlated to the proportion of bacterial DNA as found by qPCR (Pearson’s correlation test, p<0.001, cor=0.70, Extended Data Fig. 2). The assembly of the classified reads with metaSPAdes yielded a mean number of contigs of 10,444 (range 3,087–18,513, Supplementary Table 2), for a mean total number of base pairs of 8.3M (range 2.9M–16.5M, Supplementary Table 2). The total number of base pairs of contigs was higher in polymicrobial samples than in the monomicrobial ones (respectively 9.7M vs. 5.5M, t test p<0.05, Extended Data Fig. 3). The mean size of the contigs was 805bp (median 369bp, maximum 445,300bp, Extended Data Fig. 4 panel A). It was higher in polymicrobial samples than in monomicrobial samples (respectively 865bp vs 654bp, p<0.001, Extended Data Fig. 4 panel B). Besides, the genome coverage of the bacteria found in culture was higher in monomicrobial samples than in polymicrobial samples (respectively 78.7% [98.8% when not considering the sample 46 which is actually a polymicrobial infection according to metagenomic findings] vs. 22.6%, p<0.001, Extended Data Fig. 3 panel F).
In monomicrobial samples (n=8, Table 2), the sensitivity of metagenomic sequencing was 100% (8/8) at both genus and species levels. The pathogens identified by culture were found by MetaPhlAn2 at very high relative abundances (over 94.6%) with the exception of sample 46 in which Streptococcus anginosus was found at a relative abundance of 2.2% (Fig. 2). In polymicrobial samples (n=16, Table 2 ), the sensitivity of metagenomic sequencing was 58.2% (32/55) at the species level. At the genus level, the sensitivity increased to 74.5% (41/55). Besides, metagenomic sequencing could identify all bacteria found by culture in a given sample in 11/24 (45.8%) samples at the species level, including 3/16 (18.8%) for polymicrobial infections. At the genus level, 15/24 (62.5%) samples were in agreement with cultures, including 7/16 (43.8%) samples with polymicrobial infections.
Apart from the bacteria that were found in culture (n=63) in the 24 positive samples, a total of 273 bacteria, not found in culture, were identified by MetaPhlAn2 (Fig. 3). From the metagenomic sequencing of two negative controls (duplicates), we identified 10 bacterial species (Extended Data Fig. 5). In both negative extraction controls, Propionibacerium acnes was the most abundant species identified. Consistently, P. acnes was found in 20/24 clinical samples (Extended Data Figs 6 and 7). Moreover, we observed that the relative abundance of P. acnes in samples was negatively correlated to their total DNA concentration (Extended Data Fig. 6), in consistence with P. acnes being a contaminant in this study17, 18. For other species, such correlation could not be tested because of their low occurrence in samples. Beyond the 10 contaminants identified in the negative controls (found in 29 occurrences), we identified 23 putative contaminants found in 37 occurrences (Fig. 3, Supplementary Table 3), some already reported as reagents contaminants19. Others were more unexpected such as Borrelia sp. (samples 103, 104, 108 and 110) or Rickettsia japonica (sample 117). Still, the taxonomic assignment of the contigs did not confirm the presence of those species and manual Blastn of reads against the NCBI nr database supported that they were likely in silico contaminants19, 20. In summary, the mean percentage of reads assigned to bacteria that were considered as contaminants was 0.7% (range 0.01–5.4%). Besides, we identified 25 species that could be due to a misclassification of reads to closely related bacteria, such as in samples 184 (where Corynebacterium striatum was found in culture, and some metagenomic reads were identified as Dermabacter sp. and Corynebacterium pyruviciproducens), samples from patient C (19, 103, 104 where Acinetobacter baumannii and Achromobacter xylosoxidans were found in culture, and some reads were identified as from other Acinetobacter spp., Achromobacter spp., or Advenella kashmirensis, a bacterium close to Achromobacter) (Supplementary Table 3). Hence, a total of 183 bacteria not recovered in culture and not acknowledged as contaminants were identified in metagenomic sequencing. For one sample that was monomicrobial in culture (sample 46, that yielded S. anginosus), 38 other species were identified by metagenomics. Interestingly, these species appeared to be commonly found in the oropharyngeal microbiota, which was consistent with the site of the infection (mandible). In polymicrobial samples such as samples 4 and 140 (patient B), 90 and 158 (patient H), 108 and 181 (patient I), several anaerobic bacteria were identified (range 3–40, see Supplementary Table 3) in consistence with the sporadic isolation of anaerobic bacteria in culture (Table 2 and supplementary Table 3). In both samples 4 and 140 from patient B, the most abundant species was Propionibacterium propionicum (respective abundances of 71.5% and 43.2%) that was not found in culture. Arguments in contradiction with P. propionicum being a contaminant in these samples are that P. propionicum was not found in the negative controls, that the only Propionibacterium species found in other samples was P. acnes, and that the abundance of P. propionicum was high (Extended Data Fig. 7) whereas the abundance of P. acnes was low in the samples where it was identified.
Taking advantage of the depth of analysis of metagenomic sequencing, we addressed whether within the dominant species identified in the samples, more than one clone could be identified (see Extended Data). We assumed that in case of multiple clones within one species, the single nucleotide variants (SNVs) would be homogeneously distributed along the contigs (Extended Data Fig. 8, panels A–B) unlike when only one clone would be present (Extended Data Fig. 8, panels C-D). Accordingly, we found polyclonal populations for 29 of the 74 (39.2%) bacterial species that were tested. Among the bacteria that were found in culture and that were tested for polyclonal populations (n=32), 8 (25%) displayed a polyclonal population: Morganella morganii (samples 4 and 140), Streptococcus agalactiae (samples 103 and 117), Staphylococcus aureus (samples 28 and 110), S. anginosus (sample 158) and Pseudomonas aeruginosa (sample 128). Moreover, we observed for M. morganii (samples 4 and 140) no mutations on the topoisomerases were found in the sample 140, while in sample 4 the Ser83Ile and Ser84Ile were found the in GyrA and ParC, respectively. This suggests that one population of M. morganii was susceptible to fluoroquinolones and the other was not. In culture though, only the fluoroquinolone resistant clone was found.
A total of 151 ARDs (61 unique) were identified from the 24 samples (range 2–22, Table 3). The most frequent ARD families were beta-lactamases (n=30), Tet(M) (n=26), Erm (n=18) and Dfr (n=16). For monomicrobial samples, we assumed that the ARDs identified by metagenomics were expressed by the bacterium that was recovered in culture. Considering together (i) the antibiotic class the ARDs usually confer resistance to, (ii) the intrinsic antibiotic susceptibility profile of the species and (iii) the analysis of the sequence of specific genes (gyrA, parC, rpoB), we could infer a in silico susceptibility in agreement with the phenotypic susceptibility in 94.1% (111/118) cases (Fig. 4 and Supplementary Table 3), a case being defined as the susceptibility testing of one antibiotic for one sample. Of note, the six major errors (overprediction of resistance as compared to culture) originated from sample 46 where anaerobic bacteria and likely associated ARDs were found in metagenomic sequencing only. We then attempted to link the ARDs with their host through their respective depth of sequencing. Our original hypothesis was that resistant bacteria would carry at least one ARD-encoding gene copy per genome. As a consequence, one given ARD could not be sequenced less than the median depth of sequencing per contig of its host. Hence, we plotted the respective median depth of sequencing of ARDs and bacterial species (Extended Data Fig. 9). In contradiction with our hypothesis, the analysis of monomicrobial samples (in which ARDs were assumed to be carried by the identified bacteria) showed that the depth of sequencing of the contigs assigned to the pathogen could be higher than that of the ARDs (Extended Data Fig. 9). Accordingly for polymicrobial samples, we did not attempt to link ARD and their host and we separately considered the ARDs and the bacteria found in the sample (Supplementary Table 3). This way, we inferred a correct susceptibility in 76.5% (192/251) cases. Very major errors mostly occurred because some bacteria with specific resistance patterns were not detected in sequencing (Supplementary Table 3 ). Conversely and along with the observations with monomicrobial samples, most major errors occurred because some bacteria and ARDs were found in sequencing but not in culture. Of note, the prediction of susceptibility to fluoroquinolones (a pivotal antibiotic families for BJI treatment) was correct in 100% (24/24) samples.
In order to assess the potential clinical impact of the identification of bacteria by metagenomic sequencing that were not identified in culture, we retrospectively reviewed the antibiotic treatments received by the 14 patients in this study and the clinical outcome of their infection. A total of 60 different antibiotic regimens were administered, among which 6 (10.0%) were putatively not being active against bacteria found in metagenomic sequencing only (Supplementary Table 4). However as of January 2017 (median follow-up was 36 weeks, range 16–66 weeks), only one relapse was observed that involved Enterococcus faecalis. In this patient, the definitive treatment included amoxicillin, which is not active against some bacteria (Supplementary Table 4) found in metagenomic sequencing only. Still, the connection with the relapse remains speculative.
We ran the same pipeline analysis onto the 24 samples downsized at 2M paired reads. We observed that the taxonomic distribution did not apparently change for the most abundant species (Extended Data Fig. 10), but the mean genome coverage of the pathogen(s) (found in culture) was slightly lower in the downsized group than in the full-reads group (respectively 25.0% vs. 30.2%, Student paired test p<0.001, Extended Data Fig. 11). Also, only 86 ARDs were found after downsizing while 151 were detected before (Student paired test p<0.001, Extended Data Fig. 11). Of note, the impact of downsizing was observed in both monomicrobial and polymicrobial samples (Extended Data Figs 12 and 13).
The main result of this study is that we showed that metagenomic sequencing could be a potential tool in the diagnostic of BJI. Indeed for monomicrobial infections, the pathogen was identified in 100% (8/8) samples and the antibiotic susceptibility prediction was successful in 94.1% (111/128) cases. In case of polymicrobial samples, the high abundance of several bacteria (mostly anaerobes) did occasionally prevent from the correct identification of the pathogens and their antibiotic susceptibility profiles. Accordingly, our findings support that currently, metagenomic sequencing of BJI samples could not replace conventional methods based on culture due to the limitations encountered when several bacterial species are present in the samples, but rather be performed in support.
Interestingly, metagenomic sequencing yielded in some ways more information than culture. First, metagenomic sequencing identified many more bacterial species than culture. Besides likely contaminants, some bacteria that were not detected by culture were probably true positive and may not have been targeted by the selected antibiotic regimen, as we observed in 10% cases. Second, we could identify multiple clonal populations within some species, which could differ in their susceptibility to antibiotics as we observed for fluoroquinolones in M. morganii. Sequencing multiple clones obtained from the culture of BJI samples would validate this finding and shall be considered for further studies. In all, using metagenomic data could help to tailor the antibiotic regimen for the treatment of BJI, and the added-value of clinical metagenomics in BJI should now be assessed.
However, there are several hindrances to the application of metagenomic sequencing to BJI samples. First, we could only sequence 24 out of 179 samples, due to a low amount of bacterial DNA recovered from the samples. This is the main limitation of this study as it reduced the diversity of clinical situations that we could address. Nonetheless, the samples from this study have been frozen and thawed, which decays some bacteria and releases DNA. As the DNA extraction method we used eliminates free DNA after lysing eukaryotic cells, it is likely that we could have sequenced more samples if they would not have been frozen. This said, recovering enough bacterial DNA (in terms of quantity and proportion with respect to human DNA) remains challenging. Also, the high costs of NGS (at least several hundreds of USD per sample) together with an unascertained clinical impact are obstacles to its reimbursement by health agencies. We tested the impact of a lower depth of sequencing (that could be achieved at a lower cost per sample) and showed that despite the taxonomic profiles of the bacterial populations were similar, the inference of antibiotic susceptibility was less accurate due to a lower recovery of genes involved in antibiotic resistance. Our results suggest that clinical metagenomics should indeed benefit from the highest depth of sequencing, despite the high cost. Eventually, identifying contaminants remains challenging. We used laboratory negative controls that did not includes all the putative contaminants that were identified in samples, suggesting that contamination may also occur during the sampling process. In this perspective, a clinical metagenomics negative control should be taken at the sampling stage.
Besides, our observations suggest that clinical metagenomics will soon require, as for clinical microbiology, a specific expertise combining clinical, biological and bioinformatic skills in order to infer clinically relevant results from metagenomic data. In this perspective, the development of clinical metagenomics will need the definition of quality standards, e.g. what is the sufficient genome coverage for a given bacterium to consider that its antibiotic susceptibility profile can be likely inferred. In the long term, algorithms should be built to provide clinicians with clear data and robust algorithms to support clinical decisions.
In conclusion, we showed that metagenomic sequencing of BJI samples was a potential tool to support conventional methods. In this perspective, its main limitations (DNA extraction, cost and data management) should be tackled, and the clinical benefit provided by clinical metagenomics should now be assessed in a prospective fashion.
We initially included 179 per-operative samples recovered from 47 patients (range 1–8 samples per patient, Supplementary Table 1). All but 2 (swabs) were solid specimens. The quantity of material for each non-swab sample (n=177) was macroscopically estimated: less than 1mL (n=100), from 1 to 10mL (n=60) and more than 10mL (n=17). The samples were collected from September 2015 to January 2016 in the orthopedic departments of the CRIOAc (Regional Reference Center for Complex Osteo-Orticular Infections), Lyon, France (https://www.crioac-lyon.fr) and stored at −80°C until shipment in dry ice to the Genomic Research Laboratory in Geneva on April 13, 2016. The samples had previously been cultured on (i) one sheep blood agar plate (bioMérieux, Marcy l’Etoile, France) incubated under aerobic conditions at 36±1°C incubated day 1 and 2 for microbial growth, then discarded; (ii) two chocolate blood agar plates (bioMérieux, Marcy l′Etoile, France) kept under a 5% CO2-enriched aerobic atmosphere at 36±1°C: one plate was observed on day 1 and 2 then discarded; the second plate was conserved opened and read only after ten days; (iii) two Schaedler agar plates (bioMérieux, Marcy l′Etoile, France) incubated in anaerobic conditions: one plate was observed on day 2 and 5 then discarded; the second was opened only after 15 days and read; (iv) one Schaedler broth supplemented with Vitamin K1 (BD Diagnostic Systems, Le Pont-de-Claix, France) was kept 15 days in anaerobiosis and examined every day in search of a blur; if positive, Gram staining was performed and the broth was sub-cultured on agar plates (i.e., one chocolate agar and one Schaedler agar in anaerobic condition) incubated for 3 days; subculture of the broth was systematically done on day 15 even in the absence of a blur. When cultures were positive, each microorganism was identified using matrix-assisted laser desorption ionization-time of flight mass spectrometry (MALDI-TOF® MS, Bruker Daltonics, Bremen, Germany) and antimicrobial susceptibility was tested according to CA-SFM recommendations (Vitek, bioMérieux, Marcy l′Etoile, France). The quantities of bacteria were visually counted up to 100 colony-forming units (CFUs) then beyond expressed as>102. As true numbers were required to determine the proportion of bacteria in culture, we considered the value of 101 when>102 CFUs were counted on the plate. No negative sample in culture was included in this study. A single bacteria or yeast was recovered for 104 out of 179 samples (58.1%), the remaining yielding 2 (24/179, 13.4%), 3 (26/179, 14.5%), 4 (14/179, 7.8%) or 5 (11/179, 6.1%) bacteria and yeasts. Eventually, we sequenced two negative controls (the same process as samples but without any biological material). This study involved already existing, anonymous samples for which a further use was authorized by the Ethical Committee of the Lyon University Hospital (September 25, 2014). According to the French guidelines, as the exploitation of the samples and associated data is not performed in an interventional way, the consent of the patient is not needed.
Tissue samples were cut into small pieces on a disposable Petri dish support using a scalpel while the swabs were thoroughly vortexed in 1mL physiological solution. DNA was extracted from 50–100mm3 shredded sample using the Ultra-Deep Microbiome Prep kit (Molzym, Bremen, Germany) according to the manufacturer’s instructions (Version 2.0) for tissue samples. This method aims at decreasing the eukaryotes DNA through a differential lysis of human cells, followed by the elimination of free DNA. The concentration of bacterial and human DNA was determined by qPCR experiments as described previously21, using 16S rRNA and beta-actin reference genes, respectively. The reference curves for bacterial and human DNA quantitation were generated using known concentrations of E. coli DH5α genomic DNA and human genomic DNA from the TaqMan beta-Actin Detection Reagent kit (Applied Biosystems, Framingham, MA), respectively. About 3 ng of DNA (the sum of bacterial and human DNA determined by qPCR) in a 50µL volume were sent to Fasteris (Plan-les-Ouates, Switzerland) for DNA purification (resulting in a 10µL purified solution) and subsequent sequencing. The Nextera XT DNA Sample Preparation Kit was used according to the Illumina (San Diego, CA) instructions except that 16 (instead of 12) PCR enrichment cycles were used. For samples with less than 3 ng total DNA in 50µL of non-purified extracts, the volume of Nextera tagmentase was reduced from 5 to 3µL. The libraries were sequenced in Rapid Run mode for 2×250+8 cycles on an Illumina HiSeq 2500 instrument (with a HiSeq Rapid Flow Cell v2) at Fasteris using HiSeq Rapid SBS Kit v2 and HiSeq Control Software 2.2.58. Demultiplexed fastq files were generated with CASAVA-1.8.2 from on-instrument base-calling by Real-Time Analysis (RTA) software 220.127.116.11. The Trimmomatic package22 was used to remove bases that correspond to the standard Illumina adapters.
The quality of the reads was assessed by FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). The reads (fastq format) were processed with MetaPhlAn2 to get the taxonomic profile of the microbial community23, 24. MetaPhlAn2 uses a database of clade-specific marker genes, and normalizes the abundances of the species with respect to the genome size of the bacteria to which reads are assigned. In order to filter out the human reads, we used Kraken25 with default parameters and the –classified-out option on the miniKraken database, that embeds the genomes of bacteria, archaea and viruses, but noticeably not those of eukaryotes. Together with the seqtk subset command (https://github.com/lh3/seqtk), we were able to recover the pairs of reads that were found as from bacteria, archaea and viruses, and to assemble them using metaSPAdes26 with default parameters. The resulting contigs were annotated using PROKKA27. The genome coverage was assessed by processing the metagenomic contigs by QUAST28 against the reference genome of the species downloaded from the RefSeq database (http://www.ncbi.nlm.nih.gov/genome/). The total quality-filtered reads (fastq format) were mapped using Bowtie229 (using the–local argument for accepting soft clipping) onto a database made of the ARDs from the ResFinder database30 (downloaded in July 2016) and ARDs from functional metagenomic studies31–33, eventually clustered at a 95% nucleic acid identity level (using CD-HIT34). The mapped reads were then assembled with metaSPAdes with default parameters. The open reading frames (ORFs) and the amino acid sequences were obtained by Prodigal35. The identification of ARDs was performed by blastp36 (with a 10−30 e-value) using the aforementioned ARD database, using a 80% amino acid identity threshold over 80% of the reference ARD sequence. To get the depth of sequencing of the bacterial species and of the ARDs in samples, we separately mapped the reads using Bowtie2 against the contigs assigned to one given species and against the ARDs identified in this sample. The depth of sequencing (expressed in×, i.e. the number of time each nucleotide was sequenced) was calculated using Samtools37 and BBMap (https://sourceforge.net/projects/bbmap). The single nucleotide variants (SNVs) were identified using Samtools, Freebayes38 and Tablet39.
The high cost for metagenomic sequencing is a major hindrance to its application in the clinical setting. To address this issue, multiplexing samples offers the possibility to decrease the cost but with a lower number of reads per sample. To assess the possibility of using a lower throughput for our samples, we downsized our samples to 2M quality-filtered reads and applied the same pipeline as for full-reads samples.The same pipeline was applied to all samples after downsizing to 2M paired reads using the seqtk sample command (https://github.com/lh3/seqtk). The figures were drawn using R v3.2.3 and the ggplot2 package40, with colors occasionally chosen from colorbrewer2.org. The raw reads (just adaptor-trimmed) are available under the Bioproject PRJNA382079.
Besides the identification of species that were not recovered in culture, we tested whether more than one clone could be identified within the species. In this perspective, we considered the contigs>1000bp and assigned them a taxonomy with Kraken. For the main species identified in the sample (species with>0.1% total bacterial reads with Kraken and that were not considered contaminants (see section above), with at least 10 contigs available, n=74), we mapped the reads against their respective contigs and analyzed the distribution of single nucleotide variants (SNVs). For all species, we found some regions containing SNVs, but we assumed that whether more than one clone would be present within one species, the number of SNVs per contigs would be homogeneously distributed on the contigs, yielding a positive correlation between the number of SNVs per contigs and their size (Extended Data Fig. 8).
T.F., F.L. and W.M. collected the samples and the associated data. M.G. performed the DNA extractions. E.R., V.L and J.S. performed the analysis. E.R. wrote the manuscript.
The authors declare that they have no competing interests.
Electronic supplementary material
Supplementary information accompanies this paper at doi:10.1038/s41598-017-07546-5
Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.