Genotypic variation in floral volatiles influences floral microbiome more strongly than interactions with herbivores and mycorrhizae in strawberries

The floral microbiome is of significant relevance to plant reproduction and crop productivity. While plant genotype is key to floral microbiome assembly, whether and how genotypic variation in floral traits and plant-level mutualistic and antagonistic interactions at the rhizosphere and phyllosphere influence the microbiome in the anthosphere remain little known. Using a factorial field experiment that manipulated biotic interactions belowground (mycorrhizae treatments) and aboveground (herbivory treatments) in three strawberry genotypes, we assessed how genotypic variation in flower abundance and size and plant-level biotic interactions influence the bidirectional relationships between floral volatile organic compounds (VOCs) and the floral microbiome using structural equation modeling. We found that plant genotype played a stronger role, overall, in shaping the floral microbiome than biotic interactions with mycorrhizae and herbivores. Genotypic variation in flower abundance and size influenced the emission of floral VOCs, especially terpenes (e.g. αand β-pinene, ocimene isomers) and benzenoids (e.g. p-anisaldehyde, benzaldehyde), which in turn affected floral bacterial and fungal communities. While the effects of biotic interactions on floral traits including VOCs were weak, mycorrhizae treatments (mycorrhizae and herbivory + mycorrhizae) affected the fungal community composition in flowers. These findings improve our understanding of the mechanisms by which plant genotype influences floral microbiome assembly and provide the first evidence that biotic interactions in the rhizosphere and phyllosphere can influence the floral microbiome, and offer important insights into agricultural microbiomes.


