A Proposed New Species Complex within the Cosmopolitan Ring Nematode Criconema annuliferum (de Man, 1921) Micoletzky, 1925

Ring nematodes are obligate ectoparasites on cultivated and wild herbaceous and woody plants, inhabiting many types of soil, but particularly sandy soils. This study explored the morphometrical and molecular diversity of ring nematodes resembling Criconema annuliferum in 222 soil samples from fruit crops in Spain, including almond, apricot, peach and plum, as well as populations from cultivated and wild olives, and common yew. Ring nematodes of the genus Criconema were detected in 12 samples from under Prunus spp. (5.5%), showing a low to moderate nematode soil densities in several localities from southeastern and northeastern Spain. The soil population densities of Criconema associated with Prunus spp. ranged from 1 nematode/500 cm3 of soil in apricot at Sástago (Zaragoza province) to 7950 and 42,491 nematodes/500 cm3 of soil in peach at Ricla and Calasparra (Murcia province), respectively. The integrative taxonomical analyses reveal the presence of two cryptic species identified using females, males (when available), and juveniles with detailed morphology, morphometry, and molecular markers (D2-D3, ITS, 18S, and COI), described herein as Criconema paraannuliferum sp. nov. and Criconema plesioannuliferum sp. nov. All molecular markers from each species were obtained from the same individuals, and these individuals were also used for morphological and morphometric analyses. Criconema paraannuliferum sp. nov. was found in a high soil density in two peach fields (7950 and 42,491 nematodes/500 cm3 of soil) showing the possibility of being pathogenic in some circumstances.


Introduction
Plant-parasitic nematodes of the family Criconematidae Taylor, 1936 [1] received the common name of 'ring nematodes' because of the body cuticle shows wide and prominent annuli. Ring nematodes are extensively, albeit not uniformly, distributed throughout the world [2]. They are obligate ectoparasites on cultivated and wild herbaceous and woody plants, inhabiting many types of soil but particularly sandy soils [3]. Several species of ring nematodes have been reported as parasitic and important pest of crops, causing damage to roots [4,5], but experimental evidence about the potential damage for crops induced by some species of criconematids is still lacking [6]. The genus Criconema Hofmänner & ("-") Not obtained or not performed.
Plants 2022, 11, 1977 4 of 36 Jaén and Cádiz provinces (southern Spain). These samples showed population densities from 12 to 37 specimens/500 cm 3 of soil. Identification of specimens of Criconema from soil below wild and cultivated olives resulted in species identical to those found in soil below Prunus spp., whereas two species morphologically closely related were detected in the soil sample below common yew, (about 50% each in soil numbers), but clearly separated by molecular and morphometric analyses, and described herein as Criconema paraannuliferum sp. nov. and Criconema plesioannuliferum sp. nov. (Table 1).

Species Delimitation Using Morphometry by Principal Component Analysis
In the maximum likelihood factor analysis (FA), the first three components (sum of squares (SS) loadings > 1) accounted for 65% of the total variance in the morphometric characteristics of the C. annuliferum-complex ( Table 2). The eigenvalues for each character were used to interpret the biological meaning of the factors. First, the maximum likelihood component 1 (MLC1) was dominated by number of annuli between the posterior end of the body and the vulva (RV), number of annuli on the tail (Ran) and for the c' ratio with high positive correlations (eigenvalues = 0.86, 0.90 and 0.73, respectively). This component was, therefore, related with the shape and size of the posterior part of nematode and tail. The ML2 was dominated by a high positive correlation for the nematode body length (L) (eigenvalue = 0.89), thereby relating this component to the overall nematode size. Finally, the ML3 was dominated by a high positive correlation for the c' ratio (eigenvalue = 0.67). This component was then related with tail shape. Overall, these results suggest that all the extracted components were related to the overall size and shape of nematodes in each population, especially the size and shape of the posterior part of the nematode. The results of the factor analyses were represented graphically in Cartesian plots where specimens from populations of the C. annuliferum-complex were projected on the plane of the x-and y-axes, respectively, as pairwise combinations of components 1 to 3 ( Figure 2). The specimens of the C. annuliferum-complex of all species were projected in the graphic representation showing an expanded distribution along the plane for all the projected combinations of the components owing to their wide morphometric variation within species and/or populations. This was more pronounced for C. paraannuliferum sp. nov. and C. annuliferum, where a high number of populations were considered [16,19,27,39,40]. Consequently, we did not detect a clear separation among C. paraannuliferum sp. nov. and C. annuliferum, with all the specimens belonging to these species being projected at random for all the projected combinations. However, and except for the projection on the plane of MLC2 and MLC3, where specimens of all species were randomly plotted, we observed that most specimens belonging to C. plesioannuliferum sp. nov. were spatially separated from the rest of species within the C. annuliferum-complex ( Figure 2). This spatial distribution was dominated by MLC1 accounting for 27% of the total of variance (42% of the proportion explained). This spatial separation was mainly dominated by MLC1 grouped species according to the number of annuli between posterior end of body and vulva (RV), number of annuli on Plants 2022, 11,1977 5 of 36 tail (Ran) and for the c' ratio (Table 2, Figure 2). Thus, specimens of C. plesioannuliferum sp. nov. having a higher number of annuli between the posterior end of the body and the vulva (RV), higher number of annuli on the tail (Ran) and higher values for the c' ratio were located at the right side, and on the opposite side were C. paraannuliferum sp. nov., while those of C. annuliferum were characterized by lower values for these diagnostic characters ( Figure 2) [16,27,39,40]. However, specimens of C. paraannuliferum sp. nov. and C. annuliferum having similar values for RV, Ran and the c' ratio were grouped among them, not showing spatial separation between the species (Figure 2). A minimum spanning tree (MST) superimposed on the plot of the first three principal components showed the same patterns observed with factor analysis. That is, a separation of C. plesioannuliferum sp. nov. but no clear separation between C. paraannuliferum sp. nov. and C. annuliferum within the C. annuliferum-complex ( Figure 2).

