Combining whole genome shotgun sequencing and rDNA amplicon analyses to improve detection of microbe-microbe interaction networks in plant leaves

Microorganisms from all domains of life establish associations with plants. Although some harm the plant, others antagonize pathogens or prime the plant immune system, acquire nutrients, tune plant hormone levels, or perform additional services. Most culture-independent plant microbiome research has focused on amplicon sequencing of 16S rDNA and/or the internal transcribed spacer (ITS) of rDNA loci, but the decreasing cost of high-throughput sequencing has made shotgun metagenome sequencing increasingly accessible. Here, we describe shotgun sequencing of 275 wild Arabidopsis thaliana leaf microbiomes from southwest Germany, with additional bacterial 16S rDNA and eukaryotic ITS1 amplicon data from 176 of these samples. The shotgun data were dominated by bacterial sequences, with eukaryotes contributing only a minority of reads. For shotgun and amplicon data, microbial membership showed weak associations with both site of origin and plant genotype, both of which were highly confounded in this dataset. There was large variation among microbiomes, with one extreme comprising samples of low complexity and a high load of microorganisms typical of infected plants, and the other extreme being samples of high complexity and a low microbial load. We use the metagenome data, which captures the ratio of bacterial to plant DNA in leaves of wild plants, to scale the 16S rDNA amplicon data such that they reflect absolute bacterial abundance. We show that this cost-effective hybrid strategy overcomes compositionality problems in amplicon data and leads to fundamentally different conclusions about microbiome community assembly.


Introduction
Microorganisms affect many important plant traits. Opportunistic microbes can slow down growth or kill a plant, while beneficial ones can prime the plant immune system [1] , directly antagonize pathogens [2] , or indirectly inhibit pathogens by contributing to a suppressive environment [3] . Microbes may adjust plant hormone levels [4] and participate in nutrient acquisition [5,6] , among other processes [7] . Research in this area has revealed that most of the organisms present on and in healthy plant leaves are typically bacteria [8] . With the exception of pathogenic strains, fungi, archaea, and protists such as oomycetes generally have lower abundance on leaves than bacteria, and have also received less attention because they are less easily cultured and are genetically more complex. Due to the difficulty of adequately capturing microbial complexity and diversity, we still lack a good understanding of the composition and dynamics of leaf microbial communities, their abundances on the plant, and how they relate to other aspects of the host biology such as genotype, location, or environmental conditions. Amplicon sequencing, in which a specific locus common to a target group of organisms (usually 16S rDNA for bacteria and ITS for fungi) is amplified and sequenced, has been the tool of choice for revealing the taxonomic composition of a microbiome. Albeit usually very informative, amplicon sequencing relies on oligonucleotide primers that disfavor or exclude some organisms, with the consequence that some taxa are systematically ignored. Not only are there then no quantitative estimates of the excluded taxa, but there is also no information on the absolute abundances of included taxa. Furthermore, one gene cannot reliably predict the other genes and genetic and metabolic functions in a microbe beyond some conserved features [9,10] . Whole metagenome shotgun sequencing of DNA extracts has become an attractive tool for dissecting complex microbial communities, and databases and algorithms that support such efforts are being developed [11,12] . Shotgun sequencing supplies information on the total DNA content of microorganisms as opposed to a specific locus, which in principle enables us to ask many new types of questions [8,13,14] .
We used metagenome sequencing to characterize the leaf-associated (phyllosphere) microbiome of 275 wild Arabidopsis thaliana individuals from around Tübingen in southwest Germany, at four different timepoints between 2014 and 2016. Of these, we subjected 176 to 16S rDNA and ITS1 sequencing. We achieved low (<100 Mb) to high (>1 Gb) depths of microbe-associated metagenomic sequences per sample, which we mapped to reference databases. Unsurprisingly, the wild A. thaliana microbiota was highly variable between individuals, but there were some clear patterns in the dominant microbes; in particular Pseudomonas dominated in some sites and Sphingomonas in others. We estimated absolute microbial load as the ratio of microbial reads to plant chromosomal reads in each sample, and observed that the microbial load varied across samples from less than 1% to up to 77% of plant reads. After producing a load-corrected table, we often observed that intertaxa abundance correlations changed in sign compared to compositional amplicon data and metagenome data, in which microbial reads had been normalized by total sum scaling. Finally, we document the relative abundance of eukaryotic and archaeal microbes, revealing a small but noteworthy presence of fungi and oomycetes.