Introduction
The f loral microbiome that consists of thousands or more taxa of bacteria and fungi mediates processes that are essential for plant reproduction and crop productivity [1][2][3][4][5]. While we are just beginning to understand the principles and functions of microbiome assembly in f lowers [2,3], studies have revealed that the floral microbiome is governed by plant genotype [1], and it overlaps with the leaf microbiome harboring a subset of the microbiota hosted by the root microbiome following a source-sink gradient [5]. Thus, plant genotypic variation in floral traits and species interactions that dictate the microbiomes in the rhizosphere (roots) and phyllosphere (leaves) are expected to inf luence the microbiome in the anthosphere (f lowers).
Floral traits can inf luence the f loral microbiome via distinct mechanisms. For instance, f lower size is hypothesized to affect the niche availability for microbiota [6]. Flower abundance, on the other hand, can influence floral microbiome assembly by governing the source pool of colonizing microbes [1,7] and by enhancing pollinatormediated microbial dispersal [1]. In particular, floral volatile organic compounds (VOCs) have been shown to mediate pollinator foraging [8], which can potentially influence microbial dispersal [3]. In addition to the pollinator-mediated effect of floral VOCs on the floral microbiome, VOCs are known for their antimicrobial properties (e.g. terpenes, benzenoids, aliphatics) [9][10][11][12] and their roles as carbon sources (e.g. methanol) [9], and thus can potentially reduce pathogenic microbes [10,13], modulate microbe-microbe interactions [14] and thereby affect microbial communities [12]. Microbiota, in turn, can influence plant VOCs by directly contributing to VOC emissions [3] or indirectly by metabolizing plant VOCs and/or affecting plant physiology that can influence plant VOC emissions [9,10,13]. Thus, the effect of microbiota on floral VOCs is determined, in part, by the quantity and quality of plant metabolites [9] in the anthosphere. Such bidirectional relationships between floral VOCs and floral microbiome are likely to be influenced by biotic interactions experienced by the whole plant. For instance, plant interactions with rhizosphere mutualists (e.g. mycorrhizae) [15] and phyllosphere antagonists (e.g. herbivores) [16,17] may affect the source-sink dynamics of microbes that colonize the anthosphere and/or affect floral traits including VOCs [13,15], which can inf luence floral microbiome assembly. Despite the connections between whole-plant biotic interactions, f loral VOCs and the floral microbiome, no study has evaluated the relative importance of genotypic variation in f loral traits and plant-level interactions in shaping the f loral microbiome, nor in affecting the bidirectional relationships between floral VOCs and f loral microbiome.
Here we leveraged a factorial field experiment [18] that manipulated biotic interactions belowground (mycorrhizae treatments) and aboveground (herbivory treatments) in three day-neutral strawberry cultivars (Fragaria ×ananassa "Seascape", "Tribute" and "Wasatch"; Fig. 1). As clonal plants of each cultivar were used in this study, we referred to cultivars as genotypes. We first examined how epiphytic bacterial and fungal communities in flowers responded to plant genotype, mycorrhizae treatments and herbivory treatments. We then assessed how plant genotypes varied in f loral traits (f lower abundance, size and VOCs) that may shape the f loral microbiome. Lastly, we used structural equation modeling to link genotypic variation in f lower abundance and size and biotic treatments to the bidirectional relationships between floral VOCs and f loral microbiome (bacterial and fungal α-and β-diversity).
Similar to α-diversity, the β-diversity of floral bacterial communities was influenced by plant genotype (Fig. 2c

SEMs linking floral phenotype and treatment to the floral microbiome
Given the significant genotypic differences in both floral phenotype and microbiome, we then used SEMs to link genotypic variation in flower abundance, size and VOCs, as well as treatments, to the α-and β-diversity of floral bacterial and fungal communities (Figs. 4a and 5a). The SEMs focused on VOCs at three different levels: (1) individual VOCs that differed between genotypes as identified by the random forest classification (see above), and that predicted the α-or β-diversity of bacterial and fungal communities based on backward predictor selection (see Materials and methods; Supplementary Tables 4 and 5); (2) the major classes of VOCs (terpenes, aliphatics, benzenoids, S-containing compounds and C5 branchedchain compounds); and (3) VOC profile (VOC cPCoA1 and cPCoA2; Fig. 3c). The SEMs were examined in parallel  Table 2). Significant contrasts of least-squares means are denoted: * P < 0.05; * * * P < 0.001. Error bars represent 1 SE. c, The constrained principal coordinates analysis (cPCoA) revealed significant differences in f loral volatile organic compounds (VOCs) among genotypes (detailed statistics in Supplementary Table 3).
The β-diversity of floral bacterial communities was linked with a different set of VOCs relative to the α-diversity (Fig. 4), as revealed by the SEMs based on overall genotypic variation. Specifically, p-anisaldehyde was positively linked with the first axis of bacterial community composition (cPCoA1, r = 0.24, P = 0.019; Fig. 4d), which was dominated by Lactobacillus micheneri (found in both flowers [1] and pollinators [1,26]) and Methylobacterium and Acinetobacter nectaris (nectar taxa [21]) on the positive axis and Delftia (common taxa in strawberry flowers [1]) on the negative axis (  Table 1). The cPCoA2 of bacterial communities was also influenced by flower size indirectly via the enhanced emission of VOCs (ocimene isomers, r = 0.28, P = 0.021, Fig. 4e; aliphatics, r = 0.29, P = 0.008, Fig. 4f), but was not influenced by flower abundance or treatments (Fig. 4). These relationships based on overall genotypic variation were also detected within genotypes ( Supplementary Fig. 1).
The β-diversity of fungal communities, on the other hand, were influenced by treatments and flower abundance ( Fig. 5c-f), as revealed by the SEMs based on overall genotypic variation. Relative to the control, mycorrhizae treatment positively affected the cPCoA1 of fungal communities (all P < 0.10; Fig. 5c, e, f and Supplementary Table 5), which was dominated by Pseudopithomyces chartarum (potentially pathogenic [25]) and Candelaria concolor (lichenized fungal taxon [27]) on the positive axis and pathogenic Alternaria alternata [28] and Cladosporium delicatulum [25] on the negative axis ( Fig. 2d; Supplementary Table 1). In addition, herbivory The f loral bacterial microbiome was characterized using α-diversity (Shannon diversity) and β-diversity (the first two axes of constrained principal coordinates analysis, cPCoA1 and cPCoA2). VOCs were evaluated at three levels (see Materials and methods): individual VOCs, major classes of VOCs, and VOC profile (VOC cPCoA1 and cPCoA2). The bidirectional links between VOCs and bacterial microbiome ref lected potential reciprocal inf luence on each other. Treatments were coded using the control treatment as the reference level. SEMs were fitted based on overall genotypic variation using lavaan [19]. b-g, SEMs that contained notable positive (solid black) and negative (dashed black) links between VOCs and bacterial microbiome, as well as other variables, were presented (see SEM details in Supplementary Table 4): † P < 0.10; * P < 0.05; * * P < 0.01; * * * P < 0.001. Grey links indicate P > 0.10. Link width and the numbers adjacent to links indicate standardized path coefficients. + mycorrhizae treatment positively affected the cPCoA2 of fungal communities (r = 0.48, P = 0.003; Fig. 5g), which was dominated by C. concolor and A. alternata on the positive axis and Rhodotorula graminis (nectar yeast [21]) and Naganishia randhawae (=Cryptococcus randhawai [29], yeast also found in strawberry f lowers [1]) on the negative axis ( Fig. 2d; Supplementary Table 1). In addition to treatments, f lower abundance affected the cPCoA1 of fungal communities indirectly via the reduced emission of VOCs per unit f lower dry weight (Fig. 5c-f), including benzenoids (r = −0.24, P = 0.014; Fig. 5e) such as benzaldehyde (r = −0.21, P = 0.024; Fig. 5c) and VOC profile (VOC cPCoA1, r = −0.26, P = 0.027; Fig. 5f) that also primarily ref lected the variation in benzaldehyde (Fig. 3c), as well as ocimene isomers (r = −0.24, P = 0.053; Fig. 5d). For the cPCoA2 of fungal communities, f lower size exhibited an indirect effect via the enhanced emission of terpenes (r = 0.39, P < 0.001; Fig. 5g). While most of these relationships were also detected within genotypes ( Supplementary Fig. 2), the links between fungal cPCoA1 and benzaldehyde (and related VOCs including benzenoids and VOC cPCoA1 that was dominated by benzaldehyde; Fig. 5c, e, f) as well as mycorrhizae treatment showed otherwise, indicating the role of among-genotype variation in driving some of the relationships.

