Understanding Cultivar-Specificity and Soil Determinants of the Cannabis Microbiome

Understanding microbial partnerships with the medicinally and economically important crop Cannabis has the potential to affect agricultural practice by improving plant fitness and production yield. Furthermore, Cannabis presents an interesting model to explore plant-microbiome interactions as it produces numerous secondary metabolic compounds. Here we present the first description of the endorhiza-, rhizosphere-, and bulk soil-associated microbiome of five distinct Cannabis cultivars. Bacterial communities of the endorhiza showed significant cultivar-specificity. When controlling cultivar and soil type the microbial community structure was significantly different between plant cultivars, soil types, and between the endorhiza, rhizosphere and soil. The influence of soil type, plant cultivar and sample type differentiation on the microbial community structure provides support for a previously published two-tier selection model, whereby community composition across sample types is determined mainly by soil type, while community structure within endorhiza samples is determined mainly by host cultivar.


Introduction
Soil microbes play a major role in plant ecology by providing a variety of benefits such as nitrogen fixation, production of growth stimulants, improved water retention, and suppression of root diseases [1][2][3][4]. These vital microbial processes occur predominantly within the rhizosphere and rhizoplane, and are heavily influenced by fungal saprotrophs and plant-mutualists such as endomycorrhizal and ectomycorrhizal fungi [5,6]. Despite the economic and medicinal importance of Cannabis spp., little is known about its soil-based microbial associations [7,8].
Microbial composition in soil depends on complex interactions between the soil type, root zone location, and plant species [9][10][11]. Rhizosphere microbiota are highly dynamic [12], and the composition of bacterial communities can fluctuate in response to seasonal and diel temperature changes [13], water content [14], pH [15], CO 2 concentration, and O 2 levels [16]. Although evidence has been found for significant effects of plant cultivar on rhizosphere communities [17][18][19] and endomycorrhizal fungal communities [20], some work suggests that these effects are minimal compared to edaphic factors (particularly pH) or plant growth stage [21,22].
Rhizosphere bacteria not only colonize the rhizosphere and/or the rhizoplane soil, but can also colonize plant tissues. Bacteria that have colonized root tissue-more specifically known as the endorhiza [23]-have been reported to support plant growth and suppress plant diseases by providing phytohormones, low molecular weight compounds or enzymes involved in regulating growth and metabolism [24][25][26]. In addition, endorhiza bacteria assist their host plants in tolerating the phytotoxic effects of environmental toxicants [27,28]. Endorhiza communities tend to be more plant-specific, and are often shaped by the compounds or proteins produced by their host [29]. Both endophytes and epiphytes may also play a role in localized 'flavor' or terroir for crop plants, as has been shown recently for wines [30][31][32].
A growing body of work has united the colonization of both the rhizosphere and plant tissues under the two-tier selection model, where soil type defines the composition of rhizosphere and rootinhabiting bacterial communities [33][34][35]. Under this model, edaphic factors determine the structure of the local soil microbiota, which become the source for the first bacterial community shift into the nutrient rich environment of the rhizosphere. Following this first shift, migration from the rhizosphere into the plant tissues is based on plant genotype-dependent selection of the endorhiza environment [33]. Along with the prediction that rhizosphere and endorhiza microbiota should be soil-derived, the two-tier selection model predicts several broad changes in phylum-level taxon abundance associated with the shifting microbiota, such as dramatic reduction in Acidobacteria within the endosphere.
This study aims to characterize bacterial diversity in the root and soil systems of five strains of Cannabis in order to explore how soil microbiota and plant strain affect the endorhiza microbial community of this commercially important crop. We hypothesize that different cultivars maintain significantly different microbial communities, and that these differences diminish from endorhiza to rhizosphere to bulk soil.

