Sequence Analysis
Demultiplexed sequences were trimmed to remove adapter and primer binding sites using Cutadapt (Martin 2011). Amplicon sequence variants (ASVs) were inferred using DADA2 (Callahan et al. 2016). Errant ASVs due to chimerization were detected using the “consensus” method in DADA2 and discarded. Any ASV with a sequence length of greater than 256 bp or less than 250 bp were discarded. Taxonomic classifications were assigned to each ASV using DADA2’s assignTaxonomy () function using the Silva reference database (version 132, Quast et al. 2013).
Microbiome data was analyzed in R software (R Core Team 2020). ASVs assigned to mitochondrial and chloroplast lineages were discarded prior to normalization. For principal coordinates analysis (PCoA) and phylum level abundance statistics, raw counts were normalized to account for differences in sequencing depth between samples by dividing each ASV count by sequencing depth of a particular sample and multiplying by 1000 to place the counts on a per mille scale. Principal coordinate analyses were conducted using the capscale () function in the package Vegan (Oksanen et al. 2020). Bray Curtis dissimilarity on log2 transformed abundances was used for all PCoAs unless otherwise noted. Alpha diversity was calculated using Shannon Entropy from the diversity () function in Vegan. Differential abundance of aggregated phylum abundances was performed using linear models on log2 transformed abundances. Differential abundance of ASVs between conditions was conducted using DESeq2 on raw counts (Love et al. 2014).
Plant trait data from all parental replicates was analyzed to test genotypic and microbial treatment effects on plant morphological traits. We fit factorial linear mixed models using PROC MIXED in SAS (Littell et al. 1996) consisting of Ecotype, Treatment, and Ecotype x Treatment interactions as fixed effects. Preliminary analysis did not show any significant differences (in all cases, P > 0.113 between MAI and MCI treatments for parents and RILs), thus the average between them was used for this and all subsequent analyses (hereafter referred to as the Mock Inoculated (MI) treatment).
To explore the impact of the microbiome on the quantitative genetic architecture of our measured traits, we fit linear mixed models testing for GxE using the sommer package (Covarrubias-Pazaran 2018) in R based on the additive and epistatic relationship matrix determined from the genotypic data of the RIL (full approach described in Appendix S1). We calculated broad-senseH 2 as (Va+ Vaa)/Vp and present variance components and model comparisons.
The observation of different QTL effects under different treatment conditions provides evidence for QTL x environment interactions. To detect QTL present in the AI, CI and MI treatments, we completed QTL mapping on RIL values in R using the R/qtl package (Broman & Sen 2009) in each environment separately (full description in Appendix S1). We further tested for QTL x environment interactions in a full linear model incorporating data from the three treatments using the PROC MIXED procedure of SAS. Marker x treatment interaction indicates QTL x environment interaction, marker x marker interaction represents epistasis averaged over the environments, and marker x marker x treatment interaction indicates environment specific epistasis. We performed this analysis to test GxE interaction effects by contrasting the AI, CI and MI microbiomes (e.g. potentially identifying different soil or residual microbiome impacts). To test the significance of individual marker alleles at each treatment, we used the slice function in SAS as tests of simple effects (Winer 1971) for all significant marker x treatment and marker x marker x treatment interactions (full approach described in Appendix S1).
Results