Sequencing and analysis approach
To obtain an unbiased picture of microbial diversity and microbial load in wild A. thaliana plants, we took a metagenomic shotgun sequencing approach to analyze entire wild phyllospheres. To determine host DNA content, quality-filtered sequencing reads were mapped to the A. thaliana Col-0 TAIR10 reference genome [16] with bwa-mem using standard parameters [17] . Sequences that did not map to the reference genome were then translated in silico in all six reading frames and aligned against NCBI's non redundant protein database using DIAMOND [18] with standard parameters, and alignments were processed with MEGAN [19,20] .
We reasoned that the number of reads from the plant nuclear genome was highly correlated with diploid cell equivalents and thus fresh weight [21,22] . We therefore scaled the non-plant read counts to the number of reads that could be mapped to any of the five A.
thaliana chromosomes. This use of plant chromosomal DNA is analogous to studies that use internal 'spike-in' controls calibrated to sample weight or volume [23][24][25] ; our 'spike-in' is inherent to our samples (Supplementary Figure 1).
About half of non-plant reads in each sample could be assigned to microbial taxa. That the number of non-classifiable reads was positively correlated with the number of microbial reads suggested that these unclassified reads were mostly from portions of microbial genomes not present in the NCBI nr database, rather than being plant sequences not found in the A.
thaliana reference genome (Supplementary Figure 2). To gauge how likely the unclassified reads were to contain sequences missing from the A. thaliana reference genome, we mapped these to additional A. thaliana genomes assembled from long-read data [26] including five genomes available in house. The number of unclassified reads that mapped to additional plant genomes was unrelated to the quantity of unclassified reads in the sample; even in samples with up to 21% unclassified reads, the fraction of reads that mapped to the additional reference genomes was less than 1% of the total classifiable plant reads (Supplementary Figure 2). In other words, across all samples, only a small but consistent percentage of unclassified reads was likely to come from the plant. The rest most likely reflects additional microbial sequences. These sequences may belong to noncoding regions or genes from known taxa that have not been assembled and incorporated into the database. Currently, we cannot easily know how many of the nonclassified, but putative microbial reads reflect highly variable sequences of accessory genomes from known taxa, nor how many reflect the presence of microbial taxa that have not yet had their genomes sequenced. Overall, our results were reminiscent of efforts to classify metagenomic reads from soil and human gut, where more than 50% of reads could not be annotated against known databases [27][28][29] .
To further investigate species-level identification of microorganisms and search for microbial functions, we attempted metagenome assembly of all samples. We assembled short reads with MEGAHIT [30] (meta-sensitive preset parameter), filtered out contigs shorter than 200 bp, and assessed standard assembly metrics such as N50, N90, mean contig length, and total assembly size (Supplementary Figure 3). Additionally, we mapped reads back to their corresponding assemblies to determine what fraction of each library was effectively being incorporated into contigs. The N50 of contigs that were at least 200 bp long ranged from 500 to 800 bp, with the maximum length of individual contigs in the different assemblies ranging from 6 to 12 kb and the sum of their lengths ranging from 4 to 160 Mb. The mapping rate of short reads to their respective assemblies was only 25%. We made a further attempt at assembling reads into contigs using metaSPAdes [31] . We tried assembling individual samples separately as well as pooling them by sampling location to increase coverage. This yielded only modest improvements (Supplementary Figure 4). This difficulty to assemble long contigs was an apparent consequence of high sample diversity (Supplementary Figure 5). It paralleled the limited success in assembling metagenomes of deeply sequenced soil, where despite of having over 300 Gb of data, 80% of sequences could not be assembled because of low coverage of individual taxa [27] .
To evaluate the reproducibility of our approach and to infer potential biases, we prepared independent libraries for 42 samples. We used the same DNA input for 18 samples,

