Transcriptome Analysis and miRNA Target Profiling at Various Stages of Root-Knot Nematode Meloidogyne incognita Development for Identification of Potential Regulatory Networks

Root-knot nematodes (RKNs) are a group of plant-parasitic nematodes that cause damage to various plant species and extensive economical losses. In this study, we performed integrated analysis of miRNA and mRNA expression data to explore the regulation of miRNA and mRNA in RKNs. In particular, we aimed to elucidate the mRNA targets of Meloidogyne incognita miRNAs and variations of the RKN transcriptome during five stages of its life cycle. Stage-wise RNA sequencing of M. incognita resulted in clean read numbers of 56,902,902, 50,762,456, 40,968,532, 47,309,223, and 51,730,234 for the egg, J2, J3, J4, and female stages, respectively. Overall, stage-dependent mRNA sequencing revealed that 17,423 genes were expressed in the transcriptome of M. incognita. The egg stage showed the maximum number of transcripts, and 12,803 gene transcripts were expressed in all stages. Functional Gene Ontology (GO) analysis resulted in three main GO classes: biological process, cellular components, and molecular function; the detected sequences were longer than sequences in the reference genome. Stage-wise selected fragments per kilobase of transcript per million mapped reads (FPKM) values of the top 10 stage-specific and common mRNAs were used to construct expression profiles, and 20 mRNAs were validated through quantitative real-time PCR (qRT-PCR). Next, we used three target prediction programs (miRanda, RNAhybrid, and PITA) to obtain 2431 potential target miRNA genes in RKNs, which regulate 8331 mRNAs. The predicted potential targets of miRNA were generally involved in cellular and metabolic processes, binding of molecules in the cell, membranes, and organelles. Stage-wise miRNA target analysis revealed that the egg stage contains heat shock proteins, transcriptional factors, and DNA repair proteins, whereas J2 includes DNA replication, heat shock, and ubiquitin-conjugating pathway-related proteins; the J3 and J4 stages are represented by the major sperm protein domain and translation-related proteins, respectively. In the female stage, we found proteins related to the maintenance of molybdopterin-binding domain-containing proteins and ubiquitin-mediated protein degradation. In total, 29 highly expressed stage-specific mRNA-targeting miRNAs were analyzed using qRT-PCR to validate the sequence analysis data. Overall, our findings provide new insights into the molecular mechanisms occurring at various developmental stages of the RKN life cycle, thus aiding in the identification of potential control strategies.


Introduction
Plant-parasitic nematodes are the most destructive group of plant pathogens. Plantparasitic nematodes are extremely challenging to control and cause global agricultural losses of approximately USD 157 billion annually. Meloidogyne incognita is widely distributed in temperate and tropical regions worldwide; it can infect thousands of plant

