Genetically related cowpea genotypes show similar bacterial community in the rhizosphere
The experiment was conducted in a greenhouse of the Department of Soil Sciences of the Federal University of Piauí, Brazil. The climate is tropical dry with an average precipitation of 1,300 mm per year.-1 (with rainfall from December to May) and an average annual temperature of 28°C. The soil used in this study is classified as Fluvisol, the type of soil where cowpea is commonly grown in northeast Brazil, and it was collected at a depth of 0-20 cm. The soil samples were taken from an area where cowpea has grown for the past ten years. Soil chemical properties were estimated according to Tedesco et al.44, being: pH – 6.1; Organic C – 9.2 g kg-1; P – 4.9mg kg-1; K – 32.5mg kg-1; base saturation (V) – 58%.
In this study, polyvinyl chloride pots (diameter 18 cm, length 16 cm) were filled with 5 kg of soil and the experiment was organized according to a completely random scheme with three repetitions. Two African lines (IT85F-2687 and IT82D-60) and two Brazilian cultivars (BRS-Guariba and BRS-Tumucumaque) of cowpea were selected after generation advance and genetic selection (Table 1; Fig. 6). In particular, the two modern cultivars are the most used by cowpea producers in Brazil.45. All genotypes are identified and deposited in the cowpea genetic resource bank of “Empresa Brasileira de Pesquisa Agropecuária – EMBRAPA”, in Teresina (Embrapa Meio-Norte), Piauí, Brazil. Permission to use these genotypes was obtained from Embrapa Meio-Norte. Five seeds of each genotype (modern lines and cultivars) were sown per pot and ten days after germination the plants were thinned, leaving one individual plant per pot. Pots were irrigated daily with sterilized water to maintain soil moisture at 80% of field capacity. The greenhouse temperature was controlled at 30°C.
Sampling of the rhizosphere
The plants were harvested 40 days after emergence (corresponding to the flowering period). Sampling of the rhizospheric soil was carried out as follows: the roots and adherent soil of each plant were placed on a 1 mm mesh screen and washed thoroughly with a gentle stream of tap water to remove soil. All soil samples were stored at -20°C prior to DNA extraction and analysis.
DNA extraction and sequencing
DNA was extracted from 0.5 g (total wet weight) of soil using the Power Lyzer Power Soil DNA Isolation Kit (MoBIO Laboratories, Carlsbad, CA, USA) according to instructions from the manufacturer. DNA extraction was performed in triplicate for each soil sample. The quality and concentration of the extracted DNA were determined using the NanoDrop 2000 spectrophotometer (Thermo Scientific, Waltham, USA).
The V4 region of the 16S rRNA gene was amplified with region-specific primers (515F/806R)46. Each 25 μL PCR reaction contained the following: 12.25 μL nuclease-free water (Certified Nuclease-free, Promega, Madison, WI, USA), 5.0 μL 5x buffer solution (MgCl2 2Mm), 0.75μL of dNTP solution (10mM), 0.75μL of each primer (515 YF 40 μM and 806 R 10 μM), 1.0 units Platinum Taq High Fidelity polymerase at a concentration of 0.5 μL (Invitrogen, Carlsbad, CA, USA) and 2.0 μL template DNA. Additionally, a control reaction was performed by adding water instead of DNA. The PCR conditions were: 95°C for 3 min to denature the DNA, with 35 cycles at 98°C for 20s, 55°C for 20s and 72°C for 30s, with a final extension 3 min at 72°C to ensure complete elongation.
After indexing, PCR products were cleaned using Agencourt AMPure XP – PCR purification beads (Beckman Coulter, Brea, CA, USA), according to the manufacturer’s manual, and quantified using the assay kit dsDNA BR (Invitrogen, Carlsbad, CA, USA) on a Qubit 2.0 fluorimeter (Invitrogen, Carlsbad, CA, USA). Once quantified, the equimolar concentrations of each library were pooled into a single tube. After quantification, the molarity of the pool was determined and diluted to 2 nM, denatured, then diluted to a final concentration of 8.0 pM with a PhiX peak (Illumina, San Diego, CA, USA) at 20% for loading into the Illumina MiSeq sequencing machine (Illumina, San Diego, CA, USA).
Sequence data was processed using QIIME 2 version 2019.10. First, sequences were demultiplexed and quality control was performed using DADA247, using the consensus method to remove all remaining chimeric and low-quality sequences. After filtering, approximately 355,000 high quality sequences were obtained, with an average of approximately 19,600 sequences per sample. Singletons and doubletons were removed and samples were rarefied to 8600 sequences, following the lowest sample number. Taxonomic affiliation was achieved at 97% similarity using the Silva v. 13248, and the generated matrix was then used for statistical analyses. The sequences were submitted to the NCBI Sequence Read Archive under the identification PRJNA751574.
To compare bacterial community structure between treatments, we used principal component analysis. First, the data matrix was initially analyzed using trendless correspondence analysis (DCA), indicating a linear distribution of the data and the best-fitting mathematical model was PCA. To test whether different niches (i.e., loose soil and rhizosphere) and genotypes of cowpea harbored significantly different active bacterial communities, we performed two-way permutational multivariate analysis of variance (PERMANOVA). PCA analysis was performed using Canoco 4.5 (Biometrics, Wageningen, The Netherlands) and PERMANOVA using Past v.449. We also used Past to calculate the observed OTU richness and Shannon’s diversity index, and the comparison was based on Tukey’s HSD test. A Venn diagram was constructed to verify the proportion of exclusive and shared groups between treatments using the web tool Venny 2.150. To determine the differential composition of bacterial groups between treatments, we used the STAMP software51. P values were calculated based on the two-tailed Welch t-test and correction using the Benjamini-Hochberg FDR. Additionally, to investigate the correlation between bacterial groups and performance properties of cowpea plants, we used Spearman’s rank correlation coefficient using R52. We then used network analysis to assess the complexity of interactions between bacterial taxa in each treatment. For this, a non-random co-occurrence analysis was performed using the Python module ‘SparCC’53. For this, an OTU frequency table was used for the analysis, and SparCC correlations were calculated, and only strong (>0.9 or p 54. Finally, to better predict relevant potential functions performed by the bacterial community, we performed functional annotation using FAPROTAX53, a database that maps prokaryotic clades (eg, genera or species) with established metabolic functions or other ecologically relevant functions, based on cultured strains. For this, a genus-level taxa frequency table was used as input and converted into a putative functional table and comparisons were made between bulk soil and rhizosphere of different cowpea genotypes.