Large-scale gene gains and losses molded the NLR defense arsenal during the Cucurbita evolution

Main conclusion Genome-wide annotation reveals that the gene birth–death process of the Cucurbita R family is associated with a species-specific diversification of TNL and CNL protein classes. Abstract The Cucurbitaceae family includes nearly 1000 plant species known universally as cucurbits. Cucurbita genus includes many economically important worldwide crops vulnerable to more than 200 pathogens. Therefore, the identification of pathogen-recognition genes is of utmost importance for this genus. The major class of plant-resistance (R) genes encodes nucleotide-binding site and leucine-rich repeat (NLR) proteins, and is divided into three sub-classes namely, TIR-NB-LRR (TNL), CC-NB-LRR (CNL) and RPW8-NB-LRR (RNL). Although the characterization of the NLR gene family has been carried out in important Cucurbita species, this information is still linked to the availability of sequenced genomes. In this study, we analyzed 40 de novo transcriptomes and 5 genome assemblies, which were explored to investigate the Cucurbita expressed-NLR (eNLR) and NLR repertoires using an ad hoc gene annotation approach. Over 1850 NLR-encoding genes were identified, finely characterized and compared to 96 well-characterized plant R-genes. The maximum likelihood analyses revealed an unusual diversification of CNL/TNL genes and a strong RNL conservation. Indeed, several gene gain and loss events have shaped the Cucurbita NLR family. Finally, to provide a first validation step Cucurbita, eNLRs were explored by real-time PCR analysis. The NLR repertories of the 12 Cucurbita species presented in this paper will be useful to discover novel R-genes. Supplementary Information The online version contains supplementary material available at 10.1007/s00425-021-03717-x.


