Continued Organic Fertigation after Basal Manure Application Does Not Impact Soil Fungal Communities, Tomato Yield or Soil Fertility

There is currently a limited understanding of the complex response of fungal microbiota diversity to organic fertigation. In this work, a 2-year field trial with organic tomato crops in a soil previously amended with fresh sheep manure was conducted. Two hypotheses were compared: (i) fertigation with organic liquid fertilizers versus (ii) irrigation with water. At the end of both years, soils were analyzed for physical–chemical parameters and mycobiome variables. Plate culture and DNA metabarcoding methods were performed in order to obtain a detailed understanding of soil fungal communities. Fertigation did not increase any of the physical–chemical parameters. Concerning soil fungal communities, differences were only found regarding the identification of biomarkers. The class Leotiomycetes and the family Myxotrichaceae were identified as biomarkers in the soil fungal community analyzed by means of DNA metabarcoding of the “fertigation” treatment at the end of Year 1. The Mortierella genus was detected as a biomarker in the “water” treatment, and Mucor was identified in the “fertigation” treatment in the cultivable soil fungi at the end of Year 2. In both years, tomato yield and fruit quality did not consistently differ between treatments, despite the high cost of the fertilizers added through fertigation.


Introduction
Improving soil health and reducing agricultural inputs are crucial factors in ensuring sustainable and profitable farming practices. These are key elements to reversing the degradation of agrosystems and, at the same time, to addressing the substantial threat of climate change and population growth on a global scale. Hence, the implementation of efficient agricultural techniques aimed at reducing the impacts generated by productive systems in recent decades has emerged as a critical priority on the international political agenda. In this context, organic farming (OF) has been identified as playing a key role in realizing the Sustainable Development Goals outlined in the United Nations' Agenda 2030, as well as in European Union policies such as the EU Green Deal's Farm to Fork initiative or those contained in the EU Soil Mission program.
OF has experienced significant growth in recent years, as reflected in the constant increase in certified farming areas in accordance with Regulation (EU) No 2018/848. This development is evident throughout the European Union (EU), where approximately 14.7 million hectares of land assigned to OF was recorded in 2020 (including areas in conversion), representing a 56% increase between 2012 and 2020, and accounting for 9.1% of the total utilized agricultural area (UAA). This increase is also evident in Spain, the EU country with the second-largest area dedicated to OF, after France. In 2020, Spain had over 2.43 million hectares of organic production, representing nearly 10% of the country's UAA,

Plant Material, Cropping Details and Experimental Design
Two subsequent winter cycles of tomato (Solanum lycopersicum L.) "Valenciano type" plants, grafted onto Armstrong ® rootstock (Syngenta, Basel, Switzerland), were grown in the same greenhouse in the seasons 2019/20 and 2020/21 (i.e., Year 1 and Year 2, respectively). Four-week-old tomato plants were planted on 17 September 2019 and 25 September 2020, in Years 1 and 2, respectively. In both years, prior to planting, in July, fresh sheep manure was buried uniformly across all crop rows in the greenhouse at a rate of 4 kg m −2 . The incorporation was performed by removing the top layer of sand and mixing the manure with the soil using a rotavator. The area was then covered with the same sand that was previously removed. Next, drip lines were installed and, with the aim of favoring the manure's decomposition, the soil was covered with transparent polyethylene film (30 µm TIF Desinfección DS ® , Sotrafa, Almería, Spain) for a period of two months, after a single irrigation application to reach saturation at a 15 cm depth, which is known as the soil biosolarization technique. Planting took place two days after removing the film. Crop growth occurred on two axes, considering each as an individual plant for sampling purposes, thus resulting in an overall density of 2 plants m −2 . Tomato vines were vertically trained with polypropylene strings and pruned and managed according to established local practices. For correct and optimal pollination, bumblebees (Bombus terrestris) were used. Crop management and pest control were guaranteed by adhering to Regulation (EU) 2018/848 on organic production. As there were no soil-borne pathogens present in the soil, no soil treatments were necessary during the crop periods.
The study considered two treatments (Tables 1 and 2): (i) "water" received only the initial application of manure and was irrigated solely with water throughout the crop cycles and (ii) "fertigation" supplemented the manure with fertilizers listed in the Commission Implementing Regulation (EU) 2021/1165 on products and substances for use in organic production, which were applied through fertigation. The experiment was performed as a single-factor design with three replications (n = 3). The trial consisted of a total of 6 experimental plots, each being 96 m 2 in size. Table 1. Treatments considered in the study.

Water
Sheep manure at a rate of 4 kg m −2 . Irrigated solely with water throughout the crop cycle.

Fertigation
Sheep manure at a rate of 4 kg m −2 + a fertilization plan using fertilizers listed in the Commission Implementing Regulation (EU) 2021/1165 on products and substances for use in organic production which was applied through fertigation throughout tomato crop cycle ( Table 2).

Soil Sampling
Soil samples were taken at three different sampling times. To assess the starting situation, soil sampling was carried out in three areas of the greenhouse immediately before planting the tomato seedlings. Thus, the soil samples were taken after the biosolarization treatment with fresh sheep manure, two days before planting (i.e., initial, 16 September 2019). Similarly, in both years, at the end of the tomato crop season (i.e., end of Season 1, 3 April 2020, and end of Season 2, 23 April 2021, in Year 1 and Year 2, respectively), all of the experimental plots were sampled. After removing the top layer of sand, the samples were taken with an auger at a depth of 0-30 cm. Three subsamples were randomly taken from cultivation lines and then mixed and homogenized to ensure representativeness in each sample. Subsamples were taken in the center of the cultivation lines to avoid a possible edge effect. For each soil sample, one subsample was frozen at −80 • C for DNA extraction, another subsample was air-dried and sieved (<2 mm) for chemical analysis, and another subsample was air-dried and sieved (<0.2 mm) for microbiological analysis through the plate culture method.

