Reclaimed desert habitats favor entomopathogenic nematode and microarthropod abundance compared to ancient farmlands in the Nile Basin

Characterizing entomopathogenic nematode (EPN) biogeography with a goal of augmentation and conservation biological control requires fine-scale taxonomic resolution, because closely related EPN species can exhibit divergent phenotypes for key properties such as habitat adaptation and insect host specificity. Consequently, we employed high throughput genome sequencing (HTS) to identify and compare EPNs and natural enemies of EPNs in 58 citrus orchards in 2 ecoregions in Egypt (El Beheira and Al Qalyubia governorates). We designed improved primers targeting the ITS2 rDNA to discriminate EPN species and used pre-reported primers targeting D2-D3 region for soil microarthropods. Five EPN species (Heterorhabditis bacteriophora, H. indica, H. taysearae, Steinernema glaseri, and S. scapterisci) and one steinernematid not represented in Genbank databases were detected. This is the first report of S. scapterisci and possibly the unknown (perhaps undescribed) species in Egypt. Only heterorhabditid species, dominated by H. indica, were detected in the reclaimed, sandy desert soils of El Beheira governorate. In the fine textured, ancient farming lands of the Nile delta all six species were detected, but at lower frequency and abundance. Microarthropod family richness (P = 0.01) and abundance (P = 0.001) was higher in the reclaimed lands than in the Nile Delta. Soil clay content, pH and elevation explained significant variation in the mite community structure. Population density of H. indica, the only EPN found consistently and at high abundance in El-Beheira, was inversely related to abundance of species in the nematophagous mite family Rhodacaridae.

