The four populations utilized in this experiment were derived from two initial crosses made in the winter of 2015. The parents used in crosses paired a Modern, Elite Snap Bean To An Older Cultivar With Good Performance In Organic And Low-Input Conditions, As Described In Supplementary Table 1. All the parental populations were determinate (bush) beans and produced medium sized pods. The elite parent Hystyle has resistance to Bean Common Mosaic Virus (BCMV), Bacterial Brown Spot, and Curly Top Virus, while the elite parent OR5630 has resistance to BCMV, alone (Cornell University, 2012). The older cultivar Provider is resistant to BCMV, while the older cultivar Black Valentine was recently reported to have some resistance to the root rot complex found in western Oregon (Huster et al., 2021).

The two F₁ generations were grown at the Vegetable Research Farm in Corvallis, OR in the summer of 2015, and the four research populations were created from 500 seeds of bulked F₂ seed. Beginning in the summer of 2016 the populations were split into either organic or conventional management and followed a series of generation advances with only natural selection, described in Supplementary Table 2.

Snap beans were planted on the Lewis Brown Research Farm at latitude of N44.555881 and longitude of W123.214569 into a Chehalis silty clay loam. This field has been certified organic for over a decade and was previously planted to barley (Hordeum vulgare) in 2019. To minimize variation in the seed microbiome, all seeds used for this study were obtained from plants grown in the Oregon State University greenhouses in January of 2020. Treatment in the greenhouse included fertilizer applications of organic and conventional soils and fertilizers, in alignment with the breeding history of the accessions.

A subset of 40 F₇ snap bean accessions derived from the four research populations along with the four parents were grown. The snap beans were planted on May 21, 2020 and root samples were collected the first week of July, which correlated with flowering and early pod set. The field was managed organically, regardless of the breeding history of the individual accessions. Beans were planted in 3m (10 ft) plots, with 30.5 cm (12 in) spacing between plots and 76 cm (30 in) spacing between rows. Plots were fertilized on June 1, 2020, with Nutri-Rich fertilizer (8-2-4 N-P-K) at a rate of 1,159 kg ha^-1. Following the fertilizer application, plots were thinned to approximately 15 cm spacing between plants. Mechanical cultivation was used to control weeds throughout the growing season. Plots were irrigated approximately once a week with the application of about 2.5 cm water through solid-set overhead sprinklers.

Soil samples were collected in April of 2020, prior to fertilizer applications. These samples were sent to the Soil Health Laboratory at Oregon State University and tested for general soil quality and nutrient composition. A second set of soil samples was collected one week prior to root sample collection. These soil samples were tested for general quality (i.e., macronutrient and micronutrient content, organic matter, etc.), and soil microbial assessment.

In July 2020 root samples, along with the associated rhizosphere soil, were collected from the subset of 40 RILs along with the four parents. Two days prior to sample collection the field was irrigated to allow for easier harvest of root material. Samples were collected following a modified harvest protocol designed for the study of perennial grass roots (McPherson et al., 2018). Three plants were randomly selected for harvest from each plot. Each plant, including roots, was collected by digging a soil core with a radius of ~20 cm around the sample plant stem. Loose soil was shaken off in the field. Plants were cut at the hypocotyl and plant tissue from above the hypocotyl was discarded.

Using pruning scissors sterilized in 70% EtOH, 4-6 basal roots were clipped from the main root structure. This included the largest basal root in all samples. Samples from each plot were harvested in succession and each plant constituted a separate replicate from a given plot. The replicates from each plot were considered separate samples, and were analyzed without pooling. After each plot was harvested, the workstation was sterilized using 70% EtOH. The root clippings were placed in 50 mL Falcon tubes containing 35 mL autoclaved phosphate buffer and a surfactant, Tween. Falcon tubes were placed on ice and returned to the lab for further processing.

In the lab, the Falcon tubes were vortexed for 2 minutes to release rhizosphere soil. Using sterilized forceps, the roots were removed from the tubes and blotted dry on paper towels. The remaining rhizosphere soil was vortexed for several seconds to resuspend, before pelletizing via centrifugation. A secondary wash with sterile phosphate buffer was done, before pipetting the suspended soil samples into clean 2 mL microfuge tubes. The samples were centrifuged into pellets again, and the supernatant buffer was removed using a pipette. Soil samples were immediately stored at -25°C, until further processing.

