Introduction

Many bacteria function as social groups of cells in which intercellular cooperative behaviors provide fitness benefits [1,2,3]. To define a social group, cells must maintain a coherent and cooperative social structure, which often involves excluding cells outside the group. These social barriers, in turn, lead to genetic isolation and diversification of groups. However, the genetic mechanisms by which social barriers arise are not clear. To address this question, we used wild Myxococcus xanthus isolates to investigate how they evolved social barriers within their natural habitat.

M. xanthus is a model organism that exhibits highly cooperative behaviors, ranging from social motility and coordinated predation to the formation of multicellular spore-filled fruiting bodies in response to starvation [4,5,6]. These behaviors require resource sharing among kin. However, prior studies of natural soil isolates revealed that many M. xanthus strains that reside in close proximity, form social groups that are incompatible with other groups, including between groups that are highly related (e.g. >99.99% DNA sequence identity). That is, their colonies will not merge on agar (colony-merger incompatibility, CMI) [7, 8] or they are unable to develop cooperatively in mixed culture [9]. Such intraspecific social incompatibilities can arise from variations at numerous loci that act in a complex and combinatorial manner to establish strain identity [10,11,12]. In these cases, genetic determinants of kin discrimination are so numerous that the causal genetic factors by which once-identical parental strains transitioned to socially incompatible descendants are difficult to identify. In this study, we investigated very closely related M. xanthus isolates to identify the genetic determinants and molecular mechanisms that underlie their social incompatibilities. By studying recently diverged social groups with minimal sequence dissimilarity, we were able to use comparative genomics to pinpoint the causal genetic loci of their social incompatibilities.

Bacteria use a variety of kin discrimination mechanisms to maintain social identities, including bacteriocins and secretion systems, such as the type VI secretion system (T6SS) [1, 13,14,15,16]. In addition, cellular receptors can help establish social groups by kin recognition. For instance, myxobacteria identify kin through homotypic interactions of a variable surface receptor called TraA [17,18,19,20,21]. When cells expressing identical TraA receptors make contact, they engage in outer membrane exchange (OME), whereby their outer membranes (OMs) transiently fuse and their content, including OM proteins and lipids, is bilaterally exchanged. While this process aids in membrane repair and homeostasis [22, 23], it also serves in kin discrimination by transferring a suite of diverse lipoprotein toxins called SitA [24, 25]. SitA effectors belong to six families that are not homologous in sequence but share analogous domain architectures and are delivered by OME [26]. Within a family, the C-terminal toxin modules are polymorphic whereas the N-terminal domain is conserved. A downstream allele-specific immunity gene, called sitI, whose product is not transferred, always accompanies the sitA genes. Thus, following TraA–TraA recognition and OME, the transfer of SitA toxins provides a second and potent level of discrimination. Consequently, OME between SitA-incompatible strains is antagonistic and results in CMI [25, 27]. A recent comprehensive analysis of sitAI loci distribution in myxobacteria found a strikingly high number of diverse toxins [26]. On average, genomes contain 38 distinct sitAI loci and up to 83, where we propose they serve as “self-identity barcodes” that establish highly specific strain identity and CMI. How these complex barcodes arise, however, is still unclear.

CMI in M. xanthus can also emerge as a consequence of T6SS-dependent delivery of toxins [8, 15]. T6SSs are a widely distributed nanoweapon in Gram-negative bacteria that functions as a contractile “spear” to inject an array of effectors into neighboring cells [14, 28]. T6SS effectors belong to different families, including PAAR and “evolved” VgrG proteins, where these N-terminal VgrG domains involved in delivery are fused with polymorphic C-terminal toxins [29]. Given that the expression of unique OME or T6SS toxins in one strain can abolish social compatibility with its parent [15, 24,25,26], we hypothesized that horizontal transfer of toxin-immunity loci may be the causal mechanism that establishes new social incompatibilities. To test this, we explored the genetic basis of social incompatibility in natural M. xanthus isolates.

In the mentioned study, 78 environmental M. xanthus strains were isolated from a 16 × 16 cm forest soil sampling grid (Fig. 1a) [7, 30]. Initially, their phylogenomic relationship was determined by multi-locus sequence typing (MLST) and social compatibilities were determined by CMI. Here, ≥45 of the 78 isolates formed distinct social groups [7]. Notably, incompatibilities did not simply correlate with the genetic relatedness, as some strains with identical MLST genotypes were socially incompatible [7]. In a follow-up study, 22 isolates from the two largest clades, I and V, were each sequenced and their maximum likelihood phylogeny was determined (Fig. 1b) [31]. Members within clades were highly related with very few single-nucleotide polymorphisms (SNPs) across millions of bases of aligned sequence. In contrast, members between clades were divergent and contained tens of thousands of SNPs between clade members. Despite high relatedness within clades, these 22 isolates formed 11 compatibility types (CT) based on CMI. It was previously noted that the CT groupings correlated with sequence identity in a 100–150 kb polymorphic region [31]. In our prior work [24, 25], we found that a 100–300 kb polymorphic region within lab strains, containing imperfect repeats of a prophage-like element called Mx-alpha, harbors sitAI toxin-immunity cassettes involved in CMI. Mx-alpha elements are also found in other myxobacterial genomes. Here, we reasoned that discrete mobile and polymorphic elements that carry toxin effectors, such as Mx-alpha, may be responsible for the recent incompatibilities that evolved between otherwise nearly identical strains.

