A set of 96 whole-genome sequenced accessions of A. thaliana collected in the TOU-A local population (Toulon-sur-Arroux, Burgundy, France, 46°38′53.80″N–4° 7′22.65″E) was used for the purpose of this study. As previously described in Frachon et al. (2017), the TOU-A population is highly polymorphic at both the phenotypic and genomic levels. Importantly, the very short linkage disequilibrium (LD) observed in this population (r0.52∼18bp) allows the fine mapping of genomic regions associated with natural variation of phenotypic traits down to the gene level (Brachi et al., 2013; Huard-Chauveau et al., 2013; Frachon et al., 2017; Aoun et al., 2020; Lonjon et al., 2020). These 96 accessions were chosen among a set of 195 accessions to be representative of the genetic diversity observed in the TOU-A population (Platt et al., 2010).

The maternal effects of the 96 accessions were reduced by growing one plant of each family for one generation under controlled greenhouse conditions (16-h photoperiod and 20°C) in early 2011 at the University of Lille, Lille, France. Given an estimated selfing rate of ∼94% in this population (Platt et al., 2010), the 96 accessions were considered as mostly homozygous along the genome.

In this study, we used three neighboring species commonly associated with A. thaliana in natural plant communities in France and detected in the TOU-A plant community (F. Roux, personal observation). These species are the meadow grass P. annua (Poaceae) with a low spreading growth form, the chickweed S. media (Caryophyllaceae), and the speedwell V. arvensis (Scrophulariaceae) both with a crawling growth form. The seeds for these three species were obtained from the Herbiseeds company.¹

The effects of population structure on phenotype-genotype associations have been demonstrated to be limited in the TOU-A population (Brachi et al., 2013; Baron et al., 2015). Nevertheless, GWA mapping was run using a mixed-model approach implemented in the software genome-wide efficient mixed model association (GEMMA, Zhou and Stephens, 2012). This model includes a genetic kinship matrix as a covariate to control for the effect of the demographic history of the TOU-A population. In this study, we discarded SNPs with more than 17 missing values across the 91 accessions. In addition, because rare alleles may lead to an inflation of low p-values (Atwell et al., 2010; Brachi et al., 2010), we only considered the SNPs with a minor allele relative frequency (MARF) > 10%, leaving us with 630,234 SNPs. GEMMA was run independently on H1F, DIAM, HD, and BIOMASS for each of the 12 treatments. H1F, DIAM, HD, and BIOMASS were genetically correlated with FLO (mean value ± SEs across the 12 treatments, BIOMASS: Pearson’s r = 0.768 ± 0.043, DIAM: Pearson’s r = 0.360 ± 0.133, H1F: Pearson’s r = − 0.434 ± 0.085, HD: Pearson’s r = − 0.625 ± 0.015). The FLO genotypic values were thus included as a covariate in the GEMMA analyses, thereby allowing controlling for the putative effects of a confounding developmental stage (Uffelmann et al., 2021).

Thereafter, we implemented a local score approach on the set of p-values provided by GEMMA. The local score allows accumulating the statistical signals from contiguous genetic markers to detect significant genomic regions associated with phenotypic natural variation (Fariello et al., 2017). In each QTL region, the association signal, through the p-values, will cumulate locally due to linkage disequilibrium between the SNPs, which will then increase the local score (Bonhomme et al., 2019). The tuning parameter ξ was fixed at 2 (Bonhomme et al., 2019). Significant SNP-phenotype associations were identified by estimating a chromosome-wide significance threshold for each chromosome (Bonhomme et al., 2019).

To identify the candidate genes associated with natural variation of plant-plant interactions, we retrieved genes around the significant zones identified by the local score approach. Using a custom script, we selected all genes inside the zones as well as the first gene upstream and the first gene downstream of these zones. The TAIR 10 database² was used as our reference, leaving us with 4,524 unique candidate genes (N_BIOMASS = 1,538; N_DIAM = 1,008; N_H1F = 1,559; and N_HD = 1,572; Supplementary Figure 2).

In order to identify significantly over-represented biological processes (P < 0.05) in response to the different neighborhoods for the different traits, we submitted the candidate gene lists to the classification superviewer tool (³ Provart and Zhu, 2003) using the MAPMAN classification.

Highly significant genetic variation was found across the 12 treatments (Figure 1) for the five phenotypic traits scored on the focal A. thaliana plants (Table 1 and Supplementary Table 1). Similar results were observed without considering the treatments where A. thaliana was grown alone or with clones (Supplementary Tables 2, 3). The narrow-sense heritability SNP-based estimates were significant for 96.7% of the 60 “phenotypic trait × treatment” combinations (Supplementary Table 4) and ranged from 0.36 to 0.99 (mean = 0.72 and median = 0.76), suggesting that a large fraction of the phenotypic variation observed within each treatment was driven by genetic differences among the local A. thaliana accessions (Figure 1).