Overview of microbial taxa in the metagenomes
Overall, we found large variability in the fraction of assignable microbial reads, ranging from 3% to 45% of total read counts in each sample (Supplementary Figure 7). We had chosen an arbitrary depth of sequencing for our effort, and we therefore wanted to learn how much information would be lost by reducing the number of sequencing reads per sample. We made use of the replicated individuals to this end. We downsampled reads from replicated libraries and compared their taxonomic profiles. We estimated that approximately 300,000 non-plant reads constitute a lower bound for robust description of taxonomic profiles (Supplementary Figure 18). This agrees well with similar estimates for human gut microbiome samples [11] , and translates into 7.5 million total reads, or just under 1.12 Gb total sequencing reads per sample for 90% of the dataset. As a counterpoint to downsampling reads, we were curious how much could be gained by having much deeper sequence coverage from a single plant. Therefore, we processed a single plant that was visibly infected with white rust ( Albugo spp.) and downy mildew ( Hyaloperonospora arabidopsidis ), and that we in addition left unwashed to further potentially increase the fraction of microbial reads. We sequenced this plant to high depth (~20 Gb), which was 5 to 20 fold more coverage than the other samples. Fewer than 40% of reads from this sample mapped to the A. thaliana reference genome. Similar to our other samples, about half of the remaining reads could be assigned to microbial taxa, with over 90% coming from bacteria (Supplementary Figure 12). In addition to many Albugo spp. and H. arabidopsidis reads, we found many of the bacterial taxa already detected in the other samples, and in similar proportions.

Influence of site, season and host genetics
A common way to compare composition of microbiomes is based on the Bray-Curtis dissimilarity measure. However, a true distance metric is better suited than a dissimilarity measure for other downstream analyses such as principal component analysis [32] . For distance/dissimilarity measures weighted by taxa abundance, highly abundant taxa can strongly skew results while low abundance taxa contribute relatively little information. To evaluate how the microbiomes of our samples relate to each other, we used, instead of Bray-Curtis dissimilarity, pairwise Euclidean distances and PCA. We first transformed the data by taking the fourth root of the family-level abundance table, including all bacterial, fungal, and oomycete taxa. This transformation corrects for positive skewness in count distribution common in ecological datasets [33] , and also decreases the influence of high abundance microbes. No single metadata variable could clearly explain the distributions along the main axes in PCA, although collection site seemed to do best (Figure 2a, 2c, Supplementary Figure   13). Clustering of samples was most clearly driven by the most abundant taxa in each sample, a feature that correlated with collection site. Separation by Pseudomonadaceae or Sphingomonadaceae was apparent when comparing PC1 to PC2, whereas separation by less abundant taxa could be seen when comparing PC2 to PC3 (Supplementary Figure 14).
Finally, we used the abundant plant reads for host genotyping, using FreeBayes [34] to call close to 1 million SNPs from reads with high quality mapping to the TAIR10 reference genome. There were several clear clusters of plant genotypes (Fig. 2c) correlated with sampling site, in agreement with stands of A. thaliana in southwest Germany normally hosting only a limited number of genotypes [15] . It is well known that host genotype can influence the composition of the leaf microbiome [35,36] , but as genotype is strongly linked to site in wild populations, both variables are confounded and would require additional data in order to separate the effects of each.
In an orthogonal analysis, we first classified reads broadly as bacteria, fungi, plants, or unclassified, and compared overall sequence similarity between samples in each class using MASH [37] , which measures similarities in k-mer abundance. MASH does not consider the taxonomic classification of sequences, but because samples containing different taxa include different sequences and hence different k-mers, this classification-independent analysis captured many of the same patterns in the data. PCoA on MASH distances between bacterial, fungal, plant, or unclassified reads also led to some degree of clustering of samples by collection site (Supplementary Figure 15).

