Leaves were collected from S. pinnatisectum, and a total of 56 Gb of the 150-bp paired-end DNA reads was obtained after adapter trimming and quality filtering (Supplementary Table: Sequence and Analysis Data). The survey analysis estimated the S. pinnatisectum genome size to be 664–668 Mb with a heterozygosity of 1.49% (See Supplementary Table S1; Supplementary Figure S1). Regarding the genome assembly process, given the relatively high degree of heterozygosity present within the S. pinnatisectum genome, we implemented a filtration step to eliminate redundant contigs and those that aligned with the mitochondrial or chloroplast sequences from the nucleotide sequence database (NT). This approach ensures a more accurate and refined representation of the genomic data. The resulting genome size was 664 Mb, with the contig N50 value being 9 Mb (Table 1). We first evaluated the quality of the assembly using benchmarking universal single-copy orthologs (BUSCO), whose results demonstrated that the assembled S. pinnatisectum genome exhibited a high level of completeness, as evidenced by a gene set completeness rate of 99.38% (Table 1; Supplementary Table S3). To further validate the accuracy of our assembly, we mapped the Illumina paired-end reads back to the S. pinnatisectum genome; the mapping rate achieved was 99.82%, and the genome coverage with a read depth of at least 5× reached 98.10% (Table 1; Supplementary Table S2). Ultimately, employing high-resolution chromosome conformation capture (Hi-C) technology allowed us to anchor the genomic contigs onto 12 chromosomes, resulting in a total length of 600.10 Mb and an impressive loading rate of 90.36% (Figure 1F; Supplementary Table S4). Notably, two of these chromosomes (chr02 and chr11) were each composed of a single contiguous contig (Supplementary Table S4). Given that both NGS and third-generation sequencing (TGS) data were used, it is expected that there should not be any homozygous variants reported from the NGS data because of the inherent heterozygosity within the S. pinnatisectum genome. Consequently, any detected homozygous variants were considered errors rather than genuine biological events. In conclusion, the quality assessment revealed that the genome accuracy was Q40, indicating a precision greater than 99.99% (Table 1; Supplementary Table S4), thereby underscoring the reliability and high fidelity of the assembled S. pinnatisectum genome.

Genomic repeats were analyzed to assess the estimated genome quality that the S. pinnatisectum genome harbors 437.47 Mb of repetitive sequences, which account for a substantial proportion of 65.87% of the entire genomic content. Furthermore, long terminal repeats (LTRs) were found to be dominant among these genomic repeats, constituting 47.52% of the total content (Table 1; Supplementary Table S5). We proceeded to predict 34,245 protein-coding genes in S. pinnatisectum using a combined approach involving ab initio, transcriptome-based, and homology-based predictions. Of these predicted genes, an impressive 99.34% (equating to 34,019 genes) were successfully annotated with functional information across five distinct databases (refer to Table 1; Supplementary Table S6; Supplementary Figure S2). The BUSCO assessment revealed a high level of completeness of 97.96% for the protein sequences within the S. pinnatisectum genome (also detailed in Table 1; Supplementary Table S3). Additionally, our examination identified a comprehensive set of non-coding RNAs within the S. pinnatisectum genome, including 880 microRNAs, 1,036 transfer RNAs, 223 ribosomal RNAs, and 24 regulatory RNAs (the specifics can be found in Supplementary Table S7). The annotated results of the genes and repeat sequences are shown in Figure 1E. The density of repeats and genes also indicates that the genome sequence and annotated result are of high quality.

For the gene family analysis, three potato species were considered. There are 13,765 single-copy families in the genes of S. tuberosum L., S. chacoense, and S. pinnatisectum. These sets of single-copy family genes separately constitute 41.90%, 41.75%, and 40.20% of the total genes. A total of 5,416 multiple-copy families were identified in each of the species, with corresponding counts of 11,466 genes in S. tuberosum L., 12,749 genes in S. chacoense, and 11,651 genes in S. pinnatisectum. For the remaining families, S. tuberosum L. accounted for 3,875, S. chacoense accounted for 3,888, and S. pinnatisectum accounted for 2,555 family groups. There were 205, 492, and 427 unique family groups specific to S. tuberosum L., S. chacoense, and S. pinnatisectum, respectively. Moreover, there were certain genes that were not affiliated with any recognized gene families; thus, there were 1,664 non-family genes in S. tuberosum L., 3,712 genes in S. chacoense, and 1,078 genes in S. pinnatisectum (as shown in Figure 2A and detailed in Supplementary Table S8).