Stage-Wise Transcriptome Profiling of M. incognita
To acquire high-throughput transcriptome data for M. incognita, we conducted Illumina-based next-generation sequencing. In this study, we performed RNA sequence analysis of all five developmental stages of M. incognita. We used three different strategies: RNA sequencing of M. incognita at all developmental stages, targeting of mRNA using a bioinformatics pipeline based on previously identified miRNA, and functional annotation using Blast2GO. The workflow of the current study is illustrated in Figure 1. Figure 1. Flowchart representing the overall process of the current study. In silico analysis of miRNA-targeted mRNA prediction from the reference genome (ASM18041V1a) [2]; further combinations of miRNA-targeted mRNAs were elucidated. Illumina sequence data were used as the raw data; these data were cleaned using multiple software programs (RNAhybrid, miRanda, and PITA) and differentially expressed mRNAs were analyzed with a Venn diagram and heatmap. Expression levels of miRNA and mRNA based on sequencing results were normalized using DESeq2 and the fragments per kilobase of transcript per million mapped reads (FPKM) methods, respectively. Further miRNA-targeted mRNAs were validated through qRT-PCR for all developmental stages of M. incognita. Sequencing results were annotated with BLAST2GO using data for specific, common, and unique miRNA-targeted mRNAs.
Deep RNA sequencing achieved the following numbers of raw reads: egg, 70,310,943; J2, 61,150,827; J3, 51,342,127; J4, 58,016,948; and female stage, 64,618,511. From the raw reads, further exclusion of non-informative sequences was conducted using CASAVA (1.8.2) and CLC genome cell (4.0). Finally, we obtained the following numbers of highquality reads: egg, 56,902,920; J2, 50,762,456; J3, 40,968,532; J4, 47,309,223; and female stage, 51,730,234. Furthermore, we used a bioinformatics pipeline for miRNA target prediction with a reference genome. Initial preprocessing followed by mapping of sequences to the reference genome indicated that 51.31-61.12% of the high-quality sequences could be mapped to the currently available reference genome for M. incognita (Table 1). Through further analysis, we generated a dataset supporting transcription of 17,423 (85.55%) of the predicted gene models (Table S1). In silico analysis of miRNA-targeted mRNA prediction from the reference genome (ASM18041V1a) [2]; further combinations of miRNA-targeted mRNAs were elucidated. Illumina sequence data were used as the raw data; these data were cleaned using multiple software programs (RNAhybrid, miRanda, and PITA) and differentially expressed mRNAs were analyzed with a Venn diagram and heatmap. Expression levels of miRNA and mRNA based on sequencing results were normalized using DESeq2 and the fragments per kilobase of transcript per million mapped reads (FPKM) methods, respectively. Further miRNA-targeted mRNAs were validated through qRT-PCR for all developmental stages of M. incognita. Sequencing results were annotated with BLAST2GO using data for specific, common, and unique miRNA-targeted mRNAs.
Deep RNA sequencing achieved the following numbers of raw reads: egg, 70,310,943; J2, 61,150,827; J3, 51,342,127; J4, 58,016,948; and female stage, 64,618,511. From the raw reads, further exclusion of non-informative sequences was conducted using CASAVA (1.8.2) and CLC genome cell (4.0). Finally, we obtained the following numbers of highquality reads: egg, 56,902,920; J2, 50,762,456; J3, 40,968,532; J4, 47,309,223; and female stage, 51,730,234. Furthermore, we used a bioinformatics pipeline for miRNA target prediction with a reference genome. Initial preprocessing followed by mapping of sequences to the reference genome indicated that 51.31-61.12% of the high-quality sequences could be mapped to the currently available reference genome for M. incognita (Table 1). Through further analysis, we generated a dataset supporting transcription of 17,423 (85.55%) of the predicted gene models (Table S1).  [20].
The stage-wise expression patterns of gene transcripts with FPKM ≥ 0.3 were assessed in the egg, J2, J3, J4 and female stages; 15,798 15,564, 14,908, 15,226, and 14,519 transcripts were identified, respectively ( Table 2). The predicted numbers of genes were similar across all developmental stages. In terms of mRNA expression, the maximum stage-specific transcript number observed was 390 in the egg stage; 12,803 gene transcripts were commonly expressed across all stages. In total, 546 gene transcripts, the highest number shared among any four stages, were expressed in the four stages other than female stage. In the J2, J3, J4, and female stages, 326, 141, 115, and 120 transcripts were exclusively expressed, respectively. The remaining gene transcripts were shared among combinations of two, three, or four stages (e.g., egg + J2, egg + J3, egg + J4, and egg + female). The greatest number of mRNAs shared between two stages was 290 for the egg and J2 stages. No mRNAs were shared exclusively between the egg and J4 or the J2 and female stages. The numbers of gene transcripts shared among various stages are presented in a Venn diagram ( Figure 2 Table 2). The predicted numbers of genes were similar across all developmental stages. In terms of mRNA expression, the maximum stage-specific transcript number observed was 390 in the egg stage; 12,803 gene transcripts were commonly expressed across all stages. In total, 546 gene transcripts, the highest number shared among any four stages, were expressed in the four stages other than female stage. In the J2, J3, J4, and female stages, 326, 141, 115, and 120 transcripts were exclusively expressed, respectively. The remaining gene transcripts were shared among combinations of two, three, or four stages (e.g., egg + J2, egg + J3, egg + J4, and egg + female). The greatest number of mRNAs shared between two stages was 290 for the egg and J2 stages. No mRNAs were shared exclusively between the egg and J4 or the J2 and female stages. The numbers of gene transcripts shared among various stages are presented in a Venn diagram ( Figure 2).

Functional Annotation and Comparison Among Stages
To obtain functional annotations, we have utilized the Blast2GO tool suite-a comprehensive bioinformatics tool for functional annotation of sequencing data. BLAST2GO