Analysis of Physical-Chemical Characteristics of Soil Samples
Soil samples (n = 3) were analyzed for pH, electrical conductivity (EC), interchangeable cations (Ca +2 ; Na + ; Mg +2 ; K + ), active limestone, phosphorus Olsen (P Olsen), nitric nitrogen, organic matter, total carbonates (HCO 3 − ), total nitrogen, and C/N ratio, following standard soil testing procedures as described in Order 5/12/1975 [24] by an external laboratory (Eurofins, El Ejido, Spain). The results for most physical and chemical variables of soil samples were reported in different units for Years 1 and 2. Comparisons among the treatments were made for each year separately.

Study of Soil Fungi in Soil Samples
The study of soil fungal microbiota ("mycobiome") was conducted using two evaluation methods: (i) the plate culture method (culturable fungal microbiota) and (ii) a molecular method based on DNA metabarcoding.

•
Preparation of soil samples Soil samples were subjected to a drying, grinding, and sieving process, following Tello et al. [25]. Drying was conducted at room temperature (20-25 • C) for 7-10 days, until soil humidity was constant and homogenous. A porcelain mortar was used for grinding, and a 200 µm mesh-size sieve was used for sieving. The mortar and the sieve were washed and disinfected between samples by covering them with alcohol and lighting it.

•
Analytical method The soil culturable fungal microbiota were analyzed by means of the serial dilution method [25]. The culture medium used was potato dextrose agar (PDA). Five subreplicates (i.e., Petri dishes) of each soil sample (n = 3) were created at 10 −3 and 10 −4 dilutions. The Petri dishes (9 cm diameter) were incubated at 25 • C for 4-7 days. Subsequently, total colony forming units (CFU) of fungi were quantified and morphological identification at the genus level was performed [26,27], eventually expressing the results as CFU/g dry soil. The isolates unidentified by morphological means were analyzed by means of polymerase chain reaction (PCR) sequencing of amplicons of the internal transcribed spacer region (ITS) rDNA using the primers ITS4 and ITS5 [28] and subsequent database searches using BLASTN software, based on the consensus sequences created by aligning the forward and reverse sequences of the target isolates. The PCR conditions were 5 min at 94 • C, 35 cycles of 1 min at 94 • C, 1 min at 55 • C, 2 min at 72 • C, and a final elongation of 7 min at 72 • C. Purified PCR products were sequenced using the Sanger sequencing method (Instituto de Biología Molecular y Celular de Plantas, Valencia, Spain) and edited via base calling on the BioEdit program.

•
Soil DNA Isolation and Quantification DNA from soil samples (n = 3) was isolated using the DNeasy PowerSoil Pro DNA isolation kit (Qiagen, Hilden, Germany), strictly following the manufacturer's instructions. An extraction blank for cross-contamination was included. DNA was quantified using the Qubit High Sensitivity dsDNA Assay (Thermo Fisher Scientific, Waltham, MA, USA).

• Library Preparation and Sequencing
For fungal library preparation, a fragment of the ITS2 genomic region (of about 350 bp) was amplified using the primers ITS86F (5 GTG AAT CAT CGA ATC TTT GAA 3 ) [29] and ITS4R (5 TCC TCC GCT TAT TGA TAT GC 3 ) [28]. Illumina sequencing primers were attached to these primers at their 5 ends.
PCRs were carried out in a final volume of 12.5 µL, containing 2.5 µL of template DNA, 0.5 µM of primers, 6.25 µL of Supreme NZYTaq 2× Green Master Mix (NZYTech, Lisbon, Portugal), 1× CES [30], and ultrapure water at a volume up to 12.5 µL. The reaction mixture was incubated as follows: an initial denaturation step at 95 • C for 5 min, followed by 35 cycles of 95 • C for 30 s, 49 • C for 45 s, 72 • C for 45 s, and a final extension step at 72 • C for 7 min. A negative control that contained no DNA (BPCR) was included in every PCR round to check for contamination during library preparation. The libraries were run on a 2% agarose gel stained with GreenSafe (NZYTech) and imaged under UV light to verify the library size. Libraries were purified using the Mag-Bind RXNPure Plus magnetic beads (Omega Biotek, Beijing, China), following the instructions provided by the manufacturer. Then, libraries were pooled in equimolar amounts according to the quantification data provided by the Qubit dsDNA HS Assay (Thermo Fisher Scientific). This pool also contained a testimonial amount (1 µL) of the PCR blank and DNA extraction blank (Bex). The pool was sequenced in the MiSeq PE300 run (Illumina). DNA metabarcoding analyses were carried out by AllGenetics & Biology S.L. (La Coruña, Spain; www.allgenetics.eu).

•
Quality Control and Processing of Sequencing Data Illumina paired-end raw forward (R1) and reverse (R2) FASTQ reads were stored in separate files. The quality of the FASTQ files was assessed with the software (v0.11.9) FastQC [31] and the output was summarized using MultiQC [32]. The obtained amplicon reads were processed using QIIME 2 (release 2022.2) [33]. Specifically, the tool DADA2 [34] (implemented in QIIME 2) was used to remove the PCR primers, filter low-quality reads, denoise, merge the forward and reverse reads, remove the chimeric reads, and cluster the resulting sequences into amplicon sequence variants (ASVs). Cutadapt [35] was firstly used to remove primer and/or adapter sequences. The resulting output of the DADA2 pipeline was a table containing the number of occurrences of every observed ASV in each sample. The taxonomic assignment of ASVs was conducted using a pre-trained classifier of the UNITE reference database [36] (updated in May 2021). Specifically, the feature-classifier classify-sklearn approach implemented in QIIME 2 [37] was employed. The following ASVs were excluded: singletons (i.e., ASVs containing only one member sequence in the whole dataset), ASVs occurring at a frequency below 0.01%, unassigned sequences, and sequences assigned only at the kingdom level. The final filtered ASV table was converted into a Biological Observation Matrix (biom) file that was directly imported into R 3.6.1 (R-Team-Core, 2019) using the package phyloseq 1.24.2 [38]. The final filtered ASV and the biom file were used for further analysis and plotting.