DNA was extracted from the rhizosphere soil samples using the E.Z.N.A.® Soil DNA Kit (Omega Bio-tek, Inc., Norcross, GA, USA), following the manufacturer’s instructions. DNA samples were quantified at 260/280 nm wavelength on the BioTek Synergy 2 Microplate reader to ensure DNA was within the normal range for absorbance and exceeded a quantity of 2ng/uL. Following quantification, the samples were stored at -20°C until further processing.

The samples were sent to the Center for Quantitative Life Sciences (CQLS) at Oregon State University for amplification and sequencing. We used the Earth Microbiome Project protocols and primers for amplification of both the 16s V4 and ITS regions and library prep (Caporaso et al., 2012; Walters et al., 2016; Ul-Hasan et al., 2019). PCR product size was determined by visualization with the TapeStation HS DNA tape and quantified using the protocol outlined in the Illumina Library qPCR Quantification Guide. PCR product quantity was double-checked using the Qubit fluorometer, and samples that failed to show quantities more than the blank used in the standard curve (0.3051ng uL^-1) were excluded from sequencing steps. Libraries were normalized to the lowest concentration >1ng uL^-1, then pooled. Any samples<1ng uL^-1 were added to the pool without normalization. The libraries were then sequenced at the CQLS on Illumina MiSeq (300bp paired-end). Ultimately this process produced two groups of sequence data, as described in Supplementary Table 3.

Raw reads were processed for each group by assessing Phred scores, trimming forward and paired reads, denoising data, and creating contigs. Amplicon sequence variants (ASVs) were then generated by removing chimeras and classifying sequences against a reference training database using the RDP Classifier algorithm within the DADA2 package for R (Callahan et al., 2016a). DADA2 was selected for use due to the high specificity and resolution compared to other methods (Prodan et al., 2020). ASVs were classified in January 2022 for both 16s and ITS sequences. We classified 16s sequences using the Silva rRNA database with a bootstrap confidence level of 60, and ITS sequences using the UNITE database with a bootstrap confidence level of 60 (Quast et al., 2012; Nilsson et al., 2019). Phylogenetic trees were created for each group of sample sequences using the phangorn 2.8.1 package, and were compiled with the other ASV and sample data using the phyloseq package (McMurdie & Holmes, 2013; Schliep et al., 2017). Sequences determined to be outliers were removed from the phyloseq object, and sequences that were sourced from parents and positive or negative controls were combined in a separate phyloseq object for use in downstream analyses. Chloroplast and mitochondrial DNA were removed from the sequence data before further analysis.

Prior to community composition analysis, the samples underwent taxonomic filtering using the phyloseq package and low prevalence reads were removed. These steps are intended to remove rare or spurious ASVs that may result from contamination or sequencing errors (Callahan et al., 2016b). Although this process may diminish the magnitude of some α-diversity comparisons, significant relationships are retained and the results are likely to be more reproducible and comparable across studies (Cao et al., 2021). To explore community diversity and richness within each sample population, we generated Shannon and Simpson metrics on non-rarefied data. This initial exploration of α-diversity was supplemented with a secondary analysis of similar metrics, classified as Richness and Inverse Simpson, using a bootstrapping method to normalize the data (Berry et al., 2017).

All analyses of β-diversity were completed using the vegan package, unless otherwise noted (Oksanen et al., 2012). β-diversity, which encompasses the differentiation between two communities or populations, was explored first via unconstrained and constrained ordinations. For this we utilized a principal correspondence analysis (PCoA) and a canonical analysis of principal coordinates (CAP) to generate the unconstrained and constrained ordinations, respectively, for both breeding history and parentage. ANOVAs were run on the ordinations of each sample type constrained under each set of explanatory variables. Unconstrained ordinations, though valuable in initial pattern assessment, are limited in power by data variability and structure, and while constrained analysis provides a clearer view of group-based variation, it does not capture comprehensive patterns in multivariate data (Anderson and Willis, 2003). Neither approach accounts for the influence of total counts on patterns within the data, thus a secondary analysis of β-diversity was run utilizing centered-log ratio (CLR) transformations. The CLR-transformed data was used to create a Euclidean distance matrix on which permutational analyses of variance (PERMANOVAs) were run to test the significance of various explanatory variables. The homogeneity of variance was tested for each of the explanatory variables using the betadisper function.