Experiments
The data for this paper were collected in two experiments: First, an experiment to identify variation in the microbial communities, and second, an experiment designed to understand the nature and strength of cultivar-specificity. The first experiment was composed of bulk soil, rhizosphere, and endorhiza samples taken from nine plants of the three different Cannabis spp. tested strains-Burmese, BooKoo Kush, and Sour Diesel. Soil physicochemical data was taken for all bulk soil samples in the first experiment, however there was minimal edaphic variation. The second experiment sought to understand the effect of strain with more significant edaphic variation, and was accomplished using two different strains-White Widow and Maui Wowie-and two different soil types. Four plants of the two strains were grown in the same soil, and then two plants of White Widow were grown in a completely distinct soil type. Triplicate samples were taken from each plant for both the rhizosphere and endorhiza, as well as for each of the two soil types.

Cultivars
Different cultivars were used for each one of the experiments. For the first experiment, we used Sour Diesel, Bookoo Kush, and Burmese cultivars. Sour Diesel is a cultivar of Cannabis sativa, associated with a high tetrahydrocannabinol (THC) to cannabidiol (CBD) ratio. Bookoo Kush is a sativa-dominant hybrid cultivar of Cannabis sativa and Cannabis indica, associated with a moderately high THC to CBD ratio. Burmese is a balanced hybrid cultivar of both Cannabis sativa and Cannabis indica, associated with a moderate THC to CBD ratio. For the second experiment, we used Maui Wowie and White Widow cultivars. Maui Wowie is a cultivar of Cannabis sativa, associated with a high THC to CBD ratio. White Widow is a balanced hybrid cultivar of both Cannabis sativa and Cannabis indica, known to have a more moderate THC to CBD ratio.

Sample Collection
Endorhiza, rhizosphere soil, and bulk soil samples for the first experiment were taken from 9 organically-grown Cannabis plants of three different strains (Burmese, Bookoo Kush, Sour Diesel) in Vista, California, in November, 2011, for a total of 27 samples. Therefore, the triplicate DNA extracts were acquired for endorhiza, rhizosphere and bulk-soil for each of the 3 Cannabis spp. strains, resulting in a single endorhiza, rhizosphere, and bulk soil sample for each plant. The plants were grown in locally composted soil. Eight weeks following the harvesting of the Cannabis flowering bud and foliage from each plant, a 50 g bulk soil sample was taken 10 cm from the stem of each of the nine plants at a depth of 20 cm, as well as a larger sample of soil for testing edaphic factors ( Table 1). The bulk soil sample was immediately capped and transported to a 4uC refrigerator. In addition, endorhiza samples were taken from the root ball of each of the six plants. The soil that remained adhered to the roots after removal from the ground was used to produce the rhizosphere soil samples. The rhizosphere soil was removed from the roots by shaking the root into a whirlpak bag. All samples were immediately transferred to storage at 4uC for shipping back to the laboratory for processing (approximately 4 hours). All root samples were rinsed with alcohol and sterile water before the extraction. DNA was isolated from 0.25 g of soil or root per extraction using standard protocol for PowerSoil DNA Isolation Kit (MO BIO, USA), with the modification of heating the extraction at 65uC for 10 minutes prior to the initial vortex step. The soil physicochemical data was generated by Fruit Growers Laboratory (Santa Paula, CA), including total carbon and nitrogen concentrations, pH, salinity, and water content for all samples.
Endorhiza, rhizosphere, and bulk soil samples for the second experiment were taken from 6 organically-grown Cannabis plants of two different strains (White Widow and Maui Wowie) from two locations in August, 2012: Vista and Orange County, California. Triplicate samples were taken from each of the six plants (18 samples) and surrounding rhizosphere (18 samples), as well as from each of the two bulk soils used in the different locations (6 samples), totaling 42 samples. In contrast to the first experiment, all samples were taken two weeks prior to harvest. Additionally, triplicate samples from the second experiment were taken from different roots on the same plant (pseudoreplicates). Cannabinoid data was taken from the buds of three White Widow plants and one Mauie Wowie plant ( Table S1). All cannabinoid data was processed at Delta-9-Technologies, LLC (Santa Ana, California). Otherwise, sampling procedure matched the first experiment.