Fungal Community Description and Data Procesing
Four classical indexes of biodiversity were selected to describe the structure of the fungal community of each sample: Margalef genus richness (d) [39]; the Shannon-Wiener diversity index (H ) log basis [40]; Pielou's evenness (J) [41] of the distribution of individuals among genera; and Simpson's dominance (D) [42]. All calculations were carried out with the software PRIMER version 6.0 (Primer-E Ltd. 239 Plymouth, UK) for Windows [43]. Principal coordinate analysis (PCoA) was used to visualize the potential differences in the soil fungal community structure between treatments or within treatments. Fungal community data (at the genus and operational taxonomic unit (OTU) levels, for plate culture and DNA metabarcoding evaluation methods, respectively) were compared between treatments using permutational analysis of variance (PERMANOVA). This test calculates Pseudo-F statistics to obtain p-values [44]. All PERMANOVA tests used 9999 permutations from unrestricted permutations of raw data. Where there were fewer than 99 unique permutations for a meaningful test, approximate Monte Carlo p-values (p(MC)) were obtained from an asymptotic permutation distribution [45]. Similarity percentage (SIMPER) analysis was used to calculate the average contribution of the abundance of individual culturable fungal genera and OTU to the average dissimilarity of the soil fungal community of the two treatments. For all PERMANOVA and SIMPER analyses, data were square root transformed. The calculation of similarity matrices in the PERMANOVA analysis was based on Bray-Curtis distance. PCoA, PERMANOVA, and SIMPER analyses were performed using the multivariate statistical software package PRIMER-6 with the PERMANOVA+ extension (Primer-E Ltd., Plymouth, UK) for Windows.
The linear discriminant analysis effect size (LEfSe) algorithm [46] was used to identify and compare unique fungal taxa significantly related to any of the two treatments (i.e., fertigation and water) for the different seasons. The threshold for the logarithmic linear discriminant analysis (LDA) score was set at 2.0 and the Wilcoxon p-value was set at 0.05. LEfSe was conducted for the results obtained through the two evaluation methods (plate culture method and DNA metabarcoding method).
The FUNGuild database [47] (accessed on 22 March 2023) was used to classify fungal ASVs into functional groups. The FUNGuild software annotates taxonomic data within the OTU table with corresponding data on its online database. The annotations include the guild, trophic mode, and growth morphology, and only confidence scores of "Probable" and "Highly Probable" were considered. The ASVs with confidence scores of "Possible" and those that ended up without functional assignment were placed into the "Unclassified" category.
2.6. Tomato Crop Analyzed Variables 2.6.1. Tomato Crop Production The production was measured for all of the harvests in the 6 experimental plots (2 treatments with 3 replications each). In Year 1, the first harvest was undertaken on 27 December 2019 (101 days after planting (DAP)) and the last was performed on 25 March 2020 (191 DAP). In Year 2, the first harvest was on 30 December 2020 (96 DAP) and the last was performed on 13 April 2021 (200 DAP). In total, 14 and 15 harvests were undertaken in Year 1 and Year 2, respectively. The cumulative tomato production was measured (kg m −2 ) using an electronic balance with an accuracy of ±0.01 kg.

Tomato Fruit Quality
The attributes of tomato fruit quality were measured at 2 different times during each crop cycle (i.e., Season 1 and Season 2): on 8 January 2020 and 12 February 2020 in Season 1, and on 19 February 2021 and 20 April 2021 in Season 2. Three and five fruits, in Season 1 and Season 2, respectively, were taken from each experimental plot and harvested at two color stages according to the USDA ripening classes of tomatoes [48,49]: (i) color stage 3 "Turning", over 10% but not more than 30% red, pink, or orange-yellow and (ii) color stage 5 "Light red", over 60% but not more than 90% red (Figure 1). The attributes evaluated were as follows. The content of soluble solids ( • Brix) was determined with a Smart-1 refractometer (ATAGO, Tokyo, Japen) to ± 0.1 • Brix. Fruit pH and titratable acidity (TA) were determined using a pH electrode and an 862 Compac Titrosampler (Metrohm, Herisau, Switzerland) with accuracy ± 0.01. The TA was expressed as % citric acid after titration with NaOH 0.1 N [50].
1, x FOR PEER REVIEW 8 of 28 Figure 1. Information about the tomato fruit quality attributes assessed.

Statistical Analysis
Student's t-test was used to determine statistically significant effects and differences among treatments (two levels: fertigation and water) for each of the variables evaluated in the study (physical and chemical variables of soil, soil fungi variables, tomato crop production, and fruit quality variables) for both years. Previously, normality and homoscedasticity were tested using the Shapiro-Wilk and Levene tests, respectively. The Kruskal-Wallis one-way non-parametric test (p = 0.05) was performed when normality or homo-

Cost of Fertilizers Applied in Fertigated Plots
To estimate the amount of fertilizers used during the tomato crop cycles and determine the economic cost, the fertilization plan applied to the "fertigation" treatment was documented for the two crop cycles considered in this study. The cost per unit of area was calculated based on the consumption of each fertilizer and the price paid to the supplying company for the fertilizers.

Statistical Analysis
Student's t-test was used to determine statistically significant effects and differences among treatments (two levels: fertigation and water) for each of the variables evaluated in the study (physical and chemical variables of soil, soil fungi variables, tomato crop production, and fruit quality variables) for both years. Previously, normality and homoscedasticity were tested using the Shapiro-Wilk and Levene tests, respectively. The Kruskal-Wallis one-way non-parametric test (p = 0.05) was performed when normality or homoscedasticity of data was not evident (p < 0.05 Shapiro-Wilk or Levene tests, respectively). Arcsine square root transformation was applied to percentages before analyses. These statistical analyses were carried out using the statistical software package Statgraphics Centurion XIX (Statgraphics Technologies, Inc., The Plains, VA, USA) for Windows (Microsoft Corporation, Redmond, WA, USA).

Physical-Chemical Characteristics of Soil Samples
Out of all of the physical and chemical variables evaluated, only differences in Mg +2 concentration (p < 0.01, Table 3) were reported at the end of Season 2, being higher in the soil of plots irrigated only with water during the crop season compared to in the soil of those that were fertigated. Overall, the organic fertigation program did not increase any of the physical and chemical variables of soil samples. Table 3. Soil physical and chemical variables at the beginning of the study (start of Season 1), and at the end of the two growing seasons, in plots irrigated with water and through the fertigation program during the tomato crop cycles. Values (mean ± standard deviation; n = 3).

