Transcriptome Analysis of Resistant Cotton Germplasm Responding to Reniform Nematodes

Reniform nematode (Rotylenchulus reniformis) is an important microparasite for Upland cotton (Gossypium hirsutum L.) production. Growing resistant cultivars is the most economical management method, but only a few G. barbadense genotypes and some diploid Gossypium species confer high levels of resistance. This study conducted a transcriptome analysis of resistant genotypes to identify genes involved in host plant defense. Seedlings of G. arboreum accessions PI 529728 (A2-100) and PI 615699 (A2-190), and G. barbadense genotypes PI 608139 (GB 713) and PI 163608 (TX 110), were inoculated with the reniform nematode population MSRR04 and root samples were collected on the fifth (D5) and ninth (D9) day after inoculation. Differentially expressed genes (DEGs) were identified by comparing root transcriptomes from inoculated plants with those from non-inoculated plants. Accessions A2-100 and A2-190 showed 52 and 29 DEGs on D5, respectively, with 14 DEGs in common, and 18 DEGs for A2-100 and 11 DEGs for A2-190 on chromosome 5. On D9, four DEGs were found in A2-100 and two DEGs in A2-190. For GB 713, 52 and 43 DEGs were found, and for TX 110, 29 and 117 DEGs were observed on D5 and D9, respectively. Six DEGs were common at the two sampling times for these genotypes. Some DEGs were identified as Meloidogyne-induced cotton (MIC) 3 and 4, resistance gene analogs, or receptor-like proteins. Other DEGs have potential roles in plant defense, such as peroxidases, programmed cell death, pathogenesis related proteins, and systemic acquired resistance. Further research on these DEGs will aid in understanding the mechanisms of resistance to explore new applications for the development of resistant cultivars.


