Skip to main content
Download PDF

Time-course RNA-seq analysis of upland cotton (Gossypium hirsutum L.) responses to Southern root-knot nematode (Meloidogyne incognita) during compatible and incompatible interactions

BMC Genomics volume 26, Article number: 183 (2025) Cite this article

Abstract

Background

The Southern root-knot nematode (Meloidogyne incognita) poses a substantial threat to cotton (Gossypium hirsutum L.) by causing significant agricultural losses. Host plant resistance is the most plausible approach for minimizing these losses. QTL mapping and early transcriptomic studies have identified candidate genes within the QTL regions on chromosome 11 (qMi-C11) and chromosome 14 (qMi-C14). Although these QTL regions have been fine-mapped and candidate genes identified, expression profiling of Meloidogyne-Gossypium interactions across different stages of infection could further refine the list of candidate genes. This study advances our understanding of the molecular mechanisms underlying the resistance conferred by qMi-C11 and qMi-C14 against Southern root-knot nematode.

Results

Using time-course RNA-seq analyses across nematode developmental phases, we uncovered transcriptomic events—both genome-wide and within QTL intervals—underlying defense responses during compatible interactions (with Cocker 201, a susceptible line) and incompatible interactions (with M-120 RNR, a resistant line). Basal defense responses were observed in both compatible and incompatible interactions, with stronger expression in the incompatible interaction. Nematode-responsive genes associated with defense pathways showed distinct dynamics, characterized by repression during compatible interactions and early induction, greater diversity, and heightened upregulation during incompatible interactions. This study uncovers a broad repertoire of disease resistance and putative resistance genes, as well as pathogenesis-related genes, ligands, and receptors, that are differentially expressed in response to nematode parasitism. Mapping of these genes across the cotton genome identified promising candidates, including Gh_A11G3090 (PUB21) and Gh_A11G2836 (RPPL1) within the chromosome 11 QTL region, andGh_D02G0257 (RLP12) and Gh_D02G0259 (RLP12) within the chromosome 14 QTL region.

Conclusions

The findings of this study deepen our understanding of host-nematode interactions, identify candidate genes for downstream applications, and contribute to advancements in resistance breeding and sustainable nematode management strategies.

Peer Review reports

Background