Intermicrobial correlation networks
Shotgun sequencing provides a minimally-biased estimation of the true abundance of microbes in a microbiome sample. We examined microbial abundances across all samples and under different data transformations to understand colonization patterns and potential intermicrobial relationships. We first made a map of pairwise linear correlations between all bacterial families that passed filtering thresholds (1,000 assigned reads per family in at least 10 samples) using plant-scaled data (equal plant chromosomal reads), relative abundance data, and fourth root of plant-scaled data. Only cooccurrences with a Pearson correlation coefficient greater than ±0.2 | and with p-value lower than 0.05 after Student's t-test were used (Supplementary Figure 16) . | On average, any taxon was positively correlated with 13 other taxa, but this was heavily skewed toward microbes present in many samples, as correlations between taxa only seen in a handful of plants were usually not significant. In addition, bacterial load varied widely across samples and it was positively associated with the abundance of sequences associated with each taxon (Supplementary Figure 17). This resulted in only positive correlations between taxa (figure 3a). That is, a high bacterial load meant a higher abundance of nearly all taxa, although some taxa such as Pseudomonas contributed more to microbial load than others. When applying fourth root transformation to the plant-scaled microbial counts, the same trend as in untransformed data was observed for all taxa, the only difference being increased correlation values and an increased number of correlated taxa pairs. If bacterial data in all samples are transformed to relative abundance prior to analysis -a necessity for data without internal standards such as most amplicon data -taxa abundance estimations become constrained because the sum of all taxa is constant, greatly confounding the directionality of their correlation. When such a transformation was applied to our dataset, many of the positive correlations between taxa either disappeared or became negative, including the one between Pseudomonadaceae and Sphingomonadaceae ( Figure 3). If we did not have data on bacterial load, it would be tempting to jump on this as evidence for widespread antagonism between these families in the wild. Indeed, antagonism between Pseudomonas and Sphingomonas is known to occur in laboratory conditions [38] . However, if widespread antagonism exists between members of these families, our data did not reveal it. Spurious conclusions due to compositional data are a common and well-documented problem that several analysis methods, such as the centered-log-ratio transformation, attempt to overcome [39][40][41][42][43] .

Concordance between metagenome and amplicon data
In order to contrast information from metagenome and amplicon sequencing, we focused on the largest batch, batch 3, with 176 samples. We PCR amplified and sequenced the V4 region of bacterial 16S rDNA and the fungal ITS1 region for these samples. Because the V4 16S rDNA sequence of A. thaliana -associated cyanobacteria is indistinguishable from that of chloroplasts, reads with cyanobacteria assignments were ignored and cyanobacteria reads were also removed from the metagenome dataset for a fairer comparison. The agreement between assignment of bacterial families based on 16S rDNA amplicons and metagenomes was very high (Fig. 4A), with an overall Pearson coefficient of correlation R 2 of 0.78 on fourth root transformed data ( Figure 4B). Among the top taxa, compared to metagenomics estimates, 16S rDNA estimates were slightly lower for Pseudomonadaceae, and slightly higher for Sphingomonadaceae, Sphingobacteriaceae, and Oxalobacteraceae (Fig. 4A,B). In a complementary comparison, we extracted only the 16S rDNA sequences from the metagenome reads and classified them using phyloFlash [44] . When plotted against 16S rDNA amplicons that had been subsampled to match the metagenome 16S rDNA read counts, the overall correlation and overestimation/underestimation trends were the same as for the comparison of amplicons with all metagenome reads (Fig. 4E, compare to Fig. 4B).
The concordance between fungal families deduced either from ITS1 amplicons or metagenomes was weaker than for bacterial families (Fig. 4F,G), with a Pearson coefficient of correlation R 2 of 0.14. Several factors could explain this difference. First, fungi are less abundant overall, meaning their quantification is based on fewer sequences and therefore noisier. In agreement, the Pearson correlation coefficient R 2 of metagenome versus amplicon data for the most abundant fungal family, Ceratobasidiaceae, was much higher than the average, with R 2 = 0.88. Among other families, the Helotiaceae were especially poorly correlated, and were far more abundant in the ITS1 data (Fig. 4F,G). While this could be due to a bias of the ITS1 primers for this family at the exclusion of others, it could also be that metagenome sequences from this family are more often erroneously assigned to other, sequence related families, deflating counts of Helotiaceae. Another source of noise between the datasets could come from the fact that fungal genomes vary more widely than bacteria in size, and the fact that rDNA copies are not well correlated with fungal genome size [45] . This can introduce biases because organisms with larger genomes would appear to have higher abundances in metagenomes due to more mapped reads. Because these larger genomes may not have more rDNA copies, ITS1 amplicon counts are less affected by differences in genome size.
Additionally, because a much smaller fraction of fungal genomes -as compared to bacterial genomes -codes for proteins, and because gene number in fungi varies much less than genome size, many fungal sequences in the metagenome may not be represented in the protein databases used for classification and quantification. The close concordance between metagenome and 16S rDNA relative abundances enables the scaling of 16S rDNA amplicon data based on bacterial load obtained from metagenome data. Using amplicons, abundant plant host DNA is much more easily blocked using either PNA oligomers [46] or blocking primers [47] , meaning that for samples with high plant DNA content, amplicon sequencing can sample many more microbes at a lower cost than metagenome sequencing. Entries in amplicon databases currently represent the taxonomic breadth of microbes more evenly than whole genome or protein databases, and therefore may provide more consistent classification [48] . In addition, we have shown above that relatively shallow metagenome sequencing is sufficient to find enough classifiable reads with which to estimate the bacterial load of a sample (Supplementary Figure 18). Low-cost shotgun library preparation methods [49] in particular m ake a hybrid approach, in which amplicon sequencing is combined with low-depth shotgun sequencing, attractive. We used the bacterial load as calculated from plant-scaled metagenome data to adjust the abundances of 16S rDNA data to reflect estimated loads (Fig. 3C). As would be expected from the close correlation of 16S rDNA and metagenome data, the adjusted 16S rDNA data accurately captured the slightly positive correlation between Pseudomonadaceae and Sphingomonadaceae across the dataset (Fig. 3D).

