1,4-NTQ (purity ≥96.5% w/w, Sigma-Aldrich) was dissolved in an aqueous solution of Tween^® 80 (suitable for cell culture, Sigma-Aldrich) at 100 ppm, to get a final 1,4-NTQ concentration of 20 ppm. This concentration has a sub-lethal effect on M. luci J2, as previously demonstrated by Maleita et al. (2022). Approximately 18,000 J2 were incubated at 22°C for 3 days in 20 ppm 1,4-NTQ (N treatment). Water (H treatment) and Tween 80^® 100 ppm (T treatment) were included as controls. Each of the three treatments was replicated three times. After N and T treatments, J2 were centrifuged to remove 1,4-NTQ and Tween^® 80, washed three times with RNase-free water, and stored at −80°C until RNA extraction. Nematodes exposed to H treatment were also subjected to the same procedure.

Nematode RNA was extracted from the nine treated samples by Trizol/chloroform (Invitrogen) extraction and DirectZol clean up (Zymo-research) and treated with the TURBO^™ DNase (Invitrogen) to remove the remaining DNA, according to the manufacturer’s instructions. Total RNA quality was assessed by RIN (RNA Integrity Number) determination with the Agilent RNA 6000 Pico Assay Kit on the Bioanalyzer 2100 (Agilent Technologies) and quantified by Qubit RNA HS Assay kit on the Qubit 2.0 Fluorometer (Thermo Fisher Scientific). Extracted RNA was stored at −80°C until sequencing or RT-qPCR.

A poly(A) enriched strand-specific library was generated with the TruSeq Stranded mRNA sample preparation kit (Illumina) from approximately 0.5 µg of high-quality total RNA of each sample. Transcriptome sequencing was performed on a NextSeq 550 Illumina sequencer with the NextSeq 500/550 High Output v2.5 Reagent Kit (150 cycles), according to the manufacturer’s instructions (Illumina) at Genoinseq (Cantanhede, Portugal).

Sequence data were processed at Genoinseq. Raw reads were extracted from the Illumina Nextseq^® System in fastq format with bcl2fasta version 2.20.0.422 (Illumina) and quality-filtered with fastp version 0.20.0 (Chen et al., 2018) to trim bases with an average quality lower than Q25 in a window of 5 bases and to remove polyadenylated tails above 10 bases and reads with less than 50 bases. The ribosomal RNA content was estimated with SortMeRNA version 2.1 (Chen et al., 2018). High-quality reads were de novo assembled into transcripts using Trinity version 2.2.0 (Grabherr et al., 2011). Reads from each sample were aligned to the de novo transcriptome assembly with Trinity’s Bowtie module.

Transcripts were translated into amino acid sequences with Transdecoder version 5.5.0 (TransDecoder, 2021), with parameters set to a minimum length of 50 amino acids and no strand specificity. Translated sequences were then annotated using Diamond version 2.0.8 (Buchfink et al., 2015) against SwissProt and RefSeq with an e-value of 0.001.

Transcript abundance and gene count matrix were estimated using Trinity’s Transcript Quantification module with RSEM (Li and Dewey, 2011). Pre-processing and exploratory data analysis (Heatmap, K-means and PCA) of gene-level read count matrix were performed in IDEP.95 (Ge et al., 2018) with default parameters. Differentially expressed genes (DEGs) between conditions were identified using the DESeq2 package in iDEP.95, considering a fold change |log₂(FC)| ≥ 1 with a false discovery rate (FDR) ≤ 0.05 (adjusted p-value).

The gene ontology (GO) annotations of translated transcripts corresponding to DEGs were achieved with OmicsBox (2019) based on the BLAST against the non-redundant protein database NCBI and InterPro database using the default settings in each step (Götz et al., 2008). The GO analysis is done in three categories: (i) molecular function, which defines molecular activities of gene products; (ii) cellular component, which describes where gene products are active; and (iii) biological process, which clarifies the pathways and larger processes made up of the activities of multiple gene products. GO enrichment analysis was conducted for the selected DEGs in a condition against the total of DEGs among the three conditions, in OmicsBox with the statistical Fisher’s Exact Test, using a p-value of 0.05 as cutoff. Functional annotation of selected DEGs was complemented with analysis on the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database (Kanehisa and Goto, 2000) in OmicsBox.

