Population genomics supports clonal reproduction and multiple independent gains and losses of parasitic abilities in the most devastating nematode pest

Abstract The root‐knot nematodes are the most devastating worms to worldwide agriculture with Meloidogyne incognita being the most widely distributed and damaging species. This parasitic and ecological success seems surprising given its supposed obligatory clonal reproduction. Clonal reproduction has been suspected based on cytological observations but, so far, never confirmed by population genomics data. As a species, M. incognita is highly polyphagous with thousands of host plants. However, different M. incognita isolates present distinct and overlapping patterns of host compatibilities. Historically, four “host races” had been defined as a function of ranges of compatible and incompatible plants. In this study, we used population genomics to assess whether (a) reproduction is actually clonal in this species, (b) the host races follow an underlying phylogenetic signal or, rather represent multiple independent transitions, and (c) how genome variations associate with other important biological traits such as the affected crops and geographical distribution. We sequenced the genomes of 11 M. incognita isolates across Brazil that covered the four host races in replicates. By aligning the genomic reads of these isolates to the M. incognita reference genome assembly, we identified point variations. Analysis of linkage disequilibrium and 4‐gametes test showed no evidence for recombination, corroborating the clonal reproduction of M. incognita. The few point variations between the isolates showed no significant association with the host races, the geographical origin of the samples, or the crop on which they have been collected. Addition of isolates from other locations around the world confirmed this lack of underlying phylogenetic signal. This suggests multiple gains and losses of parasitic abilities and adaptations to different environments account for the broad host spectrum and wide geographical distribution of M. incognita and thus to its high economic impact. This surprising adaptability without sex poses both evolutionary and agro‐economic challenges.


