Plastic multicellular development of Myxococcus xanthus: genotype–environment interactions in a physical gradient

In order to investigate the contribution of the physical environment to variation in multicellular development of Myxococcus xanthus, phenotypes developed by different genotypes in a gradient of substrate stiffness conditions were quantitatively characterized. Statistical analysis showed that plastic phenotypes result from the genotype, the substrate conditions and the interaction between them. Also, phenotypes were expressed in two distinguishable scales, the individual and the population levels, and the interaction with the environment showed scale and trait specificity. Overall, our results highlight the constructive role of the physical context in the development of microbial multicellularity, with both ecological and evolutionary implications.


Introduction
Developmental patterns, and phenotype in general, have often been considered as univocal outcomes of the genome. Thus, research has focused-conceptually and experimentally-on studying genetic mechanisms in invariable environmental conditions. However, phenotype has been shown to involve complex, bidirectional interactions between organisms and their environment. For example, light, moisture and nutrient availability in plants, as well as temperature, population density and predator presence in animals, partially determine developmental trajectories [1][2][3][4][5][6]. The phenotypic repertoire of a given genotype across a range of environmental conditions is called reaction norm [7] and its characterization constitutes the most common approximation to the study of organism-environment interactions. Since development occurs at changing and often organism-modified environmental conditions, increased survival, adaptation and phenotypic innovation can result from flexible phenotypic responses (phenotypic plasticity). Plasticity itself contributes as a cause and not just as a consequence of development and phenotypic transitions in evolution [1,5,[8][9][10].
The origin of multicellularity is a major evolutionary transition in organismal development. Plant and animal models have guided research in organism-environment interactions as well as multicellular development and evolution. However, microorganisms are largely missing in these integrative efforts, despite their ubiquity across environments and the fact that some of them can develop stereotypical multicellular structures. Microbial groups can provide new insights regarding the contribution of the organism-environment interactions to the development and evolution of multicellularity since they develop in a scale and environment similar to those in which multicellularity might have emerged [11]. Importantly, physical processes associated with environments at the microscale happen to be relevant as both driving and restrictive forces in the development of microbial multicellular phenotypes [12][13][14]. These forces and their effect on living matter vary across scales, leading to reconsideration of the meaning of environment and phenotype at different organization levels.
Myxococcus xanthus is a cosmopolitan soil bacterium that has a multicellular life history and moves by gliding across surfaces. In nutrient-rich substrates, cells divide and swarm outwards expanding the colony radially. When nutrients become scarce, M. xanthus cells come together and develop into three-dimensional multicellular structures called fruiting bodies (FBs) that contain differentiated cells (spores) [15]. The genetics underlying this process have been studied intensively [15,16]. In addition to the multicellular organization at the level of a single FB, the development of M. xanthus involves a stereotypical, but largely unexplored, spatial arrangement of FBs [17]. Given their characteristic gliding motility, cell-tosubstrate interaction is an important aspect of organism-environment interaction in myxobacteria and is susceptible to variation by mechanical properties of the substrate [18][19][20]. We therefore hypothesize that changes in the substrate stiffness may affect velocity and aggregate formation, uncovering phenotypic plasticity in different scales and traits of the multicellular structures of M. xanthus.
In the present study, we determine reaction norms of different M. xanthus genotypes in multicellular development under varying substrate stiffness environments. Our results show that phenotypic plasticity for multicellularity in M. xanthus involves multiple levels of biological organization: single FBs and FB groups.

Strains, growth and developmental conditions
To evaluate the contribution of the genotypic and environmental context to the phenotypes, five genotypes (strains) were assayed for development over five different substrate stiffness conditions, which were modified by varying agar concentrations ( parental strain: DZF1; in-frame deletion mutants: DmkapC, DmkapA/DmkapC, DmkapA, DpktC2/DpktD1, kindly provided by S. Inouye; agar concentrations: 0.5%, 1.0%, 1.5%, 2.0%, 2.5%) [21]. DZF1 is a standard laboratory M. xanthus strain. DmkapC, DmkapA/DmkapC, DmkapA and DpktC2/DpktD1 are DZF1-derived in-frame deletion mutants of the network of PSTKs and the associated scaffold proteins Mkaps, which participate during the development of FBs [16,22]. An important feature of these mutants is that they do not arrest development and thus allow tracking of multicellular phenotypic changes at this scale. Previous studies suggest that deletion of pktC2, pktD1 and other components of this network ( pktA2, pktD9) impacts development and FB formation [17,[23][24][25][26]. MkapA and MkapC are scaffold proteins in this network, but there are not previous reports regarding the phenotypic consequences of their deletion.
Following the protocol described in Yang & Higgs [15], strains were taken from frozen stocks by spotting 50 ml of each onto a CYE agar plate and incubated at 328C for 2 days. Cells from the resulting royalsocietypublishing.org/journal/rsos R. Soc. open sci. 6: 181730 colonies were transferred to 25 ml of CYE liquid medium and incubated at 328C, shaking at 250 r.p.m. overnight. Culture dilutions of each strain were grown from 0.2 OD550 until they reached 0.7 OD550 (nutrient-rich liquid culture). Prior to the bacteria development assays, cells were harvested by spinning them at 8000 r.p.m. for 5 min. The resulting pellet was washed twice with TPM solution and resuspended in 1/10th of the original volume. Fifteen microlitres were spotted onto the different agar concentration TPM plates. TPM plates were prepared by filling each with 30 ml TPM/agar media, at the appropriate agar concentration, and storing them overnight at 328C before use. After the spots dried, the plates were incubated at 328C for 96 h, until aggregates fully developed into matured FBs, recognized by their invariant shape, size and complete darkening [17]. In order to measure the motility of each strain, 15 ml of nutrient-rich liquid culture (CYE, described above) was spotted onto the different agar concentration CYE plates. After the spots dried, plates were incubated at 328 for 96 h. Each strain/substrate condition of TPM and CYE plates were conducted in triplicate for statistical support.

Measurement of phenotypic traits
Micrographs of the resulting FB development on nutrient-poor substrates, or growth on nutrient-rich substrate, were taken at 370.8 pixels mm 21 using a LEICA m50 stereomicroscope with an ACHRO 0.63Â objective lens and a Canon-EOS Rebel T3i camera. For image processing, micrographs were binarized into black/white images and phenotypic traits were measured using FIJI (ImageJ) software v. 2.0.0 [27]. For each FB image, phenotypic traits were classified and quantified at two scales: individual and population level. Single FB scale traits included Area, and Circularity, Aspect Ratio and Roundness as shape descriptors. Population scale traits, from all the FBs resulting from a single drop, included ring formation at the edge of the drop (Ring), standard distance between FBs (Distance; note that this variable might change in unexpected directions when the FBs exhibit an uneven distribution), number of FBs (Num. FB), maturation of FBs considered as complete darkening of clearly defined aggregates (Maturation), presence of vestiges of not clearly defined aggregates (Scars), position in the sequence of maturation speed, being 1 the fastest agar % and 5 the slowest agar % for each genotype (Dev. time), and mature FBs beyond the edge of the drop (FB out). Outer FBs are formed during development from cells that migrated beyond the initial drop edge. In order to analyse the images systematically, we cropped the FBs placed outside the drop, but recorded their presence under the categorical variable 'FB out' (figure 1; electronic supplementary material,  [28]. The edge of the drop was defined at time 0 h, when it is clearly seen, but it did not expand or change its position during development on nutrient-poor media. For assays of growth on nutrient-rich substrate, the diameter of colonies was measured in micrographs taken from 0 to 96 h, every 24 h. We calculated the growth rate as the percentage that the colony diameter increased over 96 h. For this study, this rate was used as a proxy for motility. These measurements were performed using FIJI (ImageJ) software v. 2.0.0 [27].

Data analysis
In order to test contributions of genotype, substrate and their interaction, a PERMANOVA analysis of a representative sample was performed using the adonis function in the R package vegan (10% of the total data, N ¼ 2400) (table 1) [29]. This is a multivariate test suitable to distinguish sources of variation in nonparametric datasets. To investigate the phenotypic differences in micrographs of figure 2, phenotypic traits were grouped by genotype and by agar concentration, and a factorial analysis of mixed data (FAMD) was performed for each group through the FAMD function in the R package FactoMineR [30]. This multivariate analysis allows us to analyse numerical and categorical variables and generates a multidimensional space where the first and second axes explain the largest proportion of variances as a combination of phenotypic traits (figure 2b,c; electronic supplementary material, table S2).
To characterize the trait-specific phenotypic variation, reaction norms were constructed connecting median values for each phenotypic trait. A linear regression fit was superimposed to visualize tendencies, although p-values were not significant in all cases (figure 3a; electronic supplementary material, figure S2 and table S3). A second set of reaction norms based on the coordinates of the first and second axes of the FAMD analyses as a statistical synthesis of phenotypic trait variation was also performed (figure 3b).
Finally, to test if the motility was altered in response to agar concentration or in response to genotype, and whether such variability might contribute to the observed developmental plasticity, we assessed colony growth over nutrient-rich substrates measuring the increase in its diameter and plotted in function of agar concentration (figure 4a) and of time ( figure 4b). These values were normalized according to the colony diameter at 0 h. Then, we plotted the motility rate considering the proportion in which the colony diameter increases after 96 h and performed a two-way ANOVA to test the statistical significance of genotype and substrate conditions in such rate (figure 4c). Note that FBs do not develop in the presence of nutrients.

Results
We observed that most of the phenotypic variation revealed in micrographs of figure 2a is explained by genotype (adj R 2 ¼ 0.571), substrate stiffness (adj R 2 ¼ 0.166) and their interaction (adj R 2 ¼ 0.185). Statistical significance for each factor is strongly supported by the PERMANOVA analysis: pGenotype ¼ 9.9 Â 10 24 , pEnvironment ¼ 9.9 Â 10 24 , pInteraction ¼ 9.9 Â 10 24 (table 1). The FAMD analyses provided clarity in discriminating phenotypes. Specifically, by fixing each substrate condition (agar %) and assessing the differences among the genotypes, phenotypic differences can be distinguished (figure 2b). For example, all the genotypes (ellipses) are well distinguished at the lowest agar concentration (0.5%), but tend to overlap as agar concentration increases. However, the parental strain ellipse (DZF1) is clearly apart from that of the other phenotypic variants, except for 0.5%. Looking across agar concentrations, the arrangement of the genotypes follows equal directions along the X-axis (axis of major variation) in all agar concentrations, except for 0.5%, with larger contribution from the number of FBs and the standard distance between them (electronic supplementary material, table S2 and figure S3). Finally, while the DmkapA/C and DmkapC phenotypes overlap in all the cases, correlating with their similar phenotypic profiles shown in the micrographs, there is a non-additive effect of the double mutant, as it does not resemble the DmkapA phenotype. By contrast, we observed more overlap when assessing agar concentrations (ellipses) for each genotype (figure 2c). There is not a clear sequence in the ellipse disposition and overlapping is present between almost all of the ellipses, except for 0.5% ellipse, which is well distinguished in all genotypes. We also did not observe a pattern of traits contributing to X-and Y-axes within each genotype, although ring formation is predominant in all cases (electronic supplementary material, table S2 and figure S3). Altogether, our results show that at 0.5% agar concentration, there are clear phenotypic differences within and among genotypes, and that phenotypic identity is blurred as agar concentration increases. Furthermore, the consistency of variables explaining the phenotypic variation at X-axis (Dim. 1) and the constancy in the sequence of genotype ellipses when considering each agar concentration suggest that the expressions of the genotype are enabled in a given environmental range. These patterns are notably robust across replicates, which is evident from the marginal   figure S3). Reaction norms show that there is substantial phenotypic variation among genotypes for individual traits across agar concentrations (figure 3a). Each traitÂgenotype combination has its own pattern, for both the single FB and population scales. For example, in DmkapA, 'Circularity' and 'Distance' are almost linear but with opposed slope signs, while 'Area' does not show a linear behaviour.  Remarkably, threshold variation is observed for some traits with particular condition and strain combinations, such as 'Ring' at the 0.5% condition. These effects are not general, especially for the parental strain (DZF1), in which the agar concentration is less strong in affecting 'Ring' and has no effect in FB traits and 'Dev. time', except at 0.5% agar. Finally, the large effect of substrate conditions, especially at the extreme ones (0.5% and 2.5%), is evident when plotting the FAMD coordinates as a statistical synthesis of phenotypic trait variation (figure 3b).
Regarding the colony growth rate in nutrient-rich substrates, there are strong interaction effects of agar concentration. At the lowest and highest agar concentrations, there was almost no difference among genotypes for colony growth rate ( figure 4b,c). We observed a much lower motility on 0.5% compared to the rest of the agar concentrations (figure 4b,c), in line with previous reports [34]. Overall, when we track the colony growth, we observe that the strains slightly diverge according to their growth and that, as development proceeds, strains form two groups (DZF1 and DmkapA form one group and the rest of the strains another one ( figure 4a,c)). However, we found that in such nutrient-rich substrates, the agar concentration is the most significant factor variable for the colony growth rate (see ANOVA results in figure 4c).

Discussion
We experimentally measured reaction norms of bacterial multicellular development across substrates with increasing stiffness. Phenotypic plasticity in M. xanthus development was revealed, involving changes in cell-to-substrate interaction due to increasing agar concentration in culture media plates. In this study, we focus on the potential for changes in morphology during development. These include changes in development time, spore maturation and FB size and number, which are important fitnessassociated traits upon which selection may subsequently act. Previous studies with the same regulatory network reported statistically significant phenotypic differences when comparing mutants  DmkapC with parental strain (DZF1), including differences in paradigmatic functional traits, such as spore count and viability [17]. Our current results, and those from previous studies, suggest that the joint effect of the physical environment and genetic background on morphologic and developmental traits is likely to be of evolutionary relevance. The experimental design allowed for statistical discrimination of phenotypic changes at two scales: 'FB scale' and 'population scale'. In general, for each combination of genotype and substrate condition, FBs conforming the population exhibit little variation in their shape and size (figure 2a). Nevertheless, we found genotypic and substrate effects on phenotype at the 'population scale', which have been overlooked or dismissed when considering variation only at 'FB scale'. Remarkably, phenotypic identity over the genotypic and substrate conditions is not observed when considering only the FB scale (electronic supplementary material, figure S1; versus figure 2b,c).
To our knowledge, mechanisms behind the organization at the FB population scale have not been explored, although it is known that mechanical forces and living matter interact in a bidirectional way giving rise to robust patterns [13,14]. Based on our results and taking into consideration the influence of mechanisms such as substrate fibres alignment in colony formation [18] and the surface-tensiondriven coarsening process of aggregation [35], we can suggest that the mechanical properties of the medium constitute a key element when thinking about environment at microscales. Also, substrate changes could determine the individual and collective motility, which could in turn drive aggregation. In this case, the developmental process occurring in nutrient-poor substrates would be mainly determined by the agar concentration. We found that the growth rate is largely determined by the agar concentration, with slight differences among strains ( figure 4a,b). This is especially clear for the strain-independent low and high growth rate at the 0.5% and 2.5% conditions, respectively. Nevertheless, overall, the multicellular phenotypes at the FB and population scales are largely explained by agar concentration, genotype and the genotype -substrate interaction (table 1), which reflects the complexity of the developmental process under study.
Our results also demonstrate the value of investigating phenotypes from multi-level perspectives [11]. Phenotypic variation existed at both the level of single multicellular structures and the collective level of groups of these structures. It is important to note that specific environmental conditions have trait-and scale-specific effects. For instance, the substrate seems to play a much stronger role in certain ranges of agar concentration and for particular traits. We report a strain-independent ring formation, at population phenotype, surrounding the initial drop area at 0.5%. This could be explained by a higher cellular density at the edge of the drop forming when it dries at time 0 h, as in the coffee-ring formation physical phenomenon [36]. It is possible that at this agar concentration, cells cannot glide inwards easily, developing a ring of FBs at the edge. This hypothesis, and the precise mechanisms behind different phenotypes, remain to be tested.
In this work, we considered several points within the widest possible range of substrate modification which enabled a good approximation to the shape of reaction norms. We initially included a broader range of agar concentrations (electronic supplementary material, figure S4), but excluded those values where there was no development of FBs. This allowed us to notice that phenotypic change is usually greater at the extremes of the agar concentration range (figure 3b; in good agreement with previous theoretical proposals [37]). However, some limitations to this work and to the current study of plasticity in general should be considered. For instance, we modified stiffness by varying agar concentration, but other variables affecting the substrate properties such as water availability could be correlated. Indeed, reaction norms have provided useful perspectives regarding the paired relationship between a trait and an environmental factor for a specific moment in the developmental trajectory. However, new approaches considering the dynamical interaction among numerous phenotypic and environmental variables, across time and in ecologically meaningful ranges, are required. Moreover, plasticity is often studied in laboratory-adapted strains, and at least in this work, the parental laboratory-strain is not representative of the plastic responses for the rest of the genotypes (figures 1 and 2). Completely unexpected phenotypes may arise when considering natural environments or non-domesticated strains.
Overall, our study highlights the importance of the physical environment in the development of robust, yet plastic, aggregative microbial structures in different genotypic contexts. This in turn calls for systematic approaches to integrate the different factors and mechanisms (genetic and environmental) behind phenotypic plasticity and its consequent role in the emergence and evolution of multicellularity.