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