Research on entomopathogenic nematodes (EPN) in Egypt started in the 1970s, and focused heavily on imported, non-indigenous species (Abd-Elgawad, 2017). Surveys to isolate and identify indigenous EPNs began two decades later (Shamseldean and Abd-Elgawad, 1994). Inconsistent efficacy by expensive EPN products hinders their use by the Egyptian farmers, suggesting a need for further exploration to identify species which are adapted to North African conditions and best suited to infect local insect pests (Campbell and Gaugler, 1993;Koppenhöfer et al., 1998;Simões and Rosa, 1996). To this end, a comprehensive survey that employed soil baiting with sentinel insects to recover EPN from 1,000 samples from the Nile Delta and Mediterranean Sea coast, Nile valley, Upper Egypt, and Sinai Peninsula revealed just three species, Heterorhabditis indica Poinar (Rhabditida: Heterorhabditidae), Karunakar & David, Steinernema abbasi (Rhabditida: Steinernematidae) Elawad, Abbas & Hague, and S. carpocapsae Weiser (Rhabditida: Steinernematidae) (Abd-Elbary et al., 2012). However, limitations of this common EPN survey method are reflected by results of other surveys that, to date, have identified 10 described EPN species in Egypt (Table 1).
Soil baiting can fail to detect EPN for any number of reasons including incompatible host status of the sentinel (Nguyen and Smart, 1991), competition between EPN and other organisms (Duncan et al., 2007;Wu et al., 2018), and phased infectivity exhibited by many EPN populations (Baiocchi et al., 2017;Shields, 2015). Direct observation through microscopy is more likely than baiting to detect EPN but is time consuming. Both methods require considerable expertise to identify species and suffer from a lack of definitive morphological features. By contrast, quantitative PCR is both sensitive and accurate and is used increasingly to detect and identify EPN (Campos-Herrera et al., 2012). The primary limitation of qPCR is that it detects only those species matching the primer/probes that are used. Thus, for surveys, metagenomic methods provide the most reliable tool, regardless even of whether a species is represented in databases such as Genbank (Dritsoulas et al., 2020).
Citriculture has an increasing socio-economic importance in Egypt, but is subjected to considerable yield loss caused by insect pests such as such as the Mediterranean fruit fly Ceratitis capitata, Wiedemann (Diptera: Tephritidae) and hairy rose beetle Tropinota squalida Tropinota squalida, Scopoli (Coleoptera: Scarabaeidae) (Abd-Elgawad, 2020). Isolating EPN strains from the main citrus producing governorates in Egypt for practical use in biocontrol programs could provide more effective matching of nematodehost-environment in the crop (Abd-Elgawad, 2017, 2020. Al Qalyubia and El-Beheira governorates are major citriculture regions; however, from 64 sites in these governorates, just one sample from orange groves and eight from non-citrus groves were positive for unidentified EPN species (Shamseldean and Abd-Elgawad, 1994). Two strains of Heterorhabditis bacteriophora Poinar (Rhabditida: Heterorhabditidae) were also recovered from mango groves in El-Beheira and Al Qalyubia (Abd-Elgawad and Nguyen, 2007). The apparent depauperate state of EPN in Al Qalyubia and El-Beheira orchards contrasts with several surveys in other countries where citrus and other fruit orchards tended to support more abundant and diverse EPN populations than other agro-ecosystems and natural areas (Campos-Herrera et al., 2013Steyn et al., 2017;Tarasco et al., 2015). Hence, we employed metabarcoding to detect EPN in citrus orchards in Al Qalyubia and El-Beheira. We included primers in some samples as a first approach to identify microarthropods that might differentially modulate the abundance of EPNs in the two regions.

Materials and methods
Soil sampling and chemical analysis  . A total of 60 composite samples were collected for nematode analysis. All samples were kept in polyethylene bags, labeled, and transferred to the laboratory for nematodes extraction. Samples were gently sieved through a 4 mm aperture screen to remove gravel and to separate roots from soil. For each sample, 300 ml of the sieved soil were put in a stainless steel bowl, covered with tap water, and agitated before pouring the soil suspension through a 20-mesh (850 µm aperture) sieve held over a second bowl. Nematodes were extracted from each soil sample by sucrose (545 g sugar L −1 ) centrifugation after processing through 325-and then 500-mesh sieves (Jenkins, 1964). Nematode suspensions were concentrated and 100% alcohol was added to each sample. The nematodes could settle in test tubes overnight at 4°C. Thereafter, most of the water with alcohol was evaporated and samples were transferred to 1.5 mL Eppendorf tubes and stored at −20°C until sent to the laboratory in Florida for DNA extraction (Campos-Herrera et al., 2011).

DNA extraction
Based on our experience that samples in ethanol yield low DNA, following centrifugation and aspiration of excess ethanol, the tubes were refilled with 1xPBS (phosphate buffer saline) and incubated overnight at 4°C. After a second centrifugation and aspiration of excess PBS, DNA was extracted with the DNeasy® PowerSoil Kit (Qiagen).

Library preparation
Prior to library preparation process, DNA concentrations from each sample were measured using the Qubit ® dsDNA High Sensitivity Assay Kit (ThermoFisher Scientific, USA). Libraries were constructed for two groups, (i) 58 libraries targeting nematodes and (ii) 16 libraries targeting soil microarthropods. Microarthropod libraries were from eight randomly selected samples from each governorate.
According to Illumina protocols, library preparation comprises four parts: (i) amplicon PCR, (ii) amplicon PCR cleanup, (iii) index PCR, and (iv) index PCR cleanup. Samples were standardized at 5 ng/ml DNA concentration. For the EPN, samples were amplified with the following conditions: initial denaturation 95°C for 3 min, 25 cycles of denaturation at 98°C for 30 s, annealing at 56°C for 30 s, elongation at 72°C for 60 s, and terminal elongation at 72°C for 10 min. For microarthropods the conditions were initial denaturation 95°C for 3 min, 25 cycles of denaturation at 98°C for 30 s, annealing at 55°C for 30 s, elongation at 72°C for 60 s, and terminal elongation at 72°C for 10 min. For all libraries a single 25 μ L PCR reaction containing 2.5 μ L of template of 5 ng/μ L (12.5 ng total), 12.5 μ L of 2x KAPA HiFi HotStart ReadyMix (KAPA biosystems), 1 μ L of each 10 μ M overhang primer, 8 μ L of 10 mM Tris pH 8.5. Validity and reliability of PCR reactions was tested with positive controls using DNA extracted from laboratory culture of the nematodes Steinernema feltiae and Heterorhabditis bacteriophora while nega tive controls included nuclease-free water instead of DNA template. PCR products were verified on 1% agarose gels while all PCR products were purified with 1.0 × Agencourt AMPure XP beads (Beckman Coulter, Brea, CA) and eluted in 50 μ L of 10 mM Tris pH 8.5. For index PCR, purified amplicons were used as template for an eight-cycle amplification to add dual-index bar codes, P5 and P7 Illumina sequencing adapters using the Nextera XT Index Kit [FC-131-1004] for EPN and XT Index Kit [FC-131-1001] for microarthropods (Illumina, San Diego, CA, USA). The index PCR conditions were initial denaturation at 95°C for 3 min, 8 cycles of denaturation at 98°C for 30 s, annealing at 55°C for 30 s, and elongation at 72°C for 30 sec and a terminal elongation at 72°C for 10 min. Each 50 μ L PCR reaction tube contained 5 μ L of template, 25 μ L of 2x KAPA HiFi HotStart ReadyMix (KAPA biosystems), 5 μ L of Index Primers (N7XX), 5 μ L of Index 2 Primers (S5XX).
Index PCR products purified with 1.1× magnetic beads, eluted in 25 μ L, and quantified using a fluorometric quantification method that uses dsDNA binding dyes. Concentrated final libraries diluted using 10 mM Tris pH 8.5 to 4 nM. All the 5 μ L aliquots of diluted DNA from each library was mixed in a single pooling library. The pooling library was sequenced using MiSeq platform 2 × 300 bp paired-end Illumina at the Interdisciplinary Center for Biotechnology Research (ICBR) of University of Florida.

Bioinformatics
ICBR delivered raw data in fastq format which were demultiplexed and separated into respective sample identification codes. FASTQC v0.11 (Andrews et al., 2015) was used for quality assessment of each read, and then all the quality information was combined into a single document using MULTIQC (Ewels et al., 2016). ITS 2 and LSU D3-D5 amplicons of ribosomal DNA was used for the nematode and microarthropods identification respectively. The R1 and R2 reads combined and de-replicated with the ASVbased approach, using DADA2 algorithms, through QIIME2 v2019.4 (Callahan et al., 2016). ended up to a length of 350-430 bp for nematodes and 500 bp for microarthropods. Count tables were generated by mapping ASVs and assigning taxonomy. All the non-redundant nucleotide sequences from NCBI GenBank were combined to generate a standalone database for taxonomy assignment (ftp://ftp.ncbi. nlm.nih.gov/blast/db/nr.gz) through a command-line tool (BLAST +) which was integrated directly into the workflow to run BLAST.

Genetic analysis of identified ASV
Phylogenetic analysis was carried out in MEGA 10.0.5 software. Each of the identified EPN species was evaluated using a unique tree derived from the metabarcoding process. Prior to constructing trees, sequences were aligned using ClustalW alignment method on default settings. The evolutionary history was inferred by using the Maximum Likelihood method and Tamura-Nei model while the robustness of clades of the trees was assessed using 1000 bootstrap replications.

Statistical analysis
Regional differences in soil properties and differences in sites with or without EPN were evaluated by t-test. Relationships between soil properties and H. indica were measured by Spearman rank correlations. H. indica abundance in sites with or without specific mite families were compared using Kruskal-Wallis Test. Multivariate analyses were performed using the software R (R Development Core Team, 'Vegan' package). Principal component analysis (PCA) was used to reveal acari mite and soil properties that contribute to the total spatial variability of the two ecoregions. Detrended Canonical Correspondence Analysis (DCCA) was employed to determine the heterogeneity of the system. A value below 3.0 suggested that the community is homogeneous; therefore, Redundancy Analysis was used as the most appropriate constrained analysis, applying Monte Carlo permutation for significant environmental variables at the 0.05 level. RDA results were graphically presented with bi-plot scaling (R, 'Vegan' Package). Correlations between Acari mites against soil properties were performed in R studio using the corrplot function. Non-parametric test Kruskal-Wallis was employed to evaluate differences in occurrence of all microarthropod families detected in the two ecoregions.

Results
The high-throughput sequencing produced two datasets, one based on ITS2 targeting nematodes and the other on the D3-D5 region of 28 S rDNA targeting microarthropods. The ITS2 revealed 8,208,384 reads of which 43% (3,522,968) passed the quality filters and were denoised, merged and characterized as non-chimeric. 44,483 unique amplicon sequence variants (ASVs) were recovered from forty-nine phyla, with 2.1% (946 ASVs) assigned to phylum Nematoda. By setting a threshold of 80% coverage, sixty seven percent (639) of the unique ASVs belong to 19 nematode families. Finally, 22 ASVs were identified as entomopathogenic nematode. The D3-D5 dataset yielded 9,637,729 reads that reduced to 3,627,990 (37%) after filtering, denoising, merging and chimera removal. The total number of unique ASVs was 2324 from 27 phyla, of which 483 belong to arthropods. Setting a threshold of 70% coverage, revealed 221 ASVs from 28 families in the class Arachnida.
Entomopathogenic nematode isolates of Hetero rhabditis and Steinernema were detected in 41% of the samples (Fig. 3). In all, 18 of 30 sites in El-Beheira (60%) were positive for EPN (all heterorhabditids), while EPN were detected in 8 of 28 sites (29%) in Al-Qalyubia. Steinernematid ASVs identified as Steinernema glaseri, Steinernema scapterisci, and Steinernema sp. were detected in 2 sites each in Al-Qalyubia (Fig. 3b). Heterorhabditis indica was identified in 14 sites (24%), Heterorhabditis taysearae in 4 sites and Heterorhabditis bacteriophora in 3 sites (Fig. 3b). Each heterorhabditid species was found in both regions. H. indica was the only species which occurred more frequently in one of the governorates -43% of sites in El-Beheira vs 4% in Al-Qalyubia (P = 0.000). Al-Qalyubia supported all six detected EPN species, but at low population density (average EPN ASV reads per positive site = 548) compared to El-Beira (9284) (P = 0.02), due primarily to large numbers of H. indica at El Beheira.
The soils in El-Beheira governorates were consistently coarser textured with less organic matter than those in Al-Qalyubia (Table 2). Soil properties did not differ significantly in sites that were either positive or negative for EPN. The abundance of H. indica was greatest in sites with coarse soil texture and higher pH, but lower electroconductivity. However, in El-Beheira (where all but one site positive for H. indica occurred), there were no significant relationships between H. indica and soil properties.
Microarthropods comprising 28 families were identified from the eight sites in each governorate with 13 detected in Al-Qualubiya and 23 in El-Beheira. Family richness in Beheira (9.75) was more than twice (P = 0.01) that in Qualubiya (4.5). The overall microarthropod abundance was greatest in El-Beheira according to paired t-test of the ASV from the 28 families (P = 0.001). Kruskal-Wallis tests detected greater abundance in Beheira than in Qualibiya of the Ascidae, Tydeidae, Rhodacaridae, Ologamasidae, Oehserchestidae, Ereynetidae, Eupodidae. A principal component analysis of soil properties showed two very different ecoregions with respect acari mite communities (Fig. 4A). Redundancy analysis (RDA) identified clay and OM as associated with the communities of acari mites (P < 0.1) (Fig. 4B). Fewer H. indica were detected in sites with abundant Rhodacaridae species (P = 0.006). Spearman's correlations suggest that clay content, pH and elevation may modulate the occurrence of these species (Fig. 5).

Discussion
The factors regulating EPN species occurrence and abundance remain poorly understood despite an ever-expanding catalogue of EPN biogeography (Campos-Herrera et al., 2013, 2019aGarcia Del Pino and Palomo, 1996;Mráček et al., 1999;Tarasco et al., 2015;Valadas et al., 2014). While habitat biological and abiotic complexity obscures key processes affecting EPN, enhanced accessibility of metagenomic tools has vastly increased the resolution of soil food web characterization. The inventory created here of EPNs in the citrus orchards of two Egyptian ecoregions, demonstrates the enhanced capacity of metagenomic tools to detect, identify and, to some extent, quantify EPN across habits (Dritsoulas et al., 2020). The detection frequency in these samples (41%) far exceeds those of previous Egyptian surveys where EPNs were recovered from 9.5% of 661 soil samples (Shamseldean and Abd-Elgawad, 1994) and from 16% of 1,120 soil samples (Abd-Elbary et al., 2012). Critically, metabarcoding detected known and undescribed species whose relevant sequences are not registered in the GenBank databases. When blasting the sequencing output, each ASV is identified as an organism regardless of the proximity of the query to the reference sequence. The query sequences are not always identical to the reference sequences, but may differ by a few nucleotides, a phenomenon described by Porazinska et al. (2010) as generally conforming to a head (many identical) Figure 1: Phylogenetic relationships of ASVs identified in the genus Steinernema based on sequencing reads of the ITS-2 region as inferred by the Maximum Likelihood method and Tamura-Nei model. The percentage of replicate trees in which the associated taxa clustered together in the bootstrap test (1,000 replicates) are shown next to the branches. All reference sequences are indicated by the NCBI accession numbers. The reference sequences were employed to confirm the right identification of each ASV in the phylogram. C. elegans was used as a global outgroup. tail (few variants) pattern. While most sequences detected in our samples followed this pattern, we also detected ASVs identified with low affinity (< 79.5%) as Steinernema neocurtillae. Phylogenetic analysis and the derived phylogenetic tree positioned these ASV within the genus Steinernema. These findings support the need of further investigation at the indicated sites in order to isolate and describe EPN conforming to these sequences.
The constructed phylogenetic topology was critical also because it showed that multiple ASV characterized as EPN in Genbank databases actually belong to distant, unrelated nematode families (data not shown). Moreover, four ASVs of the sequencing output were identified as H. sonorensis and one was identical to both H. sonorensis and H. taysearae species. Alignment of those two reference sequences (EF043443 and FJ477730 respectively), differed by just one nucleotide in ITS1 and ITS2 region of 869 bp length. There is no difference in ITS2 region that was targeted by the primers used here. H. taysearae has been reported from Egypt (Shamseldean et al., 1996) and other African countries such as Benin (Godjo et al., 2018), Kenya (Nyasani et al., 2008) and South Africa (Nyasani et al., 2008) (Table 1). S. scapterisci is an autochthonous species of South America and was imported to other countries as a biological control agent against orthopterans, primarily in USA for application in turf grass and pastures to control Scapteriscus spp. mole crickets. It is now commonly Figure 2: Phylogenetic relationships of ASVs identified in the genus Heterorhabditis based on sequencing reads of the ITS_2 region as inferred by the Maximum Likelihood method and Tamura-Nei model. The percentage of replicate trees in which the associated taxa clustered together in the bootstrap test (1,000 replicates) are shown next to the branches. All reference sequences are indicated by the NCBI accession numbers. The reference sequences were employed to confirm the right identification of each ASV in the phylogram. C. elegans was used as a global outgroup.  Notes: Mean ± standard error mean, maximum and minimum values and the significant regional differences are represented for each component (df = 58). Hi column represents significant differences in soil properties between the sites that Heterorhabitis indica = Hi occurs and did not occur (df=12). Differences between in all cases were evaluated by non-parametric Wilcoxon Test. *P < 0.01; **P < 0.001; ***n.s. = non-significant.
detected in citrus orchards (Campos-Herrera et al., 2013 representing long-distance dispersal of the nematode in infected crickets (Parkman et al., 1993).
To date, there are no records of biological control programs in Egypt using these nematodes: however, there are local institutions which maintain S. scapterisci in culture, supporting the feasibility of its detection in Al-Qalyubia. The limited occurrence and low abundance of EPN in this survey precluded detecting relationships between most EPN and the different habitat/foodweb properties of the ancient agricultural fields of the Nile Delta compared to those in land reclaimed from desert just since the early 1950s, some in this survey as recently as 30 years ago. Only H. indica was found to occur more frequently and abundantly in the sandy soils of El-Beheira and only heterorhaditids were detected in its orchards. A similar survey detected no EPN in the natural Negev Desert soils, but recovered heterorhabditids in the rhizospheres of irrigated fruit trees there (Glazer et al., 1991). A number of other studies have reported an association of heterorhabditids with sandy and/or coastal soils in both temperate and tropical regions worldwide: Griffin et al. (2000) in Indonesia (Amarasinghe et al., 2012), in Sri Lanka, , in Britain, (Hara et al., 1991) in Hawaii, with other reports showing the association for H. indica, specifically, in Guadeloupe (Auléon et al., 2006), Japan (Yoshida et al., 1998), andFlorida (Campos-Herrera et al., 2013).
The occurrence of soil microarthropods also differed between the two ecoregions, with El-Beheira richer in microarthropods than Al-Qualubiya. Clay and elevation had the strongest relationship to microarthropod communities here and in some previous reports (Benckiser, 1997;Maraun et al., 2013;Marian et al., 2018). Although the relatively small changes in altitude in this survey suggest a relationship with a hidden variable such as water table depth (Campos-Herrera et al., 2013), soil microarthropods depend strongly on soil texture as they need pore space for all of their activities. Among the seven families found to be significantly more abundant in Beheira (Ascidae, Tydeidae, Rhodacaridae, Ologamasidae, Oehserchestidae, Ereynetidae, Eupodidae), most are nematophagous (Epsky et al., 1988;Karagoz et al., 2007;Santos and Whitford, 1981) including species of Rhodacaridae shown to be strongly inversely related to H. indica here.
This project has detected for the first time evidence of two new EPN species in Egypt, one described and represented in Genbank, the other not in Genbank  . P < 0.1; *P < 0.05; **P < 0.01; ***P < 0.001. Figure 5: Non-parametric Spearman's correlations ( N = 16) between abiotic factors and those acari families with most prominent differences between the two regions (Table 3). Positive correlations are displayed in blue and negative correlations in red. Color intensity and the size of the circle are proportional to the correlation coefficients. Significant P-values are shown.
and potentially undescribed. It more fully characterized the distribution and abundance of EPN in the citrus orchards of the Nile Delta and the reclaimed desert regions using the most sensitive methods currently available. As such the data are more comparable to those in other recent surveys employing molecular methods, rather than sentinel baiting, to survey EPN in citrus orchards in different parts of the world (Campos-Herrera et al., 2013, 2019aDritsoulas, 2020;Dritsoulas et al., 2020). Only by standardizing methodology will EPN biogeography and the mechanisms regulating occurrence and abundance be accurately revealed. Recent work has indicated that sucrose centrifugation combined with molecular identification and quantitation is highly efficient for characterizing food web components such as nematophagous fungi, parasitic bacteria, and microarthropod predators capable of modulating EPN populations (Dritsoulas and Duncan, 2020;Dritsoulas et al., 2020;Pathak et al., 2017). Improved understanding of how food webs function in different habitats is necessary to discover cultural practices that can enhance biological control by introduced or naturally occurring EPN.