Introduction
Reniform nematode (Rotylenchulus reniformis Linford and Oliveira) is an obligate parasite that has a wide host range, which hinders management of this pest in several important crop species [1].The nematode is a sedentary, semi-endoparasite, and has become one of the most important nematode species present in soils of the cotton-producing regions of the United States [1,2].In the southeast states, such as Mississippi, Alabama, and Louisiana, reniform nematode has replaced Meloidogyne incognita (root-knot nematode) as the major nematode species infecting cotton [1,3].Female vermiform nematodes can infect cotton roots throughout the growing season [4] and during this process, the nematode will modify host cells at the feeding site to form a syncytium, which provides the sole nutrient source for nematode development and reproduction [5].Feeding on the root system interferes with the uptake of water and nutrients by the plant, resulting in stunting, delayed maturity, and yield reductions [4,6].Moreover, root damage from nematode feeding makes cotton plants more vulnerable to soilborne diseases such as Fusarium wilt, which can compound yield losses [7].In 2022, it was estimated that yield losses caused by reniform nematodes were 36.67 thousand metric tons, which was 1.2% of the U.S. total cotton production, and represents a loss of USD 74 M [8].Yield losses in specific regions of the country were, however, significantly higher than the nationwide average with losses as high as 50% reported [9,10].
Limited approaches are available to manage reniform nematodes in Upland cotton (Gossypium hirsutum L.) production.Nematicide soil fumigation, seed treatment, in-furrow and foliar applications have been used to mitigate yield losses [1,10].However, only a few nematicides are currently available and their efficacy can vary based on growing conditions [11][12][13][14].Nematicides are effective in providing early season protection, but the nematode population can quickly and dramatically increase during the growing season due to the nematode's short life cycle and high reproduction rate [1].Additionally, nematicides are expensive and not all nematicide applications are profitable [14,15].Human health and environmental concerns also need to be considered with the use of nematicides.Nematode management using crop rotation with non-host crops, such as corn and peanuts, or with resistant host cultivars, such as those available for soybean, has been recommended to reduce nematode populations [11,16].This approach will require long rotational cycles to significantly lower the nematode population and therefore crop rotation is not always feasible due to the economic and resource constraints associated with cotton production [1,11].Growing resistant Upland cotton cultivars or combining the application of nematicides with the use of resistant cultivars has been reported as an effective management strategy [1,11,[17][18][19].McCulloch et al. (2021) reported that resistant cultivars significantly suppress the reniform nematode population resulting in a 26% seedcotton yield increase compared to susceptible controls [20], whereas Koebernick et al. (2021) reported 8-20% yield increase using a resistant cultivar with a nematicide application during planting [17].
The development of resistant Upland cotton cultivars for reniform nematode management will require the use of multiple sources of resistance to limit the ability of the nematode to overcome a single source of resistance.Upland cotton accounts for 97% of cotton production in the U.S. (https://www.ers.usda.gov)(accessed on 1 October 2023); however, extensive screening of the Upland cotton germplasm collection only identified a few sources with moderate resistance [21][22][23][24].These sources of resistance showed an inconsistent response [21,25] and have not been employed in commercial cultivar development.More desirable sources of resistance have been identified from tetraploid species G. barbadense L. [21,22] and several diploid Gossypium species, including G. aridum (Rose & Standl.)Skovst, G. arboreum L., G. herbaceum L., and G. longicalyx Hutch.& Lee [21,[26][27][28].
The diploid cotton species are also a potential source of novel resistance genes; however, transferring resistance to Upland cotton is technically demanding and laborious.G. longicalyx was the only cotton species reported to show an immune response to reniform nematode infection [21].This source of resistance was successfully transferred to Upland cotton using the bridging lines (G.hirsutum × G. longicalyx) and (G.hirsutum × G. herbaceum) with the breeding lines LONREN-1 (Reg.No. GP-977, PI 669509) and LONREN-2 (Reg.No. GP978, PI 669510) being developed [37][38][39].The G. longicalyx resistance gene Ren lon was mapped on chromosome AD_ch21_Dt.11[40].Resistance derived from G. longicalyx was associated with a hypersensitive response resulting in root necrosis and severe plant stunting, therefore this resistance source has not been used for commercial cultivar development [41,42].Bhandari et al. (2015) [43] reported hypersensitive reactions in G. hirsutum germplasm lines LONREN-1 and LONREN-2 with G. longicalyx resistance, and G. arboreum accession A2-190, to two isolates from Louisiana based on plant height reductions; however, reductions in root growth for A2-190 was not reported.Severe root necrosis associated with resistance has not been reported for other diploid cotton species showing high levels of resistance, although, limited research has been conducted on these resistant sources.Resistance from G. aridum and G. arboreum were transferred to Upland cotton using the (G.hirsutum × G. aridum) bridging line [27,44].The resistance QTL Ren ari derived from G. aridum was mapped to AD_ch21_Dt.11 in Upland cotton [27] and at least two resistance genes were reported for the G. arboreum source of resistance [44].Limited breeding research has been conducted on these resistance sources and no Upland cotton breeding lines have been released.
The G. arboreum germplasm collection includes more than 1600 accessions and resistant genotypes have been frequently observed in this collection [28].A few G. arboreum accessions have been genetically characterized for nematode resistance with resistance being conferred by one or a few genes [45][46][47].A genome-wide association study of 246 accessions revealed 15 SNPs significantly associated with reniform nematode resistance [48].The lack of genomic tools has hindered the transfer of resistance from G. arboreum to Upland cotton.To enhance the utilization of resistance sources, in the present study, the transcriptomes of G. arboreum and G. barbadense resistant genotypes responding to reniform nematode were evaluated to identify differentially expressed genes and genomic regions potentially associated with nematode resistance.

Sequencing Results
RNA sequences were obtained for 47 of the 48 treatments due to library construction failing for 1 sample of the non-inoculated GB 713 for the D9 treatment.The number of raw reads for each sample ranged from 19.64 to 38.04 M for single-end sequencing and 61.08 to 74.75 M for paired-end sequencing, with an average of 28.61 M raw reads.The number of bases generated from each sample ranged from 1.98 to 3.84 giga bases (Gb) for single-end sequencing and 6.17 to 7.55 Gb for paired-end sequencing, with an average of 2.89 Gb.Over 96% of bases had quality scores Q20 or higher (error p ≤ 0.01) and 93% of bases had quality scores Q30 or higher (error p ≤ 0.001).The average undetermined bases (N) were 0.03%.

Differential Expression of Genes in A2-100
For the D5 treatment, 29,168 genes were expressed with 57 genes differentially expressed between the inoculated and non-inoculated treatments.Expression was upregulated for 56 of these genes.The expression levels for 52 genes changed by at least two-fold.At least 1 differentially expressed gene (DEG) was observed on each of the 13 chromosomes, with 18 DEGs recorded for chromosome A05 (Table 1).Five DEGs located on chromosome A05 and two on chromosome A12 had uncharacterized functions.Disease resistance related genes were observed, including 11 genes encoding receptor-like proteins and 2 genes encoding Meloidogyne-induced cotton (MIC) 3 and 4 proteins (Table 2).Genes that play important roles in plant immune and defense reactions were also recorded, such as plant programmed cell death (one metacaspase-9-like, two ervatamin-B-like, and one self-pruning), disease resistance reactions (three peroxidase, two hevamine-A-like genes, one glucan endo-1,3-beta-glucosidase 7-like), and systemic acquired resistance (one protein NIM1-interacting 1).In addition, four transcript factors were differentially expressed.Data from the D9 treatment showed 4 out of the 29,072 genes were differentially expressed.The expression levels were 4 to 32 times higher in the inoculated treatment compared to the non-inoculated control and these DEGs were upregulated.Two of the DEGs were associated with chromosome A05 (Table 1) and they also showed differential expression for the D5 treatment.These four DEGs showed functions associated with plant defense (Table 2).

Differential Expression of Genes in A2-190
A total of 29,239 genes were expressed for the D5 treatment with 56 genes differentially expressed between the inoculated and non-inoculated treatments, and 44 of these DEGs were upregulated.The expression levels for 29 DEGs showed more than a two-fold change with 25 DEGs showing upregulation.The 29 DEGs were distributed on 9 chromosomes, including 11 genes located on chromosome A05 (Table 1).Two DEGs associated with resistance to nematodes, Meloidogyne-induced cotton genes (MIC-3 and MIC-4), were observed (Table 2).Several other DEGs have functions associated with plant defense and disease resistance, including two resistance gene analogs (DMR6-like), three receptor-like kinase genes, peroxidase-like, hevamine-A-like, jasmonate-induced protein, and transcription factor MYC2 (XM_017774272.2).Six DEGs on chromosome A05 and two on chromosome A10 were uncharacterized for gene function.
From the 29,096 genes expressed for the D9 treatment, 2 genes, residing on chromosomes A05 and A12, were differentially expressed (Table 1).Expression levels increased 5 and 11 times in the inoculated treatment compared to the non-inoculated control.These DEGs had functions associated with plant defense (Table 2).

Differential Gene Expression for TX 110
For the D5 treatment, expression was detected for 65,249 genes.A total of 31 of these genes were differentially expressed between the inoculated and non-inoculated treatments with 27 DEGs upregulated.The change in expression ranged from 2-to 120-fold for 29 DEGs and these genes were distributed on 18 chromosomes (Table 3).DEGs were identified on chromosomes of the A and D subgenomes with a single DEG per chromosome more frequently recorded.Five of these DEGs showed gene functions associated with plant defense (Table 4) and three genes were transcription factors.No gene function was recorded for the majority of the DEGs.
A total of 64,499 genes showed expression for the D9 treatment with 119 genes differentially expressed between inoculated and non-inoculated treatments, including 82 genes upregulated.At least a two-fold change in expression was recorded for 117 DEGs.The 117 genes were distributed on all 26 chromosomes (Table 3).Multiple DEGs were detected on each chromosome with nine genes observed on chromosomes AD_ch09_At.09 and AD_ch26_Dt.12.Twenty-nine DEGs had functions associated with plant defense (Table 4).Additionally, several transcription factors were identified.However, the gene function for the majority of the DEGs was uncharacterized.

Differential Gene Expression for GB 713
Gene expression analysis identified 65,316 genes for the D5 treatment using the genome of Pima 90 as reference.Overall, 66 genes were differentially expressed with 61 DEGs upregulated.Fifty-two of these DEGs showed at least a two-fold change in expression levels between inoculation treatments.The 52 DEGs were distributed across the 26 chromosomes with a relatively small number located on each chromosome (Table 3).The function for the majority of the DEGs was uncharacterized.Four DEGs showed a function associated with plant defense (Table 5) and transcription factors, such as Myb, were also observed.For the D9 treatment, 43 genes were differentially expressed from a total of 64,073 genes showing expression.Expression levels of these DEGs varied from 2-to 48-fold with 28 genes upregulated.DEGs were distributed on 19 chromosomes, including 6 genes on AD_ch05_At.05 and 5 on chromosome AD_ch11_At.11(Table 3).Eight DEGs showed functions association with plant defense, including one resistance gene analog (Table 5).Two transcription factors were also differentially expressed, but most DEGs were uncharacterized.
Table 5. Potential plant defense related DEGs and chromosomal locations for the G. barbadense genotype GB 713 evaluated on the fifth and ninth day after inoculation (DAI) with reniform nematodes and analyzed using three reference genomes.
Of the 60,529 genes observed for the D5 treatment using BARBREN713-32 genome as reference, 41 genes were differentially expressed with 38 genes upregulated.These genes were distributed on 19 chromosomes, including 5 on chromosome AD_ch07_At.07and 4 on chromosome AD_ch05_At.05(Table 3).The functions of 17 genes remained unknown.One transcription factor and two genes involved in abiotic responses were identified; however, no DEGs associated with plant defense were identified.For the D9 treatment, 25 genes showed differential expression from a total of 44,266 genes evaluated with expression for 16 DEGs upregulated.These 25 DEGs were associated with 15 chromosomes with three located on chromosome AD_ch05_At.05(Table 3).Three DEGs showed functions associated with plant defense (Table 5).

Comparisons of Differential Expression within and between Gossypium Species
A total of 14 genes were differentially expressed for the D5 treatment in A2-100 and A2-190 with 7 DEGs mapped to chromosome A05 (Table 6).Six DEGs showed functions associated with disease resistance, including Meloidogyne-induced cotton genes MIC-3 (XM_017750577) and MIC-4 (XM_017750341), peroxidase (XM_017747283), hevamine-A-like (XM_017781096), metacaspase-9-like (XM_053030089), and early-nodulin-75-like (XM_17753691) genes with three DEGs mapped on chromosomes A05.One gene (XM_017773121.2) showed differential expression for the D5 and D9 treatments in A2-100 and for the D9 treatment in A2-190 and was mapped on chromosome A05.This aspartyl protease gene has an important role in plant defense.Twelve genes were differentially expressed in TX 110 and GB 713 using Pima 90 as a reference genome (Table 7).Three genes (GbM_A02G0072, GbM_A01G0832, and GbM_D01G0797) showed differential expression across the genotypes on D5 with two genes (GbM_A04G0525 and GbM_A13G1710) showing differential expression on D9.On D5 and D9, three genes showed differential expression for TX 110, but only showed differential expression on one of the sampling treatments for GB 713.Two genes (GbM_A11G3523 and GbM_A10G2857) were differentially expressed across the genotypes on D5 and D9.Using the functional annotation of the DEGs as a comparison across species and genotypes, some similarity in disease resistance related genes were observed.For example, receptor-like proteins, peroxidase, hevamine-A-like, and glucan endo-1,3-beta-glucosidase were differentially expressed across the four genotypes.Metacaspase-9-like and aspartyl protease were differentially expressed for G. barbardense TX 110 and for G. arboreum A2-100 and A2-190, whereas the Meloidogyne-induced cotton gene MIC-3 was differentially expressed for G. barbardense GB 713 and for A2-100 and A2-190.

Discussion
The evolutionary and commercial breeding history for Upland cotton has resulted in a narrow germplasm pool, which makes current cultivars vulnerable to changes in pest populations.Although reniform nematode was first identified in cotton in 1940, release of resistant commercial cultivars has only occurred within the past five years.The development of resistant cultivars has been hindered by the lack of resistance genes within Upland cotton germplasm, the difficulty in screening for reniform nematode resistance, and the complex inheritance of resistance.Related cotton species are a valuable resource of resistance genes; however, genetic and genomic evaluations are essential to identify genes associated with resistance to develop marker-assisted selection strategies.
Gene expression studies are one approach to identify genes associated with resistance.Transcriptome analysis of G. hirsutum cultivars that are susceptible (DP90 and SG747), resistant (BARBREN-713) and hypersensitive (LONREN-1) to reniform nematodes revealed many DEGs associated with resistance, such as cell wall architecture, hormone metabolism and signaling, ROS levels, cell death pathways, and pathogenesis [49].In the present study, DEGs with a range of functions were identified from the four cotton genotypes in response to reniform nematode infection.Resistant gene analogs and non-specific resistance genes, such as chitinase, glucan endo-1,3-beta-glucosidase, and pathogenesis-related proteins, were differentially expressed and may have an important role in response to reniform nematode infection.The DEGs were more often upregulated than downregulated, higher expression of certain genes may help to fight the parasites and mitigate the damage of nematodes.Li et al. (2015) [49] reported that more DEGs were downregulated for GB 713, whereas the other two evaluated genotypes show slightly more DEGs being upregulated.
Meloidogyne induced cotton genes MIC-3 and MIC-4 were differentially expressed in A2-100 and A2-190, and MIC-3 gene was differentially expressed for GB 713, and they were upregulated in these genotypes.MIC-3 was found to be specifically expressed in cotton resistant genotypes when inoculated with M. incognita [50].Overexpression of MIC-3 in G. hirsutum cv.Coker 312 reduced egg production of M. incognita by 60-75%, but no effect on reniform nematode reproduction was reported [51].These genes showed expression in A2-100 and A2-190 on D5, which may indicate that expression of these genes occurs early in the infection process and may influence the establishment of the feeding site.In contrast, MIC-3 expression was observed in GB 713 on D9, suggesting a possible role in hindering nematode development.These DEGs were associated with chromosome A05 in both species.Several cotton resistance genes/QTLs to biotic and abiotic stresses have been mapped on chromosome AD-ch05_At.05[52][53][54], which may result from natural selections in cotton evolution.The physical locations on cotton genomes for these genes/QTLs were unknown, thus their distance to the resistance genes identified in this study are undetermined.
Genes for receptor-like proteins showed expression for A2-100 and A2-190 on D5 and in GB 713 on D9.These DEGs were upregulated for these genotypes.Receptor-like proteins are cell surface receptors composed of several distinct domains, including signal peptide, extracellular leucine-rich repeat (LRR) region, and transmembrane domain, with these proteins playing an important role in plant development as well as disease resistance [55].DEGs of receptor-like proteins were associated with chromosome A01 across the cotton genotypes in the present study, which may suggest a conservation response.Additionally, Li et al. (2018) [48] reported one SNP significantly associated with reniform nematode resistance for chromosome A01 in G. arboreum, which is about 9 Mb from the receptor-like gene XM_053030958.1 (Table 2).
Hevamine-A-like protein is a type of plant chitinase and lysozyme that are important for plant defense against pathogens [56].Chitin is an essential component of the nematode eggshell and pharynx, and the disturbance of chitin synthesis or hydrolysis could lead to nematode embryonic lethality, defective egg laying, or molting failure [57].Upregulated differential expression for hevamine-A-like protein was recorded for the D5 treatment across the four genotypes and for the D9 treatment in A2-190 and TX 110.This DEG was commonly associated with chromosome A12 across the four genotypes.A genome-wide association study of G. arboreum genotypes also identified candidate genes for reniform nematode resistance on chromosome A12 [48].
Several other genes associated with plant defense were also recorded for the genotypes evaluated in the current study.Multiple peroxidase genes were differentially expressed for A2-100, A2-190, and GB 713, which were upregulated.Peroxidase genes are critical for cell wall stiffening [58].Li et al. (2015) [49] reported peroxidase genes were differentially expressed in response to reniform nematode infection across the four genotypes that showed different infection responses.Aspartyl protease family protein At5g10770-like showed differential expression and was upregulated for A2-100, A2-190, and TX 110.This protein could degrade pathogen effectors contributing to reducing virulence and degrade pathogenesis-related proteins that would induce the expression of genes involved in stress and defense responses, innate immunity, and systemic acquired resistance (SAR) [59].A metacaspase-9-like gene was recorded for A2-100 and TX 110, and two ervatamin-B-like genes and one self-pruning gene were observed for A2-100.These genes have critical roles in programmed cell death during plant development and defense responses [60][61][62].Two genes differentially expressed for A2-100 and one gene for TX 110 encode early-nodulin like proteins.These proteins were reported to have an important role for enhancing fitness of the pathogen during host colonization [63].In addition, these proteins have functions associated with the transportation of nutrients and plant development [63].One DEG (XM_017774501.2) recorded for A2-100 encodes the NIM1-interacting protein.This gene has an important role in systemic acquired resistance [64].Two Kunitz trypsin inhibitor 5 genes were differentially expressed for TX 110.These genes were reported to have functions associated with the protection of plants from insect predators [65].One nematocidal crystal cry1Ag gene observed for TX 110 may have resulted from Bacillus thuringiensis contamination [66].
Mapping reniform nematode resistance genes will be critical for the development of molecular markers to enhance the introgression of these genes from related Gossypium species for the development of resistant Upland cotton cultivars.Fifteen SNPs significantly associated with reniform nematode resistance were mapped to eight chromosomes for G. arboreum [48].DEGs reported in the current study were widely distributed across the genome with DEGs observed on 13 chromosomes for A2-100 and 9 chromosomes for A2-190; however, the majority of the DEGs were associated with chromosome A05 for the 2 genotypes.Thus, DNA markers developed from chromosome A05 should be selected for screening segregating populations for nematode resistance derived from G. arboreum sources.DEGs associated with chromosome A12 may also have a role in nematode resistance for the G. arboreum genotypes.Some similarities were observed for the DEGs across the genotypes and these data may suggest some shared resistance genes; however, identifying G. arboreum genotypes with unique resistance genes will be necessary to increase genetic diversity.
DEGs were also observed on chromosome A05 of the A-genome and the homoeologous chromosome 19 of the D-genome for TX 110 and GB 713 with more DEGs frequently associated with these chromosomes for the D9 treatment.DEGs were distributed across most of the chromosomes for these two genotypes, but the frequency of DEGs associated with individual chromosomes varied across genotypes and across inoculation treatments.The occurrence of DEGs across the genome may suggest expression of genes related to root tissue damage and not associated with nematode resistance.The resistance QTLs transferred from GB 713 were mapped to AD_ch21_Dt.11 and AD_ch18_Dt.13 in the genome assembly of BARBREN-713-32 [34].Transcriptome analysis of GB 713 and TX 110 in the current research identified DEGs associated with these two chromosomes and the DEGs were in or adjacent to the intervals identified for BARBREN-713-32.Two DEGs for TX 110 and one for GB 713 that were mapped onto these chromosomes had functions associated with plant defense.A few DEGs were comparable across the two genotypes and may indicated some similarity in the mechanism of resistance.
Sampling time is very critical for the transcriptome analysis of cotton responding to reniform nematodes.While inspecting nematode infection on a small subset of plants two DAI, it was found that not all plants were infected, and the G. arboreum accessions had little root development.All sampled plants were infected on D5, and plant root mass was enough for RNA extraction.The parasitizing nematodes were in the gravid stage on D9.Therefore, root samples were collected on D5 and D9 in this study.A comparison across the two inoculation treatments showed different gene expression patterns among the genotypes evaluated in the present study.The number of DEGs recorded for the D9 treatment was dramatically reduced compared to the D5 treatment for the G. arboreum genotypes.These data could suggest resistance is associated with the early stages of nematode infection.However, fewer DEGs were observed for A2-190 for the two inoculation treatments compared to A2-100.This may suggest that resistance associated with A2-190 could hinder the establishment of the feeding site for the nematode, resulting in the low number of nematodes typically observed on the root system for this genotype.Collecting root samples prior to D5 may aid in the identification of additional DEGs associated with resistance.However, transcript evaluations prior to D5 may require a different sampling approach because of the few nematodes infecting the root system for highly resistant genotypes and low root weights characteristic of G. arboreum accessions.In contrast, resistance associated with A2-100 may occur later in the infection process, leading to a greater number of nematodes infecting the root system [67]; thus, contributing to the moderately resistant response recorded for this genotype and the higher number of DEGs observed.
The difference in the number of DEGs across inoculation treatments was reduced for the G. barbadense genotype GB 713; although, fewer DEGs were recorded for the D9 treatment.These data could suggest a different mode of resistance compared to the G. arboreum genotypes.This trend was reversed for TX 110 with more DEGs recorded for the D9 treatment than the D5 treatment.A moderately resistant response to reniform nematode infection has been reported for TX 110 [22].When infected with the MSRR04 reniform nematode population, TX 110 and GB 713 showed similar nematode development for two days after inoculation; then, more rapid progression of nematode development was observed for TX 110 compared to GB 713 [67].The expression of many non-specific pathogenesis-related genes was also observed for TX 110 for the D9 treatment in the current study, which could be associated with the moderately resistant response.The availability of data on nematode development for the two genotypes could be used in future studies to select additional time points to evaluate changes in gene expression to assess the pathways contributing to the resistant response.
Reference genomes are critical for transcriptome analysis.ShiXiYa 1 (PI 615743) was reported to have a moderately resistance response to reniform nematode infection [28].Using this reference genome, many receptor-like DEGs were identified from the two G. arboreum genotypes evaluated in the current study.While using the Pima 90 genome as a reference, relatively less receptor-like DEGs were identified from the two G. barbadense genotypes; this could result from the lack of resistance in Pima 90 to reniform nematode (its response to reniform nematode needs to be tested).Three reference genomes were used for gene identification and feature count when analyzing the transcriptomes for GB 713.The number of non-specific plant defense related DEGs recorded was reduced using the G. hirsutum genotypes BARBREN 713 and BARBREN 713-32 resistant to reniform nematodes.Since their resistance was derived from GB 713, it is possible that these DEGs of GB 713 confer resistance to reniform nematode.Three DEGs were found in each of the G. hirsutum reference genomes, suggesting some resistance genes are common in these lines.The discrepancy in the numbers of DEGs identified may result from the selections in breeding or the difference of genome assemblies and annotations.
DEGs associated with disease resistance genes identified from the current research will be further tested using other resistant and susceptible cotton genotypes to determine their role in suppressing reniform nematode infection and development.More than 100 resistant genotypes were identified from the G. arboreum germplasm collection [28] and would be useful for these evaluations.Accessions from this collection showed a wide range of variations for infection response, which could be useful for identifying additional DEGs associated with nematode resistance.Some accessions showed very low numbers of nematodes infecting the root system and resistance could be associated with the establishment of the feeding site, whereas hindering nematode egg production could be associated with host resistance for other accessions.Increasing root sampling times would also be useful to identify additional DEGs associated with this variation in infection response although constitutively expressed resistance genes cannot be identified using transcriptome analysis.Molecular markers for resistance genes can be developed using segregating populations, which would be beneficial for the introgression of resistance in Upland cotton breeding programs.

Plant Materials, Reniform Nematode Inoculation, and RNA Sequencing
Two G. arboreum accessions A2-100 (PI 529728) and A2-190 (PI 615699), and two G. barbadense accessions TX 110 (PI 163608) and GB 713 (PI 608139) were used for this study.These accessions were originally obtained from USDA-ARS Germplasm Resource Information Network (GRIN) and have been self-pollinated over multiple generations at Stoneville, MS or at the Cotton Winter Nursery, Tecoman, Mexico to develop homogenous breeding lines.The reniform nematode response for these genotypes has been evaluated across multiple experiments.The A2-100 genotype showed a moderately resistant reaction [45], TX 110 showed moderate to high resistance when tested with different reniform nematode populations [21,22,68], whereas genotypes A2-190 and GB 713 showed a highly resistant response [22,44].A susceptible G. hirsutum cultivar Deltapine 16 (PI 529251) was also included solely for checking inoculation.
These genotypes were inoculated with the Mississippi reniform nematode population MSRR04 [69] following the procedure described by Stetina et al. (2014) [70].Briefly, a single seed was planted in a conical plastic pot (Ray Leach SL-10 Cone-tainer, Stuewe & Sons, Inc., Tangent, OR, USA) containing 120 cm 3 of a steam-sterilized soil mixture of one-part sandy loam soil and two-parts sand.For this study, 60 seeds were planted for each genotype and the experiment was conducted in 2 growth chambers.One chamber was used for the inoculation treatment, whereas the second chamber was used for the non-inoculated control treatment.Growth conditions for the experiment were set to a constant temperature of 28 • C with a 16 h photoperiod.Genotypes were arranged in a completely randomized design with 3 replications and 10 pots per replication.Seven days after planting, the nematode inoculation treatments were conducted by infesting the soil in each pot with approximately 1000 nematodes (mixed vermiform life stages).On the fifth (D5) and ninth (D9) day after inoculation (DAI), three plants of each replication for each genotype were randomly selected to harvest root samples.The corresponding control samples were also harvested at each time interval.Briefly, the seedlings were removed from the pots and the soil was washed from the roots.The shoot was removed from the root system below the soil line and the roots were further washed with sterilized reverse osmosis water.The root systems from the three selected plants were combined and ground to a fine powder in liquid nitrogen and then transferred to 25 mL tube for storage at −70 • C. RNA extraction was conducted using the Spectrum Plant Total RNA kit (Sigma-Aldrich Inc., St. Louis, MO, USA) and quantified using a spectrophotometer.Libraries were constructed from total RNA using the Illumina TruSeq RNA Sample Preparation Kit v2 (Illumina, Inc., San Diego, CA, USA) with set B indexed adapter sequences.The transcriptomes were sequenced with Illumina Hiseq 1000 platform at MacroGene, Seoul, Korea.Two flowcells were used for sequencing with 42 single-end sequenced samples and 6 paired-end sequenced samples.

Data Analysis
Data analysis was conducted on the SCINet high performance computing platform (scinet.usda.gov).The raw RNA reads were checked with BBduk (sourceforge.net/projects/bbmap/, version 39.03) to remove adapters and low-quality bases (p = 0.1), then reads were aligned to reference genomes using the software Spliced Transcripts Alignment to a Reference (STAR, version 2.7.11b) [71] to identify expressed genes and count the features following the software manual.The A2 (G.arboreum) reference genome of ShiXiYa 1 (PI 615743) was used to analyze the transcriptomes for A2-100 and A2-190 [72].The reference genome for Pima90 (G.barbadense) was used to analyze the transcriptomes for TX 110 and GB 713 [73].In addition, the reference genomes of G. hirsutum germplasm lines BARBREN-713 (B713) and BARBREN-713-32 (Bar32), whose reniform nematode resistance was derived from GB 713, were used to analyze the transcriptomes of GB 713 [34].The genomic sequences and transcriptomes of these reference genomes were downloaded from CottonGen (https://www.cottongen.org)(accessed on 15 June 2023) [74].The transcripts aligned to the reference genomes and their feature counts derived from STAR analysis were further used to identify DEGs with the R package DESeq2 (adjusted p < 0.05) [75].DEGs with at least a two-fold change of expression levels were further analyzed.The sequences and chromosomal locations of the DEGs were obtained from the relevant transcript sequences of the reference genomes and the potential functions of the DEGs were determine by searching against the nr database (https://www.ncbi.nlm.nih.gov)(accessed on 31 August 2023) with BLASTX (p < 10 −6 ).

Conclusions
In this study, the transcriptome analysis of four resistant Gossypium genotypes challenged with reniform nematodes revealed differentially expressed genes and their chro-mosomal locations.Plant defense related genes, including nematode resistant genes, were identified.Some plant defense genes were common within and between Gossypium species.Further study of the identified resistance genes may help to understand the mechanism of resistance and to develop resistant cotton cultivars.

Table 1 .
Number of DEGs and their chromosome locations for G. arboreum genotypes A2-100 and A2-190 evaluated on the fifth (D5) and ninth (D9) day after inoculation with reniform nematodes.

Table 2 .
Potential plant defense related DEGs and associated chromosomal location for G. arboreum genotypes A2-100 and A2-190 evaluated on the fifth and ninth day after inoculation (DAI) with reniform nematodes.

Table 4 .
Potential plant defense related DEGs and associated chromosomal location for the G. barbadense genotype TX 110 assessed on the fifth and ninth day after inoculation (DAI) with reniform nematodes.

Table 6 .
Comparison of DEGs and chromosomal locations for Gossypium arboreum genotypes A2-100 and A2-190 on the fifth (D5) and ninth (D9) day after inoculation with reniform nematodes using ShiXiYa 1 as the reference genome.
X: gene differentially expressed.

Table 7 .
Comparison of DEGs and chromosomal locations for Gossypium barbadense genotypes TX 110 and GB 713 on the fifth (D5) and ninth (D9) day after inoculation with reniform nematodes using Pima 90 as the reference genome.
X: gene differentially expressed.