Functional Annotation and Comparison among Stages
To obtain functional annotations, we have utilized the Blast2GO tool suite-a comprehensive bioinformatics tool for functional annotation of sequencing data. BLAST2GO analysis revealed that the results of our mRNA sequencing experiment were similar to the available mRNA sequences for M. incognita (Project ID: PRJEA28837). We compared the top Gene Ontology (GO) distributions in three GO categories (biological process, molecular function, and cellular component) to the reference mRNA sequence, and the most common comparison sequence hits between ours and the reference study are as follows: In the biological process category, the cellular process (6058 vs. 4610), single-organism process (3645 vs. 5117), and metabolic process (4123 vs. 5301) were most strongly represented ( Figure 3A). In the molecular function category, binding (5139 vs. 4368) and catalytic activity (4478 vs. 3455) were most strongly represented ( Figure 3B); in the cellular component category, cell (5125 vs. 3705), organelle (3326 vs. 2619), and membrane (2569 vs. 1714) were predominant ( Figure 3C). The fewest sequence hits were for electron carrier (43 vs. 29) and antioxidant activity (58 vs. 33), both in the molecular function category ( Figure 3B). Comparison of enzymes encoded by the identified mRNAs indicated that similar families of enzymes were expressed, except for higher numbers of oxidoreductases (342 vs. 203) and transferases (653 vs. 565), compared with the reference mRNA sequence (Project ID: PRJEA28837). Among hydrolases, 934 enzymes were expressed in the reference study, but this number decreased to 583 in the present study. Overall, the families of enzymes expressed were similar between studies, but more transcripts were obtained in this study ( Figure S1). We annotated the sequences using the available InterProScan database as a reference to identify common and specific targets. Overall, InterProScan results indicated that fewer InterProScan (516) and GO (282) data were available for our contigs, compared with the reference dataset (Table S2). Gene functional annotation analyses of common and stage-wise expression, which are summarized in Figure S2, indicated that most of the genes were uncharacterized proteins from diverse species, including C. elegans, Caenorhabditis briggsae, Aedes aegypti, Herpetosiphon aurantiacus, and Vitis vinifera. Comparison between the reference dataset and our study showed that other functional proteins were also shared. Furthermore, the selected highly expressed mRNAs represented an uncharacterized protein, FMRF-amide-like peptide, cuticle collagen, and an uncharacterized protein in the egg, J2, J3, J4, and female stages, respectively. Highly expressed mRNAs from each stage of M. incognita are listed in Table 3. analysis revealed that the results of our mRNA sequencing experiment were similar to the available mRNA sequences for M. incognita (Project ID: PRJEA28837). We compared the top Gene Ontology (GO) distributions in three GO categories (biological process, molecular function, and cellular component) to the reference mRNA sequence, and the most common comparison sequence hits between ours and the reference study are as follows: In the biological process category, the cellular process (6058 vs. 4610), single-organism process (3645 vs. 5117), and metabolic process (4123 vs. 5301) were most strongly represented ( Figure 3A). In the molecular function category, binding (5139 vs. 4368) and catalytic activity (4478 vs. 3455) were most strongly represented ( Figure 3B); in the cellular component category, cell (5125 vs. 3705), organelle (3326 vs. 2619), and membrane (2569 vs. 1714) were predominant ( Figure 3C). The fewest sequence hits were for electron carrier (43 vs. 29) and antioxidant activity (58 vs. 33), both in the molecular function category ( Figure 3B). Comparison of enzymes encoded by the identified mRNAs indicated that similar families of enzymes were expressed, except for higher numbers of oxidoreductases (342 vs. 203) and transferases (653 vs. 565), compared with the reference mRNA sequence (Project ID: PRJEA28837). Among hydrolases, 934 enzymes were expressed in the reference study, but this number decreased to 583 in the present study. Overall, the families of enzymes expressed were similar between studies, but more transcripts were obtained in this study ( Figure S1). We annotated the sequences using the available InterProScan database as a reference to identify common and specific targets. Overall, InterProScan results indicated that fewer InterProScan (516) and GO (282) data were available for our contigs, compared with the reference dataset (Table S2). Gene functional annotation analyses of common and stage-wise expression, which are summarized in Figure S2, indicated that most of the genes were uncharacterized proteins from diverse species, including C. elegans, Caenorhabditis briggsae, Aedes aegypti, Herpetosiphon aurantiacus, and Vitis vinifera.
Comparison between the reference dataset and our study showed that other functional proteins were also shared. Furthermore, the selected highly expressed mRNAs represented an uncharacterized protein, FMRF-amide-like peptide, cuticle collagen, and an uncharacterized protein in the egg, J2, J3, J4, and female stages, respectively. Highly expressed mRNAs from each stage of M. incognita are listed in Table 3.  In total, 12,803 gene transcripts were commonly expressed at all stages analyzed; highly expressed common genes included a transport protein similar to SEC-2, the cell structure proteins actin and collagen (COL-1), and ribosomal and DNA-related highmobility proteins (Table 4).

Heatmap Expression and Functional Analysis
All FPKM values obtained through stage-specific analysis were subjected to two clustering processes based on the most abundant sequences and overall expression profiles ( Figure 4 and Figure S3, respectively) using a heatmap. The first clustering analysis included 60 highly expressed genes, with 10 genes specific to each stage and 10 genes commonly expressed at all stages. The second clustering method, based on overall gene expression profiles at all stages, revealed grouping of mRNA from the J3 and J4 stages, which formed a subclade similar to the egg-stage mRNA profile. Overall, mRNA from the female stage was distinct from other stages of nematode development.

Heatmap Expression and Functional Analysis
All FPKM values obtained through stage-specific analysis were subjected to two clustering processes based on the most abundant sequences and overall expression profiles ( Figure 4 and Figure S3, respectively) using a heatmap. The first clustering analysis included 60 highly expressed genes, with 10 genes specific to each stage and 10 genes commonly expressed at all stages. The second clustering method, based on overall gene expression profiles at all stages, revealed grouping of mRNA from the J3 and J4 stages, which formed a subclade similar to the egg-stage mRNA profile. Overall, mRNA from the female stage was distinct from other stages of nematode development.

miRNA Target Prediction and Regulation of mRNA Expression
To determine the targets of the miRNAs identified in our previous study, which reported 3107 total predicted miRNAs [3], the targets of these miRNAs were predicted from the 3 untranslated regions (UTRs) of M. incognita using in silico analysis. M. incognita 3 UTRs were more abundant across all ranges of sequence lengths, compared with the other nematode species C. elegans and Brugia malayi. However, for sequences 100-200 bp in length, sequences from the free-living nematode C. elegans were slightly more abundant. Furthermore, we performed sequence comparisons with mammals (human, mouse, and rat) and other organisms (fruit fly, zebrafish, and starlet sea anemone). Notably, except for sequences 0-100 bp in length, sequences from mammals and other groups were considerably more abundant at all sequence lengths ( Figure 5, Table S3). In M. incognita, a total number of 20,365 genes were predicted in the 3 UTR region and further, we selected 17,177 genes with N (any base) percentage less than 1 and sequences longer than 10 bp (Table S3).