Discussion
Describing the phyllosphere-associated microbial community in the context of natural or field cultivated plant populations is of fundamental importance for understanding and designing microbial interventions in conservation and agriculture. For years, as in studies of other microbial communities, this has been approached via isolation and culture of specific leaf microbes [50][51][52] . With the advent of high-throughput amplicon and short read sequencing, it has become easier to address the larger community of taxa that interacts with its host. Here, we investigated the advantages as well as the limits of whole metagenome shotgun sequencing to study the leaf microbiota.
A first finding was that de novo assembly produced few longer sequences, and thus did not provide advantages over directly mapping short reads to reference databases for taxonomic assignment. Using taxonomic assignments inferred from short reads directly, we found that the relative abundance of bacterial taxa among all bacteria was highly consistent with 16S rDNA amplicon measurements, with the highest correlation seen for the most abundant taxa. The relative abundance of fungi in the metagenome correlated less well with the ITS1 amplicons, which could be explained at least in part by their lower abundance as well as more complex genomes of fungi compared to bacteria. Nevertheless, the metagenome data clearly showed that fungi are ubiquitous in A. thaliana leaves, even though they are usually only a minor part of the overall A. thaliana phyllosphere microbiome.
We used shotgun sequencing to estimate microbial load across samples, which varied substantially. The absolute estimates allowed us to reveal the extent to which normalizing bacteria or fungi to a common value via rarefaction or by total sum scaling, as is common practice, may mislead researchers to equate an increase or decrease in relative abundances with a change in absolute abundances.
We used load corrected bacterial taxonomic profiles to explore similarities and differences among the microbiomes. We did not detect a strong individual influence of either environment or host genetics on the structure of leaf communities, although some combination of both contributed. We found that the most abundant taxa in a sample predicted the community structure of the other microbes in the sample. In our natural host populations, there was not enough genetic diversity at each site to allow us to disentangle the effects of site and host genotype on microbiome composition. Whether variation between genetically identical hosts reflects only stochastic effects, or variation in microenvironment, needs to be determined.
An important goal of microbial ecology is to uncover specific interactions between community members, which can point to key taxa that have a major effect on community composition [53,54] . Correlation networks present a valuable tool for such investigations; we have demonstrated that correlations based on relative abundances can lead to very different networks than correlations based on absolute abundances. In our specific case, relative abundance data had suggested that Pseudomonadaceae and Sphingomonadaceae, two of the most common bacterial families typically found in a leaf microbiome, were negatively correlated, when in fact they were positively correlated, as demonstrated with the absolute abundance data. We observed such patterns for several other taxa pairs as well. It is striking that statistically significant correlations were always positive in our dataset when using load corrected data. This could simply reflect generally stable relationships between community members, such that they tended to succeed or fail together as they colonized the plant.
In conclusion, we have demonstrated the advantages of using metagenome shotgun sequencing either alone or in combination with 16S rDNA and ITS1 amplicon sequencing for measuring microbial communities in A. thaliana leaves. Modest read depth, as few as 300,000 reads per sample, is sufficient to enable quantitative taxonomic assignment that is comparable to amplicon sequencing. In addition, it turns out once more that the small genome of A. thaliana is a substantial advantage, as it may currently be cost prohibitive to extend our direct metagenomic approach to other species. This is yet another reason to use A. thaliana (or other species with relatively small genomes) for microbiome studies in ecological settings. On the other hand, the hybrid approach of using lower coverage metagenome data to estimate microbial load and to use this information to scale amplicon data may cost effectively support endophytic analysis of plants with larger genomes as well.

