Characterization of the phenotypic and genotypic tolerance to abiotic stresses of natural populations of Heterorhabditis bacteriophora

Entomopathogenic nematodes are effective biocontrol agents against arthropod pests. However, their efficacy is limited due to sensitivity to environmental extremes. The objective of the present study was to establish a foundation of genetic-based selection tools for beneficial traits of heat and desiccation tolerance in entomopathogenic nematodes. Screening of natural populations enabled us to create a diverse genetic and phenotypic pool. Gene expression patterns and genomic variation were studied in natural isolates. Heterorhabditis isolates were phenotyped by heat- and desiccation-stress bioassays to determine their survival rates compared to a commercial line. Transcriptomic study was carried out for the commercial line, a high heat-tolerant strain, and for the natural, low heat-tolerant isolate. The results revealed a higher number of upregulated vs. downregulated transcripts in both isolates vs. their respective controls. Functional annotation of the differentially expressed transcripts revealed several known stress-related genes and pathways uniquely expressed. Genome sequencing of isolates with varied degrees of stress tolerance indicated variation among the isolates regardless of their phenotypic characterization. The obtained data lays the groundwork for future studies aimed at identifying genes and molecular markers as genetic selection tools for enhancement of entomopathogenic nematodes ability to withstand environmental stress conditions.

www.nature.com/scientificreports www.nature.com/scientificreports/ of protecting the nematodes from harmful environmental factors by development of formulation technologies aimed at stabilizing the nematodes during production, distribution and storage 24 . The second approach is screening of wild populations for resistant strains expressing the beneficial traits. An EPNs survey revealed the presence of natural populations of the genus Heterorhabditis in arid and semiarid regions in Israel, characterized by high temperatures and dry soils 25,26 . Isolation of natural populations from arid and warm environments indicates the enhanced abilities of some populations to survive harsh environmental conditions 25 . Based on the screening of natural populations, the leading approach to EPN enhancement seems to be a genetic one, through selective breeding of strains with beneficial traits, genetic manipulation (mutagenesis) and genetic engineering 18,27,28 . Advanced genetic techniques, such as gene transformation and induction of mutations using RNA silencing have been studied, and enhancement of beneficial traits for environmental stress tolerance has been reported 25 . However, the genetic study of EPNs requires more advanced information and tools. In a previous study we analyzed the transcriptome expression in Steinernema spp. with varied stress-tolerance capabilities and revealed a number of genes and metabolic pathways that are differentially expressed under stress conditions 29 . Still, further transcriptomic and genomic studies under stress conditions are expected to reveal specific genes involved in the stress response, and enable the establishment of expression markers for the enhancement or degradation of beneficial traits.
The main objective of the current study was to establish an infrastructure for the identification of molecular markers associated with heat and desiccation stress in these organisms. This was done by characterizing natural EPN populations for their tolerance to heat and desiccation stress and comparing them with a well-studied commercial line. Based on the phenotypic characterization, isolates with diverse tolerance to the studied stress conditions were chosen for further examination of gene-expression patterns and genomic variations. We hypothesize that: 1. Natural populations of EPNs differ in their tolerance to heat and desiccation stresses. This difference may be related to the characteristics of the isolation site. 2. Gene-expression pattern and genomic sequence differ significantly between isolates with high tolerance vs. low tolerance to environmental stress. To address these hypothesis we followed the following objectives: 1. Isolation of natural populations of EPNs from different habitats and climatic regions in Israel. 2. Phenotyping of heat and desiccation tolerance in these natural isolates, compared to a characterized commercial line. 3. Characterization and comparison of gene-expression patterns in high-tolerance isolates and low-tolerance isolates, using a commercial line as a reference strain for their tolerance capabilities. 4. Identification of genomic variation between high-and low-tolerance isolates, compared to a commercial line.

EPN isolation and habitat characterization. The occurrence of natural populations of EPN was
assessed as recovery frequency (number of positive samples/number total samples) and abundance (number of positive sites/number total sites) expressed as percentage 30 . EPNs were recovered from 8 out of 136 soil samples, 5.8% recovery frequency, from 34 sites, 23.5% abundance. In addition to the overall recovery and abundance of EPN, the relation with soil type, soil function, climate, soil water content and soil temperature was studied. None of the obtained habitat parameters were significant in relation to EPN recovery (χ 2 (1) =0.99 (P = 0.32), χ 2 (1) =0.75 (P = 0.38), χ 2 (2) =1.34 (P = 0.5), χ 2 (1) =0.22 (P = 0.63) and χ 2 (1) =0.01 (P = 0.91) respectively). EPN recovery in relation to soil type, soil function and climate conditions is described as positive samples out of total number of samples in Table 1. Six out of eight positive samples were isolated from cultivated soils of groves and orchards. The two isolates from non-cultivated soils were identified as Steinernema feltiae, whereas, all other six isolates were identified as Heterorhabditis spp. (Table 1, Fig. 1C). Soil profiling revealed that 50% of the positive samples were from sandy soils and 50% were from clay soils (Table 1). Climate conditions in the different sites, varied from Mediterranean in the northern and coastal regions, and arid or semi-arid in the southern regions ( Table 2, Fig. 1A). Relating to three main parameters; climate, soil type and soil function, resulted in eight habitat types. The relative frequency of EPN recovery in relation to those types is described in Fig. 1B. The highest relative recovery frequency (16.6%) was recorded in cultivated habitat with sandy soil and arid climate. Followed by non-cultivated habitat with clay soils and Mediterranean climate (8.3%), which was the only record of EPNs recovered from non-cultivated soils in the present study. Habitat with semi-arid climate and sandy soils did not yielded EPN recovery neither in cultivated soils nor in non-cultivated soils. The water content in the soil samples ranged from 89% in samples collected during January and February, and 10% in samples collected during May. The average water content in soil samples was 21.7% ± 1.8%. EPNs were isolated from soils with minimal water content of 10% and maximal water content of 47%. The average temperature in the soil samples was In addition to the isolates described above, four isolates from previous surveys conducted in Israel 21,25,31 and two non-endemic isolates, from USA and Germany were studied. Characterization of collection sites from the present study and from previous studies including isolation sites of all Israeli isolates is presented in Table 2.
Phylogenetic characterization of studied EPN isolates. Fifteen nematode isolates were characterized at the species level based on the ribosomal ITS region. Electrophoresis revealed 600-800 bp PCR product. BLAST comparison revealed that 5 out of 12 isolates are H. bacteriophora, 5 are H. indica and 2 are S. feltiae. Phylogenetic characterization and relations between the different isolates presented in Fig. 1C and Table 1. A total of 13 Heterorhabditis isolates were used in the research, including the commercial line, EN-01. The two isolates belonging to the genus Steinernema were not further studied in the present study. The studied nematodes isolates were cultured in the lab during the research period, and freshly emerged IJs were used for all bioassays and molecular studies.