Discussion
Our results demonstrated that strawberry genotype played a stronger role, overall, in shaping the floral microbiome than whole-plant biotic interactions both belowground (with mycorrhizae) and aboveground (with herbivores). In particular, we found that genotypic variation in flower abundance and size influenced floral VOCs, especially terpenes and benzenoids, which in turn The f loral fungal microbiome was characterized using α-diversity (Shannon diversity) and β-diversity (the first two axes of constrained principal coordinates analysis, cPCoA1 and cPCoA2). VOCs were evaluated at three levels (see Materials and methods): individual VOCs, major classes of VOCs, and VOC profile (VOC cPCoA1 and cPCoA2). The bidirectional links between VOCs and fungal microbiome ref lected potential reciprocal inf luence on each other. Treatments were coded using the control treatment as the reference level. SEMs were fitted based on overall genotypic variation using lavaan [19]. b-g, SEMs that contained notable positive (solid black) and negative (dashed black) links between VOCs and fungal microbiome, as well as other variables, were presented (see SEM details in Supplementary Table 5): † P < 0.10; * P < 0.05; * * P < 0.01; * * * P < 0.001. Grey links indicate P > 0.10. Link width and the numbers adjacent to links indicate standardized path coefficients. affected the f loral bacterial and fungal communities. While the effects of biotic interactions on f loral phenotype including VOCs were weak (Figs. 3-5), mycorrhizae treatments (mycorrhizae and herbivory + mycorrhizae) directly affected the fungal communities in flowers relative to the control (Fig. 5). By linking genotypic variation in f loral traits and biotic interactions to floral bacterial and fungal communities, our study revealed the mechanisms by which plant genotype inf luences floral microbiome assembly [1] and also provided the first evidence that plant-level biotic interactions at the rhizosphere and phyllosphere can inf luence the anthosphere microbiome.
In line with previous studies on the f loral microbiome [1], genotypic variation in f lower abundance inf luenced the f loral microbiome. Flower abundance has been shown to govern f loral microbiome assembly both directly by inf luencing the source pool of microbes that colonize individual f lowers [1,7] and indirectly by influencing pollinator visits [30] that can mediate microbial dispersal [1]. Our results showed that flower abundance affected fungal α-and β-diversity via the bidirectional links between VOCs and the fungal communities among genotypes but not within genotypes. In addition to flower abundance, flower size that is hypothesized to influence the floral microbiome has, nevertheless, rarely been tested [6]. In this study, we found that flower size affected bacterial and fungal communities via the bidirectional links between VOCs and the floral microbiome, which was supported by SEMs based on overall genotypic variation and withingenotype variation as well. While the causation of these relationships awaits experimental validation, it provides a new insight into the mechanisms by which flower abundance and size drive microbiome assembly.
Our results provided evidence for the hypothesized bidirectional relationships [9] between plant VOCs and microbiome in the anthosphere. Specifically, monoterpenes, namely α-and β-pinene and ocimene isomers contributed most to the variation in floral VOC profile and f loral bacterial communities. Bacterial Shannon diversity was negatively associated with α-and β-pinene, whereas bacterial community composition was associated with p-anisaldehyde, ocimene isomers, aliphatics and hexenyl tiglate (the only C5 branchedchain compound in this study). These bidirectional relationships were supported by overall genotypic variation and within-genotype variation. On the other hand, fungal Shannon diversity was positively linked with benzothiazole, the only S-containing compound detected in strawberry volatile collections here, whereas fungal community composition was associated with terpenes (especially ocimene isomers) and benzenoids (especially benzaldehyde). Yet, the bidirectional relationship between benzaldehyde and fungal community composition was only supported among genotypes. To further verify that the relationship is not caused by unmeasured genotypic differences, experimental investigation that directly tests the reciprocal effects is needed. Overall, despite the known antimicrobial effects of f loral VOCs such as terpenes, benzenoids, aliphatics and S-containing compounds [9][10][11][12], our results detected microbial taxa both negatively and those positively associated with the emission of these VOCs. This finding suggests the potential dual roles of floral VOCs as antimicrobial agents and carbon sources for some microbial taxa [9,10] and/or their role in mediating complex microbe-microbe interactions [14] where the suppression of some taxa can potentially increase other taxa via direct and indirect interactions through the microbial network [1]. While how dynamic the relationships are between f loral VOCs and the floral microbiome remains an open question, our results are expected to hold despite the time interval between VOC and microbiome collections, due to stable VOC composition (i.e. relative abundances of individual VOCs) among these genotypes detected over a longer time period (4 weeks) in this field experiment [18]. Taken together, our findings highlight the needs for future research on unraveling the specific mechanisms that underlie the relationships between f loral VOCs and the floral microbiome as well as their temporal dynamics.
While plant interactions with herbivores at the phyllosphere have been found to inf luence the leaf microbiome [16,17], our results showed, for the first time, that plant interactions at the rhizosphere and phyllosphere can inf luence the f loral microbiome. Such effects were detected on shaping the floral fungal communities. Specifically, the mycorrhizae treatment caused the deviation of fungal community composition from the control (Fig. 5) by reducing pathogenic taxa (e.g. A. alternata, Cladosporium delicatulum) in f lowers, whereas the herbivory + mycorrhizae treatment inf luenced fungal community composition by reducing f loral yeasts (e.g. R. graminis, Naganishia randhawae) and increasing pathogenic taxa (e.g. A. alternata). While the effect of mycorrhizae treatment occurred among genotypes, the effect of herbivory + mycorrhizae treatment was supported by SEMs based on overall genotypic variation and within-genotype variation as well. These effects of biotic interactions on the floral microbiome were, nevertheless, not via floral VOCs, although plant mutualists and antagonists have been shown previously to alter floral VOCs [15,31,32]. In this study, the mycorrhizae treatments and herbivory treatments exhibited weak effects on the VOC profile. Yet, we could not rule out the possibility that the effects of the biotic treatments on floral VOCs can be dynamic and thus attenuate over time, and plant genotypes and their floral VOCs might respond differently to natural mycorrhizal inocula compared to a single mycorrhizal inoculum used in this study. Whether the effects of mycorrhizae and herbivory + mycorrhizae treatments on floral fungal microbiome are caused by altered source pools of microbes from the rhizosphere and phyllosphere that can affect the source-sink dynamics or altered plant traits not measured here that can affect floral fungal taxa merits further investigation.
To conclude, plant genotype governs floral microbiome assembly via both direct and indirect mechanisms mediated by diverse floral traits, especially the bidirectional relationships between floral VOCs and bacterial and fungal communities. In addition to the intrinsic characteristics of plant genotypes, interactions with nonfocal mutualists and antagonists at the rhizosphere and phyllosphere influence the anthosphere microbiome that can profoundly affect plant reproduction [1,2,4]. These findings should not be unique to strawberries but may be generalizable to other crops and wild plants, because the VOCs of important relevance to the floral microbiome are common among plant species [9][10][11][12][13], and mycorrhizal partners and herbivore pests are ubiquitous in natural and agricultural ecosystems [33,34]. We emphasize that our findings are based on plants with pre-existing endophytic microbiomes; as such, the extent to which endophytic microbiomes influence the detected relationships remain to be determined. Nevertheless, our findings based on strawberries improve the understanding of microbiome assembly in flowers and offer important insights into agricultural microbiomes that can influence plant reproduction and crop productivity.