miRNA Target Prediction and Regulation of mRNA Expression
To determine the targets of the miRNAs identified in our previous study, which reported 3107 total predicted miRNAs [3], the targets of these miRNAs were predicted from the 3′ untranslated regions (UTRs) of M. incognita using in silico analysis. M. incognita 3′ UTRs were more abundant across all ranges of sequence lengths, compared with the other nematode species C. elegans and Brugia malayi. However, for sequences 100-200 bp in length, sequences from the free-living nematode C. elegans were slightly more abundant. Furthermore, we performed sequence comparisons with mammals (human, mouse, and rat) and other organisms (fruit fly, zebrafish, and starlet sea anemone). Notably, except for sequences 0-100 bp in length, sequences from mammals and other groups were considerably more abundant at all sequence lengths ( Figure 5, Table S3). In M. incognita, a total number of 20,365 genes were predicted in the 3′ UTR region and further, we selected 17,177 genes with N (any base) percentage less than 1 and sequences longer than 10 bp (Table S3). Pooling of common miRNAs and their targets was conducted using three software programs: miRanda, RNAhybrid, and PITA. Our analysis revealed 2431 miRNAs that could potentially target 8331 mRNAs. Selected stage-specific miRNA-targeted mRNA transcript IDs and hit descriptions are listed in Table 5. Analysis of stage-specific miRNA targets revealed that heat shock, DNA repair, domain-containing proteins, and transcription factors were targeted in the egg stage, whereas cell cycle regulation, DNA replication, heat shock, and ubiquitin conjugation pathway-related proteins were targeted in the J2 stage. In the J3 stage, the major sperm protein (MSP) domain was identified among the highly targeted proteins. This protein has been linked to the motility of nematode sperm, as well as signaling functions [21]. Ras-related and myocardial muscle degradation indicator proteins were also targeted in the J3 stage. In the J4 stage, targets included translation-related proteins (RNA polymerase and small ribosomal subunit 40S). In the female stage, targets included molybdopterin-binding domain and sequences associated with ubiquitin-mediated protein degradation. Regulatory proteins such as adenyl cyclases were also targeted at the female stage. Pooling of common miRNAs and their targets was conducted using three software programs: miRanda, RNAhybrid, and PITA. Our analysis revealed 2431 miRNAs that could potentially target 8331 mRNAs. Selected stage-specific miRNA-targeted mRNA transcript IDs and hit descriptions are listed in Table 5. Analysis of stage-specific miRNA targets revealed that heat shock, DNA repair, domain-containing proteins, and transcription factors were targeted in the egg stage, whereas cell cycle regulation, DNA replication, heat shock, and ubiquitin conjugation pathway-related proteins were targeted in the J2 stage. In the J3 stage, the major sperm protein (MSP) domain was identified among the highly targeted proteins. This protein has been linked to the motility of nematode sperm, as well as signaling functions [21]. Ras-related and myocardial muscle degradation indicator proteins were also targeted in the J3 stage. In the J4 stage, targets included translationrelated proteins (RNA polymerase and small ribosomal subunit 40S). In the female stage, targets included molybdopterin-binding domain and sequences associated with ubiquitinmediated protein degradation. Regulatory proteins such as adenyl cyclases were also targeted at the female stage.
Further analysis was conducted to elucidate the interactions of miRNAs with mRNAs using the whole-genome sequence of M. incognita from the NCBI (ASM18041v1) and Worm-Base assembly (ASM18041v1a), which resulted in 1471 miRNAs (1218 known, 263 novel) targeting 7888 mRNAs. These miRNAs had a varying number of targets; one miRNA could target one to hundreds of mRNA. More than 40% of miRNAs targeted 2-10 mRNAs, while roughly 20 to 40% of miRNAs targeted 11-50 mRNAs among all stages ( Figure 6). Table 5. Stage-specific targets of miRNA obtained using miRanda, RNAhybrid and PITA. Prediction of miRNA targets was conducted only for miRNA sequences reported in our pervious study [3] with read count > 10.   Further analysis was conducted to elucidate the interactions of miRNAs with mRNAs using the whole-genome sequence of M. incognita from the NCBI (ASM18041v1) and WormBase assembly (ASM18041v1a), which resulted in 1471 miRNAs (1218 known, 263 novel) targeting 7888 mRNAs. These miRNAs had a varying number of targets; one miRNA could target one to hundreds of mRNA. More than 40% of miRNAs targeted 2-10 mRNAs, while roughly 20 to 40% of miRNAs targeted 11-50 mRNAs among all stages ( Figure 6). GO analysis indicated that the reference (PRJE28837) mRNA-targeted biological processes included cellular processes, while the molecular function mRNAs were involved in binding and catalytic activity. Among cellular components, mRNA targets included the organelle and cell classes ( Figure 7A-C). Similarly, common targets were found in the biological processes of single organism, metabolic, and cellular processes. Common molecular function mRNAs were involved in binding and catalytic activity. Among cellular components, membranes, organelles, and cells were most important ( Figure 7D-F). Biological processes of specific (known and novel) miRNA-targeted mRNAs were associated with metabolic processes. Additionally, molecular function mRNAs were involved in transporter, binding, and catalytic activities. Cellular component targets were membrane localized, with commonly expressed miRNAs and cellular mRNAs as targets of specific miRNAs ( Figure 7G-I). Overall, all miRNAs observed here targeted mRNAs localized to the cell, organelles, membranes, and macromolecular complexes, with roles in cellular and metabolic processes as well as catalytic, binding, transporter, and molecular transducer activities. At the molecular level, very few targets were related to the nucleic acid-binding transcription factor, molecular function regulator, protein-binding transcription factor, and antioxidant activity terms. At the cellular level, the membrane-enclosed lumen, extracellular region, cell junction, synapse, and virion were represented. In the biological process category, most targets were related to localization, developmental process, GO analysis indicated that the reference (PRJE28837) mRNA-targeted biological processes included cellular processes, while the molecular function mRNAs were involved in binding and catalytic activity. Among cellular components, mRNA targets included the organelle and cell classes ( Figure 7A-C). Similarly, common targets were found in the biological processes of single organism, metabolic, and cellular processes. Common molecular function mRNAs were involved in binding and catalytic activity. Among cellular components, membranes, organelles, and cells were most important ( Figure 7D-F). Biological processes of specific (known and novel) miRNA-targeted mRNAs were associated with metabolic processes. Additionally, molecular function mRNAs were involved in transporter, binding, and catalytic activities. Cellular component targets were membrane localized, with commonly expressed miRNAs and cellular mRNAs as targets of specific miRNAs ( Figure 7G-I). Overall, all miRNAs observed here targeted mRNAs localized to the cell, organelles, membranes, and macromolecular complexes, with roles in cellular and metabolic processes as well as catalytic, binding, transporter, and molecular transducer activities. At the molecular level, very few targets were related to the nucleic acid-binding transcription factor, molecular function regulator, protein-binding transcription factor, and antioxidant activity terms. At the cellular level, the membrane-enclosed lumen, extracellular region, cell junction, synapse, and virion were represented. In the biological process category, most targets were related to localization, developmental process, multicellular organismal process, response to stimulus, signaling, and cellular component organization or biogenesis (Figure 7G-I).