Focusing on S. pinnatisectum, we selected its unique family genes (427) and non-family genes (1,078) for functional analysis, which are highlighted in yellow in Figure 2A and enumerated in Supplementary Table S8. These genes represent the species-specific repertoire compared to the two other cultivated species. To gain deeper insights into their functional significance, we conducted enrichment analyses using KEGG pathways and GO terms for these species-specific genes in S. pinnatisectum.

The top 20 significantly enriched KEGG pathways are presented in Figure 2B; here, aside from the essential cellular processes such as ribosomes, protein processing in the endoplasmic reticulum, and the MAPK signaling pathway specific to plants, we observed that the plant–pathogen interactions and phagosome pathways showed significant enrichment of species-specific genes in the wild-type S. pinnatisectum. The most significantly enriched GO terms are visually depicted in Figure 2C, where the regulation of cellular respiration, defense responses, and negative regulation of cell population functions show significant enrichment in species-specific genes in the wild type S. pinnatisectum. The objective of this analysis was to examine the specific roles and biological pathways potentially engaged by these distinct non-family genes within S. pinnatisectum.

Based on the time-series materials, a time-series differential expression analysis was conducted using data collected at 0 h, 6 h, 9 h, and 12 h. This process resulted in the identification of 330 differentially expressed genes over time (refer to Supplementary Table: Cluster Analysis). Initially, these genes were categorized into four distinct groups (Figure 3A); the initial gene expressions in Clusters 1 and 4 were significantly lower than their counterparts in Clusters 2 and 3, which displayed notably heightened expressions. However, with the passage of time, an intriguing shift was observed in the expression dynamics; the transcriptional activities in Clusters 1 and 4 increased steadily, while there were concurrent reductions in the expression levels in Clusters 2 and 3. Subsequently, these four categories could be further divided into two main clusters (Figure 3B). Upon closer examination, it was observed that genes within Cluster 1 exhibited an ascending expression pattern, starting from low levels at 0 h and increasing to higher levels by 12 h. In contrast, genes in Cluster 2 demonstrated a descending trend, with initial high expression levels that gradually decreased over the same time period from 0 to 12 h. The KEGG and GO functional enrichment analyses disclosed that Cluster 1 manifested substantial overrepresentation in multiple biological pathways. Notably, these included the KEGG pathways for sesquiterpenoid, triterpenoid, flavonoid, and phenylpropanoid biosyntheses, as depicted in Figure 3C. Regarding the GO terms, Cluster 1 exhibited remarkable enrichment in terpene synthase activity, diterpenoid biosynthesis processes, and monooxygenase activities, emphasizing its unique biological roles and molecular functions. Additional details on this gene enrichment may be found in the Supplementary Table.

In the genome of S. pinnatisectum, 303 genes containing the motifs of the NB-ARC (NBS) domain were identified (Supplementary Table: R gene statistics). These protein sequences were aligned to construct the gene tree (Figure 4). The gene FPKM of the R genes based on 12 RNAseq are displayed in the outer ring of the gene tree. The differential expression ranges from −2 to 2 (see the FPKM label in Figure 4). The non-expression label (gray color) means that these genes (39 numbers) are not expressed in the 12 RNAseq (Figure 4). The target label (pink) means that these genes (68 numbers) show ascending expression patterns from low to high with time lapse in the 12 RNAseq (Figure 4). Conversely, the non-target label applies to 196 genes (indicated in dark blue), whose expression patterns do not follow the low-to-high trend in the 12 RNAseq dataset. For the motifs in the R gene, different motifs in the sequences were identified, and the number of motif genes is 10 at most (Supplementary Figure S3). These genes have different locations in the chromosome; 299 out of the 303 genes are located in the chromosome and the remaining 4 are located in the dispersed contigs (Supplementary Table: R gene statistics). Genes located on the chromosome are depicted in the chromosome diagram, and some genes such as chr4, chr5, and chr11 together constitute 16.17% (49 genes), 10.23% (31 genes), and 13.86% (42 genes) of the total number of genes (Figure 5).

