Spatial patterns in phage-Rhizobium coevolutionary interactions across regions of common bean domestication

Bacteriophages play significant roles in the composition, diversity, and evolution of bacterial communities. Despite their importance, it remains unclear how phage diversity and phage-host interactions are spatially structured. Local adaptation may play a key role. Nitrogen-fixing symbiotic bacteria, known as rhizobia, have been shown to locally adapt to domesticated common bean at its Mesoamerican and Andean sites of origin. This may affect phage-rhizobium interactions. However, knowledge about the diversity and coevolution of phages with their respective Rhizobium populations is lacking. Here, through the study of four phage-Rhizobium communities in Mexico and Argentina, we show that both phage and host diversity is spatially structured. Cross-infection experiments demonstrated that phage infection rates were higher overall in sympatric rhizobia than in allopatric rhizobia except for one Argentinean community, indicating phage local adaptation and host maladaptation. Phage-host interactions were shaped by the genetic identity and geographic origin of both the phage and the host. The phages ranged from specialists to generalists, revealing a nested network of interactions. Our results suggest a key role of local adaptation to resident host bacterial communities in shaping the phage genetic and phenotypic composition, following a similar spatial pattern of diversity and coevolution to that in the host.


Introduction
Bacterial viruses or bacteriophages are the most diverse and abundant biological entities on earth [1,2]. They play a significant role in bacterial ecology and evolution, enabling horizontal gene transfer, and influencing bacterial diversity through their lytic or lysogenic cycles [3][4][5][6]. Nevertheless, phage biogeography, including patterns of dispersal, establishment, and community assembly, is very poorly understood [7], and its study could contribute to the advancement of microbiome engineering in agricultural and medical settings.
Bacteriophages tend to be locally adapted to sympatric bacteria [25][26][27][28]. It is commonly predicted that local adaptation is a significant underlying factor of compositional differences across phage communities [16,24,26,27,29]. Local adaptation is a process that results in a local population of a given species exhibiting higher fitness in its local environment than in allopatric populations and takes place when environmentally driven selection is more robust than migration [30,31]. The implied lower fitness experienced by immigrants may limit gene flow and increase genetic differentiation across populations (i.e., "isolation-by-adaptation," or more general "isolation-by-environment") [32][33][34][35]. Host bacteria are a key environmental driver of phage local adaptation, or more precisely, "local coadaptation", as phage adaptation (e.g., increased host range) readily invokes the counter-or coadaptation of bacteria (e.g., increased resistance) [25][26][27][36][37][38].
While bacteriophages coadapt with their bacterial hosts, symbiotic bacteria coadapt with their eukaryotic hosts [39][40][41]. For instance, it is well established that rhizobia coadapt with legumes and form a mutualistic relationship in which they provide legumes with a steady supply of plantusable nitrogen via nitrogen fixation [42,43]. Aguilar et al. [44] found that Mesoamerican and Andean common bean (Phaseolus vulgaris) genotypes were preferentially associated with Mesoamerican and Andean rhizobia, respectively, indicating local coadaptation [44] across the two domesticated common bean gene pools [45,46]. Furthermore, the Rhizobium communities were genetically differentiated, as the relative abundance of different types of the symbiotic nodC gene varied across the two bean gene pools [44]. In addition, in other legume-rhizobium systems, host legume population differentiation has led to rhizobium population differentiation [47,48] and even local adaptation [49,50].
Symbiotic organisms have significant reciprocal evolutionary effects on each other, which in turn affect third-party interactions [51][52][53][54][55][56][57][58][59]. Bacterial adaptation to plant hosts may affect bacteria-phage interactions (e.g., by altering the expression of surface receptors for plant interactions, which serve as anchor points of phages [60][61][62]). Here, we predicted that the genetic differentiation and local adaptation experienced by rhizobia across common bean gene pools shape the interaction and genetic differentiation of communities of phages infecting common bean-nodulating rhizobia. Even under the assumption of high phage dispersal capabilities, we expected communities of phages to be locally adapted, showing higher infection rates for sympatric rhizobia than for allopatric rhizobia. Hence, we aimed to elucidate whether the genomic identities and host ranges of phages infecting common bean-nodulating rhizobia are geographically structured. Our approach was to sample Rhizobium strains and associated bacteriophages from two common bean fields in Mexico and Argentina, corresponding to independent areas of common bean domestication [45,46]. Host species identity and bacteriophage genomic types were determined by Sanger sequencing and by both Sanger sequencing and wholegenome sequencing, respectively. Host range assessment results were analyzed for biogeographic signals and local adaptation. More specifically, we aimed to (i) assess bacteriophage diversity associated with common beannodulating rhizobia; (ii) compare the bacteriophage community composition across common bean gene pools; (iii) determine how bacteriophage host range pertains to bacteriophage provenance and the host species community composition; and (iv) obtain evidence of local adaptation of bacteriophages.

Soil and Rhizobium sampling
The sites of rhizobia and phage sampling were common bean (Phaseolus vulgaris) agriculture fields located in Mexico (Tepoztlán and Yautepec, Morelos) and Argentina (Chicoana and Salta, Salta) (Fig. S1a). The bean fields were named after the municipalities in which they were located except for the 'Salta' field, which was located in Cerrillos near the border with the city of Salta. Sampling was done during spring in both regions to minimalize effects of seasonal variation: Tepoztlán was sampled on March 2016 and Yautepec on May 2016 (Northern hemisphere spring); Argentinean fields were sampled on October 2016 (Southern hemisphere spring). The distance between Tepoztlán and Yautepec is c. 7.4 km; the distance between Salta and Chicoana is c. 20.8 km. In each bean field, we collected >1 L of rhizosphere soil from three plants separated by 10 m. Each soil sample was mixed and split into two aliquots, one used for phage isolation and one used to "trap" rhizobia in the root nodules of bean plants. In Tepoztlán, three bean plants (P. vulgaris var. Negro Veracruz) were collected directly from the field, and their nodules were used to isolate rhizobia. To trap rhizobia from Yautepec, we used P. vulgaris var. Negro Veracruz, while rhizobia from Argentina were caught using P. vulgaris var. Alubia Cerrillos. These cultivars belong to the Mesoamerican and Andean gene pools of domesticated common bean, respectively. Previously, axenically germinated seedlings were planted in presterilized vermiculite pots and inoculated with a 100 mL soil sample. Nodules were harvested after 8-12 weeks of growth in a greenhouse under natural environmental conditions regarding light and temperature (Supplementary Method S1). Rhizobium strains were obtained from surface-sterilized nodules via the squashing method and streaked on PY Nal agar plates (Supplementary method S2; [63]). Three subculture steps were used to purify all isolated rhizobial strains. The isolates were stored at −80°C in 50% glycerol for long-term storage.

Isolation and purification of bacteriophages
Bacteriophages were isolated from soil samples via the enrichment protocol described previously [64] using both Rhizobium isolates from local sites (local collection, or LC) and Rhizobium from the laboratory collection (standard collection, or SC; Table S1). Briefly, dry sieved soil was suspended in PY-Nal medium at a ratio of 1:2 and incubated overnight at 30°C with shaking (250 rpm). Subsequently, the soil solution was centrifuged (10,000 g, 10 min, 4°C) to remove large particles, and chloroform was added to eliminate the remaining bacteria. Chloroform was removed by centrifugation, and the solution was filtered through a 0.22 µm Millipore filter. The soil filtrate was used to inoculate one pair-member of each Rhizobium strain. The other pair-member served as an uninoculated control. LC rhizobia were inoculated with the filtrate from the soil of origin, while SC rhizobia were inoculated with soil filtrate that was pooled from each of the three soil samples per bean field. Cultures of 1 ml in PY Nal were incubated at 30°C (250 rpm) to an optical density of 0.2 at 620 nm (OD 620 ; Beckman DU650 spectrophotometer) in 96-deep-well plates. After 20 h of incubation, the cultures were centrifuged, and the supernatant was used to reinoculate new cultures of the same pairs of rhizobial strains. This enrichment process was repeated five times; each time, the OD 620 value was compared between the control and the inoculated pair member. A decrease in the OD 620 was indicative of cell lysis. The presence of phages was confirmed by plaquing a dilution series of the filtrate mixed with 2.5 mL of molten soft PY (0.65% agar) on lawns of rhizobia over solid PY Nal plates, as in Carlson (2005) [65]. Lytic plaques were picked and inoculated again in the respective bacterial cultures. Two additional dilution series were performed to ensure the purity of the bacteriophage stocks. Finally, the phage dilution was treated with chloroform to remove the remaining bacteria and stored at 4°C.

Phylogeny and taxonomic Rhizobium species identification
Two chromosomal genes, dnaB and recA, were partially amplified by colony PCR. The primers and PCR protocols used are described in Table S2. The PCR products were purified using the Exo-SAP cleanup protocol [66]. The purified products were sent to Macrogen for sequencing (Macrogen Inc, Seoul, Korea). Sanger sequences were edited and assembled in Genius Pro v. 6.1.2. Multiple alignments were performed using MUSCLE [67], followed by manual correction to remove ambiguously aligned regions. Phylogenetic trees were reconstructed and edited with MEGA 7 using the maximum likelihood (ML) method based on the GTR + G + I model and 1000 bootstrap replicates. Rhizobia were clustered into sequence types (STs) based on 100% sequence identity of dnaB and recA sequences.

Genome sequencing
We obtained the genome sequences of 100 phages from the collection chosen according to their host range differences. Phage genomic DNA was purified from phage stocks propagated in the corresponding host Rhizobium strains. Host DNA and RNA were eliminated using DNase and RNase. Subsequently, phages were precipitated using a PEG-8000/ NaCl solution. After centrifugation (10,000 g, 20 min, 4°C), the pellet was suspended in Tris-EDTA buffer. Proteins were hydrolyzed using 4% SDS and proteinase K and precipitated by adding 3 M potassium acetate. Phage genomic DNA was precipitated with 100% isopropanol and washed with 70% ethanol twice.
The remaining 96 phages were identified via PCR and Sanger sequencing of phage marker genes used to distinguish between phage genomic types (PGTs). Each of these genes was shared among members of the corresponding PGT, and showed little to none nucleotide sequence similarity to genes of other PGTs; this was assessed by the presence-absence profile obtained with the BPGA software [72]. The extracted sequences were used to design primers with Primer3 [73]. The selected genes, primers, and PCR protocol are described in Table S2.

Comparative genomics
The average nucleotide identity (ANI) of all pairs of phage genomes was calculated with pyani v.0.2.9 using the ANIm MuMmer method [74,75]. PGTs were clustered based on 80% nucleotide identity and 60% coverage of the smaller genome. The assignation of PGTs to the phage morphological families of tailed bacterial viruses (Siphoviridae, Myoviridae, and Podoviridae) was performed using the VirFam server (http://biodev.cea.fr/virfam/) [76]. PGTs assigned to Microviridae were recognized by their short genome length and through BLASTn searches against the virus database of the NCBI (https://www.ncbi.nlm.nih.gov/ genome/viruses/).

Rhizobium susceptibility and host range assessment
The infectivity of all isolated phages was evaluated by the spot test procedure [77] in the rhizobia isolated from all four bean fields. Although this method of assessing the spotting host range could overestimate the lytic properties of phages and underestimate infections due to low phage titers, different adsorption efficiencies, or alternative infection modes (e.g., pseudolysogeny) [78][79][80][81], the spotting of phages on lawns is an affordable method for testing a large number of pairwise interactions [82]. Rhizobium 'lawns' in double-agar-layer plates were spotted with 5 µl of each bacteriophage solution prepared from the respective phage Rhizobium lysates. After overnight incubation at 30°C, the plates were assessed for lysis. At least three replicates of each Rhizobium-phage combination were performed to ensure reproducibility, and at least two replicates with the same lytic or resistant phenotype were considered to indicate a positive or negative result. Spots resulting from lysis with a translucent appearance rather than a transparent appearance were recorded as 'partial lysis' but were treated equivalently to transparent spots for the statistical and BiMat analyses. The binary interaction matrix is available from Github (github.com/jvancau/interactiondata).
With this information, a matrix was constructed based on Bray-Curtis dissimilarity calculated with the vegdist function of the vegan package in R [83]. To compare the phenotypic composition with the genetic composition, we clustered the bacterial hosts showing >80% similarity in terms of their susceptibility to phages into Rhizobium phenotype groups (RPGs , Table S3). Then, phages that infected a common range of hosts (Bray-Curtis similarity >80%) formed phage phenotype groups (PPGs , Table S3).

Network structure of phage-bacterium interactions
To analyze the modularity and nestedness properties of the phage-bacterium bipartite network, we employed the BiMat program [84], which maximizes the similarities in the bacteria-phage lytic interaction matrix. The program ran in the MATLAB environment and was used according to the author's start guide [84] (https://www.github.com/cesar7f/ BiMat). Modularity was tested using the Adaptive Brim algorithm, and nestedness was tested using the NODF (Nestedness metric based on Overlap and Decreasing Fill) and NTC (Nestedness Temperature Calculator) algorithms. Statistical analysis was performed using 1000 replicates and the equiprobable null model.

Local adaptation
Rhizobium susceptibility rates were calculated for each Rhizobium strain as the number of phages able to infect the strain divided by the total number of tested phages. Similarly, phage infection rates were calculated as the proportion of rhizobia that a given phage could infect among the total number of strains tested. These values were calculated for sympatric and allopatric phage-strain combinations. The statistical significance of the differences between sympatric and allopatric susceptibility/ infection rates was tested using a generalized linear model in R with a quasi-binomial model [85]. Lower susceptibility rates with sympatric phages than with allopatric phages indicate rhizobium local adaptation; higher infection rates on sympatric rhizobia than on allopatric rhizobia indicate phage local adaptation.
Phage local adaptation was also calculated according to Vos et al. [27] as the mean difference in the mean of the rates of sympatric (S) and allopatric (A) phage infections.

Statistical analyses
To test the significance of compositional differences in Rhizobium STs and PGTs across common bean fields, we performed PERMANOVA tests based on Jaccard and Bray-Curtis distances with 999 permutations using the ecodist and vegan packages [83]. Principal coordinates of neighbor matrices (PCNM), which are orthogonal spatial variables derived from a spatial distance matrix, were calculated from the geographical coordinates using the 'pcnm' function of the vegan package for R [86]. The first PCNM value was used as a proxy for spatially related variation across the two regions and was fit on a principal coordinates analysis (PCoA) using envfit ('vegan') with 9999 permutations. This was approach was employed to assess the significance of the distance between the two regions in explaining the compositional differences among bean fields. For each distance matrix type, the correlations among all datasets of the genetic and phenotypic composition across bean fields (ST, PGT, RPG and PPG) were assessed using Mantel tests (999 permutations, vegan package).
A Mantel test was also used to correlate Rhizobium and phage genetic distances with the corresponding susceptibility range and host range distances across all isolated rhizobia and phages. PERMANOVAs were used to test the significance of the effects of the origin and genetic or taxonomic identity of rhizobia and phages on their susceptibility range and host range, respectively.

Rhizobium and phage community sample composition
We studied 229 Rhizobium strains isolated from four agricultural plots in Central Mexico (Tepoztlán and Yautepec) and Northwest Argentina (Salta and Chicoana) ( Fig. S1a-b). The rhizobial strains from each site (LC, local collections) were employed to trap phages from the soil of the same locality by the enrichment method [64]. In parallel, the same soils were pooled for each plot and used to search for phages using a SC of 94 Rhizobium strains of diverse geographic origins maintained in our laboratory (Table S1). A total of 196 phages were obtained with this protocol, 110 from LC and 86 from SC (Fig. S1b).

Genetic differentiation of Rhizobium populations between agricultural fields
Our first aim was to investigate whether the collected Rhizobium strains represent geographically structured populations. Phylogenetic analysis of the partial recA-dnaB sequences of 229 rhizobial strains identified R. etli as the predominant species at Mexican sampling sites (74.6%), while R. phaseoli was dominant in Argentina (80%) ( Fig. 1;  Fig. S2a). According to the nucleotide variations in recA and dnaB, all of the isolated rhizobia were grouped into 41 chromosomal STs (Table S4). PERMANOVA tests employing two distance measures (Jaccard and Bray-Curtis) indicated that the rhizobial communities differed significantly in terms of genetic composition and the relative abundance of genotypes (STs) among different bean fields and regions ( Fig Table S4-S5). The geographic distance between the regions, represented by a PCNM vector, was correlated with the ST composition among the rhizobia isolated from the four common bean fields (Table S5). At the species level, R. etli and R. phaseoli STs differed significantly among the sites of origin of the common bean fields. Between regions, only the R. etli populations differed in terms of the ST composition, whereas the R. phaseoli ST composition changed only marginally across regions (Table S5). The abundance of ST-5 in R. phaseoli, the only ST of the 41 total STs found across the four common bean fields, might explain this last result ( Fig. 2C; Table S4).

Diversity of phage communities
To assess phage diversity, we obtained the complete genome sequences of 100 out of 196 phages (Table S6). The phage genomes displayed a wide range of lengths (from 4.8 to 207.6 kb; median 54.4) and GC contents (from 41 to 61%; median 57%). They were clustered into 29 PGTs (defined by ANIm, see Methods), 18 of which had two or more individual genomes, and 11 were singletons. Within PGTs, the genomes exhibited nucleotide variation ranging from 85.8 to 99.7%, with a coverage of~64.1 to 100% (Fig. S3). In addition, 90 phages among the remaining 96 phages were assigned to PGTs by the PCR identification of specific phage genes conserved within the members of the 29 PGTs, followed by Sanger sequencing (see Methods) (Table S2).

Phage population differentiation across agricultural sites
Following the differences in the rhizobial ST composition per sampling site, we found that the composition of PGTs also differed significantly among common bean fields and regions ( Fig. 2 D Table S7). The differences were significantly correlated with geographic distance (Table S5). The phage communities were also significantly different between the Mexican bean fields within regions based on Jaccard distances (F 1,5 = 3.7, p = 0.024), but they were not significantly different between Mexican bean fields based on Bray-Curtis distances (F 1,5 = 3.7, p = 0.055) or between Argentinian bean fields (Jaccard: F 1,5 = 1.1, p = 0.384; Bray-Curtis: F 1,5 = 2.2, p = 0.149). Mantel tests showed that the genetic composition differences among phage communities were significantly correlated with the differences in the Rhizobium community genetic composition among bean fields (Table S5). Some PGTs coexisted at two of the sampling sites, whereas PGTs rarely coexisted at three sites and never at four (Fig. 2F). Moreover, 52% of the 29 PGTs occurred solely in one bean field, and 69% were restricted to a particular region. A spatial pattern distinction was also shown by the ANI values, as the average ANI of allopatric phages belonging to the same PGTs was 88%. In comparison, the average ANI of sympatric phages belonging to the same PGT was 96%.

Structure of the phage-bacterium interaction network
To define the rhizobium-phage interactions within and between the four communities, we tested the infectivity of 196 phages against 229 Rhizobium strains by the spot assay (see methods). We registered the following phenotypes: complete lysis (transparent spots), partial lysis (translucent spots), and resistance (absence of lysis), in three independent experiments.
A total of 44,884 interactions were examined in triplicate experiments; 19,474 plaques showed full or partial lytic phenotypes, recorded as positive interactions (1 in the bipartite network of Fig. 3) [84,87]. The Rhizobium-resistant phenotypes were recorded as 0 in the corresponding binary matrix (Fig. 3).
Overall, the BiMat network showed high connectance (0.43), indicating that there were~4/10 effective phagebacterium interactions in the community, and weak modular differentiation (Qb = 0.21; Fig. 3A). Three large modules, with high internal connectance in comparison with the outside modules, were detected (Fig. 3A). This means that to a certain extend the rhizobia and phages in these modules preferentially interact with each other rather than with members the other modules. Phages from Mexico (Tepoztlán and Yautepec) generally interacted best with Rhizobium isolates from Mexico. It was less common for these phages to infect Rhizobium from Argentina. These results suggest large-scale modularity dominated by sympatric interactions. However, some exceptions were observed; for instance, phages from Chicoana were very promiscuous, since a to the outermost circles: taxonomic classification of rhizobia, Rhizobium chromosomal STs (sequence types; no STs were assigned to SC strains), field of origin of the strains, and squares indicating the origin of the phages isolated using the corresponding strain. Bootstrap values are shown with proportionally-sized gray circles on the tree branches; smallest circles equal 75% and the largest circles equal a value of 100%. significant number of Mexican strains were cross-infected by these Argentinian phages (Fig. 3A). This may in part explain the overall nested structure shown by the BiMat matrix, quantified by the independent NODF (0.68) and NTC (0.73) algorithms. This suggests that the phage communities consist of a range of specialist to generalist phages (Flores et al. 2013). This is visualized in Fig. 3B, where phages on the left side are more generalist, while further to the right phages become increasingly specialist.

Rhizobium susceptibility and phage host range
Although all Rhizobium strains were closely related in the recA-dnaB phylogenetic tree, they varied widely in their susceptibility range, with some rhizobia being infected bỹ 10.2-74.0% of the phages (average rate of infection = 43.4%). Hence, we sought to investigate whether the susceptibility of rhizobia was associated with geographic origin, taxonomic affiliation or genetic identity. When we compared the variation in the susceptibility range between individual rhizobia, we found that the susceptibility of rhizobia was significantly different not only between species (F 1,217 = 10.3, p = 0.001) (Fig. 4A) but also among bean fields (F 3,217 = 31.0, p = 0.001), and we observed the interaction of the two factors (F 3,217 = 6.4, p = 0.001; Fig. 4A). The last finding indicates that the effect of Rhizobium species identity on the susceptibility range depends on the geographic origin of the rhizobia. When species were split into STs, the Rhizobium susceptibility range was found to differ significantly among bean fields (F 3,178 = 53.5, p = 0.001) and STs (F 40,178 = 6.1, p = 0.001), although their interaction was not significant (F 6,178 = 0.9, p = 0.542) due to the limited geographic spread of most STs. Susceptible strains were clustered in 48 RPGs (Rhizobium Phenotype Groups; see methods). Twenty-nine RPGs were singletons, whereas 19 were formed by more than two strains (maximum of 40) (Table S3). Thirty-two RPGs belonged to R. etli, twelve RPGs corresponded to R. phaseoli (RPG 4,13,19,20,22,23,24,27,31,32,33 & 37), and the other four RPGs belonged to both species (RPG 1, 2, 6, and 11). The  most abundant RPGs (1-4) were composed of the frequently identified STs of R. phaseoli (ST-5) and R. etli (ST10 and ST34). We found that the RPG composition was strongly correlated with the ST composition (Table S8). In addition, Rhizobium genetic distance and susceptibility range similarity were significantly correlated (r = 0.31, p = 0.001).
We examined the host range of phages by the spot method, to identify those with similar infection spectra. Phage infection rates, or the proportion of hosts that a given phage isolate could infect, varied considerably from 2.2 to 92.6% (average infection rate 43.4%). All but three phages were able to infect both R. etli and R. phaseoli; however, the phages were able to infect only 51% of STs on average. The host range of the phages was significantly affected by the bean field of origin (F 3,142 = 23.0, p = 0.001), phage genomic type (PGT; F 27,142 = 6.9, p = 0.001), and the interplay between these factors (F 14,142 = 2.9, p = 0.001). The last finding indicates that the effect of PGT on the host range depends on the geographic origin of the phage. Similar results were obtained when the phages were grouped by taxonomic family: phage host range was significantly affected by the bean field of origin (F 3,172 = 14.2, p = 0.001), family (F 3,172 = 8.4, p = 0.001), and the interplay between them (F 8,172 = 3.9, p = 0.001) (Fig. 4B).
Based on Bray-Curtis dissimilarity <20%, 139 phages (70% of the total) were clustered into 24 PPGs; the remaining phages (57 phages) exhibited a unique PPG. Within the PPGs, two profiles accounted for 53% of the phages, whereas 22 profiles exhibited fewer than ten phages from similar PPGs. The PPG composition across common bean fields was significantly correlated with the RPG composition based on Bray-Curtis distances, but not based on Jaccard distances (Table S8). The PPG composition was also significantly correlated with PGT composition (Table S8). Accordingly, host range similarity was considerably correlated with phage ANI (r = 0.29, p = 0.001).
Although the Mexican phages were locally adapted overall, phage cross-infection was detected at similar rates between the communities of Tepoztlán and Yautepec (Mexico), indicating nonadaptive phage differentiation between these communities (Fig. 5E). Indeed, the infection rates of Mexican phages in Argentinian rhizobia were significantly lower, suggesting an effect of geographic distance on local adaptation ( Fig. 5B; Fig. 5E). On average, the rhizobia were more susceptible to sympatric phages (mean 0.54 ± 0.01; CV = 36.9%; σ 2 = 0.039) than to allopatric phages (mean 0.40 ± 0.01; CV = 49.7%; σ 2 = 0.040), indicating that the bacterial populations were maladapted to the local phages (F 1,456 = 46.7, p = 2.72 e-11; Fig. 5C). Similar results were obtained regardless of whether susceptibility rates were calculated for each Rhizobium isolate using only infections by phages isolated by using the SC (SC phages; F 1,456 = 38.9, p = 9.99e-10) or phages isolated by using local rhizobia (LC phages; F 1,456 = 45.9, p = 3.85e-11). All of these results suggest that rhizobia are maladapted to coexisting phages. However, the high susceptibility of rhizobia to local phages appeared to hold only for Argentinian Rhizobium communities (0.57 ± 0.02 versus 0.30 ± 0.02; F 1,212 = 64,3, p = 7.20e-14; SC: F 1,212 = 102.5, p < 2.2e-16; LC: F 1,212 = 42.3, p = 5.64e-10), since the Mexican Rhizobium isolates showed no significant differences in susceptibility to sympatric and allopatric phage infection (0.51 ± 0.02 versus 0.49 ± 0.01; F 1,242 = 1.2, p = 0.268; SC: F 1,242 = 0.4, p = 0.508; LC: F 1,242 = 2.0, p = 0.159) (Fig. 5D; Fig. 5F). When we considered the susceptibility of Rhizobium to sympatric phages or allopatric phages by site, the Chicoana Rhizobium community was found to be very susceptible to its own phages but resistant to phages from the other sites (Fig. 5F). In contrast, the rhizobial strains from Tepoztlán, Yautepec, and Salta were very susceptible to Chicoana phages (Fig. 5F), which partly explains the lack of significant differences in susceptibility to sympatric and allopatric phage infection by Mexican Rhizobium isolates.  In most cases, the difference between the mean sympatric rates of phage infection (S) and the mean allopatric rates of infection (A) indicated phage local adaptation in the four rhizobium-phage communities (Fig. 5G). In contrast, resident rhizobia were more susceptible to sympatric than to allopatric phages, except for the rhizobia from Yautepec that were also efficiently infected by Tepoztlán phages (Fig. 5H). The four populations of rhizobia were susceptible to the Chicoana phages, but the Chicoana rhizobia were mostly resistant to allopatric phages (Fig. 5H).

Discussion
Overall, our data indicate that the genetic and phenotypic diversity of phages and their Rhizobium hosts is spatially structured and that phages are adapted to their local host communities. Previous research has shown that rhizobia are spatially structured and can locally adapt to their legume hosts and other local environmental factors [49,50,88,89]. For instance, across the areas of common bean domestication, rhizobia receive a greater competitive benefit when nodulating sympatric common beans (Phaseolus vulgaris) than nonnative common bean varieties [44]. In turn, our results show that sympatric phage communities have locally adapted to these rhizobia.
We found that each of the four phage communities analyzed was dominated by a particular taxonomic family, prominently Microviridae in Tepoztlán, Myoviridae in Salta, and Siphoviridae in Chicoana. Moreover, phage genome sequencing revealed high genetic diversity within taxonomic families. Phage genetic diversity varied considerably within and between communities, and phages were clustered into PGTs. Approximately 52% of the 29 PGTs occurred solely in a single phage community (PGT). Similarly, Scola et al. [90] found that 66.4% of Namib Desert soil phage OTUs were exclusive to a single sampling site. Other PGTs (31%) occurred in both Mexican and Argentinian bean fields, with an average ANI of 88% across regions. Moreover, phages of the same PGT exhibited an 8% higher ANI within regions than between regions on average. Highly similar phages have previously been found across multiple distant aquatic ecosystems around the world, in human virome samples [91,92] and even in more distinct ecosystems [11]. The results indicate an emerging pattern in which a higher fraction of PGT members present a limited geographic range, while a significant minority of relatively closely related phages are distributed globally [16,90,93]. It remains unclear to what extent such observations are due to dispersal limitation [8]. The Rhizobium communities also showed spatial structure. R. etli was the predominant species nodulating common bean in the analyzed agricultural fields in Mexico, and R. phaseoli predominated in the Argentinian fields, although both species were mainly characterized by spatially restricted STs.
Rhizobia from Tepoztlan were isolated from nodules collected directly from the bean field, while the rhizobia from the other bean fields were obtained by soil inoculation of growing beans in the greenhouse (see methods). This may have contributed to the differences observed in the composition of the rhizobia-phage communities. However, the latter is a common method for trapping rhizobia in rhizobium research to accurately sample rhizobium diversity [94,95]. The specific technique used here was a slight variation of the method used in Van Cauwenberghe et al. [96]. Briefly, they found that rhizobial diversity and composition, was mostly identical to what was found using nodules collected directly from the same sampling sites in previous studies [97,98]. Although the time of sampling was during the spring of 2016, the Yautepec samples were obtained 3 months later than those from Tepoztlan. This may have caused a confounding effect regarding to the compositional differences between these two communities. However, our analyses found only a significant difference in presence-absence (based on Jaccard dissimilarity) of phages between these two fields, which is probably also affected by various environmental differences (e.g., a 430 m altitude difference) and local historical contingencies. Moreover, we did not found significant differences in phage relative abundance (Bray-Curtis) and rhizobium presence-absence and relative abundance between the Rhizobium communities of Tepoztlan and Yautepec. Since chloroform was used during the phage isolation process, it may have excluded the membrane-containing phages or the recently described Autolykiviridae [82]. Although the results presented here suggest a minor role of these potential sources of variation, further specific experiments addressing them should be done. In particular, to extend the model to chloroform sensitive phages, that were not considered in this work.
Phage community (PGT) spatial patterns were correlated with the compositional differences among Rhizobium (ST) communities. Similar correlations between host-phage communities have been seen in aquatic systems [16,17,22,[99][100][101], which seems to verify the common assumption that the relative abundance of phages within a community depends largely on host abundance and susceptibility [8,23]. However, the presence of susceptible rhizobium lineages did not imply the presence of specific phages or vice versa. For example, the omnipresent ST-5 was susceptible to all but two members of the phage lineages, but most phage lineages were spatially limited. Similarly, F06 phages could infect members of all STs, yet their presence was limited to Mexican bean fields. Although our study provides detailed genetic information on phages, the Rhizobium lineages were broadly defined by housekeeping gene markers (recA, dnaB). It is expected that much of the still-uncharacterized phenotypic and genetic microdiversity within bacterial species [102] may explain the spatial heterogeneity of the PGT composition and host range patterns better than the presence or absence of a suitable host (defined by species or ST).
Host range breadth varied considerably among phages, from generalists to specialists, resulting in a nested structure of the inferred whole cross-interaction network. Based on the extensive analysis of experimental cross-infection data, Flores et al. [103] concluded that nestedness is the characteristic profile in most cases. A modular network structure may be significant at large phylogenetic scales [87], but genotype-to-genotype interactions are most frequent within narrow phylogenetic ranges and result from coevolutionary processes of susceptibility and resistance. High host genetic similarity may underlie the nested structure of our network.
Nevertheless, we found that phage-rhizobium interactions were significantly affected by the genetic identity of both phages and rhizobia as well as their geographic origin. This may explain the detection of three large modules with high internal connectivity and suggests ongoing local adaptation. Indeed, we showed that phages infecting common bean-nodulating rhizobia experienced higher infection rates in sympatric rhizobia than in allopatric rhizobia; hence, they were locally adapted. Furthermore, sympatric phages showed more similar host ranges than allopatric phages, and sympatric rhizobia shared similar susceptibility ranges. Only a few field studies have provided evidence that phages locally adapt to their bacterial hosts in nature [27,87,104]. Although the Rhizobium communities were generally maladapted to the local phages, they may be adapted to their local environment as a whole. A nodule may provide an isolated niche for rhizobia where they may survive competition with and the antagonistic effects of other bacteria or, more directly, phages. In free-living conditions, depredation by phages may change the population structure of rhizobia. Phages are usually considered to be slightly ahead of their host in their coevolutionary arms race due to the higher selective pressure they experience and their greater evolutionary potential [27,105]. Although local bacterial adaptation to phages has been described multiple times in coevolutionary in vitro experiments [36,106,107], evidence in nature appears to be lacking [26]. This discrepancy is probably due to the relatively high availability of resources to hosts in vitro, which sways the arms race to the benefit of the host [108].
In our model, the degree of local adaptation was spatially inconsistent. Argentinean phages (mainly from Chicoana) infected approximately as many local as nonlocal rhizobia, while Mexican phages were more infectious in local rhizobia than in nonlocal rhizobia. Spatial asymmetry in phage local adaptation is believed to be the result of the effects of nutrients on phage-host encounter rates, mutation rates and the cost of resistance [38,109] or the local mode of coevolutionary dynamics (i.e., arms race or fluctuating selection [110]). Although we did not detect local phage maladaptation and we assume that the differences in productivity across the sampled actively cultivated bean fields might be relatively minor compared to those in media based in vitro experiments, these studies show how environmental differences create spatially different intensities of phage local adaption. Indeed, the spatial heterogeneity of environmental factors results in a geographic mosaic of different evolutionary pressures [111]. Local adaptation to various environmental conditions can undermine the colonization success of allopatric individuals and limit gene flow (i.e., "isolation-by-adaptation," or more general "isolation-by-environment" [33,35,109,111]). Zhang & Buckling [34] found that host bacteria grown in the presence of phages in heterogeneous environments were limited in their ability to migrate across environments as a result of maladaptation. The limiting effect of local adaptation on phage migration has not been tested explicitly, although it has long been predicted [24,109].
The spatial structure in the genetic composition of phage communities is probably due to the interplay of a variety of factors (e.g., historical contingencies, abiotic selection, genetic drift, and, potentially, dispersal limitation [8,15]. Our results indicate that the presence of suitable hosts may play a role in shaping phage biogeography and that suitability is determined not only by the genetic identity of the host but also by local adaptation. The spatial patterns are analogous to those observed in Rhizobium-common bean interactions and suggest that the local adaptation of rhizobia to common bean may have shaped the spatial differences in the phage-rhizobium interactions. Through isolation-byadaptation, local adaptation may reinforce spatial patterns in the PGT composition. Strong local adaptation of phages has been found across much shorter distances than in our present study [27,104], and it is as yet unclear to what extent phage local adaptation leads to limited migration and at which scale this may occur. At smaller scales, spatial heterogeneity is probably under greater pressure due to higher viral migrant densities. However, across broad scales, local adaptation may be a significant barrier to successful longdistance phage migration. reviewers whose suggestions helped to improve this paper. Special thanks go to Joshua Weitz for his advice on the BiMat application.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons. org/licenses/by/4.0/.