Louca et al. Microbiome (2018) 6:41 https://doi.org/10.1186/s40168-018-0420-9SHORT REPORT Open AccessCorrecting for 16S rRNA gene copynumbers in microbiome surveys remains anunsolved problemStilianos Louca1,2* , Michael Doebeli1,2,3 and Laura Wegener Parfrey1,2,4AbstractThe 16S ribosomal RNA gene is the most widely used marker gene in microbial ecology. Counts of 16S sequencevariants, often in PCR amplicons, are used to estimate proportions of bacterial and archaeal taxa in microbialcommunities. Because different organisms contain different 16S gene copy numbers (GCNs), sequence variant countsare biased towards clades with greater GCNs. Several tools have recently been developed for predicting GCNs usingphylogenetic methods and based on sequenced genomes, in order to correct for these biases. However, the accuracyof those predictions has not been independently assessed. Here, we systematically evaluate the predictability of 16SGCNs across bacterial and archaeal clades, based on ∼ 6,800 public sequenced genomes and using severalphylogenetic methods. Further, we assess the accuracy of GCNs predicted by three recently published tools (PICRUSt,CopyRighter, and PAPRICA) over a wide range of taxa and for 635 microbial communities from varied environments.We find that regardless of the phylogenetic method tested, 16S GCNs could only be accurately predicted for a limitedfraction of taxa, namely taxa with closely to moderately related representatives (15% divergence in the 16S rRNAgene). Consistent with this observation, we find that all considered tools exhibit low predictive accuracy whenevaluated against completely sequenced genomes, in some cases explaining less than 10% of the variance.Substantial disagreement was also observed between tools(R2 < 0.5)for the majority of tested microbial communities.The nearest sequenced taxon index (NSTI) of microbial communities, i.e., the average distance to a sequencedgenome, was a strong predictor for the agreement between GCN prediction tools on non-animal-associated samples,but only a moderate predictor for animal-associated samples. We recommend against correcting for 16S GCNs inmicrobiome surveys by default, unless OTUs are sufficiently closely related to sequenced genomes or unless a needfor true OTU proportions warrants the additional noise introduced, so that community profiles remain interpretableand comparable between studies.Keywords: 16S rRNA, Gene copy number, Microbiome surveys, Phylogenetic reconstructionIntroductionAmplicon sequencing of the 16S ribosomal RNA (rRNA)gene is widely used for estimating the composition ofbacterial and archaeal communities. Global microbialdiversity initiatives, including the Human MicrobiomeProject [1], the EarthMicrobiome Project [2], and the TaraOceans global ocean survey [3], use the 16S rRNA geneto determine which microbes are present by matching 16S*Correspondence: stilianos.louca@googlemail.com1Biodiversity Research Centre, University of British Columbia, Vancouver,Canada2Department of Zoology, University of British Columbia, Vancouver, CanadaFull list of author information is available at the end of the articlerRNA sequence variants to reference databases like SILVA[4] and estimate the proportions of taxa based on rel-ative read counts. Many bacteria and archaea, however,have more than one copy of the 16S gene, which leads tobiased cell count estimates when the latter are estimatedsolely based on 16S rRNA read counts [5]. This has led toefforts to predict the distribution of 16S gene copy num-bers (GCNs) across clades based on available sequencedgenomes, in order to then correct 16S rRNA read countsto account for variable 16SGCNs in cells [5–8]. These cor-rections can substantially affect community profiles anddiversity patterns, since some clades have over 10 copies ofthe 16S rRNA gene [5, 7]. It is thus important to carefully© The Author(s). 2018 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, andreproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to theCreative 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.Louca et al. Microbiome (2018) 6:41 Page 2 of 12evaluate the accuracy [9] of predicted 16S GCNs acrossthe wide range of microbial clades encountered in micro-biome surveys. Inaccurate prediction of 16S GCNs canintroduce substantial noise to community profiles, whichcan be worse than the original GCN-related biases, par-ticularly when prediction methods differ between studies.An accurate prediction of 16S GCNs relies heavilyon the assumption that 16S GCNs are sufficiently phy-logenetically conserved. That is, 16S GCNs must beautocorrelated among related taxa at least across phylo-genetic distances typically covered by available sequencedgenomes [10]. Kembel et al. [5] found that 16S GCNexhibits a strong phylogenetic signal, as measured byBlomberg’s K statistic [11], and concluded that 16S GCNsmay be predictable based on phylogenetic placement withrespect to genomes with known 16S GCN. A similar con-clusion was reached independently by Angly et al. [7],based on a strong phylogenetic signal as measured byPagel’s λ [12]. However, neither Blomberg’s K nor Pagel’sλ make any statement about time scales (nor phylogeneticscales) over which traits vary. While 16S GCN variationis relatively rare within species, variation increases withtaxonomic distance [13] and this may lead to inaccuratepredictions for the many clades which are distant fromsequenced genomes. To date, no independent evaluationof existing 16S GCN prediction tools has been published.To resolve these uncertainties, we assessed the phy-logenetic autocorrelation of 16S GCNs across bacteriaand archaea (prokaryota) in a phylogenetic tree com-prising ∼ 570,000 OTUs (99% similarity in 16S rRNA),based on ∼ 6800 quality-checked complete sequencedgenomes. The tree was constructed from sequences inSILVA and partly constrained using SILVA’s taxonomicannotations. We predicted 16S GCNs using several com-mon phylogenetic reconstruction methods and exam-ined the accuracy achieved by each method for OTUsin the SILVA-derived tree. We assessed the predictiveaccuracy as a function of an OTU’s nearest-sequenced-taxon-distance (NSTD), that is, the minimum phyloge-netic distance (mean nucleotide substitutions per site) ofthe OTU to the nearest sequenced genome. We note thatthe average NSTD for a particular microbial community,weighted by OTU frequencies, is known as its nearestsequenced taxon index (NSTI; [6]). Further, we system-atically assessed the predictive accuracy of three recenttools for correcting 16S GCNs in microbiome surveys,PICRUSt [6], CopyRighter [7], and PAPRICA [8], whichtogether have been cited over 1000 times. While PICRUStand PAPRICA were mainly designed to predict commu-nity gene content based on 16S amplicon sequences, theyautomatically include an intermediate step for predictingand correcting for 16S GCNs. We evaluate the accuracyof these tools using the known GCNs of the aforemen-tioned sequenced genomes, as a function of a genome’sNSTD. To further evaluate these tools under more real-istic scenarios, we also compare all tools to each otherfor OTUs in 635 prokaryotic communities sampled fromdiverse natural environments, including the ocean, lakes,hot springs, soil, bioreactors, and animal guts. We findthat 16S GCNs are moderately phylogenetically conservedand that prediction of 16S GCNs for the large numberof clades without sequenced genomes from close relativeswill generally be inaccurate. This conclusion is verified byour finding of low predictive accuracies by CopyRighter,PICRUSt, and PAPRICA, both for the sequenced genomesas well as when compared to each other on microbiomes.Using phylogenetically predicted 16SGCNs to correct 16Sread counts in microbiome surveys, as previously sug-gested [5–7], worked well only for a small number ofmicrobiomes.Results and discussionHow predictable are 16S GCNs from phylogeny?We found that the autocorrelation function of 16SGCNs, that is the correlation between the GCNs of tworandomly picked OTUs at a certain phylogenetic dis-tance, decays moderately with increasing phylogeneticdistance (Fig. 1a), dropping below 0.5 at a phylogeneticdistance of ∼ 15% and to zero at a phylogenetic dis-tance of ∼ 30% (nucleotide substitutions per site in the16S gene). Hence, predictions of 16S GCNs are expectedto be inaccurate for clades with an NSTD greater thanabout 15% and close to random for clades with an NSTDgreater than about 30%. To explicitly test this conclu-sion, we predicted 16S GCNs for randomly chosen tips ofour SILVA-derived tree and compared these predictionsto the GCNs known from complete sequenced genomes,where possible. We considered the following commonancestral state reconstruction algorithms for predictingGCNs: Sankoff ’s maximum-parsimony with various tran-sition costs [14], maximum-likelihood of Mk models withrerooting (equal rates model), weighted-squared-changeparsimony [15], phylogenetic independent contrasts (PIC)[16], and subtree averaging (arithmetic average of GCNsacross descending tips). CopyRighter and PICRUSt usePIC, while PAPRICA uses subtree averaging. We mea-sured the accuracy of each method using the cross-validated coefficient of determination(R2cv)[17]. The R2cvcorresponds to the fraction of variance explained by areconstruction algorithm, when tested against a sepa-rate set of randomly chosen sequenced genomes (“testset”) than those used for state reconstruction (“trainingset”). We assessed the R2cv depending on the NSTD ofthe test set, that is, the phylogenetic distance betweenthe test set and the training set. We observed that allprediction methods only achieved high accuracies (R2 0.6) for NSTDs below about 15–30% depending on themethod (Fig. 1c), consistent with our expectations basedLouca et al. Microbiome (2018) 6:41 Page 3 of 12bnumber of OTUsNSTD (% subst. per site)aphylogenetic distance (% subst. per site)cmin. NSTD (% subst. per site)autocorrelation of 16S GCN prediction accuracyfrequencies of NSTDsautocorrelation0 20 40 60 800100002000030000400005000060000Fig. 1 Phylogenetic signal of 16S gene copy numbers (SILVA-derived tree). a Pearson autocorrelation function of 16S GCNs depending onphylogenetic distance between tip pairs, estimated based on ∼ 6,800 sequenced genomes. b Distances of tips in the SILVA-derived tree to thenearest sequenced genome. Each bar spans an NSTD interval of 2%. c Cross-validated coefficients of determination (R2cv) for 16S GCNs predicted onthe SILVA-derived tree and depending on the minimum NSTD of the tips tested, for various ancestral state reconstruction algorithms (PIC:phylogenetic independent contrasts, WSCP: weighted squared-change parsimony, SA: subtree averaging, MPR: maximum parsimonyreconstruction, Mk: continuous-time Markov chain model with equal-rates transition matrix). MPR transition costs either increased exponentiallywith transition size (“exp”), proportionally to transition size (“pr”), or were equal for all transitions (“ae”). For analogous results using the original SILVAtree, see Additional file 1: Figure S1on the autocorrelation function. At NSTDs greater than∼ 40%, the R2cv drops below zero for all methods. Max-imum parsimony with exponentially weighted transitioncosts (“MPR.exp”) was generally the best performingmethod, while Mk model maximum-likelihood was theworst method (Fig. 1c).Within the SILVA-derived tree, about 49% of OTUs havean NSTD greater than 15% and about 30% of OTUs havean NSTD greater than 30% (Fig. 1b). We note that naturalmicrobial communities are generally not a purely randomsubsample of SILVA, since SILVA overrepresents organ-isms of clinical or industrial interest, and these organismsare generally expected to have low NSTDs. Further, itis likely that a much larger number of prokaryotes notyet included in SILVA, such as from recently discoveredor yet undiscovered phyla [18, 19], has NSTDs greaterthan 30%. Consequently, predictions of 16S GCNs basedon sequenced genomes alone are expected to be inaccuratefor themajority of extant prokaryotic clades in natural envi-ronments. We note that, in principle, errors in tree topol-ogy and branch lengths could contribute to a reducedpredictive accuracy of phylogenetic reconstruction tools(Fig. 1c). As we show below, however, our expectations onthe limited predictability of GCNs are verified in severaladditional and independent analyses, as well as using theoriginal SILVA tree (Additional file 1: Figure S1).Assessment of 3rd party prediction toolsThe preceding analysis suggests that phylogenetic predic-tion of 16S GCNs based on available sequenced genomesis bound to be inaccurate for a substantial numberof prokaryotic clades, especially those with only a fewsequenced representatives. This finding casts doubts overclaims that 16S GCNs can be accurately predicted fortypical environmental 16S sequences [5] and that 16SGCN corrections should be applied systematically toevery microbiome survey [7]. We thus tested the predic-tive accuracy of three recently published tools, PICRUStv1.1.1 [6], CopyRighter v0.46 [7], and PAPRICA v0.4.0b[8]. We performed two types of tests. In the first test,we compared the GCNs of the aforementioned sequencedgenomes to the GCNs predicted by each tool based on agenome’s 16S sequence. Because many of these genomeswere also used as input to CopyRighter, PICRUSt, andPAPRICA for model calibration in the original publica-tions (“calibration genomes”), or are closely related tothose calibration genomes, we evaluate the predictiveaccuracy of each tool depending on a genome’s distance(NSTD) from the tool’s calibration genomes. In the sec-ond test, we compared the predictions of each tool tothose of the other two tools, for all OTUs in the Green-genes 16S rRNA database [20] as well as for prokaryoticOTUs found in 635 microbiomes from a diverse range ofenvironments. For each tool, a slightly different approachwas taken depending on the tool’s particular design. ForPICRUSt and CopyRighter, we used their precomputedlookup tables listing predicted 16S GCNs for entries inGreengenes and mapped genomes (first test) as well asOTUs (second test) to Greengenes entries (at ≥ 99% sim-ilarity) to obtain the corresponding GCN predictions.For PAPRICA, we used the 16S rRNA sequences of thegenomes or OTUs as input to predict their GCNs throughthe PAPRICA pipeline.We find that the predictive accuracy of all three tools,evaluated on the sequenced genomes and measured interms of the fraction of explained variance of true GCNsLouca et al. Microbiome (2018) 6:41 Page 4 of 12(R2), generally decreases with a genome’s NSTD (Fig. 2).Specifically, while accuracy is moderate to high at lowNSTDs (R2 > 0.6 for NSTDs < 10%), the R2 drops below 0.5for genomes with NSTDs above 30%. In fact, for PICRUStand PAPRICA, the R2 even becomes negative for NSTDsabove 30%, which is worse than if the average GCN overall genomes had been used as prediction. We also findpoor agreement between the predictions of different tools,when compared to each other across entries in the Green-genes database. When evaluated over the entire Green-genes database, GCNs predicted by any tool explainedat most 25% of the variance in the predictions of othertools (R2 < 0.25; Fig. 3a–c). It is noteworthy that Copy-Righter and PICRUSt use the same set of input genomes(∼ 3000 genomes from the IntegratedMicrobial Genomesdatabase; [21]) and similar reference trees (Greengenesreleases October 2012 and May 2013, respectively), andyet, GCN predictions differ substantially between Copy-Righter and PICRUSt (R2 = 0.23; Fig. 3a). When we con-sidered the agreement between tools depending on anOTU’s NSTD (Fig. 3d–f), we found that the R2 decreasesrapidly with increasing NSTD and becomes negative atNSTDs below 20%. We also found that the frequencydistributions of 16S GCNs predicted across Greengenes(Additional file 1: Figure S1) as well as across the genomes(Additional file 1: Figure S4) differ substantially betweentools. For example, CopyRighter, PICRUSt, and PAPRICApredict that the most common GCN across genomes is 1,3, and 2, respectively. As seen in Fig. 3a, c, CopyRighterindeed appears to underestimate GCNs when comparedto PICRUSt and PAPRICA.When we compared CopyRighter, PICRUSt, andPAPRICA for OTUs detected in any of the 635 samples,we found that the tools only agreed moderately to poorlywith each other for the majority of the samples. Specifi-cally, for any given pair of tools (CopyRighter vs. PICRUSt,PICRUSt vs. PAPRICA, or CopyRighter vs. PAPRICA),the fraction of variance in predictions of the 1st tool thatwas explained by predictions of the 2nd tool (R2) wasbelow 0.5 for over 84% of the samples and below 0.1 forover 55% of the samples (Fig. 4). In many cases, the agree-ment between tools was even worse than if predictionswere uncorrelated between tools(R2 < 0). A negative R2may be indicative of “overfitting” during extrapolationof GCNs to OTUs with large NSTDs. The worst agree-ment was found between PICRUSt and PAPRICA (meanR2 = − 0.70), while the best (but still bad) agreement wasfound between CopyRighter and PICRUSt (mean R2 =− 0.41). Even when we only considered animal-associatedsamples (e.g., from human guts or skin), which are consid-ered better studied than other environments and generallyhave lower NSTIs (weighted mean NSTD of consideredOTUs), we found frequent bad agreements between tools.One explanation is that even in human-associated micro-biomes, many OTUs had large NSTDs and were drivingoverall predictive accuracy down. In fact, we find thatthe poor agreement between tools in most samples is notdriven solely by a few outlier OTUs but is a reflectionof moderate to poor agreement for a large number ofOTUs in each sample (Additional file 1: Figures S5, S8,and S9). The agreement between tools generally decreasedfor larger NSTIs, although this trend became much morepronounced when we considered animal-associated sam-ples separately from non-animal-associated samples. Thestrongest trend was observed in non-animal-associatedsamples, where the R2 and NSTI exhibited significantCopyRighterNSTD (% subst. per site)aPICRUStNSTD (% subst. per site)bPAPRICAcNSTD (% subst. per site)Fig. 2 Evaluation of GCN prediction tools on genomes with known GCNs. Accuracy of GCN predictions by CopyRighter (a; [7]), PICRUSt (b; [6]), andPAPRICA (c; [8]) for sequenced genomes, as a function of the genome’s NSTD. NSTDs were calculated separately for each tool, based on the set ofgenomes used to calibrate the tool by its authors. Accuracy was measured in terms of the coefficient of determination, i.e. the fraction of variance intrue GCNs explained by each tool (R2). Genomes were binned into equally sized NSTD intervals (i.e., 0–10%, 10–20% etc.), and the R2 was calculatedseparately for genomes in each bin (one plotted point per bin). Only bins with at least 10 genomes are shownLouca et al. Microbiome (2018) 6:41 Page 5 of 12Fig. 3 Comparisons of 16S GCN predictions between tools across Greengenes. a 16S GCNs predicted by CopyRighter (vertical axis; [7]) and PICRUSt(horizontal axis; [6]) across OTUs (99% similarity) in the Greengenes 16S rRNA reference database (release May 2013; [20]). One point per OTU.b Comparison of predicted 16S GCNs by PICRUSt and PAPRICA, similarly to (a). c Comparison of predicted 16S GCNs by CopyRighter and PAPRICA,similarly to (a). Diagonal lines are shown for reference. Fractions of explained variance (R2, X-axis explaining Y-axis) and the number of consideredOTUs (n) are written in each figure. d–f Fractions of explained variance (R2) as a function of an OTU’s NSTD, for each compared pair of tools in a–c.OTUs were binned into equally sized NSTD intervals (i.e., 0–5%, 5–10% etc.), and the R2 was calculated separately for OTUs in each bin (one plottedpoint per bin). Only bins with at least 10 OTUs are shown(P < 0.05) Pearson correlations between − 0.35 and− 0.52, depending on the tools compared (Fig. 4a–c).In animal-associated samples, the R2 and NSTI exhib-ited significant Pearson correlations between − 0.17 and− 0.32 (Fig. 4d–f). The above observations are consistentwith our expectation that GCN predictions will only beaccurate for a small fraction of naturally occurring micro-biomes, namely microbiomes with low NSTIs ( 15%for the samples examined here), although tools occasion-ally disagreed substantially even on samples with lowNSTIs. We note that here, we recovered OTUs by closed-reference clustering to SILVA, thereby omitting cladesnot represented in SILVA at all. It is likely that many ofthese omitted clades, especially those from poorly stud-ied phyla, had even greater NSTDs than typical closed-reference OTUs. This realization further strengthens ourconclusions that existing GCN prediction tools performpoorly for many of those samples.Previous studies have used mock communities totest the predictability of 16S GCNs, demonstrating thatcorrecting for GCNs improves estimates of microbialcommunity composition, provided that corrections areaccurate [6, 7, 22]. While mock communities using cul-tured and sequenced strains are convenient test cases(since GCNs are known for each member), they lead to abiased assessment of predictive accuracy because strainsused in mock communities are likely to be found among(or closely related to) sequenced genomes. In other words,use of mock communities — instead of natural com-munities as in the present study — can yield the falseimpression that GCNs can be well predicted for typicalnatural microbial communities.ImplicationsAccurate correction of 16S sequence variant counts forGCNs in microbiome surveys would undoubtedly reduceLouca et al. Microbiome (2018) 6:41 Page 6 of 12bromeliad tanks (Brazil)Hot Lake sediments (Italy)Indian Ocean watersoil (Papua New Guinea)hot spring (India)hypersaline lake water (Romania)soil (Russia)Tibetan Plateau lake waterAtlantic transect watergrassland soil (USA)Yanga Lake water (Australia)hypersaline mat (Bahamas)marine sediments (Hangzhou Bay)lake water (Yosemite National Park)toluene degrading bioreactor (Belgium)restinga ponds and lagoons (Brazil)restinga soil (Brazil)biogas bioreactor (Denmark)a b cNSTI (% subst. per site)RNSTI (% subst. per site) NSTI (% subst. per site)CopyRighter vs PAPRICACopyRighter vs PICRUSt PICRUSt vs PAPRICAr = -0.52 r = -0.35 r = -0.50d e fRr = -0.32 r = -0.17 r = -0.23katydid guts (Anabrus sp.)mosquito guts (Anopheles sp.)human child guts (China)women skin (China)honeybee gutshuman cadaver gutshuman guts (HIV experiment)lizard gutsporcine gutsbird guts (Venezuela)animal-associatednon-animal associatedFig. 4 Agreement of GCN prediction tools in microbial communities, depending on the NSTI. a Agreement between 16S GCNs predicted byCopyRighter and PICRUSt (in terms of the fraction of variance in the former explained by the latter, R2) for non-animal-associated microbialcommunities, compared to the nearest sequenced taxon index (NSTI) of each community. Each point represents the R2 and the NSTI of onemicrobial community sample. b, c Similar to a, but comparing PICRUSt to PAPRICA (b) and CopyRighter to PAPRICA (c). d–f Similar to a–c, butshowing animal-associated samples. In all figures, linear regression lines are shown for reference. Pearson correlations between R2 and NSTI (r2,written in each figure) were statistically significant (P < 0.05) in all cases. Points are shaped and colored according to the original study, as listed inthe legend. Note the negative relationship between a community’s NSTI and the pairwise agreement of GCN prediction tools for that community.For a similar figure showing the spread of NSTDs in each sample, see Additional file 1: Figures S5. For detailed comparisons between tools onindividual samples see Additional file 1: Figures S6 and S7biases in cell-count estimates. As we have shown, how-ever, predicting 16S GCNs can come at a cost of sub-stantial additional errors (“noise”) when affected clades donot have close relatives with sequenced genomes. Theseerrors can even vary strongly between tools (Figs. 3 and 4).In principle, GCN corrections may be applied selec-tively to only those taxa with a sufficiently low NSTD,although this complicates the interpretation of micro-bial community profiles that include taxa with no GCNcorrection. Adoption of GCN corrections, dependent orindependent of NSTDs, as suggested by other authors[5, 7], could thus compromise the comparability betweenstudies. We recommend a careful consideration of thesecaveats before correcting for GCNs in typical microbiomestudies, until coverage by sequenced genomes is substan-tially improved. A similar conclusion was recently reachedby Edgar [22], based on tests of a specific GCN correctionalgorithm on mock microbial communities. For example,if the detection of spatiotemporal variation in commu-nity composition is the sole objective in a study, then thisvariation could be described merely in terms of 16S genevariants without the need for normalizing by GCNs. Anotable exception are cases that necessitate knowledge oftrueOTUproportions in a community, such as for biogeo-chemical modeling or for estimating gene proportions in acommunity using tools such as PICRUSt [6] or PAPRICA[8]. In these cases, it may be reasonable to correct for16S GCNs despite the high additional noise, although theeffects of this noise on estimated gene abundances remainto be investigated. Recent tools such as “16Stimator” couldhelp extend coverage to de novo assembled draft genomes[23], for which GCN counts have been hard to estimatein the past due to misassembly of ribosomal genes. Ourfindings also point to the need to explore alternativegenes with more conserved copy numbers for phyloge-netic community profiling, such as recA or rpoB [24, 25].Larger reference databases for phylogenetic identificationof sequence variants of these genes are needed, however,to make these genes a more widely adopted alternative to16S rRNA.More generally, our work demonstrates the importanceof a cautious interpretation of evolutionary statistics (inthis case, Blomberg’s K and Pagel’s λ [5, 7]) to avoidhasty conclusions about the predictability of a trait usingLouca et al. Microbiome (2018) 6:41 Page 7 of 12phylogenetic methods. Indeed, two of the most importantfactors influencing the predictability of a trait is the phy-logenetic distance to the nearest clade with known traitvalue (in our case, a sequenced genome) and the typicaldepth at which the trait is conserved [10, 26]. Unitlessstatistics such as K and λ indicate that a phylogeneticsignal exists but say little about the two factors that arecrucial for accurately predicting traits based on phylogenyin practice. Similar considerations are warranted whenphylogenetic extrapolation tools such as PICRUSt [6] orPAPRICA [8] and interpolation tools such as FAPROTAX[27] are used to estimate metabolic traits from phylogeny.Indeed, Langille et al. [6] emphasize that the accuracy ofPICRUSt depends on the availability of closely sequencedrelatives. We point out that phylogenetic conservatismvaries strongly across traits [28, 29], and thus, our abil-ity to accurately predict these traits using phylogeneticmethods also varies considerably.ConclusionsHere, we have assessed the phylogenetic conservatism of16S GCNs and examined the predictability of GCNs usingseveral common phylogenetic reconstruction algorithms,as a function of a clade’s nearest sequenced taxon distance.Our findings suggest that GCNsmay currently not be pre-dictable for a substantial fraction of extant prokaryoticclades. Further, we independently evaluated the accuracyof available 16S GCN prediction tools [6–8] on a set ofcompletely sequenced genomes, as well as for OTUs inthe Greengenes 16S database and in microbial commu-nities from a wide range of environments. Our analysisrevealed that existing tools perform poorly for a largefraction of the genomes and OTUs tested. For over 85%of the examined microbial communities, GCN predic-tions differed strongly between any two tools compared(R2 < 0.5). Thus, contrary to common assumption, 16SGCN predictions are currently bound to be inaccurate fora substantial fraction of extant prokaryotic diversity dueto insufficient coverage by sequenced genomes. We there-fore recommend that 16S GCNs should only be correctedfor in surveys of microbial communities with a low NSTI( 15%), unless the high additional noise introduced isjustified by a need to estimate true cell proportions.Materials andmethodsConstruction of SILVA-derived treeWhile the original SILVA tree is well curated taxonom-ically, it is mostly meant to be used as a guide tree, andre-calculation of branch lengths is generally advisedfor downstream phylogenetic analyses [30]. Here, toconstruct a phylogenetic tree with more meaningfulbranch lengths using OTUs in the SILVA non-redundant(NR99) 16S database (release 128; [4]), we proceeded asfollows. Aligned representative SSU sequences in SILVAwere reduced by first removing nucleotide positionswith > 95% gaps and then removing the top 5% mostentropic nucleotide positions. Taxonomic identitiesprovided by SILVA for OTUs at the domain, phylum,and class level were used to create split constraints forFastTree [31], by constraining each taxon to be on asingle side of a split and monophyletic. Taxa with fewerthan 10 OTUs were omitted from the constraints. Atotal of 354 constraints were thus defined. Using thetaxonomically generated constraints and taking theoriginal SILVA tree as a starting tree, we constructeda phylogenetic tree from the reduced alignmentswith FastTree v2.1.10 (options “-spr 4 -gamma-fastest -no2nd -constraintWeight 100”).The phylogenetic tree was rerooted so that bacteria andArchaea are split at the root. Our SILVA-derived tree isprovided as Additional file 2. For all downstream anal-yses, chloroplasts, mitochondria, and Eukaryota wereomitted from the tree. In the main article, we describeour analyses using this SILVA-derived tree (Fig. 1); anal-ogous results for the original SILVA tree are shown inAdditional file 1: Figure S1.Phylogenetic distribution of 16S GCNsTo examine how 16S GCNs are distributed phyloge-netically and to assess their general predictability usingvarious phylogenetic methods, we proceeded as fol-lows. A total of 8,767 annotated bacterial and archaealgenomes with completion status “Complete Genome”were downloaded from the NCBI RefSeq database onJanuary 4, 2018. Downloaded genomes were checked forpotential contamination using checkM 1.0.6 [32] (option“reduced_tree”), which is based on the detection ofconserved marker genes (assembly and checkM sum-maries in Additional file 3). Genomes found to exhibit acontamination level above 1% or a strain heterogeneityabove 1% were discarded, leaving us with 6,868 completegenomes for downstream analysis (Additional file 4).For each genome, 16S GCNs were determined usingtwo approaches: First, we counted the number ofannotated 16S rRNA sequences in the NCBI anno-tations (files rna_from_genomic.fna). Second, weused covariance models with the program cmsearch(as part of INFERNAL version 1.1.2, options “--noali--cut_nc”) to search for 16S rRNA sequences withinthe assembled genomes (files genomic.fna). Separatecovariance models for archaeal and bacterial 16S rRNAgenes were obtained from the Rfam database [33] (acces-sions RF00177 and RF01959). A table listing GCNs calcu-lated using both methods is provided as Additional file 5.Only genomes for which the two methods yielded thesame 16S GCNs were considered for subsequent analy-sis, yielding 16S GCNs for 6,780 genomes (“high-qualitygenomes,” Additional file 6). The accuracy of these GCNsLouca et al. Microbiome (2018) 6:41 Page 8 of 12was further verified through comparison to the RibosomalRNA Operon Copy Number Database (rrnDB, accessedon June 7, 2017; [34]) whenever a genome assembly acces-sion was present in the rrnDB (rrnDB attribute “Datasource record id”). Across 5,616 high-quality genomestested, we found an almost-perfect agreement with therrnDB (R2 > 0.999; Additional file 1: Figure S2). checkMquality summaries for the high-quality genome set areprovided as Additional file 7.Tips on the SILVA-derived tree were mapped tohigh-quality genomes, whenever possible, as follows:First, representative 16S sequences of SILVA OTUswere aligned to the longest 16S rRNA sequencefrom each genome using vsearch 2.3.4 [35] at maxi-mum (100%) similarity (vsearch options “--strandboth --usearch_global --maxaccepts 0 --top_hits_only --iddef 0 --id 1.0”). If an OTU aligned tomultiple genomes, all genomes were initially kept. Next,for each aligned OTU-genome pair, we compared theNCBI taxon ID (“taxid”) of the OTU to that of the genome.OTU taxids were obtained from a lookup table providedby SILVA (https://www.arb-silva.de/fileadmin/silva_databases/release_128/Exports/taxonomy/taxmap_embl_ssu_ref_128.txt.gz). Genome taxids were obtained fromlookup tables provided by NCBI (ftp://ftp.ncbi.nlm.nih.gov/genomes/genbank/*/assembly_summary.txt, where“*” is either “bacteria” or “archaea”). Any alignedOTU-genome pair with non-identical taxids was omit-ted. Of the remaining OTU-genome pairs with identicaltaxids, we only kept the first aligned genome for eachOTU. A total of 9,395 OTUs could thus be mappedto one of the genomes. For each mapped OTU, weassumed a GCN equal to the GCN counted for the cor-responding genome. For all other OTUs, we assumed anunknown GCN.All phylogenetic analyses were performed using theR package castor [36], available at The Comprehen-sive R Archive Network (CRAN). NSTDs for all tipswith respect to tips mapped to a sequenced genome(Fig. 1b) were calculated using the castor functionfind_nearest_tips. The phylogenetic autocorrela-tion function (ACF) of known 16S GCNs across theSILVA-derived tree (Fig. 1a) was calculated using thecastor function get_trait_acf based on 108 tippairs (options “Npairs=1e8, Nbins=100”), chosenrandomly among tips with known GCN. The functionget_trait_acf randomly picks OTU pairs on thetree, bins them into one of many intervals of phy-logenetic distance, and calculates the Pearson auto-correlation between GCNs of the OTU pairs withineach bin. Note that this analysis does not assume thatGCNs scale linearly with phylogenetic distance. Instead,the ACF merely measures the statistical correlationbetween GCNs on distinct tips, conditional upon the tipsbeing within a certain phylogenetic distance from eachother.GCNs were reconstructed on the SILVA-derivedtree using Sankoff ’s maximum-parsimony (function hsp_max_parsimony, with option transition_costseither set to “exponential,” “proportional,” or“all_equal”), phylogenetic independent contrasts(functionhsp_independent_contrasts), weighted-squared-change parsimony (function hsp_squared_change_parsimony), subtree averaging (functionhsp_subtree_averaging), and maximum-likelihoodof Mk models with rerooting (function hsp_mk_model_rerooting with options root_prior=‘empirical’, optimization_algorithm=‘nlminb’,Ntrials=5, rate_model=‘ER’).To calculate the cross-validated fraction of variance pre-dicted by (aka. cross-validated coefficient of determina-tion of) each method (R2cv; [17]) as a function of the NSTD(Fig. 1c), we proceeded as follows. We randomly chose 2%of the tips with known 16S GCN to be excluded from theinput to the reconstructions and to be used as an indepen-dent “test set” afterwards. Depending on the NSTD cutoffconsidered (for example 10% substitutions per site), wealso excluded all tips whose phylogenetic distance to thetest set was below the NSTD cutoff. The remaining tipswith known GCNs (“training set”) were used as input tothe reconstructions, and the GCNs predicted for the testset were then compared to the known GCNs of the testset. This process was repeated three times and the result-ing R2 was averaged over all repeats, yielding an R2cv foreach considered NSTD cutoff. The R script for analyzingand reconstructing 16S GCNs across the SILVA-derivedtree is available as Additional file 8. For comparison, all ofthe above analyses were also performed using the originalSILVA guide tree (Additional file 1: Figure S1).Evaluation of 3rd party GCN prediction tools on sequencedgenomesTo test the predictive accuracy of CopyRighter [7],PICRUSt [6], and PAPRICA [8] for genomes with knownGCNs, we compared their predictions with the GCNscounted in the (high-quality) sequenced genomes. Toevaluate the predictive accuracy of CopyRighter [7]on the genomes, we proceeded as follows: We firstdownloaded the precomputed lookup table listing Copy-Righter’s predictions for the Greengenes 16S rRNAdatabase (release October 2012, “GG2012”; [20]), fromthe project’s Github on June 6, 2017 (v0.46): https://github.com/fangly/AmpliCopyRighter (CopyRighter-0.46/data/201210/ssu_img40_gg201210.txt). We thenaligned the longest 16S rRNA sequence of eachgenome to OTUs (clustered at 99% similarity) inthe Greengenes database using vsearch (vsearchoptions “--strand both --usearch_globalLouca et al. Microbiome (2018) 6:41 Page 9 of 12--maxhits 1 --maxaccepts 10 --top_hits_only”), always choosing the best match in Greengenesand keeping only genomes that mapped to a Greengenesentry by at least 99% similarity (5688 genomes mapped).For each mapped genome, we took the GCN predicted byCopyRighter for the corresponding Greengenes entry asCopyRighter’s prediction for the genome. This predictionwas then compared to the GCN counted from the genomesequence. A histogram of CopyRighter’s predictionsacross mapped genomes is shown in Additional file 1:Figure S4B. The predictive accuracy of CopyRighter wasmeasured in terms of the fraction of explained variance(R2), as a function of a genome’s NSTD (Fig. 1a). NSTDsof genomes were calculated as described in a separatesection below.A similar approach was used for PICRUSt [6]: Theprecomputed lookup table listing PICRUSt’s predic-tions for the Greengenes database (release May 2013;“GG2013”) was downloaded from the project’s web-site on June 6, 2017 (v1.1.1): https://picrust.github.io/picrust/picrust_precalculated_files.html (16S_13_5_precalculated.tab.gz). A total of 5,708 high-quality genomes could be mapped to an OTU (99%similarity) in GG2013. A histogram of PICRUSt’s predic-tions across all mapped genomes is shown in Additionalfile 1: Figure S4C. The predictive accuracy of PICRUSt wasmeasured in terms of the R2 as a function of a genome’sNSTD (Fig. 1b), similarly to CopyRighter.To evaluate the predictive accuracy of PAPRICA [8]on the genomes, we proceeded as follows: We firstdownloaded and installed PAPRICA from the project’sGithub on June 6, 2017 (v0.4.0b): https://github.com/bowmanjeffs/paprica. This release includes precomputedreference trees (one for archaea and one for bacteria) andtables listing 16S GCNs for the tool’s calibration genomesrepresented in the reference trees. We used the longest16S rRNA sequence from each genome as an input tothe PAPRICA pipeline (command “paprica-run.sh”),separately for archaea and bacteria. The pipeline pro-duces, among others, a table listing the uncorrected abun-dance of each unique input sequence (this can be greaterthan 1 if multiple genomes share the same 16S rRNAsequence) and the corresponding corrected abundance(after dividing by the predicted 16S GCN). We used thistable to obtain the 16S GCNs predicted by PAPRICAfor the unique 16S sequences (representing 3473 16Ssequences), by dividing the uncorrected by the correctedabundance. We then compared these predicted GCNs tothe GCNs counted in the genome sequences, similarly toabove. A histogram of PAPRICA’s predictions across allrepresented genomes is shown in Additional file 1: FigureS4D. The predictive accuracy of PAPRICA was measuredin terms of the R2 as a function of a genome’s NSTD(Fig. 1a), similarly to CopyRighter.Comparison of 3rd party GCN prediction tools acrossGreengenesTo compare the predictions by CopyRighter to thoseby PICRUSt across all OTUs in Greengenes (Fig. 3a),we first mapped all OTUs in GG2013 to OTUs inGG2012 using vsearch (with options “--strand both--usearch_global”). We only kept matches at 100%similarity (153,375 out of 203,452 OTUs in GG2013). Toeach mapped OTU in GG2013, we compared the corre-sponding GCN predicted by PICRUSt to the GCN pre-dicted by CopyRighter for the matched OTU in GG2012.To calculate the frequency distributions of GCNs pre-dicted by CopyRighter and PICRUSt across all OTUsin Greengenes (histograms in Additional file 1: FigureS3A,B), we used the GCNs listed in their precomputedlookup tables.To compare PAPRICA to PICRUSt across Green-genes (Fig. 3b), we proceeded as follows: Representativesequences of OTUs in GG2013 were split into archaealand bacterial sequences. Each resulting fasta file was usedas input to the PAPRICA pipeline to predict the corre-sponding 16S GCN, as described above for genomes. Thisyielded a predicted GCN for all Greengenes entries. Thesepredictions were compared to the precomputed GCN val-ues provided by PICRUSt. These predictions were alsoused to calculate the frequency distribution of GCNspredicted by PAPRICA across Greengenes (Additionalfile 1: Figure S3C). To compare CopyRighter to PAPRICA(Fig. 3c), we proceeded as described above for the com-parison of CopyRighter to PICRUSt.Comparison of 3rd party GCN prediction tools acrossmicrobial communitiesTo compare CopyRighter, PICRUSt, and PAPRICAacross OTUs in various microbial communities, we pro-ceeded as follows. Publicly available 16S rRNA ampliconsequence data from various environmental samples weredownloaded from the European Nucleotide Archive(http://www.ebi.ac.uk/ena). Only Illumina sequence datafrom amplicons obtained using bacteria- and/or archaea-sensitive primers were considered. Samples were chosento cover a wide range of environments, including theocean, marine and lake sediments, soil, saline and hyper-saline lakes, hydrothermal vents, hot springs, bioreactors,and animal-associated microbiomes. All sequencingdata were processed in a similar way, where possible,as follows. Overlapping paired-end reads were mergedusing flash v1.2.11 [37] (options -min-overlap=20-max-overlap=300 -max-mismatch-density0.25 -phred-offset=33 -allow-outies), andnon-overlapping paired-end reads were omitted.Single-end reads were kept unchanged. All single-endreads and merged paired-end reads were then qualityfiltered using vsearch v2.4.3 [35] (options -fastq_Louca et al. Microbiome (2018) 6:41 Page 10 of 12ascii 33 -fastq_minlen 120 -fastq_qmin 0-fastq_maxee 1 -fastq_truncee 1 -fastq_maxee_rate 0.005 -fastq_stripleft 7). Sam-ples with more than 20,000 quality-filtered reads wererarefied down to 20,000 reads to reduce computationtime, by randomly picking reads without replacement.Quality-filtered sequences were clustered into operationaltaxonomic units (OTUs; at 97% similarity) by closed-reference global aligning against the non-redundant(NR99) SILVA SSU reference database (release 128; [4]),using vsearch. Both strands were considered for align-ment (vsearch option --strand both). Sequences notmatching any database entry at 97% similarity or higherwere discarded. Note that OTUs were thus representedby SILVA entries, namely the ones used to seed theclusters. Chloroplasts, mitochondria, and any Eukaryotawere omitted. OTUs represented by fewer than five readsacross all samples were omitted. Finally, any sampleswith fewer than 2,000 reads accounted for by OTUs wereomitted. This yielded an OTU table with 635 samplesand 65,673 OTUs represented by 4,827,748 reads (onaverage 734 OTUs per sample). Sample accession num-bers, coordinates, sampling dates, original publications,sequencing platforms, quality-filtered read lengths, andread counts and covered primer regions (where available)are provided in Additional file 9.To predict GCNs for OTUs in each sample usingCopyRighter, we used the same approach as forgenomes: Representative 16S sequences of OTUswere aligned to GG2012 using vsearch (options“--strand both --usearch_global --iddef0 --id 0.99 --maxhits 1 --maxaccepts 10--top_hits_only”), omitting any OTUs not matchedto a Greengenes entry by at least 99% similarity. Foreach OTU kept, the GCN listed by CopyRighter for thematched Greengenes entry was taken as CopyRighter’sprediction. For PICRUSt, we proceeded in an analogousway, using GG2013 instead of GG2012. For PAPRICA,we proceeded in an analogous way, using PAPRICA’sGCN predictions computed previously for GG2013 (seeprevious section).To compare any two given tools (CopyRighter vs.PICRUSt, PICRUSt vs. PAPRICA, or CopyRighter vs.PAPRICA) for a specific sample, only OTUs with at leastone read in the sample and having a GCN prediction fromboth tools were considered. We measured the agreementbetween two tools in terms of the fraction of variance inpredictions of the 1st tool that was explained by predic-tions of the 2nd tool (R2).We calculated the sample’s NSTI(nearest sequenced taxon index) according to [6], i.e., asthe arithmetic average NSTD over all OTUs considered inthe comparison and weighted by relative OTU frequen-cies. Details on how NSTDs were calculated are providedin the section below. For each pair of tools compared, wethus obtained 635 NSTIs and 635 R2s across 635 sam-ples, shown in Fig. 4. Pearson correlation coefficients(r2)between NSTIs and R2 were calculated for each pair oftools, separately for animal-associated and non-animal-associated samples. Statistical significances (P values) ofcorrelation coefficients were estimated using a permuta-tion test with 1000 permutations. Additional file 1: FiguresS6 and S7 show GCNs predicted by each tool for variousmicrobial communities. We also show relative deviationsbetween tools (|A − B| /((A + B)/2), where A and B areGCNs predicted by two tools for the same OTU) andNSTDs for OTUs in various samples (Additional file 1:Figure S8).Evaluation and comparison of GCN prediction toolsdepending on NSTDTo examine the predictive accuracy of CopyRighter,PICRUSt, and PAPRICA as a function of an OTU’sor genome’s NSTD, we proceeded as follows. Foreach OTU in SILVA, and separately for each tool,we calculated the NSTD as the phylogenetic dis-tance to the nearest sequenced genome used by thetool to make predictions (“calibration genomes”). ForPAPRICA, a list of 5,628 calibration genomes wasobtained from PAPRICA’s precomputed files (PAPRICA/ref_genome_database/*/genome_data.final.csv, where “*” is either bacteria or archaea). Cali-bration genomes were matched to SILVA OTUs via globalalignment of the 16S gene at a similarity threshold of 99%,using vsearch. Matched OTUs were assumed to have anNSTD equal to zero, and for all other SILVA OTUs, theNSTD was calculated based on the SILVA-derived treeand using the R package castor [36]. An approximatematching of genomes to OTUs (i.e., at 99% similarity)was chosen to ensure that as many of the calibrationgenomes are included as possible; note that SILVA OTUsare themselves clustered at that similarity and that theerror potentially introduced to the NSTDs and NSTIs isnegligible (< 1% nucleotide substitutions per site). ForPICRUSt, a table was downloaded from the project’s web-site listing IMG (Integrated Microbial Genomes) IDs for2,887 calibration genomes (https://github.com/picrust/picrust/tree/master/tutorials/picrust_starting_files.zip,file GG_to_IMGv350.txt). IMG IDs were translatedto GG2013 sequence IDs using the gg_13_5_img.txtlookup table downloaded from the Greengeneswebsite (http://greengenes.secondgenome.com/downloads).Matched GG2013 IDs were then mapped to SILVA OTUsvia global 16S sequence alignment with vsearch, at a simi-larity threshold of 99%. NSTDs of SILVA OTUs were thencalculated in the same way as for PAPRICA. For Copy-Righter, a lookup table was downloaded from the project’sGithub page that maps calibration genomes to GG2012sequences (https://github.com/fangly/AmpliCopyrighter,Louca et al. Microbiome (2018) 6:41 Page 11 of 12file AmpliCopyrighter-0.46/preprocessing/data/img_to_gg.txt). GG2012 sequences listed inthat table were mapped to SILVA OTUs, and NSTDs werecalculated for all SILVA OTUs, in a similar way as forPICRUSt. To determine the NSTDs for genomes exam-ined in this study (separately for CopyRighter, PICRUSt,and PAPRICA), genomes were mapped to SILVA OTUsvia global alignment of their longest available 16Ssequence at 99% similarity. For each genome, the NSTDof the most closely matched SILVA OTU was taken as thegenome’s NSTD. To determine NSTDs for all GreengenesOTUs, we mapped Greengenes OTUs to SILVA OTUs viaglobal alignment at 99% similarity. To determine NSTDsfor OTUs recovered from the sampled microbial commu-nities, we directly used the NSTDs of SILVA OTUs usedas seeds during closed-reference OTU picking. Whencomparing two GCN prediction tools on an OTU (e.g.,Figs. 3 and 4 and Additional file 1: Figure S8), in caseswhere the two NSTDs differed, we used their arithmeticaverage. To calculate the R2 between any two GCN pre-diction tools, or between a GCN prediction tool and the“true GCNs,” as a function of the NSTD (Figs. 2 and 3d–f),we binned the OTUs or genomes used in the comparisoninto equally sized NSTD intervals and calculated the R2separately for each interval. Only NSTD intervals with atleast 10 OTUs or genomes were considered.Additional filesAdditional file 1: Supplementary Information. (PDF 4320 kb)Additional file 2: SILVA-derived 16S tree (NR99). (TRE 23800 kb)Additional file 3: Assembly stats for all genomes. (TSV 13200 kb)Additional file 4: Assembly stats for quality-filtered genomes. (TSV 910 kb)Additional file 5: 16S GCNs for quality-filtered genomes using bothmethods. (TSV 135 kb)Additional file 6: 16S GCNs counted for high-quality genomes. (TSV 120 kb)Additional file 7: Assembly stats for high-quality genomes. (TSV 898 kb)Additional file 8: R script for phylogenetic analysis of 16S GCNs onSILVA. (R 34.3 kb)Additional file 9: Metadata for microbiome samples. (TSV 127 kb)Additional file 10: 16S GCNs for SILVA OTUs matched to high-qualitygenomes. (TSV 228 kb)Additional file 11: 16S GCNs predicted for SILVA 16S guide tree (NR99)via MPR.exp. (TSV 15200 kb)Additional file 12: 16S GCNs predicted for SILVA-derived 16S tree (NR99)via MPR.exp. (TSV 14700 kb)AbbreviationsGCN: Gene copy number; NSTI: Nearest sequenced taxon index; NSTD: Nearestsequenced taxon distanceAcknowledgementsWe thank Matthew Pennell for the comments on the manuscript. We thankAria S. Hahn for advice on analyzing the genomes. We thank two anonymousreviewers for helping us improve the manuscript.FundingS.L. was supported by an NSERC grant and a postdoctoral fellowship from theBiodiversity Research Centre, UBC. M.D. and L.W.P. acknowledge the supportof NSERC.Availability of data andmaterialsCorrespondence and requests for materials should be addressed to S.L.Supporting figures and tables, cited in the text, are provided as SupplementaryMaterial. Also provided is the list of 16S GCNs counted for the high-qualitygenome set (Additional file 6), as well as the list of GCNs assigned to the subsetof matched SILVA tips (Additional file 10). Calculated NSTDs and predicted 16SGCNs for all non-chloroplast, non-mitochondrial bacterial and archaeal OTUsin the SILVA guide tree and the SILVA-derived tree (method “MPR.exp”) areprovided as Additional files 11 and 12. All genomes are publicly available atthe NCBI RefSeq genome repository (ftp://ftp.ncbi.nlm.nih.gov/genomes). All16S rRNA amplicon reads of the 635 microbial communities considered arepublicly available on the NCBI Sequence Read Archive (https://www.ncbi.nlm.nih.gov/sra) under the accession numbers listed in Additional file 9. The Rscript used for analyzing the phylogenetic distribution of 16S GCNs on theSILVA tree is available as Additional file 8.Authors’ contributionsSL conceived the project, wrote the computer code, and performed all theanalyses. MD and LWP supervised the project and helped interpret the results.All authors contributed to the writing of the manuscript. All authors read andapproved the final manuscript.Ethics approval and consent to participateNot applicable.Consent for publicationNot applicable.Competing interestsThe authors declare that they have no competing interests.Publisher’s NoteSpringer Nature remains neutral with regard to jurisdictional claims inpublished maps and institutional affiliations.Author details1Biodiversity Research Centre, University of British Columbia, Vancouver,Canada. 2Department of Zoology, University of British Columbia, Vancouver,Canada. 3Department of Mathematics, University of British Columbia,Vancouver, Canada. 4Department of Botany, University of British Columbia,Vancouver, Canada.Received: 25 October 2017 Accepted: 30 January 2018References1. Turnbaugh PJ, Ley RE, Hamady M, Fraser-Liggett CM, Knight R,Gordon JI. The human microbiome project. Nature. 2007;449(7164):804–10.2. Gilbert JA, Jansson JK, Knight R. The Earth Microbiome project: successesand aspirations. BMC Biol. 2014;12(1):69. https://doi.org/10.1186/s12915-014-0069-1.3. Lima-Mendez G, Faust K, Henry N, Decelle J, Colin S, Carcillo F, Chaffron S,Ignacio-Espinosa JC, Roux S, Vincent F, Bittner L, Darzi Y, Wang J, Audic S,Berline L, Bontempi G, Cabello AM, Coppola L, Cornejo-Castillo FM,d’Ovidio F, De Meester L, Ferrera I, Garet-Delmas MJ,Guidi L, Lara E, Pesant S, Royo-Llonch M, Salazar G, Sánchez P,Sebastian M, Souffreau C, Dimier C, Picheral M, Searson S,Kandels-Lewis S, coordinators TO, Gorsky G, Not F, Ogata H, Speich S,Stemmann L, Weissenbach J, Wincker P, Acinas SG, Sunagawa S, Bork P,Sullivan MB, Karsenti E, Bowler C, de Vargas C, Raes J. Determinants ofcommunity structure in the global plankton interactome. Science.2015;348(6237):1262073. https://doi.org/10.1126/science.1262073.4. Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, Peplies J,Glöckner FO. The SILVA ribosomal RNA gene database project: improvedLouca et al. Microbiome (2018) 6:41 Page 12 of 12data processing and web-based tools. Nucleic Acids Res. 2013;41(D1):590–6. https://doi.org/10.1093/nar/gks1219.5. Kembel SW, Wu M, Eisen JA, Green JL. Incorporating 16S gene copynumber information improves estimates of microbial diversity andabundance. PLOS Comput Biol. 2012;8(10):1–11. https://doi.org/10.1371/journal.pcbi.1002743.6. Langille MG, Zaneveld J, Caporaso JG, McDonald D, Knights D, Reyes JA,Clemente JC, Burkepile DE, Thurber RLV, Knight R, et al. Predictivefunctional profiling of microbial communities using 16S rRNA markergene sequences. Nat Biotechnol. 2013;31(9):814–21.7. Angly FE, Dennis PG, Skarshewski A, Vanwonterghem I, Hugenholtz P,Tyson GW. CopyRighter: a rapid tool for improving the accuracy ofmicrobial community profiles through lineage-specific gene copynumber correction. Microbiome. 2014;2(1):11. https://doi.org/10.1186/2049-2618-2-11.8. Bowman JS, Ducklow HW. Microbial communities can be described bymetabolic structure: a general framework and application to a seasonallyvariable, depth-stratified microbial community from the coastal westantarctic peninsula. PLoS ONE. 2015;10(8):1–18. https://doi.org/10.1371/journal.pone.0135868.9. ISO 5725-1. Accuracy (trueness and precision) of measurement methodsand results - part 1: general principles and definitions. Technical report,International Organization for Standardization. 1994.10. Zaneveld JRR, Thurber RLV. Hidden state prediction: a modification ofclassic ancestral state reconstruction algorithms helps unravel complexsymbioses. Front Microbiol. 2014;5:431. https://doi.org/10.3389/fmicb.2014.00431.11. Blomberg SP, Garland Jr T, Ives AR, Crespi B. Testing for phylogeneticsignal in comparative data: behavioral traits are more labile. Evolution.2003;57(4):717–45.12. Pagel M. Inferring the historical patterns of biological evolution. Nature.1999;401(6756):877–84. https://doi.org/10.1038/44766.13. Vetrovsky T, Baldrian P. The variability of the 16S rRNA gene in bacterialgenomes and its consequences for bacterial community analyses. PLOSONE. 2013;8(2):1–10. https://doi.org/10.1371/journal.pone.0057923.14. Sankoff D. Minimal mutation trees of sequences. SIAM J Appl Math.1975;28(1):35–42. https://doi.org/10.1137/0128004.15. Maddison WP. Squared-change parsimony reconstructions of ancestralstates for continuous-valued characters on a phylogenetic tree. Syst Biol.1991;40(3):304–14.16. Felsenstein J. Phylogenies and the comparative method. Am Nat.1985;125(1):1–15.17. Shao J. Linear model selection by cross-validation. J Am Stat Assoc.1993;88(422):486–94.18. Hug LA, Baker BJ, Anantharaman K, Brown CT, Probst AJ, Castelle CJ,Butterfield CN, Hernsdorf AW, Amano Y, Ise K, Suzuki Y, Dudek N,Relman DA, Finstad KM, Amundson R, Thomas BC, Banfield JF. A newview of the tree of life. Nat Microbiol. 2016;1:16048. https://doi.org/10.1038/nmicrobiol.2016.48.19. Parks DH, Rinke C, Chuvochina M, Chaumeil PA, Woodcroft BJ,Evans PN, Hugenholtz P, Tyson GW. Recovery of nearly 8,000metagenome-assembled genomes substantially expands the tree of life.Nat Microbiol. 2017. https://doi.org/10.1038/s41564-017-0012-7.20. McDonald D, Price MN, Goodrich J, Nawrocki EP, DeSantis TZ, Probst A,Andersen GL, Knight R, Hugenholtz P. An improved greengenestaxonomy with explicit ranks for ecological and evolutionary analyses ofbacteria and archaea. ISME J. 2012;6(3):610–8.21. Markowitz VM, Chen I-MA, Chu K, Szeto E, Palaniappan K, Grechkin Y,Ratner A, Jacob B, Pati A, Huntemann M, Liolios K, Pagani I, Anderson I,Mavromatis K, Ivanova NN, Kyrpides NC. IMG: the integratedmetagenome data management and comparative analysis system.Nucleic Acids Res. 2012;40(D1):123–9.22. Edgar RC. UNBIAS: An attempt to correct abundance bias in 16S sequencing,with limited success. bioRxiv. 2017. https://doi.org/10.1101/124149.23. Perisin M, Vetter M, Gilbert JA, Bergelson J. 16Stimator: statisticalestimation of ribosomal gene copynumbers from draft genome assemblies.ISME J. 2016;10(4):1020–4. https://doi.org/10.1038/ismej.2015.161.24. Vos M, Quince C, Pijl AS, de Hollander M, Kowalchuk GA. A comparisonof rpoB and 16S rRNA as markers in pyrosequencing studies of bacterialdiversity. PLoS ONE. 2012;7(2):30600. https://doi.org/10.1371/journal.pone.0030600.25. McNair K, Edwards RA. Genomepeek—an online tool for prokaryoticgenome and metagenome analysis. PeerJ. 2015;3:1025. https://doi.org/10.7717/peerj.1025.26. Goberna M, Verdu M. Predicting microbial traits with phylogenies. ISME J.2016;10(4):959–67. https://doi.org/10.1038/ismej.2015.171.27. Louca S, Parfrey LW, Doebeli M. Decoupling function and taxonomy inthe global ocean microbiome. Science. 2016;353(6305):1272–7. https://doi.org/10.1126/science.aaf4507.28. Martiny AC, Treseder K, Pusch G. Phylogenetic conservatism offunctional traits in microorganisms. ISME J. 2013;7(4):830–8. https://doi.org/10.1038/ismej.2012.160.29. Martiny JBH, Jones SE, Lennon JT, Martiny AC. Microbiomes in light oftraits: A phylogenetic perspective. Science. 2015;350(6261):. https://doi.org/10.1126/science.aac9323.30. Glöckner FO, Yilmaz P, Quast C, Gerken J, Beccati A, Ciuprina A, Bruns G,Yarza P, Peplies J, Westram R, et al. 25 years of serving the communitywith ribosomal rna gene reference databases and tools. J Biotechnol.2017;261(169–176):. https://doi.org/10.1016/j.jbiotec.2017.06.1198.31. Price MN, Dehal PS, Arkin AP. Fasttree: Computing large minimumevolution trees with profiles instead of a distance matrix. Mol Biol Evol.2009;26(7):1641–50. https://doi.org/10.1093/molbev/msp077.32. Parks DH, Imelfort M, Skennerton CT, Hugenholtz P, Tyson GW.Assessing the quality of microbial genomes recovered from isolates,single cells, and metagenomes. Genome Res. 2014;25:1043–55. https://doi.org/10.1101/gr.186072.114.33. Griffiths-Jones S, Bateman A, Marshall M, Khanna A, Eddy SR. Rfam: anRNA family database. Nucleic Acids Res. 2003;31(1):439–41.34. Stoddard SF, Smith BJ, Hein R, Roller BR, Schmidt TM. rrnDB: improvedtools for interpreting rRNA gene abundance in bacteria and archaea anda new foundation for future development. Nucleic Acids Res. 2014;43:593–8.35. Rognes T, Flouri T, Nichols B, Quince C, Mahé F. VSEARCH: a versatileopen source tool for metagenomics. PeerJ. 2016;4:2584. https://doi.org/10.7717/peerj.2584.36. Louca S, Doebeli M. Efficient comparative phylogenetics on large trees.Bioinformatics. 2017. https://doi.org/10.1093/bioinformatics/btx701.37. Magoc T, Salzberg SL. Flash: fast length adjustment of short reads toimprove genome assemblies. Bioinformatics. 2011;27(21):2957–963.• We accept pre-submission inquiries • Our selector tool helps you to find the most relevant journal• We provide round the clock customer support • Convenient online submission• Thorough peer review• Inclusion in PubMed and all major indexing services • Maximum visibility for your researchSubmit your manuscript atwww.biomedcentral.com/submitSubmit your next manuscript to BioMed Central and we will help you at every step: