Influence of flooding and soil properties on the genetic diversity and distribution of indigenous soybean-nodulating bradyrhizobia in the Philippines.

One of the strategies that is commonly used in the Philippines to improve the production of soybean is by inoculation. However, this technique often fails mainly due to the lack of information about the indigenous soybean rhizobia in the Philippines soil. In this study, the diversity of indigenous bradyrhizobia collected from the non-flooded and flooded soil conditions at 11 locations in the country was investigated using a local soybean cultivar as the host plant. The genetic variation among the 424 isolates was detected through Polymerase Chain Reaction-Restriction Fragment Length Polymorphism (PCR-RFLP) treatment and sequence analysis for 16S rRNA gene, 16S-23S rRNA internal transcribed spacer (ITS) region and rpoB housekeeping gene. All the isolates were classified under the Bradyrhizobium species namely B. elkanii, B. diazoefficiens, B. japonicum, B. yuanmingense and a considerable proportion of the isolates were clustered under Bradyrhizobium sp. The isolates which were classified under Bradyrhizobium sp. were thought to be endemic to Philippines soil as evidenced by their nucleotide divergence against the known rhizobia and the historical absence of rhizobia inoculation in the collection sites. The major influence on the distribution and diversity of soybean bradyrhizobia is attributed to the difference in the flooding period, followed by soil properties such as pH, soil type, and nutrient content. As determined, it is proposed that the major micro-symbiont of soybean in the Philippines are B. elkanii for non-flooded soils, then B. diazoefficiens and B. japonicum for flooded soils.

Introduction

The genus Bradyrhizobium is a major micro-symbiont of soybean (Glycine max [L.] Merrill) that could form nodule with legumes through symbiosis. The symbiotic relationship between rhizobia and legume is a highly specific, complicated and energy-exhaustive process (Wang et al., 2012) which may lead to the improvement of crop's yield. It is considered that Bradyrhizobium is the predominant genus of rhizobia in the tropics (Delamuta et al., 2015). Tropical regions have diverse environmental gradients that could influence the diversity of organisms. In this regard, a high diversity of soybean-nodulating rhizobia were reported to exist in tropical regions (Delamuta et al., 2012) than temperate regions. By far, diverse species of bradyrhizobia were described as the micro-symbionts of soybean such as Bradyrhizobium japonicum, B. elkanii, B. yuanmingense, B. liaoningense, B. huanghuaihaiense and B. diazoefficiens (Jordan, 1982; Kuykendall et al., 1992; Xu et al., 1995; Saeki et al., 2006; Appunu et al., 2008; Zhang et al., 2012; Delamuta et al., 2013). Accordingly, many reports stated that soil pH, salinity, climate, nutrients, and cultural management have great influence on the diversity and distribution of soybean bradyrhizobia (Loureiro et al., 2007; Grossman et al., 2011; Adhikari et al., 2012; Shiro et al., 2013; Yan et al., 2014; Chibeba et al., 2017; Mason et al., 2017; Saeki et al., 2017). A high diversity of microorganisms in the soil plays an important role in maintaining soil health so, a vast and accurate information about this could lead to a better crop productivity.

The Philippines is a tropical country that is characterized by high temperature, humidity and abundant rainfall. It has two pronounced seasons based mainly from rainfall: (1) dry season – from December to May and (2) rainy season – from June to November. Aside from Baguio City, the average temperature in the country is almost similar at 25.5–27.5 °C throughout the year. The country's agricultural production system has been heavily dependent on chemical inputs with rice or corn mono-cropping that it has rendered the soil acidic and unproductive. Although rice is the major agricultural product in the Philippines, a recent trend recognized the role of soybean in nutrition and soil fertility restoration. This prompted its Government to create a Research and Development Roadmap for Soybean to increase the area and volume of production. This also paved the way for researchers to improve the inoculation techniques to support the government's endeavor.

Yet, inoculation does not always succeed due to several reasons such as the incompatibility between the macro and micro-symbiont (Yamakawa et al., 2003; Hayashi et al., 2012; Yamakawa and Saeki, 2013), competition between the indigenous and inoculated rhizobia (Schumpp and Deakin, 2010; Ji et al., 2017) and the agro-environmental factors stated above that influence the distribution and diversity of rhizobia. It is then deemed necessary that prior to inoculation, information on the distribution and diversity of indigenous rhizobia in the soil should be obtained.

