Acid-tolerant (pH 5.2) rhizobia strains isolated from nodules of Medicago sativa grown in acid soils in Uruguay ( Table 1 ), were obtained from the rhizobia collection of the laboratory of Biochemistry at the School of Agronomy, Universidad de la República, Uruguay. Or191 from Oregon, USA (Eardly and David, 1985) and LPU83 from Argentina (Torres Tejerizo et al., 2016) were used as reference strains of Rhizobium favelukesii. Ensifer meliloti U143 (=MCH3) strain, the inoculant for alfalfa in Uruguay (Altier et al., 2013) was supplied by Microbiology Laboratory (Ministerio de Ganadería Agricultura y Pesca). Rhizobia and E. meliloti were grown at 28°C in YEM (Vincent, 1970). Escherichia coli strain S17-1 ʎpir containing the plasmid pCAM131 with the transposon mTn5SSgusA31 (Wilson et al., 1995) grown at 37°C in LB (Miller, 1972) supplemented with Spectinomycin (Sp) 100 ug/mL and Streptomycin (St) 50 µg/mL. In liquid media the cultures were shaken at 120 rpm. Rhizobia and E. coli strains are stored at 4°C in YEM agar or LB agar Sp 100 µg/mL, respectively.

Medicago sativa cv. Chaná seeds were superficially sterilized according to (Irisarri et al., 2019). Plants grew in control condition at 23/20°C (day/night), with a 16/8 h photoperiod and 220 µE.m^-2.s^-1 photosynthetic photon flux density.

Time nodulation assay was performed in tubes with 20 mL Jensen medium (Jensen, 1942) N-free. Competitiveness, symbiotic efficiency and relative symbiotic effectiveness assays were performed in 350 mL pots with sterile vermiculite - sand (1:1), watered alternatively with distilled water and with Fahraeus medium (Vincent, 1970). Distilled water pH and Fahraeus medium were adjusted to 5.5 with MES 10 mM. Inoculation was carried out as described below in each assay.

The conjugation was carried out using the rhizobium strain ORY1 as receptor and E. coli strain S17-1 λ-pir containing pCAM131 that carries transposon mTn5SSgusA31 as donor, generally following the procedure used by Irisarri et al. (2019). The mixture remained at 28°C overnight in plate with TY (Beringer, 1974). Transconjugants were selected in plate with YEM Sp 100 µg/mL and Nitrofurantoin (Ntf) 20 µg/mL to inhibit E. coli growth.

Time nodulation assay of three transconjugants respect to the parental strain was evaluated according to time of appearance of the first nodule and total nodules at 21 days. The pre-germinated seeds, sown at a rate of one per tube with Jensen medium and 8 repetitions per treatment, were inoculated with 100 µL of a bacterial suspension with OD₆₂₀ = 0.9 (10⁸ ufc/mL). The ORY1::gusA clone was selected for the competition assay.

Nodulated roots were washed with distilled water and with 0.1 M phosphate buffer pH=7, incubated 16 h at 37°C in darkness in solution type containing 1% SDS, EDTA 0.5 M pH 8 and 1 mM 5-bromo-4-cloro-3indolil-β-D-glucuronide (X-gluc) in 50 mM sodium phosphate buffer pH 7.5 (Wilson et al., 1995). Nodules occupied by ORY1::gusA were visually identified by blue staining, in contrast to those infected by non-tagged strains which remained unstained (Irisarri et al., 2019).

Competitiveness of U143 and ORY1::gusA strains was carried out in pots with sand-vermiculite (1:1), with 5 plants per pot, and 4 replicates per treatment. Three assays with a n = 4 were conducted and the 12 resulting replicates were pooled for the statistical analysis. Pots were rinsed with water or Fahraeus medium at pH 5.5, as indicated above. Strains grew in YEM medium up to an OD₆₂₀ = 0.9. Competitiveness experiments included ORY1::gusA and U143 in a 99:1 and 1:99 ratio (named treatment 99% and 1% respectively). Bacterial concentration was determined according to Riviezzi et al. (2020). Plants were harvested 21 days after inoculation and washed root nodules were stained as described above. To assess symbiotic efficiency, biomass production was determined as shoot dry weight at 60°C for 72 h, including two controls: U143 strain and uninoculated treatment.