Stage-Wise Analysis of Highly Expressed mRNAs
To examine the abundance of target gene transcripts based on the stage-wise mRNA FPKM data resulting from high-throughput sequencing, 10 stage-specific highly expressed mRNAs were selected and validated through qRT-PCR. We used M. incognita β-actin (BE225475.1) primers as an internal control for normalization of gene expression. Upon experimental verification using qRT-PCR, we observed that all selected stage-wise mRNAs were expressed in M. incognita. Of the 10 stage-specific mRNAs, the egg stage included three mRNAs with unknown functions (Minc11396, Minc03738, and Minc10990) ( Figure 8A); the J2 stage included four mRNAs with unknown functions (Minc01592, Minc05278, Minc01593, and Minc18141) and two mRNAs with beta-1,4 endoglucanase (Minc19090) and cellulase (Minc09446) functions ( Figure 8B); one mRNA in the J3 stage possessed a double-stranded RNA-binding motif (Minc16891) and the other two (Minc00205 and Minc15344) had unknown functions ( Figure 8C); one mRNA in the J4 stage had an unknown function ( Figure 8D); and among seven mRNAs in the female stage, one mRNA encoded the secreted protein ASP-2 (Minc08697), two mRNAs were uncharacterized proteins (Minc18856 and Minc18984), and the remaining four mRNAs had unknown functions (Minc02270, Minc07322, Minc02526, and Minc02512) ( Figure 8E). Overall, the female stage had the highest level of mRNA expression (7), followed by the J2 (6), J3 (3), and egg (3) stages. In the J4 stage, we detected the fewest mRNAs (1) among the stage-specific mRNAs identified in our sequencing data. qRT-PCR results for the selected stage-specific mRNAs are shown in Figure 8. Furthermore, commonly expressed mRNAs were selected and validated; they were expressed at all stages ( Figure S4).

Stage-Wise Analysis of Highly Expressed mRNAs
To examine the abundance of target gene transcripts based on the stage-wise mRNA FPKM data resulting from high-throughput sequencing, 10 stage-specific highly expressed mRNAs were selected and validated through qRT-PCR. We used M. incognita βactin (BE225475.1) primers as an internal control for normalization of gene expression. Upon experimental verification using qRT-PCR, we observed that all selected stage-wise mRNAs were expressed in M. incognita. Of the 10 stage-specific mRNAs, the egg stage included three mRNAs with unknown functions (Minc11396, Minc03738, and Minc10990) ( Figure 8A); the J2 stage included four mRNAs with unknown functions (Minc01592, Minc05278, Minc01593, and Minc18141) and two mRNAs with beta-1,4 endoglucanase (Minc19090) and cellulase (Minc09446) functions ( Figure 8B); one mRNA in the J3 stage possessed a double-stranded RNA-binding motif (Minc16891) and the other two (Minc00205 and Minc15344) had unknown functions ( Figure 8C); one mRNA in the J4 stage had an unknown function ( Figure 8D); and among seven mRNAs in the female stage, one mRNA encoded the secreted protein ASP-2 (Minc08697), two mRNAs were uncharacterized proteins (Minc18856 and Minc18984), and the remaining four mRNAs had unknown functions (Minc02270, Minc07322, Minc02526, and Minc02512) ( Figure 8E). Overall, the female stage had the highest level of mRNA expression (7), followed by the J2 (6), J3 (3), and egg (3) stages. In the J4 stage, we detected the fewest mRNAs (1) among the stage-specific mRNAs identified in our sequencing data. qRT-PCR results for the selected stage-specific mRNAs are shown in Figure 8. Furthermore, commonly expressed mRNAs were selected and validated; they were expressed at all stages ( Figure S4).