Effects of the Treatments on Soil Fungal Diversity Alpha Diversity of Soil Fungal Community
According to the plate culture evaluation method (culturable fungal microbiota), the initial state of the soil before starting the treatments showed an average CFU of 1800 (standard deviation, s.d.: 2779) ( Table 4). By the end of Season 1, the average number of CFU was 37,333 and 24,867 in the "fertigation" and "water" soil plots, respectively, while at the end of Season 2, the CFU numbers were 6667 and 10,467, respectively. In all cases, no statistical differences (p > 0.05) were detected between the two groups. There were also no differences in the number of identified genera at the end of both seasons, with average values ranging from 5.7 to 6.3, while at the beginning of the study, an average value of 1.7 (s.d.: 1.5) was identified. Additionally, there was no evidence that the fertigation treatment, which used fertilizers listed in Commission Implementing Regulation (EU) 2021/1165, had an impact on the alpha diversity indices of soil culturable fungi (Table 4). Table 4. Alpha diversity indices of soil fungal community at the beginning of the study (start of Season 1), and at the end of the two crop seasons, in plots irrigated with water and through the organic fertigation program during the tomato crop cycles. Values (mean ± standard deviation; n = 3) are shown for the two analytical methods performed in the study of soil fungal microbiota (plate culture method and DNA metabarcoding). According to the molecular evaluation method based on DNA metabarcoding, at the beginning of the study, the average sequencing depth for the ITS library was 49,649 valid reads (s.d.: 3307) after quality trimming (Table 4). By the end of Season 1, the average number of valid reads was 54,212 and 49,679 in the "fertigation" and "water" soil plots, respectively, while at the end of Season 2, the number of valid reads was 90,645 and 93,058, respectively. In all cases, no statistical differences (p > 0.05) were detected between the two groups. There were also no differences in the number of OTUs at the end of both seasons, with average values ranging from 22.0 to 56.6, while at the beginning of the study, an average of 41.3 (s.d.: 8.3) was identified. Additionally, there was no evidence that the fertigation treatment, which used fertilizers listed in Commission Implementing Regulation (EU) 2021/1165, had an impact on the alpha diversity of soil fungi DNA (Table 4).

Comparison of Soil Fungal Community Structure between the Two Treatments Beta Diversity of Soil Fungal Community
The results of PCoA showed that the fungal communities in the samples of the two treatments had relatively discrete distribution at the end of the two seasons regardless of the analytical method ( Figure 2), with relatively large distances between the samples of the two treatments. This is more evident with data from DNA metabarcoding analysis and mainly at the end of Season 1 ( Figure 2B), since at the end of Season 2 ( Figure 2D), two samples, one from each treatment, are located relatively close to each other.
According to SIMPER analysis, the fungal community composition identified at the genera level through the plate culture evaluation method (culturable fungal microbiota) in the fertigation and water treatments showed an average dissimilarity of 63.98 and 61.12 at the end of Season 1 and Season 2, respectively ( Table 5). The contribution percentage of each culturable fungus genus to the dissimilarity between two treatments (water and fertigation), as well as the cumulative contribution percentage of that genus, is shown in Table 5. The table displays the 13 and 12 genera that were identified at the end of Season 1 and Season 2, respectively, along with their average abundance (note that data were square root transformed) in both the fertigation and water treatments. Although their contribution is different in each season, Acremonium, Penicillium, Mortierella, and Aspergillus are the top four genera that contribute the most to the dissimilarity between water and fertigation at the end of both seasons. These top four genera contributed to over 65% of the dissimilarity between the two treatments.

Comparison of Soil Fungal Community Structure between the Two Treatments Beta Diversity of Soil Fungal Community
The results of PCoA showed that the fungal communities in the samples of the two treatments had relatively discrete distribution at the end of the two seasons regardless of the analytical method ( Figure 2), with relatively large distances between the samples of the two treatments. This is more evident with data from DNA metabarcoding analysis and mainly at the end of Season 1 ( Figure 2B), since at the end of Season 2 ( Figure 2D), two samples, one from each treatment, are located relatively close to each other.   Table 5. Contribution (Contrib%) and cumulative contribution (Cum%) of the SIMPER analysis of all of the fungal genera (culturable fungi) identified through the plate culture method that contributed to the differences in the microbial community between soils irrigated with water and through the fertigation program during the tomato crop cycles at the end of the two crop seasons. On the other hand, the fungal community composition identified at the OTU level through DNA metabarcoding in the fertigation and water treatments showed an average dissimilarity of 79.46 and 76.43 at the end of Season 1 and Season 2, respectively ( Table 6). The 30 main OTUs that contributed the most to the differences in the microbial community between soils irrigated with water and through the fertigation program during the tomato crop cycles at the end of the two crop seasons contributed to over 92% and 93% of the dissimilarity between the two treatments at the end of Season 1 and Season 2, respectively. At the end of Season 1, the values range from 0.65% to 14.88%, with the most dissimilar OTU being from the phylum Ascomycota. The second and third most dissimilar OTUs are from the phylum Rozellomycota and the genus Plectosphaerella, respectively. It is worth noting that these three OTUs together contribute to more than 42% of the dissimilarity observed between the two groups. This is followed by the phylum Basidiomycota and the genus Aspergillus, both contributing more than 4%. At the end of Season 2, the results show that again the phylum Ascomycota had the highest average abundance in both treatments, and the highest contribution to the differences between the two treatments (14.18%). The second and third highest contributors were the fungal genera Aspergillus (7.81%) and Thermomyces (6.66%), respectively. This is followed by the phylum Rozellomycota and the genus Mortierella, with a contribution of more than 6% in both cases. These five OTUs contribute to more than 41% of the dissimilarity observed between the two groups at the end of Season 2. Table 6. Contribution (Contrib%) and cumulative contribution (Cum%) of the SIMPER analysis of the 30 main OTUs identified through DNA metabarcoding that contributed the most to the differences in the microbial community between soils irrigated with water and through the fertigation program during the tomato crop cycles at the end of the two crop seasons. However, PERMANOVA analysis showed that the composition of the fungal communities did not differ significantly between the treatments at the end of the two seasons, neither for cultivable fungal communities

Analysis of Differences in Fungal Genera Dominance
The relative abundance of the ten most abundant fungal genera detected in soil through the two analytical methods (plate culture method and DNA metabarcoding) at the end of the two crop seasons in plots irrigated with water and through the organic fertigation program during the tomato crop cycles is shown in Figure 3. Significant differences were only exhibited for two cultivable fungal genera at the end of Season 2 ( Figure 3C), while in no case were differences exhibited in the relative abundance of fungal genera identified through DNA metabarcoding ( Figure 3B,C). The cultivable fungal genera that showed differences were Mortierella, which exhibited higher relative abundance in the "water" treatment soils, and Mucor, which was more abundant in the "fertigation" treatment soils. These differences do not appear to be linked to the inoculum of cultivable fungal genera incorporated through irrigation in both treatments ( Figure S1).

Analysis of Differences in Fungal Genera Dominance
The relative abundance of the ten most abundant fungal genera detected in soil through the two analytical methods (plate culture method and DNA metabarcoding) at the end of the two crop seasons in plots irrigated with water and through the organic fertigation program during the tomato crop cycles is shown in Figure 3. Significant differences were only exhibited for two cultivable fungal genera at the end of Season 2 ( Figure  3C), while in no case were differences exhibited in the relative abundance of fungal genera identified through DNA metabarcoding ( Figure 3B,C). The cultivable fungal genera that showed differences were Mortierella, which exhibited higher relative abundance in the "water" treatment soils, and Mucor, which was more abundant in the "fertigation" treatment soils. These differences do not appear to be linked to the inoculum of cultivable fungal genera incorporated through irrigation in both treatments ( Figure S1).

Taxonomic Groups Associated to Fertilization Practices
LEfSe analysis identified biomarkers that caused significant differences between the two treatments at the end of the two seasons, although they did not match between the two analytical methods (Figure 4). Thus, at the end of Season 1, only two biomarkers were identified in the soil fungal community analyzed by means of DNA metabarcoding in the "fertigation" treatment ( Figure 4B) with taxa at different taxonomic levels (the class Leotiomycetes and the family Myxotrichaceae). At the end of Season 1, no biomarkers were detected in the cultivable soil fungi genera ( Figure 4A). However, contrary to the previous results, at the end of Season 2, only two biomarkers were detected in the cultivable soil fungi genera ( Figure 4C), but no biomarkers were detected in the soil fungal community analyzed by means of DNA metabarcoding ( Figure 4D). In this case, the genus Mortierella was detected as a biomarker of the "water" treatment, while the phylum Ascomycota and the genus Mucor was detected as a biomarker of the "fertigation" treatment. All identified biomarkers had LDA scores greater than 5.

Taxonomic Groups Associated to Fertilization Practices
LEfSe analysis identified biomarkers that caused significant differences between the two treatments at the end of the two seasons, although they did not match between the two analytical methods (Figure 4). Thus, at the end of Season 1, only two biomarkers were identified in the soil fungal community analyzed by means of DNA metabarcoding in the "fertigation" treatment ( Figure 4B) with taxa at different taxonomic levels (the class Leotiomycetes and the family Myxotrichaceae). At the end of Season 1, no biomarkers were detected in the cultivable soil fungi genera ( Figure 4A). However, contrary to the previous results, at the end of Season 2, only two biomarkers were detected in the cultivable soil fungi genera ( Figure 4C), but no biomarkers were detected in the soil fungal community analyzed by means of DNA metabarcoding ( Figure 4D). In this case, the genus Mortierella was detected as a biomarker of the "water" treatment, while the phylum Ascomycota and the genus Mucor was detected as a biomarker of the "fertigation" treatment. All identified biomarkers had LDA scores greater than 5.

