A reniform nematode (Rotylenchulus reniformis) culture was maintained on the susceptible cotton line M8 in a growth chamber. Eggs used for inoculum in experiments were collected from infected M8 roots by sodium hypochlorite washing according to the method of Hussey and Barker (1973). Nearly-isogenic lines (NILs) fully susceptible to R. reniformis or having the GB713-derived Ren^barb2 resistance QTL had previously been developed and characterized by our research group (Wubben et al., 2017).

Seeds of resistant and susceptible NILs were scarified, imbibed in water at 30°C for four hours, and then germinated in paper towels overnight at 30°C. Germinated seeds having ~ 0.5 cm radicles were sown into 20 cm Cone-tainers containing a 1:7 mix of autoclaved silty loam soil:sand. Inoculated Cone-tainers received ~ 5,000 R. reniformis eggs one-day after planting whereas nothing was added to control plants. Cone-tainers were arranged in rows of three within racks such that one row equaled a single biological replicate. Upon root tissue harvest, the root systems of three plants within each row were pooled for RNA extraction. Three biological replicates of each NIL (resistant, susceptible) × timepoint (5-, 9-, 13-days after inoculation) × treatment (control, inoculated) were arranged in racks in a completely randomized design. Plants were grown in a Percival PGC-9/2 growth chamber at 30°C on a 16 hr. day/8 hr. night schedule. Extra root samples were taken at each timepoint and stained with acid fuchsin to visualize the extent of R. reniformis infection and stages of development. Root staining with acid fuchsin was accomplished according to Wubben et al. (2020). Pictures of infected roots were taken with a digital camera mounted on a Nikon SMZ1500 stereomicroscope.

Root tissues were washed free of soil, blotted dry, wrapped in aluminum pouches, flash-frozen in liquid nitrogen, and stored at -80°C until RNA extraction. Total RNA was isolated by grinding the root samples under liquid nitrogen with a mortar and pestle and extracting the RNA via the hot borate method as previously described (Wubben et al., 2019). Purity and concentration of total RNA was determined using a NanoDrop 1000 spectrophotometer. A total of (24) RNA samples representing two of the three biological replicates each of NIL × timepoint × treatment was submitted to the Georgia Genomic Facility (University of Georgia). Sequencing libraries were made using the Kapa Stranded RNA-Seq Kit (Roche Inc., Indianapolis, IN) and 150-bp paired-end reads were generated using a NextSeq PE150 High Output flow cell platform (Illumina Inc., San Diego, CA). Sequence quality was assessed using FastQC (Andrews, 2010). Low-quality bases and adapter sequences were trimmed from paired reads using Trimmomatic v0.30 (Bolger et al., 2014).

The Gossypium hirsutum (AD1) Genome-Texas Interim release UTX-JGI v1.1 (available via www.cottongen.org) was used for sequence read alignment. Cotton transcriptome assembly was accomplished using the Galaxy web platform via public server at usegalaxy.org (Afgan et al., 2016) using the following built-in functions: sequence alignment was done using HISAT2 v2.1.0 (Kim et al., 2015) and transcriptome assembly using Stringtie v1.3.3 (Pertea et al., 2015). A total of 201,260,466 (~200 million) paired-end reads were generated, of which 159,566,283 (~160 million; ~79%) were concordantly mapped (using HISAT2) to the reference genome ( Supplementary Figure S1 , Supplementary Table S1 ). Based on overall alignment (concordant, discordant, or one of a mate-pair match), approximately 91% of sequences aligned to the reference genome. All catalogued exons (396,608) and loci (66,522) in the reference assembly/annotation were represented in the unified transcriptome (high sensitivity) reconstructed from 24 transcriptome assemblies (using StringTie’s merge), while 179,260 (27.2%) exons corresponding to 89,797 (56.9%) loci were specific to this dataset.