Fig. 1: Spatial distribution and phylogeny of 22 M. xanthus natural isolates (figure adapted from [31]).
figure 1

a 100 plugs were taken from 16 × 16 cm sampling grid of forest soil from Tübingen, Germany, and 78 total M. xanthus strains (circles) were isolated (no more than one isolate per plug retained) and numbered from left to right (A00-A99). Isolates belonging to the largest clades, I (blue) and V (green), are marked. b Maximum likelihood phylogeny of 22 M. xanthus environmental isolates. For tree construction details see [31]. CTs were determined by CMI [7]. Isolates in bold used in this study. The average degree of sequence dissimilarity between clade I and V was ~0.76% substitutions per site, while genetic relatedness within clades was high, i.e. 0.005% and 0.0005% substitutions/site for clade I and V, respectively. TraA compatibility groups indicated [18, 19].

Full size image

In this study, we show that intra-clade incompatibilities indeed arise from polymorphic toxin effectors from three delivery systems: OME, T6SS and an understudied system called rearrangement hotspot (RHS) [32, 33]. Although clade members contain many toxin-immunity loci, the key discriminatory loci were found on prophage-like elements. We thus propose that rapid diversification of these M. xanthus social groups occurred by horizontal gene transfer (HGT) of mobile elements that contain unique discriminatory toxin effectors. In addition, long-established social identities were abruptly reprogrammed by recombination at the traA locus. Our results therefore provide clear evidence for how microbes rapidly evolve social barriers that lead to further diversification.

Results

Social CTs correlate with cell–cell antagonism

Figure 1 summarizes how local and related M. xanthus isolates form distinct CTs based on CMI assays [7, 31]. Previously we showed that one mechanism that triggers CMI was antagonism [25, 27], and here we tested whether the natural isolates that belong to different CTs actually antagonize one another. To do this, we used one representative strain from each CT for analysis (11 strains characterized in this study are bold in Fig. 1b). In our first assay, interactions between the 11 isolates were evaluated by mixing them in all possible pair-wise combinations and observing mixed colony growth on rich media. The rationale behind this simple assay is that when strains antagonize one another, colony growth is blocked compared with monoculture growth. This analysis found that interstrain antagonism was indeed widely prevalent and seemed to be more severe between clades (Fig. S1).

Secondly, to directly measure antagonism, we conducted competition experiments to see if one strain could outcompete another, and if cell death was apparent. In order to distinguish strains, one strain was labeled with a cytoplasmic fluorescent dye (see “Materials and Methods”), which did not affect fitness outcomes (Fig. S2). Figure 2a shows one result from a 24 h competition between two clade V isolates, in which the A44–A62 strain ratio decreased ~100-fold, demonstrating that A44 was inhibited by A62. We conducted competition experiments between all intra-clade isolates and in each case, strong antagonism was found, typically with one isolate markedly outcompeting the other (Fig. 2b). In nearly all cases, the presence of abnormal cellular morphologies, e.g. filamentation and cell debris, was observed (see below). Together our findings confirm that the 11 CT isolates all antagonize one another, even when intra-clade members have extremely similar genomes (Fig. 1b). We conclude that interstrain warfare was responsible for the 11 CTs originally determined by CMI.

Fig. 2: Social incompatibility between M. xanthus environmental isolates is caused by interstrain antagonism.
figure 2

Isolates belonging to the same clade, but different CT were mixed 1:1 and incubated. To distinguish populations, one strain was labeled with CFDA dye and cell types were enumerated by microscopy. a Competition between two clade V members and strain ratios plotted at indicated times. Three biological replicates, error bars indicate SD. b–d Summary tables of intra-clade interactions after 24 h incubations. Outcomes were classified as antagonistic (red squares; ≥fivefold change in strain ratio) or nonantagonistic (dark green squares; ≤fivefold change in strain ratio). “Not available”, strain construction attempts were unsuccessful. b Competition between WT isolates. c Competition among mutants where T6SS and OME was inactivated. d Same as c, expect mutations in RHS genes were also generated in A07 and A64.

Full size image

The 11 CTs correlate with distribution of unique OME and T6SS toxin-immunity loci