Nutritional and Functional Analyses of Fungal Communities identified through DNA Metabarcoding
FUNGuild was used to predict the nutritional and functional groups of the fungal communities with different treatments that were identified through DNA metabarcoding at the end of the two seasons. At the end of Season 1, a total of 97 ASVs were identified, of which 43.3% did not receive functional assignment and 6.2% received "Possible" confidence scores, so the results are presented for 50.5% of the identified ASVs. This represents 53.8% (s.d. 16.8) and 39.7% (s.d. 10.5) of the total "fertigation" and "water" sequencing reads, respectively. At the end of Season 2, a total of 197 ASVs were identified, of which 48.2% did not receive functional assignment and 12.7% received "Possible" confidence scores, so the results are presented for 39.1% of the identified ASVs. This represents 28.6% (s.d. 35.9) and 33.2% (s.d. 27.0) of the total "fertigation" and "water" sequencing reads, respectively.
The results show that the fungi community was divided into seven trophic mode groups. The fungal community was screened for 39 and 47 identifiable species, at the end of Season 1 and Season 2, respectively, to which FUNGuild assigned confidence scores of "Probable" and "Highly Probable". They included saprotrophs (46.2% and 57.4% at the end of Season 1 and Season 2, respectively), pathotrophs (23.1% and 17.0%), pathotrophsymbiotrophs (10.3% and 2.1%), pathotroph-saprotroph-symbiotrophs (10.3% and 2.1%), saprotroph-symbiotrophs (5.1% and 12.8%), pathotroph-saprotrophs (2.6% and 6.4%), and symbiotrophs (5.1% and 2.1%), representing the general abundance of each nutrition method in the identified community ( Figure 5). Among them, pathotrophs, saprotrophs, and symbiotrophs were the major components, so the species with various trophic modes assigned were considered for all of them. The pathotrophic mode groups primarily consisted of plant pathogens (not representing a specific specie) and animal pathogens, which showed different relative abundances depending on the season, but no differences were identified between treatments ( Figure 6(A1,B1)).
communities with different treatments that were identified through DNA metabarcoding at the end of the two seasons. At the end of Season 1, a total of 97 ASVs were identified, of which 43.3% did not receive functional assignment and 6.2% received "Possible" confidence scores, so the results are presented for 50.5% of the identified ASVs. This represents 53.8% (s.d. 16.8) and 39.7% (s.d. 10.5) of the total "fertigation" and "water" sequencing reads, respectively. At the end of Season 2, a total of 197 ASVs were identified, of which 48.2% did not receive functional assignment and 12.7% received "Possible" confidence scores, so the results are presented for 39.1% of the identified ASVs. This represents 28.6% (s.d. 35.9) and 33.2% (s.d. 27.0) of the total "fertigation" and "water" sequencing reads, respectively.
The results show that the fungi community was divided into seven trophic mode groups. The fungal community was screened for 39 and 47 identifiable species, at the end of Season 1 and Season 2, respectively, to which FUNGuild assigned confidence scores of "Probable" and "Highly Probable". They included saprotrophs (46.2% and 57.4% at the end of Season 1 and Season 2, respectively), pathotrophs (23.1% and 17.0%), pathotrophsymbiotrophs (10.3% and 2.1%), pathotroph-saprotroph-symbiotrophs (10.3% and 2.1%), saprotroph-symbiotrophs (5.1% and 12.8%), pathotroph-saprotrophs (2.6% and 6.4%), and symbiotrophs (5.1% and 2.1%), representing the general abundance of each nutrition method in the identified community ( Figure 5). Among them, pathotrophs, saprotrophs, and symbiotrophs were the major components, so the species with various trophic modes assigned were considered for all of them. The pathotrophic mode groups primarily consisted of plant pathogens (not representing a specific specie) and animal pathogens, which showed different relative abundances depending on the season, but no differences were identified between treatments ( Figure 6A1,B1).
In terms of symbiotrophic mode groups, at the end of Season 1, endophytes, epiphytes, arbuscular mycorrhizae, and ectomycorrhizae were identified, while endophytes and ericoid mycorrhizae were identified at the end of Season 2. However, no significant differences in the relative abundances were found between treatments ( Figure 6A2,B2). The mode groups of saprotrophs included undefined saprotrophs, dung saprotrophs, litter saprotrophs, wood saprotrophs, and nematophagous fungi. The relative abundance of undefined saprotrophs and nematophagous fungi was higher (p < 0.05) in the "fertigation" treatment compared to the "water" treatment at the end of Season 1 and Season 2, respectively. Nevertheless, the relative abundance of dung saprotrophs was higher (p < 0.05) in the "water" treatment at the end of Season 2 ( Figure 6A3,B3).   In terms of symbiotrophic mode groups, at the end of Season 1, endophytes, epiphytes, arbuscular mycorrhizae, and ectomycorrhizae were identified, while endophytes and ericoid mycorrhizae were identified at the end of Season 2. However, no significant differences in the relative abundances were found between treatments (Figure 6(A2,B2)). The mode groups of saprotrophs included undefined saprotrophs, dung saprotrophs, litter saprotrophs, wood saprotrophs, and nematophagous fungi. The relative abundance of undefined saprotrophs and nematophagous fungi was higher (p < 0.05) in the "fertigation" treatment compared to the "water" treatment at the end of Season 1 and Season 2, respectively. Nevertheless, the relative abundance of dung saprotrophs was higher (p < 0.05) in the "water" treatment at the end of Season 2 ( Figure 6(A3,B3)).

Tomato Crop Production
There were no significant differences between the two treatments in any of the harvests nor in the accumulated marketable tomato production in either of the two years of study (Figure 7). In plots that only received basal dressing with fresh sheep manure (i.e., "Water" treatment), the total marketable tomato production was 13.10 and 15.27 kg m −2 in Season 1 and Season 2, respectively. In plots that also received nutrient inputs through fertigation (i.e., "fertigation" treatment), the production was 14.49 and 16.09 kg m −2 in Season 1 and Season 2, respectively. The higher production obtained in the second year was mainly due to a longer crop cycle than the first.

Tomato Fruit Quality
Only the Brix degrees of the tomatoes at the "Turning" color stage showed significant differences between treatments in one of the sampling times, being higher in those from the fertigation treatment compared to those from the water treatment (Table 7). In this case, the higher EC of the irrigation water from the fertigation treatment (Table S1) may have influenced the Brix degrees of the fruits evaluated. None of the rest of the tomato fruit quality variables tested showed significant differences between treatments.

Cost of Fertilizers Applied in the Two Treatments
The cost of fertilizers applied during the crop cycles in the fertigated plots (Table 2)

Tomato Fruit Quality
Only the Brix degrees of the tomatoes at the "Turning" color stage showed significant differences between treatments in one of the sampling times, being higher in those from the fertigation treatment compared to those from the water treatment (Table 7). In this case, the higher EC of the irrigation water from the fertigation treatment (Table S1) may have influenced the Brix degrees of the fruits evaluated. None of the rest of the tomato fruit quality variables tested showed significant differences between treatments. Table 7. Tomato fruit quality variables in plots irrigated with water and through the organic fertigation program during the tomato crop cycles. Assessments performed at two different times during each crop cycle (i.e., Season 1 and Season 2). Values (mean ± standard deviation; n = 3 in Season 1 and n = 5 in Season 2).

Cost of Fertilizers Applied in the Two Treatments
The cost of fertilizers applied during the crop cycles in the fertigated plots (Table 2), along with the corresponding cost of the sheep manure incorporated prior to the start of the crop seasons (i.e., EUR 4200.0 ha −1 , not including labor), amounted to EUR 17,768.7 and EUR 24,245.7 ha −1 (including VAT) in Season 1 and Season 2, respectively. The incorporation of sheep manure as a basal dressing without subsequent fertigation resulted in a cost savings of 76.4-82.7% compared to the cost of organic fertigated plots.