Factorial field experiment
The factorial field experiment [18] was conducted during March-September, 2019, at the Oakland University Student Farm, Rochester, Michigan (42.66073 • N, 83.19558 • W). Three strawberry genotypes ("Seascape", "Tribute" and "Wasatch") with 48 bare-root nursery stock plants per genotype were included in the experiment. The fresh biomass was recorded for each bare-root plant (i.e. initial plant weight) to account for potential influence of variation in initial plant size on the experiment. These bare-root plants (N = 144 total) were grown in 1-gallon pots filled with 2:1:1 mixture of sphagnum peat moss (Lambert, Rivière-Ouelle, Québec, Canada), washed play sand (Kolorscape, Atlanta, GA) and Turface MVP (Turface, Buffalo Grove, IL) during March 10-16, 2019, and re-potted into 2-gallon white fabric pots (247Garden, Montebello, CA) filled with the same media during May 24-29, 2019. These plants were subjected to four treatments (control, herbivory, mycorrhizae, and herbivory + mycorrhizae). Within each treatment, a single plant from each genotype was randomly selected (Fig. 1). The four treatments were then grouped into and replicated across 12 blocks (Fig. 1). For mycorrhizae-present treatments (mycorrhizae and herbivory + mycorrhizae), plants were inoculated with Rhizophagus irregularis spores (Elite 91 Myco Jordan, Murietta, CA) upon initial potting by dusting roots with inocula and also filling the 1-gallon pots half full of media and adding a thin layer of 5 mL inocula at the point of contact with roots. For herbivory-present treatments (herbivory and herbivory + mycorrhizae), plants experienced both natural herbivory from Vanessa cardui caterpillars (Carolina Biological Supply, Charlotte, NC; with three to five 2nd and 3rd instar larvae feeding on a single leaf per plant for six days) and simulated herbivory (by clipping the distal half each leaf from half of the leaves per plant and then applying 1 mM jasmonic acid solution in water [35] to cover all upper leaf surface) during June 17-July 6, 2019. The combination of natural and simulated herbivory is commonly used [36] to ensure that plants received the stimulus of a live herbivore and associated cues (e.g. enzymes in caterpillar saliva) while standardizing for the leaf area damaged per unit time.