Our next goal was to identify the genetic determinants that cause interstrain antagonism. Previously we discovered that myxobacteria harbor an expansive array of OME-dependent polymorphic toxins that we predicted played a key role in kin discrimination between strains with compatible TraA receptors [25, 26]. We thus identified and categorized the repertoire of SitA toxins in the genomes of all 22 isolates to test for possible correlations with the 11 CTs. A wide array of sitAI loci belonging to all six families were found in all 22 isolates (Supplemental Data File). Strains that belong to the same clade mostly contain the same sitAI loci. However, each strain or CT also encoded 1–5 sitAI loci not shared with other clade members (Fig. 3a). As described below, these unique loci were found exclusively on multiple, discrete genomic islands that are polymorphic between strains that possess them. Strikingly, these unique sitAI loci strongly correlate with CTs. To illustrate this, Fig. 3b shows predicted CTs based on the SitAI repertoire from each isolate compared with CMI. In nearly all cases, predictions were in agreement with empiric CT designations. The only exception was clade V isolate A72 that shares CT02 with A30 and A44, which, according to our predictions, forms a distinct CT by itself. However, our prediction was in agreement with a more recent study, where using a modified CMI assay, A72 was in fact shown to be incompatible with A30 and A44, thus forming a unique CT [34]. In conclusion, our analysis found a perfect correlation between unique sitAI loci distribution and experimentally determined CTs, indicating these toxins may have caused social incompatibilities. However, the repertoire of sitAI loci cannot explain different CTs between isolates with incompatible TraA receptors (Fig. 3) and, accordingly, there must be another genetic factor(s) to explain those CMIs.

Fig. 3: Incompatibility between environmental isolates correlates with genomic distribution of unique polymorphic toxins.
figure 3

a Distribution of SitA toxins, belonging to six families, among 22 isolates. Black boxes, toxins unique to a CT or a few strains used in this study and ×2 indicates two copies at different loci. b Social compatibility based on empirically determined CT (left) compared with predicted CT based on toxin distribution (right columns). Red line indicates A72 forms a distinct CT from other CT02 isolates based on SitA predicted compatibility, which is also consistent with recent empirical findings ([34]; see text for details). Highlighted isolates belong to TraA-D, while others belong to TraA-B recognition group. c Distribution of T6SS toxins, belonging to two families, among the 22 isolates.

Full size image

To investigate whether T6SS plays a role in antagonism, both between and within clades, we probed their genomes for known AHH T6SS effectors (pfam14412; HNH/ENDO VII nuclease superfamily) [15, 33]. Indeed, each strain contained multiple AHH immunity/immunity loci (Supplemental Data File and Fig. S4). However, most loci were shared by all members, and were thus not responsible for the observed antagonisms, or were clade-specific (Fig. 3c). In contrast, we discovered a second putative T6SS toxin family on Mx-alpha we named tsaE (type six secretion Mx-alpha effector; MXAN_1813) (Figs. S3 and S4b). Importantly, tsaE alleles, which all contained a conserved PAAR domain (DUF4150) and a polymorphic C-terminus (Fig. S3A), varied between genomes. As we found for SitA toxins, unique sequence identity at eight tsaE loci and one AHH locus perfectly matched CT groupings (Fig. 3b, c). Moreover, all of these unique T6SS loci were located on Mx-alpha elements (see below). We therefore hypothesized that both OME and T6SS toxins were causal factors causing social incompatibilities between strains.

Kin discrimination within clades is mediated by OME (SitA) and T6SS (TsaE/AHH) toxins

The above analyses found perfect correlations between unique toxin distribution and experimentally determined CTs (Fig. 3b). This correlation included the uncharacterized TsaE toxin family, which we indeed found acts as T6SS toxins (Supplementary Information and Fig. S3). To determine whether these toxins actually contribute to interstrain antagonism, we created mutants that abolish their delivery mechanisms and conducted competition assays. For these studies, we assayed interstrain antagonism in three scenarios: (i) absence of T6SS function (both strains t6ss−), (ii) absence of OME function (one or both strains traA−), or (iii) absence of T6SS and OME function (both strains t6ss− and one or both traA−).

First, as a single test case, we competed A44 with A62 under each of the three scenarios. Here, wild-type A62 annihilated A44 and cell debris from A44 lysis was apparent (Fig. 4a, b, upper panel). When T6SS was knocked-out in both strains, the outcome was similar to the WT control. Similarly, when OME was inactivated, the competitive outcome was the same. Remarkably, however, when both T6SS and OME were simultaneously disabled, antagonism was eliminated (Fig. 4a, b, lower panel). We conclude that OME and T6SS play equally important roles in antagonism and abrogation of these systems was sufficient to eliminate social incompatibilities.

Fig. 4: T6SS and OME are major antagonistic determinants for intra-clade social incompatibilities.
figure 4

a Competition experiments between clade V isolates A44 and A62 and their corresponding T6SS and OME mutants. Bars from left to right correspond to following strain combinations: A44 WT/A62 WT; A44 OME−/A62 OME−; A44 T6SS−/A62 T6SS−; A44 OME− T6SS−/A62 OME− T6SS−; A44 OME−/A62 tsaE_5− and A44 tsaI_5/A62 OME− (see text for details). Experiments done in biological triplicate and error bars represent SD. b Micrographs show changes in A44 WT (green; phase contrast left panels) cell morphology after co-incubation with A62 WT (upper panels) and the absence of changes when corresponding T6SS−/OME− double mutants were co-incubated (lower panel). Bar, 10 µm.