Importantly, as evidenced by highly significant “treatment × accession” interactions, a strong genetic variation of reaction norms was found for the five phenotypic traits, with or without considering the treatments where A. thaliana was grown alone or with clones (Table 1 and Supplementary Tables 2,3). For the phenotypic traits DIAM, H1F, HD, and BIOMASS, the across-environment genetic correlations ranged from 0.06 to 0.97 (mean = 0.69 and median = 0.73; Figure 2), indicating that the ranking of accessions for a given phenotypic trait largely differed among the 12 treatments. This pattern of crossing reaction norms is well illustrated for the height from the soil to the first flower for which the mean value was similar among the 12 treatments (Figure 1, Table 1, and Supplementary Table 1). In contrast, the across-environment genetic correlations for FLO were very close to unity (min = 0.91, max = 0.97, mean = 0.95, and median = 0.95; Figure 2), suggesting that the ranking of flowering time was similar among the 91 accessions across the 12 treatments. FLO was therefore discarded from further analyses.

For each treatment, the genetic correlations among the four remaining phenotypic traits were significantly different from the unity (absolute values of Pearson’s r: min = 0.02, max = 0.89, mean = 0.53, and median = 0.62), suggesting that the traits scored in this study partly behave independently (Supplementary Figure 3).

To compare the genetic architecture underlying the competitive response of A. thaliana between the bispecific and plurispecific neighborhoods, we combined GWA mapping based on 630,234 SNPs (i.e., 1 SNP every 189 bp) with a local score approach (Bonhomme et al., 2019). This combination of approaches was demonstrated to be efficient in the TOU-A population, allowing a fine description of the genetic architecture (in particular, the detection of QTLs with small effects) associated with quantitative disease resistance to the bacterial pathogen Ralstonia solanacearum and the cloning of four of the detected QTLs (Aoun et al., 2020; Demirjian et al., 2021).

The genetic architecture associated with natural variation of the traits DIAM, H1F, HD, and BIOMASS was polygenic with a number of detected QTLs ranging from 13 for BIOMASS in the SVV treatment to 72 for H1F in the A treatment (mean = 32 and median = 28) (Supplementary Table 5). Across the 12 treatments, the mean number of detected QTLs was not significantly different among DIAM, H1F, HD, and BIOMASS (Supplementary Table 5).

For each phenotypic trait, the genetic architecture of the competitive response of A. thaliana was highly flexible among the 12 treatments. First, the genetic architecture of the competitive response of A. thaliana in the context of bispecific interactions was highly dependent on the identity of the neighboring species (Figure 3). For example, no significant QTL region was shared between the treatments PPP and SSS for H1F or between the treatments PPP and SSS for HD (Figure 3). Second, the genetic architecture of the competitive response of A. thaliana was highly dependent on the number of species co-occurring with A. thaliana. For instance, two very neat peaks of association were identified for HD in the middle of chromosome 3 in the treatment PSV (Figure 4A). While the first peak was not identified in the treatments PPP, SSS, and VVV, the second peak was identified in the treatment PPP but with a much lower association score (PSV: Lindley process value = 214 and PPP: Lindley process value = 9.5). Third, the genetic architecture of the competitive response of A. thaliana to the presence of a specific pair of neighboring species largely differed not only from the genetic architecture identified in the corresponding bispecific interaction treatments but also between the two related assemblages. As an illustration, we detected a neat association peak for H1F at the beginning of chromosome 1 in the treatment PSS (i.e., one P. annua individual + two S. media individuals), but neither in the treatment PPS (i.e., two P. annua individuals + one S. media individual) nor in the treatments PPP and SSS (Figure 4B).

From these QTLs, we retrieved 4,524 unique genes located within and in the vicinity of 1 kb of the significant QTL regions detected for each “phenotypic trait × treatment” (Supplementary Dataset 2). In agreement with the values of genetic correlations between these four traits (Supplementary Figure 3), a large fraction of candidate genes was specific to a given treatment (BIOMASS: 77.4%, DIAM: 63.4%, H1F: 77.2%, and HD: 70.6%), suggesting that the genetic bases are largely distinct among the 12 treatments (Figure 5A and Supplementary Figure 4). Therefore, in this study, the genetic architecture underlying intraspecific variability of the competitive response largely depends on both the composition and assemblage of the neighborhood of A. thaliana.

Across the 48 “phenotypic trait × treatment” combinations, we identified 139 significantly over-represented gene ontology (GO, MapMan functional annotation) terms, which belong to 25 unique GO processes (Supplementary Dataset 3). Interestingly, the RNA related processes correspond to the most abundant GO term over-represented found in 28 “phenotypic trait × treatment” combinations, with eight treatments for BIOMASS, seven treatments for DIAM and HD, and six treatments for H1F (Figure 5, Supplementary Figure 5, and Supplementary Dataset 3). The second most represented GO terms are transport and minor CHO metabolism processes with 13 and 12 “phenotypic trait × treatment” combinations, respectively (Figure 5, Supplementary Figure 5, and Supplementary Dataset 3), but with different distribution across the traits. On one hand, the transport processes were represented in the four traits, while minor CHO metabolism processes were not found for BIOMASS in any of the 12 treatments. Conversely, the C1-metabolism processes were mostly enriched for BIOMASS in six treatments out of the seven “phenotypic trait × treatment” combinations identified for these processes (Supplementary Figure 5). The remaining enriched MapMan processes are related to tetrapyrrole synthesis, miscellaneous, protein, redox, secondary metabolism, stress, lipid metabolism, hormone metabolism, cell, photosystem, signaling, and not assigned genes (Figure 5, Supplementary Figure 5, and Supplementary Dataset 3).