| INTRODUC TI ON
Nematodes annually cause severe damages to the world agricultural production, and the root-knot nematodes (RKN, genus Meloidogyne) are the most economically harmful species in all temperate and tropical producing areas (Moens et al., 2009;Jones et al., 2013).
Curiously, the most polyphagous RKN species, able to parasitize the vast majority of flowering plants on Earth (Trudgill & Blok, 2001), are described as mitotic parthenogenetic based on cytogenetics comparisons with outcrossing relatives (Triantaphyllou, 1981(Triantaphyllou, , 1985. This would imply absence of meiosis and obligatory asexual reproduction. Among these mitotic parthenogenetic RKN, Meloidogyne incognita is the most widespread species and is present, at least, in all the countries where the lowest temperature exceeds 3°C. Greenhouses over the world also extend its geographical distribution (Sasser, Eisenback, Carter, & Triantaphyllou, 1983). Meloidogyne incognita is so widely distributed that it is not even included on the list of regulated pests (Singh, Hodda, & Ash, 2013). Due to its worldwide distribution and extremely large range of hosts, M. incognita has been deemed the most damaging species of crop pest worldwide (Trudgill & Blok, 2001).
However, it has become evident that the whole broad host range of M. incognita, and other major RKN species, is not present in all the individuals within the species but that different "isolates" have different and overlapping ranges of compatible hosts (Moens et al., 2009). Variations regarding host range within one given species gave rise to the concept of "host race" as soon as 1952 (Sasser, 1952). Although RKN species can be differentiated based on morphological descriptions (Eisenback & Hunt, 2009), isozyme phenotypes (Carneiro, Almeida, & Quénéhervé, 2000;Esbenshade & Triantaphyllou, 1985) and molecular analysis (Blok & Powers, 2009), this is not the case of host races within a species (Triantaphyllou, 1985). Consequently, the pattern of compatibility/ incompatibility of the nematode interaction with a set of different host plants was standardized into the North Carolina Differential Host Test (NCDHT, (Hartman & Sasser, 1985)) to differentiate races within Meloidogyne spp. In M. incognita, all the populations originally tested reproduced on tomato, watermelon and pepper and none infected peanut, but they differed in their response to tobacco and cotton defining four distinct host races (Hartman & Sasser, 1985) (Table S1). Whether some genetic characteristics are associated with M. incognita races remains unknown. Indeed, no correlation between phylogeny and host races was found using RAPD and ISSR markers (Baum, 1994;Cenis, 1993;Santos et al., 2012). A different attempt to differentiate host races was also proposed based on repeated sequence sets in the mitochondrial genome (Okimoto, Chamberlin, Macfarlane, & Wolstenholme, 1991). Although the pattern of repeats allowed differentiating one isolate of race 1, one of race 2 and one of race 4; the study encompassed only one isolate per race, and thus, the segregation could be due to differences between isolates unrelated to the host-race status itself.
Hence, no clear genetic determinant underlying the phenotypic diversity of M. incognita isolates in terms of host compatibility patterns has been identified so far (Castagnone-Sereno, 2006). This lack of phylogenetic signal underlying the host races is surprising because it would suggest multiple independent gains and losses of host compatibly patterns despite clonal reproduction. Theoretically, animal clones have poorer adaptability because the efficiency of selection is impaired, advantageous alleles from different individuals cannot be combined, and deleterious mutations are predicted to progressively accumulate in an irreversible ratchet-like mode (Glémin & Galtier, 2012;Hill & Robertson, 1966;Kondrashov, 1988;Muller, 1964).
For these reasons, the parasitic success of M. incognita has long been described as a surprising evolutionary paradox (Castagnone-Sereno & Danchin, 2014). However, this apparent paradox would hold true only if this species actually reproduces without sex and meiosis while presenting substantial adaptability. So far, no whole-genome level study conclusively supports these tenets.
A first version of the genome of M. incognita was initially published in 2008 (Abad et al., 2008) and resequenced at higher resolution in 2017, providing the most complete M. incognita reference genome available to date (Blanc-Mathieu et al., 2017). This study showed that the genome is triploid with high average divergence between the three genome copies most likely because of hybridization events. Due to the high divergence between the homoeologous genome copies, and the supposed lack of meiosis, it was assumed that the genome was effectively haploid. The genome structure itself showed synteny breakpoints between the homoeologous regions and some of them formed tandem repeats and palindromes. These same structures were also described in the genome of the bdelloid rotifer Adineta vaga and considered as incompatible with meiosis (Blanc-Mathieu et al., 2017;Flot et al., 2013). However, whether these structures represent a biological reality or genome assembly artefacts remains to be clarified. Indeed, both genomes have been assembled using the same techniques and no independent biological validation for these structures has been performed. Hence, so far M. incognita and thus to its high economic impact. This surprising adaptability without sex poses both evolutionary and agro-economic challenges.

K E Y W O R D S
agricultural pest, clonal evolution, host races, Meloidogyne incognita, parallel adaptation, population genomics no conclusive evidence at the genome level supports the absence of meiosis.
Furthermore, because the reference genome was obtained from the offspring of one single female (originally from Morelos, Mexico), no information about the genomic variability between different isolates was available. Recently, a comparative genomics analysis, including different strains of M. incognita, showed little variation at the protein-coding level between strains collected across different geographical locations (Szitenberg et al., 2017), confirming previous observations with RAPD and ISSR markers (Baum, 1994;Cenis, 1993;Santos et al., 2012). However, no attempt was made to associate these variations with biological traits such as the host-race or geographical distribution. Moreover, the whole variability between isolates, including at the noncoding level, was so far never studied.
In the present study, we used population genomics analyses to investigate (a) whether the supposed absence of meiosis is supported by the properties of genomewide single-nucleotide variant (SNV) markers between isolates, (b) the level of variation between isolates at the whole-genome level and (c) whether these variations follow a phylogenetic signal underlying life-history traits such as the host compatibility patterns, the geographical distribution or the current host crop plant.
To address these questions, we have sequenced the genomes of 11 isolates covering the four M. incognita host races in replicates from populations parasitizing six crops across different locations in Brazil ( Figure 1). We used isozyme profiles, SCAR markers and the NCDHT to characterize the biological materials and then proceeded with DNA extraction and high-coverage genome sequencing. We identified short-scale variations at the whole-genome level by comparing the M. incognita isolates to the reference genome (Morelos strain from Mexico). We conducted several SNV-based genetic analyses to test for evidences (or lack thereof) of recombination. Using two different approaches, we classified the M. incognita isolates according to their SNV patterns and investigated whether the classification was associated with the following economically important biological traits: host race, geographical localization and current host plant.
Our population genomics analysis allowed addressing key evolutionary questions such as the nature of asexual reproduction in this animal species. We also clarified the adaptive potential of this devastating plant pest in relation to its mode of reproduction. In particular, we determined whether there is a phylogenetic signal underlying variations in biological traits of agro-economic importance such as the patterns of host compatibility (host races). While association between phylogenetic signal and patterns of host compatibilities would tend to show stable inheritance from ancestral states, the nonassociation would support multiple gains and losses of parasitic abilities and substantial adaptability. Abbreviations: ct, cotton "Deltapine 61"; pr, pepper "Early California Wonder,"; pt, peanut "Florunner,"; tb, tobacco "NC95,"; tm, tomato "Rutgers,"; wm, watermelon "Charleston Gray,". Host range results for the North Carolina Differential Host Test (NCDHT), numbers represent gall index with 0 = no galls; 1 = 1-2; 2 = 3-10; 3 = 11-30; 4 = 31-100; and 5 = more than 100 galls. This resolution has important agricultural and economic applications since crop rotation and other control strategies should take into account the adaptive potential of this nematode pest.

| Purification and characterization of M. incognita isolates
The M. incognita isolates (Table 1) originate from populations collected from different crops and geographically distant origins in Brazil ( Figure 1). For each isolate, one single female and its associated egg mass were retrieved as explained in Carneiro and Almeida (2001). To determine the species (here M. incognita), we used esterase isozyme patterns on the female (Carneiro et al., 2000). Two months after inoculation, the root system was rinsed with tap water, and egg masses were stained with Phloxine B (Hartman & Sasser, 1985) to count the number of galls and eggs masses separated for each root system. We assigned a rating index number according to the scale: 0 = no galls or egg masses; 1 = 1-2 galls or egg masses; 2 = 3-10 galls or egg masses; 3 = 11-30 galls or egg masses; 4 = 31-100 galls or egg masses; and 5 > 100 galls or egg masses per root system (Table 1) We categorized M. incognita host races based on their ability to parasitize tobacco and cotton (Table 1). Classically, the index for Rutgers tomato (susceptible control) is higher than 4 (+) (Hartman & Sasser, 1985). The rest of the population was kept for multiplication on tomato plants to produce enough nematodes for sequencing (typically >1 million individuals pooled together).

| DNA preparation and SCAR test
For each characterized nematode isolate, we extracted and purified the genomic DNA from pooled eggs with the supplement protocol for nematodes of the QIAGEN Gentra ® Puregene ® Tissue Kit with the following modifications: incubation at 65°C in the cell lysis buffer for 30 min and incubation at 55°C with proteinase K for 4 hr. We verified DNA integrity on 0.8% agarose gel and the DNA quantifica-

| Sequencing library preparation
We assessed input gDNA quantity using Qubit and normalized the samples to 20 ng/μl as described in TruSeq ® DNA PCR-Free Library Prep Reference Guide (#FC-121-3001, Illumina) prior fragmentation to 350 bp with Covaris S2. We assessed the quality of fragments after size selection and size control of the final libraries using High Sensitivity DNA Labchip kit on an Agilent 2100 Bioanalyzer.

| Whole-genome sequencing
We quantified sample libraries with KAPA library quantification kit (#7960298001, Roche) twice with two independent dilutions at 1:10,000 and 1:20,000. We calculated the average concentration of undiluted libraries and normalized them to 2 nM each then pooled them for sequencing step.
We generated high-coverage genomic data for the 11 M. incognita isolates by 2 × 150 bp paired-end Illumina NextSeq 500 sequencing with High Output Flow Cell Cartridge V2 (#15065973, Illumina) and High Output Reagent Cartridge V2 300 cycles (#15057929, Illumina) on the UCA Genomix sequencing platform, in Sophia-Antipolis, France. We performed two runs to balance the read's representation among the isolates and obtain homogeneity of coverage for the different samples (Table S2).

| Variant detection
We trimmed and cleaned the reads from each library with cutadapt (Martin, 2011)  ). Hence, the genome was considered in this analysis as mostly haploid. However, the distribution of per-base coverage on the genome assembly presented a 2-peaks distribution with a second minor peak at ~twice the coverage of the main peak ( Figure S2). Genome regions of double coverage most likely represent cases where two of the three homoeologous loci have been collapsed during the assembly, probably due to lower divergence.
Such regions will systematically be responsible for "artefactual" 0/1 SNV (presenting variations within isolates) as the reads from the two homoeologous copies will map a single collapsed region in the reference genomes. To avoid confusion between SNV representing true variations between individuals within isolates from those being artefacts due to collapsed homoeologous regions, 0/1 SNV were discarded from the analysis and only 1/1 SNV fixed within isolates were considered.
We used SAMtools (Li et al., 2009) to filter alignments with MAPQ lower than 20, sort the alignment file by reference position and remove multimapped alignments.
We used the FreeBayes variant detection tool (Garrison & Marth, 2012) to call SNV and small-scale insertions/deletions, incorporating all the library alignment files simultaneously and produced a variant call file (VCF). We filtered the resulting VCF file with the vcffilter function of vcflib (Anon, 2018), retaining the positions that had more than 20 Phred-scaled probability (QUAL) and a coverage depth (DP)> 10. As a comparison with an outcrossing meiotic diploid nematode, we conducted the same analyses on the genome of Globodera pallida (Eves-van den Akker et al., 2016). We first phased the SNV to haplotypes using WhatsHap (Martin et al., 2016) because the genome assembly mainly consists of collapsed paternal and maternal haplotypes.

| Genetic tests for detection of recombination
We used custom-made scripts (cf. Data Accessibility section) to calculate the proportion of fixed markers passing the 4-gametes test and linkage disequilibrium (LD) r 2 values as a function of intermarker distance along the M. incognita and G. pallida genome scaffolds.

| Genetic diversity between isolates, clusters and efficacy of purifying selection
We used bpppopstats from the Bio++ libraries (Guéguen et al., 2013) to estimate the nucleotide variability at nonsynonymous and synonymous sites as well as efficacy of purifying selection (π N , π S and π N /π S ) using a multiple alignment of the coding regions. We calculated fixation index (F ST ) for the three clusters using vcftools (Danecek et al., 2011).

| Principal component analysis
We performed a principal component analysis (PCA) to classify the isolates according to their SNV patterns and mapped the race characteristics, geographical location or current host plants on this classification. We used the filtered VCF file as input in the statistical package SNPRelate (Zheng et al., 2012) to perform the PCA with default parameters.

| Test for association between biological traits and genetic clusters
We used Fisher's exact test in R to assess whether there was a significant association between the SNV-based clusters and the host races or the crop species from which the isolates were originally collected.
We also conducted an isolation-by-distance (IBD) analysis using the adegenet R package (Jombart & Ahmed, 2011) to check how well the genetic distances correlate with geographical distances between the sampling points of the isolates. Geographical distances were calculated from exact sampling locations, when available, or centre points if the region was known but not the exact sampling location. Sample R3-4 was excluded from this analysis since it was a mix of samples pooled together from different geographical locations. L27 was also excluded since the sampling location was unknown.

| Mitochondrial genome analysis
We subsampled genomic clean reads to 1% of the total library for each M. incognita isolate. Then, we assembled them independently using the plasmidSPAdes assembly pipeline (Antipov et al., 2016).
We extracted the mitochondrial contigs based on similarity to the M. incognita reference mitochondrial genome sequence (NCBI ID: NC_024097). In all cases, the mitochondrion was assembled in one single contig. We identified the two repeated regions (63 bp repeat and 102 bp repeat), described in Okimoto et al. (1991), and we calculated the number of each repeat present in these regions.

| The M. incognita genome is mostly haploid and shows few short-scale variations
We collected 11 M. incognita populations from six different states across Brazil and from six different crops (soybean, cotton, coffee, cucumber, tobacco, watermelon) ( Figure 1) We examined the distribution of base coverage of SNV fixed within isolates (1/1 fixed SNV) and SNV that presented variations within at least one isolate (0/1 SNV). We observed that the 0/1 SNV, which were variable within isolates, showed a peak of distribution at ~ twice the coverage of the peak for fixed 1/1 SNV in the 11 isolates ( Figure S1). This parallels the distribution of base coverage in the M. incognita reference genome scaffolds which shows a major peak at ~65X and a second minor peak at ~130X (twice the coverage; Figure S2). These genome regions at double coverage were considered as representing highly similar pairs of homoeologous genome copies that were collapsed during the assembly (Blanc-Mathieu et al., 2017). Although these regions are minority in the genome assembly, they seem to be responsible for many 0/1 SNV (presenting within isolate variations). The SNV in these minority regions of double coverage probably results from genome reads of two homoeologous regions mapped to a single collapsed region in the reference assembly. Hence, most of these 0/1 SNV might not represent variations between individuals within an isolate but between the few collapsed homoeologous regions. Because the reference genome is mostly assembled in haploid status ( Figure S2), and the nature of 0/1 SNV is ambiguous, we will utilize only 1/1 SNV fixed within isolates as markers for all downstream analyses. Although this precludes analyses of variations between individuals within isolates, this allows a comparison of variations between isolates based on >66,000 solid fixed markers.

| No evidence for meiotic recombination in M. incognita
Based on cytogenetics observation, M. incognita and other tropical root-knot nematodes have been described as mitotic parthenogenetic species (Triantaphyllou, 1981(Triantaphyllou, , 1985. However, this evolutionary important claim has never been confirmed by genome-wide analyses so far. Using the fixed SNV markers at the whole-genome scale, we conducted linkage disequilibrium (LD) analysis as well as 4-gametes test to search for evidence for recombination (or lack thereof). In an outcrossing species, physically close markers should be in high LD, with LD substantially decreasing as distance between the markers increases, because of recombination, and eventually reach absence of LD similarly to markers present on different chromosomes. In clonal species, however, in the absence of recombination, the LD between markers should remain high and not decrease with increasing distance between markers. By conducting an analysis of LD, we did not find any trend for a decrease of LD between markers as a function of their physical distance (Figure 2a).
In contrast, the LD values remained high regardless the distance and oscillated between 0.85 and 0.94. Hence, we did not observe the expected characteristic signatures of meiosis. An inversely con- F I G U R E 2 Linkage Disequilibrium and 4-gametes test of M. incognita (a) and G. rostochiensis (b) isolates. The r 2 correlation between markers, indicating linkage disequilibrium (LD), is given as a function of the physical distance between the SNV markers (red line). The proportion of pairs of two-state markers that pass the 4-gamete test is given as a function of the distance between the markers (blue line). Because the G. rostochiensis genome assembly mostly consists of merged paternal and maternal haplotypes, we had to phase the SNVs before conducting LD and 4-gametes tests. The results were totally contrasted between M. incognita and G. rostochiensis (Figure 2b). In G. rostochiensis, the LD and 4-gametes curves started at relatively lower (<0.7) and higher (>0.15) values, respectively. Furthermore, we observed a rapid exponential decrease of r 2 in the first kb for LD. At an intermarker distance of 3 kb, the r 2 value was <0.37. In parallel, we observed a concomitant rapid and exponential increase in the proportion of markers passing the 4gametes test, which was >0.38 at the same intermarker distance.
Hence, while G. rostochiensis appears to display all the expected characteristics of meiotic recombination, this was not the case for M. incognita. This validates at a whole-genome scale the lack of evidence for meiosis previously proposed based on cytological observations in M. incognita. From the SNV falling in coding regions, we constructed a multiple alignment and measured nucleotide diversity at synonymous (π s ) and nonsynonymous (π n ) sites for the 11 isolates as well as the π n /π s ratio as a measure of the efficiency of selection. Consistent with the overall low number of SNV, the π s value across the 11 isolates was low (1.29 × 10 −03 ). This is one order of magnitude lower than the values measured for two outcrosser nematodes from the Caenorhabditis genus (Romiguier et al., 2014), C. doughertyi (formerly

| There is no significant association between the short-scale variants and the biological traits
Using principal component analysis (PCA) on the whole set of fixed SNV, we showed that the eleven M. incognita isolates formed three distinct clusters, which we named A, B and C (Figure 4). Cluster A is represented by isolate R1-2 alone, which has the highest number of variants. Cluster B is constituted by R3-2 and R4-4. The rest of the isolates fall in a single dense cluster C of overall low variation. There was no significant association between the clusters and F I G U R E 4 PCA of the Brazilian M. incognita isolates groups them into three clusters (A, B, and C). The geographical origins are associated with coloured shapes: black circle: Paraná, orange diamond: Santa Catarina, green square: São Paulo, red triangle: Mato Grosso, blue star: pool. Host plant representative pictures are displayed next to the isolates: soybean pod (R1-2 and R3-2); cotton flower (R3-1, R3-4, R4-4, and R4-1); coffee grain (R2-6); cucumber vegetable (R1-3); tobacco leaves (R1-6 and R2-1); and watermelon fruit slice (R4-3) the host-race status (Fisher's exact test p-value = 1, Appendix S1, Table S4). This implies that isolates of the same host race are not genetically more similar to each other than isolates of different host races. There was also no significant association between the SNVbased clusters and the crop plant from which nematodes were collected (Fisher's exact test p-value = 0.69, Appendix S1, Table S4).
Interestingly, the four different host races are all represented in one single cluster (C). Within this cluster, the total number of variable positions was 29,597, meaning that the whole range of host-race diversity is present in a cluster that represents only 44% of the total existing genomic variation. We also conducted an isolation-by-distance (IBD) analysis, which showed no correlation between the genetic distance and the geographical distance ( Figure S3).
To assess the levels of separation versus past genetic exchanges between these clusters, we calculated fixation index values (F ST ).
Weighted F ST values between clusters were all >0.83, suggesting a lack of genetic connections between the clusters (Table S5).
Using the mean F ST values, in contrast, while we observed a mean F ST > 0.98 between clusters A and B, indicating a lack of genetic connection between R1-2 and cluster B, the F ST values were much lower between A and C (0.35) and between B and C (0.52). This would suggest isolates from clusters A and B both result from a past bottlenecked dispersal and propagation from some isolates in cluster C.

| Phylogenetic networks confirm the lack of association of SNV with biological traits and support clonal evolution
Using a phylogenetic network analysis based on SNV present in coding regions, we could confirm the same three clusters ( Figure 5). This further supports the absence of phylogenetic signal underlying the host races (patterns of host compatibilities). Interestingly, this network analysis based on fixed SNV yielded a bifurcating tree and not a network. This result further supports a lack of genetic exchanges between the isolates and clonal reproduction. To confirm this result, we conducted separate phylogenetic analyses for each of the 14 longest scaffolds with sufficient number of phylogenetically informative variable positions and the mitochondrial genome. All these analyses showed a clear separation between the same three clusters (A, B and C) with some minor polytomies within cluster C ( Figure S5 and Figure S6).
According to the two classification methods (PCA and phylogenetic network), isolate R1-2 seemed to be the most divergent from the rest of isolates, which is consistent with its higher total number of SNV and number of isolate-specific SNV. Then, a small cluster was composed of isolates R3-2 and R4-3 (equivalent to cluster B of the PCA). Finally, a cluster (equivalent to PCA cluster C) grouped the rest of the eight isolates and covered all the defined host races as well as 5 of the 6 different host plants. Consistent with the PCA and phylogenetic network analysis, we also did not observe significant association between the number of repeats in the two repeat regions in the mitochondrial genome (63R and 102R) and races, geographical origin or host plant of origin (Table 2).

| Addition of further geographical isolates does not increase the genomic diversity and confirms the lack of association between genetic distance and biological traits
To investigate more widely the diversity of M. incognita isolates in relation to their mode of reproduction and other biological traits, we included whole-genome sequencing data from additional geographical isolates (Szitenberg et al., 2017). These genome data included one isolate from Ivory Coast, one from Libya, one from Guadeloupe such as host races, nature of the host of origin and geographical distribution (Appendix S1, Figure S7).
F I G U R E 5 Phylogenetic network for M. incognita isolates based on SNV present in coding sequences. The phylogenetic network based only on changes in coding sequences produced a bifurcating tree and shows the same grouping than the PCA analysis, into 3 distinct clusters

| D ISCUSS I ON
Is the parasitic success of M. incognita an evolutionary paradox? This proposition would be true only if M. incognita is adaptive despite having a fully parthenogenetic reproduction. Our results support these two properties.
The lack of sexual reproduction in M. incognita was so far only assumed based upon cytogenetic observations (Triantaphyllou, 1981(Triantaphyllou, , 1985 but never further supported at a molecular level. Here, the different analyses we performed at the population genomics level converge in supporting the lack of recombination and genetic exchanges in M. incognita. The phylogenetic network analysis based on fixed SNVs returned a bifurcating tree that separated the different isolates and not a network. This suggests a lack of genetic exchange between the isolates. In sexual "recombining" species, the mitochondrial genome accumulates mutations much faster than the nuclear genome. This is also true in the model nematode C. elegans where the mitochondrial mutation rate is at least two orders of magnitude higher than the nuclear mutation rate (Denver et al., 2009;Denver, Morris, Lynch, & Thomas, 2004). The higher mitochondrial accumulation of mutations is supposed to be the combined result of extremely rare or total lack of recombination, the low effective population size and the effectively haploid inheritance in mitochondria (Neiman and Taylor, 2009). In M. incognita, as opposed to C. elegans, we found that the percentage of variable positions in the mitochondrial genome is only one order of magnitude higher than in the nuclear genome. This suggests the nuclear and mitochondrial genomes evolve at a comparable rate and reinforces the idea that the nuclear genome is mostly effectively haploid and nonrecombining. Theoretically, the efficacy of selection should be lower in nonrecombining species than recombining ones.
We showed that the ratio of diversity at nonsynonymous sites/ diversity at synonymous sites (π n /π s ) was indeed one order of magnitude higher in M. incognita than in two outcrossing Caenorhabditis species. Finally, under recombination, the proportion of markers passing the 4-gametes test should increase exponentially with physical distance while linkage disequilibrium should decrease exponentially.
This was not observed in M. incognita, whereas a rapid exponential decrease of linkage disequilibrium was recently observed and considered an evidence for recombination in the bdelloid rotifer Adineta vaga (Vakhrusheva et al., 2018). Collectively, these results strongly suggest absence (or extremely rare) recombination and support the mitotic parthenogenetic reproduction of M. incognita.
Despite its clonal reproduction, it was already evident that  Figure S8). This reconstruction showed that the two hypotheses concerning the host range status of the last common ancestor were If the identity of a population is unknown, the crop selected for use in a management scheme may cause dramatic increases in nematode populations (Hartman & Sasser, 1985). However, the adaptability of M. incognita casts serious doubts on the durability of such strategies and must be taken into account in rotation schemes.
Furthermore, the biological reality of host races themselves is challenged by the lack of underlying genetic signal. Actually, the initial host-race concept has never been universally accepted, in part because it covered only a small portion of the whole potential variation in parasitic ability (Moens et al., 2009). Although M. incognita was already known to parasitize hundreds of host plants, only six different host standards were used to characterize four races. New host races might be defined in the future when including additional hosts in the differential test. Furthermore, using the same six initial host plant species, two additional M. incognita races that did not fit into the previously published race scheme have already been described (Robertson et al., 2009). Although the terminology "races" of Meloidogyne spp. has been recommended not to be used since 2009 (Moens et al., 2009), several papers related to M. incognita diversity of host compatibility or selection of resistant cultivars are still using this term, including on coffee (Lima et al., 2015;Peres et al., 2017); cotton (Mota et al., 2013;da Silva et al., 2014); or soybean (Fourie, McDonald, & Waele, 2006). This reflects the practical importance to differentiate M. incognita populations according to their different ranges of host compatibilities.
However, because these variations in host ranges are not monophyletic and thus do not follow shared common genetic ancestry, we recommend abandoning the term "race". Terms like "pathotype" or "biotype" that only refer to a phenotype and do not imply an underlying phylogenetic signal should be preferred (Sturhan, 1985).
What level of intraspecific genome polymorphism is required to cover the different ranges of host compatibilities in M. incognita and their ability to survive in different environments, despite their clonal reproduction? In this study, we found that the cumulative fixed divergence across the eleven isolates from Brazil and the reference genome (sampled initially from Mexico) reached ~ 0.02% of the nucleotides. Addition of isolates from Africa, the French Antilles and the United States did not increase the maximal divergence. This relatively low divergence is rather surprising, considering the variability in terms of distinct compatible host spectra (host races). Host-specific SNV were found only for Race 2 and no functional consequence for these SNV could be found, as they did not fall in predicted coding or evident regulatory regions. Furthermore, the existence of race-specific SNV themselves is even questionable as addition of other isolates might disqualify the few Race 2-spe- showing a surprisingly low number of nucleotide variation (only ~400 SNV on a ~3 Gb genome representing a proportion of variable positions of 1.3 × 10 −7 only). This also led the authors to suggest that mechanisms other than acquisition of point mutations and selection of the fittest haplotype must be involved (Gutekunst et al., 2018).
Previously, we have shown that the genome structure of M. incognita itself could participate in its versatility. Indeed, being allopolyploid, M. incognita has >90% of its genes in multiple copies. The majority of these gene copies showed diverged expression patterns one from the other and signs of positive selection between the gene copies have been identified (Blanc-Mathieu et al., 2017). How the expression patterns of these gene copies vary across different geographical isolates with different host compatibilities would be interesting to explore in the future. We thank the SPIBOC Bioinformatics platform for data hosting. The preprint version of this paper has been peer-reviewed and recommended by Peer Community in Evolutionary Biology (https ://doi. org/10.24072/ pci.evolb iol.100077)

CO N FLI C T O F I NTE R E S T
The authors declare that they have no financial conflict of interest with the content of this article.