Floral traits and VOC collection
In this study, we recorded the number of open flowers produced over a two-week period (August 4-17, 2019) prior to and during the f loral microbiome collection. We sampled f loral volatile emissions eight weeks after herbivory treatments, and 5-15 days prior to f loral microbiome collection (July 30-August 10, 2019) to minimize the potential inf luence of the activities of VOC measurements that can cause artificial microbial arrival and destructive f lower collection on f loral microbiome. Floral volatiles were collected from new fully expanded flowers between 11 AM and 3 PM on sunny to partly sunny days, with no to little chance of rain. During the 4-h sampling period each day, 16 f loral and two blank samples were collected using a dynamic headspace sampling method [31], by enclosing f lowers in 12 oz. polyethylene terephthalate cups with dome lids (Comfy Package, Rikkel Corp, NY) and pulling air through a semi-open system at a f low rate of 200 mL/min through volatile traps with 30 mg of HayeSep Q adsorbent (VAS, Rensselaer, NY). We used either a PVAS22 portable volatile assay system (VAS, Rensselaer, NY) or an Air Lite lowflow air sample pump (SKC, Eighty Four, PA), which were shown functionally equivalent in volatile collection during initial trials. Volatiles were analyzed on an Agilent 7890A gas chromatograph (GC) with an Agilent 5977B mass spectrometer (MS) and separated on a HP-5 column (30 m × 250 μm × 0.25 μm) with helium as a carrier gas. We measured flower diameter and dry weight for all flowers sampled. Floral VOC emissions (per hour) were standardized by flower dry weight.

