Study area and inclusion criteria

The research was conducted in the Alto Chicamocha irrigation and drainage district, which includes the municipalities of Paipa, Duitama, Tibasosa, Firavitoba, Santa Rosa de Viterbo, Sogamoso and Nobsa in the provinces of Tundama and Sugamuxi located in the department of Boyaca (Fig. 1). The average climatic conditions in the provinces are: Average temperature of 14°C, average precipitation of 778 mm, relative humidity of 70%, and an average altitude of 2,560 m a.s.l. (Suppl. Tab. 1, https://revistas.uptc.edu.co/index.php/ciencias_horticolas/article/view/15702/13005) (^{Martínez et al., 2014}). Most of the land is used for agriculture, with some dairy cattle ranching and, to a lesser extent, dual-purpose cattle ranching (^{GISSAT-UPTC, 2012}).

Figure 1. Location of the sampling zones in the study. 

Sample selection criteria

The farms chosen for this study had the following selection criteria:

Sampling was carried out on farms with active production in August and September 2021, located in the municipalities of Duitama and Tibasosa.

Field and laboratory phases

In each zone sampled, a soil sample and an irrigation water sample associated with the crop were obtained. Each sampled zone received a paired coding for the soil and irrigation water sample as T and the farm number (i.e., T1) (Suppl. Tab. 1). Soil samples were taken in a 50×50 m quadrat on each farm. In each quadrat a zigzag path was made and 20 random sampling points were taken at a depth of 15 cm, for a total of 200 g of sieved soil sample and 2 soil clods (maximum allowed by the crop owners). Each of the soil samples was stored in dark bags and sent to the laboratory for physicochemical analysis in 15 mL dark falcon tubes in cold chain (-4°C) for transport and subsequent DNA extraction. The irrigation water samples were collected in sterile plastic bottles to avoid contamination by salts or other substances. From each sampling zone, 3 L of irrigation water were collected, stored, and sent to the laboratory for the corresponding analyses. The analyses of the physicochemical properties of the soil and water samples were carried out in the soil and water laboratory of the Faculty of Agricultural Sciences of the UPTC (Universidad Pedagógica y Tecnológica de Colombia), which uses the methodologies proposed by the Instituto Geográfico Agustín Codazzi (^{Colombia IGAC, 2006}). The physicochemical properties and methods implemented are shown in table 1.

The postharvest crop attributes were collected from the information provided by each grower. The attributes were size and weight of the fruit, the productivity in tons and the age of the crop, to classify the crops in the following categories: Q1 = tomato type "1^st Quality: Crop with the fruit size/weight ratio in quantile 1"; Q2= tomato type "2^nd Quality: Crop with the fruit size/weight ratio in quantile 2"; Q3= tomato type "3^rd Quality: Crop with the ratio fruit size/weight in the quantile 3". The attributes analyzed were: production size of the crop quantified in tons per hectare (t/ha), age of the crop quantified in years, the geographical altitude and the planting density in plants per hectare (plants/ha).

Nucleic acid extraction, library amplification, and sequencing

DNA was extracted in triplicate from each soil and irrigation water sample (technical replicate) using the Qiagen soil DNA extraction kit and the Qiagen water DNA extraction kit. The DNA concentration was quantified by spectrophotometry (Epoch™ Microplate Spectrophotometer, BioTek). Once the DNA extraction of the three technical replicas for each sampling point was carried out, these were integrated into a single sample to carry out the sequencing process. Samples were sequenced on the Miseq sequencing platform (2×300 bp paired-end protocol with a depth of 100 k reads) using the universal primers 341F 5′- CCTACGGGGNGGCWGCAG-3′ and 805R 5′-GACTACHVGGGTATCTAATCC-3′ to target to the V3 and V4 regions of the 16S rRNA gene (~460 bp) (^{Klindworth et al., 2013}; ^{Ren et al. 2017}).

Bioinformatic analysis of 16S rRNA gene amplicon library sequences

The raw sequence data were processed using the QIIME 2 v. 2019.7 packages (^{Caporaso et al., 2010}). Quality filtering was performed on the merged sequences and sequences that did not meet the following criteria were discarded: sequence length < 200 bp, no ambiguous bases, and average quality score ≥ 20. The sequences were then clustered using Uclust (^{Edgar, 2010}) for subsequent identification and removal of chimeric sequences using Uchime (^{Edgar et al., 2011}). Effective sequences were clustered and classified into Amplicon Sequences Variants (ASVs) with a > 97% similarity threshold using Vsearch and Classify-sklearn (^{Caporaso et al., 2010}) taxonomy assignment was performed using Silva v. 138 as the reference database (^{Quast et al., 2013}).

Data availability

The data presented in the study are deposited in the SRA repository, accession number PRJNA936423.

Statistical analysis

The variability of the physicochemical properties of the soil and water was evaluated through descriptive statistics. Contingency tables were constructed to analyze quantitative and qualitative variables. Multivariate exploratory data analyses and Pearson correlations between physicochemical properties and crop attributes were executed using the FactoMineR and Hmisc package (^{Husson et al., 2014}; ^{Harrell and Dupont, 2023}). Non-parametric multidimensional scaling (NMDS) was performed using the Vegan package (^{Oksanen et al., 2017}) to observe which soil and water properties could present possible groupings or relationships in each of the sampled areas. A canonical correlation analysis (CCA) was performed using the CCA package (^{González and Déjean, 2021}) to find the relationship between the sets of physicochemical properties and crop attributes for each of the sampled areas. The graphical visualization of each of the procedures was performed using the GGvegan package (^{Simpson, 2019}). The diversity of the bacterial community was evaluated by the index of observed ASVs. Likewise, the relationships between the composition of the bacterial communities present in the soil and irrigation water samples and the different crop quality traits were evaluated by principal coordinate analysis (PCoA) of the Jaccard index. Distances and clustering were evaluated by the ANOSIM test, using the ggplot2 package v. 3.3.1 (^{Wickham, 2011}). Finally, to simplify the visualization of the bacterial community related to each quality of the Chonto tomato crop, the bacterial core was established at the phylum and genus level using the Microbiome package (^{Lahti and Shetty, 2019}). All the packages were run on the free software Rstudio v. 4.2.1.

Physicochemical properties

The lowest pH value (5.42) was given by sample T4 and the highest (7.73) by sample T10. The lowest value of EC (dS cm^-1) was (0.81) given by sample T15 and the highest (8.98) by sample T2. The OM% values remained in a range of (1.49-6.15) with sample T15 being the lowest and sample T14 the highest. Regarding nutrition, the highest values of Mg, K, and Na (cmol₊ kg^-1), were 3.67, 2.75, and 5.25 corresponding to the same sample (T2), while the lowest values were 0.52 Mg, 0.37 K, and 0.33 Na belonging to samples T8, T6, and T10 respectively. Finally, there was a high proportion of soil texture, clayey (Ar), followed by sandy clay loam (F-Ar-A) (Tab. 2). The alterations in the physical and chemical properties of the soil are caused by long-term saline irrigation, seen in plots with alkaline-saline soils, reduced structure index, low water infiltration rate, high water content, poor aeration of the root zone, superficial crusting and pH ≥8.0 (^{Bouksila et al., 2013}; ^{Cabrera et al., 2016}; ^{Bonachela et al., 2018}).

In the physicochemical values of the irrigation water, the pH did not significantly change at increasing salinity of the irrigation water or of the soil, and was similar to that reported by ^{Maggio et al. (2004)}. A pattern was observed with sample T2, which, in this case, is the sample that presented the lowest values of EC (dSm^-1), hardness, Mg, Na, Cl (cmol₊ kg^-1) and HCO₃ (Meq L^-1). Likewise, we highlight that for samples T3 and T7, extremely low values of SO₄ (Meq L^-1) were obtained. Finally, high hardness and electrical conductivity were found in the T13 sample, with hardness at 2.95 and EC at 851 dS cm^-1, which could indicate a possible salinity problem in the irrigation water (Tab. 3).

The relationship between the physicochemical variables of the soil and the different characteristics of the tomato crop was studied using the Pearson coefficient. Among the most striking results was the high inverse correlation between pH value and productivity (r=-0.82; P<0.05) and the mean inverse correlation between EC and Ca values and crop altitude (r=-0.44; P<0.05). The only major positive correlation was found between soil bulk density with crop altitude (r=0.38; P<0.05) (Tab. 4).

The values obtained by Pearson correlation in the physicochemical properties of irrigation water showed an inverse mean correlation between the variables EC and HCO₃ with crop productivity (r=-0.31; P<0.05). The highest positive mean correlation was seen between K and crop planting density (r=0.38; P<0.05), followed by turbidity and productivity (r=0.33; P<0.05) (Tab. 5).

Salinized soils over time can undergo permanent changes in their physicochemical properties, further aggravating the damage to crop plants caused by excess salts. Hyperosmotic stress, hypoxia, and nutritional imbalance often coexist in soils exposed to salinization and persist throughout the growing season (^{De Pascale et al., 2003}). Tomato is moderately tolerant to salinity (1.3 dS cm^-1< electrical conductivity of the saturated soil < 6 dS cm^-1) and is typically cultivated in regions that are exposed to soil salinization (^{Maggio et al., 2004}). The best position for monitoring the concentration of nutrients and salts appears to be the center of the wet bulb, because this position might respond faster to changes in the nutrient solution supplied or the root activity, especially for very mobile elements, such as nitrate (^{Bonachela et al., 2018}). The bulk density values reported in this study range from 0.75-1.17, classifying them in the order of entisols, similar to that reported by ^{Alvarado and Forsythe (2005)}. The world's soils are grouped into nine orders ranging from 0.14-2.00 (^{Davis, 1959}; ^{Alvarado et al., 2001}). ^{Shaheb et al. (2021)} observed that in loam-textured soils the bulk density increased due to the intensive use of agricultural machinery. In addition, the bulk density in clayey soils was lower than in loam soils, which is explained by the greater porosity that clayey soils generally have (^{Jaurixje et al., 2013}).

The non-metric multidimensional scale showed medium-to-low relationships between the physicochemical properties of the soil (Fig. 2A) and the irrigation water for tomato cultivation (Fig. 2B), showing highly conservative patterns within the multidimensional space represented. In general, the samples are distributed along the physicochemical conditions, without forming any type of grouping corresponding to geographic distribution or any macrogeographic aspect. It should be noted that samples T2 (Fig 2A) and T4 (Fig. 2B) present a more than considerable distance with respect to the other samples. This could be due to extreme values in Na nutritional inputs (Fig. 2A) and irrigation water turbidity (Fig. 2B). Even so, in general, all these variables show an extremely low variability within the multidimensional space. This could indicate the ability of greenhouse crops to maintain at least supremely homogeneous soil conditions along a geographical gradient.

Figure 2. The non-metric multidimensional scale (NMDS) relates the sampling zones to the physicochemical properties of soil (A) and irrigation water (B) associated with the tomato crop. 

The relationships between the physicochemical variables of soil and irrigation water with the attributes of crop productivity and quality and their geographic association were analyzed by canonical correspondence analysis (Fig. 3). In general, it was observed that for soil physicochemical properties, fruit quality is highly related to soil P and pH values (P<0.05). Likewise, high P and pH values, which are related to fruit quality, are inversely proportional to the crop productivity (Fig. 3A). This interaction could be related to the fact that a living biological system with smaller or poorer quality fruits could optimize the nutritional resources available in the soil and, therefore, obtain a greater number of fruits (^{Li et al., 2019}; ^{An et al., 2023}). For the physicochemical values of irrigation water, there was low variability except for hardness values, which are significantly related to productivity and crop age (P<0.05) (Fig. 3B). In addition to this and considering that the low variability of soil physicochemical values may be linked to the type of crop (greenhouse), the low variability in irrigation water samples could be related to mechanical conditions. This is a response to the combination of irrigation water sources and local adaptations, resulting in differential conditions in the irrigation water samples. Thus, and considering the results, we can support the idea that greenhouse crops could have extremely high effects on soil nutritional loads, because the microclimatic conditions that are generated in one way or another homogenize the soil conditions in these crops that are of high nutritional demand (^{Usero et al.,2021}).

Figure 3. Tri-plot of the Canonical Correspondence Analysis (CCA) between sampling zones and the tomato crop productivity attributes. A) soil physicochemical properties and the squares listed correspond to: 1 = P mg kg^-1, 2 = Mg cmol⁺ kg^-1, 3 = bulk density g cm^-3, 4 = K cmol⁺ kg^-1, 5 = EC dS cm^-1, 6 = Na cmol⁺ kg^-1, 7 = Ca cmol⁺ kg^-1, 8 = OM %, 9 = pH. B) irrigation water properties and the squares listed correspond to: 1 = hardness (°F), 2 = turbidity NTU, 3 = EC dS cm^-1, 4 = Ca meq L^-1, 5 = Mg meq L^-1, 6 = SO₄ meq L^-1, 7 = Cl meq L^-1, 8 = K meq L^-1, 9 = Na meq L^-1, 10 = pH. 

Bacterial composition associated with the crop

Of the 15 samples sequenced, 14 viable irrigation water samples and 12 viable soil samples associated with the Chonto tomato crop were obtained (Suppl. Tab. 2, https://revistas.uptc.edu.co/index.php/ciencias_horticolas/article/view/15702/13006). Likewise, 1,333 consensus amplicon sequence variants (ASVs) were obtained for all irrigation water samples and 1,718 ASVs for the Chonto tomato crop microbiota profile soil samples associated at the species level. After rarefaction and collapse of each sample type, the ASV table was assigned to 310 genera for the irrigation water microbiota profile and 340 genera for the soil microbiota profile.

To evaluate the relationship between bacterial diversity/composition and the quality of the fruit of the Chonto tomato crop, the samples were grouped by quality and alpha and beta diversity analyses were performed. The comparison of alpha diversity based on the three observed qualities showed significant variations (P<0.05) in the bacterial diversity evaluated. It was observed that quality Q1 (1^st Quality: Crop with the fruit size/weight ratio in quantile 1), which is the best fruit quality, has a lower number of bacterial ASV associated with irrigation water with respect to quality Q2 (2^nd Quality: Crop with fruit size/weight ratio in quantile 2) and Q3 (3^rd Quality: Crop with fruit size/weight in the quantile 3) (Fig. 4A). For β-diversity, PCoA analysis was performed based on the Jaccard index in combination with ANOSIM. Bacterial composition was significantly different (P<0.001) between fruit quality Q1 and Q2 vs. Q3. ANOSIM results showed that the greatest variation in bacterial composition was always against Q3 (R=0.688, P<0.001), while comparing Q1 and Q2 fruit quality the bacterial composition was more like each other (R=0.245, P>0.05), but not statistically significant (Fig. 4B).

While the comparisons evaluated for the associated microbiota for the irrigation water samples, both alpha and beta diversity comparisons were found to be non-significant (P>0.05) (Suppl. Fig. S1, https://revistas.uptc.edu.co/index.php/ciencias_horticolas/article/view/15702/13007).

Figure 4. Alpha (observed ASVs) and beta (Jaccard index) diversity analysis of the fruit quality microbiota of the Chonto tomato crop. (A) Boxplot plots visualize the number of observed ASVs for fruit qualities Q1, Q2 and Q3. Significant differences (P≤0.05) were assessed by Kruskal Wallis paired t-test (Pairwise group A, B and C). (B) PCoA plots showing distances between bacterial community composition of fruit qualities Q1, Q2 and Q3. Significant differences in bacterial composition were tested using the ANOSIM pairwise test. Each fruit quality type is associated with a color bar (A) and dots (B). 

At the level of bacterial composition, the samples associated with irrigation water presented a similar composition of ASVs for the three quality groups evaluated (Q1, Q2 and Q3). The bacterial ASVs associated with the water samples do not offer a homogeneous taxonomic resolution at the genus level, however, the bacterial groups observed are highly conserved and in abundance and dominance in the three groups of samples evaluated (Fig. 5A).

Figure 5. Relative abundance of bacterial community at genus level associated with fruit quality in the Chonto tomato crop. Bacterial genera in irrigation water samples (A) and soil samples (B). Only taxa with occurrence greater than 0.1% of total relative abundance are presented. 

The composition of the tomato crop microbiota was related to fruit quality and planting soil. The abundance of the genus Sphingomonas spp. (~38%), which is more dominant in Q1 and Q2 quality soil, was substantially reduced in Q3 quality soil (Fig. 5B). Likewise, while the genus Sphingomonas spp. showed differential abundance between Q1 - Q2 versus Q3 at the soil level (ANOVA: P<0.05), the genus Sphingobium spp. was highly dominant and abundant (42%) in the samples related to the planting soil for fruit quality Q3 (Fig. 5B).

Currently, it has been demonstrated that the bacterial microbiota associated with tomato cultivation plays an important role in plant growth and development through nutrient uptake and environmental adaptation (^{Schmitz et al., 2022}). In this study, we identified that the soil bacterial microbiota has a direct relationship with fruit quality of the Chonto tomato crop under greenhouse conditions. Bacterial groups such as Sphingomonas spp., previously reported in studies of other tomato varieties have a relevance in all phenological stages of the plant, mainly as protection against diseases affecting the fruit (^{Khan et al., 2012}; ^{Anzalone et al., 2022}; ^{Zhang et al., 2022}). Likewise, several studies have shown that soil-associated bacteria show greater fluctuation in their composition and abundance (^{Manzotti et al., 2020}). This confirms that the structuring of the crop soil-associated microbiota community is primarily determined by the niche of the planting site and the host species, while location (and site-related soil and climate factors) had little effect on effect on the microbiota.

The impact of rhizobiales is influenced by root exudates, which also play an important role in signaling and recognition processes between plants and microorganisms (^{Venturi and Keel, 2016}). The tomato crops evaluated presented high levels of salinity that do not allow the plant to have the optimal medium to grow, but these crops showed mostly good fruit growth, which is explained by the symbiotic relationship of Rhizobial bacteria found in the metagenomic analysis, as these bacteria help the plant to withstand the stress caused by high temperatures, drought, salinity and infection by pathogens (^{Foyer et al., 2016}; ^{Gouda et al., 2018}). Rhizobacteria have been studied more in legumes and in terms of biological nitrogen fixation. However, studies suggest that most rhizobacteria have the potential to promote plant growth by mechanisms other than biological nitrogen fixation or biocontrol of certain microorganisms harmful to plants or to regulate plant development under stress conditions, both in legumes and non-legumes (^{Singh et al., 2010}).

Groups of rhizobacteria capable of solubilizing minerals such as potassium have been identified and have been associated with increased productivity (^{Parmar and Sindhu, 2013}), where the bacterial genus Sphingomonas has been reported as a rhizobacterium capable of solubilizing potassium (^{Etesami et al., 2017}). Likewise, bacterial groups described within the genera Sphingomonas and Sphingobium capable of degrading these compounds have been identified and isolated from soils contaminated with different petroleum and coal hydrocarbons (^{Lal et al., 2006}; ^{Prakash and Lal, 2006}; ^{Dadhwal et al., 2009}). The few studies carried out on these genera have been associated with wastewater and soils contaminated with hydrocarbons. The need to further characterize the microbiota present in crops and to carry out applied research to better understand the relationship or presence of certain bacterial groups present in soil or water is emphasized.

Finally, the little variation observed in the structure of the bacterial community associated with the irrigation water samples of the Chonto tomato crop can be attributed mainly to the effects of the pretreatments carried out within the irrigation infrastructure of the association, as previously stated. This has been demonstrated before in other studies of tomato cultivation (^{Ma et al., 2020}), suggesting and supporting the hypothesis that the bacterial communities present in irrigation water are not the main source of variation in the quality of the fruit of the crop.