Full size image

In the match-up where OME was inactivated and antagonism was driven by T6SS, A62 strongly outcompeted A44. Our hypothesis that mobilization of toxin effectors was responsible for strain incompatibility suggested that the one tsaE cluster present in A62 but absent in A44 was the determinant of T6SS-dependent antagonism (TsaE_5, Fig. 3c). To test this, we inactivated tsaE_5 in A62. Indeed, this A62 tsaE_5 mutant was unable to antagonize A44 (in the absence of OME) (Fig. 4a). In a complementary approach, expression of the tsaI_5, the predicted immunity gene located downstream of the tsaE_5 toxin, rescued A44 from antagonism by A62. We also note that A44 has a T6SS toxin (TsaE_9), to which A62 does not appear to encode immunity. Therefore, we expected that the A62 tsaE_5 mutant would be outcompeted by A44. However, A44 did not outcompete A62 tsaE_5, suggesting that TsaE_9 was not active in A44. Regardless, in the absence of OME, our results show that a single T6SS effector located on a mobile genomic island, caused strain incompatibility. Our result also rules out the possibility that other, overlooked T6SS-dependent toxins contributed to the phenotype.

Since inactivation of the T6SS and OME abolished antagonism, we tested the impact of such mutations on the CMI phenotype, to determine if interstrain antagonism causes CMI. We found when both systems were inactivated in A06 versus A92, antagonism was abolished (Fig. 2c) and the two colonies freely merged, whereas the WT parents formed a demarcation, as did the T6SS and OME single mutants (Fig. S5). Other intra-clade interactions that we tested also did not show a demarcation when OME and T6SS were inactivated, although some strains exhibited distinct morphologies that made interpreting the CMI assay difficult (Fig. S5). We conclude that toxin delivery by OME and T6SS are the primary causes of colony demarcations (CMI) and antagonism among most intra-clade strains.

Given the above results, we systematically tested each pair-wise intra-clade strain mixtures for antagonism. This analysis excluded A15 and A31, because, after repeated attempts, we were unable to make the desired mutants. Importantly, for both clades, we again found that inactivating either T6SS or OME was not sufficient to abolish interstrain antagonism. However, inactivating both systems was sufficient to eliminate antagonism in seven of the nine tested strains (Fig. 2c). The only exceptions were two clade I isolates, A07 and A64, which continued to antagonize other clade I isolates when OME and T6SS systems were disabled, suggesting these isolates have additional kin discrimination mechanisms.

RHS systems are the third and final determinant of intra-clade antagonism

Our findings involving A64 and A07 suggest these strains use a third kin discrimination system apart from the T6SS and OME toxins (Fig. 2c). To identify this third system, we directed our attention to annotated features on polymorphic islands that make up the only significant variation between these otherwise highly related strains. In this effort, we found that A07 and A64 each have a large and conspicuous RHS-repeat ORF on their Mx-alpha elements (Figs. 6 and S6), which are not present in the other clade V isolate genomes. We named these ORFs rhs4 and rhs5, respectively. Proteins that contain RHS repeats are widely distributed in Gram-negative bacteria. They are typically associated with toxin domains, and are determinants of interstrain antagonism [33, 35]. Like sitA and tsaE (Fig. 3), rhs distributions correlated to CT groupings (Fig. 5b). Thus, we hypothesized that these rhs genes were the remaining determinants of antagonism in A64 and A07.

Fig. 5: RHS plays a role in kin discrimination.
figure 5

a Mx-alpha 3 from DK1622 (MXAN_1800 to MXAN_1900) with indicated ORFs. In strains A07 and A64, 38 ORFs from Mx-alpha 3 are substituted with the two large genes (see Fig. S6 for additional details). Figure not to scale. b Table shows RHS protein distribution in 22 isolates where black boxes indicate unique toxins. Note, the rhs3 and rhs9 are in tandem on clade 1 genomes. c Domain architecture of RHS toxins. SpvB (pfam03534), VCBS-repeats (pfam13517), TcdB (pfam12256), RhsA (CO63209), RHS core (TIGR03696), TM (transmembrane helices), SP (signal peptide), PIN-like (IPR038765), CP-carboxypeptidase-like (pfam13620), TG-transglutaminase-like (pfam01841), Gln_amidase (pfam15644), and TolB-like (IPR011042).

Full size image

To investigate a possible role of Rhs4 and Rhs5 in kin discrimination, we disrupted the corresponding genes in A07 and A64 and competition experiments within clade I were repeated. Remarkably, with the available strains, all antagonism between A07 and A64 and other clade I members was eliminated when all three systems (T6SS, OME, and RHS) were disabled (Fig. 2d). This result shows that along with SitA and TsaE, Rhs4 and 5 also contribute to intra-clade social incompatibilities. We conclude that shuffling of these toxic effectors in unique combinations is what spawned social incompatibilities.