To identify the shared and unique genes among the resistance genes and Cluster 1, gene annotation results were obtained (Supplementary Figure S4); the genes were exclusively assigned to a single group of 273 (constituting 90.10% of the total) in the R-resistance category, 119 (representing 93.70%) in Cluster 1, and 64 (comprising 62.75%) in Function. Additionally, there may exist genes that simultaneously belong to more than one subgroup. Our analysis showed that 30 genes bear functional significance and are related to R-resistance. Only eight genes (StuW01G003151, StuW01G003152, StuW01G003155, StuW02G002548, StuW04G001654, StuW05G001553, StuW10G002376, and StuW10G002377) were specifically found to have significant functionality within Cluster 1, as depicted graphically. It is worth noting that no genes were found to be common to both the R-resistance category and Cluster 1.

For genome sequencing, a single tissue culture seedling of S. pinnatisectum was collected from Jiangxi Key Laboratory of Crop Growth and Development Regulation (28^oN 115^oE) (Figures 1A–D). The total genomic DNA from the young fresh leaves of one plant was extracted using the CTAB method (Doyle and Doyle, 1987). Approximately 10 µg of the DNA was sheared into 10–50 kb fragments, followed by size selection on the BluePippin instrument. Approximately 5 μg of the recovered DNA was retrieved for library construction using the Ligation Sequencing 1D kit (SQK-LSK109, ONT, United Kingdom) according to manufacturer instructions, and the final library was sequenced on the Oxford Nanopore PromethION platform (ONT, United Kingdom) (Lu et al., 2016) at the Genome Center of Grandomics (Wuhan, China). Reads with mean quality scores higher than 7 were retained. For the Illumina library construction, the extracted DNA was fragmented and fractionated from the same source and was subjected to paired-end library construction, which was subsequently sequenced on the Illumina NovaSeq 6000 platform (Illumina Inc., CA, USA). Furthermore, leaves were collected from the same S. pinnatisectum source, and RNAseq reads were generated for genome annotation using the Illumina platform.