Determination of Mean Tolerated Temperature (MT°C50) of commercial line.
Since the commercial line EN-01 served as the reference line throughout the study, its mean tolerated temperature was determined and further used for comparison of the studied isolates. Survival rates of IJs of the commercial line, EN-01, decreased with increasing temperature. Survival rates of the control group at 25 °C were 99% ± 0.4. Gradient treatments between 37-39.4 °C had survival rates between 94% ± 0.8% -9% ± 3%. Tukey's HSD test revealed significant differences in the survival rates in temperatures higher than 37.7 °C (P = 0.0001). Temperature above 38.2 °C caused more than 50% ± 5% mortality (Fig. 3A). Probit analyses determined the MT 50°C of EN-01 as 38.8 °C.
Determination of Mean tolerated Water Activity (MW50) of commercial line. Since the commercial line EN-01 served as the reference line throughout the study, its mean tolerated active water was determined and further used for comparison of the studied isolates. Survival rates of IJs of EN-01 decreased with increasing PEG concentrations at the range of 20-40% v/v equivalent to active water of 0.98-0.96 a w . Survival rates of the control group in water were 99% ± 0.3. Tukey's HSD test revealed significant differences in the survival rates in treatments with PEG percentage higher than 33% v/v (P = 0.03), equivalent to active water lower than 0.97 a w . The www.nature.com/scientificreports www.nature.com/scientificreports/ survival rates in treatments of 0-30% v/v PEG were between 99% ± 0.3%-70% ± 9%, whereas survival rates in treatment of 35% v/v PEG were 41% ± 10% (Fig. 3B). MW 50 was determined as 33.88% v/v PEG.
Comparative heat tolerance in the studied isolates. The mean tolerated temperature determined for the commercial line, EN-01 was used as a standard to compare the survival rates of the isolates in this study. The survival percentages of IJs of the different species and isolates following exposure to heat stress at 38.8 °C is presented in Fig. 4. The results show significant differences (F 11, 117 = 7.06, P < 0.0001) in the survival rates between the species and among isolates belonging to the same specie. The mean survival percentage was 60% ± 8% and 38% ± 5% among H. indica and H. bacteriophora isolates respectively. Among H. indica, Mop, RA and ENL isolates displayed the highest tolerance, and Gilat isolate displayed the lowest tolerance (Fig. 4A). Among H. bacteriophora, the commercial line, EN-01 together with the wild type isolate, KH, displayed the highest tolerance. The other wild-type isolates Grofit, Magen and Porat displayed the lowest tolerance (Fig. 4B).  www.nature.com/scientificreports www.nature.com/scientificreports/ Comparative desiccation tolerance in the studied isolates. The mean tolerated active water determined for the commercial line, EN-01, was used as a standard to compare the survival rates of the isolates in this study. The survival percentages of IJs of the different species and isolates following exposure to desiccation stress of 0.97 a w is presented in Fig. 4. The results show significant differences between the two species and among isolates belonging to the same specie (F 11, 132 = 5.63, P < 0.0001). The mean survival rates of H. indica and H. bacteriophora were 40% ± 7% and 34% ± 13% respectively. Among H. indica, RA displayed the highest tolerance, mire and HPL displayed the lowest tolerance (Fig. 4C). Among H. bacteriophora, KH and Porat displayed the highest tolerance, whereas Grofit and Magen showed the lowest tolerance (Fig. 4D). The commercial line, EN-01 did not differ significantly from the highest tolerance isolates (F 4, 55 = 5.96, P = 0.14, P = 0.08 for differences with KH and Porat respectively).  Survival rates under increasing gradient temperature. IJs were exposed to adaptation phase of 3 hours at 35 °C, recovery phase of 1 hour at 25 °C and stress phase of 4 hours at gradient temperatures. Survival rates were determined following recovery phase at 25 °C over-night. (B) Survival rates under different PEG concentrations (v/v %). IJs were exposed to adaptation phase of 72 hours at 10% PEG solution, followed by stress exposure of 16 hours at gradient PEG dilutions. Error bars represent standard error for each mean and different letters express significantly different means (P < 0.05). Error bars represent standard error for each mean and different letters express significantly different means (P < 0.05).
www.nature.com/scientificreports www.nature.com/scientificreports/ Transcriptome mapping. In order to study gene expression patterns under stress condition of heat compared to control conditions, RNA samples were obtained from H. bacteriophora isolates with different heat tolerance phenotype. Transcriptome studies using RNA-Seq of the libraries of two isolates of H. bacteriophora IJs at early stages of exposure to heat stress and parallel control representing four libraries with three biological replications. The studied isolates were chosen according to the phenotype of survival under heat stress. EN-01 was chosen as a high heat-tolerance isolate and a commercial line and Grofit was chosen as a low heat-tolerance isolate (Fig. 4B). RNA-seq produced approximately 503 million raw reads. The reads were mapped to the reference genome of H. bacteriophora, IL3, homogenous inbreed line of the commercial line, EN-01. The average mapping percentage was 87% (Table 3).
Principle component analysis (PCA) revealed four clusters representing heat treatment and control for the two studied isolates (Fig. 1S). Each cluster is comprised of three biological replicates of the specific treatment. Results show distinctive clusters of the two heat groups (    www.nature.com/scientificreports www.nature.com/scientificreports/ Differential expression patterns. During early stages of exposure to heat stress of 35 °C compared to control at 25 °C, a large fraction of transcripts were up-regulated among high tolerance strain (EN-01) and low tolerance isolate (Grofit) (Fig. 5A). 378 transcripts were up-regulated and 231 transcripts were down-regulated in the high tolerance strain, the commercial line EN-01. In the low tolerance isolate, Grofit, 425 transcripts were up-regulated and 161 transcripts were down-regulated. 230 transcripts were commonly up-regulated and 50 transcripts were commonly down-regulated in both isolates (Fig. 5B). As mentioned according to the PCA plot (Fig. 1S), comparison of both isolates under heat treatment, showed similar expression patterns, however, 45 transcripts were identified as differentially expressed (2FoldChange&padj<0.05).
Putative functional classification using gene ontology (GO enrichment) and KEGG pathway analysis. All the differentially expressed transcripts (2FoldChange&padj<0.05) in the comparisons of heat and control treatments within the isolates, and heat treatments between the isolates were analyzed for GO enrichment. Differentially expressed transcripts between heat treatment and parallel control revealed 154 GO terms in the high tolerance strain, EN-01 and 115 GO terms in the low tolerance isolate, Grofit (Table S1). Further analysis revealed 29 enriched GO terms reported as related to stress response in nematodes, including response to topologically incorrect protein, response to heat, membrane organization, programmed cell death, phosphorylation, heat-shock protein activity, protein folding, fatty acid metabolism, response to abiotic stimulus and response to unfolded protein (Fig. 6). The significant difference of the enrichment of those stress-related GO terms between the high tolerance strain, EN-01 and the low tolerance isolate, Grofit is presented in Fig. 6 as -log10 (P value). Overrepresentation of significant (padj < 0.05) enriched GO terms that are related to stress was higher in the high tolerance strain, EN-01. Moreover, enriched GO terms of response to stress, detection of abiotic stimulus, response to temperature stimulus and thermos-sensoring were more significantly differential expressed in the high tolerance strain, EN-01, than the low tolerance isolate, Grofit. Since other genes might be related to stress, yet have not been characterized, further analyses was carried out using the KEGG mapper search tool in order to identify the key metabolic pathways and processes of which the differentially expressed transcripts were mapped to. The functional annotation revealed that 15 of the commonly up-regulated transcripts in both isolates are heat-shock proteins (HSP) related to longevity regulating pathways (Table S2). Two transcripts of the specific down-regulated transcripts of EN-01 were annotated as known heat stress-related proteins, glycerol kinase (GK) and fatty acid desaturase (FAD). Among the up-regulated specific transcripts of EN-01, three were annotated as heat-shock proteins and six were annotated as zinc finger proteins (ZFP), which are also related to stress response (Table 4). Among the specifically down regulated transcripts of the low thermos tolerance isolate, Grofit, two genes were annotated as trehalose coding gene (TRE) as part of the starch and sucrose metabolism (Table 4). Functional annotation of the up regulated transcripts of Grofit did not reveal any stress related pathways and proteins. GO enrichment of the differentially expressed transcripts between Grofit and EN-01 under heat treatment, did not yield significant stress-related GO terms for further evaluation.
Genome sequencing results and SNPs calling. In order to study genomic variation between isolates with varied degrees of stress-tolerance, three H. bacteriophora isolates were chosen for genomic sequencing and detection of SNPs. List of the isolates and their determined phenotype under heat and desiccation stress is presented in Table 5. Genomic DNA was sequenced successfully and resulted in 81,269,768, 105,659,854 and 77,291,447 paired end reads for Grofit, KH and Magen respectively. Sequences were mapped to the reference genome of H. bacteriophora homogenous inbreed of the commercial line EN-01, IL-3 genome. Mapping www.nature.com/scientificreports www.nature.com/scientificreports/ percentages are presented in Table 6. The variations among the isolates were identified in comparison with the reference genome, IL-3 regarding homozygous SNPs. There were 6188, 6340 and 8810 homozygous SNPs for Grofit, Magen and KH respectively. Statistics of shared and unique SNPs of each isolate are shown in Table 4, indicating that 2384 SNPs from the Grofit, 2408 SNPs from KH and 3304 SNPs from Magen are located in genes or up/downstream of genes (1000 bp). Additional annotation data of SNPs from each isolate are presented in Table 6. Similarly, 8 SNPs were detected in 5 out of 6 of the differentially expressed stress-related genes up/downstream of genes (1000 bp) ( Table 4).

Analysis of variation between isolates compared to the reference genome.
In order to compare the genomic variation between the studied isolates, homozygous SNPs positions of each isolate compared to the reference genome, were analyzed for common and unique SNPs. Results revealed that Magen isolate had the highest number of unique SNPs position (2880 SNPs, 30.7% of total positions). Followed by KH with 352 unique SNPs and Grofit with the lowest number of unique SNPs (46 SNPs, 0.5% of total positions) (Fig. 7A). Further analysis of common and unique genes containing SNPs revealed similar trend, as Magen had the highest number of unique SNPs in genes (819 genes), followed by KH with 23 unique SNPs in genes and Grofit with no unique SNPs located in genes (Fig. 7B).

Discussion
The results stem from both transcriptomic data for differential expression between isolates with high vs. low stress tolerance, and analysis of genomic data for SNPs among specific isolates with different degrees of stress tolerance. The obtained information contributes to our understanding of the differences between isolates with varied degrees of abiotic stress tolerance and to a future establishment of molecular markers to breed strains with improved performance.

Isolation and characterization of natural EPN populations.
To enrich the genetic pool, natural EPN populations were collected from different habitats across Israel. Abiotic factors are known to have strong effects on EPN occurrence, persistence and species variation [32][33][34] . However, previous studies conducted worldwide have shown varied correlations between effects of habitat parameters and EPN occurrence [35][36][37] . Initially, the recovery frequency and abundance of EPNs were evaluated according to the number of positive samples out of the total number of samples and collection sites. The results (5.8% recovery frequency and 23.5% abundance) were similar to recovery rates reported in previous surveys conducted in Israel and in other Mediterranean areas 30,31,[38][39][40] . The sampling period (January to May) was considered part of the rainy season in Israel, which extends from October to early May. However, rainfall peaks are common in December through February and rainfall varies considerably by region from north to south. To maximize the potential for EPN isolation, sampling was performed close to rainfall events and under tree canopies. To investigate the effect of habitat characteristics on the occurrence of EPNs, the collection sites were deliberately selected in different climatic regions, and in each region, samples were collected from cultivated and non-cultivated sites. Soil profiling of sand, silt and clay content revealed two main types of soil related to the climatic region: high clay soils were found only in Mediterranean regions, whereas high sand soils were found in arid and Mediterranean regions. Therefore, to study the relationship with EPN occurrence, the parameters of climate and soil type were integrated, producing four habitat types. Additional www.nature.com/scientificreports www.nature.com/scientificreports/ integration with the soil function parameter resulted in eight habitat types according to three parameters: climate, soil type and soil function. The occurrence of natural EPN populations was evaluated in relation to these habitat types, but none of the parameters were found to be significant. Nevertheless, certain trends corresponding with previous EPN surveys were recorded. Overall, the number of positive samples was higher in soils from cultivated sites compared to non-cultivated sites, in agreement with surveys reporting high recovery rates of EPNs from orchards and groves in the United States, Portugal and Israel [36][37][38]40 . Moreover, species variation was found among cultivated and non-cultivated sites. Isolates belonging to the species S. feltiae were only recovered from non-cultivated sites, whereas Heterorhabditis species were recovered only from cultivated soils. These findings might be explained by the foraging strategy of Heterorhabditis species, which are known to disperse deeper in the soil compared to 'ambusher' species of the genus Steinernema; they are therefore more influenced by tillage and agronomic practices 18,41 . Another interesting trend was that Steinernema isolates were only recovered from Mediterranean habitats, whereas Heterorhabditis isolates were recovered from arid as well as Mediterranean habitats. Nematodes of the genus Heterorhabditis are known to exist in habitats with a warm and arid climate, compared to Steinernema which are more abundant in humid and high-altitude habitats, usually characterized by lower temperatures 18,37,42 . Additional habitat parameters, such as soil temperature and altitude, showed minor differences between collection sites, and were therefore found irrelevant for the analysis of EPN occurrence.
In the current study we provide a taxonomic classification of naturally occurring Israeli EPNs to the species level, using molecular tools based on the ribosomal ITS region. Six Heterorhabditis isolates and two Steinernema isolates were identified. Another four isolates from previous studies were molecularly identified. Phylogenetic analyses revealed that the isolates belong to three different species: H. bacteriophora, H. indica and S. feltiae. This is the first confirmed report of H. indica in Israel, a species that was first reported in India and is widely distributed across all continents except for Antarctica. It has been reported as a dominant species in citrus groves in Florida, as well as in coastal regions in other Caribbean locations 38 32,34,37,38,43 . This might be due to aggregated spatial distributions of the species, which is a well-known distribution pattern for EPNs 18,33,35 . The patchy distribution might also provide a reasonable explanation for the lack of association of habitat parameters with EPN occurrence 30 . For ecological research purposes, this must be taken into account by randomly collecting a larger number of samples from each site.    www.nature.com/scientificreports www.nature.com/scientificreports/

Characterizing heat and desiccation tolerance of natural isolates of Heterorhabditis compared to a commercial line. Application of EPNs for biocontrol is restricted to habitats with favorable microenvi-
ronments, such as soil and cryptic habitats that provide protection from extreme environmental conditions 21,44 , unless applied with protective formulations. To increase the use of EPNs as biocontrol agents, the main objective is to overcome their intolerance to heat and desiccation stress. In the present study, the survival properties of natural isolates under extreme heat and desiccation conditions was evaluated in comparison to the commercial line EN-01. As the tolerance properties of this latter strain have been studied, it was chosen as a reference strain for the characterization of the tolerance capabilities of native populations from Israel. For the present study, we determined the MT 50°C and MW 50 of EN-01. The MT 50 was 38.8 °C, and the MW 50 was 0.97 aw , indicating high heat tolerance and low to medium desiccation tolerance according to previous studies 20, 45 . Phenotyping of the natural H. bacteriophora isolates under heat stress revealed three isolates with significantly lower tolerance and one isolate with similar tolerance compared to the commercial line EN-01. Under desiccation stress, none of the natural isolates displayed significant differences in tolerance compared to EN-01.
As the main objective of the research was to establish genetic selection tools based on the phenotypic characterization of natural EPN populations, we needed to obtain a wide range of phenotypic characteristics. We hypothesized that natural isolates will differ in their tolerance to heat and desiccation stresses. This hypothesis was confirmed by the variations in survival rates among species and among isolates within species. In general, H. indica isolates displayed higher survival rates than H. bacteriophora isolates. Such variations in tolerance among isolates and species have been reported previously; moreover, H. indica has been reported to be a high heat-and desiccation-tolerant 21,27,46,47 . It should be noted that the phenotypic characterization of stress tolerance in the present study was determined according to a single parameter-survival rate under stress-determined by stereoscopic observation. Further evaluation of morphological characteristics, virulence and reproduction of the surviving nematodes may provide valuable information on their tolerance capabilities.

Gene-expression patterns under heat stress in H. bacteriophora isolates.
RNASeq was only carried out for H. bacteriophora heat-exposed and control samples. This enabled comparing the transcriptomic and genomic data to the existing data on the commercial line EN-01. Isolates with significantly different survival rates under heat stress were chosen for RNASeq analysis in search of variation corresponding to the phenotypic differences. We hypothesized that gene expression would differ significantly between isolates with varied degrees of tolerance to environmental stress. The genes that varied the most were selected as candidate markers for the stress-resistance trait. As already noted, the commercial line EN-01 is characterized by high heat tolerance. This is the first study to provide transcriptomic data for this strain. The Israeli isolate, Grofit, was isolated in the Arava region and was phenotyped as having low heat tolerance in the present study. A comparison of gene expression in these two strains under heat stress of 35 °C vs. control conditions of 25 °C revealed a higher number of upregulated transcripts than downregulated transcripts in both strains. However, a comparison of the numbers of expressed transcripts revealed different expression levels between the strains. The high heat-tolerant strain EN-01 expressed both a lower number of upregulated transcripts and a higher number of downregulated transcripts, compared to the low heat-tolerant isolate Grofit. Previous studies have shown a similar relationship between gene-expression patterns and stress-tolerance capabilities, i.e., downregulation of transcripts in strains with high www.nature.com/scientificreports www.nature.com/scientificreports/ tolerance and upregulation of transcripts in those with low tolerance 28,29 . Although the overall expression pattern of both strains in the present study was upregulation, the results correspond with those previous findings.
Specific stress-related genes. The analysis of gene ontology and pathways of the differentially expressed transcripts revealed common expression of certain heat-shock proteins (HSPs) and related pathways of HSP activity, regardless of the strains' survival rates, supported by reports of HSP synthesis as a common response to stress 28 . Transcripts expressed only in the high heat-tolerant strain EN-01 were annotated as known stress-related genes, encoding glycerol kinase (GK), fatty acid desaturase (FAD) and zinc finger protein (ZFP). These three genes have been previously reported with similar expression patterns under heat stress, suggesting high suitability as candidate expression markers for heat tolerance 14,29,48,49 . Functional annotation of the uniquely downregulated transcripts of the low heat-tolerant isolate Grofit revealed the known stress-related gene encoding trehalose (TRE). High temperature is known to trigger the accumulation of trehalose in EPNs and the accumulated trehalose correlates with enhanced heat and desiccation tolerance 14 . These findings suggest downregulation of the trehalose-encoding gene as a candidate expression marker for the identification of low heat-tolerant strains.
Genomic variation. Identification of genomic variation between high and low stress-tolerant isolates might be used as a foundation for selection markers based on polymorphism in the genome. Natural EPNs populations have highly diverse phenotypes and genotypes 20,47 . In the present study, the genomes of three natural isolates of H. bacteriophora were sequenced to explore the genomic differences between them. We hypothesized that the genomic sequence would differ significantly between isolates with high tolerance vs. low tolerance to environmental stress. However, the results only partially supported this hypothesis. Two of the isolates, Magen and Grofit, were characterized as having low tolerance to heat and desiccation. The third isolate, KH, was characterized as having high tolerance to heat and desiccation, and the reference strain had high heat tolerance and low to medium desiccation tolerance. Magen was found to be the most distant isolate from the reference, whereas Grofit was most similar to the reference line EN-01. As Magen and Grofit have the same low-stress-tolerance phenotype, this result did not support the hypothesis. However, the distance of Magen (low tolerance) from KH (high tolerance), as well its distance from the reference strain, indicated distinct genomic differences between isolates with varied degrees of stress tolerance. Since the H. bacteriophora genome is still in its assembly version represented as scaffolds, we could not create a circus map to compare SNPs positions among isolates 50 . But, this comparison enabled the identification of common SNPs in the natural isolates which can be further studied for their association with phenotypic characterizations of EPNs. We identified 8 SNPs in the candidate genes for heat stress tolerance. Part of these SNPs were polymorphic between homozygous isolates. It is important to note that genotype-phenotype associations in natural populations using the single phenotypic parameter of survival is not ideal. The main objective of this work was to establish a basis for further research by screening natural populations and identifying various positions in their genomes that are potentially suitable for use as molecular markers for beneficial traits.
Further research and use of the obtained knowledge. The present study provides a phenotypic characterization of the survival rates of natural populations of EPNs belonging to two species: H. bacteriophora and H. indica. The phenotypic variations led to variations in the nematodes' survival rates. To identify the genotypic variations, the isolates that varied the most in their tolerance to each of the studied stress conditions were chosen for transcriptomic comparison. The results presented here are for two H. bacteriophora strains under heat stress. Several differentially expressed genes between the highly heat-tolerant strain EN-01 and the low heat-tolerant isolate Grofit were identified as stress-related genes reported in previous studies. These genes were chosen as candidate marker genes. Previous studies of EPNs under heat and desiccation stress have revealed that IJs use a conserved stress-tolerance response that may be triggered by changes in temperature or desiccation regime. Once this stress response is induced, the nematodes acquire resistance to multiple factors, including heat, desiccation and ultraviolet radiation 51 . Therefore, further study of the transcriptome of H. bacteriophora isolates under desiccation stress and H. indica isolates under heat and desiccation stress will provide additional important information regarding the expression patterns of stress-related genes under stress conditions and their suitability for use as expression markers for beneficial traits of multiple-stress resistance in EPNs. Eventually, the candidate genes will be further validated for use as expression markers using RT-qPCR under stress conditions. By studying the genomic sequence variations between natural isolates with varied degrees of stress tolerance, specific SNPs were found that may be associated with beneficial traits. Future study of the specific physical and functional positions of those SNPs along the H. bacteriophora genome will provide important information on their relevance as markers.
conclusion Natural populations of EPNs were isolated from different habitats across Israel and characterized at the species level, phenotypic level (degree of tolerance to heat and desiccation stress) and genotypic level (variations in gene expression and genome sequence). The phylogenetic analysis revealed three EPN species, S. feltiae, H. bacteriophora and H. indica, the latter reported for the first time in Israel. The phenotypic characterization of the different isolates under stress conditions revealed different degrees of tolerance in terms of survival after a defined period of stress. Pursuant to the main objective of the study, strains with the most different survival properties were chosen for genotypic characterization. This characterization revealed differences in the expression of specific known stress-related genes that need to be further validated. Furthermore, genomic sequence variations between natural isolates with varied degrees of stress tolerance were identified. Among these variations, specific SNPs were found which may be associated with beneficial traits. This study is the first to provide a phylogenetic characterization and genomic information on natural EPN populations in Israel, in comparison to a well-studied commercial line.
www.nature.com/scientificreports www.nature.com/scientificreports/ The obtained knowledge lays the groundwork for future studies aimed at identifying molecular markers as genetic selection tools for enhancement of EPNs' ability to withstand environmental stress conditions.

Soil samples collection and characterization.
During the rainy season of January to May 2018, 136 soil samples were collected from 34 sites across Israel (Fig. 6). Each site was of 100 m 2 cultivated or non-cultivated land. At each site, four samples were collected into pots from a distance of 30-40 cm from the tree trunks, under the canopy, at 15-20 cm depth using hand shovel. Each sample contained 500 gram of soil. Three of the samples were subjected for EPN isolation and the fourth sample was subjected to water content analysis and soil profiling. Water content analysis was performed according to the gravimetric method 52 . Soil samples were air-dried for one week in room temperature. The differences between the weight of wet soil and the weight of dry soil (g) were calculated as the mass of water in the sample. The gravimetric water content was calculated as the mass of water per mass of dry soil, presented as percentages. Soil profiling of sand, silt and clay content was done at the Field Service Lab, Neve Ya'ar. Soil temperature was measured at each site using a thermometer probe inserted to the soil at 10 cm depth. Soil function (cultivated or not) and the type of vegetation were recorded. Climate characteristics were obtained from the Israel Meteorological Service (IMS) (http://www.ims.gov.il), according to the annual precipitation.
Isolation of EPNs from soil. Isolation of nematodes was done selectively for entomopathogenic nematodes using a baiting technique as described before by Galleria mellonella (L.) larvae obtained from a lab colony at Volcani center 31,53 . One trap with 5 larvae was inserted to each soil sample and samples were incubated at 25 °C for one week. At the end of incubation, mortality rates of larvae in the traps were recorded. Dead larvae were further incubated in white traps at 25 °C in order to collect emerging IJs from the insect cadavers, as described by 54 . Emerging IJs were collected in water and stored in a Flask tissue culture at 12-14 °C.
Nematode culture. Recovered IJs from the white traps were re-exposed to last-instar G. mellonella larvae in a 50 cm diameter petri dish lined with moist filter paper. New generation of emerging IJs were collected in Flask tissue culture and kept in 12-14 °C for further research. Nematodes rearing during the research was done similarly, by infection of last-instar G. mellonella larvae with freshly emerged IJs every two weeks as described by 54  Amplified products at size of 600-800 bp were sequenced by Sanger method. Sequences were compared with ITS sequences in the GenBank by means of BLAST in NCBI and WormBase Parasite database. The ITS sequences of the various isolates were compared using the MAFFT software (Version 7, http://mafft.cbrc.jp/alignment/ server/). Bootstrap values were calculated by 100 repetitions as described by 55 . Phylogenetic tree, describing the phylogenetic relations between the various nematodes isolates was designed using TOL software (https://itol. embl.de/).
Host invasion rate assay. The infectivity rate of the different EPN isolates was determined by invasion host assay as described by Glazer and Lewis (2000). Five last-instar G. mellonella larvae were placed in a 50-cm-diameter petri dish lined with moist filter paper. A thousand IJs of each isolate were added to each petri dish and the dishes were incubated at 25 °C for 72 hours to ensure sufficient time for nematodes invasion and adult hermaphrodite development inside the larvae cadavers. Subsequently, larvae were washed in tap water in order to remove any nematodes from the outer surface. Larvae cadavers were dissected by scalpel under a stereomicroscope and tap water were added to the plates in order to encourage movement of the invading nematodes, for better visualization. The invasion rate of each nematode isolate was determined according to the number of nematode counted in each plate and presented as percentage of invasion out of the initial number of nematodes applied in each plate (1000 IJs per 5 larvae). Three replicates were carried out for each nematode isolate (one petri dish = one replicate) and the experiment was repeated three times.
Heat tolerance bioassay -Gradient temperature. Gradient 47 . Control tubes kept at 25 °C during the heat stress phase. The experiment was performed in PCR tubes with ~200 IJs per tube. All treatments were performed with a single batch of infective juveniles and repeated in three biological replicates. The experiment was repeated three times. Nematodes survival rate was determined for each treatment at the end of the experiment by counting the number of live and dead nematodes under a stereomicroscope (See scheme Fig. 8).
www.nature.com/scientificreports www.nature.com/scientificreports/ Desiccation tolerance bioassay -Gradient active water. Polyethylene glycol 600 (PEG 600, Sigma-Aldrich) dilutions were used to create gradient desiccation conditions by the hygroscopic method. PEG is a clear, non-toxic and non-ionic, hygroscopic solution 46 . Desiccation stress was measured as water activity (a w value) in each PEG dilution. Water activity is defined as the relative proportion of unbound water molecules in a sample. The lower the water activity of a solution is, the higher the desiccation is. IJs of EN-01 were subjected to the different PEG dilutions in order to determine their survival under increasing desiccation conditions. Initially, 3000 IJs were exposed to adaptation phase in 15 ml tubes containing 10% v/v PEG (0.987 a w ). The tubes shaken at 155 rpm for 72 hours at 25 °C. After adaptation, tubes were centrifuge in 5000 rpm for 5 minutes in order to precipitate the nematodes. The supernatant was discarded and IJs were re-suspend in 2 ml ringer buffer (9 gl-1 NaCl, 4.42 gl-1 KCl, 0.37 gl-1 CaCl2 ×2 H2O, 0.2 gl-1 NaHCO3) in order to wash the PEG solution. Wash step was performed three times. In the third time the supernatant was discarded, control IJs were re-suspended in 2 ml water, and the treated ones were suspended in 2 ml of 20% v/v (0.98 a w ), 30% v/v (0.975 a w ), 33% v/v (0.971 a w ), 35% v/v (0.968 a w ), 37% v/v (0.965 a w ), 40% v/v (0.96 a w ) and 50% v/v (0.93 a w ) PEG solution for 16 hours at 25 °C, shaken at 155 rpm. Nematode survival percentage was determined for each treatment by counting the number of live and dead nematodes under a stereomicroscope at the end of the experiment (See scheme Fig. 8). All treatments were performed with a single batch of infective juveniles and repeated in three biological replicates. The experiment was repeated three times.
Comparative heat tolerance bioassay -constant temperature. Two hundreds IJs from each nematode isolate and the commercial line, EN-01, were subjected to heat tolerance bioassay in PCR tubes. PCR program was set as follows: adaptation phase of 3 hours at 35 °C, recovery phase of 1 hours at 25 °C, stress phase of 4 hours at 38.8 °C (the determined MT °C 50 of EN-01), recovery phase at 25 °C overnight. Control tubes were kept at 25 °C during the heat stress phase. Each isolate had three biological replicates and the experiment was repeated three times. Nematodes survival percentage was determined for each isolate at the end of the experiment by counting the number of live and dead nematodes under a stereomicroscope.
Comparative desiccation tolerance bioassay -constant active water. Three thousands IJs from each nematode isolate and the commercial line, EN-01, were subjected to 15 ml tubes in a 10% v/v PEG solution for 72 hours in order to apply adaptation for desiccation stress. After adaptation, centrifugation and wash were carried out similarly to the gradient desiccation experiment. In the third wash, the supernatant was discarded, the control IJs were re-suspended in 2 ml water, and the treatment ones were re-suspend in 2 ml of 33.8% v/v PEG (the determined MW 50 of EN-01). Tubes were kept for 16 hours at 25 °C, shaken at 155 rpm. Each isolate had three biological replicates and the experiment was repeated three times. Nematodes survival percentage was determined for each isolate at the end of the experiment.
Preparation of total RNA. Gene expression patterns were studied in isolates with contrasting phenotype of heat stress tolerance. The commercial line, EN-01 (e-nema, GMBH, Schwentinental, Germany) was used as a reference strain for studying the gene expression patterns under stress conditions. Samples were taken at early stages of exposure to stress (Fig. 8) as the main events of gene expression take place at early stages of detection of stress by the nematodes 29 . Heat stress samples were taken after exposure to 35 °C for 3 hours. Three biological replicates performed for each condition. The control consisted of non-heated nematodes from the same batch heated counterparts, kept in distilled water at 25 °C. After the exposure period, samples were immediately frozen in liquid nitrogen and dried in a lyophilizer. Freeze-dried samples were ground to a fine powder and total RNA extracted from the desiccated, heated and control samples using the Plant/Fungi Total RNA purification Kit www.nature.com/scientificreports www.nature.com/scientificreports/ (Norgen Biotek Corp.) according to the manufacturer's instructions. Analysis was done in the Core Genomics Facility University of Illinois at Chicago. From the total RNA, the mRNA was isolated using selective Poly (A) column. cDNA libraries were enriched by PCR using NuGen Universal Plus mRNA chemistry kit according to the manufacturer's instruction. Cluster generation and sequencing of the cDNA transcripts were performed on a single lane of NovaSeq SP flow cell. Transcriptome analysis. Homogenous inbreed line of EN-01, IL-3, was used as a reference genome in the present study (Unpublished genome). Reads were aligned to IL-3 genome using Tophat2 software (v2.1). Gene abundance estimation was performed using Cufflinks (v2.2) combined with gene annotations from previous study 56 .
Differential expression analysis was completed using the DESeq. 2 R package. Genes that varied from the control more than twofold, with an adjusted P-value of no more than 0.05, were considered differentially expressed. Venn diagrams were calculated using "Venny" tool 57 . The transcriptome was used for a search of the NCBI non-redundant (nr) protein database, employing the DIAMOND program 58 . The results were exported to Blast2GO version 4.0 59 for gene ontology (GO) assignments. The eggNOG-mapper v2 tool) http:// eggnog-mapper.embl.de/ (was used for functional annotation. The KAAS tool (Kegg Automatic Annotation Server; http://www.genome.jp/tools/kaas/) was used for KEGG ontology and KEGG pathway assignments. Gene ontology enrichment analyses was carried out using Blast2GO 59 software based on Fisher's Exact Test 60 . The ReviGO web server was used for visualization of the GO terms in a semantic similarity-based scatterplot [http:// revigo.irb.hr] 61 .
Genome sequencing. Isolates of the specie H. bacteriophora with varying degrees of tolerance to each of the stress conditions were chosen for genome sequencing in order to search for polymorphism in the genome in relation to the tolerance capabilities. Two-hundred μl (~50 ng/μl) of genomic DNA of each isolate were prepared and sent for sequencing. Macrogen Inc. (South Korea) performed the genome sequencing using TruSeq DNA PCR free for library construction and NovaSeq. 6000 for the sequencing.
SNPs calling on genomic sequences. Homogenous inbreed line of EN-01, IL-3, was used as a reference genome in the present study (Unpublished genome). Reads were aligned to IL-3 genome using Tophat2 software (v2.1) using default parameters. Gene abundance estimation was performed using Cufflinks suite (v2.2) by the Cuffquant combined with gene annotations from previous study 57, and then Cuffnorm was used for the normalization. Paired-end reads of the 3 samples were mapped to the reference genome (IL-3) using the BWA mem program with default parameters 62 . The resulting mapping file was processed using Picard tool (http://broadinstitute. github.io/picard/; version 1.95) for adding read group information, sorting, marking duplicates, and indexing. Then, the local re-alignment process for locally re-aligning reads such that the number of mismatching bases is minimized across all the reads was performed using the Realigner Target Creator of the Genome Analysis Toolkit version 3.4-0 [GATK; version http://www.broadinstitute.org/gatk/] 63 . Finally, the variant calling procedure was performed using Haplotype Caller of the GATK toolkit, for the detection of SNPs between the variants and the reference IL-3 genome. Only homozygous SNPs were further analyzed. An in-house Perl script was used to define SNPs in genes or upstream/downstream to genes (1000 bp) and the define stop codon or non-synonymous substitutions.
Statistical analysis. All the statistical analysis was done using JMP, Version 14. SAS Institute Inc., Cary, NC, 1989-2019. Results are presented as mean ± SE of replicate analysis and are either representative of or include at least three independent experiments. Means of replicates were subjected to statistical analysis and considered significant when P < 0.05.
To determine Mean Tolerated Temperature (MT°C 50 ) and Mean Water Activity (MW 50 ) The mean percentages of survival in each treatment of heat or desiccation submitted to a Probit test with inverse prediction probability of 0.5 in order to find the predicted temperature and predicted active water value in which the survival of the commercial line, EN-01, is 50%, referring to MT °C 50 and MW 50 respectively. For comparative analysis of the heat and desiccation tolerance, the square roots of the proportion data were arcsine transformed and then subjected to ANOVA test to determine the percentage survival by strain of nematode. A multiple comparison, Dunnett with EN-01 as control was done in order to compare the survival rates of the different isolates with the commercial line.