Floral microbiome collection and sequencing
The floral microbiome was collected from individual strawberry plants that flowered on August 16, 2019 (N = 93 out of 144 plants). A single flower (without pedicel) per plant was collected into a sterile 15 mL centrifuge tube using ethanol-rinsed forceps. These flower samples were transported on dry ice to the University of Pittsburgh and stored at −20 • C within 12 h. Following the previous protocol [1], epiphytic microbes were collected by sonicating flowers in 3 mL phosphate-buffered saline at 40 kHz for 10 min and vortexing for 3 min. The microbes were pelleted by centrifuging at 10000 x g for 5 min, and then subject to DNA extraction using Quick-DNA Fecal/Soil Microbe Kits (Zymo Research, Irvine, CA). Samples and one negative control (without flower but in the process of microbial isolation and DNA extraction) were sent to Argonne National Laboratory for bacterial (16S rRNA V5-V6 region, 799f-1115r primer pair) and fungal (ITS1f-ITS2) library preparation [1]. Because the negative control failed in library preparation, the 93 samples were sequenced on a 1/2 lane of Illumina MiSeq (paired-end 250 bp).
Bacterial and fungal ASV tables were further filtered separately before conversion into microbial community matrices using package phyloseq [40]. First, we removed non-focal ASVs (Archaea, chloroplasts and mitochondria). Second, we filtered out samples of low reads (<200 reads). Third, we normalized per-sample reads to the same number (i.e. the median reads, 4114 and 10 377, bacterial and fungal data set, respectively) [1,5]. Lastly, we removed low-frequency ASVs (<0.001% of total observations). The final community matrices consisted of 1500 and 976 ASVs for bacteria (N = 86 samples) and fungi (N = 87), respectively.

Floral microbial α-and β-diversity
To evaluate how f loral microbial α-diversity was influenced by plant genotype, herbivory treatments (present: herbivory, herbivory + mycorrhizae; absent: control, mycorrhizae), and mycorrhizae treatments (present: mycorrhizae, herbivory + mycorrhizae; absent: control, herbivory), we analyzed the bacterial and fungal communities separately using general linear mixed models (LMMs) in package lme4 [41]. The response variable (Shannon diversity) was calculated using package vegan [42], and power transformed to improve normality with the optimal power parameter determined using the Box-Cox method in package car [43] (i.e. power parameter = 1 and 2 for bacterial and fungal Shannon diversity, respectively). The predictors included genotype, herbivory treatments and mycorrhizae treatments as well as their two-way interactions and three-way interaction, with initial plant weight as a covariate. The split plot with randomized complete block design in this experiment (Fig. 1) required nested random effects (i.e. treatment plots nested within blocks). However, due to only a subset of the 144 samples (N = 86 and 87 for bacterial and fungal data sets, respectively) for the LMMs, nested random effects led to difficulties in model convergence and thus were not used. Instead, the LMMs considered the random effect of blocks alone. Statistical significance (type III sums of squares) and least-squares means (LS-means) of predictors were assessed using packages lmerTest [44] and emmeans [45].
Microbial β-diversity (Bray-Curtis dissimilarity) was evaluated using PERMANOVA and constrained principal coordinates analysis (cPCoA) in vegan for bacterial and fungal communities separately. To assess the significance of the main effects, PERMANOVA and cPCoA included genotype, herbivory treatments and mycorrhizae treatments, with initial plant weight as a covariate and block as the random effect. To assess the significance of two-way interactions, PERMANOVA and cPCoA included both the main effects and their two-way interaction terms. Similarly, to assess the significance of the three-way interaction, the main effects, two-way interactions and three-way interaction were included. Once significant genotypic effect was identified, we further conducted the same PERMANOVAs and cPCoAs between pairwise genotypes to detect which specific genotypes were divergent from each other in microbial community composition.