rhs4 and 5 are the only RHS toxins that are located on Mx-alpha and are distributed in unique CTs within clade I. However, along with rhs4 and 5, at least nine distinct rhs loci with diverse architectures were found in the 22 isolates (Fig. 5b, c). rhs1 is present in all 22 strains and rhs2 is present in all clade V strains. rhs2 also has a close homolog in laboratory strains (MXAN_6679; 95% identity, except C-terminal toxin domains are divergent). rhs6, 7, and 8 are found in all clade V isolates, while rhs3 and 9 are shared by all clade I isolates. Since these RHS toxins tend to show clades specific distribution they likely contribute to inter-clade antagonism.

Mx-alpha and other mobile elements are engines of social diversification

Our analysis revealed that Mx-alpha and other mobile elements played key roles in generating new social groups by harboring kin discrimination loci. To determine the distribution of the described toxins on chromosomes and their relation to mobile elements, we used the fully assembled DK1622 genome as a reference to align the draft genomes of the 11 isolates used in this study with progressiveMauve [36] (Fig. 6). As implied above, this analysis found that all discriminatory T6SS toxin loci (i.e. TsaE_3 to 10, AHH_7; Fig. 3c) were located on Mx-alpha elements (Fig. 6). Similarly, the two unique rhs loci within clade I, rhs4 and rhs5, resided on Mx-alphas linked to tsaEI and sitAI loci. For the discriminatory sitAI genes, 8 of 18 resided on Mx-alphas, while the other ten resided on three different families of prophage-like islands or on a family of polymorphic insertion elements (Fig. 6 and Supplementary Information). Some of the discriminatory genes were found in more than one isolate. For example, tsaE7 was present in A06 (clade I) and A15 (clade V) on Mx-alphas, but the linked sitAI loci were different (Fig. 6). Similarly, sitA3_1 was found in four genomes (including DK1622) on four Mx-alphas with linked t6ss or rhs loci that were all different. Unlike other isolates, A92 (and all CT-11 members) lack Mx-alpha and, accordingly, Mx-alpha associated toxins (Figs. 3 and 7; see Supplementary Information and Fig. S6 for additional details on the five genomic island families). We conclude these mobile islands drove social diversification among isolates.

Fig. 6: Unique polymorphic toxins are associated with mobile genetic elements.
figure 6

Each line represents a genome, where core genome sequences within clades are >99.99% identical, triangles depict genetic islands. Scale of genetic islands is 10× larger than of genomes. Colors represent mobile elements (islands) that belong to different families (legend; see Supplementary Information for details), as identified by BLAST homology searches of characteristic ORFs to each element, as well as conserved gene order. Numbers after underscores indicates different toxin alleles, i.e. unique specificity (see Figs. 3 and 5). Strain A64 has a number of duplicate Mx-alpha genes and, therefore, may contain two elements, demarcated by a white line. Note, SitA5_8 (not shown here) is present in all clade V strains, but A62 (shown) has a second copy located at a distinct locus. Contigs from isolates were aligned relative to the DK1622 reference genome using progressiveMauve. DK1622 contains islands with sitA (MXAN_1231 and MXAN_1899) and T6SS (MXAN_1813) toxins.

Full size image
Fig. 7: Inter-clade social compatibility reprogramed by HGT.
figure 7

a Schematic of inter-clade HGT by phage carrying the traA locus and surrounding region between natural isolates. HGT resulted in the traA allele from clade I being swapped with the traA allele from clade V in an ancestral strain. A subset of clade I members now carry the traAcladeV allele (Fig. 1b). Sequence alignment of a clade V strain (A15) compared with a clade I HGT descendent (A60) and a clade I strain (A07) with the original clade I region. Transferred region shown wherein the traA locus is 2.1 kbp. MXAN locus tags denote the relative boundaries of recombination that corresponds to the DK1622 reference genome. MXAN_6895 is the traA locus tag. Matching green and blue lines indicate identical DNA sequence between respective chromosomal regions. SNPs densities for three regions given. Figure not to scale. b Antagonism between A44 (clade V) and A06 (clade I) is neutralized by T6SS mutations, while engineered OME between these strains reintroduces social conflict. Bars from left to right correspond to strain combinations: A44 WT/ A06 WT; A44 WT/A06 T6SS−; A44 T6SS−/A06 T6SS−; A44 T6SS−/A06 T6SS− traAcladeV (strains compatible for OME). Experiments done in biological triplicate and error bars represent SD.

Full size image

T6SS is a major inter-clade antagonistic determinant