In agreement with the strong crossing reaction norms observed among the three bispecific interaction treatments, a GWA mapping approach revealed that the genetic architecture of A. thaliana competitive response in the bispecific neighborhoods was highly dependent on the identity of the neighboring species. Similar results were observed (i) in A. thaliana challenged with four neighboring species (including the three species used in this study) in field conditions (Baron et al., 2015) and (ii) in Oryza sativa challenged with the three weed species Echinochloa oryzicola, Monochoria korsakowii, and Schoenoplectus juncoides in greenhouse conditions (Onishi et al., 2018). Altogether, these results suggest that the effect of the identity of the neighboring species on the genetic architecture of competitive response is conserved among the focal plant species and across diverse phenotyping environments. The genetic architecture of A. thaliana competitive response was also highly flexible between the bispecific and plurispecific neighborhoods, suggesting that the genetic response of a particular accession of A. thaliana to a plurispecific neighborhood can be hardly predicted from the observation of its genetic responses in the corresponding bispecific neighborhoods. These results confirm our previous conclusions obtained from a GEA approach performed on 145 whole-genome sequenced natural populations of A. thaliana located in the southwest of France and described for various plant community descriptors (Frachon et al., 2019). Our results are also in line with previous studies reporting in A. thaliana (i) the limited overlap between QTLs underlying resistance to combined stress responses, such as drought plus the specialist herbivore Pieris rapae, and QTLs underlying resistance to the corresponding single stress response (Davila-Olivas et al., 2017; Thoen et al., 2017) and (ii) the distinct transcriptomic responses for single and combined stresses in diverse plant species, as illustrated for pathogen attack combined with heat stress (Suzuki et al., 2014; Desaint et al., 2021).

In this study, three non-exclusive and interconnected hypotheses can be proposed to explain the emergence of new QTLs in plurispecific neighborhoods. They are based on (i) the putative generation of new signals or elimination/modification of the pre-existing signals (e.g., light, aerial volatile organic compounds, root exudates, and nutrient availability) in the context of plurispecific interactions and (ii) their specific perception and subsequent signaling in A. thaliana. First, the amount of signals produced by the neighboring species A is strongly reduced in a plurispecific neighborhood, thereby leading to the strong reduction of the perception/signaling events activated in the context of bispecific interactions. Second, the production of new or modified signals in a given neighboring species is triggered by the presence of another neighboring species. Third, the simultaneous presence of different signals produced by the neighboring species A and B leads to the generation of new signaling processes. The two latter cases correspond to the production of new signals that emerged from higher-order interactions among the neighboring species (Levine et al., 2017).

While our experimental design identified a large set of genes (4,524 genes), we clearly observed an RNA process as a significant enriched biological process for most of the “phenotypic trait × treatment” combinations investigated, such as RNA processing, RNA binding, and a large variety of transcription factors, with a more pronounced proportion of BHLH and homeobox genes. These transcription factors are known to regulate multiple developmental processes, such as stress and light response (Bürglin and Affolter, 2016; Guo et al., 2021), and might exert a crucial role in promoting plant adaptation to diverse plant environments.

The second major over-represented biological process is related to the transport functions. Interestingly, 14 out of the 79 transport related proteins, identified in the seven treatments for ratio HD corresponding to the height of the first flower divided by rosette diameter, correspond to ABC (ATP-binding cassette) transporter proteins (Supplementary Dataset 2). Although none of them have been functionally characterized, these transporters have been recently proposed as key players of plant adaptation to their abiotic or biotic environment (Hwang et al., 2016). Interestingly, different ABC transporter genes were identified in the case of the bispecific interactions (ABC-B2 and G13) as compared with the plurispecific interactions (ABC-C3, C6, C7, and G26). Because of their diverse substrate specificities, they might constitute the essential components of perception/signaling pathways activated during plurispecific plant-plant interactions. Among the transport processes, nutrient foraging might also be a major response strategy in the context of plurispecific plant-plant interactions (Pierik et al., 2013).

Our observations, if considering all the enriched biological processes found in this study, are partially in agreement with the previously reported seven categories of functions identified in the artificial environments simulating plant-plant interactions: photosynthesis, hormones, cell wall modification and degradation, defense against pathogens, ABC transporters, histone modification, and meristem identity/life history traits (Subrahmaniam et al., 2018). However, transcriptional regulation appears here as a major previously unidentified process, in line with a highly flexible genetic reprogramming in response to the different plant assemblage environments.