High‐Resolution Community Profiling of Active Bacteria and Eukaryotes in Replant‐Diseased Blueberry Farm Soils From New Jersey, USA

TL;DR

This study identifies differences in soil microbes between low- and high-yield blueberry farms, which could help diagnose and improve replant disease.

Contribution

The study introduces a method combining long-read sequencing and SIP to profile active microbes in blueberry farm soils.

Findings

Low-productivity soils showed active Bacillus species, while high-productivity soils had active Burkholderia and Paraburkholderia.

Eukaryotic species like Candida blankii were enriched in low-yield soils.

The method can differentiate soil microbiomes and may help detect harmful microbes.

Abstract

Highbush blueberry (
Vaccinium corymbosum
L.) fields can remain productive for decades. However, some older fields decline in plant health and exhibit lower yields. After re‐planting with new stock, the yields continue to suffer. This condition is termed ‘Replant Disease’. The causative agent(s) in replant disease in New Jersey blueberry fields are unknown. To assess if low‐ and high‐yield blueberry farm soils from two separate farms contained different microbiomes, we coupled long‐read bacterial and eukaryotic ribosomal rRNA operon sequencing using the Oxford Nanopore MinION with stable isotope probing (SIP) to detect 13C/15N‐utilising soil microbial communities. The results indicate multiple Bacillus species were active on 13C/15N‐growth media (predominantly amino acids and glucose) in low‐productivity soils from both farms, while high‐productivity and adjacent forest soils contained active Burkholderia and Paraburkholderia species. Eukaryotic community profiling indicated Candida blankii, Nadsonia starkeyi and Sugiyamaella chiloensis were slightly enriched and active in low‐productivity soils compared with high‐productivity soils. This approach differentiates low‐ and high‐productivity blueberry farm soils by ribosomal RNA operon profiling and SIP. The findings also suggest a diagnostic test of blueberry replant affected soils is feasible and may ultimately be used to improve productivity and potentially detect the responsible pathogenic agent(s) or other deleterious microbes.

Long read bacterial and eukaryotic ribosomal operon sequencing combined with stable isotope probing (SIP) allowed the differentiation of microbiomes in low‐ and high‐yield blueberry farm soils. This report may ultimately be used to improve productivity and potentially detect the responsible pathogenic agent(s) or other deleterious microbes.

Figures

Bubble plot comparing differences in relative abundance for the top 10 13C–15N‐active bacterial taxa identified by DESeq2 analysis for low‐ and high‐productivity soils. Classification was done using the bacterial ribosomal operon database (rOPDB). Red boxes indicate active taxa from Farm 1 and Farm 2.

Bubble plot comparing differences in relative abundance for the top 10 13C–15N‐active eukaryotic taxa identified by DESeq2 analysis for low‐ and high‐productivity soils. Classification was done using the eukaryotic ribosomal operon database (GoldenRod). Red boxes indicate active taxa from Farm 1 and Farm 2.

Bubble plot comparing differences in relative abundance for the top 10 12C–14N‐resident bacterial taxa identified by DESeq2 analysis for low‐ and high‐productivity soils. Classification was done using the bacterial ribosomal operon database (rOPDB). Red boxes indicate resident taxa from Farm 1 and Farm 2.

Bubble plot comparing differences in relative abundance for the top 10 12C–14N‐resident eukaryotic taxa identified by DESeq2 analysis for low‐ and high‐productivity soils. Classification was done using the eukaryotic ribosomal operon database (GoldenRod). Red boxes indicate resident taxa from Farm 1 and Farm 2.

Bubble plot comparing differences in relative abundance for the top 10 bacterial taxa identified by DESeq2 analysis. Both 12C–14N‐resident and 13C–15N‐active fractions are indicated in soils from Forest 1 and Forest 2. Classification was done using the bacterial ribosomal operon database (rOPDB).

Bubble plot comparing differences in relative abundance for the top 10 eukaryotic taxa identified by DESeq2 analysis. Both 12C–14N‐resident and 13C–15N‐active fractions are indicated in soils from Forest 1 and Forest 2. Classification was done using the eukaryotic ribosomal operon database (GoldenRod).

Full text

Introduction

1

Highbush blueberry ( Vaccinium corymbosum L.) is a long‐lived woody perennial native to North America which grows in organic, well‐aerated acidic sandy soils (pH of 4–5) (Hayes 1988). The United States produced over 725 million pounds of blueberries in 2024 (National Association Blueberry Council report 2024) with a production value of > $800 million, due to their richness in natural antioxidants and other bioactive compounds (Prior et al. 1998). Highbush blueberry fields can remain productive for decades, but older fields tend to decline in plant health and blueberry yield (Figure S1). Replanting often does not solve this problem, suggesting the issue is soil‐associated (Jagdale et al. 2013; Noe et al. 2014). This condition, termed replant disease, is common in woody perennial crops and is well studied in fruit trees (Mazzola and Manici 2012; Waschkies et al. 1994). However, replant disease in blueberry and its possible association with the soil microbiome is not well studied. Plant‐parasitic nematodes have been implicated in replant disease in Georgia and North Carolina (Jagdale et al. 2013; Noe et al. 2014). Phytophthora cinnamomi is known to cause root rot in blueberry and cranberry ( Vaccinium

macrocarpon); however, Phytophthora spp. are not associated with replant disease in New Jersey blueberries (Kawash et al. 2023).