Illumina sequencing of the V4 region of the 16S rRNA gene
We utilized Illumina 16S rRNA sequencing to analyze samples of the endorhiza, the rhizosphere, and the bulk soil of three different strains of Cannabis in the first study (27 samples), and two different strains of Cannabis in the second study (42 samples), for a total of 69 samples. The V4 region of the 16S rRNA gene was amplified and sequenced using the primers specified in Caporaso et al. (2012) following the Earth Microbiome Project's standard pipeline (http://www.earthmicrobiome.org/emp-standardprotocols/) [36]. The 291 bp length V4 region amplification was performed using the 515F primer and the 806R Golay-barcoded reverse primers (for a full list of these primers visit http://www. earthmicrobiome.org/emp-standard-protocols/). Each 25 mL PCR reaction contained 12 mL of MO BIO PCR Water (Certified DNA-Free), 10 mL of 5 Prime HotMasterMix (1x), 1 mL of Forward Primer (5 mM concentration, 200 pM final), 1 mL Golay Barcode Tagged Reverse Primer (5 mM concentration, 200 pM final), and 1 mL of template DNA. The conditions for PCR are as follows: 94uC for 3 minutes to denature the DNA, with 35 cycles at 94uC for 45 s, 50uC for 60 s, and 72uC for 90 s, with a final extension of 10 min at 72uC to ensure complete amplification. PCR was completed in triplicate and products were pooled. Each pool was then quantified using Invitrogen's PicoGreen and a plate reader. Once quantified, different volumes of each of the products were pooled into a single tube so an equal amount (ng) of DNA was in the pool, and cleaned using the UltraClean PCR Clean-Up Kit (MO BIO). After quantification, the molarity of the pool is determined and diluted down to 2 nM, denatured, and then diluted to a final concentration of 6.1 pM with a 30% PhiX spike for sequencing on the Illumina MiSeq. A 151 bp612 bp6151 bp MiSeq run was performed using the custom sequencing primers and procedures described in the supplementary methods in Bioinformatic analysis of the 16S rRNA V4 sequence data All sequence analysis was done using QIIME 1.7.0 [38]. QIIME defaults were used for quality filtering of raw Illumina data. In the second study, both closed and open reference OTU-picking methods were employed. In the first study, OTUs were picked against the Greengenes [39] database pre-clustered at 97% identity, and sequences that did not hit the reference collection were clustered de novo (i.e. open reference). Representative sequences were aligned to the Greengenes core set with PyNAST [38]. All sequences that failed to align were discarded. A phylogenetic tree was built from the alignment using FastTree [40], and taxonomy was assigned to each sequence using the RDP classifier [41] retrained on Greengenes. Samples for the first experiment were rarified to an even depth of 3,000 sequences. Four samples were discarded due to insufficient sequence coverage. For the second experiment, samples were rarified to an even depth of 45,000 sequences. One sample was discarded due to insufficient coverage. Alpha, and beta-diversity metrics were produced using QIIME [38]. Relationships between samples were visualized and evaluated using redundancy analysis (RDA) and principal coordinate analyses (PCoA) calculated from pairwise sample distances (weighted and unweighted UniFrac metrics) [42]. Significance tests were run using the compare_categories.py (ANOSIM, ADONIS, ANOVA, and RDA) and compare_distan-ce_matrices.py (Mantel) scripts in QIIME [38]. To evaluate the most important abiotic factors in structuring the communities, a Best Subset of Environmental Variables with Maximum (Rank) Correlation with Community Dissimilarities (BEST) analysis was run in QIIME (see vegan::bioenv) [43].