Discussion
The sustainability of agricultural systems relies on the type and intensity of the cultural practices carried out. Basal dressing fertilization and continuous fertigation during the crop season are common practices for vegetable growers in intensive systems, including greenhouse organic farming. Both types of fertilization, using organic matter, can influence crop profitability and soil fertility, including its impact on soil fungal communities. In the present study, the incidence and impact of continued fertigation with fertilizers listed in the Commission Implementing Regulation (EU) 2021/1165 on products and substances for use in organic production on the main physical and chemical variables of soil and its fungal communities, as well as on the production and quality of greenhouse tomato fruit, were evaluated for two consecutive years after a basal fresh sheep manure application. To obtain a more detailed understanding of the impact on soil fungal communities, the study was conducted using both plate culture (culturable fungi) and DNA metabarcoding methods. Additionally, a study of the cost of fertilizers applied during the two crop cycles in the fertigated plots was conducted. At the end of the two crop cycles, none of the physicochemical variables were increased in the fertigated plots compared to those where soils only received water during the crop cycles and the only added fertilizer was the basal fresh sheep manure. Concerning soil fungal communities, there were no significant differences between the treatments in the alpha diversity indices of the soil fungal community or in the fungal community composition, and differences were only found regarding the identification of biomarkers. The fertilizers added through fertigation did not consistently improve the tomato yield and fruit quality in any of the two years studied, despite their high cost.
In terms of soil physicochemical properties, in this study, we only detected differences in the exchangeable magnesium (Mg +2 ) content at the end of the second year. It was found in higher concentrations in the soils irrigated with water during the tomato crop cycles, which had been previously amended with fresh sheep manure. These differences could be influenced by the extraction of this fraction of the mineral element from the soil by the metabolic products of a living organism. For instance, the importance of Mg in the metabolism of fungi in the Aspergillus genus (e.g., A. niger) has been widely reported, and several authors have successfully demonstrated the ability of the fungi to extract Mg from the soil for its own growth [51][52][53]. In this regard, a technique based on the correlation between the weight of mycelium of the fungus and the available magnesium content of the soil has been successfully used [52,53]. In our study, although we did not detect significant differences in the relative abundance of Aspergillus in soils due to the treatments, this fungal genus was identified through the two analytical methods as one of the main OTUs and/or genera that contributed the most to the differences in the microbial community between soils irrigated with water and those fertigated during the tomato crop cycles at the end of the two crop seasons. The average abundance of Aspergillus was higher in the fertigated soils; thus, the extraction of Mg +2 could also be higher, which could be the reason why the content of Mg +2 in the soil was lower. This is only a possible hypothesis, and it should be noted that while the role of soil fungi, such as mycorrhizae, in N and P uptake and balance has been intensively studied, the effects on metal nutrients such as Mg are less clear and more inconsistent in cultivated soils for vegetable crops [54].
The availability of an adequate nutrient supply is crucial for ensuring optimal growth and productivity in plants [55]. Inadequate synchronization between nutrient availability and plant uptake has been recognized as a significant limitation to the productivity of organic plant production systems [56,57]. In our study, the basal fresh sheep manure application was sufficient for the correct development of a greenhouse tomato crop, without showing any deficiencies due to a lack of nutrient availability during the two-year study period. In this sense, the incorporated manure could have promoted the prevalence of soil microbial communities, including the fungal community, capable of meeting the crop's needs throughout the growing season. In this regard, while standardized methods for evaluating soil physicochemical parameters have been in place for several decades, the same cannot be said for the study of fungal microbiota in soil. Continuous efforts have been devoted to the standardization of the numerous methods that have been developed over the last few decades to enable a better understanding of the structure and ecology of soil microbial communities and their roles in soil functioning [58][59][60]. In this respect, the emergence of high-throughput sequencing (including DNA metabarcoding) in recent years has allowed scientists worldwide to identify soil fungi. The majority of published studies utilize Illumina MiSeq to characterize fungi, with a focus on targeting ITS regions, especially ITS2, which is also considered the optimal technique for comparative studies [10]. On the other hand, the plate culture method makes it possible to detect and identify culturable fungal species. This method is increasingly less widely used mainly because it is laborious, time-consuming, and requires the expertise of taxonomists for accurate identification. In addition, the microbial diversity of soil microbiome is highly underestimated, as it is estimated that only 1-5% of total soil microorganisms have been cultured by current methodologies [61,62]. However, while through high-throughput DNA sequencing both the live (including active and inactive fractions) and dead fungi are identified, the culturable fungi identified from a soil sample are considered active or potentially active fungal groups. The latter exhibit specific functions or trophic requirements with a certain degree of certainty, allowing for studies on their functionality [61]. In our study, we used both methods for studying the soil mycobiome. Using the plate culture method, a total of 13 and 12 fungal genera were identified at the end of Seasons 1 and 2, respectively, resulting in a total of 16 different fungal genera identified throughout the entire study. Through DNA metabarcoding, 40 and 50 fungal genera were identified at the end of Seasons 1 and 2, respectively, resulting in a total of 74 different fungal genera identified throughout the entire study. Although the results we obtained must be considered complementary, for the purpose of comparing the evaluated treatments in relation to soil fungal communities, both techniques led us to the same conclusions. Thus, the analysis of fungi alpha diversity showed no differences in any of the calculated indices due to the treatments, regardless of the evaluation method for soil fungi assessment. Similarly, the PERMANOVA analysis did not show significant differences in the composition of the fungal communities (fungal beta diversity), although the SIMPER analyses showed high dissimilarity values that were reflected in the PCoA plots. Nevertheless, the LEfSe analysis identified biomarkers that caused significant differences between the two treatments at the end of the two seasons, although they did not match between the two analytical methods. The class Leotiomycetes and the family Myxotrichaceae (included in the class Leotiomycetes) were identified as biomarkers in the soil fungal community analyzed by means of DNA metabarcoding of the "fertigation" treatment at the end of Season 1. The genus Mortierella was detected as a biomarker of the "water" treatment, and Mucor was identified as a biomarker of the "fertigation" treatment in the culturable soil fungi at the end of Season 2. Leotiomycetes (inoperculate discomycetes) is a class of Ascomycota which includes powdery mildews that are classified in the order Erysiphales and the family Erysiphaceae [63]. Likewise, it includes the family Myxotrichaceae, a cellulolytic group that encompasses several species. In this regard, Myxotrichum deflexum was exclusively identified in the 'fertigation' treatment at the end of Season 1. Gymnascella nodulosa, on the other hand, was identified in both treatments but only at the end of Season 2. Additionally, Oidiodendron periconioides and Oidiodendron echinulatum, which are Oidiodendron anamorphs, were exclusively identified in the 'fertigation' treatment at the end of Season 1. Moreover, Pseudogymnoascus roseus was exclusively identified in the 'fertigation' treatment at the end of Season 1. Lastly, Gymnascella dankaliensis, a Gymnoascus anamorph, was identified in both treatments but only at the end of Season 1. These are fungal species that have been shown to be able to form typical ericoid mycorrhizas [64]. Ericoid mycorrhizal fungi have saprotrophic and biotrophic lifestyles and are able to biodegrade SOM and establish symbiotic relationships with plants in the family Ericaceae [65]. Therefore, these types of mycorrhizae are not classified as beneficial for tomato plants, which belong to the Solanaceae family. Furthermore, this result was not consistent after two years of treatments (i.e., end of Season 2). Meanwhile, Mortierella is a fungal genus known for its ability to decompose organic matter and participate in nutrient cycles in the soil, and some species of fungi in the Mortierella genus demonstrate the ability to solubilize phosphorus (P) [66]. Recently, the work of [67] revealed that substituting manure for mineral P fertilizer significantly increases P availability in greenhouse soil. The authors suggested that this could largely be attributed to the enrichment of P-mineralizing microbes and indicated that manure is a feasible substitute for mineral P fertilizers in greenhouse farming. In addition, species from the Mortierella genus have also been recognized as plant-growth-promoting fungi [68] and as potential antagonistic agents against various plant pathogens [69][70][71]. On the other hand, the Mucor genus is a fungal genus classified in the Mucorales, the most prominent order of zygospore-forming fungi that formerly constituted the Zygomycota, a phylum which is not currently accepted due to polyphyly. While Mucor is primarily a saprobic fungal genus, some species have been identified as plant endophytes, and others have been found to be opportunistic pathogens of animals and humans. It is worth noting that no species of Mucor have been reported to have mutualistic associations with plants [72].
Our results also show a significant presence of fungi from other fungal genera, such as Acremonium, Aspergillus, and Penicillium, which were identified through the two analytical methods and were present at all sampling time points regardless of the treatment. These fungi have been recognized as phosphorus solubilizers [73][74][75][76][77] as well as plant-growthpromoting fungi [76,[78][79][80]. In order to obtain additional information, in the present study, the FUNGuild database was used to predict the nutritional and functional groups of fungal communities identified through DNA metabarcoding. No differences were identified between treatments in both pathotrophic and symbiotrophic modes, but differences were observed in the saprotrophic mode. Specifically, the relative abundance of undefined saprotrophs and nematophagous fungi was significantly higher in the "fertigation" treatment compared to the "water" treatment at the end of Season 1 and Season 2, respectively. However, the relative abundance of dung saprotrophs was higher in the "water" treatment at the end of Season 2. When interpreting these results, it is crucial to exercise caution due to their limited consistency, as they were not consistently observed throughout all sampling time points. Additionally, their low representativeness should be taken into account, considering the relatively small proportion of sequencing reads they encompass. Additional in-depth studies on soil fungal functional groups are also needed to further investigate the function of the soil fungal community as a function of the type of sustainable fertilization and to obtain a more comprehensive understanding of the role of fungi in soil agroecosystems. On the other hand, an issue that caught our attention is that while Exophiala appeared as the fungal genus with the highest relative abundance both in the water and in the fertilizer solution applied in the plots, this yeast was not identified in the soils by any of the evaluation methods in any case. Soil microbial communities are highly complex and diverse and are typically resistant to invasion by applied microorganisms [81,82].
Numerous studies have evaluated the impact of fertilizers on soil fertility, including the soil fungal community. In recent years, there have been several studies comparing the impacts of conventional fertilizers (i.e., chemical synthetic fertilizers) versus organic fertilizers [20,[83][84][85][86][87]. There are also studies comparing different organic fertilizers applied as basal amendments before crop installation, organic fertilizers applied through fertigation, or both in combination [21,86], and even studies comparing combinations of chemical fertilizers with organic fertilizers [88][89][90][91][92]. According to these studies, these practices generally have a significant impact on soil fungal communities. For instance, several studies have demonstrated that the use of inorganic fertilizers leads to a decrease in both fungal diversity and biomass [84,93]. Conversely, the application of organic fertilizers can alter the composition and abundance of fungal communities, potentially resulting in an increase in fungal diversity in some cases [93,94]. However, there have been reports of an increase in fungal diversity when inorganic fertilizers were applied, both with and without basal manure application [86,87,90]. In addition, the organic fertilizer type and composition seem to be relevant for determining the diversity of soil fungi communities [91]. Thus, a recent meta-analysis of 37 studies on the effect of organic and mineral fertilizers on soil microbial diversity suggests that a large amount of further research is required to fully understand the influence of fertilizer regimes on microbial diversity and ecosystem function [83]. Therefore, the extent of the impact on soil fungal communities varies, and it remains uncertain whether results obtained from studies conducted on specific soils and agroecosystems can be used to infer trends at a broader scale. Furthermore, irrigation and fertilization levels have also shown significant effects on fungal community diversity [95].
To the best of our knowledge, no research to date has examined the impact of organic fertilizers incorporated through fertigation following basal fresh sheep manure application on soil fertility, including the soil fungal community, and on tomato crop yield in intensive greenhouse vegetable production systems. This practice has a significant impact on the sustainability and profitability of organic vegetable production systems in greenhouses, where improving soil health and reducing agricultural inputs are essential. Thus, this study is equally relevant to farmers who decide to initiate the necessary conversion process to certify their greenhouse area as "organic" according to Regulation (EU) No 2018/848 on organic production. The study presented here was conducted for two consecutive years, the required timeframe for horticultural crops. According to the obtained results, the continued organic fertigation after basal manure application does not impact soil fungal communities, tomato yield, or soil fertility. Therefore, incorporating manure as a basal fertilizer in intensive tomato crops is proposed as a crucial practice to achieve more sustainable and profitable greenhouse agriculture. Consequently, it emerges as a key practice to overcome the primary challenge of shifting soil modeling from a chemical-based approach to an agroecological approach [4], and thus enables most farmers to make this transition.

Data Availability Statement:
The data presented in this study are available within the article and Supplementary Materials.