Plants were grown in a repeated measures design in a climate-controlled greenhouse with mean daytime temperature of 28°C and night time of 22°C at the University of Hawai'i at Mānoa. The Lahaina series (Very-fine, kaolinitic, isohyperthermic Rhodic Eutrustox, pH 7.9) was used as the growth medium, chosen for its naturally low fertility but agricultural importance in Hawai'i. This Oxisol is comparable to other Oxisols in tropical regions of the world (38), although the alkalinity, prior use for sugarcane production, and anecdotal evidence suggests it had been previously limed. The soil was collected from the 0 to 30 cm depth from an uncultivated field directly adjacent to an organically managed agricultural system in Waialua, HI (21.555 °N, −158.117 °W). This location had a homogenous vegetation cover of Guinea grass (Megathyrsus maximus) that helped to limit the amount of variation in the microbial community due to plant diversity (39, 40). Collected soil was air-dried, homogenized, and sieved to 2 mm. Plastic pots (7.6 L) with drainage holes, covered with a fiberglass screen mesh, were filled with 1 kg of acid-washed sand, followed by 3.5 kg of the homogenized soil mixed with each fertilizer treatment (urea, compost, and a no fertilizer control) at an application rate of 100 kg/ha plant-available nitrogen, which was determined as the sum of inorganic N content and estimated mineralizable N. We used a 3% total N EcoSan compost (3–3.7–1.8, N-P-K), a product of aerobic, thermophilic composting process from human feces and sugarcane bagasse (41, 42). An estimated N mineralization rate of 7% was used for compost based on a recent review of N mineralization in composted biosolids (43). As such, we added 215.4 g of compost (dry weight) to each pot (equivalent to 47.6 metric tons/ha). We used urea (46-0-0, N-P-K) as a contrasting synthetic fertilizer. All N in urea was considered plant-available. Fertilizers were applied only once, immediately prior to the first planting. Each treatment was replicated three times (n = 3) for a total of nine pots (N = 9). Pots were arranged in a randomized control block design to account for possible environmental heterogeneity in the greenhouse. Because the soil had been dried, pots were watered to field capacity and allowed to equilibrate for 2 days prior to planting. Nine organic radish seeds (Raphanus raphanistrum subsp. sativus ‘Cherry Belle', Burpee) were planted in each pot and grown for seven days; they were then thinned to three plants per pot. We chose radishes because of their fast growth from sowing to harvesting that allows for multiple growing cycles within a short period of time. Plants were hand-irrigated with deionized water (100–300 ml/day depending on the weather and stage of growth) to maintain approximate field capacity. Each cycle of radishes was grown for 36 days from seed to maturity, at which point they were harvested. The plants were uprooted, washed with deionized water, and divided at the crown for below and aboveground biomass. All three plants per pot were summed for a total biomass per replicate. Aboveground (shoots) and belowground (taproots) samples were put into paper bags and dried in an oven at 65°C for 2 and 7 days, respectively, for dry weight measurements. Large radish roots were sliced in half to facilitate the drying process. At the end of each harvest, pots were watered and allowed to equilibrate for 2 days prior to the sowing of the next cycle. These plantings and harvesting cycles were repeated for a total of 6 crop cycles, or 227 days, beginning on August 22, 2017.

Soil cores were collected immediately after treatment application (C0) and immediately prior to harvest on the 36th day for each of the six crop cycles (C1–C6). From each pot, three cores (1 cm diameter x 10 cm deep, total 63 samples) were taken about 2 cm from each radish taproot. After collection, these cores were composited and homogenized in a clean plastic bag. A subsample was put into a 10 mL transfer tube and stored in a−20°C freezer for microbial analysis. The remaining samples were air-dried and analyzed for pH using a slurry method with a 2 to 1 ratio of deionized water:soil. Soil samples collected after the sixth crop cycle (C6) were analyzed for total carbon, total nitrogen, and exchangeable cations (Table 1). Soil carbon and nitrogen concentrations were analyzed by combustion on an elemental analyzer (Costech 4100 Elemental Analyzer) at the University of Hawai‘i at Hilo Analytical Laboratory. Exchangeable calcium, magnesium, potassium, and sodium were extracted from soils using the ammonium acetate method buffered at pH 7.0 (44) and analyzed on a Thermo iCAP DUO 7400 ICP-OES. Cation exchange capacity was calculated as the sum of base cations.

DNA was extracted from 0.25 g of frozen soil using the DNeasy PowerSoil kit (QIAGEN, Germany) following the manufacturer's standard protocol. Fungal and bacterial mock communities served as positive controls (45). Primer design and sample barcoding followed a two-step amplification, dual barcoding system using a combination of a P5/P7 Illumina Adapter, 8 bp barcodes attached to each adapter, and partial P5/P7 overhangs (Supplementary Figure 1). For bacteria and archaea, the 16S rRNA gene (V4 region) was targeted using updated Earth Microbiome Project primer pairs 515F (46) and 806RB (47). For fungi, the ITS1 region was targeted using primer pairs ITS1F (48) and ITS2 (49). High-fidelity, hot-start polymerase mastermixes were used for the PCR reaction with iProof (Bio-Rad Laboratories, USA) for 16S and Phusion (ThermoFisher Scientific Inc., USA) for ITS. Internal testing showed that these polymerases amplified each respective gene more efficiently. In the first PCR, the targeted loci (16S, ITS) were amplified under the following thermocycling conditions using 3-5 ng of DNA: For 16S, polymerase activation at 98°C for 30 s, followed by 20 cycles of 98°C for 10 s, 55°C for 15 s, 72°C for 10 s, and a final extension at 72°C for 7 min; and for ITS, polymerase activation at 98°C for 30 s, followed by 20 cycles of 98°C for 10 s, 53°C for 15 s, 72°C for 10 s, and a final extension at 72°C for 7 min. PCR products were assessed using gel electrophoresis and successfully amplified samples were cleaned using AMPure SPRI beads (Beckman Coulter, USA). In the second PCR step, 1 μl of the cleaned products was used as templates and amplified using a second set of primers that included P7/P5 overhangs, barcodes, and Illumina adapters (Supplementary Figure 1) under the following thermocycling conditions: polymerase activation at 98°C for 30 s, followed by 14 cycles of 98°C for 10 s, 52°C for 15 s, 72°C for 10 s, and final extension at 72°C for 7 min. PCR products were assessed using gel electrophoresis, and cleaned using SPRI beads as above. Barcoded amplicon libraries were quantified using a Qubit 3 Fluorometer (ThermoFisher Scientific Inc., USA). Each sample for a gene library (16S, ITS) was combined at equimolar concentration, including negative controls according to Nguyen et al. (45). Gene libraries were then combined at a ratio of 3:7 (by mass) ITS to 16S, spiked with 11.75% PhiX, and sequenced together on a single lane of Illumina (MiSeq) PE250 at the University of California, Davis Genome Center. Sequence data was deposited in the Sequence Read Archive under BioProject #PRJNA551045.

Sequence data processing and quality control were performed using the QIIME2 v2018.11 workflow and available plugins (50). Raw sequences were demultiplexed and primers and adapters removed. For 16S sequences, ends of sequences with base quality of q < 25 were truncated, and reads were paired and denoised using the DADA2 plugin (51). DADA2 Amplicon Sequence Variants (ASVs) were used for further analyses. For ITS sequences, conservative regions (18S and 5.8S) that flanked the ITS1 gene were removed using ITSXpress (52), followed by pairing and denoising as described above. This step was essential to accurately classify ITS sequences as verified using our mock community data. Furthermore, DADA2 ASVs overinflated fungal mock community richness (multiple OTUs per species) so it was necessary to further cluster these OTUs into 97% similarity OTUs using open reference clustering via the VSEARCH plugin (53). This clustering was essential to accurately recover the expected mock community diversity in the dataset (45). Hereafter both 16S DADA2 ASVs and ITS 97% clustered OTUs will be referred to simply as OTUs. The Naïve Bayes Classifier was used to classify OTUs using the Greengenes Database specifically trimmed for the primer pairs 515F-806R (gg-13-8-99-515-806-nb, (54) for prokaryotic 16S, and the UNITE database version 7.2 (55) for fungal ITS. For 16S, sequences classified as “unassigned”, “mitochondria”, or “chloroplast” were removed. For ITS, only sequences that aligned with at least 70% of reference sequences were kept. Internal validation showed that sequences that do not currently match to 70% of the UNITE v7.2 database are typically non-fungal. Sequences labeled as “unassigned”, “rhizaria, “Protista”, and “Metazoa” were removed. Negative PCR controls showed a maximum occurrence of 3 OTUs for fungi and 6 OTUs for bacteria and archaea. As a quality control measure, these maximum occurrence numbers from the negative control were then used as minimum occurrence thresholds for the rest of the data. As a result, only OTUs that had more than 7 sequences were kept for the 16S dataset and only OTUs that had more than 4 sequences were kept for ITS. The 16S dataset was rarefied to 2831 sequences, and the ITS dataset was rarefied to 514 sequences. Sampling saturation was variable for datasets (Supplementary Figure 2). Trophic guild data for fungi was identified using FUNGuild (56).

Statistical analyses were conducted in R version 3.5.1 (57) and QIIME2 version 2018.11 (50). A repeated measures approach was used to compare changes in above and below plant biomass as well as microbial community richness and diversity as response variables to the different fertilizer treatments and control. A generalized least squares (GLS) model was chosen due to its flexibility when dealing with observations across time (58), which is common for repeated measures experiments with autocorrelated residuals (59, 60). The GLS model y_i= X _iβ+ε_i, ε_i~(V_i), was created using the “lme” function in the package nlme (58), where y_i denotes the outcome variable for the i-th group (our response), X_i is the design matrix for the fixed effects β (Treatment + Time; Treatment & Time interactions) and V_i, the variance-covariance matrix for the error terms. V_i was specified to CorAR1, which confined errors within each individual replicate (pot). For the biomass data, a box-cox transformation was performed to meet the assumption of normality and equal variance. An Analysis of Variance (ANOVA) was performed using the “anova” function in base R. Pairwise comparisons of least-squares means were calculated for sample groupings by both treatment and time. These were calculated using functions “lsmeans” and “cld” via the lsmeans and emmeans (61, 62) packages. Least-squared mean comparisons were averaged over the levels of treatment using the Satterthwaite method (63) at a confidence level of 0.95. For the biomass data only, the “glht” function in the multcomp package was used (64). Pairwise comparisons were corrected for multiple comparisons with the Tukey p-adjustment.

Sequence reads obtained through high-throughput sequencing are compositional (65). Hence, in our community composition analysis, we followed protocol as suggested for compositional datasets. The first steps of data wrangling were conducted in R using the vegan package (66). First, unrarefied abundance matrices were normalized using the function “decostand” followed by the Hellinger square root method (67). We found this was sufficient in balancing our data. Next, Euclidean distances were created using the function “vegdist” on the transformed data matrices. We chose Euclidean distances as our metric for dissimilarity because of their appropriateness in Eigenvector-type analyses such as Principal Components Analysis (68–70), which we used for visualization. Community composition was compared using the function “adonis2”, which conducted permutational analysis of variances (PERMANOVAs) on the Euclidean distances as grouped by treatment. These PERMANOVAs were performed at the first crop cycle (1) and then included subsequent crop cycle accumulations (crop cycles 1–2, 1–3, 1–4, 1–5, and 1–6). This allowed us to compare how samples grouped in “species-space” at different points throughout the study. Comparisons were conducted using both an un-nested approach that looked at interactive effects and a nested approach that limited errors within “strata” of time within each crop cycle. Pairwise comparisons of treatment groupings were performed using the Wilks statistic with FDR corrections for multiple testing (59) via the function “pairwise.perm.manova” from the package RVAideMemoire version 0.9–72 (71). Results were visualized using ggplot2 (72) and sjPlot (73).

To determine how different fertilizer types affected the abundance of any specific OTU, we conducted an Analysis of Composition of Microbiomes [ANCOM, (74)], using the “Composition” plug-in for the QIIME2 platform. ANCOM is a robust analysis tool that uses compositional constraints to reduce false discoveries and can improve identification of differentially abundant microbes in complex datasets (75). Differences in relative abundances across treatments were compared at each crop cycle for both fungi and prokaryotes. To differentiate changes that could be confounded by time, relative abundances of OTUs were also compared across crop cycles. This occurred for crop cycles (1–6) for fungi, and because of computational power constraints even at the supercomputer level, we minimized analysis of prokaryotes at intervals of crop cycles (1–4 and 4–6). In addition, we performed a similar test, “log2FoldChange”, for changes in the microbiome using the package DESeq2 (76) to detect significant changes in OTU relative abundance between urea and compost treatments for cycle 1 and cycle 6.

All microbial association networks were created using the R package SPIECEASI (77). The SPIECEASI method was chosen because it was developed for compositional data analysis and assumes taxon-taxon associations scale linearly with the number of measured taxa. It also features greater precision when compared to other methods such as CoNet, gCoda, SparCC, and Spearman network inferences (34). Best-practices in microbial association network construction require large sample sizes to increase the precision of networks (34, 77). Although the sample size in our study was small (n = 3 for each treatment), we found that time across the 6 cropping cycles, although significant, had a much smaller effect in shaping microbial community structure when compared to the effect of the fertilization treatments. The effect of time across cropping cycles had an R² value ranging from 0.022 to 0.067 for fungi and 0.024 to 0.062 for bacteria & archaea, as compared to the effect of fertilizer type with R² ranging from 0.298 to 0.387 for fungi and 0.041 to 0.117 for bacteria & archaea (Table 2). Because time across cropping cycles explained a much smaller portion of the variation compared to treatment, we combined samples from all crop cycles (C1–C6) for each fertilizer treatment, which allowed for network inference with a pseudo n = 18 (3 pots x 6 cycles). Using these combined data, we created prokaryote networks for bacteria & archaea alone, fungi alone, and a combined network of bacteria, archaea, and fungi. We recognize that the prokaryotic network is technically a cross-domain network, but for ease of discussion, will refer to the combined network as the cross-domain network henceforth.

To reduce the likelihood of false positive associations within networks, OTU counts were first filtered for minimum occurrence across ~33% of samples using the R package Phyloseq (78). OTUs with 0 abundance were filtered out of any individual sample, but their sums across the whole dataset were kept so as not to change the overall sample counts. SPIECEASI networks were then assembled via the “spiec.easi” function using the Meinshausen and Bühlmann or “MB” method (79) with a pulsar parameter threshold of 0.05 and lambda minimum ratio of 1e⁻² following basic user guidelines (77, 80). The function “multi.spiec.easi” was used for creating cross-amplicon networks under the same parameters. Resulting networks were visualized in Cytoscape V-3.7.1 (81) and network topological properties were assessed using NetworkAnalyzer (82). Subnetworks were created by selecting individual hub OTUs from the CytoScape interface, selecting nearest-neighbors of those OTUs, and then creating a new network from the combined selections.

Network hub taxa were selected based on node degree (the number of edges linked to that node) and closeness centrality (how close a node is to all other nodes) as measurements of overall network connectivity using NetworkAnalyzer in Cytoscape. Currently, there is no consensus for what determines statistically significant network hub taxa, although previous studies determined good candidates by selecting those with higher node degree, betweenness centrality, or closeness centrality based on normal distribution fit with p < 0.1 and correlation cut-offs (83, 84). The SPIECEASI method of network inference is a conservative method that although results in fewer connections than other inference tools, can help to reduce false positives (77, 80). Therefore, we took a simpler approach for determination of hub taxa that did not dismiss non-outliers by choosing the top 5 OTUs that contained the highest degree and closeness centrality for each treatment. We limited our network hub taxa analysis to only the cross-domain networks as these were the most comprehensive network inferences.

At the end of the experiment, fertilization treatments did not significantly affect soil pH, which averaged 7.8 across all treatments (Table 1). There was a significant treatment effect on soil total N (p = 0.0002) and a marginally significant effect on soil carbon (p = 0.0696). Compost increased concentrations of soil nitrogen and carbon relative to both urea and control. There was no significant treatment effect on cation exchange capacity by the end of the sixth crop cycle. Calcium was the dominant cation across all treatments. Compost significantly increased concentrations of potassium (p = 0.0031), magnesium (p = 0.0069), and sodium (p = 0.0005).

Plant growth responded differently to treatments in the six consecutive cropping cycles following a single initial application of urea or compost (Figure 1). When compared to the control, compost fertilization significantly increased both aboveground and belowground plant biomass for all crop cycles (p < 0.01). Urea fertilization, in contrast, only significantly increased above and belowground plant biomass in the first crop cycle (p < 0.001). When the treatments were compared to each other, compost had greater biomass relative to the urea treatment for crop cycles 2–6 (p < 0.001).

Fertilizing soils with either compost or urea affected microbial richness and diversity differently (Figure 2). Fertilizer type did not affect prokaryote richness (p = 0.850), Shannon's Index (p = 0.984) nor Faith's Phylogenetic Diversity (p = 0.635). These metrics remained stable across the six cropping cycles (p > 0.4). These patterns were somewhat different for fungi where fertilization did not affect overall richness (p = 0.326) and Shannon's Index (p = 0.628) when compared to unfertilized soils, but over the six cropping cycles richness (p = 0.001) and diversity (p = 0.004) significantly decreased. Fungal richness and diversity were highest immediately following fertilization of field-collected soil prior to first sowing of seeds, then continued on a decreasing trend until the sixth cycle. These effects were observed for both urea and compost fertilizer sources.

We compared changes in microbial communities of different fertilizers across the six cropping cycles relative to the first cycle and found that fertilizer type significantly changed soil microbial community composition, but these effects were latent in prokaryotic communities (Figure 3, Tables 2, 3). Fungal communities in compost-fetilized soil differed significantly from unfertilized controls (p < 0.013) and urea (p < 0.013), but did not differ between urea-fertilized and unfertilized soils (p > 0.327) (Table 3). The variation in the data explained (R²) across the six cycles ranged from 0.297 to 0.386. Prokaryotic communities did not respond significantly to fertilization until the third crop cycle (p = 0.025) with increasing statistical significance until the sixth cycle. Similar to the fungal communities, prokaryotic communities from compost-fertilized soils differed from unfertilized control (p < 0.048) and urea (p < 0.036), while communities in control and urea-fertilized soils were not statistically different (p > 0.27). The variation in the data explained (R²) across the six cycles ranged from 0.249 at the first crop cycle, and consistently decreasing to 0.041 in cycle 6.

Fertilizer type significantly changed the relative abundance of certain individual fungal OTUs based on ANCOM analysis, but trends were more prominent for fungi than for prokaryotes across the 6 cycles (Supplementary Table 1). Two fungal OTUs, Thermomyces lanuginosus and a Myceliophthora sp. were more abundant in compost-fertilized soils. Thermomyces lanuginosus was more abundant for crop cycles 1 (w-score = 142), and crop cycle 6 (w = 146), while Myceliophthora sp. was more abundant in the second (w = 283), third (w = 245), fourth (w = 227), and fifth (w = 205) crop cycles. A Pezizaceae sp. was more abundant in the control and urea treated soils, but only at crop cycle 4. Only one bacterial OTU, an unidentified Acidomicrobiales, was more abundant in the control treatment at the first crop cycle (w = 3952). No differences in the relative abundance of prokaryotic OTUs were detected across all crop cycles.

Fold change analysis of OTU relative abundance in urea vs. compost treatments showed that prokaryotic taxa responded more strongly in crop cycle 1, whereas fungi responded more strongly in crop cycle 6. We detected 214 OTUs within 42 taxa (mostly at the genus level but some could be not assigned) that responded either strongly to urea or compost fertilization (Figure 4). Of these, OTUs of Actinomadura, Catellatospora, Iamia, Phytohabitans, Parasegitibacter, unidentified Chloroflexi, Nitrospira, Gemmata, Planctomyces, Arthrospira, Dechloromonas, Devosia, Luteibacter, Massilia, Tepidimonas, Opitutus, and Prosthecobacter significantly increased in abundance in urea-fertilized soils. In contrast, Kibdelosporangium, Mycobacterium, Rubrobacter, Streptomyces, Flavisolibacter, an unidentified Nitrospiraceae, unidentified Pirellulaceae, Geobacter, Hyphomicrobium, and Reyranella significantly increased in abundance in compost-fertilized soils. Other genera contain OTUs that show different preferences for the different types of fertilizers. No differences in the relative abundance of prokaryotic OTUs were detected at crop cycle 6. We found only one OTU of fungi, Thermomyces lanuginosus, that significantly favored compost at crop cycle 1, but six OTUs at crop cycle 6. These were Acremonium dichromosporum, T. lanuginosus, a Myceliophthora sp., and three unidentified Ascomycota in the orders Sordariales and Eurotiales. General visual comparisons of all OTUs aggregated at the phylum level across all 6 cycles can be found in Supplementary Figures 3, 4.

Soils treated with fertilizers altered network topography for the individual prokaryotic network, the fungal network, and the cross-domain network (Figure 5). Network complexity, as measured by the number of nodes and edges, increased when compost was added to the soil, and decreased when urea was added, relative to the control. For the prokaryotic network, the control contained 134 nodes and 150 edges, the compost-fertilized network was higher with 141 nodes and 217 edges, and the urea-fertilized network was fewer with 117 nodes and 154 edges. We observed a similar trend for the fungal networks. The control network had 14 nodes and seven edges, the compost-fertilized network was higher with 44 nodes and 31 edges. The urea-fertilized network, however, had 18 edges and 11 nodes, both of which were higher than the control. The cross-domain networks followed this trend as well where the control network had 229 nodes and 372 edges, the compost-fertilized network had the highest number of 248 nodes and 408 edges, and the urea network had the lowest number of 207 nodes and 307 edges.

Within the networks, fertilizer type had a strong effect on network hub taxa. Not only did fertilizer type affect which OTUs became hub taxa, but it affected how they were connected in the overall networks (Figure 6). The compost-fertilized network contained three bacterial (a Nostocaceae sp., Steroidobacter sp., Bradyrhizobiaceae sp.) and two fungal (Lasidioploidia lignicola, and Tetracladium furcatum) hub taxa. The unfertilized control network contained three bacterial (a Nostocaceae sp., Rhizobiales sp., Micrococales sp.) and two fungal (Hypocreales sp., and Aspergillus purpureus) hub taxa. The urea-fertilized network contained one bacterial taxon (a Solirubrobacteriales sp.) and four fungal (Ascomycota sp., Aspergillus sp., Trechispora sp., and an unidentified fungus) hub taxa. Only one OTU, a member of the Nostocaceae, was present as a hub taxon in both the unfertilized soil and compost-fertilized network, but not the urea-fertilized network. Each hub taxon in the compost-fertilized network was directly connected to other hub taxa, and indirectly connected to 2 or more other hub taxa. In stark contrast, the hub taxa in the unfertilized network shared only indirect connections to other hub taxa, and averaged less than two connections to other hub taxa. Hub taxa in the urea-fertilized network similarly had fewer direct and indirect connections, including one hub taxa that was not connected to any other hub taxa by nearest neighbors.

A comparison of OTUs found in the ANCOM analysis with those from the hub taxa analysis showed that the taxon Thermomyces lanuginosus (significantly more abundant in compost-fertilized soil) was directly connected to a hub taxon, a Steroidobacter sp. (Supplementary Figure 5). The Myceliophthora sp. was found in a separate node cluster and was neither a hub taxon nor connected to one (data not shown). No other taxa from the ANCOM analysis were found to interact closely with network hub taxa of any treatment group. The same T. lanuginosus and Myceliophthora sp. were also detected in the fold change analysis, favoring compost-fertilized soil.