As noted, the 11 isolates belong to two TraA recognition groups (TraA-D and B), which precludes OME of SitA toxins between them. In contrast, T6SS does not discriminate between target strains. Although inter-clade relatedness is much lower than strains within a clade, we asked if inter-clade antagonisms were primarily driven by the T6SS. To this end, we conducted competition experiments between clade I and clade V T6SS mutants. In one example, where WT A06 (clade I) outcompetes WT A44 (clade V), disabling T6SS in A06 was sufficient to switch the competitive outcomes (Figs. 7b and S7). Moreover, when T6SS was inactivated in both isolates, antagonism was abolished. In three other inter-clade competitions tested, inactivating T6SS reproducibly reversed the fitness outcomes (Fig. S7). That is, the winning strain became the losing strain when its T6SS was disabled. In general, when both strains were T6SS–, antagonism was abolished in our assay. However, in some cases, signs of antagonism, such as cell lysis and decreased colony density from mixed colonies and/or a CMI persisted (Fig. S5c). Combined, these results suggest that, while T6SS is a major antagonistic determinant between clades, other toxins play a role in inter-clade antagonism, which likely includes RHS toxins (Fig. 5b and also 2d). This conclusion is not surprising, given the more distant relationship between inter-clade members (Fig. 1b), which allows more time for them to obtain divergent kin discriminatory systems.

Social compatibilities are reprogrammed by HGT of traA loci

Clade I CT-11 is a curious group because despite otherwise high relatedness, it contains isolates that belong to two TraA compatibility groups, TraA-D (A92) and TraA-B (all other isolates; Fig. 1b). Sequence comparisons revealed, strikingly, that an HGT event occurred between a clade V member and an ancestor cell that gave rise to a subset of CT-11 members. Here, the traA locus and flanking DNA from clade V replaced the corresponding sequence from clade I (Fig. 7a). DNA alignment from A07 (clade I), A15 (clade V), and A60 (clade I, CT-11), found this ~33 kb region of A60 perfectly aligns with clade V sequence; i.e. 415 SNPs plus two indels of 27 and 1 bp, were identical to clade V sequence, while only a single SNP matches to clade I DNA. Interestingly, the highest concentration of SNPs was found in traA (154 SNPs for a 2.1 kb ORF), illustrating the polymorphic nature of this gene (Fig. 7a).

We next sought to test our prediction that changing TraA recognition influences social interactions. To do this we expressed a traAB allele (clade V) in a clade I strain (A06) lacking T6SS function. By doing so, OME of the SitA toxin repertoires occurs between members belonging to different clades. In this example, where the A44 T6SS− and A06 T6SS− parent strains interacted without antagonism, the expression of TraA-B in A06 now allowed it to outcompete A44 (Fig. 7b). This result again shows that TraA is a major determinant of self-identity and social fitness [19, 25]. Moreover, among isolates, this natural HGT event had a significant impact on fitness, because from that one ancestral cell its resulting progeny now constitute 42% of clade I measurable members (5 out 12 isolates, Fig. 1). Since traA is, by far, the most polymorphic locus in this region and TraA plays key roles in social interactions, we suggest that changes in TraA recognition was the cause for these fitness gains.

Discussion

The ability of organisms to recognize and associate with kin and to discriminate against non-kin is essential when conducting cooperative behaviors. Myxobacteria participate in elaborate and resource-intensive social behaviors, such as fruiting body formation and OME, and rely on kin recognition to ensure that benefits are retained within highly related siblings. Here we investigated a key question, how diversification between highly related M. xanthus isolates evolved within a fine-scale soil habitat. Specifically, we identified the molecular mechanisms underlying social incompatibilities and how they arose. We showed that social incompatibilities emerged as a consequence of antagonism and that OME and T6SS each play major roles in intra-clade conflict. For two isolates, RHS also played important roles. For inter-clade antagonism, T6SS plays a key role and we suspect RHS is also important. OME can direct inter-clade conflict when isolates share compatible TraA receptors. For instance, a change in TraA compatibility naturally occurred within CT-11 and influenced their competitive fitness. Importantly, social diversification was driven through HGT of mobile elements harboring polymorphic toxin-immunity cassettes.

Previously, we showed that expansive families of SitA toxins govern discriminatory behaviors between strains expressing compatible TraA receptors, but not between strains with incompatible receptors [24, 25]. We postulated that these combinatorial sets of effectors act as self-identity barcodes that allow one to predict social compatibility among strains with identical receptors. Here we tested our predictions with wild isolates and demonstrated that OME indeed acts as a major kin discrimination determinant that allows social interactions to be predicted (Fig. 3). In contrast, others only used variations in traA sequences to test for correlations with CMI [34]; however none was found because this approach fails to consider the repertoire of sitAI discriminatory loci [25, 26] involved in toxin exchange between TraA compatible strains.

T6SS tsaABCEIZ gene cluster

We described a new set of myxobacteria polymorphic toxins delivered by T6SS. These clusters contain TsaE effectors with N-terminal PAAR domains and a downstream cognate TsaI immunity protein. Interestingly, most of these clusters (8/11) contained unique toxins and were found on Mx-alpha islands, suggesting they were recently acquired by HGT. In addition, we found that toxin activity depends on two upstream genes: tsaB and tsaC. TsaB contains a predicted adapter (DUF2169) that is homologous to T6SS proteins from other species [37]. Adapters act as chaperones and facilitate attachment of cognate effectors to the T6SS by binding to VgrG [38]. The gene immediately upstream of tsaE, tsaC, encodes a 3-oxoacyl-(acyl carrier) synthase protein, which may acetylate TsaE. Homologous synthases are found in Proteus mirabilis T6SS toxin clusters [39]. It is also noteworthy that the AHH effectors described here are typically found in clusters with PAAR and 3-oxoacyl-(acyl carrier) synthase genes [8, 15, 40].