Transcriptome analysis was done using GffCompare (https://github.com/gpertea/gffcompare), gene count using htSeq-count (Anders et al., 2015), and differential expression (DEG) using DESeq2 (Love et al., 2014). For DEG analysis, an adjusted P-value (FDR) of ≤ 0.05 was calculated using the approach of Benjamini and Hochberg (1995). A significant fold-change (FC) in transcript expression was determined at 2-FC (log₂FC > 1.0 or < -1.0). DEG contrasts were made between infected and control at each timepoint for resistant and susceptible cotton NILs. Gene ontology (GO) and KEGG pathway enrichment of genes was accomplished using the data fetch and enrichment tool at CottonFGD (cottonfgd.org) (Zhu et al., 2017). Hierarchical cluster analysis of DEGs was done using the ‘Expression’ tool at www.heatmapper.ca/expression with the settings average linkage for clustering method and Euclidean distance measurement method (Babicki et al., 2016).

Total RNA used for transcriptome sequencing was also used for qRT-PCR. Approximately 1 µg of total RNA was used for genomic DNA digestion and cDNA synthesis using the iScript^Tm gDNA Clear cDNA Synthesis Kit (Bio-Rad, Carlsbad, CA). cDNA reactions were diluted 1:10 to use as template in qRT-PCR reactions. Reactions were performed on a Bio-RAD CFX^TM Opus 96 Real-Time PCR Detection System. Reaction volumes were 10 µl and comprised as follows: 5.0 µl 2X iQTM SYBR^® Green Supermix (Bio-Rad), 0.25 µM primer, 1 µl cDNA template, and deionized water. Three technical reps of each primer × template combination were performed. GhUBQ14 Ct values were used for expression normalization. Two-sample t tests (P ≤ 0.05) were performed on the average normalized ΔCt value of three biological replicates of inoculated versus control samples. Mean fold-change in expression was converted to 2^-ΔΔCt for data presentation. Primers used in qRT-PCR experiments are provided in Supplementary Table S2 .

Resistant and susceptible NILs showed clear differences in R. reniformis infection levels and rate of female development over the course of the experiment and in particular at the 9-dai and 13-dai timepoints. At 5-dai, both resistant and susceptible NILs showed roughly equivalent levels of infection as evidenced by the presence of vermiform sedentary female R. reniformis protruding from the root epidermis ( Figure 1A ). By 9-dai, sedentary females had attained their characteristic kidney-shape; however, the susceptible line showed a much greater number of females compared to resistant roots, and in many instances the females on susceptible roots had begun to produce egg masses, whereas female development on resistant roots lagged behind ( Figure 1A ). By 13-dai, sedentary females on resistant and susceptible root systems had formed egg masses; however, far fewer females were observed on resistant roots versus susceptible ( Figure 1A ). These observations indicate an early resistance reaction occurring shortly after infection in the resistant NIL roots that appears to inhibit feeding site establishment and female development.

To obtain a global view of gene expression changes in resistant and susceptible plants during R. reniformis infection, RNA-Seq was performed on control and infected roots at 5-, 9-, and 13-dai. A total of 1,099 unique DEGs were identified with the vast majority of DEGs coming from resistant plant roots ( Figure 1B ). In resistant roots, up- and down-regulated DEGs increased over time with total DEGs being 129, 210, and 760 at 5-, 9-, and 13-dai, respectively, with up-regulated DEGs always being in the majority. In contrast, in susceptible roots, only 44, 10, and 87 DEGs were detected at the same timepoints with up- and down-regulated DEGs being split near equally at 13-dai ( Figure 1B ). These findings demonstrate a significantly stronger transcriptional response, both positively and negatively, occurred in resistant roots versus susceptible following R. reniformis infection.

The majority of DEGs in resistant and susceptible roots were restricted to single timepoints. Among up-regulated DEGs in resistant roots, 87/625 appeared in more than one timepoint and only 16 genes were up-regulated across all timepoints ( Figure 2A ). Similarly, among down-regulated DEGs in resistant roots, 19/341 appeared in more than one timepoint; however, no DEGs were down-regulated across all timepoints ( Figure 2A ). In susceptible roots, 5/81 up-regulated DEGs appeared in more than one timepoint and only a single gene was up-regulated across all timepoints ( Figure 2B ). This gene also happened to be up-regulated across all timepoints in resistant roots, i.e., Gohir.D10G155200, and encodes a putative early nodulin-75-like protein ( Table 1 ). A total of 52 DEGs were down-regulated in susceptible roots with only a single gene appearing at more than one timepoint ( Figure 2B ). A comparison between up- and down-regulated DEGs between resistant and susceptible roots identified some commonalities in gene expression. For example, approximately 63% of susceptible up-regulated DEGs were also found to be up-regulated in resistant roots for at least a single timepoint ( Figure 2C ). Likewise, approximately 22% of susceptible down-regulated DEGs were also down-regulated in resistant roots ( Figure 2C ). In contrast, a group of 10 genes was identified that showed up-regulation in resistant roots but down-regulation in the susceptible line ( Figure 2C ). Among these 10 genes were multiple defense-related proteins, including laccases, nematode-resistance protein-like, ERF5-like protein, phosphatase 2C, and phospholipase. Sixteen (16) genes were induced in resistant roots across all three timepoints ( Table 1 ). In many instances, genes were also induced in the susceptible genotype at one or more timepoints but not all three with the exception of two early-nodulin-75-like homologs that were induced at all timepoints in both lines. Many of the genes were defense-related such as chitinase, WRKY transcription factor, lipoxygenase, and MIC3. MIC3 was induced in resistant and susceptible plants at 5-DAI, but induction only continued in resistant plants through 9- and 13-dai. A similar pattern of expression between resistant and susceptible roots was observed for chitinase, WRKY5, lipoxygenase, and two uncharacterized proteins ( Table 1 ). These findings suggest that an initial response to R. reniformis infection occurred in the roots of the susceptible NIL but there was no sustained resistance response as in the resistant root tissues.

Hierarchical clustering of all resistant and susceptible DEGs across all timepoints was performed to provide a comprehensive perspective of DEG expression patterns ( Figure 3 ). In resistant roots, two predominant DEG expression patterns were apparent: (i) genes being induced early in the R. reniformis infection process at 5- and 9-dai but then returning to baseline expression levels or being down-regulated by 13-dai and (ii) genes staying at baseline levels early in the infection process but then being dramatically up-regulated at 13-dai ( Figure 3A ). Interestingly, these two trends were observed in susceptible roots, only with far fewer genes compared to resistant plants ( Figure 3B ). These observations indicate a biphasic response in resistant roots, i.e., an ‘early’ and ‘late’ response, that happens to be much stronger and with more genes in resistant versus susceptible plants.

GO enrichment analysis was performed on DEGs from resistant and susceptible roots to identify biological processes that were overrepresented at each timepoint. Due to the low number of DEGs in susceptible roots, only a handful of GO annotations were identified. At 5-dai, oxidation-reduction process (GO:0055114) and carbohydrate metabolism (GO:0005975) were enriched in susceptible roots with 10 and 3 genes, respectively, that were up-regulated ( Table 2 ). These genes would be part of the early response clade in Figure 3B with representatives including peroxidase, lipoxygenase, and other oxygenases with varying cellular functions.

Substantially more GO annotations were enriched in resistant versus susceptible roots ( Table 2 ). Biological processes related to oxidation-reduction and cellular redox homeostasis were heavily represented in resistant roots at 5-dai. DEGs with roles in oxidation-reduction continued to be up-regulated at 9-dai; however, down-regulation of this process also began to occur at this timepoint. By 13-dai, oxidation-reduction DEGs were wholly down-regulated in resistant roots. In contrast, DEGs associated with regulation of transcription (GO:0006355) and various protein secondary modifications, e.g., ubiquitination (GO:0016567), dephosphorylation (GO:0006470), and glycosylation (GO:0006487), were absent at early timepoints but strongly up-regulated at 13-dai. Genes involved in carbohydrate metabolism (GO:0005975) were up-regulated at each time-point with some down-regulation at 9-dai. Processes associated with cell wall fortification such as lignin catabolism (GO:0046274) were up-regulated early at 5-dai but then down-regulated by 9-dai (cell wall biogenesis; GO:0043546) and 13-dai (cell wall organization; GO:0009664). Enrichment by molecular function showed an overabundance of genes involved in DNA binding and transcription factor activity in resistant roots ( Supplementary Figure S2 ). Genes with functions in oxidoreductase activity or other catalytic activity were also heavily represented in resistant roots ( Supplementary Figure S2 ). Enrichment by cellular component further highlighted the resistance/defense nature of gene activity in resistant roots in that cellular spaces ‘outside’ of the cell or related to signaling or secretory processes were overrepresented such as the plasma membrane, apoplast, exocyst, and extracellular space ( Supplementary Figure S2 ). In susceptible roots, only the apoplast was enriched, likely representing an early, yet unsustained, response to R. reniformis infection.

The number of KEGG-enriched pathways increased over time in resistant roots with metabolic pathways (ko1100) and biosynthesis of secondary metabolites (ko1110) predominating at each individual timepoint ( Supplementary Table S3 ). Host plant resistance to nematodes many times involves the activity of compounds that have nematicidal or nemastatic activity and are the products of secondary metabolic pathways (Sato et al., 2019). Phenylpropanoid biosynthesis (ko00940) and plant hormone signal transduction (koko04075) became apparent at 9-dai and remained through 13-dai. At 13-dai the aforementioned pathways were joined by a host of others including ethylene biosynthesis (M00368), plant-pathogen interaction (ko04626), and nitrogen metabolism (ko00910), indicating a large change in transcriptional regulation of multiple signaling and biosynthetic pathways in resistant roots. In contrast, susceptible roots were enriched only at 13-dai largely categorized as metabolic pathways and secondary metabolite biosynthesis which likely reflects the nutrient-sink nature of the syncytium in susceptible plant roots and indicate a compatible interaction with the nematode.

Ninety-seven (97) unique genes having a role in GO biological processes oxidation-reduction (GO:0055114) and response to oxidative stress (GO:0006979) were differentially expressed in resistant and susceptible roots following R. reniformis infection ( Supplementary Table S4 ). This group of DEGs was largely represented, in decreasing order of abundance, by peroxidases, 2-oxoglutarate-dependent dioxygenases, laccases, Berberine bridge enzymes, ACC oxidases, and cytochrome P450 enzymes ( Supplementary Table S4 ). Two distinct clades/expression profiles emerged when hierarchical clustering of these DEG profiles is presented as a heat map across all timepoints for resistant and susceptible plants ( Figure 4A ). In Clade I, there is a general up-regulation or no change in expression across timepoints in resistant plants but an overall down-regulation in susceptible roots. In contrast, in Clade II, there is a clear trend of significant down-regulation of gene expression at 13-dai in resistant roots versus susceptible where gene expression remains largely unchanged from baseline levels. A closer examination of this 13-dai down-regulated group in resistant roots shows it is partially comprised of multiple peroxidases ( Figure 4B , Supplementary Table S4 ). A total of 14 peroxidase genes were differentially expressed in resistant roots with more than half showing significant down-regulation at 13-dai. The remaining peroxidase genes show an inverse expression pattern, being up-regulated at various points throughout infection ( Figure 4B ).

A total of 84 genes representing multiple transcription factor (TF) families were differentially expressed following R. reniformis infection with the overwhelming majority being in resistant plants ( Supplementary Table S5 ). Ethylene-responsive factors and WRKY domain-containing genes comprised more than half of all TFs identified ( Figure 5A ). NAC, auxin-related, and homeobox-related TFs followed in abundance. Numbers of genes corresponding to the ERF, WRKY, and NAC TF families increased steadily over the course of R. reniformis infection in resistant roots ( Figure 5B ).

Ethylene-responsive transcription factors (ERFs) are involved in multiple signal transduction pathways such as hormone signaling (ethylene, auxin, JA) and different aspects of the cotton defense response (Zafar et al., 2022). Twenty-one (21) genes with significant identity to various ERFs were differentially expressed in resistant roots upon R. reniformis infection ( Table 3 ). Gohir.D03G166600 was the single ERF differentially expressed in susceptible plants being repressed by > 2-fold at 5-dai, while in resistant roots this gene was induced by 4-fold at 13-dai. The remaining 20 ERFs were all induced in resistant plants at 13-dai with the exception of Gohir.A08G189900 which was induced at 9-dai. Cotton genes similar to ERF17 family members were induced at 8-16-fold in resistant plants by 13-dai ( Table 3 ). Likewise, Gohir.A13G040200, a homolog of ERF109, was induced by ~16-fold in resistant plants at 13-dai. ERF17 and ERF109 play roles in plant defense signaling including oxidative stress responses and programmed cell death ( Table 3 ).

WRKY transcription factors represented the second major TF class with (17) members being differentially expressed in resistant roots upon R. reniformis infection ( Table 4 ). WRKY genes are transcriptional regulators that are involved in activating and propagating PTI and ETI signaling pathways as well as in response to various abiotic stressors (Javed and Gao, 2023). Similar to the ERF transcription factor class, most of the WRKY genes identified here were induced at a single later timepoint with only four genes being differentially regulated at more than one timepoint. GhWRKY5 was the only gene induced across all timepoints and also the only WRKY gene induced in susceptible plants but only at 5-dai. Three GhWRKY70-like genes were identified with Gohir.D02G003400 being the most up-regulated of all WRKYs discovered at ~ 5-fold. GhWRKY70 has been shown to regulate jasmonic acid production and promote resistance to verticillium wilt (Zhang et al., 2023). Two genes sharing identity with GhWRKY40 were identified that showed induction at 9- and 13-dai. GhWRKY40 has been shown to be regulated by salicylic acid and other hormones while mediating wound- and pathogen-induced responses (Wang et al., 2014). Another regulator of JA signaling, GhWRKY33 (Ji et al., 2023), was identified in our analysis as being induced more than 4-fold at 13-dai.

A number of genes with homology to known defense-related pathways were up-regulated in resistant roots versus susceptible. Cotton genes representative of PR-1, basic and acid endochitinases, thaumatin-like proteins, and MIC-3 were induced across all three timepoints. Multiple genes with similarity to the nematode resistance protein-like HSPRO2 from Arabidopsis thaliana were up-regulated in resistant roots and, in some cases, down-regulated in susceptible roots. Gohir.A05G303900 and Gohir.A10G233600 were both induced 2-3-fold in the resistant line at 13-dai but repressed more than 2-fold in the susceptible at the same timepoint. Interestingly, the likely homoeologous copies of these genes, Gohir.D05G302900 and Gohir.D10G245900 were also induced ~ 3-fold in resistant plants but were not differentially regulated in susceptible roots.

Nine (9) defense-related genes were selected for qRT-PCR to check the validity/accuracy of the RNA-Seq data analyses. The cotton genes GhPR1 and GhMIC3 were both induced across all timepoints in resistant plant roots while both genes were induced only at 5-dai in susceptible roots ( Figure 6 ). These results are consistent with the RNA-Seq analysis with the exception that GhPR1 induction was not detected at such a high level as by qRT-PCR. This likely reflects the difference in sensitivity between the assays for this particular transcript. qRT-PCR showed that the WRKY genes GhWRKY50 and GhWRKY92 were both induced at 13-dai in resistant roots but not differentially expressed at all in the susceptible line ( Figure 6 ). Likewise, the ERF genes GhERF109 and GhERF17 showed dramatic induction in resistant roots at 13-dai, whereas no significant change in expression was observed in susceptible roots ( Figure 6 ). The WRKY gene GhWRKY70 was induced at all timepoints, but significantly so at 5- and 9-dai ( Figure 6 ). This expression profile correlates well with that determined by RNA-Seq ( Table 4 ). The peroxidase gene GhPER52 was also assessed by qRT-PCR. RNA-Seq showed early induction of this gene in resistant roots and somewhat later induction in the susceptible ( Figure 4B ), and in this assay, induction of GhPER52 was at significant levels in the resistant line at 5-dai ( Figure 6 ). GhPER52 was significantly induced in susceptible roots at 9-dai. Overall, the qRT-PCR expression profiles of the selected genes mirrored those provided by RNA-Seq analysis.

The Ren^barb2 QTL is defined by the SSR markers BNL3279 and BNL4011 on chromosome 21 (D11) (Gutierrez et al., 2011). Using the TX-JGIv1.1 genome sequence, these markers delineated a ~ 1.2 Mbp region, bp 65886069 – 67094220, that contained 59 predicted genes. Using all of the available RNA-Seq data from the resistant root samples, we detected 25 transcripts that contained SNPs compared to the TX-JGIv1.1 genome ( Supplementary Table S6 ). Furthermore, we were able to identify non-synonymous SNPs specific to the resistant line that occurred within the coding region of 16/25 genes whose expression we detected. We were also able to determine that the non-synonymous SNPs for 3/16 genes were shared within the genomic sequence of the B713 cotton line that also contains the Ren^barb2 QTL (Perkin et al., 2021). The non-synonymous SNPs within the remaining 13/16 genes were found to be present within the R. reniformis-susceptible G. barbadense genome sequence 3-79 HAUv2 and were dropped from further study. The gene IDs and SNP position data for the remaining candidate Ren^barb2 resistance genes, based on our RNA-Seq SNP analysis of resistant roots, are provided in Table 5 .

Gohir.D11G301700, having similarity to the Lr10 rust resistance kinase, showed a single non-synonymous SNP that results in a F64L change within the predicted galacturonan-binding region ( Table 5 ). A single SNP was also identified within Gohir.D11G306000, a cold-responsive protein kinase homolog, that causes a G284R change within the predicted Ser/Thr kinase catalytic domain. The third candidate gene, Gohir.D11G302300, a predicted CC-NBS-LRR resistance protein, contained six non-synonymous SNPs, with four of the SNPs causing mutations within the predicted NB-ARC domain of the protein.

In addition to the genes listed in Table 5 , gene Gohir.D11G304100, which encodes a probable LRR-RLK, contained a non-synonymous SNP absent in the susceptible 3-79 HAUv2 genome sequence; however, the SNP was also absent from the resistant B713 genome sequence. Gene Gohir.D11G304600, which lacks homology to any known protein, contained two non-synonymous SNPs but extensive BLAST analyses failed to identify a G. barbadense homolog (data not shown).

We did not detect differential expression of any of the resistance candidate genes in the RNA-Seq data; therefore, qRT-PCR was conducted to detect possible changes in expression. No significant change in transcript levels was observed for any of the candidate genes or for genes D11G304100 and D11G304600 in response to R. reniformis infection (data not shown). However, significant differences in baseline expression between resistant and susceptible control roots for some of the candidates was observed ( Figure 7 ). D11G304100 and D11G304600 both showed slight but significant down-regulation in resistant roots compared to susceptible plants; however, D11G302300, which encodes a putative CC-NBS-LRR resistance protein, showed a nearly 4-fold increase in expression in control resistant roots versus control susceptible roots ( Figure 7 ).