Except for one location in the Philippines, which is Nueva Ecija (Mason et al., 2017), there is no information about the indigenous soybean bradyrhizobia in Philippine soils. Since the Philippines is composed of more than 7,000 islands surrounded by bodies of water with considerable variation in agro-environmental conditions, it is hypothesized that a diverse species of bradyrhizobia can be collected and identified. Therefore, this study aimed to (i) identify the indigenous soybean bradyrhizobia in the Philippines soil and (ii) determine the factors that influence its distribution and diversity. The output of this study could help to improve the current status of soybean production in the country and to better understand the biogeography of tropical bradyrhizobia.

Materials and methods

Soil collection

In this study, soil samples were collected as formerly described (Mason et al., 2017) from eleven locations which were previously and/or currently planted with soybean and/or other legumes (Fig. 1). The locations that were selected have two different periods of flooding condition (flooded means the soil usually contains gravitational water for a period of at least 10 months in one year while non-flooded means the soil usually contains gravitational water for a period of at most 5 months in one year). The location and basic information on the study sites are summarized in Table 1. From the one (1) kg composite soil, a 0.5 kg was air dried and pulverized for soil analyses that include pH (1:2.5 water extraction method) and electrical conductivity (EC: 1:5 water extraction method), Total C, Total N, Bray P, K (flame spectrophotometer) and soil texture then, the remaining 0.5 kg was freshly used for soybean cultivation. For soil texture analysis, the Hydrometer method was used as described (Bouyoucos, 1962), whereas C and N were analyzed by an automatic high sensitive NC Analyzer Sumigraph NC-220F (Sumika Chemical Analysis Service. Ltd., Tokyo, Japan). Data of the annual average temperature and rainfall were obtained from Philippine Atmospheric Geophysical and Astronomical Services Administration (PAGASA) website (https://www1.pagasa.dost.gov.ph/) which were averages from the last two decades.

Fig. 1

Geographic distribution of the soybean-nodulating bradyrhizobia in the Philippines from the result of Restriction Fragment Length Polymorphism (RFLP) and sequence analysis of the 16S-23S rRNA gene ITS region.

Table 1

Information on the agro-environmental gradients of the Philippines soil used in this study.

Location	Coordinate	pH	EC (dS/m)	C (%)	N (%)	C/N	Bray P (mg/kg)	K (mg/kg)	Sand (%)	Silt (%)	Clay (%)	Ave. Temp. (°C)	Flooding Period (month)
Ilagan, Isabela (IS)	17.30°N,122.01°E	5.90	0.08	1.34	0.13	10.3	1.86	51.80	28.0	34.5	37.5	26.9	5
Gamu, Isabela (GI)	17.08°N,121.79°E	5.52	0.15	1.85	0.17	10.9	2.30	58.60	28.2	33.8	38.0	27.0	10
Baguio, Benguet (BA)	16.40°N,120.60°E	5.22	0.20	3.10	0.24	12.9	22.22	51.00	19.0	41.4	39.6	19.3	10
Nueva Ecija1 (NE1)	15.74°N,120.93°E	6.21	0.05	1.37	0.13	10.5	6.74	73.90	28.7	34.6	36.7	26.9	5
Nueva Ecija2 (NE2)	15.74°N,120.93°E	5.81	0.12	2.36	0.22	10.7	21.63	49.40	27.4	34.7	37.9	26.9	10
Irosin, Sorsogon (SO)	12.72°N,124.04°E	5.26	0.15	1.92	0.22	8.7	2.57	55.80	28.9	33.6	37.5	27.1	10
Abuyog, Leyte (LT)	10.67°N,125.04°E	5.80	0.12	1.50	0.15	10.0	6.39	174.20	29.2	32.8	38.0	27.0	10
La Carlota, Negros Occidental (NR)	10.24°N,122.59°E	5.62	0.15	0.63	0.07	9.0	20.44	74.10	28.0	34.2	37.8	27.7	5
Ubay, Bohol (BO)	9.99°N,124.45°E	5.82	0.11	0.63	0.06	10.5	2.80	47.80	29.7	33.5	36.8	27.2	5
Sultan Kudarat, Maguindanao (SK)	6.51°N,124.42°E	6.64	0.14	2.48	0.19	13.1	4.53	59.60	24.0	34.1	41.9	27.3	10
Tupi, South Cotabato (SC)	6.34°N,124.97°E	5.52	0.15	1.36	0.14	9.7	31.18	47.20	19.5	42.5	38.0	25.4	10

Soybean cultivation and isolation of indigenous rhizobia

A commonly available soybean cultivar in the Philippines, with the local name PSB-SY2, was used to isolate the indigenous soybean rhizobia. Soybean seeds were surface-sterilized by soaking in 70% ethanol and sodium hypochlorite solution as formerly described (Saeki et al., 2006) and planted in 1-liter culture pots (n = 4). Next, the culture pots were filled with vermiculite then, N-free nutrient solution was added (Saeki et al., 2004) at 40% (vol/vol) distilled water content and were autoclaved for 20 min at 121 °C. A soil sample (2–3 g) was placed on the vermiculite at a depth of 2–3 cm, the surface-sterilized seeds were then sown on the soil, and the pot was weighed. Afterwards, the plants were grown for 28 days inside the growth chamber (33 °C for 16 h, day; 28 °C for 8h, night), and were supplied weekly with sterile distilled water until the initial weight of the pot was reached.

After 28 days, 15 to 20 nodules were randomly collected from the roots of each soybean plant and sterilized with ethanol and sodium hypochlorite solution as described previously (Suzuki et al., 2008). Each nodule was crushed in sterile distilled water inside a 1.5 mL microtube by using a sterilized toothpick and streaked on the yeast extract mannitol agar (YMA; Vincent, 1970) plate. It was then incubated for about 1 week in the dark at 28 °C. To determine the genus, a single colony was picked-up by a sterile wire loop and streaked onto YMA plate containing 0.002% (wt/wt) bromothymol blue (BTB; Keyser et al., 1982) and incubated as described above. Thereafter, the pure single colonies were obtained by repeated streaking into YMA plates.

All the isolates obtained from the 11 locations were subjected to an inoculation test to determine their nodulation capability. Each isolate was cultured in YM broth culture (YMB; Vincent, 1970) for about 1 week at 28 °C, and the cultures were diluted with sterile distilled water to approximately 106 cells ml−1. Then, surface-sterilized soybean seeds were sown as described previously but without soil, and a 1-ml aliquot of each isolate per seed was inoculated with three replications. The plants were grown inside the growth chamber for 4 weeks then, nodule formation was evaluated. A control pot (un-inoculated) was prepared under similar conditions.

DNA extraction

The pure single colony isolated from the YMA plate from repeated streaking was cultured in HEPES-MES (HM) broth culture (Cole and Elkan, 1973; Sameshima et al., 2003) for 3–4 days at 28 °C with continuous agitation at 120 rpm. Thereafter, the bacteria cells cultured in the HM broth were collected by centrifugation and washed with sterile distilled water. The extraction of DNA was done by using BL buffer as from the method reported previously (Hiraishi et al., 1995) which was described by Minami et al. (2009).

PCR-RFLP analysis

All the isolates that produced nodules on the soybean plants were used for PCR amplification. The amplification of the 16S rRNA, ITS region, and rpoB gene was conducted using Ex Taq DNA polymerase (TaKaRa Bio, Otsu, Shiga, Japan) with primers and PCR cycle conditions previously used (Mason et al., 2017). Then, the RFLP analyses of the 16S rRNA and ITS region were performed using the restriction enzymes HaeIII, HhaI, MspI and XspI (TaKaRa Bio) whereas for the rpoB gene, enzymes HaeIII, MspI and AluI (TaKaRa Bio) were used. For the reference strains, B. japonicum USDA 4, 6T, 38, 122, 123, 124, 129, 135, B. diazoefficiens USDA 110 T, B. elkanii USDA 31, 46, 76T, 94, and 130 and B. liaoningense USDA 3622T (Saeki et al., 2004) were used. From the PCR product, a 5.0 μl aliquot was digested with the restriction enzymes for 16 h at 37 °C in a 20 μl reaction mixture. A submerged gel electrophoresis using 3 or 4% agarose gels in TBE buffer was conducted to separate the fragments and were visualized with ethidium bromide.

Sequence analysis

According to the OTUs obtained from the RFLP analysis of the 16S rRNA, ITS region and rpoB gene, representative isolates were selected for sequence analysis. The amplified products from the PCR reaction of the representative isolates for the ITS region and rpoB gene were purified according to the protocol of NucleoSpin® Gel and PCR Clean-up (Macherey-Nagel, Germany). Then, the DNA concentration of the purified product was confirmed with the NanoDrop 2000 Spectrophotometer (Thermo Scientific, U.S.A.). After then, samples were prepared according to the protocol for a premixed template and primer of the manufacturer (EUROFINS GENOMICS) using the previously designed sequence primers (Mason et al., 2017) and were sent to the company (EUROFINS). Based from the similarity in the ITS region-rpoB gene type of the representative isolates, a few isolates were randomly selected to confirm the sequences of 16S rRNA gene.

Sequence alignment and construction of phylogenetic trees

The Basic Local Alignment Search Tool (BLAST) program in DNA Databank of Japan (DDBJ) was used to determine the nucleotide homology. The sequences of type strains having similarity with our isolates of at least 99% for 16S rRNA, 96% for ITS region and 98% for rpoB gene were retrieved from the BLAST database. The phylogenetic trees also included previously determined sequences of the 16S rRNA and ITS region Bradyrhizobium genospecies (Saeki et al., 2000; Van Berkum and Fuhrmann, 2000). Then, the alignment of sequences obtained were performed using ClustalW. The Neighbor-Joining (Saitou and Nei, 1987) method for the 16S rRNA, ITS region and rpoB gene was used to construct the phylogenetic trees. Genetic distances were calculated using the Kimura 2-parameter model (Kimura, 1980) in the Molecular Evolutionary Genetic Analysis (MEGA v7) software (Kumar et al., 2016). Thereafter, the phylogenetic trees were bootstrapped with 1,000 replications of each sequence to evaluate the reliability of the tree topology. All the nucleotide sequences determined in this study were deposited in DDBJ at http://www.ddbj.nig.ac.jp/.

Cluster and diversity analysis of indigenous bradyrhizobia

Only those isolates with reproducible fragments longer than 50 bp in the electrophoresis gels were used for the cluster analysis. The genetic distance between pairs of isolates (D) was calculated using the equation DAB = 1 [2NAB/(NA + NB)], where NAB represents the number of RFLP bands shared by strains A and B whereas NA and NB represent the numbers of RFLP bands found only in strains A and B, respectively (Nei and Li, 1979; Sakai et al., 1998). Then, diversity analysis was performed by the Shannon-Wiener diversity index as described previously (MacArthur, 1965; Pielou, 1969; Saeki et al., 2008).

Principal component analysis (PCA)

To detect the relationship between the agro-environmental factors and the distribution of soybean bradyrhizobia, PCA was performed in R software. The variables for the principal component include soil chemical properties, soil texture and some environmental data.

Results

Isolation of indigenous bradyrhizobia

A total of 771 isolates was obtained from the study sites with a range of 63–79 isolates per location. The samples were labeled with a combination of the abbreviation of the sampling site (IS - Ilagan; GI - Gamu; BA - Baguio; NE1 - 1st location in Nueva Ecija; NE2 - 2nd location in Nueva Ecija; SO - Sorsogon; LT - Leyte; BO - Bohol; NR - Negros; SK - Sultan Kudarat; and SC - South Cotabato) and the number of the isolate (1–63 or 1–79) (e.g., for South Cotabato, SC – 1–74).

The inoculation test showed that 424 isolates were capable to form nodules on the host plant and all these isolates has produced alkaline or neutral reaction in the YMA plates with BTB thus, were considered as Bradyrhizobium (Jordan, 1982). We also confirmed that the un-inoculated pots did not produce any nodules.

RFLP and sequence analysis of the 16S rRNA, ITS region and rpoB gene

The 424 isolates were all used for the RFLP treatment of the 16S rRNA, ITS region and rpoB gene. The differences in the fragment sizes and patterns through the RFLP treatment showed 31 OTUs which were used for sequence analysis and the results are summarized in Table 2. It is evident that the genetic variations mostly occurred for Be clusters within the strains of B. elkanii 31, 46 and 76.

Table 2

Summary of the genetic variations detected from the 31 representative indigenous soybean-nodulating bradyrhizobia in the Philippines through the Restriction Fragment Length Polymorphism (RFLP) treatment and sequence analysis of the 16S rRNA gene, 16S-23S rRNA ITS region, and rpoB housekeeping gene.

Rep. isolate	Isolate (no.)	Cluster (RFLP)			Cluster (sequence analysis)			ITS- rpoB type
		16S rRNA	ITS region	rpoB gene	16S rRNA	ITS region	rpoB gene
IS-2	40	Be	Be76	Be46	B.elkanii	Be31	Be46	Be31-Be46
GI-4	30	Bj	Bj	Bj	B. japonicum	Br	Br	Br-Br
GI-8	6	Be	Be76	Be46	nd	Be31	Be46	Be31-Be46
BA-24	31	Bj	Bj6	Bj6	B. japonicum	Bj6’	Bj6	Bj6’-Bj6
BA-41	1	Bj	Bj124	Bj124	B. japonicum	Bj124’	Bj124’	Bj124’-Bj124’
BA-42	1	Be	Be76	Be76	B.elkanii	Be31	Be46	Be31-Be46
NE1-6NE1-19	495	BeBr	Be76Br	Be46Br	B.elkanii nd	Be76’Br	Be46Br	Be76’-Be46Br-Br
NE1-34	1	Bj	Bj	Bj	B. japonicum	BrNE	BrNE	BrNE-BrNE
NE1-65	1	Bj	Bj	Bj	B. diazoefficiens	Bd110	Br	Bd110-Br
NE2-1	3	Be	Be76	Be46	B.elkanii	Be76	Be46	Be76-Be46
NE2-3	5	Bj	Br	Br	nd	Br	Br	Br-Br
NE2-37	26	Br	Br	Br	B. japonicum	Br	Br	Br-Br
NE2-66	1	Bj	Bj	Bj	nd	Br	Bj124’	Br-Bj124’
SO-1	44	Bj	Bd110	Bd110	B. diazoefficiens	Bd110	Bd110	Bd110-Bd110
LT-3	42	Bj	Bd110	Bd110	B. diazoefficiens	Bd110	Bd110	Bd110-Bd110
LT-36	1	Be	Be76	Be46	nd	Be76’	Be46	Be76’-Be46
NR-1	4	Bj	Bj6	Bj6	nd	Bj6	Bj6	Bj6-Bj6
NR-2	22	Be	Be76	Be46	nd	Be76	Be46	Be76-Be46
NR-40	3	Bj	Bj124	Br	B. japonicum	Bj124’	Bj	Bj124’-Bj
NR-48	1	Bj	Bd110	Br	B. diazoefficiens	Bd110	Bd110’	Bd110-Bd110’
NR-60	1	Br	Be76	Br	nd	Be76	Br	Be76-Br
BO-4	24	Be	Be76	Be46	B.elkanii	Be31	Be46	Be31-Be46
BO-15	4	Bj	Bd110	Bd110	nd	Bd110	Bd110	Bd110-Bd110
BO-52	1	Bj	Bj	Br	nd	Br	Br	Br-Br
SK-1	4	Be	Be46	Be76	B.elkanii	Be46’	Be130	Be46’-Be130
SK-2	6	Be	Be94	Be94	B.elkanii	Be94’	Be94	Be94’-Be94
SK-5	29	Bj	Bd110	Bd110	nd	Bd110	Bd110	Bd110-Bd110
SK-12	4	Br	Br	Br	B. yuanmingense	BySK	BySK	BySK-BySK
SC-3	31	Bj	Bj6	Bj6	B. japonicum	Bj6	Bj6	Bj6-Bj6
SC-49	3	Be	Be76	Be46	nd	Be76	Be46	Be76-Be46
Total isolates	424

Note: The prime (’) symbol indicates a variation in the sequence identity of at least 1.00% between the isolates and the Bradyrhizobium USDA reference strains. Be – B. elkanii, Bd – B. diazoefficiens, Bj – B. japonicum, By – B. yuanmingense, Br – Bradyrhizobium sp., nd – not determined due to the similarity with the other isolates that possess an identical ITS-rpoB type.

Presented in Fig. 2 is the phylogenetic tree from the sequence analysis of the ITS region showing the genetic diversity of soybean bradyrhizobia across the country. The isolates were predominantly grouped under four bradyrhizobia species which are B. elkanii (37.74%), B. diazoefficiens (29.48%), B. japonicum (16.51%), and B. yuanmingense (0.94%). A Bradyrhizobium sp. clade that makes up about 15.09% of the population was also observed including two independent single isolates NE1-34 and NE2-66. On the other hand, almost similar clusters were observed from the rpoB gene phylogenetic tree (Fig. 3) wherein eight delineations were also distinguished similarly with the ITS region. Notably, some genetic variations were detected only on specific isolates between the two gene loci. For example, SK-1 belonged to Be46 cluster in the ITS region but it was reclassified into Be130 cluster in the rpoB gene. This observation mainly occurred with the isolates clustered under B. elkanii strains except for those clustered with the Be94. On the other hand, no remarkable genetic diversity was observed for the other clusters particularly for those under B. diazoefficiens and B. yuanmingense.

Fig. 2

Phylogenetic tree based on the sequence analysis of the 16S-23S rRNA gene ITS region. The tree was constructed using the Neighbor-Joining method with the Kimura 2-parameter distance correlation model and 1000 bootstrap replications in MEGA v.7 software. The accession numbers are indicated only for sequences obtained from the BLAST database. The isolates in this study are indicated with letters and number combinations, for example: IS-2 – isolate no. 2 collected from Ilagan, Isabela.


Fig. 3

Phylogenetic tree based on the sequence analysis of the rpoB housekeeping gene. The tree was constructed using the Neighbor-Joining method with the Kimura 2-parameter distance correlation model and 1000 bootstrap replications in MEGA v.7 software. The accession numbers are indicated only for sequences obtained from BLAST database. The isolates in this study are indicated with letters and number combinations, for example: IS-2 – isolate no. 2 collected from Ilagan, Isabela.

Distribution of isolates and diversity analysis

The distribution of the isolates in the locations and the diversity indices of soybean bradyrhizobia in the Philippines are presented in Table 3 and are graphically seen in Fig. 1. Those isolates that were clustered under the B. elkanii strains were present in almost all the locations, suggesting a prevalence of this species in the regions with high temperature. The highest diversity and equitability indices were obtained from Sultan Kudarat (SK) (H′ = 0.98, Eh = 0.71), followed by Negros Occidental (NR) (H′ = 0.82, Eh = 0.59). On the other hand, the lowest indices (H′ = 0.00, Eh = 0.00) were obtained from Isabela (IS) and Sorsogon (SO) where all isolates belonged to B. elkanii and B. diazoefficiens, respectively.

Table 3

Population distribution and diversity indices of the indigenous bradyrhizobia in the Philippines. The Shannon's diversity (H′) and equitability (Eh) indices were computed with the formulae (H′ = −∑Pi ln Pi; Eh = H′/ln S).

Cluster/location	IS	GI	BA	NE1	NE2	SO	LT	NR	BO	SK	SC
Be31	40	6	1	0	0	0	0	0	24	0	0
Be46	0	0	0	0	0	0	0	0	0	4	0
Be76	0	0	0	49	3	0	1	23	0	0	3
Be94	0	0	0	0	0	0	0	0	0	6	0
Bj6	0	0	31	0	0	0	0	4	0	0	31
Bj124	0	0	1	0	0	0	0	3	0	0	0
Bd110	0	0	0	5	0	44	42	1	4	29	0
By	0	0	0	0	0	0	0	0	0	4	0
Bradyrhizobium sp.	0	30	0	1	31	0	0	0	1	0	0
Independent bradyrhizobia	0	0	0	1	1	0	0	0	0	0	0
Total	40	36	33	56	35	44	43	31	29	43	34
H′	0.00	0.45	0.27	0.48	0.42	0.00	0.11	0.82	0.55	0.98	0.30
Eh	0.00	0.65	0.25	0.34	0.38	0.00	0.16	0.59	0.50	0.71	0.43

Note: The Pi is the dominance of the isolate, expressed as (ni/N), where N and ni are the total number of isolates tested and the number of isolates belonging to a particular cluster, respectively. S is the total number of clusters, indicating the taxonomic group, at each field site.

PCA of the factors affecting the diversity and distribution of bradyrhizobia

Presented in Fig. 4 is the PCA plot showing the factors that has influence on the diversity and distribution of indigenous bradyrhizobia in the country. PC1 shows 47.74% proportion that is accounted for most of the variance and indicating the correlation between the parameters considered and the distribution of bradyrhizobia in the specific locations. The abundance of the isolates clustered under the Bradyrhizobium sp. and Bd110 can be attributed to the influence of flooding period (FP), N, C, pH and clay content while the dominance of the isolates clustered under Bj6 can be attributed to the influence of high silt and P content in the soil. On the other hand, the abundance and distribution of the isolates under the B. elkanii clusters is mainly correlated with the high temperature in the country and possibly with the higher amount of sand. Since only a few isolates were identified to be related with B. yuanmingense, we could not make a definite correlation between these isolates and the agro-environmental parameters.

Fig. 4

Principal Component Analysis (PCA) plot depicting the relationship between the dominance of Bradyrhizobium species in the respective locations and the agro-environmental factors considered in this study. • Blackened circle indicates the agro-environmental factors. ♦ Blackened diamond indicates the location of the soil sampling collection. FP-period of flooding. Dotted straight line indicates the separation between the flooded and non-flooded soil condition.

Discussion

Genetic diversity of bradyrhizobia

Aside from the high species diversity as revealed from the Shannon's index, the genomic variations between the ITS region and the rpoB gene further showed a high genetic diversity. This is not a new observation as the same was reported previously (Mason et al., 2017) where genomic variations in the clusters Be76 and Be46 were detected and attributed to the temperature and pH gradient. Although the specific reason for this phenomenon cannot be explained accurately by our results alone, we hypothesize that a gene transfer might have occurred for the isolates which are classified under the Be76, Be46, and Be31 clusters. From the genome map of B. elkanii USDA76T with a genome size of 9,484,767bp, rpoB gene (6,539,500bp) is located on almost opposite locus with the ITS region (411,600–410,900bp) as reported previously (Reeve et al., 2017) so it may easily reflect genetic changes that occurred within the species, or even within the strain.

In general, the genomic variations between the two genetic loci on B. elkanii clusters are helpful for the identification of tropical bradyrhizobia and merit further investigation.

Factors affecting the diversity and distribution of bradyrhizobia

In this study, several agro-environmental factors were considered to investigate which parameter/s could impact the distribution and genetic diversity of soybean bradyrhizobia in the Philippines. The positive influence of Phosphorus content on the abundance of B. japonicum sp., which were phylogenetically clustered with Bj6 was also observed earlier (Yan et al., 2014). Another study stated that the Andisols in Japan were dominated by the isolates clustered to Bj6 (Shiina et al., 2014). This is similar to the soils of South Cotabato (SC) which are classified as Andisols with mixed alluvium and sedimentary deposits where dominant strains of Bj6 cluster was found in this study. Since B. japonicum USDA6T was reported to release N2O (Sameshima-Saito et al., 2006), which is one of the greenhouse gases, its impact in both agriculture and environment merits more attention for better understanding.

Meanwhile, the distribution of Bd110 cluster in flooded soils is supported by the previous reports where the dominance of B. diazoefficiens USDA110T was enhanced by flooding condition in the soil (Saeki et al., 2017) and that the USDA110 cluster is dominant on fine-textured soil that are affected by water status and oxidation-reduction potential in the soil (Saeki and Shiro, 2014). Additionally, it was reported that the anaerobic condition in flooded soil of alluvial origin resulted in the dominance of Bd110 cluster (Shiina et al., 2014) which is also similar to our results. The isolates clustered under Bd110 were found in areas which are usually planted with rice during wet season then, planted with rice and/or legume during dry season. In both seasons, rice farming is always done in waterlogged status. This provided an interesting insight particularly for the Philippines' agriculture since the strain B. diazoefficiens USDA110T was proven to be a highly effective and efficient inoculant for soybean (Siqueira et al., 2014; Liu et al., 2017) including its complete denitrification ability (Itakura et al., 2013; Shiina et al., 2014; Akiyama et al., 2016) that is beneficial for climate change mitigation.

In Nueva Ecija (NE), wherein the locations were almost similar in all aspect except for the soil water status, a different species of bradyrhizobia dominated each location. The NE1 isolates (dominated by B. elkanii) were collected from non-flooded condition while the NE2 isolates (dominated by Bradyrhizobium sp.) were collected from flooded condition. These results indicate that soil management could alter the population dominance of certain species of bradyrhizobia. An earlier report stated that the diversity and abundance of Bradyrhizobium species were altered by cultural management and other soil-related properties (Yan et al., 2014). We considered that these strains might be novel species in the Philippines since the stated locations have no history of rhizobia inoculation and we did not obtain any highly similar sequences from the BLAST engine for its identification.

As expected, the abundance and widespread distribution of isolates belonging to B. elkanii clusters in the Philippines is consistent with the previous findings wherein this species is distributed in areas with slight to moderate acidity and sub-tropical to tropical region (Saeki et al., 2006; Adhikari et al., 2012; Shiro et al., 2013; Mason et al., 2017). It might be one of the reasons for the low yield of soybean across the country. It was reported that B. elkanii provides lower N fixation and symbiotic efficiency (Risal et al., 2010) in comparison with B. diazoefficiens USDA110T. Additionally, B. elkanii species produce NO2-, and it is known that the interaction of nitrites to some soil components are related to environmental concerns. Thus, these species of bradyrhizobia are important to conduct research on for its role in agriculture and environment.

Conclusion

In relation to our objectives, we propose the following ideas: (i) the major micro-symbionts of soybean in the Philippines are the strains under the clusters of B. elkanii, Bd110, Bj6 and Bradyrhizobium sp., (ii) the dominance of the strains in the Philippines are classified accordingly: B. elkanii clusters for non-flooded soil with high temperature, Bd110 cluster for fine-textured flooded soils, and Bj6 cluster for flooded soils that are high in phosphorus and silt content, (iii) the isolated Bradyrhizobium sp. strains are endemic to the country and might be potential novel species, and (iv) the distribution and genetic diversity of bradyrhizobia in the country was mainly influenced by the period of flooding and soil properties (pH, soil type, nutrient). Hence, this study was able to identify the indigenous and potentially new endemic strains of bradyrhizobia, and identified the agro-environment conditions wherein they are abundant and dominant. Further studies on the strains’ characteristics, possession of denitrification genes, and symbiotic relationship with different soybean cultivars that are already adapted to the local conditions will be helpful to test its potential as an efficient and effective inoculant. This work was able to provide information that the ecology of bradyrhizobia is indeed complicated, particularly in a tropical region; and it depends on several abiotic and biotic factors.

Declarations

Author contribution statement

Maria L. T. Mason: Conceived and designed the experiments; Performed the experiments; Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper.

Baby L. C. Tabing: Performed the experiments.

Akihiro Yamamoto, Yuichi Saeki: Conceived and designed the experiments; Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper.

Funding statement

This work was supported by the JSPS KAKENHI (18K05376) and the Ministry of Education, Culture, Sports, Science and Technology (MEXT) Japanese Government Scholarship Program.

Competing interest statement

The authors declare no conflict of interest.

Additional information

Data associated with this study has been deposited at DDBJ under the accession numbers LC367064-LC367078, LC367080-LC367087, LC367091-LC367112, LC367114-LC367121, LC367125-LC367131, LC386868-LC386886, LC415435-LC415436.