Along with PAAR and DUF2169 domains, VgrG domains are frequently required for cognate toxin delivery [38, 41]. However, in our searches of clade I and V members, we only found one bona fide vgrG gene linked to a toxin gene (the previously described tsxE) [42, 43]. In contrast to the polymorphic T6SS toxins described here, this TsxE toxin was nearly invariant between all 22 isolates and the DK1622 protein. These findings suggest distinct roles of AHH and TsaE toxins versus the TsxE toxin. The former function in kin discrimination, while TsxE functions to resolve physiological conflict (heterogeneity) among genetically identical clones [42, 43].

RHS toxins

We also identified a new family, the RHS toxins, as mediators of kin discrimination. Toxins containing RHS motifs are sometimes delivered by T6SS, however, Rhs4 and 5 do not require the T6SS. Furthermore, these nine RHS proteins are large (up to 4600 amino acids) and contain signal peptides (eight of nine), suggesting they are exported to the cell surface and possibly assemble as filamentous cell surface toxins similar to CDI effectors [32, 44]. These RHS proteins show distinct domain architectures and sequence variability at their N- and C-termini. For example, Rhs2 and 4 contain C-terminal PIN-like domains characteristic of toxin-antitoxin systems with a predicted RNase activity [45], while Rhs3 contains a C-terminal amidase domain found in polymorphic toxins [33]. Future studies need to investigate the mechanism of RHS delivery, their mode of action and roles in microbial warfare.

Mobile elements as agents of social changes

Strikingly, OME, T6SS and RHS systems all contain toxin loci in highly variable Mx-alpha islands. In the case of sitAI cassettes, they were also associated with smaller prophage-like islands (not related to Mx-alpha) and insertion elements. Genome alignments showed that of the 11 isolates characterized, 10 contained variable Mx-alpha elements at a similar location found in DK1622 (between 2.1 and 2.25 Mbp) (Fig. 6). Along with the unique sitAI, tsa, and rhs loci, the Mx-alpha islands contain conserved genes homologous to bacteriophage (Fig. S6). Given their recent acquisition and prior findings that Mx-alpha transduces, without lytic growth, as small tail-less icosahedral particles (35 nm diameter), suggest they are indeed mobile [46, 47]. Since divergent myxobacterial genomes contain these elements, we propose they act as widespread agents of diversification. The resulting fluid organization of social groups may in turn help them guard against the emergence and propagation of cheater cells [48], to ensure that cooperative behaviors are maintained. In contrast, from the perspective of Mx-alpha and the other selfish elements, they can exploit their association with polymorphic toxin-immunity cassettes to provide powerful selections for their propagation and retention in populations by poisoning individuals that lack them. In essence, these cassettes act to addict and maintain the selfish elements in host genomes [25].

Prior studies revealed a striking degree of social diversity within the Tübingen soil patch, i.e. ≥45 compatibility groups among 78 isolates [7, 30, 31], and here we show these incompatibilities are driven by HGT of discrimination cassettes. This recent and striking diversification may benefit the overall fitness of M. xanthus by protecting them against exploitation and social collapse [48, 49]. However, rapid social diversification also leads to elevated levels of conflict and hence fitness loss. Future studies need to investigate the interplay between these important evolutionary forces.

Unlike intra-clade social diversification, which recently occurred, inter-clade diversification is likely multifaceted, involving mobile elements and changes in toxin-immunity specificity by mutations, deletions, and recombination between cassettes. In addition, horizontal acquisition of other discriminatory systems, such as bacteriocins and secondary metabolites, may contribute toward social incompatibility between distant strains. Accordingly, increases in spatial distances between isolates generally correlates with greater genetic diversity [50], and hence the mechanisms that govern social incompatibilities become more complex.

Environmental reservoir of mobile polymorphic elements