Species Delimitation Based on Ribosomal and Mitochondrial DNA
Species delimitation by ribosomal and mitochondrial DNA was based on four statistical parameters: intra-/inter-species variation, P ID (Liberal) values, posterior probability of clades on Bayesian analyses, and Rosenberg's P AB (Table 3). Analyses of species delimitation demonstrated that C. annuliferum from Belgium was clearly separated from C. paraannuliferum sp. nov. and C. plesioannuliferum sp. nov. The intra-and inter-species molecular variation for D2-D3 region of all three species was higher than 0.10, except for C. plesioannuliferum that was 0.09 (Table 3), suggesting that the probability of species identification with this gene is low. However, the variation for the ITS and COI genes were clearly below 0.10 ( Table 3), suggesting that the probability of species identification with these loci is high [42]. Similarly, the P ID (Liberal) values for all three species and locus were ≥0.93 suggesting that species can be adequately delimited [43,44]. Furthermore, all clade supports for the three loci were strong (PP = 1.00), and the Rosenberg's P AB values also support the monophyly of each species separately [45]. Table 3. Parameters evaluating the Criconema annuliferum-complex delimitation based on two rRNA genes (D2-D3 region of the 28S rRNA, ITS) and one mtDNA barcoding locus, COI for three Criconema species of the complex.