Seeds were sown at the rate of 5 per pots and inoculated with ORY1 and U143 strains as described previously, and harvested 35 days after inoculation. Each treatment included 5 replicates. Relative symbiotic effectiveness (%) is based on ORY1 parasitic strain shoot dry weight production relative to effective U143 strain.

Genomic DNA extraction from different rhizobial strains were extracted with DNA Extraction kit (Qiagen, Alemania). DNA quality and concentration were determined with Nanodrop 2000 (Thermo Scientific, USA), and its integrity was visualized in agarose gel 1.2%.

16S rRNA partial gene was amplified using the primers 27f and 1525r and conditions described by Lane (1991). PCR products were run at 100 V in Tris-acetate buffer pH 8.2 and visualized in 1.2% agarose gel stained with SYBR Safe DNA Gel Stain (Thermo Fisher Scientific, USA). Accuruler 1 kb DNA Ladder (Maestrogen Inc, Taiwan) was used as molecular weight marker. Sequencing was performed in Macrogen Inc. (Korea).

Total genomic DNA was amplified with BOXA1 primer (Koeuth et al., 1995). Amplification condition used was 1 cycle of 2 min at 95°C, 30 cycles of 45 seg at 93°C, 1 min at 50°C and 6 min at 65°C, and a final extension of 8 min at 68°C. PCR products were visualized as described previously.

Genomic sequencing was obtained using Illumina Technology (Macrogen, Korea) yielding 11.1 million paired-end reads with a length of 151 bp. Low quality reads were filtered out using Trimmomatic software with parameters 2:30:10 LEADING:3TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:75 (Bolger et al., 2014). SPAdes software was used to perform a de novo assembly using the remaining paired reads (Bankevich et al., 2012). The quality of the assembly was assessed with QUAST software using LPU83 genome as reference (Gurevich et al., 2013). Plasmid contigs were detected using plaSquid software (Giménez et al., 2022) with the option –minidist and using LPU83 plasmids as reference database. Genomic annotation of all the analysed genomes was done using prokka software with default parameters (Seemann, 2014).

The calculation of the Average Nucleotide Index (ANI), correlation indexes of tetra-nucleotide signatures, and tetra correlation search (TCS) were performed using the JSpeciesWSserver (Richter et al., 2016). The genome sequence data were uploaded to the Type Genome Server (TYGS), for a whole genome-based taxonomic analysis (Meier-Kolthoff and Göker, 2019).

Symbiotic gene nifH, nifD, nifK, nodA, nodB and nodD, sequences were extracted from the ORY1 genome. The obtained sequences were compared with reference strains sequences available in public National Center for Biotechnology Information (NCBI) database. Sequences were aligned with the ClustalW algorithm and Maximum likelihood phylogenomic tree (Tamura and Nei, 1993) were constructed for each multilocus sequence analysis (MLSA). Alignments and trees were constructed with MEGA 11 software (Tamura et al., 2021). Statistical support for tree nodes was evaluated by bootstrap analyses using 1000 replicates (Felsenstein, 1985).

The GenBank accession numbers for 16S rRNA sequences used in phylogenetic analysis are: Rhizobium favelukesii ORY1 (OP294988), ORY2 (OP294989), ORY3 (OP294990), ORY4 (OP294991), ORY5 (OP294992), Or191 (EU928874.1), Rhizobium anhuiense PVPR1(MT476932.1), Rhizobium indigoferae NAC94(MK872365.1), Bradyrhizobium yuanmingense VAFW14 (LC585439.1), Rhizobium croatiense 9T (MK753104.1), Rhizobium altiplani BR 10423 (NR 152084.1), Rhizobium tibeticum CCBAU85039 (NR116254.1) and Sinorhizobium meliloti 2011 (CP004140.1).

The GenBank accession numbers for genomes and plasmids from which the symbiotic genes and some partial 16S sequences were extracted are indicated in brackets: Rhizobium favelukesii ORY1 (JAIRAY000000000.1); LPU83 (HG916852.1) and Or191 (GCA_000419725.1); Rhizobium tibeticum CCBAU85039 (GCF_900108425.1); Rhizobium altiplani BR 10423 (GCA001542405.1); Rhizobium grahamii BG7 (CP043498.1); and Sinorhizobium meliloti 2011 (CP004140.1); Rhizobium tropici CIAT 899 (NC 020061.1); Rhizobium esperanzae N561 (CP013501.1); Rhizobium phaseoli R650 (CP013533.1); Rhizobium gallicum sp. R602 (CP006878.1); Rhizobium etli NXC12 (CP020907.1).