The Southern root-knot nematode (RKN), Meloidogyne incognita, is one of the most important agricultural pests of upland cotton [Gossypium hirsutum L.; 2n = 4x = 52 (subgenomes ‘A’ and ‘D’ with 26 chromosomes each) due to its widespread distribution, genetic diversity, broad host range, and polyphagous feeding habits [1, 2]. In the United States alone, M. incognita has caused losses exceeding 65 million kilograms of cotton, valued at over 123 million dollars [3, 4]. Recommended management practices include cultural controls, such as crop rotation and cover crops, biological control with nematopathogenic bioagents, and chemical control through nematicide applications [5]. However, the broad host range of M. incognita complicates cultural and biological management methods, leaving growers with limited nonhost crop options for rotation. Nematicides, while effective for early-season nematode suppression, significantly increase production costs.

Host plant resistance is the most economical, practical, and environmentally sound method for managing this subterranean pest. Currently, the resistant line ‘Auburn 623 RNR’ and the elite breeding ‘M’ lines derived from it remain the primary source of RKN resistance. Histopathological comparisons of compatible and incompatible interactions—resulting in successful infection or defensive response respectively—using M-lines or their derivatives have shown that resistance operates through a two-stage post-penetration interference. The first, or ‘early’ mechanism occurs soon after infection (8–12 days after inoculation) to prevent juvenile nematodes from developing functional feeding sites. The second, or ‘late’ mechanism, occurs later in the infection process (25–30 days after inoculation) to impede the development of nematodes into egg-laying adult females [6,7,8,9,10]. These findings were supported by molecular genetic studies, which identified two QTLs, qMi-C11 and qMi-C14, that confer RKN resistance. The qMi-C11 region suppresses root gall formation, while qMi-C14 reduces egg production with minimal impact on galling [11,12,13,14,15]. Resistance conferred by qMi-C14 has been attributed to the gene Gh_D02G0276, which encodes for hAT-like transposases with novel N- and C-terminal domains that resemble targets of known RKN effector molecules [16]. However, the causal gene underlying qMi-C11 remains unidentified.

Comparative transcriptomic analysis using time-course RNA sequencing (RNA-seq) data enables researchers to pinpoint candidate genes within QTL regions and uncover genome-wide, spatiotemporal expression dynamics of genes involved in complex regulatory and biochemical pathways. This approach has been successfully applied to identify nematode-responsive genes and pathways involved in host-nematode interactions. Time-course RNA-seq has been used to elucidate RKN-host interactions in plant families such as Fabaceae [17], Rosaceae [18], Cucurbitaceae [19, 20], and Solanaceae [21, 22]. Comparative genomics has also been used as a predictive tool to identify conserved plant genes responsive to RKN infections [23]. The resistance mechanisms observed in the RKN-resistant lines derived from Auburn 623 RNR are similar to those found in Cucumis metuliferus [20], Vigna unguiculata [24], and Medicago truncatula [25], where resistance does not prevent nematode penetration and does not involve a hypersensitive response (HR). Comparative host-RKN transcriptomics, using data from diverse plant species exposed to RKN, offers great potential for identifying candidate genes involved in host-RKN interactions. This approach is facilitated by advancements in whole genome sequencing, which have improved our understanding of synteny and collinearity among plant genomes.

In an exploratory RNA-seq study, numerous differentially expressed genes were identified in a susceptible cotton line (compatible interactions) and its near-isogenic line with nematode resistance (incompatible interactions) 12 and 30 days after inoculation [26]. Two key differences in gene regulation were observed between the early and late stages of infection. First, the resistant genotype induced the expression of genes to resist invasion and restrict the movement of RKN within the roots, whereas the susceptible genotype expressed genes that facilitated the establishment of feeding sites and the development of female RKN. Second, while the QTL loci may contribute to RKN resistance, the major contribution likely stems from the elevated expression of defense response genes common to both genotypes. An RNA-seq study conducted 10 days-after-inoculation supports this observation, revealing that defense response genes regulated by hormones such as jasmonic acid (JA) and salicylic acid (SA) are constitutively expressed but accumulate at higher levels in the roots of the resistant genotype ‘Acala NemX’ [27]. Thus, key genetic switches controlling resistance need to be investigated by examining gene expression across various stages of infection, from the establishment of feeding sites by the M. incognita second-stage juveniles (J2s) to their postinfection development into adult females.

We conducted a comprehensive time-course RNA-seq analysis using the resistant and susceptible genotypes and experimental conditions reported by Kumar et al. [26], tracking the entire life cycle of M. incognita development in cotton roots [26]. Comparative transcriptomic analysis was performed across five developmental phases (4, 8, 12, 16 and 20 days after inoculation, DAI), providing a clearer understanding of the transcriptomic changes and molecular mechanisms underlying the two modes of resistance conferred by qMi-C11 and qMi-C14 QTLs. This study lays the groundwork for narrowing down the list of candidate genes previously identified for these QTLs and deepens our understanding of the molecular mechanisms involved in host-plant interaction.

Results

Transcriptome sequencing and mapping

A total of 44 RNA-seq libraries were sequenced. After trimming adapters and removing low-quality reads, the number of reads per library ranged from 8.75 to 20.30 million, with an average of 12.98 million. Of the 574.2 million total reads, 483.99 million (84.28%) were mapped to the reference 'TM-1' genome [28] (Supplemental Figure S1), obtained from CottonGen database (www.cottongen.org) [29].

Differential expression of genes and enriched categories

A total of 2,471 differentially expressed genes (DEGs), corresponding to 1,674 unique genes, were found to be significant (FDR < 0.05; log2 FC ≤ −2 or ≥ + 2) in RKN-treated plants compared to the corresponding untreated controls at five different time points (Table 1, Supplemental Table S1). In response to nematode parasitism, incompatible interactions exhibited twice as many DEGs (1,662) as compatible interactions (809), with upregulated genes comprising 60-70% of both interaction types. The highest number of DEGs during incompatible interactions occurred at 8 DAI (724), followed by 16 DAI (493). Venn diagrams illustrate the number of DEGs across the five time points (Fig. 1).

Table 1 Number of DEGs during compatible and incompatible interactions between cotton and southern root-knot nematode
Full size table
Fig. 1
figure 1

Venn diagrams of the number of DEGs during compatible and incompatible interactions. Diagrams show compatible (susceptible, C201) and incompatible (resistant, M120) interactions between cotton and southern root-knot nematode at a. 4 days after infection (DAI), b. 8 DAI, c. 12 DAI, d. 16 DAI and e. 20 DAI. ‘C’ and ‘M’ correspond to C201 and M120, respectively. Numbers following ‘C’ and ‘M’ correspond to different DAI. ‘Up’ and ‘Down’ correspond to upregulated and downregulated

Full size image

DEGs during compatible interactions

A total of 247 downregulated and 562 upregulated DEGs were significant during compatible interactions in the C201 genotype (Supplemental Table S1). Nearly half of the downregulated genes were observed at 8 DAI, whereas the highest number of upregulated genes (40%) were found at 16 DAI. Downregulated DEGs were enriched in gene ontology (GO) categories such as stress response and oxidation-reduction (Fig. 2a). A total of 27 downregulated genes were common between 4 and 8 DAI, indicating an early transcriptomic repression response during compatible interactions (Fig. 3a). In contrast, relatively few downregulated genes (1 to 4) were recurrent at other time points. Transcriptomic repression at 12 and 16 DAI was less dynamic than at 8 DAI, with only seven and 18 genes significantly downregulated, respectively. Conversely, 59% of upregulation events were recurrent across different time points (Fig. 3b; Supplemental Table S1). Enriched GO classes included biological processes such as oxidation-reduction, response to oxidative stress, and chitin catabolism (Fig. 2a).

Fig. 2
figure 2

Enriched categories of DEGs during compatible and incompatible interactions. ‘Red’ bars correspond to upregulated and ‘blue’ downregulated categories of genes in a. C201 and b. M120

Full size image
Fig. 3
figure 3

Venn diagrams of common DEGs during compatible and incompatible interactions at different DAI. Diagrams show a. downregulated genes in C201, b. upregulated genes in C201, c. downregulated genes in M120 and d. upregulated genes in M120, each at different days after inoculation (DAI). ‘C’ and ‘M’ correspond to C201 and M120, respectively. Numbers following ‘C’ and ‘M’ correspond to different DAI. ‘Up’ and ‘Down’ correspond to upregulated and downregulated. For DEGs common between two different DAI, significance threshold was kept at (FDR < 0.05) given that at least one of the occurrences is also significant at 4-log2 FC.

Full size image

DEGs during incompatible interactions

During incompatible interactions, 655 downregulated and 1,007 upregulated DEGs were significant (Supplemental Table S1). As observed in compatible interactions, nearly half of the downregulated genes occurred at 8 DAI, while the highest proportion of upregulated genes (~ 40%) were also detected at this time point. Approximately 41% of the DEGs were recurrent across at least one other sampling time point (Fig. 3c; Supplemental Table S1). At the earliest time point (4 DAI), only four genes were significantly downregulated, and no specific enrichments were observed at the final time point (20 DAI). Enriched categories for downregulated DEGs included regulation of transcription, cell wall biogenesis, and carbohydrate metabolic process, as well as associated processes like xyloglucan and cellular glucan metabolism and hydrolase activity of O-glycosyl compounds (Fig. 2b). Fewer downregulated genes were observed at 12 and 16 DAI (Supplemental Table S1). Among the upregulated genes, approximately 41% were recurrent across at least one other time point (Fig. 3d; Supplemental Table S1). Key enriched biological processes included metabolic processes like oxidation-reduction, nicotianamine biosynthesis, and chitin and cell wall macromolecule catabolism. Additionally, defense responses to biotic stimuli, including bacteria and fungi, were significantly enriched (Fig. 2b).

Common DEGs during compatible and incompatible interactions

A total of 0, 8, 4, 7, and 9 downregulated genes and 54, 52, 45, 100, and 37 upregulated genes were common between C201 and M120 at 4, 8, 12, 16 and 20 DAI, respectively (Supplemental Figure S2). No genes showed contrasting expression patterns (i.e., opposite directions) between the two genotypes at the same time point. While none of the common downregulated genes were recurrent across time points, several common upregulated genes were recurrent across different time points. For example, 16 DEGs, including DNA-directed RNA polymerase subunit beta, were upregulated in both genotypes throughout the infection. Similarly, two ATP synthase subunit alpha (ATPA) genes were recurrently upregulated. Gh_A03G1169 was upregulated at 4 and 12 DAI during compatible interactions and at all time points except 16 DAI during incompatible interactions. A second gene, Gh_D06G0518, was upregulated at all time points in both interactions except at 20 DAI during incompatible interactions.

Heat maps

Hierarchical clustering of significant DEGs revealed two primary gene clusters (Group A and B), each further divided into three sub-clusters (Fig. 4). Sub-clusters A1 and B1 contained down- and upregulated genes at 8 DAI, whereas A2 and B3 contained those at 16 DAI, respectively – these four clusters were predominantly associated with incompatible interactions. Sub-cluster A3 contained both up- and downregulated genes involved in both types of interactions, whereas B2 was primarily composed of genes that were upregulated during both compatible and incompatible interactions.

Fig. 4
figure 4

Heatmap clustering DEGs during compatible and incompatible interactions. The map shows clustering of DEGs at five different DAI (4, 8, 12, 16 and 20 DAI) in C201 (compatible) and M120 (incompatible). Upregulated DEGs identified with different shades of Red whereas downregulated are shown with different shades of blue. Black shade indicates that the genes do not display significant differential expressions.

Full size image

Differentially expressed genes at the QTL regions

The physical locations of the QTL regions on chromosome 11 (designated as ‘A11’ in subgenome ‘A’) and chromosome 14 (designated as ‘D02’ in subgenome ‘D’) were determined by BLAST homology search against the G. hirsutum genome [28]. This search utilized clone sequences of SSR markers CIR316 and CIR069 spanning qMi-C11, and BNL3644 and NAU5467 spanning qMi-C14 [11, 14]. Additionally, we included chromosome D11, the homeolog of A11, and chromosome A02, the homeolog of D02, to catalog DEGs in the QTL homoeologs. For genotype-specific contrasts at different DAIs, a liberal threshold (FDR < 0.05; log2 FC ≤ −1 or ≥ + 1) was applied to declare DEGs. In total, 47 differentially expressed transcripts corresponding to 38 genes within the two QTLs and their homoeologous regions were identified (Table 2).

Table 2 List of DEGs at the QTL regions
Full size table

DEGs at qMi-C11 QTL region

Three DEGs were identified in the A11 QTL region during incompatible interactions: ‘putative disease resistance RPP13-like protein 1’ (RPPL1: Gh_A11G2836), ‘alpha-xylosidase 1’ (XYL1: Gh_A11G2856), and ‘IQ domain-containing protein’ (IQM3: Gh_A11G2865). RPPL1 was upregulated at 8 DAI, whereas XYL1 and IQM3 were downregulated at 8 and 16 DAI, respectively. Notably, the QTL homoeolog in D11 contained a ‘putative disease resistance protein RGA4’ (Gh_D11G3190), which was upregulated at 8 DAI, whereas the homoeologs of XYL1 and IQM3 (Gh_D11G3231 and Gh_D11G3247) were significantly downregulated at 8 DAI. No significant DEGs were found in the qMi-C11 QTL region during compatible interactions. qPCR indicated the upregulation of Gh_A11G2836 at 8 DAI in both compatible and incompatible interactions (Supplemental Figure S3c).

DEGs at qMi-C14 QTL region

Among the 15 DEGs identified in the qMi-C14 QTL region, three showed differential expression at multiple time points. Notably, ‘receptor-like protein 12’ (Gh_D02G0259) was upregulated at three consecutive time points—8, 12 and 16 DAI—during incompatible interactions. Additionally, six genes from D02 and its homoeolog (A02) were differentially expressed, such as ‘MACPF domain-containing protein’ (At1g14780: Gh_A02G0154 and Gh_D02G0191), ‘cysteine proteinase inhibitor 2’ (CYS2: Gh_A02G0165, Gh_D02G0206), ‘cannabidiolic acid synthase-like 1’ (CBDAS2: Gh_A02G0172, Gh_D02G0216), ‘cytochrome P450 CYP749A22’ (Gh_A02G0220, Gh_D02G0289), ‘selenium-binding protein 2’ (SBP2: Gh_A02G0178, Gh_D02G0241), and ‘glutamine synthetase leaf isozyme, chloroplastic’ (Gh_A02G0242, Gh_D02G0311). Both SBP2 (Gh_A02G0178) and ‘glutamine synthetase’ gene were also differentially expressed during compatible interactions. Additionally, four genes, including ‘receptor-like protein 12’ (RLP12: Gh_D02G0257), showed differential expression during compatible interactions. qPCR analysis of the two putative ‘receptor-like protein 12’ candidate genes, Gh_D02G0257 (Fig. 5a) and Gh_D02G0259 (Fig. 5b), confirmed their significant overexpression in response to nematode parasitism. Gh_D02G0259 was upregulated beginning at 8 DAI and remained so through 16 DAI during incompatible interactions, further confirming the RNA-seq findings.

Fig. 5
figure 5

Quantitative real-time PCR (qRT-PCR) of qMi-C14 candidate genes a. Gh_D02G0257 and Gh_D02G0259. Charts show mean fold-change (with standard error bars) in expression at six different time points after RKN infection. Asterisks mark significant upregulation in expression in inoculated samples compared to non-inoculated plants and as determined by t-test of ΔCt values (P ≤ 0.05) using two biological replicates.

Full size image

Specific categories of nematode-responsive genes

Stress- and defense-responsive genes

The CottonFGD search function identified defense-related genes associated with the JA and SA pathways, resistance gene analogs (RGAs), pathogenicity-related (PR) genes, and ethylene biosynthesis. Genes involved in plant defense responses, including those in the alpha-linolenic acid metabolism (JA pathway), RGAs and PR genes, were differentially regulated in both compatible and incompatible interactions (Supplemental Table S2). In the alpha-linolenic acid pathway, which is a precursor of JA, two and four genes were repressed, while four and 14 genes were overexpressed during compatible and incompatible interactions, respectively (Supplemental Figure S4).

Further, Arabidopsis orthologs were identified for the significant DEGs found during G. hirsutum-M. incognita interactions. The NEMATIC database was used to identify common interactions between Arabidopsis and cotton in response to RKN (Supplemental Table S3). Focusing on nematode-responsive genes in the ‘stress’ category led to the identification of 25 common DEGs in the stress category during RKN-host interactions (Supplemental Table S4). Of these, six genes were downregulated during compatible interactions in both Arabidopsis and cotton with RKN.

Cell wall biogenesis and remodeling

Nematode-responsive DEGs involved in cell wall biogenesis, cell wall macromolecule catabolism, cell wall modification, and cell wall organization were identified (Supplemental Table S5). Incompatible interactions exhibited greater diversity in GO categories and a higher number of DEGs in response to nematode parasitism than in compatible interactions. These DEGs were categorized into four major molecular functions: xyloglucan:xyloglucosyl transferase activity, chitinase activity, pectinesterase activity, and expansin enzyme family (Supplemental Figure S5a-d). For instance, all five genes encoding endochitinases were upregulated during compatible interaction at 16 DAI, with fold changes ranging from 2.6 to 20.2. While three of these genes were also upregulated during incompatible interactions, five other endochitinases were specific to incompatible interactions. The highest fold change was observed for CHIA1 (FC = 54.7). Additionally, 17 DE events corresponding to seven members of ‘fasciclin-like arabinogalactan proteins’ were downregulated at 8 and/or 16 DAI during incompatible interactions, whereas only one gene was downregulated at 8 DAI during compatible interactions. Furthermore, several DEGs in the phenylpropanoid biosynthetic pathway indicated changes in cellular biogenesis in response to nematode infection (Supplemental Figure S6a-b).

Cell growth and development

Several genes involved in cell proliferation were differentially regulated in response to nematode parasitism (Supplemental Table S6). For example, during incompatible interactions, ‘phytosulfokine 3’ (PKS3) growth factors, which facilitate the positive regulation of plant cell proliferation, were upregulated at 8 and 16 DAI, whereas two putative PSK6 homoeolog genes were downregulated at 12 DAI. Additionally, three types of ‘cyclin-dependent protein serine/threonine kinase regulators’ (CYCA, CYCB and CYCD), along with two ‘G2/mitotic-specific cyclin S13-6’ genes, which regulate the cell cycle, were downregulated during incompatible interactions. Identical expression patterns were observed for two classes of ‘COBRA-like proteins’ (COBL1 and COBL7), which regulate cell expansion and cellulose crystallinity. Specifically, two sets of homoeologous COBL7 and COBL1 genes were downregulated at 8 DAI and 20 DAI, whereas COBL4 was downregulated at 16 DAI.

Transcription factors

Several transcription factor (TF) families altered expression patterns following RKN infection, likely driving widespread transcriptional changes. A total of 174 TFs across 29 families displayed differential expression in response to RKN infection (Fig. 6; Supplemental Table S7). During compatible interactions, 11 TF families were downregulated, and six were upregulated, whereas 18 and 17 families were differentially expressed during incompatible interactions. Among these, WRKY was one of the most abundant and exclusively overexpressed TF families of both interaction types. Similarly, ‘trihelix’ was repressed only in compatible interactions, while ‘GATA’ was downregulated only in incompatible interactions. In addition to these TF families, the expression of key genes involved in transcriptional regulation, such as ‘DNA-directed RNA polymerase’, was also altered. Notably, two to four ‘DNA-directed RNA polymerase’ genes (rpoB, rpoB1 and rpoB2) were upregulated with an average fold change of 22 to 26 (average log2 FC > 4) across all sampling time points in both interaction types, supporting enhanced transcriptional activity in response to RKN infection.

Fig. 6
figure 6

Distribution of G. hirsutum transcription factors (TF) families during RKN infection and development. Downregulated (a and c) and upregulated (b and d) TF families during compatible (a and b) and incompatible (c and d) interactions. Corresponding number of TF genes are provided following TF family identifiers (after comma)

Full size image

Carbon and energy metabolism

A wide array of genes involved in glycolysis/gluconeogenesis, amino sugar and nucleotide sugar, and starch and sucrose synthesis pathways were differentially expressed in response to nematode parasitism during both compatible and incompatible interactions (Supplemental Table S8; Supplemental Figure S7-S10). For example, three genes corresponding to ‘basic vacuolar isoform of glucan endo-1,3-beta-glucosidase’ were upregulated (log2FC range from 4.1 to 5.3), whereas several plasma membrane-based members were moderately downregulated during incompatible interactions at 8 DAI. These genes were subsequently upregulated in both interaction types. Similarly, genes encoding ‘UDP-glucuronate 4-epimerase’ were downregulated in compatible interactions (two genes at 8 DAI) and in incompatible interactions (one gene at 16 DAI). Additionally, two ‘UDP-glucose 4-epimerase genes’ were upregulated during incompatible interactions at 16 DAI.

During compatible interactions, all DEGs associated with the glycolysis/gluconeogenesis pathway were downregulated at 4 and 8 DAI. For example, homoeologous copies of ‘alcohol dehydrogenase class-P’ (ADH1) and ‘pyruvate decarboxylase 1’ (PDC1) were downregulated at 4 and 8 DAI, whereas those of ‘ATP-dependent 8-phosphofructokinase 3’ (PSK3) were downregulated at 8 DAI. During incompatible interactions, ADH1 was downregulated at 8 DAI but upregulated at 16 and 20 DAI. Additionally, a set of homoeologous PSK3 genes was upregulated at 8 DAI during incompatible interactions. Similarly, ‘cytosolic fructose-1,6-bisphosphatase’, a key enzyme in the sucrose biosynthetic pathway, was upregulated at both 8 and 16 DAI.

Discussion

Differential gene expression analysis using global transcriptome data is an effective approach to understand the regulatory networks of genes and pathways involved in plant host-pathogen interactions. Based on histopathological evidence, the current working hypothesis suggests that resistance in the Auburn 623 RNR source is conferred by two major genes with different mechanisms of resistance. The first mechanism occurs soon after infection (8–12 DAI), preventing juvenile nematodes from developing functional feeding sites, while the second mechanism occurs later (14–20 DAI), impeding the development of nematodes into egg-laying females [6,7,8,9,10]. Mechanisms of RKN resistance are dynamic biological processes, with nematode suppression occurring at multiple stages of the life cycle. Prior RNA-seq experiments, which were limited to a single time point (10 DAI) [27] or two time points (12 and 30 DAI) [26], provided important insights but were insufficient to fully capture the dynamics of RKN resistance. Therefore, the current study conducted a more detailed time-course RNA-seq experiment, constructing and sequencing 44 libraries across five time points (4, 8, 12, 16 and 20 DAI) from RKN-treated and untreated (controls) samples of C201 (susceptible) and M120 (resistant) genotypes. The DEGs identified between genotypes with and without nematode infections, as well as key candidate genes in the QTL regions, provide a more comprehensive view of the molecular events underlying RKN resistance in cotton.

Massive transcriptional reprogramming of defense response underlies cotton-RKN interactions

A dynamic transcriptional response distinguishes incompatible from compatible interactions. At 4 DAI, more DEGs were detected in compatible interactions (Table 1), supporting the idea that the gene types, rather than expression level, influence nematode survival [30]. From 8 DAI onward, incompatible interactions displayed stronger transcriptomic responses, both genome-wide and at specific QTL regions, despite higher nematode loads in C201 (Supplemental Table S9). The near-immunity in M120 appears to be driven by enhanced basal defense gene expression and selective activation of R genes—hallmarks of RKN resistance in several plant species, including tobacco [31, 32] and alfalfa [33]. In contrast, susceptibility involves active suppression of host defenses during nematode parasitism [29].

Resistance in M120 aligns with the zig-zag model of the plant immune system [34], involving two defense layers: nematode-associated molecular patterns (NAMPs)-triggered immunity (PTI) and effector-triggered immunity (ETI). PTI, associated with pathogenesis-related (PR) genes, is initiated by NAMP recognition via pattern recognition receptors (PRRs), triggering downstream cellular signaling cascades, including the activation of ‘mitogen-activated protein kinases’ (MAPKs), the production of reactive oxygen species (ROS), and the initiation of JA or SA pathways [35]. We detected 37 Mi genes (mapped to the M. incognita genome) encoding effector proteins (Supplemental Table S9) and several Gh genes encoding ‘receptor-like kinases’ (RLKs), receptor-like proteins’ (RLPs) and ‘wall-associated receptor kinases’ (WAKs) that were differentially expressed in RKN-treated samples. Six PR genes were differentially regulated during incompatible interactions, whereas only two were upregulated in compatible interactions. Furthermore, homoeologous RLK genes (At1g72540: Gh_A05G1164 and Gh_D05G1341) were downregulated in compatible interactions at 4 and 8 DAI, consistent with the general downregulation of PRR-encoding genes in susceptible hosts, as observed in potato [22].

One of the early PTI and/or ETI cellular responses involves the rapid generation of reactive oxygen species (ROS), which activates defense genes and reinforces cell wall [36]. In this study, several key ROS-associated genes were differentially expressed. For example, three different ‘peroxidases’ were downregulated at 8 DAI, whereas five ‘peroxidase’ genes and two ‘respiratory burst oxidases’ (RBOHA) were upregulated at 16 DAI during compatible interaction. This suggests nematodes suppress ROS-mediated defense during early infection, mirroring observations in potato [22]. Furthermore, the early (8 DAI) downregulation of ROS-scavenging genes coincides with histopathological evidence of localized HR during incompatible interactions [6], suggesting that while HR is inhibited in C201, it remains active in M120, though it does not involve classic HR-associated cell death.

’Patatin-like proteins’ (PLPs) are key mediators of HR in various host-pathogens interactions [37,38,39]. In this study, two PLP1 and four PLP2 genes were upregulated at 8 DAI during incompatible interaction (average log2FC = 4.85), potentially enabling limited HR, as suggested by histopathological studies in cotton [6]. Finally, the early induction of ‘dirigent’ genes in both this study and de Deus Barbosa et al. [30] suggests that during incompatible interactions, the endodermal barrier system actively impedes nematode migration toward the vascular tissues [40]. This barrier, absent in susceptible lines, may contribute to the higher success rate of nematodes in establishing feeding sites in C201.

Several TF families show differential expression in response to M. incognita infection

The significant upregulation of the transcription catalysts (rpoB, rpoB1 and rpoB2) during cotton-RKN interactions supports widespread transcriptional reprogramming in both compatible and incompatible interactions. Transcription factors (TFs), such as WRKY genes, play a crucial role in adaptive biotic stress responses. While WRKY genes were upregulated in both compatible and incompatible interaction types, they exhibited earlier induction, greater diversity, and in certain cases, a higher degree of upregulation during incompatible interactions. Notably, among the WRKY genes, the SlyWRKY75 gene, a JA-pathway regulator [41], showed a striking 17- and 86-fold increase at 8 and 16 DAI, respectively, in RKN-treated M120, indicating a strong JA response. Additionally, six WRKY70 genes were upregulated during compatible interactions, consistent with their role in mediating crosstalk between SA and ET/JA pathways [42], as observed in soybean- M. incognita interactions [43].

Ethylene response factors (ERF) are closely associated with hormonal signal transduction pathways that are involved in plant defense responses. In this study, the predominant downregulation of ERFs during compatible interactions suggests a suppressed ethylene response. However, multiple ERFs were differentially regulated during incompatible interactions, reflecting their complex role in RKN-cotton interactions, although reports of their functions remain inconclusive [44].

MYB TFs, known to mediate phytohormones defense signaling [45], were mostly downregulated in compatible interactions but displayed mixed regulation in incompatible interactions. Two MYB108 genes were downregulated in compatible interactions, while four were upregulated in incompatible interactions, indicating enhanced downstream defense signaling responses. Similarly, three NAC TFs, which regulate cell wall remodeling in response to RKN infection [19] (Gh_A03G0024, Gh_A08G1428 and Gh_D03G1627), were significantly upregulated in both compatible and incompatible interactions but showed sustained overexpression in incompatible interactions at later stages.

During cotton-RKN interactions, hormone-modulated responses abound

Phytohormones regulate plant defense by triggering defense molecules such as phenylpropanoids, phytoalexins, and pathogenicity-like proteins. Auxin-related genes exhibit differential expression during plant-nematode interactions [46, 47]. Indole-3-acetic acid (IAA), the primary natural form of auxin, is synthesized from tryptophan (Trp) via a two-step pathway: Trp is first converted to indole-3-pyruvate (IPA) by the TAA family aminotransferases, then IPA is converted to IAA by the YUC family flavin monooxygenases [48]. In this study, a TAA gene (TAR4) was upregulated at 16 DAI in both interactions, while YUC10 genes (Gh_A08G1012 and Gh_D08G1285) were upregulated only in incompatible interaction at 8 DAI. Additionally, NGA4, a TF potentially linked to YUC overexpression [49], was upregulated at 16 DAI during incompatible interactions. Similarly, auxin-response factors ARF5 and ARF6 were overexpressed between 8 and 16 DAI in incompatible interactions, while ARF7 and ARF19 were upregulated only during compatible interactions at 20 DAI.

Nematodes also influence hormonal pathways, linked to root apical meristem maintenance, lateral root formation, and nematode feeding site development. Several SAUR and auxin-responsive genes were downregulated in compatible interactions, with mixed regulation in incompatible interactions. For example, SAUR41 homologs (Gh_A10G0478 and Gh_D10G0502), involved in Arabidopsis cell expansion and root meristem patterning [50], were downregulated at 4 and 8 DAI in compatible interactions. Similarly, three LBD genes, which encode LOG domain-containing proteins associated with auxin-regulated root formation [51], were downregulated in response to nematode parasitism.

Genes involved in Jasmonic acid (JA) and Salicylic acid (SA) biosynthesis, key pathways in host resistance, were differentially expressed, particularly at 8 and 16 DAI. For instance, two genes encoding allene oxide synthase (CYP74A), a key enzyme in JA biosynthesis, were strongly upregulated (log2FC > 4.6) in incompatible interactions at 8 and 16 DAI. Other key genes of JA biosynthesis pathway were also upregulated, including AOC4 (allene cyclase) and LOX1.5 (lipoxygenase) [52] at 16 DAI, as well as several ACX (Acyl-coenzyme A oxidase) genes at 8 DAI. Similarly, the PR-1 gene, a key SA-regulated defense gene [53], showed strong upregulation (FC = 32.6) in the incompatible interaction at 8 DAI. Additionally, several PAL genes, integral to phenylalanine-mediated SA production, were upregulated in both interactions, which is consistent with the response of C. metuliferus to RKN [20].

Hydrolases contribute to plant resistance to fungus, viruses, and plant-parasitic nematodes like Heterodera rostochiensis [54,55,56]. These enzymes, upregulated in incompatible interactions in cotton [30] and tobacco [32], are associated with hypersensitive responses (HR), which involve toxin accumulation, calcium flux, and cell death [34]. Such responses may deprive nematodes of nutrients, leading to giant cell collapse [20, 24]. In this study, three HGN1 genes (Glucan endo-1,3-beta-glucosidase, basic vacuolar isoform), which activate toxin precursors, were strongly upregulated (average FC = 27.79) during incompatible interaction at 8 DAI, with one gene also upregulated at 16 DAI.

Ethylene-related genes, known to be induced in nematode feeding sites [57], also enhance systemic expression of basic chitinases involved in nematodes’ defense [58, 59]. In this study, endochitinases were mostly upregulated in both interactions, but incompatible interactions showed earlier induction, greater diversity, and higher upregulation. For instance, endochitinase 1 expression had an average FC of 14.6 in compatible interactions at 16 DAI, but 54.7 and 27.4 in incompatible interactions at 8 and 16 DAI, respectively, aligning with soybean responses [60].

DEGs involved in metabolic/transport activities and giant cell formation

Nematodes exploit host systems in compatible interactions to divert nutrients from the giant cells to their feeding sites. During these interactions, six aminotransferases were downregulated and three were upregulated, indicating reduced catabolism and increased allocation of amino acids to nematodes. Conversely, in incompatible interactions, the trend was reversed, with six aminotransferases upregulated and two downregulated, suggesting enhanced catabolism and restricted nutrient flow to nematodes. Also, genes linked to oligopeptide transport were mostly upregulated at 8 DAI (five upregulated and one downregulated) but showed repression at 16 DAI, reducing nutrient supply to giant cells. PSK (phytosulfokine) signaling, involved in nematode-induced redifferentiation in Arabidopsis [61], was repressed. In this study, four PSK6 genes were downregulated in compatible interactions, and two in incompatible interactions. Gh_A01G1788 was recurrently downregulated at 12 and 16 DAI in both interactions and at 20 DAI in incompatible interactions. A homoeologous PSK6 gene pair (Gh_A01G1788 and Gh_D01G2029) was also downregulated at 12 DAI across both interactions. Downregulation of ‘E2F-transcription factor’ gene (E2FE) and its modulator, the retinoblastoma-related (RBR) protein, at 8 DAI, along with repression of six cyclin-encoding genes, suggest that characteristic endoreduplication in giant cells was repressed during incompatible interactions [62]. The lower galling index and smaller galls in resistant genotypes, compared to the larger galls in susceptible genotype, suggest that differential regulation of PSK3, PSKR, RBR, E2F and cyclin genes during early infection plays a crucial role, as observed in soybean [43].

Basal responses to nematode-secreted cell wall-degrading enzymes involve strengthening the cell wall barrier through the action of endochitinases, pectinesterases, expansins, fasciclin-like arabinogalactan proteins, and dirigent proteins. Many of these genes were differentially expressed here and are conserved, nematode-responsive orthogroups in resistant genotypes [23]. Further, giant cell formation requires altered expression of genes regulating cell wall remodeling. Four such gene families—xyloglucan:xyloglucosyl transferase, chitinase, pectinesterase, and expansin—showed notable changes in both interactions. Additionally, phenylpropanoid biosynthesis, essential for physical and biochemical defenses (e.g., lignin and suberin), is critical for pathogen resistance [63]. Ye et al. reported that timely and efficient activation of this pathway distinguishes incompatible from compatible interactions [20]. This study confirmed significant expression changes in phenylpropanoid-related genes, with incompatible interactions exhibiting earlier induction, greater diversity, and higher upregulation.

DEGs at the QTL regions emerge as key RKN-resistance candidates

DEGs within the QTL intervals of qMi-C11 and qMi-C14 are of particular interest in this study because transcriptional changes observed during incompatible interactions are likely driven by one or more genes that interfere with nematode parasitism [7]. A total of 22 genes within these QTL regions were differentially expressed during incompatible interactions, compared to only six during compatible interactions.

Previous findings showed that nematode survival and reproduction were markedly reduced in resistant lines, with a lower proportion of nematodes progressing to J3 (third-stage juvenile) and J4 (fourth-stage juvenile) stages in M120 compared to C201 at 8 DAI (3.4% vs 11.9%) [7]. The delay in J2s reaching the central cylinder and becoming sedentary aligns with the first two sampling time points in this study (4 and 8 DAI), where the highest number of DEGs was detected at 8 DAI. This supports the idea that the qMi-C11 QTL mediates M120-specific defense responses early in the infection cycle. During this phase, J2s in compatible interactions would have developed into J3s, and hypertrophied giant cells would already be evident [6, 7]. In contrast, in incompatible interactions, necrotic and dying cells in the central cylinder, indicative of HR, would be prominent [6].

The period between the second and third sampling time points (8 to 12 DAI) corresponds to the HR phase, ultimately leading to nematode feeding site degeneration by the fourth and fifth sampling time points (16 to 20 DAI). Interestingly, the second highest number of DEGs were observed at 16 DAI, suggesting that the qMi-C14 QTL activates defense responses that disrupt nematode development in the later stages of the infection cycle during incompatible interactions. This disruption likely prevents RKNs from developing into mature females, thereby reducing their reproductive capacity [7]. While qMi-C14 appears critical for defense responses at later stages, the precise mechanism by which it confers resistance remains unclear and warrants further investigation.

For the qMi-C11 QTL on chromosome 11, our previous study identified a pair of colocalized putative disease resistance genes (RPP13-like protein 1: Gh_A11G2835 and Gh_A11G2836) that exhibited contrasting expression patterns in RKN-treated M120 and C201 genotypes [26]. While Gh_A11G2835 did not show differential expression in the current study, two other RPPL1 genes (Gh_A11G2836 and Gh_A11G2838) were downregulated in RKN-treated M120 versus C201 contrasts at 12 and 20 DAI. However, qPCR did not validate the downregulation at these time points (Supplemental Figure S3c). Notably, Gh_A11G2836 was upregulated during incompatible interactions at 8 DAI in both RNA-seq and qPCR analyses, which suggests dynamic regulation in response to nematode pathogenicity.

Another candidate gene, ‘U-box domain-containing protein 21’ gene (PUB21: Gh_A11G3090), located within the qMi-C11 QTL region, was previously reported to be upregulated in the resistant line at both early (12 DAI) and late (30 DAI) infection stages [26]. Interestingly, several E3 ubiquitin ligases, including PUB21 orthologs in Arabidopsis, are known to respond to chitin, a well-established plant defense elicitor in fungal cell walls and the exoskeletons of insects and nematodes [64]. While RNA-seq analysis in the current study did not identify Gh_A11G3090 as differentially expressed, qPCR indicated its upregulation at four and five different time points during compatible and incompatible interactions, respectively (Supplemental Figure S3d). The inconsistencies in DEGs identified within the qMi-C11 region across RNA-seq studies or between RNA-seq and qPCR approaches may reflect subtle expression changes that initiate downstream resistance responses or suggest that post-transcriptional modifications are involved in qMi-C11-mediated resistance. Nevertheless, Gh_A11G2836 and Gh_A11G3090 remain the putative candidate genes within the qMi-C11 QTL region.

In a recent study, Ojeda-Rivera et al. employed a single sampling time point (23-day-old seedlings) for their comparative transcriptomic analysis. They concluded that the basal expression of candidate genes might be constitutively higher in ‘Acala NemX’, potentially conferring resistance without imposing a fitness penalty [27]. Notably, higher basal expression levels were observed in two homologs of an Arabidopsis NBS-LRR R protein gene (AT5G36930), namely Gohir.A11G297600 (log2FC = 1.28) and Gohir.D11G312800 (log2FC = 1.99), in mock-treated libraries of ‘Acala NemX’ compared to ‘Acala SJ2’. As such, these genes were proposed as the most promising candidates for RKN-resistance in ‘Acala NemX’, particularly as they are located within the rkn1 QTL region [27]. In the current study, these genes did not show differential expression in any of the contrasts examined, whether within genotypes (mock vs. treated) or between genotypes (mock vs. mock or treated vs. treated) across all time points. Furthermore, there were no shared DEGs between the two studies within the qMi-C11 QTL region. This disparity is not unexpected, as the resistance source in this study, ‘Auburn 623 RNR’ (qMi-C11), and the ‘Acala NemX’-sourced rkn1 QTL are likely different, even though the QTLs are mapped to the same chromosome region but in adjacent, yet distinct, marker-intervals. Further, resistance in ‘Acala NemX’ is conferred by a recessive gene, whereas in ‘Auburn 623 RNR’, it is conferred by dominant gene(s) [12, 65], strongly indicating that they have unique modes of action [12, 27]. Collectively, these results support the hypothesis that rkn1 and qMi-C11 represent unique and independent QTLs.

For the qMi-C14 QTL on chromosome 14, our prior research identified two RLP12 genes (Gh_D02G0257 and Gh_D02G0259) as promising candidates based on both RNA-seq and qPCR analyses [15]. In this study, comparative contrasts of RKN-treated genotypes revealed differential expression of these genes. Specifically, Gh_D02G0257 was relatively overexpressed in the susceptible line at 12, 16 and 20 DAI, whereas Gh_D02G0259 was consistently overexpressed in the resistant line across all time points (Supplemental Table S10). Further, RNA-seq analysis showed that Gh_D06G0257 was upregulated during compatible interaction at 16 DAI, whereas Gh_D06G0259 was upregulated during incompatible interactions at 8, 12 and 16 DAI, with highest fold change (log2FC = 2.1) at 16 DAI. These timelines align closely with histopathological observations from da Silva et al., who reported a six-fold reduction in nematodes progressing to females and a tenfold decrease in egg production in the resistant genotype (M120) at 16 and 20 DAI, respectively [7]. These findings strongly suggest that RLP12 genes within the chromosome 14 QTL region may act as plant receptors capable of detecting nematode parasitism and initiating signaling pathways that activate innate immune responses, thereby restricting RKN reproduction in the resistant genotype.

It is worth noting that neither this study nor our earlier study [26] showed significant upregulation of Gh_D02G0276, a gene previously proposed as the casual gene for the qMi-C14 [16], in response to RKN infection (Supplemental Figure S3f). This gene was identified in MAGIC RIL lines derived from M-240 RNR, whereas our studies used M-120 RNR; however, both share Auburn 623 as their source of resistance. Similarly, Ojeda-Rivera et al. reported significant upregulation of two putative candidate genes on chromosome 14, Gohir.D02G023200 (log2FC = 3.28) and Gohir.D02G023300 (log2FC = 2.73), in the susceptible genotype ‘Acala SJ2’ in response to nematode infection, whereas their expressions did not differ significantly in the resistant genotypes [27]. Although ‘Acala NemX’ does not carry the qMi-Chr14 QTL [66], its presumed source, ‘WMJJ’ [13, 14], also showed no differential expression of these genes at 10 DAI [27]. Interestingly, basal expression comparisons (control libraries) revealed a significant difference for Gohir.D02G023200 (log2FC = 4.42) between ‘Acala NemX’ and ‘Acala SJ2’, leading researchers to propose that ‘Acala NemX’ employs constitutive overexpression of these genes, rather than infection-induced changes, to suppress nematode infection [27]. These discrepancies in gene expression, observed in response to RKN infection, may reflect genotypic differences between ‘WMJJ’ and ‘M-120 RNR’.

Conclusion

Basal defense responses were activated both during compatible and incompatible interactions, although the responses were more pronounced in RKN resistant genotypes compared to the susceptible lines. In summary, nematode responsive genes related to defense pathways were often repressed during compatible interactions, whereas earlier induction, greater diversity, and a higher degree of upregulation of those genes were archetypal of incompatible interactions. A wide spectrum of disease resistance and putative resistance genes, pathogenesis-related genes, and genes corresponding to ligands and receptors were differentially expressed in response to nematode parasitism in G. hirsutum. These genes were mapped across the cotton genome and include potential candidate genes Gh_A11G3090 (PUB21) and Gh_A11G2836 (RPPL1) in chromosome 11 and Gh_D02G0257 (RLP12) and Gh_D02G0259 (RLP12) in chromosomes 14 QTL regions.

Comparative transcriptomics of a whole root system has limitations, as it may mask, dilute, and/or modify transcriptomic changes at localized interaction sites, such as specific cells or tissues [46, 67]. However, using controls at each time point and comparing near isogenic lines can justify this broader approach. Also, since nematode infection and development within plant roots are not synchronous, pooling root tissues over several days, as done by Shukla et al. [21], may better capture these dynamic interactions. Ongoing research on comparative transcriptomics of near-isogenic lines with QTLs on chromosomes 11 and 14 aims to validate these findings and identify candidate genes for future functional studies that could explain the ‘near immunity’ levels of resistance in M-120 RNR.

Methods

Plant material and greenhouse evaluation

Two Upland cotton genotypes, M-120 RNR (M120) and Coker 201 (C201), with contrasting responses to RKN parasitism, were used in this study. C201 is an obsolete cotton variety that is susceptible to M. incognita, while M120 is a near-isogenic line highly resistant to both galling and RKN egg production. It was developed by crossing Auburn 634 RNR with C201, followed by several backcrosses to C201, before being self-pollinated and publicly released as M-120 RNR [68]. Selfed progenies of these two publicly available germplasm lines are continually developed and maintained in our breeding program. The greenhouse experiments used in this study are described in the companion study [7]. Briefly, seeds of both genotypes were germinated in vermiculite. Fourteen days after germination, seedlings were individually transplanted into 10.6 cm x 10.6 cm x 12.4 cm pots filled with approximately 500 ml of steam-pasteurized soil (Tifton loamy sand). At transplant, half of the seedlings of each line were infected with 4,000 M. incognita J2s per plant, while the other half was left inoculated, serving as control. The entire root systems from inoculated and un-inoculated M120 and C201 seedlings were carefully removed from the potting-mix at 4, 8, 12, 16 and 20 days after inoculation (DAI), washed with deionized water, dried with sterilized paper towels and immediately frozen in liquid nitrogen. Root tissues of uninoculated seedlings at 0 DAI were also harvested as controls. Samples were stored at −80°C until RNA extraction. In total, 44 samples were collected, representing the resistant (M120) and susceptible genotype (C201), two treatments (RKN-treated vs. control), five time points (4, 8, 12, 16 and 20 DAI), and two biological replications for each condition. Phenotypic data from the companion study [7] corresponding to M. incognita counts at different developmental stages in the roots of C201 and M-120 RNR are presented in Supplemental Tables 10 and 11.

Library construction and sequencing

Frozen root samples from C201 and M120, collected at six different time points, were individually ground with mortar and pestle in liquid nitrogen. Total RNA was then extracted using the Spectrum™ Plant Total RNA extraction kit (Sigma-Aldrich), following the manufacturer’s protocol. Library construction was performed by the Georgia Genomics and Bioinformatics Core at the University of Georgia, Athens, GA, using the Kapa Stranded RNA-seq kit (Roche), and sequencing was conducted on a NextSeq PE150 High Output flow cell platform. The RNA sequencing data from this study have been deposited in the NCBI database under BioProject ID PRJNA1062816 and SRA accessions SRX23150901 to SRX23150944.

Transcript assembly and differential expression analyses

The quality of raw reads was checked using the FastQC program (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Low-quality bases and adapter sequences from the paired reads were trimmed using the Trimmomatic v0.30 program [69] (http://www.usadellab.org/cms/?page=trimmomatic). Trimmed reads were then mapped to the G. hirsutum reference genome [28] using HISAT2 v2.1.0 (https://daehwankimlab.github.io/hisat2/) [70]. Sequence alignment files were input into the software Stringtie v1.3.3 to assemble potential transcripts (http://ccb.jhu.edu/software/stringtie/) [71]. A Python script was used to obtain gene-level raw counts from each library for differential gene expression analysis using the R package DESeq2 (https://bioconductor.org/packages/release/bioc/html/DESeq2.html) [72].

DESeq2 package, incorporated in SARTools [73], an R-based comprehensive differential expression (DE) analysis platform, was run with default parameters for RNA-seq data analysis. Adjusted P-value of ≤ 0.05, based on Benjamini-Hochberg procedure, was used to declare significant DE genes, while ‘cooksCutoff’ was used for outlier detection and ‘independentFiltering’ was used to remove low count genes. Contrasts corresponding to DE of genes within the genotypes (C201 or M120) under different conditions (RKN treatment vs. control) and time points (4, 8, 12, 16 and 20 DAI), as well as between the genotypes (M120 vs. C201) under identical conditions and time points, were also performed. The web-based Next-Generation Clustered Heat Map (NG-CHM) (https://www.ngchm.net/) was employed using Ward’s agglomerative hierarchical clustering and Euclidean metric to cluster significant DEGs and generate heat maps [74].

Functional characterization, gene enrichment analysis, and real-time qPCR

Differentially expressed genes (DEGs) involved in plant defense response and disease resistance, in signaling (such as ligands, receptors, and transcription factors), or in cell wall and cellular functions were identified using keyword search tools available in the Cotton Functional Genomics Database CottonFGD (https://cottonfgd.net/) [75]. Gene ontology (GO) analysis, gene list enrichment, and KEGG pathway analyses were performed using the ‘Data Fetch and Enrichment’ tool in CottonFGD. Default parameters of P ≤ 0.0001 and ‘Minimum Gene Number’ of 3 were used to produce enrichment categories. Graphical pathways were obtained from KEGG Mapper’s Search&Color Pathway web tool (https://www.genome.jp/kegg/tool/map_pathway2.html). FunRich [76] and Genvenn (http://habtom.github.io/biojs-vis-genvenn/examples/index.html) were used for the graphical display of Venn diagrams. Expression patterns of common genes and orthologs of identified DEGs were further explored using the NEMATIC database [77], an Excel-based resource harboring several A. thaliana-RKN transcriptomic studies.

Details of primer design and RT-qPCR-based relative gene expression analysis (normalized against endogenous GhACT4) for putative candidate genes within the qMi-C11 and qMi-C14 QTL regions have been previously reported [26]. Briefly, two biological replications of C201 and M120 from six different time points post-inoculation (including 30 DAI, which was not sequenced) and their respective controls were analyzed following Kumar et al. [26].

Data availability

The plant materials used in this study are publicly available resources. The raw RNA-seq data supporting the conclusions of this article are available through the National Center for Biotechnology Information database (https://www.ncbi.nlm.nih.gov/), under BioProject (I.D. PRJNA1062816) titled ‘Global transcriptome of Southern root-knot nematode resistant and susceptible Upland cotton genotypes.’ Direct link to datasets: http://www.ncbi.nlm.nih.gov/bioproject/1062816. Processed data are included within the manuscript or provided as Supplemental Files.

References

  1. Moens M, Perry RN, Starr JL. Meloidogyne species – a diverse group of novel and important plant parasites. In: Perry RN, Moens M, Starr JL, editors. Root-knot nematodes. Wallingford: CABI; 2009. p. 1–17.

    Google Scholar 

  2. Faske TR, Mueller J, Becker JO, Bernard EC, Bradley C, Bond J, Desager J, Eisenback J, Grabau Z, Hu J. Summarized distribution of the Southern Root-Knot Nematode, Meloidogyne incognita, in field crops in the United States. Plant Health Progress. 2023;24(4):522–4.

    Article  Google Scholar 

  3. Cotton Price S. 2023–2024. Annual Report 105(13): 1-21. United States Department of Agriculture, Agricultural Marketing Service, Cotton and Tobacco Program, Memphis, TN. Retrieved from chrome-extension://efaidnbmnnnibpcajpcglclefindmkaj/https://www.ams.usda. gov/mnreports/cnaacps.pdf.

  4. Faske T, Sisson A. Cotton disease loss estimates from the United States —. 2023. Retrieved from https://cropprotectionnetwork.org/publications/cotton-disease-loss-estimates-from-the-united-states-2023. https://doiorg.publicaciones.saludcastillayleon.es/10.31274/cpn-20240219-0.

  5. Hajihassani A, Lawrence KS, Jagdale GB. Plant parasitic nematodes in Georgia and Alabama. In: Subbotin SA, Chitambar JJ, editors. Plant parasitic nematodes in sustainable agriculture of North America: Vol2 - Northeastern, Midwestern and Southern USA. Cham: Springer International Publishing; 2018. p. 357–91.

    Chapter  Google Scholar 

  6. Lopes CML, Suassuna ND, Cares JE, Gomes ACMM, Perina FJ, Nascimento GF, Mendonca JSF, Moita AW, Carneiro RMDG. Marker-assisted selection in Gossypium spp. for Meloidogyne incognita resistance and histopathological characterization of a near immune line. 2020;216:19. https://doiorg.publicaciones.saludcastillayleon.es/10.1007/s10681-020-2554-7.

  7. da Silva MB, Davis RF, Kumar P, Nichols RL, Chee PW. Resistance Quantitative Trait Loci qMi-C11 and qMi-C14 in Cotton Have Different Effects on the Development of Meloidogyne incognita, the Southern Root-Knot Nematode. In. United States: American Phytopathological Society; 2019. p. 853.

    Google Scholar 

  8. da Silva MB, Davis RF, Nichols RL, Kumar P, Chee PW. The effect of two QTLs for resistance to Meloidogyne incognita in cotton on nematode egression from roots. J Nematol. 2021;52:e2020-122. https://doiorg.publicaciones.saludcastillayleon.es/10.21307/jofnem-2020-122.

  9. Jenkins JN, Creech RG, Tang B, Lawrence GW, McCarty JC. Cotton resistance to root-knot nematode: II. Post-penetration development. Crop Sci. 1995;35(2):cropsci19950011183X003500020013x.

    Article  Google Scholar 

  10. Wubben MJ, Gaudin AG, McCarty JC Jr, Jenkins JN. Analysis of cotton chromosome 11 and 14 root-knot nematode resistance quantitative trait loci effects on root-knot nematode postinfection development, egg mass formation, and fecundity. Phytopathology. 2020;110(4):927–32.

    Article  CAS  PubMed  Google Scholar 

  11. Shen X, Becelaere GV, Kumar P, Davis RF, May LO, Chee PW. QTL mapping for resistance to root-knot nematodes in the M-120 RNR Upland cotton line (Gossypium hirsutum L.) of the Auburn 623 RNR source. Theoretical Appl Genet. 2006;113(8):1539–49.

    Article  CAS  Google Scholar 

  12. Shen X, He Y, Lubbers EL, Davis RF, Nichols RL, Chee PW. Fine mapping QMi-C11 a major QTL controlling root-knot nematodes resistance in Upland cotton. Theor Appl Genet. 2010;121(8):1623–31.

    Article  PubMed  Google Scholar 

  13. Gutierrez OA, Jenkins JN, McCarty JC, Wubben MJ, Hayes RW, Callahan FE. SSR markers closely associated with genes for resistance to root-knot nematode on chromosomes 11 and 14 of Upland cotton. Theor Appl Genet. 2010;121(7):1323–37.

    Article  CAS  PubMed  Google Scholar 

  14. He Y, Kumar P, Shen X, Davis RF, Van Becelaere G, May OL, Nichols RL, Chee PW. Re-evaluation of the inheritance for root-knot nematode resistance in the Upland cotton germplasm line M-120 RNR revealed two epistatic QTLs conferring resistance. Theor Appl Genet. 2014;127(6):1343–51.

    Article  CAS  PubMed  Google Scholar 

  15. Kumar P, He Y, Singh R, Davis RF, Guo H, Paterson AH, Peterson DG, Shen X, Nichols RL, Chee PW. Fine mapping and identification of candidate genes for a QTL affecting Meloidogyne incognita reproduction in Upland cotton. In, vol. 17; 2016.

  16. Wubben MJ, Thyssen GN, Callahan FE, Fang DD, Deng DD, McCarty JC, Li P, Islam MS, Jenkins JN. A novel variant of Gh_D02G0276 is required for root-knot nematode resistance on chromosome 14 (D02) in Upland cotton. Theoretical Appl Genet. 2019;132(5):1425–34.

    Article  CAS  Google Scholar 

  17. Beneventi MA, da Silva OB Jr., de Sá MEL, Firmino AAP, de Amorim RMS, Albuquerque EVS, da Silva MCM, da Silva JP, Campos MA, Lopes MJC, et al. Transcription profile of soybean-root-knot nematode interaction reveals a key role of phythormones in the resistance reaction. BMC Genomics. 2013;14:322.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Cao K, Li H, Wang Q, Zhao P, Zhu G, Fang W, Chen C, Wang X, Wang L. Comparative transcriptome analysis of genes involved in the response of resistant and susceptible peach cultivars to nematode infection. Sci Hort. 2017;215:20–7.

    Article  CAS  Google Scholar 

  19. Ling J, Mao Z, Zhai M, Zeng F, Yang Y, Xie B. Transcriptome profiling of Cucumis metuliferus infected by Meloidogyne incognita provides new insights into putative defense regulatory network in Cucurbitaceae. In, vol. 7; 2017.

  20. Ye D-Y, Qi Y-H, Cao S-F, Wei B-Q, Zhang H-S. Histopathology combined with transcriptome analyses reveals the mechanism of resistance to Meloidogyne incognita in Cucumis metuliferus. J Plant Physiol. 2017;212:115–24.

    Article  CAS  PubMed  Google Scholar 

  21. Shukla N, Yadav R, Kaur P, Rasmussen S, Goel S, Agarwal M, Jagannath A, Gupta R, Kumar A. Transcriptome analysis of root-knot nematode (Meloidogyne incognita)‐infected tomato (Solanum lycopersicum) roots reveals complex gene expression profiles and metabolic networks of both host and nematode during susceptible and resistance responses. Mol Plant Pathol. 2018;19(3):615–33.

    Article  CAS  PubMed  Google Scholar 

  22. Macharia TN, Bellieny-Rabelo D, Moleleki LN. Transcriptional profiling of potato (Solanum tuberosum L.) during a compatible interaction with the root-knot nematode. Meloidogyne javanica bioRxiv 2019:849414.

  23. Mota APZ, Fernandez D, Arraes FBM, Petitot A-S, de Melo BP, de Sa MEL, Grynberg P, Saraiva MAP, Guimaraes PM, Brasileiro ACM, et al. Evolutionarily conserved plant genes responsive to root-knot nematodes identified by comparative genomics. Mol Genet Genomics. 2020;295(4):1063–78.

    Article  CAS  PubMed  Google Scholar 

  24. Das S, DeMason DA, Ehlers JD, Close TJ, Roberts PA. Histological characterization of root-knot nematode resistance in cowpea and its relation to reactive oxygen species modulation. J Exp Bot. 2008;59(6):1305–13.

    Article  CAS  PubMed  Google Scholar 

  25. Dhandaydham M, Charles L, Zhu H, Starr JL, Huguet T, Cook DR, Prosperi JM, Opperman C. Characterization of Root-Knot Nematode Resistance in Medicago truncatula. In. United States: SOCIETY OF NEMATOLOGISTS; 2008. p. 46.

    Google Scholar 

  26. Kumar P, Khanal S, da Silva M, Singh R, Davis RF, Nichols RL, Chee PW. Transcriptome analysis of a nematode resistant and susceptible upland cotton line at two critical stages of Meloidogyne incognita infection and development. PLoS ONE. 2019;14(9):1–19.

    Article  Google Scholar 

  27. Ojeda-Rivera JO, Ulloa M, Roberts PA, Kottapalli P, Wang C, Nájera-González HR, Payton P, Lopez-Arredondo D, Herrera-Estrella L. Root-knot nematode resistance in Gossypium hirsutum determined by a constitutive defense-response transcriptional program avoiding a fitness penalty. Front Plant Sci. 2022;13:858313. https://doiorg.publicaciones.saludcastillayleon.es/10.3389/fpls.2022.858313.

  28. Zhang T, Hu Y, Jiang W, Fang L, Guan X, Chen J, Zhang J, Saski CA, Scheffler BE, Stelly DM, et al. Sequencing of allotetraploid cotton (Gossypium hirsutum L. acc. TM-1) provides a resource for fiber improvement. Nat Biotechnol. 2015;33(5):531–7.

    Article  CAS  PubMed  Google Scholar 

  29. Yu J, Jung S, Cheng C-H, Ficklin SP, Lee T, Zheng P, Jones D, Percy RG, Main D. CottonGen: a genomics, genetics and breeding database for cotton research. Nucleic Acids Res. 2013;42(D1):D1229–36.

    Article  PubMed  PubMed Central  Google Scholar 

  30. de Deus Barbosa AEA, da Rocha Fragoso R, de Lima e Souza DS, Freire É, de Oliveira Neto OB, Viana AAB, Togawa RC, Guimarães LM, Martins NF, Cia E, et al. Differentially expressed genes in cotton plant genotypes infected with Meloidogyne incognita. Plant Sci. 2009;177(5):492–7.

  31. Xing X, Li X, Zhang M, Wang Y, Liu B, Xi Q, Zhao K, Wu Y, Yang T. Transcriptome analysis of resistant and susceptible tobacco (Nicotiana tabacum) in response to root-knot nematode Meloidogyne incognita infection. In: 2017 2017; United States. Elsevier Science B.V., Amsterdam: 1114.

  32. Li X, Xing X, Tian P, Zhang M, Huo Z, Zhao K, Liu C, Duan D, He W, Yang T. Comparative transcriptome profiling reveals defense-related genes against Meloidogyne incognita invasion in tobacco. Molecules. 2018;23(8).

  33. Postnikova OA, Hult M, Shao J, Skantar A, Nemchinov LG. Transcriptome analysis of resistant and susceptible alfalfa cultivars infected with root-knot nematode Meloidogyne incognita. PLoS One. 2015;10(2):e0118269.

    Article  PubMed  PubMed Central  Google Scholar 

  34. Jones JDG, Dangl JL. The plant immune system. Nature. 2006;444(7117):323–9.

    Article  CAS  PubMed  Google Scholar 

  35. Ali MA, Anjam MS, Nawaz MA, Lam H-M, Chung G. Signal Transduction in plant-nematode interactions. Int J Mol Sci. 2018;19(6):1648.

  36. Goverse A, Smant G. The activation and suppression of plant innate immunity by parasitic nematodes. In. United States: Annual Reviews Phytopathology. 2014;243-65. https://doiorg.publicaciones.saludcastillayleon.es/10.1146/annurev-phyto-102313-050118.

  37. Camera SL, Geoffroy P, Samaha H, Ndiaye A, Rahim G, Legrand M, Heitz T. A pathogen-inducible patatin-like lipid acyl hydrolase facilitates fungal and bacterial host colonization in Arabidopsis. Volume 44. Great Britain: Blackwell Publishing Ltd; 2005. pp. 810–25.

    Google Scholar 

  38. Camera SL, Balague C, Gobel C, Geoffroy P, Legrand M, Feussner I, Roby D, Heitz T. The Arabidopsis Patatin-Like Protein 2 (PLP2) plays an essential role in cell death execution and differentially affects biosynthesis of oxylipins and resistance to pathogens. In. Great Britain: Aps Press. 2009:469.

  39. Tirumalaraju SV, Jain M, Gallo M. Differential gene expression in roots of nematode-resistant and -susceptible peanut (Arachis hypogaea) cultivars in response to early stages of peanut root-knot nematode (Meloidogyne arenaria) parasitization. J Plant Physiol. 2011;168(5):481–92.

    Article  CAS  PubMed  Google Scholar 

  40. Holbein J, Franke RB, MarhavA, Fujita S, Gorecka M, Sobczak M, Geldner N, Schreiber L, Grundler FMW, Siddique S. Root endodermal barrier system contributes to defence against plant-parasitic cyst and root-knot nematodes. Plant J. 2019(2):221.

  41. López-Galiano MJ, González-Hernández AI, Crespo-Salvador O, Rausell C, Real MD, Escamilla M, Camañes G, García-Agustín P, González-Bosch C, García-Robles I. Epigenetic regulation of the expression of WRKY75 transcription factor in response to biotic and abiotic stresses in Solanaceae plants. Plant Cell Rep. 2018;37(1):167–76.

    Article  PubMed  Google Scholar 

  42. Derksen H, Rampitsch C, Daayf F. Signaling cross-talk in plant disease resistance. Plant Sci. 2013;207:79–87.

    Article  CAS  PubMed  Google Scholar 

  43. Ibrahim HMM, Hosseini P, Alkharouf NW, Hussein EHA, Gamal El-Din AEKY, Aly MAM, Matthews BF. Analysis of Gene expression in soybean (Glycine max) roots in response to the root knot nematode Meloidogyne incognita using microarrays and KEGG pathways. BMC Genomics. 2011;12:1–16.

    Article  Google Scholar 

  44. van Loon LC, Geraats BPJ, Linthorst HJM. Ethylene as a modulator of disease resistance in plants. Trends Plant Sci. 2006;11(4):184–91.

    Article  PubMed  Google Scholar 

  45. Cheng HQ, Han LB, Yang CL, Wu XM, Zhong NQ, Wu JH, Wang FX, Wang HY, Xia GX. The cotton MYB108 forms a positive feedback regulation loop with CML11 and participates in the defense response against Verticillium dahliae infection. Journal of Experimental Botany. 2016;67(6):1935-50. https://doiorg.publicaciones.saludcastillayleon.es/10.1093/jxb/erw016.

  46. Barcala M, Garcia A, Cabrera J, Casson S, Lindsey K, Favery B, Garcia-Casado G, Solano R, Fenoll C, Escobar C. Early transcriptomic events in microdissected Arabidopsis nematode-induced giant cells. Plant J. 2010;61(4):698-712. https://doiorg.publicaciones.saludcastillayleon.es/10.1111/j.1365-313X.2009.04098.x.

  47. Santini L, Munhoz CF, Bonfim MF Jr., Brandao MM, Inomoto MM, Carneiro Vieira ML. Host transcriptional profiling at early and later stages of the compatible interaction between Phaseolus vulgaris and Meloidogyne incognita. Nematology. 2016;106:282–94. https://doiorg.publicaciones.saludcastillayleon.es/10.1094/PHYTO-07-15-0160-R.

  48. Zhao Y. Auxin biosynthesis: a simple two-step pathway converts tryptophan to indole-3-acetic acid in plants. Mol Plant. 2012;5(2):334–8.

    Article  CAS  PubMed  Google Scholar 

  49. Martínez-Fernández I, Sanchís S, Marini N, Balanzá V, Ballester P, Navarrete-Gómez M, Oliveira AC, Colombo L, Ferrándiz C. The effect of NGATHA altered activity on auxin signaling pathways within the Arabidopsis gynoecium. Front Plant Sci. 2014;5(210).

  50. Kong Y, Zhu Y, Gao C, She W, Lin W, Chen Y, Han N, Bian H, Zhu M, Wang J. Tissue-specific expression of SMALL AUXIN UP RNA41 differentially regulates cell expansion and root meristem patterning in Arabidopsis. Plant Cell Physiol. 2013;54(4):609–21.

    Article  CAS  PubMed  Google Scholar 

  51. Lee HW, Cho C, Pandey SK, Park Y, Kim M-J, Kim J. LBD16 and LBD18 acting downstream of ARF7 and ARF19 are involved in adventitious root formation in Arabidopsis. BMC Plant Biol. 2019;19(1):46.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Ozalvo R, Cabrera J, Escobar C, Christensen SA, Borrego EJ, Kolomiets MV, Castresana C, Iberkleid I, Brown Horowitz S. Two closely related members of Arabidopsis 13-lipoxygenases (13-LOXs), LOX3 and LOX4, reveal distinct functions in response to plant-parasitic nematode infection. Mol Plant Pathol. 2014(4):319.

  53. Tada Y, Spoel SH, Pajerowska-Mukhtar K, Mou Z, Song J, Wang C, Zuo J, Dong X. Plant immunity requires conformational charges of NPR1 via S-nitrosylation and thioredoxins. Science. 2008;321(5891):952–6.

    Article  CAS  PubMed  Google Scholar 

  54. Linthorst HJ, Melchers LS, Mayer A, van Roekel JS, Cornelissen BJ, Bol JF. Analysis of gene families encoding acidic and basic beta-1,3-glucanases of tobacco. Proc Natl Acad Sci USA. 1990;87(22):8756–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Okinaka Y, Mimori K, Takeo K, Kitamura S, Takeuchi Y, Yamaoka N, Yoshikawa M. A structural model for the mechanisms of elicitor release from fungal cell walls by plant [beta]-1,3-endoglucanase. Plant Physiol. 1995;109(3):839–45.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. Giebel J. Beta-glucosidase activity in potato roots and its possible role in plant tissue response to Heterodera rostochiensis. Bull Acad Pol Sci Biol. 1976;24(1):37–42.

    CAS  PubMed  Google Scholar 

  57. Cabrera J, Díaz-Manzano FE, Fenoll C, Escobar C. Developmental pathways mediated by hormones in nematode feeding sites. Adv Bot Res. 2015;73:167–88. https://doiorg.publicaciones.saludcastillayleon.es/10.1016/bs.abr.2014.12.005.

  58. Zinov’eva SV, Vasyukova NI, Ozeretskovskaya OL. Biochemical aspects of plant interactions with phytoparasitic nematodes: a review. Appl Biochem Microbiol. 2004;40(2):111–9.

    Article  CAS  Google Scholar 

  59. Gortari MC, Hours RA. Fungal chitinases and their biological role in the antagonism onto nematode eggs. A review. In. Netherlands: Springer Science + Business Media; 2008:221.

  60. Qiu J, Hallmann J, Kokalis-Burelle N, Weaver DB, Rodriguez-Kabana R, Tuzun S. Activity and Differential Induction of Chitinase Isozymes in Soybean Cultivars Resistant or Susceptible to Root-knot Nematodes. United States: Society of nematologists; 1997. p. 523.

    Google Scholar 

  61. Rodiuc N, Barlet X, Hok S, Perfus-Barbeoch L, Allasia V, Engler G, Seassau A, Marteu N, Almeida-Engler J, Panabieres F et al. Evolutionarily distant pathogens require the Arabidopsis phytosulfokine signalling pathway to establish disease. Plant Cell Environ. 2016(7):1396.

  62. de Vleesschauwer V, Burssens S, Celenza JL, Inzé D, Van Montagu M, Engler G, Gheysen G. Molecular markers and cell cycle inhibitors show the importance of cell cycle progression in nematode-induced galls and syncytia. Plant Cell. 1999;11(5):793–807.

    Article  PubMed  PubMed Central  Google Scholar 

  63. Wuyts N, Lognay G, Swennen R, Waele DD. Nematode infection and reproduction in transgenic and mutant Arabidopsis and tobacco with an altered phenylpropanoid metabolism. J Exp Bot. 2006;57(11):2825.

    Article  CAS  PubMed  Google Scholar 

  64. Libault M, Wan J, Czechowski T, Udvardi M, Stacey G. Identification of 118 Arabidopsis Transcription Factor and 30 Ubiquitin-Ligase Genes Responding to Chitin, a Plant-Defense Elicitor. Mol Plant Microbe Interact. 2007;20(8):900–11. Great Britain: APS Press.

    Article  CAS  PubMed  Google Scholar 

  65. Zhou E, Star JL, Smith CW. Inheritance of resistance to Meloidogyne incognita in the cotton cultivar Acala NemX. J Nematol. 1999;31:584–5.

    Google Scholar 

  66. Congli W, Mauricio U, Tra TD, Philip AR. QTL analysis of transgressive nematode resistance in tetraploid cotton reveals complex interactions in chromosome 11 regions. Front Plant Sci. 2017;8:1979. https://doiorg.publicaciones.saludcastillayleon.es/10.3389/fpls.2017.01979.

  67. Damiani I, Baldacci-Cresp F, Hopkins J, Andrio E, Balzergue S, Lecomte P, Puppo A, Abad P, Favery B, Hérouart D. Plant genes involved in harbouring symbiotic rhizobia or pathogenic nematodes. New Phytol. 2012;194(2):511–22.

    Article  CAS  PubMed  Google Scholar 

  68. Shepherd RL, McCarty JC Jr., Jenkins J, Parrott WL. Registration of nine cotton germplasm resistant to root-knot nematode. Crop Sci. 1996;36:820.

    Article  Google Scholar 

  69. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30(15):2114–20. https://doiorg.publicaciones.saludcastillayleon.es/10.1093/bioinformatics/btu170.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  70. Kim D, Paggi JM, Park C, Bennett C, Salzberg SL. Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nat Biotechnol. 2019;37(8):907–15.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  71. Pertea M, Pertea GM, Antonescu CM, Chang T-C, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. 2015;33(3):290–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  72. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.

    Article  PubMed  PubMed Central  Google Scholar 

  73. Varet H, Brillet-Guéguen L, Coppée J-Y, Dillies M-A. SARTools: a DESeq2- and EdgeR-based R pipeline for comprehensive differential analysis of RNA-Seq data. PLoS One. 2016;11(6):1–8.

    Article  Google Scholar 

  74. Ryan MT, Stucky M, Wakefield C, Melott JM, Akbani R, Weinstein J, Broom N. BM: Interactive Clustered Heat Map Builder: an easy web-based tool for creating sophisticated clustered heat maps. F1000Res. 2019;8:1750. https://doiorg.publicaciones.saludcastillayleon.es/10.12688/f1000research.20590.2.

  75. Zhu T, Liang C, Meng Z, Sun G, Meng Z, Guo S, Zhang R. CottonFGD: an integrated functional genomics database for cotton. BMC Plant Biol. 2017;17(1):101.

    Article  PubMed  PubMed Central  Google Scholar 

  76. Pathan M, Keerthikumar S, Ang C-S, Gangoda L, Quek CYJ, Williamson NA, Mouradov D, Sieber OM, Simpson RJ, Salim A et al. FunRich: an open access standalone functional enrichment and interaction network analysis tool. Proteomics. 2015(15):2597.

  77. Cabrera J, Bustos R, Favery B, Fenoll C, Escobar C. NEMATIC: a simple and versatile tool for the in silico analysis of plant-nematode interactions. Mol Plant Pathol. 2014;15(6):627–36.

    Article  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

We extend our heartfelt acknowledgment to the late Dr. Robert L. Nichols (1946-2020) for his invaluable supervision and profound contributions to our root knot nematode cotton manuscript. We would also like to thank Edward Lubbers and Jennifer McBlanchett for their technical help in designing and managing field experiments.

Funding

This project was funded in part by Cotton Incorporated and the Georgia Cotton Commission.

Author information

Authors and Affiliations

  1. Department of Crop and Soil Sciences and Institute of Plant Breeding, Genetics, and Genomics, University of Georgia, 2356 Rainwater Road, Tifton, GA, 31793, USA

    Sameer Khanal, Pawan Kumar, Mychele B. da Silva, Rippy Singh, Nelson Suassuna & Peng W. Chee

  2. Bayer Crop Science Division, 700 W. Chesterfield Pkwy W, Chesterfield, MO, 63017, USA

    Pawan Kumar & Rippy Singh

  3. Tidewater AREC, Virginia Tech, 6321 Holland Road, Suffolk, VA, 23437, USA

    Mychele B. da Silva

  4. Empresa Brasileira de Pesquisa Agropecuária (Embrapa), Santo Antônio de Goiás, GO, Brazil

    Nelson Suassuna

  5. Cotton Incorporated, 6399 Weston Parkway, Cary, NC, 27513, USA

    Don C. Jones

  6. U.S. Department of Agriculture – Agricultural Research Service, 115 Coastal Way, Tifton, GA, 31793, USA

    Richard F. Davis

Authors
  1. Sameer Khanal
    View author publications

    You can also search for this author inPubMed Google Scholar

  2. Pawan Kumar
    View author publications

    You can also search for this author inPubMed Google Scholar

  3. Mychele B. da Silva
    View author publications

    You can also search for this author inPubMed Google Scholar

  4. Rippy Singh
    View author publications

    You can also search for this author inPubMed Google Scholar

  5. Nelson Suassuna
    View author publications

    You can also search for this author inPubMed Google Scholar

  6. Don C. Jones
    View author publications

    You can also search for this author inPubMed Google Scholar

  7. Richard F. Davis
    View author publications

    You can also search for this author inPubMed Google Scholar

  8. Peng W. Chee
    View author publications

    You can also search for this author inPubMed Google Scholar

Contributions

PWC and RFD conceived and designed the experiments. SK, PK, MDS, RS, NS conducted the experiments. SK and PK analyzed the data. SK wrote the paper. DCJ provided supervisory support on behalf of the funding agency. All the authors have reviewed the manuscript.

Corresponding authors

Correspondence to Sameer Khanal or Peng W. Chee.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

12864_2025_11339_MOESM1_ESM.xlsx

Supplementary Material 1: Supplemental Table S1. Differentially expressed genes (DEGs) during compatible and incompatible interactions between G. hirsutum-M. incognita interactions.

Supplementary Material 2: Supplemental Table S2. Stress- and defense-related DEGs.

12864_2025_11339_MOESM3_ESM.xlsx

Supplementary Material 3: Supplemental Table S3. Significant DEGs identified during compatible and incompatible G. hirsutum-M. incognita interactions and their putative nematode-responsive homologous genes in Arabidopsis cataloged in NEMATIC.

12864_2025_11339_MOESM4_ESM.xlsx

Supplementary Material 4: Supplemental Table S4. Significant DEGs identified during compatible and incompatible G. hirsutum-M. incognita interactions and their putative nematode-responsive orthologous genes in Arabidopsis cataloged in 'stress' category in NEMATIC.

Supplementary Material 5: Supplemental Table S5. Cell wall biogenesis and remodeling related DEGs.

Supplementary Material 6: Supplemental Table S6. Cell growth and development-related DEGs.

Supplementary Material 7: Supplemental Table S7. Differentially expressed transcription factors (TFs).

Supplementary Material 8: Supplemental Table S8. DEGs related to carbon and energy metabolism.

12864_2025_11339_MOESM9_ESM.xlsx

Supplementary Material 9: Supplemental Table S9. Detailed information of annotated M. incognita genes identified in RNA-seq libraries and includes RNAi pathway proteins, carbohydrate active enzymes (CAZymes) and nematode effectors.

12864_2025_11339_MOESM10_ESM.xlsx

Supplementary Material 10: Supplemental Table S10. Differentially expressed genes (DEGs) in contrasts corresponding to compatible and incompatible interactions between G. hirsutum-M. incognita interactions.

12864_2025_11339_MOESM11_ESM.tiff

Supplementary Material 11: Supplemental Figure S1. Distribution of mapped and unmapped paired-end sequences corresponding to two different genotypes (‘C’ C201 and ‘M’ M120), six different time points (‘0’, ‘4’, ‘8’, ‘12’, ‘16’ and ‘20’), and two different treatments (‘Cont’ control and ‘RKN’ treated with RKN).

12864_2025_11339_MOESM12_ESM.tiff

Supplementary Material 12: Supplemental Figure S2. Number of common DEGs during compatible and incompatible interactions at different DAI.

Diagrams show a. Venn diagram of upregulated genes (‘red font’ are total numbers for the group), b. Venn diagram of downregulated genes, and c. frequencies and percentages of upregulated genes.

12864_2025_11339_MOESM13_ESM.tiff

Supplementary Material 13: Supplemental Figure S3. Quantitative real-time PCR (qRT-PCR) of candidate root-knot nematode (RKN) resistance genes at six different infection time points corresponding to two different genotypes (‘C’ C201 and ‘M’ M120), six different time points (‘4’, ‘8’, ‘12’, ‘16’, ‘20’, and ’30’ days after infection, DAI), and two different treatments control and treated with RKN).

β- actin gene was used as reference gene (endogenous control) for normalizing the data. Charts show mean fold-change (with standard error bars) in relative expression of target genes against the endogenous control after RKN infection. Green asterisks mark significant upregulation, and red asterisks mark significant down regulation in expression in inoculated samples compared to non-inoculated plants and as determined by t-test of ΔCt values (P ≤ 0.05) using two biological replicates each with three technical replicates per sample.

12864_2025_11339_MOESM14_ESM.tiff

Supplementary Material 14: Supplemental Figure S4. Expression profiles of the genes in alpha-linolenic acid pathway during compatible and incompatible interactions.

The heat map (a) shows DEGs from C201 (compatible) and M120 (incompatible) samples at five different time points (4, 8, 12, 16 and 20). Darker shades of blue and red bars indicate greater magnitude of down- or upregulation of DEGs, whereas black suggests indifferentiable expression. The pathway (b) shows upregulated (red highlighted) and downregulated (blue highlighted) genes during incompatible interaction at 8 DAI.

12864_2025_11339_MOESM15_ESM.tiff

Supplementary Material 15: Supplemental Figure S5. Expression profiles of the genes encoding enzymes involved in cell wall biogenesis represented by xyloglucan endotransglucosylase/hydrolase protein (a), cell wall macromolecule catabolic process dominated by chitinase and endochitinases (b), cell wall modification represented by pectinestarase (c), and cell wall organization represented by expansins (d) during compatible and incompatible interactions between G. hirsutum and M. incognita at five different time points. Darker shades of blue and red bars indicate greater magnitude of down- or upregulation of DEGs, whereas black suggests indifferentiable expression.

12864_2025_11339_MOESM16_ESM.tiff

Supplementary Material 16. Supplemental Figure S6. Expression profiles of the genes encoding enzymes involved in phenylpropanoid biosynthetic pathway during compatible and incompatible interactions between G. hirsutum and M. incognita at five different time points (a), and phenylpropanoid biosynthesis pathway at 8 DAI in the incompatible interaction (b). Heatmap: darker shades of blue and red bars indicate greater magnitude of down- or upregulation of DEGs, whereas black suggests indifferentiable expression. The pathway (b) shows upregulated (red highlighted) and downregulated (blue highlighted) genes during incompatible interaction at 8 DAI.

12864_2025_11339_MOESM17_ESM.tiff

Supplementary Material 17: Supplemental Figure S7. Expression profiles of the genes encoding enzymes involved in amino sugar nucleosugar metabolism (a), glycolysis/gluconeogenesis (b), and starch and sucrose metabolism (c) during compatible and incompatible interactions between G. hirsutum and M. incognita at five different time points. Darker shades of blue and red bars indicate greater magnitude of down- or upregulation of DEGs, whereas blacksuggests indifferentiable expression.

12864_2025_11339_MOESM18_ESM.tiff

Supplementary Material 18: Supplemental Figure S8. Expression profiles of the genes encoding enzymes in amino sugar and nucleotide sugar pathway at 8 (‘left’) and 16 (‘right’) DAI in an incompatible interaction between G. hirsutum and M. incognita. The pathways show upregulated (red highlighted) and downregulated (blue highlighted) genes.

12864_2025_11339_MOESM19_ESM.tiff

Supplementary Material 19: Supplemental Figure S9. Expression profiles of the genes encoding enzymes in glycolysis/gluconeogenesis pathway at 8 DAI in a compatible interaction (‘left’) and 8 and 16 DAI (‘center’ and ‘right’) in an incompatible interaction between G. hirsutum and M. incognita. The pathways show upregulated (red highlighted) and downregulated (blue highlighted) genes.

12864_2025_11339_MOESM20_ESM.tiff

Supplementary Material 20: Supplemental Figure S10. Expression profiles of the genes encoding enzymes in starch and sucrose metabolism pathway at 8 (‘left’) and 16 (‘right’) DAI in an incompatible interaction between G. hirsutum and M. incognita. The pathways show upregulated (red highlighted) and downregulated (blue highlighted) genes.

12864_2025_11339_MOESM21_ESM.xlsx

Supplementary Material 21: Supplemental Table S11. Number of Meloidogyne incognita (different developmental stages) in the roots of cotton near isolines that differ in M. incognita resistance genes.

12864_2025_11339_MOESM22_ESM.xlsx

Supplementary Material 22: Supplemental Table S12. Percentages of Meloidogyne incognita (different developmental stages) in the roots of cotton near isolines that differ in M. incognita resistance genes.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Khanal, S., Kumar, P., da Silva, M.B. et al. Time-course RNA-seq analysis of upland cotton (Gossypium hirsutum L.) responses to Southern root-knot nematode (Meloidogyne incognita) during compatible and incompatible interactions. BMC Genomics 26, 183 (2025). https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12864-025-11339-w

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12864-025-11339-w

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us