Stage-Specific miRNA-Targeted mRNA Expression Analysis
The abundance of target gene transcripts was examined based on miRNA; mRNA results from high-throughput sequencing, represented as read counts, were normalized using the DeSeq2 and FPKM methods. We selected candidates based on inverse expression analysis of the DeSeq2 of upregulated miRNAs and corresponding FPKM of their downregulated target mRNAs. Table S7 is a validated list of miRNA ID, 5 -3 miRNA sequence, mRNA gene ID, genome contig location, and function for the candidates. We used M. incognita β-actin (BE225475.1) primers as an internal control for normalization of gene expression. Upon experimental verification using qRT-PCR, we observed that all selected miRNAs targeting mRNAs were irreversibly expressed at all stages of nematode development. Expression profiles of the 29 miRNAs targeting mRNAs were consistent when expression was compared with the internal standard β-actin (BE225475.1). Upon experimental verification using qRT-PCR, we observed that approximately 60% of the selected stage-wise miRNAtargeted mRNAs were expressed in M. incognita. The stage-specific miRNAs targeting mRNAs provided the following results: six miRNA-targeted mRNAs were detected in the egg stage with functions of glucosyltransferase (MIN00371-Minc13648) and patched family protein (MI00171-Minc18132) ( Figure 9A). In the J2 stage, the major functions were nematode cuticle collagen (MI03464-Minc01401) and beta-tubulin (MIN00380-Minc13117) ( Figure 9B). In the J3 and J4 stages, fucosyltransferase (MI02061-Minc00091), hypothetical protein (MI02049-Minc00288), ubiquitin-conjugating protein (MI03255-Minc00661), and unknown function (MI01808-Minc02324) ( Figure 9C,D) were the major functions. In the female stage, among five mRNAs, nematode cuticle collagen (MI01064-Minc12284) and putative secretion protein (MI02712-Minc11857) ( Figure 9E) were stage-specific. Our experimental validation showed that few miRNA-targeted mRNAs were expressed during other stages. These mRNAs had their lowest expression levels during the stages in which they were characterized as stage-specific according to high-throughput sequencing data. Selected qRT-PCR results for stage-specific miRNA-targeted mRNAs are shown in Figure 9.

Stage-Specific miRNA-Targeted mRNA Expression Analysis
The abundance of target gene transcripts was examined based on miRNA; mRNA results from high-throughput sequencing, represented as read counts, were normalized using the DeSeq2 and FPKM methods. We selected candidates based on inverse expression analysis of the DeSeq2 of upregulated miRNAs and corresponding FPKM of their downregulated target mRNAs. Table S7 is a validated list of miRNA ID, 5′-3′ miRNA sequence, mRNA gene ID, genome contig location, and function for the candidates. We used M. incognita β-actin (BE225475.1) primers as an internal control for normalization of gene expression. Upon experimental verification using qRT-PCR, we observed that all selected miRNAs targeting mRNAs were irreversibly expressed at all stages of nematode development. Expression profiles of the 29 miRNAs targeting mRNAs were consistent when expression was compared with the internal standard β-actin (BE225475.1). Upon experimental verification using qRT-PCR, we observed that approximately 60% of the selected stage-wise miRNA-targeted mRNAs were expressed in M. incognita. The stage-specific miRNAs targeting mRNAs provided the following results: six miRNA-targeted mRNAs were detected in the egg stage with functions of glucosyltransferase (MIN00371-Minc13648) and patched family protein (MI00171-Minc18132) ( Figure 9A). In the J2 stage, the major functions were nematode cuticle collagen (MI03464-Minc01401) and beta-tubulin (MIN00380-Minc13117) ( Figure 9B). In the J3 and J4 stages, fucosyltransferase (MI02061-Minc00091), hypothetical protein (MI02049-Minc00288), ubiquitin-conjugating protein (MI03255-Minc00661), and unknown function (MI01808-Minc02324) ( Figure 9C, D) were the major functions. In the female stage, among five mRNAs, nematode cuticle collagen (MI01064-Minc12284) and putative secretion protein (MI02712-Minc11857) (Figure 9E) were stage-specific. Our experimental validation showed that few miRNA-targeted mRNAs were expressed during other stages. These mRNAs had their lowest expression levels during the stages in which they were characterized as stage-specific according to high-throughput sequencing data. Selected qRT-PCR results for stage-specific miRNA-targeted mRNAs are shown in Figure 9.