Floral phenotype analysis
To evaluate whether strawberry genotypes differed in flower abundance (sample size N = 93) and flower size (N = 80), we conducted LMMs with the response variables power transformed to improve normality (i.e. power parameter = 0.5, square-root transformation, for flower abundance, and 2 for flower diameter). The predictors included genotype, herbivory treatments and mycorrhizae treatments as well as their two-way interactions and three-way interaction, with initial plant weight as a covariate and block as the random effect.
To evaluate whether and how strawberry genotypes differed in VOC profile (sample size N = 66), we conducted PERMANOVA and cPCoA. The predictors included genotype, herbivory treatments and mycorrhizae treatments as well as their two-way interactions and three-way interaction, with covariates (i.e. initial plant weight, pump ID and average daily temperature during VOC collection) and random effect (block). The PERMANOVA and cPCoA were conducted as in the above-mentioned microbial analysis to assess the significance of the main effects and interaction terms. Once significant genotypic effect was detected, we further identified the specific VOCs that differed between genotypes using random forest (RF) classification in packages caret [46] and randomForest [47]. The RF classification models were run for the full data with 1000 trees, and the number of randomly selected variables (i.e. individual VOCs) at each split of a decision tree was optimized using 10-fold cross validation in caret. RF model performance was assessed using outof-bag (OOB) error. The set of important, non-redundant VOCs (in the presence of potentially correlated VOCs) were selected using backward variable elimination with package varSelRF [48].

SEMs linking floral phenotype to microbial α-and β-diversity
SEMs were conducted to link genotypic variation in flower abundance, size and VOCs, as well as treatments, to the α-diversity and β-diversity of floral bacterial and fungal communities (Figs. 4a and 5a). We focused on VOCs at three different levels: (1) individual VOCs that differed between genotypes as identified by the random forest classification (see above), and that predicted the α-diversity (Shannon diversity) or β-diversity (the first and second axis of cPCoA, cPCoA1 and cPCoA2) of bacterial and fungal communities based on backward predictor selection using packages caret and leaps [49]; (2) the major classes of VOCs (terpenes, aliphatics, benzenoids, S-containing compounds and C5 branched-chain compounds); and (3) VOC profile (VOC cPCoA1 and cPCoA2; Fig. 3c). For the SEMs, flower abundance and flower diameter were power transformed (i.e. power parameter = 0.5 and 2, respectively). For the categorical variables of treatments, we coded the control treatment as the reference level so that the effects of the other three treatments (herbivory, mycorrhizae, herbivory + mycorrhizae) were relative to the control. All the variables of VOCs were log transformed (i.e. power parameter = 0, log(x + 1)). For floral bacterial microbiome (Fig. 4a), a total of 41 SEMs were conducted (14 SEMs for Shannon diversity, 14 for cPCoA1 and 13 for cPCoA2, Supplementary Table 4).
Similarly, for f loral fungal microbiome (Fig. 5a), 41 SEMs were conducted (14 for α-diversity, 13 for cPCoA1 and 14 for cPCoA2, Supplementary Table 5). To reduce the number of models and model complexity, SEMs were only retained when notable paths (P < 0.10) were identified, and then re-fitted with only the notable paths present to estimate standardized path coefficients. The SEMs were fitted using package lavaan [19] based on overall genotypic (within-and among-genotype) variation. Model estimation used robust maximum likelihood with Satorra-Bentler scaled χ 2 that can accommodate nonnormality in lavaan, and model fit was confirmed for each retained SEM (comparative fit index, CFI > 0.9; root mean squared error of approximation, RMSEA, the lower bound of 90% confidence interval < 0.05; standardized root mean squared residual, SRMR <0.1) [50]. With these retained SEMs (Supplementary Tables 4 and 5), to further confirm whether the detected relationships based on overall genotypic variation were also supported within genotypes, we re-fitted the retained SEMs (with only the notable paths present) based on withingenotype variation using package piecewiseSEM [20], where genotype was included as the random effect using package nlme [51]. The SEMs that contained notable links (P < 0.10) between VOCs and bacterial and fungal α-or β-diversity based on overall genotypic variation were presented in Figs. 4 and 5 (six for bacterial microbiome and six for fungal microbiome), and the same set of SEMs fitted based on within-genotype variation were presented in Supplementary Figs. 1 and 2.