Ensifer meliloti USDA1963 peptidase HrrP was retrieved from NCBI reference proteins database and used to detect its domains in Pfam web server (Bateman et al., 2004). Hidden Markov models from M16 Peptidase (PF00675) and M16_C Peptidase (PF05193) were downloaded from Pfam database. To understand the presence of potential peptidases in plasmids obtained from rhizobium genus, 461 plasmids were downloaded from PLSDB database (V. 2021_06_23_v2) (Galata et al., 2019). All encoded proteins were predicted using prokka software. HMMER 3.0 software was used to scan for M16 peptidase domains using the option –cut_ga as filtering threshold (Eddy, 2011). All predicted proteins were clustered to 98% of identity, using CD-HIT software, to reduce sequence redundancy (Fu et al., 2012). R. favelukesii predicted M16 peptidases were added to the phylogenetic analysis. Sequences (available at: https://doi.org/10.6084/m9.figshare.21641816.v1) were aligned using muscle algorithm through msa R package (Bodenhofer et al., 2015). Phylogenetic trees were computed using neighbor-joining algorithm with ape R package (Paradis et al., 2004) and drawn using ggtree (Yu et al., 2017). Domain organization of each sequence was manually obtained from Pfam database search and added to the analysis.

ANOVA was performed in R software to identify statistical significance and a post hoc test (Tukey´s test) was used for pairwise comparison to determine the P-values.

To characterize the M. sativa-parasite rhizobia ORY1, ORY2, ORY3, ORY4, and ORY5 strains, isolated from Uruguayan acidic soils, we analysed BOX profile, the 16S rRNA sequences, and Relative symbiotic effectiveness.

BOX profiles for the five autochthonous isolates and the Argentinian and USA R. favelukesii strains showed that ORY1 and ORY5 strains had the same profile as the foreign strains, whereas ORY2, ORY3, and ORY4 strains had the same profile, but different to the foreign strains ( Figure 1A ).

The partial 16S rRNA sequence gene of ORY1, ORY4 and ORY5 were 100% identical to the R. favelukesii strains, LPU83 and Or191, isolated from Argentinian and USA acid soils respectively. ORY2 and ORY3 strains showed respectively a 99.9% and 99.7% identity to ORY1, LPU83, and Or191. Next, we performed a phylogenetic analysis of the different strains isolated from Uruguayan soils (ORY) among other related species. All ORY strains clustered together with R. favelukesii LPU83 and Or191 strains, and closely related to R. tibeticum ( Figure 1B ). Hereafter, ORY1 strain was used as the strain of study for all further assays. The parasitic behaviour of ORY1 was confirmed by shoot dry weight production in presence of this strain as the sole inoculant. This parameter was only 15% of such obtained with the commercial inoculant strain U143, evidencing a very low symbiotic efficiency.

The analysis of the draft genome of R. favelukesii ORY1 strain revealed that this genome is 7.44 Mb in size; has 59.7% GC content; containing 7,509 coding and 51 RNA genes. Five replicons were identified corresponding to the chromosome and four plasmids. The analysis of plasmids revealed that 151 contigs covered 96% of the total length of Rhizobium favelukesii LPU83 plasmids. This represents a total length of 3.25 Mb which encode 3,349 ORFs. Based on TCS or tetra-nucleotide signature, only R. favelukesii and R. tibeticum were found closely related with z-scores > 0.999. ANI values of > 96% consistently group genomes from strains of R. tibeticum and R. favelukesii together with ORY1, but ANI values of > 99% identify ORY1 as R. favelukesii. According to TYGS analysis ORY1 belongs to R. favelukesii. The phylogenetic analysis of R. favelukesii ORY1, Or191, and LPU83 strains showed that they clustered together, separated from R. tibeticum, when using the complete genomes ( Figure 2A ). In addition, R. favelukesii and R. tibeticum also grouped separated from R. grahamii and R. altiplani ( Figure 2A ). Likewise, R. favelukesii strains does not group with R. tibeticum when a MLSA of symbiotic genes, nifD, nifK, nifH, nodA, nodB, and nodD, was considered. In addition, these concatenated genes are related to E. meliloti symbiotic genes ( Figure 2B ). To better understand the similitude between the three R. favelukesii strains, ORY1, Or191, and LPU83, we analyse the overlap between the three complete genomes by identifying the homologous genes. We observed that most of the coding genes are common to all strains (around 87.5%, relative to ORY1), Or191 had a greater overlap with ORY1 than LPU83, and around 2.7% coding genes were unique for ORY1 ( Figure 2C ).

To evaluate the competitiveness of ORY1 strain we tagged ORY1 with the reporter gene, gusA. To confirm that the gusA insertion did not affect competitiveness and shoot biomass production, we compared the performance of ORY1::gusA strain for these parameters, to the parental strain ORY1. No differences were observed between these strains in terms of time to observe the first nodule (8 days), number of nodules per plant at 21 days, and shoot biomass.

Two treatments were used to test competitiveness and biomass production, consisting of mixes of ORY1 and the commercial inoculant, E. meliloti U143. One treatment consisted of a 99:1 ratio (ORY1::gusA/U143, hereafter 99% parasite strain treatment), whereas the other treatment was a 1:99 ratio (ORY1::gusA/U143, hereafter 1% parasite strain treatment). Nodule occupation by ORY1::gusA strain in 99% treatment was greater than in 1% treatment ( Figure 3A ). It is worth noting that even in 1% treatment, the percentage of occupied nodule exceeded the 50% of total nodules, evidencing the great competitiveness of this strain. Despite of the fairly similar occupation of both treatments, a remarkable difference was observed in terms of biomass production between 99% and 1% treatments ( Figure 3B ). In fact, the biomass production in 99% treatment was very low and not statistically different to the one obtained in the non-inoculated control, whereas the biomass production in the 1% treatment was not statistically different to the one obtained with U143 strain ( Figure 3B ). Regardless the percentage of parasite occupied nodules, the total number of nodules was much greater in the 99% treatment ( Figure 3C ). Despite of the greater number of nodules in the 99% treatment, the number of white nodules was significantly greater in the 1% treatment ( Figure 3D ).

Nodule morphology induced by ORY1::gusA and U143 strains was different ( Figure 3E ). We observed that ORY1::gusA induces the development of a great number small spherical nodules unlike U143, which induces a small number of larger undetermined nodules ( Figure 3F ). Exceptionally, a few U143-occupied-like nodules were occupied by ORY1::gusA ( Figure 3G ), a phenomena that we attribute to co-occupation of these strains, given that we never observed this morphology when plants are solely inoculated with ORY1::gusA.

Given that the bacterium genes hrrP and sapA encode for M16 peptidases that have been related to rhizobia with low symbiotic efficiency, we sought for genes coding M16 peptidase in the ORY1 genome. Bioinformatic analysis of protein domains using HMMER software, together with Pfam proteins database allowed us to evidence in ORY1 genome the presence of three genes encoding aminoacidic sequences for M16 zinc metallopeptidase domains. Two of these genes are contiguous gene located on an accessory plasmid, and the other is a chromosomic gene.

We performed a NCBI protein BLAST search and found that the 448 (ID 07295) and 532 (ID 07296) aminoacidic contiguous sequences encoded on the accessory plasmid of ORY1 strain showed an identity of 22.78 and 26.55%, and a coverage of 50 and 63%, respectively, to the HrrP peptidase of E. meliloti USDA1963 ( Table 2 ). This protein is an homodimeric M16 zinc metallopeptidase involved in the cleavage of NCR peptides. In the chromosome of ORY1 strain, we identified a DNA sequence coding a 432 aminoacid sequence presenting an 82.78% of identity and 98% of coverage to SapA peptidase from E. meliloti 1021 strain (CAC45492.2). This E. meliloti SapA is a member of a conserved M16B subfamily peptidase present in many different rhizobia.

To assess the presence of M16 genes in rhizobia, we analyzed 461plasmids available at PLSDB database. Only 30 genes encoding this type of protein domains were found, and a few of them corresponded to the same replicons. M16 domain-encoding genes were found in rhizobia plasmids of different species such as R. leguminosarum, R. etli, R. favleukesii and R. gallicum ( Figure 4 ). Additionally, clusters of sequences found in this phylogenetic tree do not show a species-based distribution.

To understand remote homology patterns in these genes we further analyzed domain architecture of the encoded peptidases in rhizobium plasmids ( Figure 4 ). We observed that there are different domain organizations but with a limited number of protein domains; specificallyM16 (PF00675) and M16_C (PF05193). Peptidase M16 followed by Peptidase M16_C is the domain architecture most present in these proteins, at least one gene for the 5 rhizobia species represented in the tree have this domain organization. A reference of an uncharacterized protein having this domain is deposited in UniProt under accession number W2S8M8 with its predicted AlphaFold structure. As a functionally characterized protein we included HrrP of E.meliloti USDA 1963, which resulted in a similar domain architecture only having an extra Peptidase M16_C domain. A reference protein having these domains is also deposited in UniProt under accession W9RIX3.

Alfalfa parasitic Oregon-like rhizobia are a problem for the implantation and production grasslands of this legume (Torres Tejerizo et al., 2011). These rhizobia had been identified in USA, Canada (Eardly et al., 1985; Bromfield et al., 2010) and in Argentina, where the novel species R. favelukesii was described (Torres Tejerizo et al., 2016). Castro-Sowinsky et al. (2002) purposed the presence of Oregon-like rhizobia in Uruguayan soils. However, the presence of this rhizobium in Uruguayan soils had not been confirmed before. In this study, based on the whole genome sequencing analysis of five isolates obtained from Uruguayan acid soils, where there are yield problems in alfalfa production, we confirmed the presence of Oregon strains ( Figure 2A ).

R. favelukesii and R. tibeticum can nodulate M. sativa (Hou et al., 2009), but R. tibeticum has multiple hosts such as Trigonella archiducis-nicolai, Medicago lupulina, Phaseolous vulgaris, among others (Abd-Alla et al., 2014; Torres Tejerizo et al., 2017). In congruence, our genome based phylogenetic analysis showed that R. favelukesii strains (ORY1, Or191, and LPU83) clustered together among them but separated from R. tibeticum strain CCBAU85039 ( Figure 2A ). Moreover, the genome comparison based on common genes revealed that ORY1 strain is more similar to Or191 isolated from USA than LPU83 isolated from Argentina ( Figure 2C ). The limited number of unique genes in ORY1 (2.7%) are mostly uncharacterized genes and we did not identify any symbiotic genes.

Our results also showed a greater diversity than that observed in Argentinian acid soils. In our study, among the five strains isolated from Uruguayan acid soils, two strains presented the same BOX profile as LPU83 and Or191, identified in Argentina and USA respectively, while the other three presented a profile different to those. Unlike this result, twelve isolated from Argentinian acid soils resulted in the very same ERIC fingerprint profile as Or191 from USA (Wegener et al., 2001), representing a very low genetic diversity. The low genetic diversity is different to what is often reported for rhizobia having a wide geographic distribution (Martínez-Romero and Caballero-Mellado, 1996), and to what we observed for clover-nodulating rhizobia (Irisarri et al., Tartaglia et al., 2019) and lotus-nodulating rhizobia (Sotelo et al., 2011; Batista et al., 2013) isolated from Uruguayan soils. Nevertheless, the twelve isolates tested by Wegener et al. (2001) were selected from an initial collection of 465 isolated from Argentinian soils, with pH between 5.5 and 6.2, by growing them at pH 5.0 (Del Papa et al., 1999). This means that the twelve isolated with identical ERIC fingerprint were selected by high tolerance to low pH, and extreme environmental conditions are known to reduce genetic diversity, whereas the pH of Uruguayan soils where our strains were isolated were only slightly acid (pH 5.4 - 5.8), explaining the greater genetic diversity observed in our collection.

The high competitiveness of the parasitic ORY1 strain was evidenced by the occupancy of nodules in the 1% treatment ( Figure 3A ). Interestingly, the biomass production was high, and not different to that for U143 ( Figure 3B ), even though more than 50% of nodules were occupied by the parasite strain ORY1.This result suggests that the presence of parasite rhizobia should not be a problem for biomass production if our design (1:99, parasite:efficient rhizobia) represents the situation occurring during the first few days after sowing, where the strain used as commercial inoculants are expected to be at greater concentration than native parasite rhizobia. On the other hand, our results show that if the load of commercial inoculant is very low, or if the load of parasite rhizobia in the soil is high, severe problems in biomass production can be found ( Figure 3B ).

At first glance, it sounds counterintuitive the fact that we observe a high biomass production when the occupation by the parasite strain, ORY1, is greater than 50%. This can be in part explained by the fact that the 1% treatment had approximately 30% more white (effective) nodules relative to the 99% treatment ( Figure 3D ). It is worth noting that the incompatible plasmid-carrying strain can induce a higher number of nodules than the effective strains, even though the resulting nodules are less beneficial to the host (Crook et al., 2012). This goes in line with our observations, where we found a greater number of total nodules in the treatment enriched in ORY1 relative to the 1% treatment ( Figure 3C ). However, the nodules produced by the parasite strain were shown to be less beneficial for the plant compared to those produced by effective strains, which were suggested to produce bigger and more dense nodules (Crook et al., 2012). Therefore, if we only consider the white nodule as nitrogen fixing active nodules, it is not surprising that the 1% treatment had a greater biomass production than the 99% treatment. Another factor that could account for the high biomass production even when the occupation of ORY1 is relatively high, is the possibility of co-occupation with efficient strains. Co-occupation has been reported in M. sativa – E. meliloti symbiosis (Checcucci et al., 2016). Hence, we hypothesize that many of the blue nodules can be also occupied with the effective strain which contributes to nitrogen fixation. In fact, we observed a few well-developed nodules that were marked with the parasite strain ( Figure 3G ), strongly suggesting the possibility of co-occupation in those nodules. Thus, our observations remark the importance of not only evaluating competitiveness as percentage of occupied nodules, as usually tested (Irisarri et al., 2019), but also in terms of total, well-developed, not developed, and co-occupied nodules, to better understand the consequence at plant level.

One of the plasmid detected in ORY1 genome encoded contagious putative M16 peptidase, which are highly conserved in the three available genomes of R. favelukesii. The phylogenetic analysis of putative M16 peptidases in a comprehensive set of Rhizobium plasmids showed that the R. favelukesii putatives peptidases do not cluster by species in the phylogenetic tree, similarly to what occurs for Ensifer sequences ( Figure 4 ). This has been proposed as a result of lateral transfer events shaping the evolutionary pathways of these plasmid peptidases (Price et al., 2015).

M16 plasmid peptidases have been previously studied in E. meliloti and could be linked to host restriction functions (Price et al., 2015). The two M16 peptidases encoded in a R. favelukesii plasmid are contiguous and both of them together mimic the domain architecture of E. meliloti HrrP protein, where one sequence encodes for M16 - M16C domains, and the next sequence for another M16C domain. Given the proximity of M16 peptidase-encoding genes we propose that they could be related to hrrP as part of a gene fusion or fission event, as it has been seen in a set of Gram + and Gram – bacterias (Pasek et al., 2006). Additionally, these proteins could be implied in a similar host restriction mechanism as the one previously described in E. meliloti species. In particular, HrrP and SapA peptidases play a role in the digestion of NCR peptides produced by the host, for instance, the indispensable peptide for bacteroid differentiation NCR169 (Horváth et al., 2015; Price et al., 2015; Benedict et al., 2021). These peptides are necessary for the correct establishment of the symbiosis E. meliloti - M. truncatula (Price et al., 2015; Benedict et al., 2021). The plasmidic peptidase HrrP degrades NCR peptides, attributing the rhizobia an active mechanism to control de host (Pan and Wang, 2017; Benedict et al., 2021), and its expression is often associated with a parasite behaviour (Price et al., 2015). The presence of putative hrrP and sapA genes in R. favelukesii Oregon strains, ORY1, Or191, and LPU83 can partially explain the parasitic behaviour of R. favelukesii evidenced by the unusual globular aspect of nodules, hyper-nodulation, and low symbiotic efficiency.