From the non-overlapping DEGs found in N–H comparison, four genes were selected from the list of enriched GO terms or most representative KEGGs to validate relative quantitative data obtained from sequencing data. The relative transcript abundance of selected genes was assessed by reverse transcription quantitative real-time PCR (RT-qPCR) with SsoAdvance Universal SYBRGreen supermix (Bio-Rad), according to standard protocols, using the CFX96 Touch™ Real-Time PCR Detection System (Bio-Rad). Extracted RNA was converted into cDNA using the iScript cDNA Synthesis kit (Bio-Rad) and used in qPCR. The amplification kinetics of each transcript was normalized with the amplification kinetics of the malate dehydrogenase (mdh) and ribosomal protein S6 (rps6) genes, chosen as endogenous controls according to Wu et al. (2019) and amplified with primers adapted to M. luci corresponding transcript sequences obtained in this study. Other primers used in qPCR were also designed based on M. luci transcript sequences ( Table S1 ). qPCRs were done at 98°C for 30 s, followed by 40 cycles of 98°C for 5 s and 60°C for 30 s. Melting curve analyses were performed and validation experiments were first carried out to ensure equivalent amplification efficiency for all transcripts. The RT-qPCRs were conducted for the three repetitions of each treatment, with three technical replicates for each qPCR. Amplification efficiencies and Ct values were determined by the CFX ManagerTM Software 2.1 (Bio-Rad) and the mean Ct values used in the REST software (Pfaffl et al., 2002) for relative transcript level and statistically significant differences analysis using the Pairwise Fixed Reallocation Randomization Test^©.

Sequencing of the nine M. luci RNA samples produced 722,024,674 high-quality reads, consisting approximately of a total of 55.5 Gbases. The sample library had an average size of 82,298,380 paired reads and the percentage of high-quality reads after quality filtering was >97% ( Table S2 ). From this, de novo assembly generated 58,042 transcripts corresponding to the de novo transcriptome of M. luci J2 available through the Short Read Archive (SRA) under the BioProject accession number PRJNA940699, which constitutes the first transcriptomic data available for this species.

From the 58,042 transcripts, 46,271 were predicted to be translated into amino acid sequences, and from these, 19,393 (33.4%) and 11,663 (20.1%) were annotated against the Swiss-Prot and RefSeq databases, respectively. The low number of annotated transcripts is in accordance with what is found in other PPN species with a significant number of candidate genes lacking annotation and a predicted function (Vieira et al., 2015; Petitot et al., 2016). Compared to genome data on other Meloidogyne species, the number of coding transcripts is close to that described for M. incognita BioProject PRJEB8714, for which a total of 43,718 coding genes were predicted (Blanc-Mathieu et al., 2017).

Considering that some of the transcripts correspond to isoforms of the same gene, a gene abundance matrix with a total of 47,435 expected counted genes was obtained, and from these, after IDEP.95 default filtering (at least 0.5 counts per million in at least one library), 28,165 genes were left. Pre-process analysis revealed homogeneity among the nine libraries, with small variation in library sizes ( Figure 1A ), small variation of transformed data between replicates ( Figure 1B ), and similar distribution of the transformed data ( Figures 1C, D ).

The top 1,000 most variable genes were ranked by their standard deviation across all samples in hierarchical clustering and consistent clusters were obtained in Heatmap and K-means analysis ( Figures 2A, B ). A PCA plot using the first and second principal components revealed a clear difference between the water (H) control treatment, 1,4-NTQ (N) treatment, and Tween 80^® (T) treatment ( Figure 2C ). Variations among replicates were minimal ( Figure 2 ).

From the 28,165 genes, 7,854 were found differentially expressed among the three conditions. A higher number of DEGs was found between N–H treatments (5,744) than between T–H (4,550) or T–N (3,060) treatments. As expected, a lower number of DEGs was found between T–N treatments, since Tween^® 80 (T treatment) was used as the 1,4-NTQ solvent (N treatment). The number of downregulated genes was much higher in both N and T treatments compared to H treatment than the number of upregulated genes in the same comparisons. In the T–N comparison, the differences in the number of up- or downregulated genes were the lowest ( Figure 3 ).

Overall, the M. luci gene expression was highly influenced by 1,4-NTQ treatment with a higher number of DEGs (5,744) than those reported for M. incognita treated with fluensulfone (1,831), fluazaindolizine (3,256), oxamyl (334), and fluopyram (166) (Wram et al., 2022b) or in M. incognita treated with P. tenuis huperzine A (1,079) or fungal filtrate (1,635) (Kassam et al., 2023).