One of the first reports of healthy microbial communities in the wild lowbush blueberry ( Vaccinium angustifolium ) fields was from Nova Scotia, Canada, which utilised short‐read sequencing (200–400 bp) of the V6–V8 region of the 16S rRNA gene for bacteria identification plus the V4 region of the 18S rRNA for fungal profiling (Yurgel et al. 2017). Subsequently, a comparison of bacterial and eukaryotes from managed and forested wild lowbush blueberry patches using the same target genes, found productivity to be negatively correlated with several fungal and protist taxa (Yurgel et al. 2018). Other studies utilised the same target SSU gene sequences for bacterial community characterisation and added the intergenic transcribed spacer (ITS) region of the eukaryotic ribosomal RNA operon for fungal identification (Li et al. 2020). Additional efforts employed a different region of the 16S rRNA gene (V3–V4) for bacterial identification and the ITS region for fungal identification (Li et al. 2020; Zhou et al. 2022). Nearly all of this prior research found Acidobacteriota, Actinomycota, Bacillota, Pseudomonadota and Verrumicrobiota as predominant bacterial phyla, while Ascomycota, Basidiomycota, Cryptomycota and Mucoromycota as predominant fungal phyla. However, the study by Yurgel et al. (2018) did not find any correlation with bacterial taxa and low fruit production, and only a small percentage of the fungal diversity could be associated with lower fruit productivity (Adonis test: R ^2^ = 0.04; p < 0.05). This inability to discern marked differences between low‐ and high‐ productivity soils could be a result of microbiome classification at taxonomic levels which are too high (i.e., grouping of SSU data by phylum/class/family level) or by not considering the active fraction of the microbiome (i.e., those microbes responding to various growth substrates vs. dormant cells, etc.).

