Measurement(s)genome • sequence_assembly • sequence feature annotationTechnology Type(s)DNA sequencing assay • sequence assembly process • sequence annotationSample Characteristic - OrganismMeloidogyne enterolobiiMachine-accessible metadata file describing the reported data: 10.6084/m9.figshare.12410363
Root-knot nematodes (genus Meloidogyne) are plant parasites causing huge economic loss in the agricultural industry and affecting severely numerous developing countries. Control methods against these plant pests are sparse, the preferred one being the deployment of plant cultivars bearing resistance genes against Meloidogyne species. However, M. enterolobii is not controlled by the resistance genes deployed in the crop plants cultivated in Europe. The recent identification of this species in Europe is thus a major concern. Here, we sequenced the genome of M. enterolobii using short and long-read technologies. The genome assembly spans 240 Mbp with contig N50 size of 143 kbp, enabling high-quality annotations of 59,773 coding genes, 4,068 non-coding genes, and 10,944 transposable elements (spanning 8.7% of the genome). We validated the genome size by flow cytometry and the structure, quality and completeness by bioinformatics metrics. This ensemble of resources will fuel future projects aiming at pinpointing the genome singularities, the origin, diversity, and adaptive potential of this emerging plant pest.
Background & SummaryThe root-knot nematode Meloidogyne enterolobii (syn. M. mayaguens1) is a polyphagous species, attacking a wide range of life-sustaining and ornamental plants2. M. enterolobii is considered one of the most pathogenic and virulent root-knot nematode (RKN) species, as it is able to develop and reproduce on host plants carrying resistance to the major tropical RKN species3–5. M. enterolobii is already widespread, damaging cotton and soybean in the United States6, causing detrimental effects for watermelon production areas in Mexico7, and potato in Africa8.In 2010, M. enterolobii was added to the European Plant Protection Organisation (EPPO) A2 list and is now recommended for regulation as a quarantine species9. Based on the widespread distribution and potential to establish in the Mediterranean and subtropical regions, several countries now have designated M. enterolobii as a quarantine pest10,11. Currently, only a few potential sources of resistance such as the Ma gene from the Myrobalan plum (Prunus cerasifera)12 or in guava (Psidium spp.) and pepper accessions were reported13,14. However, these genetic resources were only tested against a single regional M. enterolobii population and not the full range of isolates from various sources and host plants.M. enterolobii is described as reproducing clonally via mitotic parthenogenesis, similarly to the most damaging RKN to worldwide agriculture (M. incognita, M. javanica, and M. arenaria). The absence of meiosis and sexual reproduction in these species was based on cytological observations in the 80’s, where no pairing of homologous chromosomes had been observed in hundreds of isolates15. Recent population genomics analysis confirmed the absence of sexual meiotic recombination in M. incognita16. The genomes of the three above-mentioned tropical RKN revealed some interesting characteristics. All three are polyploid with highly diverged genome copies most likely resulting from hybridization events17. This peculiar genome structure was shown to provide these species with diverged homoeologous gene copies that presented different gene expression patterns. These diverged gene copies might be involved in the extreme polyphagy of these species despite their asexual reproduction. More recently, it was shown that convergent gene copy losses were associated with the breaking down of a resistance gene in tomato by M. incognita18, suggesting gene copy number variations are involved in adaptive processes. M. enterolobii belongs to the same Clade I19 as the most damaging RKN species and has the same reproductive mode. However, this species is not controlled by the resistance genes deployed in plants for protection against the other RKN. Therefore, it is of great interest to explore the genome of M. enterolobii and determine similarities and differences between the species of the same RKN clade.A first draft genome of M. enterolobii from Burkina Faso, based solely on Illumina short reads, was previously available20. However, the genome assembly was quite fragmentary.Here, we used both PacBio long reads and Illumina short read sequences to study the genome of the first M. enterolobii strain officially reported in Europe. Compared to the previous draft genome of M. enterolobii, the present genome assembly is more contiguous, reducing the number of fragments by more than 10 times (from >46,000 scaffolds to <4,500 contigs) and improving the N50 length by >15x (from 9.3 kb to 143 kb).We produced RNA-seq transcriptome data that we used as a source of evidence to predict protein-coding and non-coding genes. We also annotated transposable elements (TE) as well as other repeats, and used the annotated genome to determine the ploidy level, the degree of divergence between genome copies and TE abundance. So far, this new genome is one of the most contiguous and most complete annotated one, publicly available, for a root-knot nematode species.Providing a robust reference genome for this species constitutes an important resource for further analyses with important evolutionary and agro-economic implications. This resource will pave the way for comparative as well as population genomics towards pinpointing the genome singularities, the origin, diversity, and adaptive potential of this emerging plant pest.MethodsNematode collection and DNA/RNA extractionFor genome and transcriptome sequencing, the Swiss M. enterolobii population was originally isolated from infected tomato root-stock obtained from an organic farm21. The population was maintained in a greenhouse at 25 °C and 16 h supplemental light on the cultivar Oskar F1 carrying the Mi-1 gene, conferring resistance to other Meloidogyne species (e.g. M. incognita, M. javanica and M. arenaria). We confirmed the purity of the population and correct species identification by species-specific PCR and DNA barcoding11,22 both on the bulk population and on randomly picked individual juveniles obtained from separate egg masses.Heavily infected (i.e. galled) tomato roots were collected and carefully washed free from soil to prepare the nematode material. Through incubation in a mist chamber, freshly-hatched second-stage juveniles (J2) were collected every 2 to 3 days over two weeks and used for DNA (PacBio long reads and Illumina short-read sequencing) or RNA (Illumina short read) extraction. For eggs, three galled tomato roots were extracted using the NaOCl method23. In order to obtain pure egg suspensions, density-gradient centrifugation was used to separate eggs from organic debris, according to Schaad et al.24. All J2 and egg suspensions were checked using a microscope and impurities were removed before DNA or RNA extraction.After collecting eggs and J2s, samples were frozen in liquid nitrogen. Total DNA was extracted using the DNeasy Blood & Tissue Kit (Qiagen, The Netherlands), while total RNA was extracted using the Nucleospin RNA plant kit with the protocol modified for maximum yield (Macherey-Nagel, Oensingen, Switzerland). The RNA concentration was determined with a Qubit 3.0 (Thermofisher).Genome and transcriptome sequencingThe Swiss M. enterolobii population long-fragment DNA libraries were sequenced with the PacBio RS technology at the FGCZ (ETH, Switzerland) sequencing center and were complemented by additional PacBio and Illumina short read data at GATC Biotech company. DNA extracted from nematode egg suspensions was used to generate the long-read genomic libraries. The libraries were subsequently size-selected with the BluePippin approach and then sequenced using the P6C4 chemistry. This procedure allowed us to obtain long DNA reads of up to 35 kb in length, an average length around 7 kb, and a total of ca. 13 Gigabases of long read sequence data (Table 1). Genomic Illumina short read data were obtained from J2 pre-parasitic juveniles. Transcriptomic data from eggs and J2, separately, were also produced in order to be used as a source of evidence for gene prediction. A total of one µg RNA for each stage (egg and J2) was used for sequencing with the Illumina MiSeq reagent kit V3 with 2 x 300 bp paired-end reads (Table 2).Table 1Genomic Sequence Data Statistics.DatasetMean read lengthInsert SizeNumber of readsSum (bp)Illumina Paired-end101250382,534,94638,636,029,546Illumina Mate-pair513,000430,206,97221,940,555,572PacBio7,008NA1,854,39912,996,274,301Table 2Transcriptomic Sequence Data Statistics.DatasetMean read lengthInsert Size (range)Number of readsSum (bp)Eggs Illumina Paired-end196300 (56–336)54,424,80210,660,505,616J2 Illumina Paired-end196300 (57–347)56,468,70811,083,687,102Genome assemblyPacbio reads were corrected and trimmed using sprai with default parameters (http://zombie.cb.k.u-tokyo.ac.jp/sprai/). The trimmed reads were then used as input to Canu assembler25 for a first preliminary assembly. This assembly was then used to check for contamination with the blobtools pipeline26,27. Briefly, Illumina genomic reads from both this study and the previously produced M. enterolobii draft genome20 were mapped to the genome with bwa28, and each contig was given a taxonomy affiliation based on BLAST results against the NCBI nt database (Fig. 1). Contigs that had coverage only in one Illumina library, GC percentage outside of the range of the estimated M. enterolobii GC content, and affiliation to different taxa, were considered possible contaminants. The GC estimate was calculated to be around ~30% by taking into account the GC% of nematode contigs from the blobplot analyses, the GC% of the assembly done by Szitenberg et al.20, and the GC% estimate of Meloidogyne species as shown in Table 3. Careful investigation of each such contig was then employed to limit the false positive rate. After, only the Pacbio reads that mapped to the clean contigs were retained, and a new assembly was created with Canu. Similarly, as described above, this assembly was checked for contamination, and after these two steps of contamination removal a total of 676 contigs spanning ~73 Mbp were eliminated. The clean assembly was then corrected with pilon29 using the transcriptomic reads generated in this study (–fix snps parameter), which resulted in the final frozen assembly used for downstream analyses.Fig. 1BlobPlot of the genome assembly before removing contamination. Each circle is a contig proportionally scaled by contig length and coloured by taxonomic annotation based on BLAST similarity search results. Contigs are positioned based on the GC content (X-axis) and the coverage of PacBio reads (Y-axis). There are some contigs of Proteobacteria origin at high GC and variable coverage indicating possible contamination. These contigs were removed from the assembly.Table 3Root-knot nematode genome features.Species/FeatureAssembly size (Mb)Nuclear DNA content (Mb)% complete CEGMA (C) (copies) %Partial (P)% complete BUSCO (C) % fragment (F) % missing (M)N50 (kb)# contigs/ scaffoldsGC %M. enterolobii(Swiss)*240275 ± 19C: 94.76 (3.30)P: 96.77C:87.5%F:3.6%; M:8.9%1434,43730M. enterolobii(L30)20162.4NAC: 81.45 (2.66)P: 89.92C:79.9%F:8.9%; M:11.2%9.346,09030.2M. incognita (Morelos, 2017)17183.5189 ± 15C: 94.76 (2.93)P :96.77C:88.5%F:3.0%; M:8.5%38.612,09129.8M. incognita(W1)20122NAC: 82.66 (2.34)P :89.52C:80.2%F:7.9%; M:11.9%16.533,73530.6M. incognita (Morelos, 2008)6986189 ± 15C: 74.6 (1.77)P :78.63C:71.3%F:7.3%; M:21.4%82.82,81731.4M. javanica (Avignon)17235.8297 ± 27C: 92.74 (3.68)P :95.56C:90.1%F:2.3%; M:7.6%10.431,34130M. javanica(VW4)20142.6NAC :89.52 (2.71)P: 95.16C:87.5%F:4.3%; M:8.2%14.134,39430.2M. arenaria(Guadeloupe)17258.1304 ± 9C :94.76 (3.66)P :95.56C:87.1%F:4.3%; M:8.6%16.526,19630M. arenaria(A2-O)71284.05NAC: 94.76 (3.57)P :96.77C:87.1%F:2.6%; M:10.3%204.62,22430M. arenaria(HarA)20163.8NAC :91.53 (2.74)P :95.97C: 78.2F :12.2; M:9.610.546,50930.3M. hapla(VW9)6553.6121 ± 3C :93.55 (1.19)P: 95.56C:87.4%F:4.3%; M:8.3%83.61,52327.4M. floridensis(JB5)7299.9NAC: 56.45 (1.95)P: 76.61C:54.1%F:28.7%; M:17.2%3.581,11129.7M. floridensis(SJF1)2074.9NAC: 77.42 (1.71)P: 83.87C:76.5%F:7.6%; M:15.9%13.39,13430.2M. luci(SI-Smartno)73209.2NAC: 95.56 (2.92)P :96.77C:87.8%F:4.0%; M:8.2%1,71232730.2M. graminicola(IARI)7438.18NAC: 84.27 (1.34)P: 90.73C:73.6%F:15.2%; M:11.2%20.44,30423.05*This work.Transcriptome assemblyAdapters, low-quality regions (Phred score <30), and regions with two or more consecutive ambiguous nucleotides, were cropped using prinseq30, resulting in 56,468,708 and 54,424,802 cleaned paired-end reads for J2 and egg libraries, respectively. The clean reads were assembled using CLC Genomics Workbench 9.0 (https://www.qiagenbioinformatics.com/) with automatic bubble size and word size estimation, in simple contig mode, and minimum contig length of 200. The transcriptome assemblies consisted of 110,068 and 110,263 contigs for J2 (J2 stage transcriptome assembly of Meloidogyne enterolobii (INRA/JKI). available in Figshare31) and Egg (Egg transcriptome assembly of Meloidogyne enterolobii (INRA/JKI). available in Figshare31) libraries, respectively.Gene prediction and annotationDetection of gene models was done with the fully automated pipeline EuGene-EP version 132. EuGene has been configured to integrate similarities with known proteins from Wormpep 22133, Meloidogyne incognita predicted proteome17 and UniProtKB/Swiss-Prot database34, with the prior exclusion of proteins that were similar to those present in RepBase35.Three datasets of Meloidogyne enterolobii transcribed sequences were aligned on the genome and used by EuGene as transcription evidence: (i) the de novo assemblies of the egg and J2 transcriptomes (ii) the egg stage RNAseq clean reads and (iii) the J2 stage RNAseq clean reads. The alignments of dataset (i) on the genome, spanning 50% of the transcript length with at least 95% identity were retained. The alignments of datasets (ii) and (iii) spanning 90% of the read length with 97% identity were retained.The EuGene default configuration was edited to set the “preserve” parameter to 1 for all datasets, the “gmap_intron_filter” parameter to 1, the minimum intron length to 35 bp, and to allow the non canonical donor splice site “GC”. Finally, the Nematode specific Weight Array Method matrices were used to score the splice sites (http://eugene.toulouse.inra.fr/Downloads/WAM_nematodes_20171017.tar.gz).Functional annotation of predicted proteinsWe analysed all the predicted proteins in order to identify their putative functions and the putative secretory subset.Prediction of conserved protein domains and gene ontology assignmentWe used PfamScan36 on the whole set of M. enterolobii proteins against the Pfam-A v32 library of HMM profiles37 to identify conserved protein domains. We used the Pfam annotation as well as the set of predicted proteins as an input to BLAST2GO pro v.1.12.1138 to assign standard and ‘slim’ gene ontology terms according to the presence of Pfam domains, to BLASTp39 homology searches against the NCBI’s nr database and to Interproscan40 annotation performed internally in BLAST2GO.Secretome predictionSignal peptides for secretion were detected in the set of M. enterolobii predicted proteins using Signalp4.141 and cleaved. From the cleaved proteins, transmembrane regions were predicted using TMHMM2.0c42. We considered as possibly secreted all the proteins that featured a predicted signal peptide and no transmembrane region.Prediction and annotation of transposable elementsPrediction and annotation of transposable elements (TE) were performed with the REPET meta-pipeline, which combines TEdenovo and TEannot pipelines43.Data pre-processingTo correctly perform the all against all comparisons of the genome contigs step (see below), REPET automatically trims stretches of Ns of length 11 or more. However, this can result in a higher fragmentation of the genome with a risk of identification of spurious short repetitive matches. The M. enterolobii genome comprises few unresolved ‘N’ bases (1,401). We first created a modified version of the genome by splitting it at N stretches of length 11 or more and then trimming all N, using dbchunk.py from the REPET tools collection. To circumvent the potential problem with shorter DNA fragments, we then only kept chunks of length above 8,659 bp, defined as the L99 chunk length threshold.TE de-novo predictionThe previously pre-processed genome was used to build a TE consensus library using the TEdenovo pipeline (REPET 2.5 configuration files used for TE prediction and annotation in the M. enterolobii genome (JKI/INRA). available in Figshare31). The genome was aligned to itself using Blaster44 and High Scoring Segment Pairs (HSPs) were detected. The repetitive HSPs were clustered by Recon45 and Grouper44. A consensus sequence was created for each cluster based on the corresponding multiple alignment with MAP46. Consensus sequences were analysed by PASTEClassifier47 in order to find homology with known TE and to detect TE-related features such as specific HMM-profiles and structural characteristics. For identification of TE features, we used the curated libraries provided by the REPET development team (URGI): Repbase 20.05 (aa and nt), a concatenation of Pfam 27.0 and GyDB 2.0 hmm-profiles databanks and rRNA Eukaryota databank. The detected features were used by PASTEClassifier to classify consensus sequences in the different TE orders defined in Wicker’s classification48. Simple sequence repeats (SSR) and under-represented unclassified consensuses were filtered out.Semi-automated TE consensus filteringThe TE consensus library was filtered as follows in order to minimise false positives and redundancy. A draft annotation of the whole genome based on the consensus library was performed using TEannot (steps 1, 2, 3, 7; (REPET 2.5 configuration files used for TE prediction and annotation in the M. enterolobii genome (JKI/INRA) available in Figshare31). Only consensus sequences with at least one Full-Length Copy (FLC) annotated on the genome were retained to constitute a new filtered TE consensus library (M. enterolobii cleaned TE consensus sequence library available in Figshare31). Consensuses with at least one FLC annotated on the genome were identified with PostAnalyzeTELib.py and the corresponding sequences were extracted using GetSequencesFromAnnotations.py (REPET tools).TE whole genome annotationThe TEannot pipeline was used to annotate the whole M. enterolobii genome from the previously created ‘filtered TE consensus library’ (steps 1, 2, 3, 4, 5, 7, 8; (REPET 2.5 configuration files used for TE prediction and annotation in the M. enterolobii genome (JKI/INRA). available in Figshare31). TE consensus sequences from the ‘filtered TE consensus library’ were aligned on the genome using Blaster, Censor49, and RepeatMasker. The results of the three methods were concatenated and MATCHER44 was used to remove overlapping HSPs and make connections with the “join” procedure. SSRs were detected by TRF50, Mreps51, and RepeatMasker, and then merged. Eventually, after some redundancy removal and application of the “long join” procedure (distant fragments belonging to the same copies are joined), raw annotations (Unfiltered transposable elements annotation of M. enterolobii genome (INRA/JKI). available in Figshare31) were exported for post-processing.Annotations filtering and post-processingWe used in-house Python scripts (Python script used to filter TE annotations in the M. enterolobii genome (INRA/JKI). available in Figshare31) to extract canonical TE annotations from the whole repeatome annotation, but also increase the consistency between TE consensus sequences and their annotated copies.Only TE annotation with >85% identity to the consensus sequence, a length >250 nucleotides, and classified as retrotransposon (Wicker’s class I) or DNA-transposon (Wicker’s class II) were retained. TE-annotations covering less than 1/3 of the consensus sequence length were discarded. Each genomic sequence corresponding to TE annotations were individually blasted against the TE consensus library, and only the ones with their consensus as the best hit were retained. Eventually, overlapping annotations on the same strand were removed. This ensemble of filters yielded the final TE annotation (Filtered final transposable elements annotation of M. enterolobii genome (INRA/JKI). available in Figshare31) used to compute statistics on the percentage of the genome occupied by TE and relative repartition in TE orders.These scripts also allowed computing basic statistics and to describe TE coverage on the genome regarding Wicker’s classification.Data RecordsAll the Illumina and PacBio sequence data supporting the results of this article as well as the genome assembly and gene models have been deposited and are publicly available at the EMBL-EBI’s European Nucleotide Archive52 and at the NCBI53,54. All the processed data, including genome and transcriptome assemblies and all the annotation results have been deposited and are publicly available as a Figshare collection31. The genome assembly as well as the predicted gene models, a blast server and a genome browser are publicly available at https://meloidogyne.inrae.fr/.Technical validationDe novo genome assemblyAfter elimination of the few contaminant contigs, we assembled the M. enterolobii Swiss strain in 4,437 contigs for a total genome size of 240 Mb with a contig N50 length of 143 kb (Meloidogyne enterolobii genome assembly (V1, INRA/JKI). available in Figshare31). Overall, the genome assembly we produced constitutes a significant improvement in terms of contiguity and completeness compared to the previous M. enterolobii draft genome from a Burkina Faso isolate20. The genome size grew from 162.4 to 240 Mb. The number of contigs / supercontigs was divided by >10 (46,090 to 4,437) while the N50 length was multiplied by >15 (9.3 to 143 kb). This N50 length is in the top 3 highest for a Meloidogyne genome and the best for a publicly available annotated one (Table 3).Validation of the genome sizeWe used flow cytometry to perform accurate measurement of cells DNA contents in M. enterolobii compared to internal standards with known genome sizes: Caenorhabditis elegans strain Bristol N2 (approximately 200 Mb at diploid state) and Drosophila melanogaster strain Cantonese S. (approximately 350 Mb at diploid state) as previously described17. Briefly, nuclei were extracted from two hundred thousand J2 infective juvenile larvae as described in55 and stained with 75 µg/mL propidium iodide and 50 µg/mL DNAse-free RNAse. Flow cytometry analyses were carried out using an LSRII / Fortessa (BD Biosciences) flow cytometer operated with FACSDiva v6.1.3 (BD Biosciences) software. The DNA contents of the M. enterolobii samples were calculated by averaging the values obtained from three biological replicates. For estimation of total nuclear DNA content, we used the M. enterolobii strain Godet (from Guadeloupe France; N°75 available in Institut Sophia Agrobiotech collection).The calculated nuclear DNA content via flow cytometry experiments was 274.69 ± 18.52 Mb (Fig. 2). With 240 Mb, the de novo genome assembly size is close to the estimated total DNA content in M. enterolobii cells and also represents a substantial improvement compared to the previous genome assembly (162.4 Mb, Table 3)20.Fig. 2Relative DNA staining in nuclei of M. enterolobii. Cytogram example obtained after gating on G0/G1 nuclei (arbitrary units) from M. enterolobii when processed mixed together with C. elegans, the internal standard (diploid genome size is 200 Mb).This suggests that most of the M. enterolobii genome has been captured in this assembly. The difference between the genome assembly size and the total estimated DNA content ranges from 16 to up to 53 Mb and could be due either to duplicated and repetitive genome regions that have not been correctly separated during assembly or to differences in genome sizes between the Swiss M. enterolobii strain we have sequenced and the ‘Guadeloupe’ strain used for DNA content measurement via flow cytometry.Genome completeness assessmentTo assess the completeness of the genome assembly, we ran CEGMA 2.556 and BUSCO v3.057 with the Eukaryotic dataset (odb9) in fast mode and C. elegans as a model for Augustus predictions58. Both pipelines search in genome assemblies for genes universally or largely conserved in Eukaryotes and produce reports on the number of genes found in complete, partial, single copy or duplicated versions in the genome under consideration. Although a nematode dataset exists in BUSCO, it only includes 8 genomes from just three of the 12 described nematode clades (2, 8 and 9)59. Because the root-knot nematodes belong to Clade 12, which is not represented at all, we decided to use the Eukaryotic dataset which is more comprehensive (65 species, including the 8 nematodes).Overall, 94.76 and 87.5% of the CEGMA and eukaryotic BUSCO genes, respectively, were found in complete length.For comparison purposes, we also ran CEGMA and BUSCO with the same parameters on all the root-knot nematode genome assemblies that are publicly available (Table 3) The CEGMA and BUSCO completeness scores are both among the highest obtained for a Meloidogyne genome to date. For comparison, the model nematode C. elegans reaches CEGMA and eukaryotic BUSCO completeness scores of 96.37 and 94.7, respectively. This is higher than all the Meloidogyne genomes but it should be noted that the whole protein set of C. elegans is part of the training set for both CEGMA and BUSCO.Gene predictionUsing the Eugene-EP pipeline, 63,841 genes were predicted, of which 59,773 were protein-coding with the exons spanning 61.9 Mb (~26%) of the genome assembly length. The GC percentage was higher in the coding portion (34.3%) than in the whole genome (30.0%). Spliceosomal introns were detected in 88% of protein-coding genes, with an average of 6.2 exons per gene. Similarly to the other tropical RKN genomes17, less than 1% of splice sites have a non-canonical GG donor dinucleotide (canonical is GT).Eugene-EP also predicted 4,068 non protein-coding genes (e.g. ribosomal, tRNA, splice leader genes). None of them had predicted intron but similarly to protein-coding genes the average GC content was higher (34.1%) than for the rest of the genome.The overall statistics of the EUGENE gene predictions (Statistics of the gene predictions of M. enterolobii (genome V1, INRA/JKI).) and GFF3 annotations (GFF3 of the EUGENE gene prediction of M. enterolobii (genome V1, INRA/JKI).) as well as the fasta files of the predicted protein-coding genes (Predicted protein-coding genes of M. enterolobii (genome V1, INRA/JKI).), coding sequences (Predicted coding sequences (CDS) in M. enterolobii (genome V1, INRA/JKI).), proteins (Predicted proteins of M. enterolobii (genome V1, INRA/JKI).), messenger RNAs (Predicted messenger RNAs of M. enterolobii (genome V1, INRA/JKI).) and non-coding RNAs (Predicted non coding RNAs (ncRNA) of M. enterolobii (genome V1, INRA/JKI).) are all available in Figshare31.Validation of the predicted proteinsWe ran a BUSCOV3 analysis in protein mode with the eukaryotic (odb9) dataset of 303 highly conserved genes to assess the completeness and quality of the set of predicted proteins in M. enterolobii and compared to the previous version of the M. enterolobii predicted proteins from the Burkina Faso isolate. Overall, 94.7% (287/303) of the eukaryotic BUSCO genes were found in complete length in the M. enterolobii proteome. This is a substantial improvement compared to the previously available M. enterolobii proteome, in which only 78.2% of complete eukaryotic BUSCO genes were found (Table 4). This improvement was due to both a reduction of the number of fragmentary BUSCO genes (from 11.2% to 3.0%) and of the number of missing BUSCO genes (from 10.6% to 2.3%).Table 4Estimation of the protein set completeness with BUSCO.BUSCO genes (303)M. enterolobii (Swiss)*M. enterolobii (L30)20Complete94.7% (287)78.2% (237)Single-copy complete9.2% (28)45.9% (139)Duplicated complete85.5% (259)32.3% (98)Fragmentary3.0% (9)11.2% (34)Missing2.3% (7)10.6% (32)*This work.Predicted secreted proteinsRKN and other plant-parasitic nematodes secrete in planta, proteins called effectors, where they support parasitism by manipulating plant structure and functions60. We assessed whether the set of predicted proteins in M. enterolobii could support effector discovery. To do so, we studied whether some gene ontology terms enriched in predicted secreted proteins corresponded to obvious known effector functions. We identified 4,973 M. enterolobii proteins (~8.3% of the total) as possibly secreted by having a predicted signal peptide for secretion and no transmembrane domain (List of predicted secreted proteins of M. enterolobii (genome V1, INRA/JKI). available in Figshare31).Using BLAST2GO pro functional annotation (Functional annotation of predicted proteins in M. enterolobii (genome V1, INRA/JKI). available in Figshare31), we investigated whether some functions were enriched in the set of predicted secreted proteins. We used a gene set enrichment analysis (GSEA) as described in61 to identify significantly overrepresented gene ontology (GO) terms in the set of predicted secreted proteins in comparison to the rest of the proteins. We considered as significantly overrepresented, the GO terms that returned a false discovery rate (FDR) value < 0.05 in one-tailed Fisher’s exact test. We identified 49 overrepresented GO ‘slim’ terms in the predicted secreted proteins (Over-represented GO terms in predicted secreted proteins M. enterolobii (genome V1, INRA/JKI). available in Figshare31).In total, 20, 24 and 5 enriched GO terms were identified in the ‘biological process’ (BP), ‘molecular function’ (MF) and ‘cellular component’ (CC) ontologies, respectively. In the BP category, enriched terms encompassed several enzymatic / catabolic processes in relation with plant cell wall macromolecules (e.g. ‘cell wall macromolecule catabolic process’, ‘cellulose catabolic process’) or other macromolecules (e.g. ‘chitin metabolic process’, ‘peptidoglycan catabolic process’). These terms echo enriched terms in the MF category such as ‘pectate lyase activity’, ‘polysaccharide binding’, ‘cellulase activity’, all related to degradation / modification of the plant cell wall. Plant cell wall-degrading enzymes are one category of known effectors secreted by plant-parasitic nematodes with the clearest function62,63. The above-mentioned enriched GO terms in the BP and MF ontologies, validate the set of predicted secreted proteins an interesting resource towards discovery and description of M. enterolobii effectors.As expected for predicted secreted proteins, in the CC category, the two most significantly enriched GO terms were ‘extracellular space’ and ‘endoplasmic reticulum lumen’.Genome structure and ploidy levelThe genomes of M. incognita, M. javanica and M. arenaria, other mitotic parthenogenetic RKN from Clade I, have been described as allopolyploids17. Because M. enterolobii belongs to the same clade and is also described as a mitotic parthenogenetic species, it is important to estimate the ploidy level of the genome assembly.In the genomes of M. incognita, M. javanica and M. arenaria, respectively described as triploid and degenerate tetraploids17, the CEGMA genes were found in 2.93, 3.68 and 3.66 copies, on average, respectively (Table 3). In contrast, the genomes of M. hapla and C. elegans have an average number of CEGMA gene copies of 1.19 and 1.09, respectively, which is consistent with their homozygous diploid nature and the corresponding haploid genome assemblies. Thus, the average number of CEGMA gene copies in RKN genomes seem to recapitulate the estimated ploidy levels. In M. enterolobii, the average number of copies per CEGMA gene was 3.3, suggesting a polyploid (at least triploid) genome structure as well. To further investigate and validate the ploidy level, we performed additional analyses using gene predictions as markers.Duplication level estimation via mapping of coding sequences (CDS)We used BLAT64 to align the ensemble of predicted coding sequences (CDS) to the respective reference genomes of M. enterolobii, M. incognita Morelos17 and the meiotic facultative sexual species M. hapla VW965. We only retained genome matches that covered at least 2/3 of the CDS length at a minimum 95% identity, and counted the number of matched loci per CDS query. In the diploid meiotic species M. hapla, 90% of the CDS mapped to a unique locus in the reference genome and can be considered single-copy genes. In contrast, only 12 and 13% of M. enterolobii and M. incognita CDS map a unique locus on their respective genome assemblies. Hence, the vast majority of protein-coding genes are in multiple copies in M. enterolobii and M. incognita (Fig. 3). Furthermore, both M. enterolobii and M. incognita show a peak of CDS alignments at 3 different loci. This is consistent with the average CEGMA gene copy number of 2.93 and 3.3 for M. incognita and M. enterolobii, respectively, and reinforces the possibility of triploidy.Fig. 3Number of CDS mapping n loci in Meloidogyne genomes. Number of CDS mapping with >95% identity on >66% on their length at n loci on the M. enterolobii (red), M. incognita (green) and M. hapla (blue) genomes.Genome structure estimation by classification of gene duplicatesAlthough CEGMA gene copy numbers and alignment of CDS to the genome suggest most of the genes are triplicated, this does not necessarily imply a triploid genome. To further validate ploidy, it is important to show that gene duplicates are forming whole duplicated blocks and not only dispersed independent duplications.We used MCScanX66 to detect and classify duplications and study the genome structure (MCScanX duplication analysis results on the M. enterolobii genome (INRA/JKI). available in Figshare31). In a first step, all the predicted protein sequences were self-blasted to determine homologous relationships between them, with an e-value threshold of 1E−10. Using homology information from the all-against-all blast output and gene location information from the GFF3 annotation file, MCScanX detects duplicated protein-coding genes and classifies them in the following categories, (i) singleton when no duplicates are found in the assembly, (ii) proximal when duplicates are on the same contig and separated by 1 to 10 genes, (iii) tandem when duplicates are consecutive, (iv) whole genome duplication (WGD) or segmental when duplicates form collinear blocks with other pairs of duplicated genes, and (v) dispersed when the duplicates cannot be assigned to any of the above-mentioned categories. For detection of WGD / segmental duplications we required at least 3 collinear gene pairs.MCScanX analyses showed that 94.8% of the protein-coding genes are duplicated in M. enterolobii. Overall, the majority of the genes (61.8%) are part of the whole genome duplication category and form 2,892 collinear blocks, that span 164 Mb (ca. 68% of the genome). This reinforces the idea that the genome is polyploid. Then, 27.6% of the genes have been classified as dispersed duplicates, 3.1% and 2.2% form proximal or tandem duplicates, respectively. By comparison, in M. incognita, only 28.5% of the genes formed duplicated collinear blocks and 59.4% were dispersed copies. This highlights the gain of resolution allowed by improved genome assemblies. Indeed, the proportion of genes forming whole duplicated blocks seems to be correlated to the genomes N50 lengths. We investigated the duplication depth of genes in the collinear blocks and observed a peak at a depth of three (Fig. 4), which is consistent with the peak at three loci matched by the CDS on the genome and the average copy number of CEGMA genes. This ensemble of results strongly suggests that the genome of M. enterolobii is triploid.Fig. 4Distribution of genes based on the duplication depth in the collinear blocks. The peak of the distribution is observed at a depth of three, suggesting a triploid genome.Divergence level between duplicated genome blocksIn the mitotic tropical RKN M. incognita, M. javanica and M. arenaria, the duplicated genomic regions have a high pairwise average nucleotide divergence of ~8%17. To check whether a similar high nucleotide divergence existed between the M. enterolobii genome copies, we used nucmer67 and aligned the 2,892 duplicated blocks identified by MCScanX (WGD/segmental category). The average pairwise nucleotide identity between duplicated blocks was 92%, thus validating a similar 8% nucleotide divergence. The distribution of pairwise identities between the M. enterolobii duplicated genomic blocks shows a major peak at 89.9% identity and a second minor peak at 94.2% (Fig. 5). This distribution is different from that of M. incognita, where almost two overlapping peaks at 89.5% and 91.4% identity are observed. This suggests despite similar average divergence, different evolutionary trajectories in M. incognita and M. enterolobii and independent polyploidization events.Fig. 5Nucleotide identity computed from NUCmer alignments between homologous duplicated blocks. M. enterolobii (blue) shows one major peak at 89.9% and two smaller ones at 94.2% and 99.9% nucleotide identity, whereas M. incognita (green) shows almost two overlapping peaks at 89.5% and 91.4% of identity.Transposable and other repetitive elementsRepetitive elements span 47.62% of the M. enterolobii genome assembly size, but canonical Transposable Elements (TE) annotations occupy only 8.70% (Fig. 6, Table 5). Using Wicker’s classification48, Class-I Retrotransposons span 2.87% of the genome assembly while class-II DNA-transposons occupy a twice higher proportion (5.83%). Using the same methodology, a much higher proportion of DNA transposons compared to retrotransposons was also observed in the genomes of M. incognita and C. elegans68. Because TE annotation is highly dependent on genome assembly quality and on the methodology used, we will not make direct comparison of the percentage of the genome occupied by TE with other nematode genomes. Nonetheless, re-annotating the C. elegans genome with the same methodology as a control yielded a very similar proportion of the genome occupied by TE than reported in previous studies68.Fig. 6Transposable Elements annotation in the genome of M. enterolobii. Canonical-TE annotations distribution as a genome percentage. Annotations are grouped according to Wicker’s classification48.Table 5Summary of M. enterolobii Canonical-TE annotations statistics.ordernb. of featuressum of features length (bp)median features length (bp)median identity with consensus (%)Retro-transposonsCLASS_1_LIKE86624,8687,19798.65LTR6473,329,7014,01297.55DIRS58494,4838,323.597.9LINE3431,379,8843,79897.1SINE3636,9061,07699.5LARD108233,7971,795.593.83TRIM469794,64469097.9DNA-transposonsCLASS_2_LIKE46184,4741,193.599TIR3,9744,415,24291896.95HELITRON179150,7648,32898.9MAVERICK5465,253,8128,95497.3MITE4,452264,54459096.4Total10,94417,163,119As a further validation of the accuracy of TE annotation, we observed that all the previously reported TE orders in RKN genomes17,69–74 (i.e. class-I Retrotransposons:DIRSs, LINEs, LTRs, SINEs; and class-II DNA-transposons:HELITRONs, TIRs, MAVERICKs, and MITEs) were detected in M. enterolobii as well. Furthermore, TRIM and LARD non-autonomous retrotransposons, previously reported in RKN genomes as ‘unclassified’17, were also now correctly identified and classified in the genome of M. enterolobii.