Introduction
Cucurbita is an important genus of the Cucurbitaceae family that encloses independently domesticated species native to Central and South America (Castellanos-Morales et al. 2018). The taxonomy of this genus is still much debated, as the number of accepted species varies from 13 to 30 (https:// plants. usda. gov/ and https:// www. itis. gov/). A burst of morphological transformations along the Cucurbitaceae Communicated by Dorothea Bartels. evolutionary history was found after early genome duplication events (Guo et al. 2020). Cucurbita species can be divided into two main groups: (i) mesophytes, annual or short-lived perennial herbaceous vines; and (ii) xerophytes, perennials growing in arid zones (Kates et al. 2017;Castellanos-Morales et al. 2018). Cultivated Cucurbita species belong to the first group, of which five domesticated crops (C. argyrosperma, C. ficifolia, C. maxima, C. moschata, and C. pepo) are grown worldwide for their edible fruits and seeds (Castellanos-Morales et al. 2018), variously known as squash, pumpkin or gourd. The earliest known evidence of the domestication of Cucurbita is dated 10,000 calendar years B.P. (Renner and Schaefer 2016;Kates et al. 2017;Smith 1997), preceding the domestication of other crops (e.g., for maize is started 8,700 calendar years B.P.). At least six independent domestication events from distinct wild ancestors were identified in domesticated Cucurbita (Sanjur et al. 2002). Certainly, domestication and further breeding contributed to narrowing the genetic base of the species, making the Cucurbita crops susceptible to existing pests and diseases (Román et al. 2020). This vulnerability, along with the increasing interest in developing new, environmentally friendly Cucurbita-resistant cultivars, steers the identification and the study of structural and functional mechanisms of the evolution of pathogen-recognition genes. The wild cucurbits represent a valuable genetic resource for crop improvement (Khoury et al. 2020), in which novel disease-resistance genes could be identified and transferred into modern Cucurbita cultivars (Capuozzo et al. 2017). To date, the information on the defense arsenal of wild Cucurbita species is very limited, therefore, its exploration is an important goal for developing cultivars.
Basically, the plants have developed two pathogen-recognition layers to repel their attacks (PTI: PAMP-triggered immunity and ETI: Effector-triggered immunity) (Andolfo and Ercolano 2015). The pathogen perception mainly occurs through the recognition of effectors via plant disease-resistance (R) proteins. The major class of R-genes is the nucleotide oligomerization domain (NOD)-like receptors (NLRs), which encode nucleotide-binding site (NB) and leucine-rich repeat (LRR) gene families (Michelmore et al. 2013). NLRs are generally divided into three sub-classes, based on different N-terminus architecture, namely TIR-NB-LRR (TNL), CC-NB-LRR (CNL) and RPW8-NB-LRR (RNL) (Andolfo et al. , 2020. Recent studies demonstrated that NLRs can specifically detect pathogens using integrated domains (IDs) by monitoring the manipulation of host targets (Ortiz et al. 2017).
The number of NLRs varies dramatically between different plant genomes. Some crops such as hot-pepper, barrelclover, apple and wheat showed approximately 1000 NLRs (Andolfo et al. 2020). Diversely, the Cucurbitaceae family shows a lower number of NLR genes per lineage than other species (Sun et al. 2017;Andolfo et al. 2020). All genome-wide analyses showed that a large fraction of NLRs is arranged in clusters (Di . NLR gene cluster organization is very ancient, the first gene cluster involved 5 TNL genes and was identified in a green alga . The gene duplications have played an important role in the evolutionary history of the R-gene family, contributing to the diversification of the defense arsenal. Until now, the NLR repertory was described in few Cucurbitaceae sequenced genomes such as cucumber, melon, watermelon, zucchini and pumpkin (Garcia-Mas et al. 2012;Lin et al. 2013;Andolfo et al. 2017;Sun et al. 2017).
Recently, a set of Cucurbita genomic tools has become available for breeders and researchers, including de novo genome assemblies, transcriptomes, high-quality saturated maps, TILLING platforms and genome resequencing (Sun et al. 2017;Montero-Pau et al. 2018;Barrera-Redondo et al. 2019;Xanthopoulou et al. 2019). RNA sequencing is becoming increasingly a frequent technology orienting and driving the basic genetic research (Andolfo et al. 2014a;D'Esposito et al. 2019). Transcriptome sequencing encompasses a wide variety of applications from the prediction of protein-encoding gene function to the monitoring of dynamic gene expression regulation (Wang et al. 2009).
To date, several high-quality Cucurbita transcriptomes have been explored, most of them attributable to C. pepo cultivars (Blanca et al. 2011;Andolfo et al. 2017;Montero-Pau et al. 2018). However, still little is known about the Cucurbita species genetic diversity and even less has been done to explore their resistance gene arsenals. The availability of transcriptomic resources in Cucurbita provides an opportunity to conduct the first study on the evolutionary dynamics of the NLR-encoding gene family for a representative group of Cucurbita species. To this end, we used 40 transcriptomes and 5 genomes belonging to 12 Cucurbita species carrying out the NLR gene family annotation, characterization and diversification.

Identification of NLR-domain encoding genes
A total of 1,990,997 transcripts and 151,141 protein sequences were analyzed to identify and annotate the Cucurbita expressed-NLR (eNLR hereafter) genes in 40 de novo assembled transcriptomes and 5 Cucurbita genome assemblies, respectively. To identify the NLRdomain encoding genes, we translated the transcripts in all six possible reading frames (three in the forward direction and three in the reverse) using EMBOSS Sixpack (Madeira et al. 2019). The Multiple EM for Motif Elicitation (MEME) (Bailey et al. 2006) algorithm was used to decompose in motifs (Suppl. Table S4) the NB Pfam domain sequences (PF00931) of functionally characterized plant R proteins. The identified MEME motifs were used by Motif Alignment and Search Tool (MAST) (Bailey et al. 2006) to identify the NLR-domain encoding genes in our dataset. The Cucurbita proteomes were also scanned with the hidden Markov model (HMM) of the nucleotidebinding (NB Pfam ID: PF00931) using HMMER v3.0 (Finn et al. 2011). The analysis was carried out using the default cutoff value for statistical confidence. To detect additional Cucurbita eNLR/NLR genes that have escaped by the foregoing annotation, a reiterative process of gene identification was implemented. The annotated Cucurbita NB-encoding genes were used as queries against the 45 Cucurbita proteomes using a PSI-BLAST search (E-value cutoff of 1e-10). The protein domain architecture of identified Cucurbita eNLR/NLR proteins was annotated using IntetProScan v5 (Jones et al. 2014) and Conserved Domain Search (CD-Search) (Marchler-Bauer and Bryant 2004) with default parameters. Finally, a manual curation process was included in our annotation pipeline to discover the genes encoding single or partial NLR domains (Fig. 1). Cucurbita NLR sequences in FASTA format can be accessed at FIGSHARE with the following link (https:// figsh are. com/s/ b2808 ef8f5 a29a8 a522d). Users can download and use the data freely for research purposes only with acknowledgment to us and quoting this paper as a reference to the data.

Multiple sequence alignments and phylogenetic analysis
MAFFT (Multiple Alignment using Fast Fourier Transform) v6.814b was used to align the NB Pfam domain (PF00931) of both eNLR-/NLR-encoding genes using the L-INS-i algorithm (Katoh et al. 2002). Sequences with less than 50% of the full-length NB Pfam domain were excluded. Evolutionary analyses were conducted using MEGA 7 (Kumar et al. 2016). The phylogenetic relationships of Cucurbita NLR proteins were inferred using the maximum likelihood method based on the JTT model (Jones et al. 1992). The model with the lowest BIC score (Bayesian Information Criterion) was considered to describe the substitution pattern the best. The bootstrap consensus tree inferred from 100 replicates was taken to represent the evolutionary history of the analyzed sequences. The trees were drawn to scale, with branch lengths measured in the number of substitutions per site. Finally, the human Apaf1 gene was used to reroot the phylogenetic trees, because the split of plant and animal kingdoms predates the origin of NLR gene family groups (Andolfo et al. 2020).

Estimation of selection pressure and evolutionary divergence
Selective pressure acting on the Cucurbita NLRs was investigated by determining the nonsynonymous to synonymous nucleotide substitution (dN/dS) indicated as ω. Tests were conducted to estimate the evolution of each codon: positive (ω > 1); neutral (ω = 0); and negative (ω < 1). All positions with less than 80% site coverage were eliminated. All the NLR-coding DNA sequences were aligned using ClustalW v1.74 (Larkin et al. 2007). To depict the proportion of sites under selection, the evolutionary fingerprint analysis was carried out using the SLAC algorithm implemented in the Datamonkey server (Delport et al. 2010). The evolutionary divergence (π) was estimated as the number of base substitutions per site from averaging over all sequence pairs. Analyses were conducted using the Maximum Composite Likelihood model (Tamura et al. 2004). All positions containing gaps and missing data were eliminated. Evolutionary analyses were conducted in MEGA7 (Kumar et al. 2016).

Changes in NLR gene family size
The changes in Cucurbita NLR gene family sizes, across the dated phylogeny (Barrera-Redondo et al. 2019), were assessed using CAFE v4.2.1 with default settings (De Bie et al. 2006). To analyze significant expansion/contraction events of NLR orthogroups throughout the Cucurbita tree, we used the λ value (0.01397) estimated by Barrera-Redondo et al. (2019), which takes into account the genome assembly and annotation errors for all analyzed genomes (Han et al. 2013). The average expansion size indicates the mean number of genes gained or lost per orthogroup. C. sororia NLR gene family was excluded from evolutionary change analyses because our ortholog dataset was inadequate to correctly estimate a ʎ value including C. sororia.

Validation of eNLRs' expression by qPCR
NLRs' expression was validated by quantitative PCR (qPCR) in nine different transcriptomes of the analyzed Cucurbita accessions. Three seeds of each accession were disinfected with 10% sodium hypochlorite, rinsed with distilled water and then germinated by incubation over cotton in Petri plates at 37 °C for 48 h. Seedlings were grown in a climatic chamber at 26 °C and 60% of relative humidity, under a photoperiod of 16 h day/8 h night. Leaf tissue of plants at the third true-leaves stage was sampled and immediately frozen in liquid nitrogen. Total RNA was extracted using Extrazol (EM30) (Blirt, Gdańsk, Poland) and quantified using a NanoDrop 1000 spectrophotometer (Thermo Scientific™, San José, CA, USA). Aliquots of 1.5 µg were treated with Perfecta DNase I (Quanta Biosciences, Gaithersburg, MD, USA), and the first strands of cDNA were reverse transcribed with RevertAid synthesis Kit (Thermo Scientific™) and Oligo(dT) primer.
Four NLRs were amplified per accession using 39 primer pairs (Suppl . Table S5), which were designed by Primer 3Plus software (Untergasser et al. 2007). qPCR reactions were performed on a LightCycler ® 480 System (Roche Diagnostics, Mannheim, Germany), conducting two technical replicates and using 7.5 μL of the ready-to-use FastStart Essential DNA Green Master mix (Roche Diagnostics), 1.5 μL (1 µM) of each primer and H 2 O to a final volume of 15 μL. Reactions were conducted under the following conditions: 95 °C for 5 min, 40 cycles of 95 °C for 15 s, 60 °C for 30 s, and 72 °C for 15 s, ending up with a melting curve step (60 °C-95 °C) to check the specificity of each fragment. Relative expression of the NLRs amplified was calculated following the 2 −ΔCt method, described by Schmittgen and Livak (2008). The ubiquitin fusion protein (UFP) gene was used as an endogenous calibrator (Obrero et al. 2011).

NLR-resistance gene annotation in Cucurbita species
A total of 1603 eNLR and 256 NLR genes were identified in 40 transcriptomes and 5 genomes, respectively, using a hybrid strategy of gene annotation (Table 1 and 2; Suppl.  Tables S6 and S7). A fine R-gene motif-/domain-based search (Suppl . Table S4) was combined with a data-mining approach through a reiterative process (Fig. 1). In the first step, approximately 500 NB-encoding genes were identified, using MAST and HMMER software for sequence similarity searches. The resulting NB-encoding genes were annotated and classified based on the features of their N-terminal and C-terminal protein structures. Finally, all NLR-domain encoding sequences were used as queries against the entire 45 Cucurbita assemblies to check for additional genes that had not been previously identified (Fig. 1).
The full NLR complements ranged widely, from 87 to 25 members in C. moschata and C. argyrosperma, respectively (Table 2), a number substantially lower than in other plant species Barchi et al. 2019;Andolfo et al. 2020). Interestingly, C. sororia exhibited 51 NLR genes a number significantly (P < 0.05) higher than its domesticated counterpart, C. argyrosperma. Likewise, the analyzed transcriptomes, deriving from 40 Cucurbita accessions, come from 6 geographic regions of the world (Fig. 2), exhibited an average of 40 eNLR/NLR genes. Currently, C. pedatifolia (from Mexico) holds the highest number of eNLRs (73) among all Cucurbita spp., while, the lower number of eNLRs (10) was recorded in the 2 Spanish C. pepo accessions (Table 1). Differences in the total number of eNLRs within species were found (e.g., C. foetidissima; C. argyrosperma; C. okeechobensis) (Table 1). Moreover, only 4% (58 out of 1603) of annotated eNLRs encode the minimal domain structure (NB and LRR) necessary for a full-length NLR-encoding gene (Table 1). This finding was unlikely to reflect the true structure of 40 Cucurbitaresistance (R) gene families and might be due to the fragmented nature of the 40 transcriptomes since more than 70% (1162) of the eNLR proteins are encoded from very short transcripts (length < 700 bp). At least one full-length eNLR gene for each of the main R protein classes (CNL, RNL and TNL) was annotated in each of the analyzed Cucurbita sp. (Table 1). Generally, the NLRs belonging to the CNL class outnumbered those of the TNL class of one-third (Table 2). Interestingly, the full-length RNLs were especially numerous in the Cucurbita genomes (Table 2 and Suppl. Fig. S1).
The chromosomes 5, 6 and 14 of C. maxima and C. moschata genomes are particularly enriched of NLRs (over 50% of total genes), while NLRs are absent on chromosomes 8, 12 and 19. The 37% of C. sororia NLRs were located on chromosomes 5 and 9. Based on Richly et al. (2002) definition of a gene cluster (a spatial arrangement of NLRs by a gap of no more than 8 non-NLRs), we found that 51 (out of 87) and 12 (out of 53) NLRs resided in clusters in C. maxima and C. moschata, respectively. The largest cluster was located on chromosome 5 and includes 8 and 12 NLRs in C. maxima and C. moschata, respectively. C. sororia exhibited the biggest cluster on chromosome 9 including 8 TNLs. The 5 Cucurbita genomes hold 21 pairs and 6 triplets. Therefore, 77%, 55%, 26%, 22% and 65% of NLRs resided either in a gene cluster or in an array of two or three genes in C. moschata, C. maxima, C. argyrosperma, C. pepo and C. sororia, respectively (Table 3).
A substantial domain reshuffling was observed in the Cucurbita R-gene family, several NLRs showed atypical protein domain assemblies (Suppl. Fig. S2) Fig. S2).

Phylogenetic analysis of Cucurbita NLR genes
The evolutionary events, driving the Cucurbita NLR gene family, were inferred by analyzing the full NB domains extracted from each NLR gene annotated in four Cucurbita genome assemblies (Table 2). For comparative purposes, we included 96 well-characterized plant R-genes (Osuna-Cruz et al. 2018) and 1 out-group gene with an NB domain, the human Apaf-1 (Suppl . Table S3). A total of 106 NLRs were aligned together with the same reference R dataset ( Fig. 3a; Suppl. Fig. S3; Suppl. Table S3). The analyzed NB Pfam domains collapsed in 12 robust clades (bootstrap index > 80), of which only 5 (CNL4, CNL5, CNL6, RNL1 and RNL2) included at least 1 gene belonging to each analyzed genome ( Fig. 3a; Suppl. Fig. S3). Moreover, we found the direct orthologs only for Fom-2, Vat, ADR1 and NRG1 reference genes (Suppl. Fig. S3). The C. pepo and C. argyrosperma proteomes no-included Fom-2 homologs, while in C. maxima (8) and C moschata (13) were especially numerous (CNL3 clade in Suppl. Fig. S3). The BLASTn search of melon Fom-2 coding DNA sequence against the Cucurbita genome assemblies uncovered 64 and 76 significant exonmatches (E-value < 1e-6) in C. pepo and C. argyrosperma, respectively (Suppl . Table S8). Surprisingly, the clades CNL3 and TNL1 grouped ~ 35% of C. maxima and C moschata NLR complements. Overall, the difference in NLR arsenals is due to the gain/loss of specific gene subfamilies across Cucurbita species. To integrate and validate the above-mentioned findings, we performed a phylogenetic analysis on eNLR gene identified in 40 Cucurbita transcriptomes. A total of 243 amino acid sequences were aligned and collapsed into three groups (TNL in blue, RNL in green and CNL in red) supported by bootstrap values ≥ 50 (Fig. 4). Within the three groups, ten clades (CNL1 to 7, RNL1 and 2 and TNL1) were defined. Twenty-eight Cucurbita NLRs belonging to all domesticated species, except C. maxima, as well as to all wild types were harbored into the clade CNL1. They were identified in C. pepo accessions of the three subspecies: subsp. pepo pumpkins from Central America, subsp. ovifera acorn from USA and wild subsp. fraterna from Mexico, but not in the Spanish subsp. pepo accessions. Clade CNL2 comprises a CNL gene (CAR-CO098UC19951) annotated in the Indian accession C. argyrosperma PI-451712, the one with the highest NLR number. This gene showed very high homology to the melon gene Fom-2, which confers resistance to Fusarium oxysporum f.sp. melonis races 0 and 1 (Joobeur et al. 2004). Seventeen similar sequences were identified in 6 different Cucurbita spp. (Central America domesticate C. argyrosperma, C. ficifolia, C. pedatifolia, C. lundelliana and C. foetidissima and a Nigerian C. moschata accession) and are grouped in clade CNL3. It is interesting to note that clade CNL4 and CNL5 share a common ancestor (supported by 83% bootstrap indexes). CNL4 includes eight members from three domesticates, including C. maxima but not C. pepo, and two wild species, whereas Clade CNL5 holds genes from the economically important C. pepo species. This clade also includes genes from C. argyrosperma, C. ficifolia and C. moschata and from the wild types C. pedatifolia, C. lundelliana and C. foetidissima, all from Central America. The clade CNL6 contains nine CNLs. The direct ortholog (CFIC0097UC00660 identified in C. ficifolia) of Vat (virus aphid transmission) gene (Dogimont et al. 2014), conferring resistance to both A. gossypii and the viruses, falls into clade CNL7 (Fig. 4). Distinct from the canonical CNL genes are the RNL homologs divided into two robust clades (RNL1 and RNL2) (Fig. 4). RNL sequences are located on an ancestral position between TNL and CNL groups, and harbor homologs to the Arabidopsis ADR1 (Activated Disease Resistance 1) gene and Nicotiana benthamiana NRG1 (N requirement gene 1) (Peart et al. 2005;Bonardi et al. 2017). Overall, 27 and 21 Cucurbita RNL genes were grouped into the RNL1 (ADR1-homologs) and RNL2 (NRG1-homologs) clades, respectively. At least one ADR1-homolog was identified in each of the analyzed Cucurbita sp., while not all showed their NRG1-homolog.

Cucurbita NLR gene expansion and contraction events
The wide range of NLR family sizes and diversification prompted us to investigate the gene birth and death processes that occurred in the Cucurbita evolutionary history (Table 2). To assess the expansions and contractions of NLR family along the Cucurbita phylogenetic tree, we have used the orthology relationships indicated in the NLR phylogenetic tree ( Fig. 3a; Suppl. Fig. S3). We employed the global gene birth-death rate (λ) reported by Barrera-Redondo et al. (2019), to estimate the NLR gene gain-loss across Cucurbita phylogeny (Fig. 3b). Five NLR orthogroups (CNL1, CNL2, CNL3, TNL1 and TNL4 clades of Fig. 4a) significantly changed (P < 0.01) into the Cucurbita genus. Moreover, the TNL2 clade showed a lower, but significant (P < 0.05) level of change.
The average expansion size of C. moschata NLR family was 2, while C. argyrosperma (average expansion size = − 1) and C. pepo (average expansion size = − 0.92) showed a global contraction of NLR families. The difference in the proportion of C. maxima NLRs, across all orthogroups, was almost nil (average expansion size = 0.08). The branches leading to C. maxima and C. moschata exhibited the highest number of orthogroups The gene birth/death rate (λ) in the most recent common ancestor of the phylogeny and the whole-genome duplication (yellow star) in Cucurbita sp. were indicated in the rapid expansion (Fig. 3b). In particular, the Fom-2 homologs (CNL3 clade in Fig. 3a) were significantly expanded in C. maxima and C. moschata. The differences in Fom-2 homolog content across Cucurbita species were corroborated by several significant exon-matches identified into the Cucurbita assemblies (Suppl . Table S8). Even though few orthogroups (CNL2 and TNL4 clades in Fig. 4a) showed significantly rapid (P < 0.01) levels of change in the branch leading to C. pepo, this branch had the highest number of orthogroup (OG) changes (8 out of 12) within the entire phylogeny (Fig. 3b).

Evaluation of selection pressure acting on Cucurbita NLRs
To infer the direction and degree of natural selection acting on the Cucurbita NLRs, we estimated the variation level between the nonsynonymous substitution (dN) and synonymous substitution (dS) values. Results of the neutrality test performed for eight (out of ten) phylogenetic clades (Fig. 3) using the SLAC method (Pond and Frost 2005) are summarized in Table 4. The selection pressure expressed as dN/dS (ω) values of eight NLR-coding Clades are collapsed based on a bootstrap value over 50 and numerated. The TNL, RNL and CNL groups are drawn with blue, green and red backgrounds, respectively sequence alignments greatly varies from 0.48 to 1.2. Purifying selection (ω < 1) was observed in TNL1, RNL1-2, and in all CNL clades excepted for CNL5 (ω = 1.2). A consistent number (until 49) of codons under negative selection were identified for each clade (Table 4). A global stabilizing selection (ranging from 0.49 to 0.57) is acting against polymorphic variants in all the Cucurbita R-gene groups (Fig. 3), indicating the need to preserve their state (Table 4). The single codon analysis underlined 6, 6 and 1 positively selected sites in CNL1, RNL1 and RNL2 clades, respectively (Table 4). Five positively selected codon sites of ADR1-homologs (RNL1 clade in Fig. 3) are located on the NB Pfam domain (PF00931). Probably, the global protein structure of ADR1-like genes has been conserved, but the positive selection in specific sites of NB domains has been promoted to generate novel variability. The overall analysis of CNL and RNL phylogenetic groups underlined a conspicuous number of codons under positive (21 and 13 in CNL and RNL, respectively) and negative (98 and 67 in CNL and RNL, respectively) selection (Table 4).
Finally, the evolutionary divergence (π values in Table 4), calculated as mean base substitution per site over all sequence pairs, oscillated from 0.035 (TNL1 clade) to 0.159 (RNL2 clade). The low level of nucleotide diversity within eight clades indicated a limited genetic variation between phylogenetic homologs. While comparing the three phylogenetic groups (Fig. 3), the Cucurbita RNLs exhibited the highest nucleotide variability (π = 0.094). Although, the ADR1-/NRG1-homologs are conserved in Cucurbita spp., with respect to well-characterized plant R-genes, they exhibited the widest range of sequence diversity in coding regions.

Molecular validation of Cucurbita expressed-NLR genes
To verify that the de novo assembled eNLRs were present in Cucurbita transcriptomes a molecular analysis was carried out. A total of 36 eNLR genes selected in 9 Cucurbita accessions were molecularly validated by qPCR. Specific fragments were amplified to all the genes tested, verifying their transcription in the assayed accessions (Supp. Fig. S4).
Relative expression values (2 −ΔCt ) were consistent with the level of transcripts detected by sequencing, which demonstrates the veracity of the eNLR gene expression in the studied transcriptomes.

Discussion
A reiterative process implemented in an ad hoc R-gene annotation approach was used to identify the genes that encode domains similar to plant R proteins in 45 Cucurbita de novo assemblies (including both transcriptomes and genomes). The NLR gene families show a wide range variations into analyzed Cucurbita species. In our study, NLR gene families ranged from 87 to 27 members in C. moschata and C. argyrosperma, respectively. On average, the number of NLR loci in the Cucurbitaceae family is significantly lower than in other plant families (Lin et al. 2013;Andolfo et al. 2020). Recently, Guo et al. (2020) identified wholegenome duplication (WGD) events in the Cucurbiteae, which have surely contributed to their evolutionary diversification. The polyploidizations events provide abundant duplicated genes (Soltis and Soltis 2016), but only a small fraction is maintained (Zhang et al. 2012). Probably, most of the ancient duplicated NLRs were lost before the cucurbit speciation.
On average, a slightly higher number of CNLs compared to TNLs was identified in most Cucurbita assemblies. A similar CNL\TNL ratio was found in Solanaceae Andolfo et al. 2021), while the Brassicaceae and Fabaceae contain from twofold to sixfold more TNL than CNL genes (Guo et al. 2011;Kang et al. 2012). Surprisingly, the Cucurbita RNL (also known as helper) genes varied from ~ 18 to ~ 6% of the full NLR gene complement (Suppl. Fig. S1). It is a conspicuous fraction compared to that of many other Angiosperms (Andolfo et al. 2020). The helpers are essential in defense-related gene regulatory networks for increasing the robustness of the innate immune system (Wu et al. 2018). Cucurbita NLR gene complements include several genes carrying integrated domains (IDs), which act as "baits" enabling pathogen recognition (Sarris et al. 2016). NLR-ID proteins for functioning may require an appropriate helper or NLR partner and when transferred as a unit between plant families confer resistance to diverse pathogens (Nguyen et al. 2021).
The speed of adaptation to changing pathogen populations of plant species is crucial for their survival, therefore, the availability of a varied NLR defense arsenal may be useful (Andolfo and Ercolano 2015;Di Donato et al. 2017). Tandem duplication, exon shuffling, gene gains and losses have a decisive role in NLR genes diversification of Cucurbita genus (Lin et al. 2013;Andolfo et al. 2017). C. maxima and C. moschata preferentially tend to organize theirs defense arsenals in clusters, and most of them are derived from recent gene duplications. By contrast, C. pepo and C. argyrosperma genomes showed a random genome-wide NLR organization without clear gene clustering events. The cluster formation may facilitate the NLR gene variability by sequence recombination (Andolfo et al. 2013;Joshi and Nayak 2013), but the fitness cost can prevent R locus fixation. Therefore, NLR expansion in a genome could be limited to balance the energetic cost for R-gene transcription, translation and regulation (Tian et al. 2003). Comparative analyses of zucchini, cucumber, melon and watermelon reveal that the NLR lineages present in Cucurbitaceae family predate speciation (Harris et al. 2009;Andolfo et al. 2017;Román et al. 2020). The rare events of recent duplications influenced the expansion trajectories of cucurbits R-gene family (Wan et al. 2012). The plant species of Fabaceae and Rosaceae taxa, evolutionarily close to the cucurbits, displayed significantly larger NLR families (Jia et al. 2015), pointing out a meaningful impact of Cucurbitaceae clade divergence on the fate of its R-genes. In three Cucurbitaceae species (C. pepo, C. melo and C. sativus), several cucurbit-specific expansions of cell surface receptors (receptor-like kinase (RLK) and receptor-like protein (RLP) were identified, suggesting that the limited number of NLRs is complemented from the expansions of RLKs and RLPs .
In our dataset, the variation by multiple fold size of the eNLR family has been found not only within the genus but also among accessions belonging to the same species. The eNLR genes found in this study in some species resulted lower than Cucurbita repertories already reported Román et al. 2020). The NLRs are biotic-stress responsive genes, usually, constitutively expressed at a low level in the plants (McHale et al. 2006). All the eNLR genes molecularly tested were shown to be transcribed. The transcriptome-wide atlas of eNLR sequences was produced from RNA sequencing data of healthy leaf tissues of Cucurbita species and slight differences in eNLR expression levels may influence their annotation. In addition, the fragmented nature of the 40 transcriptomes may reduce the correct identification of eNLR genes. It is well known that the duplicated genes can collapse in assemblies (Bayer et al. 2018) and automated annotation systems produced approximately 30% of errors in NLR identification (Meyers et al. 2003;Andolfo et al. 2014b). In our study, only 4% of automated annotated genes showed full-length NLR structure. Therefore, we integrated a manual curation of data that enabled the discovery of NLR genes encoding single or partial NLR domains. Moreover, the R-genes tend to be masked by the genome annotation pipelines using public databases for transposable elements (Bayer et al. 2018). Indeed, our BLAST search identified several significant exon-matches for Fom-2 into the C. pepo and C. argyrosperma genomic regions devoid of predicted gene models. The observed variation of eNLR family sizes could be also accounted for R-genes losses in domesticated Cucurbita accessions. In addition, the differences between the BLASTn matches and the annotated gene models for Fom-2 could be the result of pseudogenization processes during domestication. This is consistent with the fact that 4 (out of 5) of analyzed genomes belong to domesticated taxa. Our R-gene annotation underlined a significant contraction of NLR gene family in C. argyrosperma compared to its wild relative C. sororia (Barrera-Redondo et al. 2021). The selective pressures imposed during domestication actively purge deleterious alleles, such as those involved in the defense mechanisms (Moreira et al. 2018).
The reconstruction of the evolutionary history of the Cucurbita eNLR and NLR gene family allowed the identification of three main branches (TNL, CNL and RNL). Comparing our Cucurbita NLR repertory with well-characterized plant R-genes, we noticed that the Cucurbita eNLRs collapsed in clades without reference R-genes. Six individual large clades (CNL1, CNL3-6 and TNL1) do not have similarity to any well-characterized reference R-gene, and might thus be potential sources of resistance to specific Cucurbita pathogen and pest elicitors (Andolfo and Ercolano 2015). Surprisingly, the drastic contraction of Cucurbita R family is associated to a species-specific diversification of CNL and TNL protein classes. Indeed, plant NLRs have been evolved in diverse families and subfamilies (Shao et al. 2016;Andolfo et al. 2020). Direct orthologs were identified only for Vat and Fom-2, two important CNL genes cloned in melon (Joobeur et al. 2004;Dogimont et al. 2014). Instead, the Cucurbita RNLs collapsed in two distinct clades together with the relative reference R-genes, ADR1 and NRG1 (Peart et al. 2005;Bonardi et al. 2017) confirming that RNLs are well conserved (Shao et al. 2016). The variations in Cucurbita-specific NLR subfamilies may depend on their involvement in successful stress responses during the host-pathogen interaction, resulting in their selection to the detriment of several NLR lineages. The ratio of nonsynonymous (Ka) to synonymous (Ks) nucleotide substitutions among coding DNA sequences of phylogenetic clades was below 1, except for clade CNL5. Heterogeneity of impact of selection on NLRs has been previously reported in plants (McHale et al. 2006). The NB domain of R-genes is reported to be more commonly subjected to purifying selection (Yang et al. 2006;Wan et al. 2012); the negatively selected codon sites were located on the NB-encoding region. Our results confirmed that global purifying selection act on Cucurbita TNL, CNL and RNL genes, but neuroses codon sites under positive/negative selection were identified.

Conclusion
In this study, we report the first Cucurbita-wide annotation and characterization of NLR and eNLR genes. The IDs identification, chromosomal locations, phylogenetic relationships and species gene gains/losses' analysis provided in this paper will contribute to the identification, isolation and characterization of R-genes. New insights into the TNL, CNL and RNL gene evolution and diversification of Cucurbitaceae family and were revealed. An important result for future applications was the reduction of R-gene family complexity from eNLR atlases, thus avoiding sequence analysis of non-expressed paralogues. The eNLR and NLR repertories identified into the 45 de novo assemblies may represent a useful resource for performing both function analysis of R-genes and genetic improvement.