Discussion
M. incognita is an economically important pathogen and one of the most important plant pathogens worldwide. Previous studies regarding the transcriptome of M. incognita were limited to juvenile stages J2 and J3 [13,14]. These studies indicated an abundance of cytoskeletal and ligand-binding proteins during the J2 stage, as well as high expression of host immune response genes and parasitic functional glutathione-S-transferases in the J3 stage. In the current study, we conducted RNA sequencing of M. incognita in all developmental stages. In general, obligate parasites have compact genomes, and members of the genus Meloidogyne reportedly have the smallest genomes among plant parasites [8,22,23]. Compaction has been proposed as the reason for reduced genome size in organisms that have small genomes with high gene density. The previously reported draft genome studies of nematode revealed 16,419, 18,074, 14,420, 18,348, and 19,212 genes from M. incognita, G. pallida, B. xylophilus, M. hapla, and B. malayi, respectively [2,8,[24][25][26]. Overall, 17,423 transcripts were commonly expressed across all stages of M. incognita in our study (Figure 2). As a result, the transcripts obtained in our results are also close to the values in other nematodes.
Alignment of the sequencing results to available databases indicated that our sequences were orthologous to the sequences of B. malayi (3729 hits), C. elegans (2489 hits), and C. briggsae (1871 hits), comprising 46.43% of transcripts (Table S2). The filarial nematode B. malayi also exhibits a parasitic lifestyle with infective and reproductive stages, which may explain the similarity in expressed genes. Annotation of the mRNA sequences indicated that a large number of mRNA sequences had no matches in the database search ( Figure S2). The same phenomenon was observed in previous transcriptional studies of M. incognita and other nematodes [8][9][10][11]13]. Elevated expression of actin and the high-mobility group of mRNAs in the egg stage suggest that cell proliferation and DNA-associated processes were active. Annotation of the RNA sequences also indicated the presence of the collagen family protein COL-45 and alpha-(1,3)-fucosyltransferase transcripts; the alpha-(1,3)-fucosyltransferase in M. incognita is similar to fut-1 in C. elegans, which is expressed in the egg stage. The J2 stage had abundant FMRF-amide-like peptides, which play key roles in neuroregulatory functions. Cellulases similar to the cellulases in Meloidogyne javanica (β-1,4-endoglucanase) and Heterodera glycines (guanylate cyclase) were also expressed at the J2 stage. During infective stages of the nematode, cellulases can aid the infection of plant tissues, and several cellulases and endoglucanases have been reported in the stylet secretions of the cyst nematode Globodera rostochiensis [27,28]. Therefore, motility and infection may be the crucial functional roles of genes expressed at the J2 stage. The J3 and J4 stages exhibited large numbers of collagen and actin-related mRNAs, which are involved in nematode growth and metamorphosis. The cytochrome P450 and cuticle collagen family proteins CYP-23A1 and COL-149, as well as poly(A) polymerases orthologous to the polymerases in C. briggsae, were expressed in the J3 stage, whereas several zinc finger proteins and zinc metalloproteinases were specific to the J4 stage. The female-stage nematode also possessed large numbers of cuticle collagen mRNAs, SEC-2-like proteins (which control vesicular transport), and nucleic acid-binding KH domains ( Figure 5, Table 4). Differential expression patterns of miRNAs and their effects on mRNA levels showed that several miRNAs were expressed at specific stages of nematode development; these miRNAs may regulate mRNA expression. We observed 1218 known miRNAs and 253 novel miRNAs that targeted mRNAs in M. incognita. By comparing the expression levels of miRNAs and mRNAs, we found 76 known and 4 novel miRNAs regulating more than fourfold changes in the expression levels of their target mRNAs.
Among the targeted mRNAs with fourfold changes, the egg and female stages had the greatest number and diversity of targets compared with all other stages. In the egg stage, the greatest numbers of targets were unknown or putative uncharacterized proteins from C. elegans and C. briggsae. Known targets in the egg stage included carboxypeptidase, a regulator of the nonsense transcripts homolog; β-1,4-N-acetylgalactosaminyltransferase, a metabolite transporter protein homolog; putative esophageal gland cell secretory, helicase-like, and cytochrome P450 family 33-related proteins; and ecdysteroid UDP-glucosyltransferase from C. briggsae. In the J2 stage, highly targeted known proteins included β-tubulin, histonelysine N-methyl transferase, collagen, ubiquitin, and parasitism-related Vap-1 homolog proteins. The J3 and J4 stages had fewer targets that showed significant changes in expression (fourfold). Fucosyltransferase family proteins and guanylate cyclase were major targets in the J3 stage, whereas Nk homeobox protein was the major target in the J4 stage. In the female stage, most target proteins were cuticle-and collagen-related proteins; other targets included the PAN domain, ion channel receptor, HAT family dimerization domain, leucine-rich repeats, RNase L inhibitors, and TRP homologous cation channel protein.
Quantitative analysis of expression using PCR validated that stage-specific miRNAs targeting mRNAs identified in the sequencing-based analysis were expressed in M. incognita. We confirmed that these miRNAs targeting mRNAs were authentic; 29 miRNAs targeting mRNAs in the egg, J2, J3, J4, and female stages exhibited stage-specific expression patterns in sequencing reads and corresponding patterns in quantitative PCR analysis. Moreover, selected mRNAs-identified as stage-specific through digital expression analysis based on RNA sequencing-were expressed at all stages. However, the functions of these miRNAs that target mRNAs in M. incognita remain unclear. In summary, the present study aimed to supplement existing miRNA data regarding M. incognita to identify targeted mRNAs at all stages of development. The data presented here will be useful for future nematologists when selecting potential regulatory networks for targeted control of nematode-induced agricultural losses.

Biological Sample Preparation
M. incognita was collected from roots of an infected tomato plant (Solanum lycopersicum var. Rutgers) grown in our laboratory greenhouse at a constant 25 • C. The method used to collect each stage of nematode samples was described in our previous study [20]. In brief, eggs were washed from the infected roots and treated with 10% NaClO for 5 min, then excess water was passed through a 25 µm mesh to collect the eggs, which were purified through 35% sucrose gradient centrifugation at room temperature. J2 samples were collected by hatching eggs at 25 • C for 5 days in autoclaved distilled water, and samples were collected using 5-7 layers of Kimtech Science Wipers on a Petri dish placed on a laboratory bench. To collect the J3, J4 and female stages, infected roots were washed, chopped, and treated with 7.7% cellulase and 15.4% pectinase, then washed and filtered through a 75-µm filter. The samples stored on the filters were re-rinsed in water, and nematodes were hand-picked using a pipette under a stereomicroscope.