A striking realization from this study is that strain diversity at the genomic and social levels is driven by a wide array of mobile, prophage-like elements. As the core genome sequences within clade isolates are virtually identical, it argues that the variable islands were recently acquired and, perhaps, in some case deleted. The former scenario raises the question where did all these elements come from? A partial answer, at least, comes from the large and diverse pool of M. xanthus strains found within this Tübingen soil patch [30]. The original study successfully cultured 78 independent isolates and the two largest clades are represented here. The remaining 56 isolates thus represent smaller and more numerous clades. Figure 6 shows 25 polymorphic islands, 22 of which are unique, within the 11 isolates and presumably some or perhaps many of these islands will be found in the other 56 genomes. In addition, this soil patch contains an unknown number of other Myxococcus strains that were not cultivated, but likely contain the mentioned polymorphic islands. Since natural competence and conjugation are not known to exist in myxobacteria, we hypothesize these elements are likely transferred by phage or transducing particles. In addition, some of the toxin alleles, e.g. sitA3_1, are found on different versions of Mx-alpha, suggesting that recombination and mixing of toxin effectors occurs between elements. The mechanism(s) by which recombination occurs is unknown. Recombination between clade V and I members has also led to strain diversification. In fact the deep division between clade I members (Fig. 1b) exists because of an HGT event (i.e. transduction) replaced a ~33 kb stretch of clade I DNA with clade V DNA. Although this was not discussed in a prior study where it was instead argued recombination between clades was negligible [31], it nevertheless occurred and was important as it reconfigured self-identities by changing traA alleles and therefore the targets of SitA toxin delivery. In addition, a more comprehensive study of all 78 isolates from the Tübingen soil patch concluded traA alleles were extensively transferred across genomic backgrounds [34]. In summary, the exchange of mobile elements and recombination between clades at the traA locus played key roles in how these strains rapidly diverged into distinct and complex social groups.

Materials and methods

Growth conditions

M. xanthus strains (Table S1) were routinely grown in the dark at 33 °C in CTT medium [1% casitone; 10 mM Tris⋅HCl (pH 7.6); 8 mM MgSO4; 1 mM KPO4]. E. coli TOP10 used for cloning was grown in LB medium at 37 °C. When required for selection or induction, media was supplemented with kanamycin (50 µg/mL), oxytetracycline (10 µg/mL), zeocin (Zeo 50 µg/ml), or IPTG (1 mM). TPM buffer (CTT without casitone) was used to wash cells. For all assays, strains were grown to mid-log phase, washed, and resuspended to appropriate densities.

Cloning and strain construction

Plasmids and primers used in this study are listed in Tables S2 and S3. Insertion mutations were created by PCR amplification of ~500 bp fragment that was internal to the gene of interest, ligated into pCR2.1 TOPO TA (kanamycin), pCR2.1 TOPO XL (kanamycin, zeocin) (Life Technologies), pDW263 (zeocin), or pSWU22 (tetracycline) vectors and electroporated into TOP10. Constructs were verified by colony PCR, restriction analysis, and/or sequencing and were then electroporated into M. xanthus. Transformants were selected by their ability to grow on medium supplemented with kanamycin, zeocin, or oxytetracycline.

Inducible gene expression constructs of Mx-alpha 1 T6SS toxin cluster and A62 tsaI_5 were created by PCR amplification of the gene (s) of interest and ligation of the amplicon into pMR3487 [51] downstream of the IPTG-inducible promoter using Gibson assembly (New England BioLabs). Constructs were verified and electroporated into M. xanthus, and transformants were selected by their ability to grow on medium supplemented with oxytetracycline.

Competition experiments

Environmental isolates were grown to mid-log phase and washed in TPM before labeling one of the competing isolates with a metabolic fluorescent dye, carboxyfluorescein diacetate succinimidyl ester (CFDA) (Thermo Fisher Scientific). To label cells, they were resuspended in CTT to ~3 × 108 CFU/ml where CFDA was added to a final concentration of 20 µM. After a 30 min incubation with shaking at 33 °C in the dark, allowing the dye to be metabolically activated and trapped in the cytoplasm, the cells were washed and resuspended in TPM to a density of ~3 × 109 CFU/ml and mixed 1:1 with unlabeled cells of the competing strain. Mixtures were spotted on CTT agar containing 2 mM CaCl2, which facilitates motility and OME. After 24 h incubation (where CFDA labeling was stable) cells were harvested, washed and counted by microscopy. Typically, ~400 cells were counted per replicate. Based on labeled/unlabeled cell counts, CI was determined as the strain ratio at 24 h normalized to the starting ratio. We note that the dynamic range of this assay was limited to 100–1000-fold changes, but in many cases, the changes were actually much larger, and typically resulted in the complete absence of the losing strain. Cell morphology was evaluated on 1% agarose CTT pads containing 2 mM CaCl2. Competition experiments using DK1622 labeled with TdTomato were done similarly, but omitting CFDA labeling. All expression constructs were induced with 0.1 mM IPTG. Mutants defective in OME were made by disrupting the traA gene, while T6SS was disabled by disrupting the tssA gene (see [42] for T6SS operon structure). To prevent T6SS activity, it was disabled in both competing strains, while disrupting traA in just one strain was sufficient to block OME.

Sequence analysis

M. xanthus isolate genomes were assembled and annotated from public reads acquired from the NCBI Short Read Archive [52] by using the Comprehensive Genome Analysis Service provided by the PATRIC web server [53]. SitA and T6SS toxins were found using BLAST against the isolate genomes with known SitA and T6SS family toxins as queries [15, 25, 26]. Multiple sequence alignments were generated in MUSCLE and visualized in Jalview [54]. ProgressiveMauve [36] was used to align 11 environmental isolates’ genomes to the reference DK1622 genome. Clustal Omega [55] with default parameters was used to construct phylogenetic trees of the T6SS toxins.