Description
Female: Body ventrally arcuate after heat relaxation, body annuli thick 7.5-10.5 µm wide, rounded with smooth edges and without anastomoses. Lip region low, with two annuli, the first annulus wide, anteriorly directed, and the second annulus narrower and forming a collar-like appearance. The SEM photographs show an elongate oral aperture, surrounded by a rounded oral plate, without submedian lobes, and six distinct hexaradiate pseudolips. Stylet long, generally straight, but slightly curved in some specimens, According to the dichotomic key of species within Criconema by Geraert [6], the new species needs to be compared with species of Criconema sharing a group of characters including mid-body annuli smooth and tail annuli not ornamented, stylet 80-130 µm long, no distinct differentiation of the lateral field, anterior vulval lip overhanging or overlapping posterior lip, RV = 8-19, R = 53-76, last tail annuli in line with the general slope of the tail, and the first lip annulus wider than following lip annuli [6]. The new species is morphologically very similar to C. annuliferum, and resembles Criconema crotaloides (Cobb, 1924) Schuurmans-Stekhoven & Teunissen, 1938 [46,47] and Criconema iranicum Azimi & Pedram, 2020 [48]. The morphological and morphometrical data from C. paraannuliferum sp. nov. for the three populations collected in soils from below peach, one in common yew, and two in cultivated and wild olives during this study are within the ranges of the original description of C. annuliferum, as well as those of Peneva et al. [40] from oak forests in Russia, and those of Etongwe et al. [16] from Belgium, and these minor differences are within the range of intraspecific variation. In addition, the new species can be separated from C. annuliferum by spicule length (30.5-37.0 vs. 51-53 µm). However, molecular characterization of the Belgian populations by ribosomal and mitochondrial genes clearly separated C. annuliferum from C. paraannuliferum sp. nov. in this study. The new species can be differentiated from C. crotaloides by body length (396-657 vs. 517-820 µm), RV (7-10 annuli from terminus vs. [11][12][13][14], and Ran (3-6 annuli from terminus vs. [6][7][8][9] [19,49]. Finally, it can be separated from C. iranicum by RV (7-10 vs. 9-11), and longer stylet (85-113 vs. 76.5-84.0 µm) [48].

Type Habitat and Locality
Criconema paraannuliferum sp. nov. was found in the rhizosphere of peach (coordinates 38 • 12 21.3" N, 1 • 42 23.9" W); the municipal district of Calasparra, Murcia province, southeastern Spain. Additional localities and host plants are reported in Table 1.

Etymology
The species epithet, paraannuliferum, refers to Gr. prep. para, alongside of and resembling, N.L. masc. n. annuliferum, because of its close resemblance to Criconema annuliferum.

Description
Female: Body almost straight to slightly ventrally arcuate after heat relaxation, body annuli thick 8.0-12.0 µm wide, rounded with smooth edges and without anastomoses. Lip region low, with two annuli, the first annulus wide, anteriorly directed, and the second annulus narrower and forming a collar-like appearance. The SEM photographs show a rounded oral aperture, surrounded by a rounded oral plate, without submedian lobes, and six distinct hexaradiate pseudolips. The SEM photographs show a type of rudimentary cord joining the body annuli from the 5th annulus to the caudal region ( Figure 10). Stylet long, generally straight but slightly curved in some specimens, representing 16.
According to the dichotomic key of species within Criconema by Geraert [6], the new species needs to be compared with species of Criconema sharing a group of characters including mid-body annuli smooth and tail annuli not ornamented, stylet 80-130 µm long, not distinct differentiation of the lateral field, anterior vulval lip overhanging or overlapping posterior lip, RV = 8-19, R = 53-76, last tail annuli in line with the general slope of the tail, and the first lip annulus wider than following lip annuli [6]. The new species is morphologically very similar to C. annuliferum and C. paraannuliferum sp. nov., and resembles Criconema crotaloides (Cobb, 1924) Schuurmans-Stekhoven & Teunissen, 1938 [46,47] and Criconema iranicum Azimi & Pedram, 2020 [48]. The morphological and morphometrical data of the population from common yew is within the ranges of C. paraannuliferum sp. nov. and C. annuliferum, except for ratios related with tail [16,40]. The new species can be separated from C. paraannuliferum sp. nov. by tail shape (abruptly narrowing to a conical shape, with 8-10 annuli, finely pointed terminus, and last annulus not lobed vs. conical with terminal annulus usually bi-lobed), RV (9-12 vs. 7-10 annuli from terminus), Ran (8-10
The COI gene sequences (ON648861-ON648884) showed high similarity values of 97.1% (differing in 20 to 21 nucleotides and 0 indel) from Criconema sp. from Denmark TSH-2020 (MN710767-MN710769) suggesting that they could be the same species. Unfortunately, no morphological data were found for Criconema sp. Denmark TSH-2020 and only COI gene sequences were available in GenBank. Species delimitation parameters based on COI supported that both species cannot be separated, although additional morphological data are needed to confirm this status.
All molecular markers studied, except for 18S rRNA, clearly separated the new species from other Criconema species. Because two closely related species were detected in the sample from common yew, the few male specimens were confirmed molecularly (100%) to belong to C. plesioannuliferum sp. nov.

Etymology
The species epithet, plesioannuliferum, refers to a compound name from the Greek word plesios = near, and annuliferum, the closet species of the genus Criconema.

Distribution of the Criconema annuliferum-Complex
In an exhaustive review of the geographical distribution of the Criconema annuliferum -complex in cultivated and natural environments in Spain and all over the world, we detected that this species complex has a wide distribution across a wide variety of herbaceous and woody hosts ( Figure 12). The Criconema annuliferum-complex is widely distributed in several European countries, but also in Africa [24], Asia [2,25], New Zealand [26], and South America [27]. It should be noted that the highest diversity seems to be in Spain, with three species in cultivated and wild environments ( Figure 12). Although these data suggest that other species in this species complex may be found in other countries after accurate integrative taxonomical identifications are conducted on them.

Distribution of the Criconema annuliferum-Complex
In an exhaustive review of the geographical distribution of the Criconema annuliferum -complex in cultivated and natural environments in Spain and all over the world, we detected that this species complex has a wide distribution across a wide variety of herbaceous and woody hosts ( Figure 12). The Criconema annuliferum-complex is widely distributed in several European countries, but also in Africa [24], Asia [2,25], New Zealand [26], and South America [27]. It should be noted that the highest diversity seems to be in Spain, with three species in cultivated and wild environments ( Figure 12). Although these data suggest that other species in this species complex may be found in other countries after accurate integrative taxonomical identifications are conducted on them.

Phylogenetic Analyses of the Criconema annuliferum-Complex
The D2-D3 domains of the 28S rRNA gene alignment (702 bp long) included 54 sequences of nine Criconema species and three outgroup species (Paratylenchus bukowinensis (MN088372), Paratylenchus enigmaticus (MZ265080), and Paratylenchus parastraeleni (MZ265065). Seventy-eight new sequences were included in this analysis. The Bayesian 50% majority rule consensus tree inferred from the D2-D3 alignment is given in Figure 13. For this region, all species that belong to the species complex Criconema annuliferum clustered together in a well-supported (PP = 1.00) clade, which was subdivided into two subclades, one of them (PP = 0.99) formed by Criconema paraannuliferum sp. nov. (ON705053-ON705072) and Criconema plesioannuliferum sp. nov. (ON705073-ON705080) and the other one (PP = 1.00) by C. annuliferum (MN783697-MN783702) and C. demani (MH828126, MH828128, MN628432 and MW938521). In this analysis, we detected that an unidentified species of Criconema from Italy (AY780952) clustered together with C. plesioannuliferum sp. nov., being molecularly identical. Thus, although no morphometrical data on this population are available, it most probably is conspecific with the new species.

Phylogenetic Analyses of the Criconema annuliferum-Complex
The D2-D3 domains of the 28S rRNA gene alignment (702 bp long) included 54 sequences of nine Criconema species and three outgroup species (Paratylenchus bukowinensis (MN088372), Paratylenchus enigmaticus (MZ265080), and Paratylenchus parastraeleni (MZ265065). Seventyeight new sequences were included in this analysis. The Bayesian 50% majority rule consensus tree inferred from the D2-D3 alignment is given in Figure 13. For this region, all species that belong to the species complex Criconema annuliferum clustered together in a well-supported (PP = 1.00) clade, which was subdivided into two subclades, one of them (PP = 0.99) formed by Criconema paraannuliferum sp. nov. (ON705053-ON705072) and Criconema plesioannuliferum sp. nov. (ON705073-ON705080) and the other one (PP = 1.00) by C. annuliferum (MN783697-MN783702) and C. demani (MH828126, MH828128, MN628432 and MW938521). In this analysis, we detected that an unidentified species of Criconema from Italy (AY780952) clustered together with C. plesioannuliferum sp. nov., being molecularly identical. Thus, although no morphometrical data on this population are available, it most probably is conspecific with the new species. The ITS rRNA gene alignment (739 bp long) included 49 sequences of eight species of Criconema and two outgroup sequences of Paratylenchus baldaccii (MW798336, MZ265015). Thirty-six new sequences were included in this analysis. The Bayesian 50% majority rule consensus tree inferred from the ITS alignment is given in Figure 14. The tree showed a well-supported subclade (PP = 1.00) with Criconema paraannuliferum sp. nov. (ON705081-ON705102) and Criconema plesioannuliferum sp. nov. (ON705103-ON705116), but clearly separated from each other. Thirty-six new sequences were included in this analysis. The Bayesian 50% majority rule consensus tree inferred from the ITS alignment is given in Figure 14. The tree showed a well-supported subclade (PP = 1.00) with Criconema paraannuliferum sp. nov. (ON705081-ON705102) and Criconema plesioannuliferum sp. nov. (ON705103-ON705116), but clearly separated from each other.
The 18S rRNA gene alignment (1694 bp long) included 48 sequences of 14 species of Criconema and two outgroup species (Tylenchocriconema alleni (KJ636364), and Paratylenchus shenzhenensis (KF668498)). Twenty-eight new sequences were included in this analysis. The Bayesian 50% majority rule consensus tree inferred from the 18S rRNA sequence alignment is given in Figure 15. The tree showed two well-supported (PP = 1.00) major clades, one of them appears in the basal part of the tree and includes only C. permistum, while the other major clade included the remaining species of Criconema. For this region, Criconema plesioannuliferum sp. nov. clustered into a poorly supported subclade with the other morphologically related species, C. crotaloides (HM116022). Criconema paraannuliferum sp. nov. formed a unique clade. The 18S rRNA gene alignment (1694 bp long) included 48 sequences of 14 species of Criconema and two outgroup species (Tylenchocriconema alleni (KJ636364), and Paratylenchus shenzhenensis (KF668498)). Twenty-eight new sequences were included in this analysis. The Bayesian 50% majority rule consensus tree inferred from the 18S rRNA sequence alignment is given in Figure 15. The tree showed two well-supported (PP = 1.00) major clades, one of them appears in the basal part of the tree and includes only C. permistum, while the other major clade included the remaining species of Criconema. For this region, Criconema plesioannuliferum sp. nov. clustered into a poorly supported subclade with the other morphologically related species, C. crotaloides (HM116022). Criconema paraannuliferum sp. nov. formed a unique clade. The COI gene alignment (674 bp long) included 110 sequences of 20 species of Criconema and three outgroup species (Paratylenchus baldaccii (MZ262220), Paratylenchus hamatus (MW797016) and Paratylenchus indalus (MW797005)). Sixty new sequences were included in this analysis. The Bayesian 50% majority rule consensus tree inferred from the COI sequence alignment is given in Figure 15. The phylogenetic relationships of Criconema spp. nov. inferred from analysis of this region were not well defined. The two new species grouped together in a low supported clade with Criconema sp. Denmark TSH-2020, which appeared forming a high supported (PP = 1.00) subclade with C. plesioannuliferum. The other morphologically related species to this species, C. crotaloides and C. annuliferum, clustered into two clearly separate well-supported clades (PP = 1.00). An unidentified species of Criconema from Ireland (MN710781) clustered within C. annuliferum (Figure 16), suggesting that it may be conspecific, but as no morphological data are available, additional studies are needed for confirmation Plants 2022, 11, x FOR PEER REVIEW 28 of 38 Figure 15. Phylogenetic relationships within the genus Criconema. Bayesian 50% majority rule consensus tree as inferred from 18S rRNA sequence alignment under the transition model with invariable sites and a gamma-shaped distribution (TIM2 + I + G). Posterior probabilities of more than 0.70 are given for appropriate clades. Newly obtained sequences in this study are in bold. The scale bar indicates expected changes per site, and the coloured boxes indicate the clade association of the Criconema annuliferum-complex.
The COI gene alignment (674 bp long) included 110 sequences of 20 species of Criconema and three outgroup species (Paratylenchus baldaccii (MZ262220), Paratylenchus hamatus (MW797016) and Paratylenchus indalus (MW797005)). Sixty new sequences were included in this analysis. The Bayesian 50% majority rule consensus tree inferred from the COI sequence alignment is given in Figure 15. The phylogenetic relationships of Criconema spp. nov. inferred from analysis of this region were not well defined. The two new species grouped together in a low supported clade with Criconema sp. Denmark TSH-2020, which appeared forming a high supported (PP = 1.00) subclade with C. plesioannuliferum. The other morphologically related species to this species, C. crotaloides and C. annuliferum, clustered into two clearly separate well-supported clades (PP = 1.00). An unidentified species of Criconema from Ireland (MN710781) clustered within C. annuliferum (Figure 16), suggesting that it may be conspecific, but as no morphological data are available, additional studies are needed for confirmation

Discussion
The primary objective of this study was to identify and molecularly characterize species belonging to the Criconema annuliferum-complex associated with cultivated and wild plants (viz. Prunus spp., wild and cultivated olive, and common yew) in southern Spain using integrative taxonomical approaches (morphological, morphometrical and molecular). Our results demonstrated that morphological studies integrated with rRNA and mitochondrial DNA molecular markers revealed the cryptic diversity of Criconema annuliferum species complex, enabling the description of two new species, C. paraannuliferum sp. nov. and C. plesioannuliferum sp. nov. Ribosomal and mitochondrial markers (D2-D3 Figure 16. Phylogenetic relationships within the genus Criconema. Bayesian 50% majority rule consensus tree as inferred from cytochrome c oxidase subunit 1 (COI) sequence alignment under the three-parameter model with invariable sites and gamma distribution model (TPM3uf + I + G). Posterior probabilities of more than 0.70 are given for appropriate clades. Newly obtained sequences in this study are in bold. The scale bar indicates expected changes per site, and the coloured boxes indicate the clade association of the Criconema annuliferum-complex. *** Criconema sp. Ireland TSH2020 should be considered to be C. annuliferum.

Discussion
The primary objective of this study was to identify and molecularly characterize species belonging to the Criconema annuliferum-complex associated with cultivated and wild plants (viz. Prunus spp., wild and cultivated olive, and common yew) in southern Spain using integrative taxonomical approaches (morphological, morphometrical and molecular). Our results demonstrated that morphological studies integrated with rRNA and mitochondrial DNA molecular markers revealed the cryptic diversity of Criconema annuliferum species complex, enabling the description of two new species, C. paraannuliferum sp. nov. and C. plesioannuliferum sp. nov. Ribosomal and mitochondrial markers (D2-D3 expansion domains of the 28S rRNA gene, ITS rRNA gene, and the mtDNA gene COI) are important tools for accurate identification of Criconema spp. In the present study, the mtDNA gene COI was used with others to separate C. plesioannuliferum sp. nov. and C. paraannuliferum sp. nov. from the Valdepeñas population, as it was the region that showed the highest interspecific variability. The broad range of cultivated and wild hosts indicated that C. paraannuliferum sp. nov. is quite generalized and lives in a wide range of environments and host plants. Multivariate morphometric analyses have proven to be useful tools for species delimitation within the genera Longidorus and Xiphinema [50,51]. Cryptic speciation has commonly been reported in criconematids, consequently these data boosted the hypothesis that criconematid nematodes are a hyper diverse group of organisms [52][53][54]. The recognition of these two new species within the C. annuliferumcomplex may have a direct influence on the real geographic distribution of these species. The low to moderate nematode population density in the common yew sample, as well as the sequencing of numerous individuals from this site, confirmed that each of the two species represented about half of the sample, suggesting that plant resources are good enough for supporting both species in the same area, and no competitive exclusion is asserted [55]. Notwithstanding this data, additional studies need to be conducted, such as life history phenology to confirm this hypothesis. Interestingly, in some cases, the levels of some of this species (C. paraannuliferum sp. nov.) are extremely high (more than 7000 nematodes/500 cm 3 of soil) in peach and presumably could affect the plant growth. However, phytopathological tests are necessary to understand if these levels are significantly affecting peach production. Nevertheless, the wide distribution and higher molecular diversity of C. paraannuliferum sp. nov. in several cultivated and natural environments suggests indirect dispersion by agricultural practices (movement of soil or plant materials) from naturally infested to uninfested soils, although additional analyses on all previous records on C. annuliferum from Spain and all over the world can support this hypothesis.
Multivariate morphometric analyses have proven to be useful tools for species delimitation within soil nematodes, especially in plant-parasitic nematodes such as those belonging to the genus Xiphinema [56][57][58] and the genus Longidorus [51,58]. The present study is, to our knowledge, the first one to apply multivariate methods to delimit and decipher species boundaries within species complexes within the genus Criconema. Our data support that the C. annuliferum-complex comprises a model example of morphostatic speciation (that is, genetic modifications not reflected in morphology and morphometry) [59], as independent approaches based on molecular analyses using ribosomal and mitochondrial sequence data clearly separated the species considered within the C. annuliferum-complex. The results of the multivariate analysis identified the shape and size of the posterior part of the nematode described by the number of annuli between the posterior end of the body and the vulva (RV), the number of annuli on tail (Ran) as well the c' ratio as key morphometric characters to differentiate some closely related species within the C. annuliferum-complex (Table 2, Figure 2). These results agree with the taxonomic statement outlining the number of annuli describing the nematode body as a fundamental feature in identifying the species within the genus Criconema [2]. Although some specimens of C. plesioannuliferum sp. nov. share similar values for most morphological characters with the other species included in this study, multivariate analysis allowed us to differentiate this species within this cryptic complex using a discrete number of characters (Table 2, Figure 2). However, multivariate analysis also supports the idea that C. paraannuliferum sp. nov. and C. annuliferum could resemble the same species (Figure 2) based in the wide morphometric variation and the similar values they share for the most diagnostic characters that identify species in the genus Criconema. This point could difficult their accurate identification. Finally, our results may also suggest that the C. annuliferum-complex comprises an endemic lineage within Criconema that has diversified in the Iberian Peninsula. However, further studies exploring the occurrence of the new taxa in other areas are necessary to confirm this hypothesis.
The lack of intraspecific variability in ribosomal and mitochondrial markers in C. plesioannuliferum sp. nov., may suggest a continuous isolation of this population under this natural environment under maintained biological (plant-hosts) and ecological characteristics (soil, temperature, etc.). Alternatively, it may suggest a recent speciation event from C. paraannuliferum sp. nov., like that inferred for other criconematids, such as M. erinaceum Powers, Mullin, Higgins, Harris & Powers, 2016 [60].
Phylogenetic analyses based on the D2-D3 region, ITS, 18S rRNA, and the COI gene using BI, mostly agree with the clustering obtained by other authors [13,16,60,61]. Ribosomal and mitochondrial based phylogenies clearly separate the C. annuliferum-complex into three separate species, which was confirmed by morphometric and molecular species delimitation analyses. Unfortunately, a concatenated analysis of the three ribosomal genes was not undertaken due to some sequences not being available for all species.
The present results confirmed previous data describing the remarkable biodiversity of several groups of plant-parasitic nematodes in southern Spain, such as species within the family Longidoridae (including virus vector nematodes of the genera Xiphinema and Longidorus) or pin nematodes of the genus Paratylenchus [51,62,63], and warrant additional sampling efforts to clarify the real biodiversity in our country.

Sampling Sites and Nematode Morphological Identification
A survey was conducted in the five principal areas of stone-fruit production (Prunus spp.) in Spain (Almanzora, Ebro, Guadalquivir, Júcar and Segura Valleys), which were established by river valleys. A total of 219 sites were sampled in the present study, and twelve showed the presence of specimens of putative Criconema annuliferum. Soil samples for nematode analysis were collected with a shovel from below four to five randomly selected trees and mixed to constitute a soil sample from each sampling site; samples came from the upper 5-40 cm depth of soil. Nematodes were extracted from a 500 cm 3 sub-sample of soil by centrifugal flotation [64]. In addition, three samples from below common yew and cultivated and wild olives in two localities of Jaen (Valdepeñas and Castillo de Locubín), and one locality from Cádiz province (Prado del Rey), respectively, were selected and studied here, because they contained specimens of Criconema closely resembling those from below Prunus spp. Nematode identification was performed using an integrative approach, combining morphological and morphometrical evaluation with molecular techniques. Morphological and morphometrical analyses were conducted using fixed individuals mounted on permanent slides. To prepare the fixed material, specimens of Criconema were killed at 70-75 • C and fixed in an aqueous solution of 4% formaldehyde + 1% glycerol, dehydrated using alcohol-saturated chamber and processed to pure glycerin using Seinhorst's method [65] as modified by De Grisse [66]. A total of 94 individuals including 79 females, 8 males and 7 juveniles were used for morphological and morphometrical analyses. Fixed mounted individuals were then examined and measurements of each nematode population, including important diagnosis characteristics (i.e., de Man indices, body length, stylet length, R, Rst, Roes, Rex, RV, Rvan, Ran, and the ratio VL/VB) [41], were performed using a Leica DM6 compound microscope with a Leica DFC7000 T digital camera.
Females and fourth-stage juveniles of each species mounted in glycerin were selected for SEM observations. The nematodes were hydrated in distilled water, dehydrated in a graded ethanol-acetone series, critical point-dried, coated with gold, and observed with a Zeiss Merlin scanning electron microscope (5 kV) (Zeiss, Oberkochen, Germany) [67].
Voucher specimens of these described species were deposited in the nematode collection of Institute for Sustainable Agriculture, IAS-CSIC, Córdoba, Spain.

DNA Extraction, PCR and Sequencing
For molecular analyses, and to avoid mistakes in case of mixed populations in the same sample, single nematodes were previously mounted in a drop of NaCl and used for molecular identification after recording morphological data. Genomic DNA extraction from single specimens was conducted as described by Palomares-Rius et al. [68]. Briefly, an individual nematode is cut using a scalpel in a drop of PCR buffer (ThermoPol ® , Biolabs, New England, USA) (20 µL) and 2 µL proteinase K (600 µg/mL) were added. The tubes were frozen at −80 • C (15 min), then incubated at 65 • C (1 h) and at 95 • C (10 min) consecutively. Tubes were centrifuged (1 min, 16,000× g) and kept at −20 • C until use in PCR; and more importantly, all three molecular markers for each population of Criconema were extracted from the same single individual in each PCR tube without any exception. In addition, male conspecificity was confirmed by single DNA extraction of males for each species.
The D2 and D3 expansion domains of the 28S rRNA were amplified using the D2A (5 -ACAAGTACCGTGAGGGAAAGTTG-3 ) and D3B (5 -TCGGAAGGAACCAGCTACTA-3 ) primers [69]. The Internal Transcribed Spacer region (ITS) was amplified by using forward primer TW81 (5 -GTTTCCGTAGGTGAACCTGC-3 ) and reverse primer AB28 (5 -ATATGCTTAAGTTCAGCGGGT-3 ) [70]. The partial 18S rRNA was amplified using the primers 988 (5 -CTCAAAGATTAAGCCATGC-3 ) and 1912R (5 -TTTACGGTCAGAACTAG GG-3 ) [71]. The COI gene was amplified using the primers COI-F5 (5 -AATWTWGGTGTTG GAACTTCTTGAAC-3 and COI-R9 (5'-CTTAAAACATAATGRAAATGWGCWACWACAT AATAAGTATC-3 ) [72]. The PCR cycling conditions for the 28S rRNA, ITS and 18S rRNA were as follows: 95 • C for 15 min, followed by 35 cycles of 94 • C for 30 s, an annealing temperature of 55 • C for 45 s, and 72 • C for 1 min, and one final cycle of 72 • C for 10 min. The PCR cycling for COI primers was as follows: 95 • C for 15 min, 39 cycles at 94 • C for 30 s, 53 • C for 30 s, and 68 • C for 1 min, followed by a final extension at 72 • C for 7 min. The PCR volumes were adapted to 25 µL for each reaction, and primer concentrations were as described in De Ley et al. [69], Subbotin et al. [13], Holterman et al. [71] and Powers et al. [72]. We used 5x HOT FIREpol Blend Master Mix (Solis Biodyne, Tartu, Estonia) in all PCR reactions. The PCR products were purified using ExoSAP-IT (Affimetrix, USB products, Kandel, Germany) and used for direct sequencing in both directions with the corresponding primers. The resulting products were analysed in a DNA multicapillary sequencer (Model 3130XL Genetic Analyzer; Applied Biosystems, Foster City, CA, USA), using the BigDye Terminator Sequencing Kit v.3.1 (Applied Bio-systems) at the Stab Vida sequencing facility (Caparica, Portugal). The sequence chromatograms of the four markers (18S rRNA, ITS, COI and D2-D3 expansion segments of 28S rRNA) were analysed using DNASTAR LASERGENE SeqMan v. 7.1.0. The Basic local alignment search tool (BLAST) at the National Center for Biotechnology Information (NCBI) was used to confirm the species identity of the DNA sequences obtained in this study [73]. The newly obtained sequences were deposited in the GenBank database under accession numbers indicated on the phylogenetic trees and in Table 1.

Recognition of Putative Species within the Criconema annuliferum-Complex and Species Delimitation Approach
Two independent strategies of species delimitation using morphometric and molecular data were used to determine species boundaries within this species complex. The recognition of the group of species included in this complex was established from large-scale taxonomic studies in the genus Criconema [2]. Specifically, morphological comparison showed that several of the diagnostic characters defining the genus Criconema were characteristic of the whole species complex, highlighting the lip region shape composed by two annuli being the first wider than the second one as well as the total number of annuli on body (R) and between posterior end of body and vulva (RV) [2]. We named the group the C. annuliferum-complex, after the oldest described species within it. The main diagnostic features of C. annuliferum were used to determine the morphologically closest species to the complex. Thus, in addition to the new taxa, C. paraannuliferum sp. nov. and C. plesioannuliferum sp. nov., the selected species included in this species complex were C. crotaloides and C. annuliferum [2]. However, the recognition of the group of species used in both species delimitation approaches as belonging to the C. annuliferum-complex (i.e., similar key morphometric characters) was also determined by the presence of accurate studies using an integrative taxonomic strategy to confirm species identity (that is, availability of molecular data linked to accurate morphology identification by morphometry in order to avoid misidentifications). For this reason, the previously described species, C. crotaloides, was not included in both strategies because of the lack of available information for it in the literature [17].
Species delineation using morphometry was conducted with exploratory factor analysis (FA) to estimate the degree of association among species within the C. annuliferumcomplex [74]. Factor analysis was based upon the following morphological characters: L (body length), stylet length, R, Rst, Roes, Rex, RV, Rvan, Ran, and the ratios a, c, c', V, VL/VB (Table 2) [41]. Several populations from natural and agricultural areas were used for some of the new taxa included in the C. annuliferum-complex (Table 1). In the case of the described species C. annuliferum, all the nematode populations were selected based on the availability of molecular data, their accurate identification using morphological and morphometrical comparison with other reported populations and closely related species, and/or the species distribution. Thus, the nematode populations selected for C. annuliferum were as follows: two populations from Netherlands including the type locality [19,39], a population from Chile [27], two populations from Russia [40], and three populations from Belgium [16]. For these populations, we used the mean value of the morphological characters mentioned above. Overall, 85 female specimens were used in a multivariate approach for the C. annuliferum-complex. Prior to the statistical analysis, diagnostic characters were tested for collinearity [75]. We used the collinearity test based on the values of the variance inflation factor (VIF) method that iteratively excludes numeric covariates showing VIF values > 10 as suggested by Montgomery et al. [76]. The FA was performed using the Maximum Likelihood algorithm through a decomposition of the data matrix between populations using the fa function implemented in the software package 'psych' [77]. This analysis produced a set of variables (factors) that are linear combinations of the original variables, independent of each other and ranked according to the amount of variation accounted for. Orthogonal varimax raw rotation was used to estimate the factor loadings and only factors with sum of squares (SS) loadings > 1 were extracted. Finally, a minimum spanning tree (MST) based on the Euclidean distance was superimposed on the scatter plot of the C. annuliferum-specimens complex against the FA axes. The MST was performed using the ComputeMST function implemented in the software package 'emstreeR' [78]. All statistical analyses were performed using the R v. 3.5.1 freeware [79].
Species delineation based on molecular data was performed using the species delimitation plugin [42] from the program Geneious Prime v2022.1.1. (Geneious, Auckland, New Zealand), and was used to calculate intra-and inter-species variation by means of the P ID liberal and the Rosenberg's P AB value. The intra-and inter-species molecular variation was determined by calculating the ratio between the average genetic distance between individuals within a species and the average genetic distance between individuals belonging to the sister species (the average pairwise tree distance among members of a putative species/the average pairwise tree distance between the members of one putative species and the members of the closest second putative species), if the ratio is less than 0.10 the probability of species identification is high [42]. The P ID (Liberal) value [44] represents the probability that a correct species identification would be made using the best sequence alignment (BLAST), closest genetic distance or placement on a tree (falling within or being sister to a monophyletic species clade). Species with P ID (Liberal) ≥ 0.93 were considered to be adequately delimited [43]. The Rosenberg's P AB represents the probability that the monophyly of a group of sequences is the result of random branching [45].

Phylogenetic Analyses
The D2-D3 expansion segments of the 28S rRNA, ITS rRNA, 18S rRNA, and COI mtDNA sequences of the seven populations of Criconema were obtained in this study. These sequences and other sequences of Criconema spp. from GenBank were used for phylogenetic analyses. Selection of outgroup taxa for each dataset were based on previously published studies [16,61]. Multiple sequence alignments of the different genes were completed using the FFT-NS-2 algorithm of MAFFT V.7.450 [80]. The BioEdit program V. 7.2.5 [81] was used for sequence alignments visualization and edited by Gblocks ver. 0.91b [82] using options for a less stringent selection (minimum number of sequences for a conserved or a flanking position: 50% of the number of sequences + 1; maximum number of contiguous nonconserved positions: 8; minimum length of a block: 5; allowed gap positions: with half). The phylogenetic analyses of the sequence datasets were based on Bayesian inference (BI) using MrBayes 3.1.2 [83]. The best-fit model of DNA evolution was achieved using JModelTest V.2.1.7 [84] with the Akaike information criterion (AIC). The best-fit model, the base frequency, the proportion of invariable sites, and the gamma distribution shape parameters and substitution rates in the AIC were then used in MrBayes for the phylogenetic analyses. The general time-reversible model with a gamma-shaped distribution (GTR + G) for the D2-D3 segments of 28S rRNA, the transversion model with a gamma-shaped distribution (TVM + G) for ITS rRNA region, the transition model with invariable sites and a gammashaped distribution (TIM2 + I + G) for the partial 18S rRNA gene, and the three-parameter model with invariable sites and gamma distribution model (TPM3uf + I + G) for COI gene, were run with four chains for 4 × 10 6 generations, respectively. A combined analysis of the three ribosomal genes was not undertaken due to some sequences not being available for all species. The sampling for Markov chains was conducted at intervals of 100 generations. For each analysis, two runs were conducted. After discarding burn-in samples of 30% and evaluating convergence, the remaining samples were retained for more in-depth analyses. The topologies were used to generate a 50% majority-rule consensus tree. On each appropriate clade, posterior probabilities (PP) were given. FigTree software version v.1.4.3 [85] was used for visualizing trees from all analyses.

Conclusions
This study proves the importance of using integrative taxonomy for the identification of species of Criconema. It establishes the existence of cryptic biodiversity within the C. annuliferum-complex, increasing and expanding the diversity of this group of nematodes in Spain. For the first time, our results demonstrate the coexistence of two closely related species of Criconema within the same sample in a natural environment, suggesting that both species could coexist without competing between themselves (similar individual numbers in the soil sample), while only C. paraannuliferum sp. nov. was found in cultivated Prunus, probably distributed within cultivated areas by several indirect processes (movements of soils, plant material, etc.). This study provides ribosomal and mitochondrial markers for precise and unequivocal diagnosis of species in the C. annuliferum-complex and suggests that other reports of C. annuliferum in Spain and around the world need to be confirmed with molecular markers.