Sampling, Processing of Plants, Metagenomic Library Preparation
Plants were sampled from previously-described sites Eyach (EY), Pfrondorf (PFN) and Jugendhaus Einsiedel (JUG) around Tübingen, Germany [15,55] , in four distinct sampling batches which also had different processing details, representing our evolving pipeline. was used for DNA extraction, using a custom protocol we previously described [55] . Briefly, the sample was subjected to bead-beating in the presence of 1.5% sodium dodecyl sulfate (SDS) and 1 mm garnet rocks, followed by SDS cleanup with ⅓ volume 5 M potassium acetate, and then SPRI beads. The library was prepared using the TruSeq Nano kit (Illumina), with DNA shearing performed with a S2 focused ultrasonicator (Covaris) as suggested in the manufacturer's protocol. Rather than Illumina adapters, we used custom adapters described in [56] . End-repair, A-tailing, and adapter ligation were performed similar to [57] following "Alternative Protocol 2" with double DNA size selection after the End-repair step, and using homemade SPRI beads instead of AMPure XP beads. Other minor modifications were that the total volume of the end repair reaction was scaled down to ¼ volume, with DNA eluted after SPRI-cleanup in 17 μ L. The total volume of the A-tailing reaction was scaled down to ½ volume. Again, custom adapters described in [56] were substituted for Illumina adapters. were removed with sterile scissors and tweezers, and washed 3x with sterile water. Two leaves were removed and independently processed to culture bacteria as previously published in [55] , and the remaining rosette was flash-frozen on dry ice and processed for metagenomic sequencing and 16S rDNA sequencing of the V4 region. The metagenomic libraries were prepared using a modification of the Nextera protocol for smaller volumes similar to [58] , as previously described [55] . As for Batch 1 and Batch 2 plants, the full set of 176 libraries was quantified and combined to one pool, size selected for 350 -700 bp final library size, and sequenced with 2x150 paired ends over multiple lanes of a HiSeq3000 instrument.

16S rDNA V4 library construction and sequencing for Batch 3 plants
Amplicon sequencing of the V4 region of the 16S rDNA gene performed using a two step PCR protocol using PNAs to block chloroplast and mitochondrial sequences, slightly modified from [46] . The first PCR step amplified the rDNA using 515F [59] and 806R [60] primers as well as short overhangs (Supplementary Table 1) and a second step primed these overhangs to added Illumina adapters. The primers differed from Lundberg et al. 2013 in two key ways. First, although the frameshifting nucleotides were kept, the molecular tagging nucleotides were removed from the primers to make the protocol simpler and robust to more variable DNA quantities. Second, the primers were modified such that the Illumina TruSeq priming sequences were used on both the forward and reverse primers, as opposed to the use of a Nextera sequence on the forward primer. Unique barcoding of samples was accomplished by use of 96 independent indexing primers in the second PCR, combined with two combinations of frameshift primers in the first PCR as explained in [46] . Half of the samples from the first PCR were amplified with 515F frameshifts 1, 3, and 5 paired with 806R frameshifts 2, 4, and 6. The other half of the samples from the first PCR paired 515F frameshits 2, 4, and 6 with 806R reverse frameshifts 1, 3, and 5. This strategy allowed up to 192 samples to be uniquely indexed.
In the first PCR, each reaction was prepared in 60 μL, which was split into three 20 μL amounts, and sequenced using a MiSeq V2 2x500 reagent kit (Illumina) which was sufficient to overlap and assemble the forward and reverse reads. The frameshifts built into the primers used in the first PCR made the addition of PhiX to increase sequence diversity unnecessary [46] .

ITS1 library construction and sequencing for Batch 3 plants
ITS1 rDNA amplicons were prepared similarly to 16S rDNA amplicons, using gene-specific primers for the first PCR and adding indexes and adapters in the second PCR. We used a protocol modified from [47] , which uses blocking primers to prevent amplification of plant sequences. Because blocking primers, unlike PNAs, result in a quantifiable PCR product, we used the cycling conditions suggested in [47] to prevent the product of the blocking product from reaching detectable levels. As with the 16S rDNA protocol, we used six frameshifting forward ITS1F primers and six frameshifting reverse ITS2R primers. The first 60 μL PCR reaction (also run as three parallel 20 μL reactions) included 6 μL of 10X ThermoPol Taq