From the 1,338 non-overlapping DEGs found in the N–H comparison ( Figure 3A ), 250 were found to be upregulated and 1,088 were found to be downregulated ( Table S3 ). These genes were further explored as their regulation could be associated with a more specific effect of the 1,4-NTQ.

In order to find which group of genes are overrepresented in the 1,338 DEGs unique in the N–H comparison, a GO enrichment analysis was done for the 250 upregulated and 1,088 downregulated genes against all DEGs. Among the upregulated genes, these analyses revealed a strong enrichment on ABC-type transporter activity, binding, and enzyme regulator activities, belonging to GO molecular function terms, and regulation of biological process and biological regulation, belonging to GO biological process terms ( Table 1 ).

From the 250 upregulated genes, 82 transcript sequences were associated with various KEGG pathways ( Table S4 ) with a higher number of transcript sequences assigned to the thyroid hormone signaling pathway and insulin signaling pathway in the Metabolism category; ABC transporters, sphingolipid signaling, AMPK signaling, PI3K-Akt signaling, Ras signaling, and MAPK signaling pathways in the Environmental Information Processing category; vitamin digestion and absorption and longevity regulating pathway-worm in the Organismal Systems category; and proteasome in the Genetic Information Processing category. A higher number of transcript sequences of upregulated genes were found to be associated with KEGG pathways belonging to the Environmental Information Processing category, namely, ABC transporters associated with nematode xenobiotic detoxification ( Figure 4 ).

On the other hand, among the downregulated genes, an enrichment of GO terms associated with molecular transducer activity and adenylyltransferase and molybdopterin molybdotransferase activities, belonging to the GO molecular function terms, and associated with the molybdopterin cofactor metabolic process and prosthetic group metabolic process, belonging to GO biological process terms, was found ( Table 2 ). Molybdopterin cofactor is essential for the catalytic activity of some enzymes such as sulfite oxidase and xanthine dehydrogenase, whose failure is associated with severe neurological abnormalities (Duran et al., 1980), and aldehyde oxidase, which plays an important role in the metabolism of several drugs due to its involvement in cellular redox stress (Garattini and Terao, 2012).

From KEGG analysis of 1,088 downregulated genes, 316 transcript sequences were associated with various KEGG pathways ( Table S5 ), and from these, a higher number of transcript sequences were associated with ribosome and spliceosome in the Genetic Information Processing KEGG category, related to translation and transcription processes, respectively ( Figure 5 ). The considerable higher number of downregulated transcripts associated with ribosome reveals a clear inhibitory effect of 1,4-NTQ on translational processes, compromising protein synthesis and consequently several core biological processes. Ribosomes are attractive targets in the antitumor, antiviral, antibacterial, antifungal, and antiparasitic therapies, with several compounds suppressing protein biosynthesis in eukaryotic cells, such as antibiotics, being widely used in medicine and also in agriculture (Lin et al., 2018; Mann et al., 2021). Downregulation of a high number of genes related to ribosome pathway suggests ribosome as one of the main targets of 1,4-NTQ in M. luci.

The biosynthesis of cofactors pathway was also highly represented in downregulated genes’ KEGG analysis and is in accordance with molybdopterin cofactor-related activities identified in GO enrichment analysis. Additionally, the thermogenesis pathway, implicated in stress response and ensuring the normal cellular and physiological function under conditions of environmental challenge, was also compromised by 1,4-NTQ, with a high number of downregulated genes associated with this pathway. The neurotoxicity of 1,4-NTQ was also evident with the highly represented GABAergic synapse pathway among the downregulated genes.

Moreover, sphingolipid metabolism is also highly represented among the downregulated genes; however, the related sphingolipid signaling pathway was also identified among the most represented KEGGs in upregulated genes. Sphingolipids are bioactive lipid molecules found in the membranes of all eukaryotic cells and have critical functions in the control of cell growth, senescence, differentiation, and programmed cell death. Different intermediates of sphingolipid pathways can have opposing effects on cell signaling (Haughey, 2010), and therefore, imbalances in sphingolipid metabolism and signaling, caused by 1,4-NTQ, can deregulate key cellular processes in M. luci.