The Illumina paired-end reads were filtered using fastp (v0.19.6) (Chen et al., 2018) with default parameters and were then applied toward genome size and heterozygosity estimations using Jellyfish (v2.2.3) (Marcais and Kingsford, 2011). Approximately 62.59 Gb of the pass reads sequenced from the Nanopore PromethION platform by Guppy (Sherathiya et al., 2021) were obtained, and the S. pinnatisectum genome was subsequently assembled using NextDenovo software (https://github.com/Nextomics/NextDenovo) (read_cutoff = 1k, seed_cutoff = 28k). To obtain a genome with greater accuracy, error correction was performed on the assembled contigs using Racon (v1.5.0) (Vaser et al., 2017) with the Nanopore long reads and NextPolish (v1.2.4) (Hu et al., 2019) with the Illumina short reads for three and four rounds, respectively. The genome redundancies were detected and removed by Redundans (Pryszcz and Gabaldón, 2016) (with --identity 0.88 and --overlap 0.88).

To evaluate the accuracy of the genome assembly, the Illumina genomic paired-end reads were mapped to the genome contig sequences using the “mem” submodule of BWA (Jung and Han, 2022; Langarita et al., 2023). The mapping identity and genome coverage of the genome assembly were calculated from the mapping results obtained with SAMtools v1.4 (Li et al., 2009) with the default parameters. Homozygous single-base variations were subsequently detected using BCFtools v1.8.0 (Narasimhan et al., 2016; Danecek and Mccarthy, 2017) with the default parameters. Furthermore, the Illumina RNAseq reads were mapped to the genome sequence using HISAT2 v2.1 (Kim et al., 2019) with the default parameters, and the mapping rate of the RNAseq reads was calculated with SAMtools (Li et al., 2009). The completeness of the conserved genes and eukaryote core gene assembly were evaluated using BUSCO v5.1.3 (Simao et al., 2015; Seppey et al., 2019) with the “embryophyta_odb10” dataset.

To further eliminate contaminated sequences of the genome that could cause potential problems during downstream analysis, the error-corrected genome contigs were aligned with the Nucleotide Sequence Database (NT) (Harger et al., 2000) using BLASTN v2.9 (Camacho et al., 2009) with the parameter “E-value 1e-5,” and the sequence alignment results were classified based on species taxonomy. The contigs aligned to taxonomies except “Viridiplantae” and “Nohit” were classified as contamination sequences and filtered out from the genome assembly.

The Hi-C library was constructed using young fresh leaves from the same S. pinnatisectum Dunal and sequenced using the Illumina platform. An improved Hi-C procedure (Lieberman-Aiden et al., 2009; Louwers et al., 2009; Rao et al., 2014) was adapted. Briefly, the leaves were fixed with 1% formaldehyde to induce crosslinking (Sigma) and were subsequently lysed to form the cohesive ends by restriction endonuclease DPN II (NEB). The digested DNA was blunt-ended by filling the nucleotides by Klenow enzyme (NEB) with biotin-14-dATP (Invitrogen), followed by ligation by T4 DNA ligase (NEB). After overnight incubation to reverse the crosslinks, the ligated DNA was sheared into 300–600 bp fragments. The DNA fragments were blunt-end repaired and A-tailed, followed by purification through biotin-streptavidin-mediated pull down. Finally, the Hi-C library was sequenced on the Illumina NovaSeq-6000 platform (Illumina Inc., CA, USA). For chromosome-level scaffolding, Hi-C paired-end reads were filtered using fastp (v0.19.4) (Chen et al., 2018) with the default parameters and were then aligned to the decontaminated genome contigs using bowtie2 (v2.3.2) (Langmead and Salzberg, 2012) with the end-to-end model (-very-sensitive -L 30). LACHESIS (Burton et al., 2013) (https://github.com/shendurelab/LACHESIS) was subsequently applied according to the agglomerative hierarchical clustering algorithm to cluster the contigs with CLUSTER MIN RE SITES, CLUSTER MAX LINK DENSITY, CLUSTER NONINFORMATIVE RATIO, ORDER MIN N RES IN TRUNK, and ORDER MIN N RES IN SHREDS set to 100, 2.5, 1.4, 60, and 60, respectively, to assemble the genome contigs into groups that were further ordered and oriented into chromosomes. Finally, the chromosome-level genome was revised manually based on the heat-map matrix of Hi-C.

Tandem repeats (TRs) across the S. pinnatisectum genome were annotated using GMATA (v2.2) (Wang and Wang, 2016) with the default parameters and Tandem Repeats Finder (TRF) (v4.07b) (Benson, 1999) (2 7 7 80 10 50 500 -f -d -h -r). The plant transposable elements (TEs) were searched separately using LTR_finder (v1.0.6) (Xu and Wang, 2007) and LTR_harvest (v1.6.5) (Ellinghaus et al., 2008) with the default parameters, and their results were applied to construct an LTR library file using LTR_retriever (Ou and Jiang, 2018). A MITE transposon library was generated using MITE-hunter (Han and Wessler, 2010) (-n 20 -P 0.2 -c 3) for plants and animals, and a de novo TE library was predicted using RepeatModeler (v1.0.11) (Bedell et al., 2000) (-engine wublast). The LTR, MITE transposon, and de novo TE libraries as well as Repbase database (Jurka et al., 2005) were combined to construct the final TE library for S. pinnatisectum, which was then used as the repeat library for RepeatMasker (Bedell et al., 2000) (v4.0.6; www.repeatmasker.org) (nolow -no_is -gff -norna -engine abblast -lib lib) to identify the TE elements in the S. pinnatisectum genome. The results of the TRs and TEs were merged and masked from the genome sequence. Finally, further repetitive sequences in the masked genome were found using RepeatProteinMask.

Gene models of the S. pinnatisectum genome were constructed by ab initio, homology-based, and transcriptome-based predictions. The RNAseq paired-end reads were mapped to the S. pinnatisectum genome using HISAT2 (v2.1.0) (Kim et al., 2019), and StringTie (v1.3.3) (Pertea et al., 2015; Pertea et al., 2016) with the default parameters was subsequently applied to assemble the transcripts that were then used as the inputs to PASA (v2.0.2) (Haas et al., 2003) (--ALIGNERS gmap, blat) for transcriptome-based gene prediction. Augustus (v3.3.1) (Stanke et al., 2006) was used for the ab initio gene prediction with default parameters. Moreover, the protein sequences of six homologous species (Nicotiana attenuata, Arabidopsis thaliana, S. aethiopicum, S. pennellii, Caosicum baccatum, and S. chacoense) were used for the homology-based prediction through GeMoMa (v1.5.3) (Keilwagen et al., 2019) with default parameters. The gene prediction results from the three methods were integrated using EVidenceModeler (EVM; v1.1.1) (Haas et al., 2008) to obtain the raw gene set. Finally, the genes whose sequences were composed of TEs were filtered using TransposonPSI (Urasaki et al., 2017) (http://transposonpsi.sourceforge.net). The completeness of the predicted genes was evaluated using BUSCO v5.1.3 (Simao et al., 2015; Seppey et al., 2019) with the “embryophyta_odb10” dataset.

The annotation of the non-coding RNA set was performed next, and the genome of S. pinnatisectum was aligned to the non-coding database Rfam (v11.0) (Griffiths-Jones et al., 2005) using INFERNAL (Nawrocki et al., 2009) to annotate the genes of the small nuclear RNAs (snRNAs) and microRNAs (miRNAs). Transfer RNAs (tRNAs) were predicted using tRNAscan-SE (v1.3.1) (Griffiths-Jones et al., 2005). Finally, the ribosome RNAs (rRNAs) were predicted using RNAmmer (v1.2) (Lagesen et al., 2007).

The biological functions of the predicted genes in the S. pinnatisectum genome were annotated using two strategies with protein sequences. First, the predicted protein sequences were aligned with the Swiss-Prot protein database (Bairoch and Apweiler, 2000; Boutet et al., 2007), non-redundant protein sequence database (NR) (Harger et al., 2000), Kyoto Encyclopedia of Gene and Genomes (KEGG) database (Kanehisa and Goto, 2000; Kanehisa et al., 2017), and Eukaryotic Orthologous Groups (KOG) of protein database (Tatusov et al., 2003) using BLASTP (v2.7.1) (Camacho et al., 2009) with parameters “-evalue 1e-5, -max_target_seqs 1.” The Gene Ontology (GO) (Ashburner et al., 2000) analysis was subsequently performed using InterProScan (Zdobnov and Apweiler, 2001) v5.32-71.0 with the default parameters and databases. For the circles, all proteins were aligned against each other using BLASTP (-e 0.01), and Python package MCScanX (Wang et al., 2012) was used to find the collinear segments based on the protein alignment files. Subsequently, the number of genes in the collinear block was found to be more than 40, which was retained in Figure 1E.

The RNAseq data of S. pinnatisectum were obtained from NCBI with the BioProject accession number PRJNA616420 (S. pinnatisectum sample accession numbers: SRX8168235, SRX8168240, SRX8168241, SRX8168242, SRX8168243, SRX8168244, SRX8168245, SRX8168246, SRX8168247, SRX8168248, SRX8168249, and SRX8168250) (Gu et al., 2020). These data were based on histological observations of infected leaf tissues 0, 6, 9, and 12 hours post inoculation (hpi) as the time points to investigate the transcriptional dynamics of S. pinnatisectum. Quality control of the RNAseq reads were performed using fastp (Chen et al., 2018) with the default parameters and mapped to the S. pinnatisectum genome sequence using HISAT2 (v2.1.0) (Pertea et al., 2016; Kim et al., 2019) with the default parameters. Read alignments for the transcripts in each sample were extracted and counted using StringTie (v1.3.3) (Pertea et al., 2015; Pertea et al., 2016). The expression level of each gene was measured in terms of the fragments per kilobase per million (FPKM) values estimated in StringTie. The read count of each gene generated by StringTie script was used for differential expression analysis. DESeq2 (Love et al., 2014) was employed in this analysis with false discovery rate (FDR) ≤0.05 and fold change ≥2. Because the sample of RNAseq is time ordered, the genes with significantly different expressions in the time series were selected to determine the gene clusters using pheatmap with kmeans. Two and four clusters (k = 2 or k = 4) were implemented, and two clusters were found to be more suitable by comparing the two results. Thus, the gene set with low to high expressions from 0 to 12 hpi was selected for enrichment analysis through GO and KEGG using the R-package of clusterProfiler (Yu et al., 2012; Wu et al., 2021). Finally, ggplot2 (Wickham, 2009) was used to draw the bubble diagram.

Diploid cultivated potatoes S. tuberosum L. (http://spuddb.uga.edu/dm_v6_1_download.shtml) and S. chacoense (http://spuddb.uga.edu/M6_v5_0_download.shtml) were downloaded for the downstream analyses. In these downloaded and S. pinnatisectum data, the longest mRNA in each gene was selected from the annotated file and translated to a protein sequence. Then, all selected protein sequences were aligned and clustered using Orthofinder v2.5.4 (Emms and Kelly, 2019) (-S diamond). Next, the clustered gene family was classified into four groups. The single-copy family showed that there was only one gene from each species in this group. Multiple families would indicate that the number of genes, which could be from any one species, is equal or greater than one and is not same as the single-copy family in this group. A unique family indicates that the number of genes is greater than one for any single species and that the others are zero. The family that does not belong to the single-copy, multiple, or unique categories is called as the other family. Meanwhile, there may still be some genes that may not belong to the above gene families. According to the above classification method, all genes were categorized into one of these five groups as four gene families and one non-family. For each species, the genes of the unique family and non-family are the species-specific genes. Thus, the species-specific genes in S. pinnatisectum were collected for gene set enrichment analysis, including GO and KEGG, using the R-package of clusterProfiler. Finally, ggplot2 was used to draw the bubble diagram.

The hidden Markov model (HMM) file of the NB-ARC (NBS) domain (PF00931: http://pfam-legacy.xfam.org/search/keyword?query=PF00931) (Pandolfi et al., 2017) was download from the pfam database (http://pfam-legacy.xfam.org/) (Mistry et al., 2021). Then, the predicted protein sequences from the S. pinnatisectum genome were first aligned against the HMM of PF00931 using hmmscan in HMMER v3.3.2 (Potter et al., 2018). Next, from the raw aligned results, a high-quality protein set (E-value < 1×e⁻²⁰) was selected to construct the species-specific HMM file using hmmbuild in HMMER v3.3.2 (Potter et al., 2018). Then, the predicted protein sequences from the S. pinnatisectum genome were aligned again with the species-specific NBS HMM using hmmscan in HMMER v3.3.2 (Potter et al., 2018). The genes with E-values less than 0.01 were obtained for the newly aligned file as the final resistance-related (R) genes (Lozano et al., 2015).

For NBS encoding proteins, the sequences were aligned using mafft v7.471 (Katoh et al., 2002) with default parameters, and the alignment file was subsequently input into Gblock v0.19b (Talavera and Castresana, 2007) (-t = p -b5 = h) to product the trimmed alignment file. Then, iqtree v2.2.0 (Nguyen et al., 2015) (-b 1000) was used to construct the gene tree with the trimmed alignment file. Based on the transcriptome data analysis, the FPKM of the R genes was obtained and then modified in the tree using iTOL (https://itol.embl.de/) (Letunic and Bork, 2007; Letunic and Bork, 2016).