Following significance tests of explanatory variables, we sought to identify biologically relevant information on the specific taxa responsible for the variation between sample groups. To do this, we utilized the sparse partial least squares-discriminant analysis (sPLS-DA), with a maximum distance method, from the mixOmics package (Rohart et al., 2017). The sPLS-DA method is an extension of a PLS regression analysis that includes a LASSO penalization to select variables, which in the context of this study is ASV taxon groups (Lê Cao et al., 2011). It is further extended to become a discriminant analysis by encoding the response matrix with class variables, such as organic and conventional breeding history. The taxon groups most responsible for shifts in centroids can then be extracted and are interpreted to be the taxa most responsible for community composition differences. Following our discriminant analysis, the 10 ASV groups most responsible for shifts between the organic and conventional classes along each principal component were compiled. The raw nucleotide sequence for each ASV group was run through the NCBI BLAST taxonomy search, and additional information generated from the sequence regarding identity and putative function were collected (Schoch et al., 2020).

A series of ANOVAs were run to compute variance components for the estimations of narrow-sense heritability. Narrow-sense heritability was estimated for several microbiome metrics, including those representing both α- and β-diversity.

Narrow-sense heritability was estimated for each category of data, using the following formula:

where σ G 2 is the estimated genotypic variance component (M. Singh et al., 1993).

Heritability was calculated for each trait on a family means basis using the dplyr and EnvStats packages for R (Wickham, 2009; Millard, 2013). An analysis wherein all model components were treated as random effects, was used to obtain the variance components.

In April of 2020 the soils in the Lewis Brown Farm research plots had a pH of 6.4, along with several other nutrient and organic matter metrics seen in Supplementary Table 4. The first test took place one month before planting, and prior to fertilizer and compost application. The second test, in July 2020, indicated a drop in pH to 6.25 along with changes in several other nutrient metrics. Notably, phosphate decreased slightly while potassium increased. Microbial respiration was 35.9 µg CO₂-C/g dry soil/day. Comparing this to other estimates of microbial respiration, this is considered very low soil activity. Microbial respiration increases with temperature, therefore the value would have increased into August, although may not have exceeded moderately low activity at that point either.

The four populations of Phaseolus vulgaris hosted a diverse community of fungi and bacteria. After the removal of outliers, controls, and parents, the analysis identified 21 bacterial phyla comprised of 809 unique ASVs and 4 fungal phyla comprised of 129 unique ASVs. Bacterial communities were dominated by Proteobacteria (47% of community), followed by Bacteroidota (18.5% of community), Acidobacteriota (6.3% of community), and Actinobacteriota (6.2% of community). Among fungal communities Basidiomycota dominated (53.5% of community), followed by Ascomycota (34.9% of community) and Mortierellomycota (9.3% of community).

Our assessment of within-group microbiome diversity (α-diversity) for the two breeding history groups identified trends across both 16s and ITS communities. Broadly, the organically bred populations displayed slightly more evenness in both the micro- and myco- biomes, as seen in Figures 1A, C). This trend was also observed in the richness values determined for the 16s sequences, Figure 1B), although it was diminished in the ITS sequence analysis, Figure 1D). While the increase in evenness and richness was consistent across metrics and sample types, the increases were not statistically significant.

Initial and exploratory unconstrained ordinations were generated using a PCoA. These indicated a general shift in the bacterial microbiome associated with breeding history, as well as a smaller shift associated with cross, as seen in Supplementary Figures 1, 2. Similar shifts were observed in the mycobiome, as seen in Supplementary Figures 3, 4, wherein both breeding history and cross were associated with community shifts. Subsequent constrained analyses accounted for these shifts within the model of the ordination, and were paired with ANOVAs to test the significance of community shifts. For the bacterial microbiome, breeding history retained an observable shift in community composition (Figure 2A) and the ANOVA indicated the shift described by breeding history was significant at α< 0.05. When the model included both breeding history and cross, a heightened shift was observed in the ordinations (Figure 2B) and the paired ANOVA found the shift to be significant at < 0.01. Although cross amplified the breeding history shift, when analyzed independently it was not found to cause a shift in ordination of any significance. Exploring shifts in the mycobiome, the constrained ordination of the data modeled on breeding history showed a very slight shift that was not significant, as seen in Figure 3A. The ordination of the data modeled on cross had a markedly larger shift (Figure 3B), which was significant at <0.05.