Furthermore, in common with most represented KEGGs in upregulated and downregulated genes was the mitogen-activated protein kinase (MAPK) signaling pathway. The MAPK signaling pathway is known to serve as a transducer of extracellular stimuli that allow cellular adaptation to changes in the environment, playing important roles in complex cellular processes like proliferation, differentiation, development, transformation, and apoptosis (Andrusiak and Jin, 2016), and consequently, deregulation of this pathway is another key adverse effect of 1,4-NTQ in M. luci.

On the other hand, the most represented KEGG pathway in upregulated genes, ABC transporters, was not identified in downregulated genes. Moreover, the most represented KEGG in downregulated genes, ribosome, was also not identified among the upregulated genes, reflecting a consistent stimulant and inhibitory effect of 1,4-NTQ on ABC transporters and ribosome pathways, respectively.

Overall, it is perceptible that the antagonist effect of 1,4-NTQ on M. luci transcriptome disturbs several fundamental cellular processes. These results support the findings in previous bioassays, where 1,4-NTQ showed nematicidal effects on M. luci mortality, hatching, root penetration, and reproduction in tomato (Maleita et al., 2022).

Nematodes use detoxification pathways to protect themselves from the effect of toxic compounds. Nematode xenobiotic metabolism includes three phases. In phase I, functional groups such as hydroxyl groups are added to xenobiotics by enzymes mainly from cytochrome P450 (CYP) and short-chain dehydrogenase/reductase (SDR) families, to increase their solubility. In a second phase, enzymes mainly from UDP-glucuronosyl transferase (UGT) and glutathione S-transferase (GST) classes catalyze the conjugation of xenobiotics with polar molecules, promoting their water solubility. Finally, the excretion of xenobiotics (phase III) is made by ATP-binding cassette (ABC) transporters and other transmembrane transporters that actively export catalyzed compounds across the cytoplasmic membrane (Hartman et al., 2021).

Transcript levels of enzymes involved in detoxification were clearly affected in M. luci J2 exposed to 1,4-NTQ ( Table 3 ). From the DEGs belonging to the main enzyme families involved in phase I of detoxification, two CYP were found to be downregulated and one SDR was found to be upregulated. Moreover, from the seven DEGs associated with the metabolism of xenobiotics by cytochrome P450 (ko00980), two were found to be upregulated (one glycosyl hydrolase and one UGT) and five were found to be downregulated (two UGTs, one GST, and two dehydrogenases), compromising phase I and phase II of the detoxification process and possibly compromising the ability to encounter this external toxic compound. In addition, and as discussed above, from GO enrichment analysis, the activity of several enzymes that could be involved in this detoxification process, such as aldehyde oxidase, could be affected not at the transcript level but by their cofactor’s availability. On the other hand, phase III of the detoxification process was stimulated with the upregulation of all transcripts related to ABC transporters and the activation of the final excretion of metabolized toxic compounds. The same tendency was found in the pinewood nematode (PWN), Bursaphelenchus xylophilus, in response to emamectin benzoate, used to control PWN infection by pine tree trunk injection, with the downregulation of several UGTs and some of the differentially expressed GSTs and the upregulation of ABC transporters (Lu et al., 2020). Also in B. xylophilus, the response to the tree host-derived α- or β-pinene in reaction to pathogen attack resulted in the downregulation of GSTs and the upregulation of ABC transporters, but UGTs were mainly upregulated (Li et al., 2019; Li et al., 2020). In M. incognita, transcripts related to xenobiotic detoxification were also revealed to be affected by several nematicides; however, the pattern of up- and downregulated expression varied among the several tested nematicides (Wram et al., 2022b).

The relative transcript level of four genes, selected from the list of non-overlapping DEGs in the N–H comparison, was determined by RT-qPCR. From these, TRINITY_DN20723_c2_g1 and TRINITY_DN21456_c0_g1 were found to be significantly upregulated in the N-treated condition compared to the H-treated condition, whereas TRINITY_DN9569_c0_g1 and TRINITY_DN8190_c0_g1 were found to be significantly downregulated ( Figure 6 ), in accordance with relative quantitative data of those genes obtained from RNA-seq ( Table S6 ). Additionally, the pattern of changes in transcript levels of these four genes, among the other two pairwise comparisons (T–H and N-T), was consistent using both qPCR and RNA-seq approaches ( Figure 6 , Table S6 ).