Stage-Wise RNA Extraction and Sequencing
Nematode samples were ground with a mortar and pestle in liquid nitrogen; an approximately 200 mg sample was used for RNA extraction for each of the five developmental stages (i.e., egg, J2, J3, J4 and female) with three biological replicates collected. Total RNA from the samples was quantified using a Nanodrop spectrophotometer (Thermo Scientific, CA, USA); quality was assessed with the RNA 6000 Nano Kit (Agilent Technolgies, Santa Clara, CA, USA) and Bioanalyzer 2100 (Agilent Technolgies, Santa Clara, CA, USA). One microgram of treated RNA was used for next-generation sequencing; libraries were generated using the TruSeq RNA Sample Prep Kit (Illumina Inc., San Diego, CA, USA), in accordance with the manufacturer's instructions. In brief, poly(A)-containing RNA molecules were extracted using a poly-T oligonucleotide attached to magnetic beads. The purified total poly(A) RNA was fragmented into small pieces using divalent cations at elevated temperature. The cleaved mRNA fragments from each stage were treated with DNase I at 37 • C for 30 min to remove residual DNA. Using random primers, the firststrand cDNA was synthesized and the fragments were purified with the QIAquick PCR extraction kit. Finally, elution buffer was added for end repair and addition of poly(A). Subsequently, the repaired short fragments were attached to sequencing adapters, and then, each library was assigned an adjoining distinct multiplex identifier tag. The resulting cDNA libraries were sequenced with the Illumina HiSeq™ 2500 system.

Assembly and Genome Mapping
Complete paired-end sequences were obtained as individual fasta files (forward and reverse) from the images using CASAVA v.1.8.2 base-calling software with an ASCII Q-score offset of 33. Raw reads sequences were trimmed by the parameters, quality trimming based on Phres quality scores (Q ≤ 20), adaptor trimming and minimum length discard (<90 bp) using the CLC genome cell (v. 4.0). To characterize the quantitative expression patterns of individual genes, the clean sequence reads from five libraries with three replicates (egg, J2, J3, J4 and female) were mapped individually to the reference genome under the guidance of gff3 using Bowtie from TopHat (v. 2.0.8) with the following parameters: mate-inner-dist (100), mate-std-dev (200), splice-mismatch (1), library-type (fr-unstranded), and microexon-search. Differential expression between pairs was calculated using Cuffdiff (v. 2.1.1) [29]. Final sequences were filtered with the following criteria: FPKM ≥ 0.3 [30], FDR (q value) < 0.01 and log 2 (FC) ≥ 1.00.

Functional Annotation and Expression Analysis
Each gene sequence was subjected to functional annotation using Pedant-Pro suite (Biomax Informatics AG, Martinsried, Germany) [31] with customized database parameters: eukaryote non-plant database; analysis type: Intronic rich genes (nematodes, C. elegans); genetic code: stranded, blastp, and Blast2GO [32] and other default parameters. Functional annotations such as GO terms and Clusters of Orthologous Groups were assigned to the gene sequences using Pedant-Pro default scores. The expression patterns of genes across different stages were compared using a heatmap generated with MeV software (v. 4.8.1).

miRNA Sequencing and Target Prediction
miRNA sequencing and classification at various stages of the M. incognita lifecycle were described in detail in our previous report [3]. In total, 3107 miRNAs were used for target prediction. To obtain the 3 UTR region of each gene, prediction was conducted using Python scripts as follows. Initially, gene coordinates were taken from the GFF3 file for M. incognita available from WormBase (ASM18041v1a) [2], and the 3 flanking region from the stop codon to a maximum of 5000 bases downstream was searched for 16 polyadenylated site motifs. The region from the stop codon to the first polyadenylated site motif plus 20 downstream bases was considered the 3 UTR for each gene. The polyadenylated site motifs were taken from the C. elegans genome [33].
To predict potential targets, we used bioinformatics software including miRanda, RNAhybrid, and PITA. First, we used two different methods, miRanda [34] and RNAhybrid [35,36], to predict miRNA targets by detecting the most energetically favorable hybridization sites for small RNA within the 3 UTR based on sequence similarity features [37]. The default parameters for miRanda were used and an output filter was applied (∆G ≥ −25.0 kcal/mol). The parameters for RNAhybrid were number of hits per target = 1, energy cutoff ≥ −25.0 kcal/mol, and maximal internal or bulge loop size per side = 4. Finally, the common mRNA targets obtained from miRanda and RNAhybrid were used for prediction with PITA, where ∆∆G values below −10 were considered representative of endogenous miRNA expression levels.

Comparison of miRNA and Complementary mRNA Expression
miRNA and mRNA expression patterns were integrated with targets through inverse expression analysis (miRNA upregulated and corresponding target mRNA downregulated) using Python scripts. miRNAs were normalized to obtain DESeq2 and mRNA expression was normalized to obtain FPKM. The miRNA and target expression levels were then compared.

Stage-Wise Gene Expression Analysis
qRT-PCR was performed for 29 miRNAs targeting mRNAs and 20 stage-specific mRNA genes selected based on stage specificity, high expression levels determined using