Results
Work for this study was accomplished in two experiments. First, we performed an experiment to identify variation in the microbial communities between roots and soil in three different Cannabis strains (Burmese, BooKoo Kush, and Sour Diesel), and second, an experiment designed to understand the nature and strength of plant cultivar-specificity between two different strains (White Widow and Maui Wowie) in two different soil types (with significant differences in edaphic variables). Triplicate samples were taken from each plant for both the rhizosphere and endorhiza, as well as for each of the two soil types.
Both endorhiza and bulk soil microbiomes were significantly distinct from other sample types, and strain level differences were only observed in the endorhiza In the first experiment, using unweighted UniFrac, betadiversity comparisons of each individual sample type against all other sample types (Fig. 1a) yielded significant clustering of endorhiza (ADONIS: R 2 = 0.26, p = 0.001) and bulk soil (ADO-NIS: R 2 = 0.14, p = 0.001) samples from the other categories, but rhizosphere samples were not significantly different (ADONIS: R 2 = 0.07, p = 0.07). Weighted UniFrac distances yielded similar results with endorhiza (ADONIS: R 2 = 0.59, p = 0.001) and bulk soil (ADONIS: R 2 = 0.29, p = 0.004) samples demonstrating significant differences from other sample types, but no significant differences for rhizosphere (ADONIS: R 2 = 0.09, p = 0.10) samples. Division of all communities via strain (Fig. 1b) was not significant for weighted (ADONIS: R 2 = 0.11, p = 0.25) or unweighted (ADONIS: R 2 = 0.11, p = 0.15) analyses, however, division of endorhiza communities via strain was significant for both weighted (ADONIS: R 2 = 0.59, p = 0.004) and unweighted (ADONIS: R 2 = 0.39, p = 0.003) analyses. The abundance of Methylophilus explained a significant portion of this difference (FDR: p = 0.012), comprising 13% of the microbial community in the endorhiza of Bookoo Kush, 0.13% in Burmese and was absent in Diesel. Despite these significant differences, all endorhiza samples maintained a core community of Pseudomonas, Cellvibrio, Oxalobacteraceae, Xanthomonadaceae, Actinomycetales, and Sphingobacteriales. With the exception of the aerobic cellulytic bacterium Cellvibrio, all prevalent members of the core endorhiza community were well known endophytic bacteria [44,45] primarily within the orders Gammaproteobacteria and Alphaproteobacteria, which supports observations from other plant systems [46,47].
Community composition across all samples was determined predominantly by soil properties, but differences in community structure (abundance) within endorhiza were driven by Cannabis cultivar In the second experiment, using unweighted UniFrac, community beta diversity was significantly different between soil types (  (Fig. 2b) yielded significant differences for endorhiza (ADONIS: R 2 = 0.10, p = 0.001) and rhizosphere (ADONIS: R 2 = 0.05, p = 0.04) samples, but no significant differences for bulk soil (ADONIS: R 2 = 0.04, p = 0.12) samples. Using weighted UniFrac, community beta diversity varied significantly by soil type (Fig. 2d) (ADONIS:  Soil followed by strain had the largest affect on OTU abundances, but strain showed no impact on OTU presence/absence For individual OTUs, both unweighted (g-test) and weighted (ANOVA) analyses showed that soil type had the strongest influence over significant OTU differences ( Table 2). While strain showed a larger effect than sample type for weighted OTU differences, there were no significant unweighted OTU differences between strains, further suggesting the importance of strain in structuring OTU abundances -rather than OTU presence/ absence.
As suggested by the two-tier model [33], our results demonstrate a decrease in abundance of Acidobacteria and an increase of Proteobacteria and Actinobacteria relative to the rhizosphere and bulk soil (Fig. 4). Furthermore, the most significant OTU abundance difference between sample types was the decrease in Acidobacteria from the order iii1-15 in endorhiza samples (Bonferroni-corrected ANOVA: p = 1.12e-7). Of the 51 OTUs significantly differentiating between sample types, the 17 OTUs which increased in abundance within the Cannabis endorhiza relative to rhizosphere were predominantly Proteobacteria, including several from the Rhizobiales order. Mean abundance of the 51 OTUs were highly correlated between bulk soil and rhizosphere samples (Pearson's rho: 0.92), versus a lower correlation between rhizosphere and Cannabis endorhiza (Pearson's rho: 0.63), and even lower between bulk soil and Cannabis endorhiza (Pearson's rho: 0.42).
Significant OTU abundance differences between strains were composed mostly of differences in Proteobacteria, notably Pseudomonadales, Burkholderiales, Sphingomonadales, and Rhizobiales. Apart from Proteobacteria, Bacteroidetes orders Sphingobacteriales and Flavobacteriales were also responsible for several significant OTU differences between Cannabis strains. Intriguingly, one of the significant OTUs between strains was the prevalence of Sphingomonas wittichii in the Maui Wowie strain, which in some contexts can metabolize phenazine-1-carboxylic acid and has been implicated in increased survival in soil environments.
Bulk soil and rhizosphere microbiomes are more similar to each other than to endorhiza microbiomes ( Fig. 5) Beta distances between rhizosphere and bulk soil communities were significantly lower than distances between rhizosphere and endorhiza communities for both unweighted and weighted analyses (unweighted: t = 24.59, p,0.001 weighted: t = 211.82, p,0.001). Beta distances between rhizosphere and bulk soil communities were significantly lower than distances between bulk soil and endorhiza communities for both unweighted and weighted analyses (unweighted: t = 25.15, p,0.001; weighted: t = 211.56, p,0.001). Beta distances between rhizosphere and endorhiza communities were not significantly different from distances between bulk soil and endorhiza communities for both unweighted and weighted analyses (unweighted: t=22.10, p = 0.109; weighted: t = 22.23, p = 0.078).
Endorhiza share more OTUs with the soil they are grown in than with another soil in which the same strain is grown Previously, a two-step process of root colonization, first fueled by rhizodeposition and followed by fine-tuning by host genotype has been posited [34,35,48]. This two-step selection model was tested by pooling samples by Cannabis strain and analyzing the core microbiome within each strain. As bulk soils are the putative source of microbes for the plant, endorhiza communities would be expected to share more OTUs with their own soil than with another. White Widow was grown in two different soils, and roots shared more OTUs with the soil they were grown in than with the different soil in which the different white widow plant was grown. The number of shared OTUs between endorhiza and their own soil (n = 45, mean = 2934) was significantly greater (t = 210.05, p = 1.209e-15) than the number of shared OTUs between endorhiza and the other soil (n = 45, mean = 2162).

Cannabinoid concentration and composition was significantly correlated to structure of endorhiza communities
Each plant in the second experiment was tested for a variety of cannabinoids, including delta-9-tetrahydrocannibinol. Cannabinoid data associated with the plants was used in Mantel tests to understand the potential biochemical associations with community composition or structure; with significant differences between strains (unweighted; r-stat: 0.863, p-value = 0.001). However, due to higher THC composition and concentration in plants from one of the soil types, THC variables were also significantly correlated to the soil edaphic variables, and as such any association between microbiota and THC is very hard to disassociate from soil physicochemical variables.
Edaphic factors were strongly linked to structure of microbial communities in rhizosphere and endorhiza communities In both experiments, edaphic data associated with the plants was used in Mantel tests to understand the effect of edaphic factors on structuring bulk soil, rhizosphere, and Cannabis endorhiza communities. For all experiments, the soil texture was defined as a sandy loam, with significant differences in clay and other edaphic unweighted UniFrac distances, and demonstrate relationship between soil type (A), sample type (B), strain (C), and the major PC axes (PC 1 = 32.06% variance, PC 2 = 11.34% variance, PC 3 = 5.67% variance). Plots for weighted analysis are based on weighted UniFrac distances, and demonstrate relationship between soil type (D), sample type (E), strain (F), and the major PC axes (PC 1 = 34.51% variance, PC 2 = 25.41% variance, PC 3 = 19.31% variance). Note that PC 1 in the unweighted analysis is dominated by variation in soil type (A), but PC 1 in weighted analysis is dominated by strain (F). Grey points (Fig. 2c, 2f)   Alpha diversity peaks in bulk soil and declines with the transitions into the rhizosphere and endorhiza microbiomes   Composition of endorhiza communities in the first experiment suggest potential root decay After analysis of the data from the first experiment yielded high abundances (greater than 10% of taxonomy assigned reads) of the known cellulytic bacterium Cellvibrio, we sought to investigate the possibility that Cellvibrio was an indication of root decay rather than its unexpected presence as a member of the endophytic core community. This was of particular interest because the samples were taken 8 weeks post-harvest. Comparisons of relative abundances of Cellvibrio between the first and second experiments yielded rather convincing results demonstrating the early stages of root decay despite significant cultivar-specificity within the samples. Specifically, the relative abundance of Cellvibrio within the endosphere of the first experiment was 16.9% (s = 13.0%, N = 9) versus 0.095% (s = 2.7%, N = 18) in the endosphere of the second experiment.

Discussion
Recent literature has suggested a two-step selection model for the endorhiza, where bulk-soil microbial communities are filtered by increased concentration of rhizodeposits, followed by convergent host genotype-dependent selection on endophytic communities [34,35,48]. Results from both experiments support many of the expectations produced by this model. Most importantly, the principal coordinate analysis (PCoA) plots for the second experiment demonstrate highly significant clustering patterns. First, soil type is the main determinant of PC1 (32.06%) for the unweighted analysis of the second experiment, revealing that soil is undoubtedly the most important factor in all samples for determining what microbes are present. Second, communities within both soil types demonstrate a similar community shift from bulk soil to endorhiza samples along PC2 (11.34%), which is dominated by differentiation between sample types. Specifically, endorhiza samples have high, positive values along PC2, rhizosphere samples have intermediate values, and bulk soil samples have more negative values. Third, Cannabis strain is the main determinant of PC1 (34.51%) for the weighted analysis of all samples in the second experiment, suggesting that convergent host genotype-dependent selection acts through controlling community structure (abundance) more than composition. PCoA results exhibit how all sample types form significantly differentiated clusters in weighted analyses but that only rhizosphere and endorhiza samples form significantly differentiated clusters in unweighted analyses, suggesting niche-filtering of microbes in rhizosphere and endorhiza samples from bulk soil. Furthermore, there were no significant segregating OTUs based on unweighted analysis between cultivars in endorhiza and rhizosphere samples in the second experiment, however there were 71 when abundance was accounted for. This differs greatly from the 657 OTUs that significantly differ between soil types in the same dataset. Testing of the two-step selection model with pairwise comparisons of shared OTUs between endorhiza and bulk soil samples also validated the hypothesis that a portion of the endophytic microbes are inherited and selected from the surrounding soil, showing significantly more OTU overlap between endorhiza and their own bulk soil compared to endorhiza and foreign bulk soil.
Given the results from the second experiment strongly suggesting that Cannabis cultivars have important structuring effects on both rhizosphere and endorhiza samples, it may seem troubling that results from the first experiment do not suggest this for the rhizosphere samples. However, differences in Cellvibrio abundance between experiments show that root decay could have diminished the rhizosphere effect, thus diminishing this potential signal. Sampling for the first experiment was done post-harvest, when plant tissues were undergoing senescence and decay, while samples for the second experiment were taken from actively growing plants. Considering the extensive work demonstrating the importance of plant growth stage on the microbiota [21,49], as well as the plant-soil feedbacks identified in structuring belowground microbial communities [50,51], the differences between the first and second experiments are unsurprising. The similarities, however, are surprising. In particular, that cultivar-specificity could be identified in the microbiota within the endorhiza samples in the first experiment without any input of cultivar-specific metabolites from the living plant for weeks.
Although we have presented several highly significant findings supporting expectations of the two-step selection model, some expectations remain to be validated. Specifically, although the mean beta-diversity distances indicate that rhizosphere and endorhiza samples are closer than bulk soil and endorhiza samples, this difference was not significant and thus provides little evidence for the first differentiation step of the two-step selection model [34,35,48].
Future work with the Cannabis microbiome should focus on elucidating the role of cultivar on rhizosphere, as well as what aspects of host genotype are producing the structure observed across Cannabis strains. Increased testing of cannabinoids and decoupling this variation from edaphic factors will improve our understanding of the importance of cannabinoid production in structuring endorhiza communities. Sampling a time series of endorhiza communities across several plants may help us to understand natural variation in the endorhiza during the reproductive cycles of Cannabis. Understanding this natural variation will help direct future mechanistic studies aimed at using microbial communities to increase plant fitness, suppress disease, or augment desired metabolite production.