To account for the skew of count data, the relationships identified via the unconstrained and constrained ordinations were further explored through a permutational analysis of variance performed on centered log transformed data. Within the 16s sequence analysis, breeding history was a significant explanatory variable, while cross was not (Supplementary Table 5). Both cross and breeding history had normal dispersion as well. For the ITS sequence analysis, both breeding history and cross were significant explanatory variables for the variation seen across the samples (Supplementary Table 6). The homogeneity of dispersion for these variables did not indicate any significant, non-normal dispersion.

Exploring these significant relationships further, a sPLS-DA model was tuned to fit the data described in (Table 1). Using the classification that was developed through this model tuning, it was possible to estimate the predicted ordination area for each class. For each of the sample types described in Table 1 a supplemental plot was created that includes the sample ordination along with the predicted background. In Figure 4 the sample ordination plots cleanly within the predicted background for the 16s sequences with no crossover across the background prediction divide. For the ITS sequences evaluated by breeding history classes, several samples from the organic class were plotted into the conventionally predicted background, as seen in Figure 5. The ITS sequences evaluated by parental cross, shown in Figure 6, had similar crossover events. In general, for all data sets, most of the samples fell into the ordination space that was predicted using the classification model. The number of crossover events across the predicted background divide in each ordination trends upwards in concert with the error rates for each of the ordination components.

From the ordination loadings it was possible to extract the ASV information that were most responsible for shifts between each of the two classes in each data set. The nucleotide sequence information associated with these loadings was compared to available GenBank sequences using NCBI BLAST. The top match identified in GenBank is reported in Tables 2–4 along with relevant taxonomic and functional information for each sequence.

Within the context of the 16s sequences (Table 2), there were notable increases in several potentially beneficial taxa in samples derived from beans with both an organic and conventional breeding history. Within conventionally bred accessions, plant-growth promoting bacteria Caulobacter hibisci and Arthrobacter pokkalii were increased along with the nitrogen-fixing bacteria Rhizobium mesoamericanum and Azospirillum canadense. The plant-growth promoter Arthrobacter pokkalii was increased in multiple instances in the organically bred accessions.

In the context of ITS sequences, both beneficial and pathogenic fungal taxa were present in samples from both breeding histories, as seen in Table 3. Among the conventionally bred samples there was an increase in the plant pathogen Rhizoctonia carotae, however there was also an increase in the potentially beneficial fungi Minimedusa polyspora and Laetisaria arvalis, the latter of which is known to have biocontrol properties over Rhizoctonia. In the organic bean samples, there were multiple fungi with increased abundance that are known to have beneficial properties in plant-association. Among these are: Laccaria moshuijun, a known ectomycorrhizal basidiomycete; Pleurocatena brevior, a mycoparasite; and Rhizophagus prolifer, an arbuscular mycorrhizal fungus. Additional species were identified in the organically bred samples that had less definitive roles in the rhizosphere mycobiome, such as Cadophora gamsii and Hormonema viticola. Increases in the plant pathogen Ceratobasidium ramicola were also observed.

Shifts in the mycobiome due to parentage or cross are presented in Table 4. Among accessions derived from the HYPR cross, the only fungal species that increased was Minimedusa polyspora, which is notable for its antibiotic properties. In accessions from the ORBV cross several plant pathogens were elevated: Ceratobasidium ramicola, a relative of Rhizoctonia; Phaeosphaeria gahniae and Helgardiomyces anguioides, both pathogens in cereal crops; and Clarireedia bennettii, a turfgrass pathogen. There were no notable increases observed among the ORBV accessions in taxa that have putatively beneficial function in the mycobiome.

In Table 5 the narrow-sense heritability estimates for microbiome community composition traits are summarized. β-diversity represents the amount of diversity between sample groups, such as between two breeding histories, and is represented by two principal components derived from an ordination of the trait. Heritability for β-diversity in the 16s data exceeds the estimates made in the ITS data, although this distinction diminishes in the second PC metric used as a trait. The two metrics evaluated for α-diversity, which considers the within sample group diversity, were the taxon richness and the evenness of taxon abundance. Regarding α-diversity, the 16s and ITS samples for each population were quite variable, ranging from low to moderately high heritability. Species richness appeared to have more stable heritability in the fungal environment.