As an alternative approach to assess microbiome differences between replant affected and healthy soils from blueberry farms, we coupled ribosomal RNA operon sequencing with stable isotope probing (SIP) for high‐resolution community profiling of bacteria and eukaryotes. The reasoning was that prior difficulties in observing microbial community differences associated with replant disease were caused by a lack of taxonomic resolution at a species‐ or strain‐level using traditional short‐read approaches. Alternatively, discerning the metabolically active fraction of microbial communities could help differentiate low‐ and high‐productivity soils when differences in microbial relative abundance were not as readily apparent. For this study, we chose to utilise the long‐read capability of the Oxford Nanopore Technologies Ltd. (ONT) MinION sequencer to sequence prokaryotic and eukaryotic rRNA operons (Benítez‐Páez and Sanz 2017; Kerkhof et al. 2017; Cusco et al. 2018), since it can provide species‐ and strain‐level identification of bacteria (Jain et al. 2015, 2018; Dowden et al. 2020; Ibironke et al. 2020; Kerkhof et al. 2022), while profiling eukaryotic rRNA operons have been used for fungal characterisation in mock communities, canine outer‐ear samples, or turf grasses (Mafune et al. 2020; D'Andreano et al. 2020; Groben et al. 2023).

As proof‐of‐concept, we also employed SIP to delineate the active members of the blueberry rhizosphere soil microbial communities using ^13^C–^15^N‐BioExpress growth media (Cambridge Isotope Laboratories Inc., Tewksbury). SIP is based on the uptake and metabolism of ^13^C and/or ^15^N substrates by microbes which then synthesise ^13^C‐ and/or ^15^N‐DNA during growth (Radajewski et al. 2000). BioExpress media is nearly the same as Luria–Bertani (LB) medium with 0.1–0.5 g glucose per litre. LB medium is reported to be rich in oligopeptides and contain all proteinogenic amino acids (predominantly Asp + Asn, Leu, Glu + Gln and Pro) (Sezonov et al. 2007). Interestingly, there is evidence that amino acids/peptides and carboxylic acids with associated aromatic rings are the dominant root exudates for Avena barbarta late in plant development before senescence (Zhalnina et al. 2018). Additionally, there is a report of l‐glutamic acid as a major exudate in blueberries during growth/development and is thought to be a major driver of the microbial community (Che et al. 2024). As such, our ^13^C–^15^N‐BioExpress‐based SIP results may reflect both utilisation of a ‘typical laboratory media’ and actual rhizosphere activity in the NJ soils. However, it should also be noted that BioExpress is a copiotrophic media, which will select for heterotrophic bacteria/eukaryotes capable of consuming organic substrates in high concentrations. The soil microbes in the studied fields may be experiencing different concentrations of these amino acid/glucose substrates or different exudates entirely. Regardless, the initial aim of our SIP approach is to first differentiate low‐ and high‐productivity soils by the analysis of microbial communities, which can ultimately be used as biomarkers for further investigation of replant diseased soils.

In this report, we reveal both the resident and the active bacterial and eukaryotic (predominantly fungal) community members which assimilate ^13^C‐^15^N‐amino acids or ^13^C‐glucose in low‐ and high‐productivity blueberry farm soils as well as forest soils harbouring wild blueberries. We hypothesised that distinct active communities are linked to different productivity levels. Our results indicate that multiple Bacillus spp. and three fungal taxa (Candida bankii, Nadsonia starkeyi‐henricii and Sugiyamaella chiloensis) actively metabolise the ^13^C–^15^N‐amino acids or ^13^C‐glucose in low‐productivity soils. In contrast, active Burkholderia and Paraburkholderia spp., as well as Mucor lusitanicus, labelled with ^13^C‐^15^N are associated with high‐productivity soils. Additionally, mycorrhiza‐forming Glomeromycota are enriched in high‐productivity soils, while Cryptomycota (also known as Rozellomycota), which are parasites of other fungi, are more prevalent in low‐productivity soils. This study demonstrates the feasibility of identifying actively growing microorganisms by sequencing near full‐length rRNA operons labelled with ^13^C/^15^Nisotopes. Our results also establish a foundation for developing molecular biomarkers to differentiate low‐ and high‐productivity soils in blueberry fields and may ultimately lead to a mechanistic understanding of replant disease in these systems.

Methods and Materials

2

Sample Collection and Stable Isotope Probing

2.1

Farm field soils (visually exhibiting low‐ and high‐productivity as in Figure S1) were hand sampled from two commercial blueberry farms, along with soils from two adjacent forests containing wild blueberries in 2019. The six sites were sampled in triplicate. Rhizosphere soil was collected, with a small shovel, from the root zone (between 30 and 40 cm from the crown) of 3–4 mature ‘Bluecrop’ blueberry plants per sample location to a depth of 16 cm and homogenised. Five hundred milligrams of homogenised soil was placed in tubes and amended with either light (^12^C/^14^N) or heavy (^13^C/^15^N) 1× (BioExpress 1000) Growth Media to fully wet the soil (e.g., 100 μL of 10× media plus 900 μL of sterile water). The amended soil samples were incubated for 2 days in a lighted growth chamber (Percival Scientific, Perry, IA) with a 12‐h day/night cycle and fluorescent illumination (80 μmol m^−2^ s^−1^) at 25°C, after which the soil was collected and processed as detailed below.

DNA Extraction

2.2

Total genomic DNA was extracted from the soil using a modified CTAB/phenol‐chloroform procedure (Männistö et al. 2013). Specifically, 50 mg of soil was wetted with 50 μL of Solution 1 (50 mM glucose, 10 mM EDTA, 25 mM Tris‐Cl (pH 8.0)) and subjected to five freeze–thaw cycles. Each sample was then incubated with 50 μL of lysozyme solution (4.0 mg/mL in Solution 1) for 10 min at room temperature, and 450 μL CTAB solution (100 mM Tris‐Cl, 20 mM EDTA, 1.4 M NaCl, 4% (w/v) CTAB and 2% β‐mercaptoethanol) was added. This aqueous mixture was extracted twice with 800 μL phenol: chloroform: isoamyl alcohol (24:25:1) to remove proteins/lipids. The DNA was ethanol precipitated with 20 μg of glycogen, resuspended in 50 μL of 10 mM Tris‐HCl and stored at −80°C.

Separation of Light and Heavy DNA

2.3

The ^12^C/^14^N‐ and ^13^C/^15^N‐labelled DNA was physically separated via caesium chloride (CsCl) gradient ultracentrifugation (225,000g for 36 h) using an Optima MAX‐TL ultracentrifuge with a TLA‐110 rotor (Beckman Coulter Inc., Indianapolis, IN, USA). To enhance the visualisation of the ^12^C/^14^N‐and ^13^C/^15^N‐labelled genomic DNA from the soils, each gradient was amended with ^12^C/^14^N‐ and ^13^C/^15^N‐labelled archaeal carrier DNA from Halobacterium salinarum (Gallagher et al. 2005) and ethidium bromide (Figure S2). This ethidium bromide addition alters the buoyant density of DNA (Gallagher et al. 2010) and nullifies any GC effects during DNA separation, allowing for easy visualisation of the ^13^C/^15^N fraction for amplification. The ‘heavy’ carrier method also significantly shortens incubation times and improves sensitivity when detecting ^13^C–^15^N bacterial or ^13^C–^15^N eukaryotic DNA (Gallagher et al. 2005).

The ^13^C–^15^N‐labelled archaeal DNA clearly defines the ‘heavy’ region within the cesium gradient and the archaeal DNA is not amplified with bacterial or eukaryotic universal PCR primers. This ^13^C–^15^N carrier approach has been used in a variety of complex systems to discern active bacteria, archaea and eukaryotes in soils/sediments, deep sea samples or at sub‐zero temperatures (Gallagher et al. 2005; Gallagher et al. 2010; Kerkhof et al. 2011; Seyler et al. 2014; Seyler et al. 2018; Seyler et al. 2019; Tuorto et al. 2014; Liu et al. 2016; Gadkari et al. 2020). The light (^12^C/^14^N‐labelled) DNA contained in the top band represents the ‘resident’ community. The heavy (^13^C/^15^N‐labelled) DNA contained in the bottom band of the CsCl gradient represents the ‘active’ community. Both top and bottom bands were collected by pipette and CsCl was removed by micro‐dialysis against 10 mM Tris‐HCl. Control SIP incubations were performed with ^12^C/^14^N‐amendments of Bioexpress and did not amplify the ^13^C/^15^N‐bottom bands using bacterial or eukaryotic rRNA operon primers (Figure S3).

PCR Amplification of Bacterial and Eukaryotic Ribosomal RNA Operons

2.4

A two‐step amplification/barcoding procedure using Hi‐Fidelity Taq polymerase (Bimake LLC, Houston, TX, USA) described in Kerkhof et al. (2017) was used to amplify target rRNA operons for sequencing rather than employing a barcoding kit provided by ONT. Specifically, universal prokaryotic 16S forward 27F (5′AGAGTTTGATCCTGGCTCAG 3′) and 23S reverse 2241R (5′ ACCGCCCCAGTHAAACT 3′) primers were first used to amplify bacterial ribosomal RNA operon fragments from each CsCl band using a 30‐cycle touchdown amplification described in further detail in Data S2. At the end of the 16th cycle, 10 μL of the PCR amplification mixture was removed and stored frozen. The PCR reaction then was allowed to proceed to 30 cycles to ensure that the top ^12^C/^14^N bands would amplify and the ^13^C/^15^N heavy DNA bands would not show any PCR amplification in the ^12^C/^14^N‐treated control samples (Figure S3). Once these control conditions had been met, the low cycle number amplicons (16 cycle) were purified using Ampure beads (Beckman Coulter Inc.) and barcoded via a second round of PCR with tailed PCR primers containing a specific barcode sequence for multiplex sequencing. A similar approach was used to generate barcoded eukaryotic rRNA operons with the primers Euk 18S F 5′ TACCTGGTTGATYCT GCCAGT 3′ and Euk 28S 3050 R 5′ TCTGACTTAGAGGCGTTCAGTCATAAT 3′ as described more fully in Data S2.

Library Preparation and MinION Sequencing

2.5

Barcoded rRNA operon amplicons were quantified by agarose gel electrophoresis and image analysis. Twelve amplicon samples (each sample having 150 ng of DNA (1800 ng in total)) were combined per multiplex DNA sequencing library, end‐repaired and dA‐tailed using NEB kits (New England Biolabs, Ipswich, MA, USA) as per manufacturer's instructions. The reactions were purified using Ampure beads, supplemented with 100 nM dNTP's and ligated to the ONT sequencing adaptor using Blunt/TA ligase (NEB) amended with 1 μL of freshly made ATP solution (~4 mg/mL). All sequencing was done using the ONT SQK‐LSK109 kit and analysed on FLO‐MIN106 R9.4.1 flow cells. The sequencing of PCR negative controls was not performed since prior research provided an estimate of 0.03% contamination from the ‘kitome’ (Ibironke et al. 2020).

Demultiplexing and BLAST Analysis

2.6

All raw reads were base‐called with the Guppy 3.2.2 software (ONT) via the high‐accuracy algorithm. Each library resulted in 2–3.5 million raw reads. Demultiplexing was performed via a custom workflow using the Geneious R11 software (Biomatters Ltd.), which included: filtering the DNA sequence reads by size (3700–5700 bp for bacterial rRNA operons; 4000–6300 bp for eukaryotic rRNA operons), truncating the assigned ONT name, annotating the reads with barcode primers (90% pairwise identity) and extracting copies of the reads containing specific barcodes into separate folders. An assessment of whether annotated sequences could be assigned to multiple barcode folders indicated a 0.375% ‘contamination’ of mis‐assigned reads per sample. This low percentage of mis‐assigned reads were deemed too low to warrant extensive detection and clean‐up.

All barcoded raw reads following QA/QC were then exported as fasta files and initially screened against the EZBioCloud 16S rRNA gene database (Yoon et al. 2017) for bacterial taxonomic assignments and the UNITE/all eukaryotes ITS database (UNITE Community 2019) for eukaryotic assignment. These databases were available in 2019. MegaBLAST parameters included: word size of 28 bp, match/mismatch of 1/−2, a linear gap penalty and maximum target sequences of 1. Subsequently, the raw rRNA operon reads were re‐classified using more recently available long‐read databases available in 2022 and 2024, specifically the bacterial (Kerkhof et al. 2022) and the eukaryotic rRNA operon databases (Krabberød et al. 2024). MegaBLAST parameters for screening against the longer databases were: word size of 65 bp, match/mismatch of 2/−3, gap open/extend of 0/4 and maximum target sequences of 1.

Statistical Analysis and Data Visualisation

2.7

The BLAST results were imported as hit tables in CSV format to Microsoft Excel. Pivot tables were generated based on the counts of species/strains for all soil samples originating from both ^12^C/^14^N (light) and ^13^C/^15^N (heavy) DNA bands and imported to R studio (2023.06.0 + 421), using R version 3.6.2. Non‐metric multidimensional scaling (NMDS) was performed on the species abundance data using the ‘vegan’ package version 2.5–6 (Oksanen et al. 2018). The number of dimensions used for the analysis (k = 3) was based on the NMDS stress plots and R ^2^ values for the model fit. The R ^2^ value for non‐metric fit of observed dissimilarities was 0.99 for the eukaryotic species‐abundance data and 0.997 for bacterial data with k = 3.

Differential Abundance Analysis

2.8

The Excel pivot tables were imported into the DESeq2 package (Love et al. 2014) to perform differential abundance analysis between resident (^12^C^14^N) and active (^13^C^15^N) taxa from low‐ and high‐productivity farm soils and adjacent forest soils in R studio. The ggplot2 package was used to generate volcano and bubble plots (Wickham 2016). Adjusted p < 0.05 and log2‐fold‐change > 2 were considered threshold values in DESeq2 for differentially abundant taxa.

Results

3

Overall, 1,794,249 bacterial rRNA operon reads were generated, while the eukaryotic rRNA profiling produced 3,202,957 reads. The average read numbers were ~50,000 per sample for bacteria and ~89,000 per sample for eukaryotes. These rRNA operon sequences were classified against both short‐read databases (EZBioCloud [bacterial] and UNITE [eukaryotic]) and long‐read databases (rOPDB [bacterial] and GoldenRod version 1.2 [eukaryotic]) by MegaBLAST as described above.

Taxonomic Classification of Bacterial rRNA Operon Reads

3.1

Screening by EZ BioCloud detected 10,897 bacterial best BLAST hits (BBH or operational taxonomic units) at the species level, belonging to 44 different bacterial phyla. The bacterial phyla most abundant in the samples included the Pseudomonadota, the Actinomycota, the Acidobacteriota, the Dormibacteraeota and the Bacillota, accounting for > 90% of the taxa in the ^12^C/^14^N labelled (resident) community (Figure S4B). The reproducibility of rRNA operon profiles and the challenge of identifying differing bacterial taxa in replant disease‐affected soils are evident in relative abundance bubble plots of dominant bacterial BBHs (Figure S5). Most samples displayed identical taxa with similar relative abundances. Only two ^12^C/^14^N‐labelled taxa in low‐productivity soil had higher abundances: Steroidobacter‐related (DQ480487_s) and Rhodospirilliaceae‐related (AY913269_s; Figure S5 [resident]). These database entries were not well characterised beyond the Order/Family level in EZBioCloud. In contrast, the ^13^C–^15^N‐labelled taxa in low‐productivity soils were enriched with Bacilli (Bacillus spp.): * B. bataviensis, B. cucumis, B. drentensis, B. fumaroli, B. niacini, B. soli, B. novalis *, B. pocheonensis and B. spp. (DQ129292_s; Figures S4B and S5).

The re‐classification of bacterial rRNA operon reads against the rOPDB near‐full length ribosomal operon database provided in a 2.6‐fold increase in query sequence/reference alignment length thus improving BLAST classification results as seen in the scatter plot of BBHs. There is a shift away from alignments to the 16S rRNA gene (1500 bp alignment) toward the full operon (4200–5000 bp) with roughly the same percent identity (70%–95%) for the SSU matches (Figure S4A) and the near full operon (Figure S8). More importantly, the number of taxa detected using the rOPDB increased from 10,897 to 12,418 along with the taxonomic resolution of the BBHs to species/strain level.

Taxonomic Classification of Eukaryotic rRNA Operon Reads

3.2

Nineteen major clades of eukaryotic organisms were identified by classification with the UNITE database, constituting 51 eukaryotic phyla and 2501 distinct BBHs (or operational taxonomic units) in the soil sample set. Soil communities were predominantly fungal, but the soils also contained numerous protist and metazoan reads (data not shown). The top resident eukaryotic phyla were Ascomycota, Basidiomycota, Glomeromycota, Mortierellomycota, Mucoromycota and Rozellomycota, while the active ^13^C–^15^N‐labelled soil communities were more enriched with Mortierellomycota species in all soils (low‐, high‐productivity and forest) (Figure S6B). In low‐productivity soils, there was an enrichment of an unclassified Rozellomycota taxa in the ^12^C–^14^N‐labelled fraction. In the ^13^C–^15^N‐labelled fraction, the differences between the low‐ and high‐productivity soils were very minor. Additionally, low‐productivity soils had the slightest enrichment in a member of the Pythiaceae family in the ^13^C–^15^N‐labelled fraction.

Re‐classifying the eukaryotic rRNA operon reads with the GoldenRod database also provided 2.5‐fold increase in alignment length, and an increase in the number of assigned taxa was observed as well (from 2504 to 3144). Likewise, the taxonomic resolution was at the species level. However, some eukaryotic taxa classified using the UNITE database (e.g., members of the Rozellomycota) were not well represented in the GoldenRod database and would likely have been misclassified.

Differential Abundance Testing of Bacteria and Eukaryotes

3.3

To identify the taxa exhibiting differential abundance, DESeq2 analysis was performed on classifications using the long‐read databases. The top 10 bacterial and eukaryotic taxa demonstrating both high community relative abundance and differential abundance (log2‐fold‐change > 2 and adjusted values p < 0.05) from the DESeq2 analysis are presented in Figures 1, 2, 3, 4, 5, 6. The long‐read analysis indicated the active bacterial community from both Farm1 and Farm2 were similar in low‐ and high‐ productivity soils. Specifically, ^13^C–^15^N‐rRNA operon reads in low‐productivity soils were related to Bacillus megaterium S2, B. mesonae H20.5, B. soli DSM15604, four different strains of Bacillus sp. (BRMEA1, OV166, S3, X1) and Streptomyces sp. SAJ15 (Figure 1). While the high‐productivity soils had the highest ^13^C–^15^N‐active taxa associated with Paraburkholderia or Streptomyces species. Likewise, there were shared reads in the active fraction of eukaryotes incorporating ^13^C‐^15^N‐BioExpress (Figure 2). Specifically, reads associated with Candida blankii, Nadsonia starkeyi‐henricii and Sugiyamaella chiloensis were detected in the active fraction from low‐productivity soils in both Farms. The relative abundance of these active eukaryotes in low‐productivity soils was small, reflecting that these taxa are probably not consuming large amounts of amino acids or glucose‐rich substrates in situ.

Bubble plot comparing differences in relative abundance for the top 10 13C–15N‐active bacterial taxa identified by DESeq2 analysis for low‐ and high‐productivity soils. Classification was done using the bacterial ribosomal operon database (rOPDB). Red boxes indicate active taxa from Farm 1 and Farm 2.

Bubble plot comparing differences in relative abundance for the top 10 13C–15N‐active eukaryotic taxa identified by DESeq2 analysis for low‐ and high‐productivity soils. Classification was done using the eukaryotic ribosomal operon database (GoldenRod). Red boxes indicate active taxa from Farm 1 and Farm 2.

Bubble plot comparing differences in relative abundance for the top 10 12C–14N‐resident bacterial taxa identified by DESeq2 analysis for low‐ and high‐productivity soils. Classification was done using the bacterial ribosomal operon database (rOPDB). Red boxes indicate resident taxa from Farm 1 and Farm 2.

Bubble plot comparing differences in relative abundance for the top 10 12C–14N‐resident eukaryotic taxa identified by DESeq2 analysis for low‐ and high‐productivity soils. Classification was done using the eukaryotic ribosomal operon database (GoldenRod). Red boxes indicate resident taxa from Farm 1 and Farm 2.

Bubble plot comparing differences in relative abundance for the top 10 bacterial taxa identified by DESeq2 analysis. Both 12C–14N‐resident and 13C–15N‐active fractions are indicated in soils from Forest 1 and Forest 2. Classification was done using the bacterial ribosomal operon database (rOPDB).

Bubble plot comparing differences in relative abundance for the top 10 eukaryotic taxa identified by DESeq2 analysis. Both 12C–14N‐resident and 13C–15N‐active fractions are indicated in soils from Forest 1 and Forest 2. Classification was done using the eukaryotic ribosomal operon database (GoldenRod).

Interestingly, the resident bacterial community from both Farm1 and Farm2 were also similar with respect to differential abundance in the low productivity soils (Figure 3). These low‐productivity soils were enriched in rRNA operon reads associated with 3 taxa associated with the Rhodospirillales (strain R5959, SYSU, URHD0088) and 2 taxa associated with Azospirillum sp. (MTCC4039 and TSO35.2). While the resident eukaryotic community displaying differential abundance only detected Stylonichia lemnae in the low‐productivity soils (Figure 4). In the high productivity soils, for Farm1 both the resident and the active fraction were enriched with Mucor lusitanicus (Figures 2 and 4). M. lusitanicus was not seen as differentially abundant in Farm2 or in the control site Forest1 (Figure 6). This taxon was a minor component of the resident and active fraction in Forest2. Another eukaryotic taxon shared between Farm1 and Forest1 was Podila verticulata, which was enriched in the resident fraction in Farm1 (Figure 2) and the active fraction of Forest1 (Figure 6).

None of the above active bacterial taxa associated with low‐productivity soils from farms were observed to be enriched or active in the control Forest soils (Figure 5). However, there were different Paraburkholderia and Burkholderia species than the ones observed to be active in high‐productivity soils, which were also active in the control Forest soils (Figure 3).

Discussion

4

This study was initiated to ascertain if high‐resolution profiling or detection of active microbiota could elucidate differences between low‐productivity versus high‐productivity soils on commercial blueberry farms. Our results demonstrate that rRNA operon profiling coupled with SIP can detect differences in soil microbiomes from replant disease affected areas. Specifically, our findings demonstrate that various Bacillus species were active in low‐productivity soils from two separate farms. It is known that many bacilli can synthesise antibacterial, antifungal and nematocidal substances (Pérez‐García et al. 2011). For instance, there are reports of B. soli DSM15604 possessing nematocidal activity (Englebrecht et al. 2020; Horak et al. 2022) B. mesonae H2O‐5 suppressing tomato wilt from Ralstonia solanacearum (Yoo et al. 2022), and B. megaterium strains promoting plant growth and possessing insecticidal properties (Migunova and Sasanelli 2021; Tarigan et al. 2022). Additionally, many Bacillus species are also reported to evoke induced systemic resistance (ISR) during various disease conditions in multiple plant species, resulting in reduced incidence or severity of disease (Choudhary and Johri 2009; Yang et al. 2023). These prior findings suggest the activity of Bacilli on ^13^C–^15^N‐ Bioexpress media might not be directly related to low productivity. Rather, this increase in Bacilli activity may be the response to plant disease and stress.

At this point we are not aware of any potential direct or indirect interactions of the Bacillus spp. identified by our rRNA operon reads with blueberries. It is important to recognise that finding a BBH in the rOPDB or GoldenRod database can often mean that novel species and strains are present in the sample rather than evidence that any specific strain from the database has been detected (Dowden et al. 2024). Indeed, the overall matches of our rRNA operon reads to B. soli (the most relatively abundant active taxa in replant disease‐affected soils), B. megaterium and B. mesonae from the rOPDB are 87%–88% over 4000–4085 bp. These lower percent identities suggest our rRNA operon reads represent novel genera or species within the Bacillota. For simplicity, we are reporting the BBH here as a means of conveying the overall character of the microbial community. However, understanding the species‐strain level character of replant disease‐affected soils is important since different bacterial strains can have very different functional properties. Furthermore, it is likely these ‘berryland soils’ are unique with respect to other major crops fields (e.g., soybean, wheat, corn, apples, etc.). Additional high‐resolution profiling coupled to SIP analyses is warranted for other diseased soils. Ultimately, it is hoped the rRNA operon profiling of active microbes (e.g., Bacillus spp. in this case) can provide critical biomarkers to aid in the isolation and characterisation of novel strains active in diseased soils and to test if the active microorganisms are involved in stress‐response to plant disease or may be responsible for the pathogenesis.

For the high‐productivity farm soils, various Paraburkholderia taxa appeared to be the most active with our SIP substrates. For example, Paraburkholderia species were found in all types of soil. However, these microorganisms were more active on ^13^C/^15^N‐amino acids or ^13^C‐glucose in forest and high‐productivity farm soils than in low‐productivity farm soils. Some Paraburkholderia representatives are thought to be beneficial plant symbionts able to fix nitrogen and promote plant health (Kaur et al. 2017). In our study, Paraburkholderia oxyphyla was found to be the most actively growing bacterial species in high‐productivity soils by EZBioCloud classification. This bacterium, first isolated from acidic forest soils, is not capable of N_2_‐fixation, unlike most of its phylogenetically closely related Paraburkholderia species. However, a novel species from the same genus, a symbiont capable of nitrogen fixation, was recently isolated from the nodules of Mimosa gymnas and was shown to share almost 99% similarity with Paraburkholderia oxyphyla at the 16S rRNA gene level (Otsuka et al. 2011; Paulitsch et al. 2019).

With respect to the short‐read classification of the resident eukaryotic community, our findings indicated high‐productivity blueberry farm soils were enriched with fungi from the Glomeromycota phylum. These fungi are an ancient group of eukaryotes and are thought to be obligatory plant symbionts, forming arbuscular mycorrhizas with the roots of many vascular plants (Oehl et al. 2011; Redecker et al. 2013; Schüßler et al. 2001). They are considered beneficial symbionts for the plants and are thought to contribute by fighting fungal pathogens/harmful nematodes, aid in the uptake of phosphate/other inorganic ions and increase the fitness of plants in polluted environments (Oehl et al. 2011; Redecker et al. 2013). It has previously been shown that multiple Glomeromycota species can enhance the height and the dry weight of blueberry seedlings and increase the levels of N and P in the plant roots (Farias et al. 2014). Glomeromycota species were also found in the rhizosphere of the wild blueberry plant and were relatively more abundant in the rhizosphere itself than in the bulk soil (Yurgel et al. 2018). Another study investigating mycorrhizal fungi in the roots of Vaccinium oldhamii Miq. found a unique OTU very similar to Rhizophagus diaphanous (Baba et al. 2016), which suggests that Glomeromycota species may, in fact, colonise the roots of ericaceous plants though previously reported to be symbionts of other vascular plants. Ericaceous plants are thought to be colonised by only ericoid mycorrhizae (Straker 1996), and formation of arbuscules has not been reported in Vaccinium spp. Thus, more research is required to understand whether these beneficial fungi are critical in the health of the highbush blueberry plant.

In contrast, short‐read classification by the UNITE database of low‐productivity soils from the two investigated farms indicated Rozellomycota species were enriched in the eukaryotic resident communities. Rozellomycota, which are mainly unicellular endoparasites of algae, amoebae, oomycetes and different fungal species. Most of the members of this phylum are uncultured environmental clones and not much information is known about their pathogenic lifestyle and potential host range (Corsaro, et al. 2014; Corsaro, Walochnik, Venditti, Steinmann, et al. 2014). Since Rozellomycota species are pathogens of other fungi, they may also present another indirect mechanism of lowering productivity by removing other beneficial fungal species. It is also suggested that Rozellomycota feed via phagocytosis in contrast to other fungi feeding via osmotrophy, thus they might have low sensitivity to SIP since they do not assimilate ^13^C–^15^N‐labelled amino acids. Furthermore, a recent study found that the abundance of Rozellomycota species was substantially higher in disease conducive soil (DCS) for tobacco cultivation (Gao et al. 2020). Our results also indicate that in high‐productivity soils, where the members of the phylum Rozellomycota have very low abundance, members of the Glomeromycota thrive instead and vice versa. Further functional genomic evidence is needed to better elucidate the potential role of Rozellomycota in blueberries affected by replant disease.

Most of the ^13^C–^15^N‐amino acid assimilation in the eukaryotic community appeared to be from an unclassified Mortierella sp. across all soil types possibly due to the presence of its trophic (vegetative) life stage primarily feeding by osmotrophy. Other fungal taxa may be present in the form of spores or not directly assimilating amino acids/glucose and remain undetected in the ‘active’ community fraction. In addition, active Oomycota (unclassified Pythiaceae spp.) species were abundant in low‐productivity soils from both farms but were almost absent in high‐productivity soils. Members of the Oomycota phylum are known pathogens of agricultural crops, nematodes, algae, zooplankton and fish. Some Pythium species can cause damping‐off disease in seedlings, root rot and other plant diseases, and are associated with replant disease in apple and ginseng species (Ontario Ministry of Agriculture, Food and Rural Affairs [OMAFRA] 2015; Somera and Mazzola 2022). They are fungi‐like organisms as well and can either form hyphae or live as very simple holocarpic thalli. Some species, though, are not predominantly pathogens; they are saprophytes but can act as opportunistic pathogens when the immune system of the host organism is compromised (Hassett et al. 2019; Klinter et al. 2019). Low‐productivity farm soils were also relatively enriched with ciliates. Although it is not clear how ciliates interact with plant rhizospheres, it has been shown that their diversity does not depend on the resident plant species unlike fungi and bacteria. Rather, the richness of ciliates in rhizosphere microbiomes may depend on the availability of inorganic N or the composition of rhizosphere bacterial communities (Acosta‐Mercado and Lynn 2004). This suggests predation may play a major role in removing beneficial microbes from the blueberry rhizosphere in low‐productivity soils.

Finally, the long‐read classification of eukaryotes indicated some surprising findings as well. Candida blankii is an emerging human pathogen with reduced susceptibility to antifungals (de Almeida Jr et al. 2018; Mirchin et al. 2022; Pinho et al. 2024). Perhaps the Candida strains we observe with our rRNA operon reads are emerging plant pathogens. Nadsonia starkeyii‐henricii has been found in forest soils (O'Boyle et al. 2018) and Sugiyamaella chiloensis was isolated from a fallen beech tree in Chile (Shi et al. 2021). None of these fungal species are reported as plant pathogens. Interestingly, Mucor lusitanicus has been reported as a biocontrol agent against fungal phytopathogens (Alejandre‐Castañeda et al. 2023) and the large signal in the resident/active community of the high‐productivity soils of Farm1 is consistent with that finding. Additionally, the rRNA operon reads associated with Podila verticillata in both the low productivity soils in Farm1 and in the active signal from Forest1 suggest that sampling at these sites may have disturbed ant mounds. Ants are common in New Jersey blueberry fields, and their occurrence sometimes indicates the presence of root‐colonising mealy bugs (Dysmicoccus vaccinii) (Pavlis 2009). It is worth highlighting that even though the long‐read classification of eukaryotes provided increased taxonomic resolution, several taxa were not detectable with this approach. This result indicates that different reference databases likely have different representative taxa, and misclassification is possible, especially for those raw reads demonstrating long BLAST alignments and lower percent identity. Currently, it is suggested to screen query rRNA operon sequences against more than one reference database for methodological validation of the taxonomic assignments of rRNA operons until more robust database entries can be generated.

In conclusion, our study has demonstrated the ability to profile ribosomal RNA operons from ^13^C/^15^N‐labelled DNA for increased resolution of active bacteria and eukaryotes associated with the rhizosphere soils of blueberry plants. Although we do not know the mechanisms by which these active groups of microorganisms affect blueberry productivity in situ, we have provided a rapid diagnostic tool for discerning high‐ and low‐productivity farm soils.

Author Contributions

J.J.P. and L.J.K. conceived and designed the experiments. S.M. performed DNA extractions, amplified the rRNA operons and created the sequence libraries with L.J.K. S.M. performed the data/statistical analysis in consultation with L.J.K. on the short‐read results. L.J.K. performed the data/statistical analysis on the long‐read results. S.M., J.J.P. and L.J.K. discussed the findings and interpreted the results. S.M. wrote the first draft. All authors read, edited and approved the final manuscript.

Funding

This work was supported by the Natural Resources Conservation Service, New Jersey Blueberry and Cranberry Research Council.

Consent

All authors of the manuscript have read and agreed to its content and are accountable for the accuracy and integrity of the manuscript.

Conflicts of Interest

The authors declare no conflicts of interest.

Supporting information

Figure S1: Example of highbush blueberry plants in rows exhibiting high productivity (normal growth) and low productivity (symptomatic of replant disease) where soil samples were collected. Figure S2: Schematic combining stable isotope probing (SIP) with sequencing of rRNA operon amplicons via the Oxford Nanopore MinION device. Figure S3: Example of QA/QC for SIP when testing both ‘light’ (^12^C/^14^N) and ‘heavy’ (^13^C/^15^N) bands for amplification depending on organic substrate amendment. Figure S4: Plot of alignment length versus pairwise identity from 10,000 randomly selected bacterial rRNA operon reads following classification using EZBioCloud 16S rRNA database (A). Overall, > 2 million reads were obtained for this study. Relative abundance of the phyla detected in this study is shown in (B). Figure S5: Bubble plot of the top bacterial taxa based on relative abundance for the ^12^C/^14^N band (resident) and the ^13^C/^15^N band (active) from classification by the EZBioCloud 16S rRNA database. Red arrows indicate taxa enriched in low productivity (replant diseased) soils. Figure S6: Plot of alignment length versus pairwise identity from 10,000 randomly selected eukaryotic rRNA operon reads following classification using the UNITE database (A). Overall, > 3 million reads were obtained for this study. Relative abundance of the phyla detected in this study is shown in (B). Figure S7: Bubble plot of the top eukaryotic taxa based on relative abundance for the ^12^C/^14^N band (resident) and the ^13^C/^15^N band (active) from classification by the UNITE ITS database. Red arrows indicate taxa enriched in low productivity (replant disease‐affected) soils. Figure S8: Plot of alignment length versus pairwise identity from 10,000 randomly selected bacterial rRNA operon reads following classification using the bacterial rRNA operon database.Figure S9. Plot of alignment length versus pairwise identity from 10,000 randomly selected bacterial rRNA operon reads following re‐classification using the eukaryotic rRNA operon database. Figure S10: Enriched taxa detected by DESeq2 analysis for the bacterial rRNA operon reads following classification using EZ BioCloud database. Those taxa exhibiting log2‐fold‐change > 2 and adjusted values p < 0.05 are indicated with red dots. Figure S11: Enriched taxa detected by DESeq2 analysis for the eukaryotic rRNA operon reads following classification using the UNITE database. Those taxa exhibiting log2‐fold‐change > 2 and adjusted values p < 0.05 are indicated with red dots.

Table S1: Bacterial (A) and eukaryotic (B) universal rRNA operon primers attached to ONT (Oxford Nanopore Technologies Ltd.) barcode sequences.

Data S1: emi70280‐sup‐0003‐Supinfo1.xlsx.

Data S2: emi70280‐sup‐0004‐Supinfo2.docx.