Amplicon quality processing, clustering, and classification
Raw sequences from both 16S and ITS1 rDNA amplicons were first demultiplexed according to their 9 bp barcodes added in the second PCR, not allowing any mismatches. All sequences were further demultiplexed by the frameshift combinations using strict regular expressions without mismatches in any part of the primer sequence ( https://github.com/derekLS1/Metagenome2019 ). Forward frameshifts 1, 3, and 5 were only allowed pairings with reverse frameshifts 2, 4, or 6. Forward frameshifts 2, 4, and 6 were only allowed pairings with reverse frameshifts 1, 3, and 5.
Forward and reverse reads from the 16S rDNA sequences were merged with FLASH

Metagenome read QC and host data removal
Sequencing libraries were subject to adapter trimming and quality control with Skewer [65] .
Reads were trimmed to a minimum length of 30 bp and minimum average Phred score of 20.
After sequencing, samples were composed of a mixture of mostly host Arabidopsis thaliana reads and microbial origin reads. In order to remove most of the plant reads, libraries were aligned against the A. thaliana reference genome [16] using the bwa mem algorithm with standard parameters [66] . After mapping, only read pairs for which neither of the mates mapped against the plant reference genome were mapped against the metagenomic reference.
Data aligned to the host was later used for host plant genotyping.

Taxonomic correlations and network computation
After computing intermicrobial linear correlations and filtering out weakly associated taxa pairs, network representation was computed with the networkx Python package using the kamada kawai layout function with standard parameters based on the correlation values.

Plant genotyping
Individual host genotypes were determined using plant associated reads from each metagenome. Reads aligned to the TAIR10 reference genome [16] were filtered to a minimum mapping quality of 20, resulting in an average genome coverage from 15x to 40x. Single nucleotide polymorphisms were called using FreeBayes [34] , and resulting VCF files were filtered using custom scripts. SNPs with a minimum alternative count above 3, minimum read depth of 6, and no more than 5% missing data across all samples were kept for downstream analysis.
For determining genotype groups, a genetic distance matrix was computed with ngsDist [67] from the alternative allele count matrix of all SNPs that passed filtering thresholds. This distance matrix was used as input in tsne [68] to visualize sample clustering.

Adjusting 16S rDNA amplicon data by bacterial load
To adjust the 16S rDNA amplicon dataset to correct for bacterial load, the abundance of each OTU from a sample in the total sum-scaled OTU table can be multiplied by a load scaling factor calculated from the metagenome data for that sample. The simplest load scaling factor is the ratio of all bacteria to plant chromosomal reads in the metagenome sample. If read depth in the metagenome allows, a more precise scaling factor can be calculated based on bacterial families detectable by both methods. Both methods yield similar results in our dataset, because the majority of sequencing reads in both methods fall into bacterial families shared by both methods. We scaled based on bacterial families shared by both methods.
To correct the 16S rDNA dataset by shared taxa in the scaled metagenome dataset, the 16S rDNA dataset was first normalized to 100% in each sample by total sum scaling. The common bacterial families that could be identified by at least a single read in both the metagenomic and 16S rDNA datasets were then identified for each sample, and the sum of read counts in falling in these common taxa was calculated for each sample in both the metagenome and 16S rDNA datasets. The sum of common taxa for each sample in the plant-chromosome scaled metagenome dataset was divided the sum of reads in these common taxa in the corresponding 16S rDNA table to yield a load scaling ratio. Metagenome reads from each sample in the Batch 3 dataset were mapped to the RDP 16S rDNA training set using phyloFlash [44] . Each metagenome sample yielded an average of 280 mapped 16S rDNA reads, with many yielding fewer than 100 reads. Because for most samples there were too few reads to compare directly to their 16S rDNA amplicon counterparts, samples were pooled to make one aggregate metagenome dataset containing 52,589 16S rDNA phyloFlash sequences. This was then compared to a corresponding 16S rDNA amplicon dataset comprising 52,589 sequences subsampled from the full 16S rDNA dataset. Each sample contributed a matching number of phyloFlash 16S rDNA or amplicon 16S rDNA reads to either the phylFlash 16S rDNA pool or the amplicon 16S rDNA pool respectively. The family level relative abundances for these pools were then plotted against each other.

Data availability
All data in this manuscript has been deposited in the European Nucleotide Archive (ENA). It

Competing Interests
OD is also an employee of Computomics GmbH. The other authors declare no competing interests.