# Cultivar-specific transcriptome and pan-transcriptome reconstruction of tetraploid potato

Marko Petek, Maja Zagorščak, Živa Ramšak, Sheri Sanders, Špela Tomaž,

Elizabeth Tseng, Mohamed Zouine, Anna Coll, Kristina Gruden

► To cite this version:

Marko Petek, Maja Zagorščak, Živa Ramšak, Sheri Sanders, Špela Tomaž, et al. Cultivar-specific transcriptome and pan-transcriptome reconstruction of tetraploid potato. Scientific Data, 2020, 7, pp.249. 10.1038/s41597-020-00581-4. hal-02913152

HAL Id: hal-02913152

https://hal.inrae.fr/hal-02913152

Submitted on 7 Aug 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
SCIENTIFIC DATA

<figure><img src="image_2.png" /></figure>

Check for updates

OPEN

DATA DESCRIPTOR

# Cultivar-specific transcriptome and pan-transcriptome reconstruction of tetraploid potato

Marko Petek<sup>1,6✉</sup>, Maja Zagorščak<sup>1,6✉</sup>, Živa Ramšak<sup>1</sup>, Sheri Sanders<sup>3</sup>, Špela Tomaž<sup>1,2</sup>, Elizabeth Tseng<sup>4</sup>, Mohamed Zouine<sup>1,5</sup>, Anna Coll<sup>1</sup> & Kristina Gruden<sup>1</sup>

Although the reference genome of *Solanum tuberosum* Group Phureja double-monoploid (DM) clone is available, knowledge on the genetic diversity of the highly heterozygous tetraploid Group Tuberosum, representing most cultivated varieties, remains largely unexplored. This lack of knowledge hinders further progress in potato research. In conducted investigation, we first merged and manually curated the two existing partially-overlapping DM genome-based gene models, creating a union of genes in Phureja scaffold. Next, we compiled available and newly generated RNA-Seq datasets (cca. 1.5 billion reads) for three tetraploid potato genotypes (cultivar Désirée, cultivar Rywal, and breeding clone PW363) with diverse breeding pedigrees. Short-read transcriptomes were assembled using several *de novo* assemblers under different settings to test for optimal outcome. For cultivar Rywal, PacBio Iso-Seq full-length transcriptome sequencing was also performed. EvidentialGene redundancy-reducing pipeline complemented with in-house developed scripts was employed to produce accurate and complete cultivar-specific transcriptomes, as well as to attain the pan-transcriptome. The generated transcriptomes and pan-transcriptome represent a valuable resource for potato gene variability exploration, high-throughput omics analyses, and breeding programmes.

## Background & Summary

At species level, genomes of individuals can differ in single nucleotide polymorphisms (SNPs), short insertions and deletions (INDELs), gene copy numbers, and presence or absence of genes<fnref n="1" />. The latter leads to the concept of species specific pan-genomes, namely the core genome present in most individuals and the dispensable genome comprised of genes present only in a subset of individuals, which results in the emergence of particular subgroup-specific phenotypes. This concept has been extended to pan-transcriptomes, where the presence or absence of variations is not bound only to gene content, but also to genetic and epigenetic regulatory elements. Pan-genomes and pan-transcriptomes have been described in the model plant species *Arabidopsis thaliana*² and several crop species, including maize³,⁴, rice⁵, wheat⁶ and soybean⁷.

While the genome of a double-monoploid clone of *Solanum tuberosum* Group Phureja (DM) is available⁸, this diploid potato group differs from the tetraploid Group Tuberosum, which includes most varieties of cultivated potato. Through domestication and modern breeding efforts, different potato cultivars also acquired genes from other closely related Solanum species or lost some ancestral genes<fnref n="1" />. Different breeding programmes have resulted in accumulation of smaller genome modifications, e.g. SNPs and INDELs. Consequently, each distinct potato cultivar harbours a unique set of transcripts, resulting in physiological and developmental differences and different responses to biotic and abiotic stress. SNP and INDEL profile differences and novel gene variants in anthocyanin pathway were identified in a comparative transcriptome analysis of two Chinese potato cultivars⁹. Unfortunately, we could not include these transcriptomes in our pan-transcriptome because the assemblies were not publicly accessible.

Based on the DM genome, the PGSC and ITAG annotation consortia⁸,¹⁰ have each independently produced potato gene models. For practical reasons, most potato researchers use only one genome annotation, either PGSC or ITAG, especially when conducting high-throughput analyses. Using an incomplete gene set can lead to false

<footnote n="1">¹Department of Biotechnology and Systems Biology, National Institute of Biology, Ljubljana, Slovenia. ²Jožef Stefan International Postgraduate School, Ljubljana, Slovenia. ³National Center for Genome Analysis and Support (NCGAS), Indiana University, Bloomington, USA. ⁴PacBio, Menlo Park, CA, USA. ⁵Laboratoire Génomique et Biotechnologie des Fruits, INRA-INP/ENSAT, Castanet-Tolosan, France. ⁶These authors contributed equally: Marko Petek, Maja Zagorščak. ✉e-mail: marko.petek@nib.si; maja.zagorscak@nib.si</footnote>
conclusions regarding gene presence or gene family diversity in potato. Using a computational pipeline followed by manual curation, we have consolidated the two publicly available Group Phureja DM gene model sets to produce a unified one.

While a combined DM gene set is useful, it is still not as useful as a pan-transcriptome that includes assemblies from cultivated potatoes. However, obtaining an optimal transcriptome from short-read RNA-Seq data is not a trivial task. Each *de novo* assembler suffers from different intrinsic error generation and no single assembler performs best on all datasets¹¹. To maximise diversity and completeness of potato cultivar transcriptomes, usage of multiple *de novo* transcriptome assemblers and various parameter combinations over the same input data was employed. Following this"over-assembly" step, we used tr2acds pipeline from EvidentialGene¹² to reduce redundancy across assemblies and obtain cultivar-specific transcriptomes. Finally, we consolidated representative cultivar-specific sequences to generate a potato pan-transcriptome (StPanTr). These transcriptomes will improve high throughput sequencing analyses, from RNA-Seq and sRNA-Seq to more specific ones like ATAC-Seq, by providing a more comprehensive and accurate mapping reference. The translated protein sequences can enhance the sensitivity of high-throughput mass-spectroscopy based proteomics. The resource is valuable also for the design of various PCR assays, e.g. quantitative PCR, where exact sequence information is required. Additionally, the knowledge generated regarding variations in transcript sequences between cultivars, such as SNPs, insertions and deletions, will be a key instrument in assisting potato breeding programmes.

## Methods

**Merging PGSC and ITAG gene models of reference genome Group Phureja.** GFF files corresponding to their respective gene models (PGSC v4.04, ITAG v1.0) were retrieved from the Spud DB (solanaceae, plantbiology.msu.edu) potato genomics resource¹³. The two models (39,431 PGSC and 35,004 ITAG) were then compared on the basis of their exact chromosomal location and orientation in order to create the most complete set of genes from both predicted genome models. Several combinations arose from the merge (Fig. 1a), those for which no action was required (singletons, i.e. sole PGSC or ITAG genes not covering any other genes); 1-to-1 or 1-to-2 combinations between PGSC and ITAG genes, which were solved programmatically; and lastly, combinations of more than 3 genes in various combination types, which continued on to manual curation. The latter (318 gene clusters; example in Fig. 1b) were considered to be non trivial merge examples (overlapping genes in two models or multiple genes in PGSC corresponding to a single gene in ITAG). This resulted in a merged DM genome GFF3 file with 49,322 chromosome position specific sequences, of which 31,442 were assigned with ITAG gene IDs and 17,880 with PGSC gene IDs¹⁴.

**Data pre-processing.** The complete bioinformatic pipeline is outlined in Fig. 2. Sequence quality assessment of raw RNA-Seq data, quality trimming, and removal of adapter sequences and polyA tails was performed using CLC Genomics Workbench v6.5-v10.0.1 (Qiagen) with maximum error probability threshold set to 0.01 (Phred quality score 20) and no ambiguous nucleotides allowed. Minimal trimmed sequence length allowed was set to 15 bp while maximum up to 1 kb. Orphaned reads were re-assigned as single-end (SE) reads. Processed reads were pooled into cultivar datasets as properly paired-end (PE) reads or SE reads per cultivar per sequencing platform. For the Velvet assembler, SOLiD reads were converted into double encoding reads using perl script "denovo_preprocessor_solid_v2.2.1.pl"¹⁵. To reduce the size of cv. Désirée and cv. Rywal datasets, digital normalization was performed using khmer from bbmap suite v37.68¹⁶ prior to conducting *de novo* assembly using Velvet and rnaSPAdes.

PacBio long reads were processed for each sample independently using Iso-Seq 3 analysis software (Pacific Biosciences). Briefly, the pipeline included Circular Consensus Sequence (CCS) generation, full-length reads identification ("classify" step), clustering isoforms ("cluster" step) and "polishing" step using Arrow consensus algorithm. Only high-quality full-length PacBio isoforms were used as input for further steps.

**PacBio Cupcake ToFU pipeline.** Cupcake ToFU (github.com/Magdoll/cDNA_Cupcake) scripts¹⁷ were used to further refine the Iso-Seq transcript set. Redundant PacBio isoforms were collapsed with "collapse_isoforms_by_sam.py" and counts were obtained with "get_abundance_post_collapse.py". Isoforms with less than two supporting counts were filtered using "filter_by_count.py" and 5'-degraded isoforms were filtered using "filter_away_subset.py". Isoforms from the two samples were combined into one non-redundant Iso-Seq transcript set using "chain Samples.py".

**De Bruijn graph based de novo assembly of short reads.** Short reads were *de novo* assembled using Trinity v2013-02-25¹⁸, Velvet/Oases v. 1.2.10¹⁹, rnaSPAdes v.3.11.1²⁰ and CLC Genomics Workbench v8.5.4-v10.1.1 (Qiagen). Illumina and SOLiD reads were assembled separately. For CLC Genomics *de novo* assemblies, combinations of three bubble sizes and 14 k-mer sizes were tested on PW363 Illumina dataset. Varying bubble size length did not influence the assembly statistics much (Supplementary Fig. 2), therefore we decided to use the length of 85 bp for Illumina datasets of the other two cultivars. Bubble size and k-mer length parameters used for Velvet and CLC are given in Table 1. The scaffolding option in CLC and Velvet was disabled. More detailed information per assembly is provided in Auxiliary Table 2²¹.

**Decreasing redundancy of assemblies and annotation.** 739 mio Désirée short reads were assembled into 3,765,661 potential transcripts, 284 mio PW363 short reads were assembled into 6,022,291 potential transcripts, and 710 mio Rywal short reads and 1.4 mio Rywal PacBio sequences were assembled into 1,912,821 potential transcripts. While generation of several transcriptomes from diverse input data and various parameter combinations increases the likelihood of capturing and accurately assembling transcripts²², redundancy reduction without loss of information and error removal from the over-assemblies is required. All cultivar-specific transcriptome assemblies, compiled into cultivar-specific transcriptome over-assembly, were initially filtered with the tr2acds pipeline (part
<figure><img src="image_3.png" /><figcaption>**Fig. 1** Merging DM Phureja PGSC and ITAG gene models. (a) Decision tree for the merging of both genome models, with 6 possible outcomes: singleton genes ('keep gene as is'), manual curation ('gene cluster to manual curation' and programmatic solution (all remaining 4 options). Green solid lines represent a satisfied condition (Y: Yes), dashed red lines, an unsatisfied condition (N: No). (b) Example of manual curation in the merged DM genome GTF, region visualisation (chr12:11405699..11418575) in the Spud DB (solanaceae.plantbiology.msu.edu) Genome Browser¹³. ITAG defined Sotub12g014200.1.1 spans three PGSC defined coding sequences (PGSC0003DMT400005728, PGSC0003DMT400005745 and PGSC0003DMT400005726). Below the gene models, RNA sequence tracks are shown, showing how these genes are expressed in various plant organs. In the concrete case, Sotub12g014200.1.1 was preferred due to RNA-Seq evidence being in concordance, and no evidence for PGSC0003DMT400005745.</figcaption></figure>

of EvidentialGene v2016.07.11¹²) which consists of four steps. First, all perfect redundant nucleotide sequences are removed using fastnrdb, part of the exonerate package²³, leaving only the transcript with the longest coding region. Next, all perfect fragments of the remaining transcripts are removed using cd-hit²⁴. These first two steps are important in reducing transcriptome redundancy, as true transcripts are expected to be assembled independently by multiple of the assembly methods. Keeping the transcripts with the longest and most complete coding region helps eliminate chimeric and misassembled transcripts, as these errors tend to occur more often in UTR regions or in a manner that causes frameshifts and long, incomplete coding regions¹²).

The third and the fourth step of the tr2acds pipeline segregate transcripts that are likely isoforms, alleles, or other variations seen at a single locus. This is done through amino acid sequence clustering, which identifies putative transcripts that vary only in silent mutations, and through reciprocal BLAST, which detects high-identity exon-sized alignments (likely isoforms). A tag is assigned to all transcripts providing detailed information on why they were discarded (e.g. perfect fragments, perfect duplicates, very high similarity, ...) or why they were marked
<figure><img src="image_4.png" /><figcaption>**Fig. 2** Bioinformatics pipeline for generation of potato transcriptomes. Software used in specific steps are given in bold. Input datasets (sequence reads) and output data (transcriptomes) are depicted as blue cylinders. Data upload steps to public repositories are shaded in orange. Abbreviations: SRA – NCBI Sequence Read Archive, PGSC – Potato Genome Sequencing Consortium, ITAG – international Tomato Annotation Group, CLC – CLC Genomics Workbench, PacBio – Pacific Biosciences Iso-Seq sequencing, Tr – transcriptome, StPanTr – potato pan-transcriptome, tr2aacds – “transcript to amino acid coding sequence” Perl script from EvidentialGene pipeline.</figcaption></figure>

as alternatives (and which sequence they are an alternative form of). The final output of the tr2aacds pipeline are three sets of data – a non-redundant set of representative sequences (i.e. main set), a set of putative alternatives mapped to the representative set (i.e. alt set), and a discarded set (i.e. drop set) of redundant sequences. It is important to note that not all dropped sequences are of poor quality or incorrect – many of them are dropped due to full or partial redundancy.

Representative and alternative sets (termed okay sets by EvidentialGene) were merged into initial cultivar reference transcriptomes and, as tr2aacds only uses internal evidence for data curation, used in further external evidence for assembly validation, filtering and annotation steps (Fig. 2). The *de novo* cultivar-specific transcript sets were first mapped to the DM reference genome by STARlong 2.6.1d²⁵ using parameters optimized for *de novo* transcriptome datasets (all scripts are deposited at FAIRDOMHub (fairdomhub.org) project home page²⁶). Aligned transcripts were analysed with MatchAnnot to identify transcripts that match the PGSC or ITAG gene models. Domains were assigned to the polypeptide dataset using InterProScan software package v5.37-71.0²⁷. For all transcripts and coding sequences, annotations using DIAMOND v0.9.24.125²⁸ were generated by que-rying UniProt (www.uniprot.org) retrieved databases (E-value cut-off 10⁻⁰.⁵ and query transcript/cds and tar-get sequence alignment coverage higher or equal to 50%; custom databases: *Solanum tuberosum*, *Solanaceae*, plants). Initially assembled transcriptomes were also screened for nucleic acid sequences that may be of vector origin (vector segment contamination) using VecScreen plus taxonomy program v.0.16²⁹ against NCBI UniVec Database (ftp.ncbi.nlm.nih.gov/pub/UniVec). Potential biological and artificial contamination was identified as up to 3.3% of sequences per cultivar, if taking into account cases when potential contaminants covered less than 1% of the sequence (number of sequences with strong, moderate and weak proof of contamination as follows: 182, 547 and 10,509 for Désirée; 48, 228 and 7,877 for PW363; 169, 179 and 4,103 for Rywal). The results from
<table><caption>Table 1. Parameters used for short read *de novo* assembly generation. NAp – not applicable. NA – not available.</caption><thead><tr><th>Genotype</th><th>Assembly ID</th><th>Read type</th><th>Assembler</th><th>Assembler version</th><th>k-mer length (word size)</th><th>Bubble size</th></tr></thead><tbody><tr><td>Désirée</td><td>CLCdnDe8</td><td>SOLiD</td><td>CLC de novo</td><td>9.1</td><td>24</td><td>50</td></tr><tr><td>Désirée</td><td>CLCdnDe1</td><td>SOLiD</td><td>CLC de novo -</td><td>10.0.1</td><td>24</td><td>50</td></tr><tr><td></td><td></td><td></td><td>- transcript discovery as reference</td><td></td><td></td><td></td></tr><tr><td>Désirée</td><td>VdnDe8, ...</td><td>SOLiD</td><td>Velvet/Oases</td><td>1.2.10</td><td>23, 33, 43</td><td>Default</td></tr><tr><td></td><td>... VdnDe10</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Désirée</td><td>CLCdnDe9, ...</td><td>Illumina</td><td>CLC de novo</td><td>9.1</td><td>21, 23, 33, 43, 53, 63</td><td>85</td></tr><tr><td></td><td>... CLCdnDe14</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Désirée</td><td>CLCdnDe2, ...</td><td>Illumina</td><td>CLC de novo -</td><td>10.0.1</td><td>21, 23, 33, 43, 53, 63</td><td>85</td></tr><tr><td></td><td>... CLCdnDe7</td><td></td><td>- transcript discovery as reference</td><td></td><td></td><td></td></tr><tr><td>Désirée</td><td>TDe</td><td>Illumina</td><td>Trinity</td><td>r2013-02-25</td><td>25</td><td>NA</td></tr><tr><td>Désirée</td><td>VdnDe1, ...</td><td>Illumina</td><td>Velvet/Oases</td><td>1.2.10</td><td>23, 33, 43, 53, 63, 73, 83</td><td>Default</td></tr><tr><td>PW363</td><td>CLCdnPW1</td><td>SOLiD</td><td>CLC de novo</td><td>8.5.4</td><td>24</td><td>50</td></tr><tr><td>PW363</td><td>CLCdnPW2</td><td>SOLiD</td><td>CLC de novo -</td><td>9.1</td><td>24</td><td>50</td></tr><tr><td></td><td></td><td></td><td>- transcript discovery as reference</td><td></td><td></td><td></td></tr><tr><td>PW363</td><td>VdnPW8, ...</td><td>SOLiD</td><td>Velvet/Oases</td><td>1.2.10</td><td>23, 33, 43</td><td>Default</td></tr><tr><td></td><td>... VdnPW10</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>PW363</td><td>CLCdnPW3, ...</td><td>Illumina</td><td>CLC de novo</td><td>8.5.4</td><td>21, 23, 24, 25, 30, 33, 35, 40, 43, 45, 50, 53, 55, 63</td><td>50, 65, 85</td></tr><tr><td></td><td>... CLCdnPW44</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>PW363</td><td>CLCdnPW45, ...</td><td>Illumina</td><td>CLC de novo -</td><td>10.0.1</td><td>21, 23, 33, 43, 53, 63</td><td>85</td></tr><tr><td></td><td>... CLCdnPW50</td><td></td><td>- transcript discovery as reference</td><td></td><td></td><td></td></tr><tr><td>PW363</td><td>SdnPW1</td><td>Illumina</td><td>rnaSPAdes</td><td>3.11.1</td><td>43</td><td>Default</td></tr><tr><td>PW363</td><td>TPW</td><td>Illumina</td><td>Trinity</td><td>r2013-02-25</td><td>25</td><td>NA</td></tr><tr><td>PW363</td><td>VdnPW1, ...</td><td>Illumina</td><td>Velvet/Oases</td><td>1.2.10</td><td>23, 33, 43, 53, 63, 73, 83</td><td>Default</td></tr><tr><td></td><td>... VdnPW7</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Rywal</td><td>PBdnRY1</td><td>PacBio Isoseq</td><td>Iso-Seq. 3, Cupcake ToFU</td><td>2017</td><td>NAp</td><td>NAp</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Rywal</td><td>CLCdnRY1, ...</td><td>Illumina</td><td>CLC de novo</td><td>9.1</td><td>21, 23, 33, 43, 53, 63</td><td>85</td></tr><tr><td></td><td>... CLCdnRY6</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Rywal</td><td>CLCdnRY7, ...</td><td>Illumina</td><td>CLC de novo -</td><td>10.1.1</td><td>21, 23, 33, 43, 53, 63, 73, 83</td><td>85</td></tr><tr><td></td><td>... CLCdnRY12</td><td></td><td>- transcript discovery as reference</td><td></td><td></td><td></td></tr><tr><td>Rywal</td><td>SdnRY1</td><td>Illumina</td><td>rnaSPAdes</td><td>3.11.1</td><td>43</td><td>Default</td></tr><tr><td>Rywal</td><td>VdnRY1, ...</td><td>Illumina</td><td>Velvet/Oases</td><td>1.2.10</td><td>23, 33, 43, 53, 63, 73, 83</td><td>Default</td></tr><tr><td></td><td>... VdnRY7</td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

MatchAnnot, InterProScan and DIAMOND were used as biological evidence in further filtering by in-house R scripts. Transcripts that did not map to the genome nor had any significant hits in either InterPro (www.ebi.ac.uk/interpro) or UniProt (www.uniprot.org) were eliminated from further analysis to obtain higher reliability of constructed transcriptomes³⁰⁻³². Pajek v5.08³³, in-house scripts, and cdhit-2d from the CD-HIT package v4.6²⁴ were used to re-assign post-filtering representative and alternative classes and to obtain finalised cultivar-specific transcriptomes.

The whole redundancy removal procedure reduced the initial transcriptome assemblies by 18-fold for Désirée, 38-fold for Rywal, and 24-fold for PW363. Completeness of each initial *de novo* assembly to cultivar-specific transcriptome was estimated with BUSCO (Supplementary Figs. 1–3).

Individual contributions by various assembly methods were investigated in light of what contributed to the final, cleaned cultivar transcriptomes. SOLiD assemblies (Supplementary Fig. 1: CLCdnDe1, CLCdnDe8, VdnDe8-10), produced by either CLC or Velvet/Oases pipelines, contributed least to transcriptomes, which can mostly be attributed to short length of the input sequences. Interestingly, increasing k-mer size in the CLC pipeline for Illumina assemblies produced more complete assemblies according to BUSCO scores and more transcripts were selected for the initial transcriptome (Supplementary Fig. 1: CLCdnDe1-7, CLCdnDe9-14). On the contrary, increasing k-mer length in Velvet/Oases pipeline lead to transcripts that were less favoured by the redundancy removal procedure (Supplementary Fig. 1: VdnDe1-7). The Trinity assembly was comparable
<figure><img src="image_5.png" /><figcaption>**Fig. 3** Structure of the potato pan-transcriptome. Stacked bar plot showing the overlap of paralogue groups in cultivar-specific transcriptomes and merged Phureja DM gene model. Only representative and alternative transcripts of the pan-transcriptome are counted (i.e. cultivar representative sequences) while disregarding additional cultivar alternative transcripts. For Phureja DM, the merged ITAG and PGSC DM gene models were counted. DM and at least one Group Tuberosum: sequences shared by Phureja DM and at least one tetraploid genotype, core: sequences shared among all genotypes in the pan-transcriptome.</figcaption></figure>

to the high k-mer CLC assemblies in transcriptome contribution and BUSCO score (Supplementary Fig. 1). It might seem that the PacBio Iso-Seq transcripts contributed less than expected to the cv. Rywal transcriptome (Supplementary Fig. 3), however it should be noted that a considerable number of PacBio transcripts was assigned to the EvidentialGene drop set because they had perfect or near-perfect CDS identity of transcripts assembled by CLC. The EvidentialGene pipeline also prioritised CLC-assembled transcripts over PacBio transcripts because the redundancy removal algorithm reorders the near-perfect duplicates by transcript name and only retains the first transcript listed (Auxiliary Table 6³⁴).

**Potato pan-transcriptome construction.** While the PGSC gene model defined transcripts as well as coding sequences, the ITAG gene model defined only coding sequences. Therefore, the potato pan-transcriptome construction was conducted at the level of CDS.

Cultivar-specific representative coding sequences (57,943 of Désirée, 43,883 of PW363 and 36,336 of Rywal) were combined with coding sequences from the merged Phureja DM gene models (17,880 and 31,442 non-redundant PGSC and ITAG genes, respectively) and subjected to the cdhit-est²⁴ algorithm (global sequence identity threshold 90%, alignment coverage for the shorter sequence 75%, bandwidth of alignment 51 nt and word length of 9) to create potato pan-transcriptome. Sequences that did not cluster using cdhit-est were separated into tetraploid and DM datasets and subjected to the cdhit-2d²⁴ algorithm (local sequence identity threshold 90%, alignment coverage for the shorter sequence 45%, bandwidth of alignment 45 nt and word length of 5).

Sequences that are shared by the DM merged gene model and *de novo* assembled cultivar-specific transcriptomes were designated as "core" transcripts, and sequences that were assembled in only one transcriptome were designated "genotype-specific". The total pan-transcriptome includes 96,886 representative, non-redundant transcripts and 90,618 alternative sequences (covering alternative splice forms, allelic isoforms and partial transcripts) for those loci (Fig. 3, Supplementary Fig. 4, Auxiliary Table 7³⁵). The core subset of the pan-transcriptome contains 68,708 sequences, among which 12% are partial sequences.

Polyploid crop pan-genomes generally consist of many cultivar-specific genes³⁶. Therefore we included all genotype-specific sequences in our potato pan-transcriptome (Fig. 3, Supplementary Fig. 4). This subset contains 64,529 sequences, among which 13% sequences are partial³⁵. Genotype-specific transcripts are generally shorter in length than the core transcripts, however they do not differ much in the percentage of complete transcripts.
<table><caption>**Table 2.** Samples used to generate the *de novo* transcriptome assemblies. <sup>a</sup>PVY, *Potato virus Y*; CPB, Colorado potato beetle. <sup>b</sup>PE, paired-end library (the number stands for read length in nt); SE, single-end library (the number stands for read length in nt); DSN-normalized, RNA-Seq library utilizing the crab duplex nuclease; CCS, circular consensus sequences. <sup>c</sup>For paired-end libraries, pairs are counted as two reads.</caption><thead><tr><th>Genotype</th><th>Sample description<sup>a</sup></th><th>Sequencing platform</th><th>Library type<sup>b</sup></th><th>Number of reads<sup>c</sup></th><th>SRA ID</th></tr></thead><tbody><tr><td rowspan="2">Désirée</td><td rowspan="2">PVY inoculated leaves</td><td rowspan="2">Illumina</td><td>DSN-normalized</td><td>~54 mio</td><td>SRR10070125</td></tr><tr><td>PE90 unstranded</td><td></td><td></td></tr><tr><td>Désirée</td><td>non-transformed and PVY-inoculated plants, non-infested and CPB infested leaves</td><td>Illumina</td><td>PE90 unstranded</td><td>~195 mio</td><td>SRR1207287, ..., SRR1207290</td></tr><tr><td>Désirée</td><td>mock and PVY inoculated leaves and stem</td><td>SOLiD</td><td>SE50 unstranded</td><td>~154 mio</td><td>SRR10065428, SRR10065429</td></tr><tr><td>Désirée</td><td>leaves</td><td>Illumina</td><td>SE50 unstranded</td><td>~172 mio</td><td>SRR3161991, SRR3161995, SRR3161999, SRR3162003, SRR3162007, SRR3162011, SRR3162015, SRR3162019, SRR3162023, SRR3162027, SRR3162031, SRR3162035</td></tr><tr><td>Désirée</td><td>seedlings</td><td>Illumina</td><td>SE100 unstranded</td><td>~80 mio</td><td>SRR4125238, ..., SRR4125247</td></tr><tr><td>Désirée</td><td>roots</td><td>Illumina</td><td>SE100 unstranded</td><td>~31 mio</td><td>SRR4125248, ..., SRR4125252</td></tr><tr><td>Désirée</td><td>mock and <i>Phytophthora infestans</i> inoculated leaves</td><td>Illumina</td><td>PE90 unstranded</td><td>~53 mio</td><td>ERR305632</td></tr><tr><td rowspan="2">Rywal</td><td rowspan="2">mock and PVY inoculated leaves</td><td rowspan="2">PacBio</td><td>Iso-Seq, 0.7–2 Kb, 2–3.5 Kb, &gt;3.5 Kb</td><td>~1.4 mio</td><td>SRR8281993, ..., SRR8282008</td></tr><tr><td>CCS</td><td></td><td></td></tr><tr><td>Rywal</td><td>mock and PVY inoculated leaves</td><td>Illumina</td><td>PE100 strand-specific</td><td>~710 mio</td><td>SRX6801457, ..., SRX6801468</td></tr><tr><td rowspan="2">PW363</td><td rowspan="2">PVY inoculated leaves</td><td rowspan="2">Illumina</td><td>DSN-normalized</td><td>~104 mio</td><td>SRR10070123, SRR10070124</td></tr><tr><td>PE90 unstranded</td><td></td><td></td></tr><tr><td>PW363</td><td>mock and PVY inoculated leaves</td><td>SOLiD</td><td>SE50 unstranded</td><td>~180 mio</td><td>SRR10065430, ..., SRR10065433</td></tr></tbody></table>

## Data records

Transcriptomic sequences of three potato genotypes, cv. Désirée, cv. Rywal and breeding clone PW363, were obtained from in-house RNA-Seq projects<sup>37–41</sup> and supplemented by publicly available cv. Désirée datasets retrieved from SRA<sup>42–44</sup> (Table 2). All generated files have been deposited to FAIRDOMhub<sup>26</sup> under project name _p_stRT_ (https://fairdomhub.org/projects/161).

The largest quantity of reads, cca. 739 mio reads of various lengths, was obtained for cv. Désirée, using Illumina and SOLiD short-read sequencing platforms. For cv. Rywal and breeding clone PW363 only mature leaf samples were available. For cv. Désirée leaf samples were augmented with samples from stems, seedlings and roots. For cv. Rywal short-read sequencing was complemented with full-length PacBio Iso-Seq sequencing of independent samples. Detailed sample information is provided in Auxiliary Table 1<sup>45</sup>.

Cv. Rywal NPR1-1 coding sequences, sequenced by the Sanger method, were deposited at NCBI GenBank under accession numbers MT210578, ..., MT210585<sup>46–53</sup>.

The GTF file with merged ITAG and PGSC gene models for *S. tuberosum* Group Phureja DM genome v4.04<sup>14</sup> is also available at FAIRDOMHub project page, as well as the cultivar-specific and pan-transcriptome assembly FASTA and annotation files<sup>54</sup>.

## Technical Validation

As a measure of assembly accuracy, the percentage of correctly assembled bases was obtained by mapping Illumina reads back to cultivar-specific initial transcripts using STAR v2.6.1d RNA-seq aligner with default parameters (Table 3). To assess the quality of transcriptomes via size-based and reference-based metrics, we ran TransRate v 1.0.1<sup>55</sup> on cultivar-specific transcriptomes (prior to and post filtering, Table 4), cultivar-specific representative transcript sequences and PGSC defined representative transcripts (Table 5). Comparative metrics for cultivar-specific coding sequences (CDS) were obtained using Conditional Reciprocal Best BLAST (CRBB)<sup>56</sup> against merged Phureja DM gene model coding sequences.

To estimate the measure of completeness and define the duplicated fraction of assembled transcriptomes (prior and post filtering cultivar-specific transcriptomes and pan-transcriptome), BUSCO v3<sup>57</sup> scores were calculated using embryophyta_odb9<sup>58</sup> lineage data (Table 6). At the cultivar-specific transcriptome level, the most diverse dataset in terms of tissues and experimental conditions, the cv. Désirée dataset, resulted in the highest BUSCO score as expected. The success in classification of representative and alternative transcripts is evident from the pan-transcriptome BUSCO scores (i.e. differences in single-copy and duplicated BUSCOs for representative and alternative dataset). The highest number of fragmented BUSCOs is observed for the breeding clone PW363, which we can probably attribute to the highest number of short-contig assemblies. Furthermore, the long-read assembly presumably contributed to the shift in favour of single-copy BUSCOs (Table 6) and uniquely mapped reads (Table 3) for cv. Rywal.

To inspect the quality of paralogue cluster assignments, multiple sequence alignments (MSA) using MAFFT v7.271<sup>59</sup> were conducted on representative and alternative sequences from paralogue clusters<sup>35</sup> containing sequences from each of the four genotypes (Désirée, PW363, Rywal and DM). Alignments were visualized using MView v1.66<sup>60</sup> (Auxiliary File 2<sup>61</sup>). These alignments were used to quality check *de novo* transcriptome assemblies
<table><caption>**Table 3.** Transcriptome quality control by RNA-seq reads remapping. Illumina paired-end reads used for generating assemblies were mapped back to the corresponding cultivar specific transcriptomes using STAR.
\*Number of reads crossing supposed splice sites. 'Initially constructed transcriptomes (prior to filtering steps).
<sup>a</sup>Relevant % of mapped reads: % of uniquely mapped reads + % of reads mapped to multiple loci.</caption><thead><tr><th>Mapping statistics/genotype</th><th>Désirée'</th><th>PW363'</th><th>Rywal'</th></tr></thead><tbody><tr><td>Number of input reads</td><td>177,149,132</td><td>52,171,015</td><td>342,767,035</td></tr><tr><td>Average input read length</td><td>178</td><td>179</td><td>199</td></tr><tr><td colspan="4"><strong>UNIQUE READS:</strong></td></tr><tr><td>Uniquely mapped reads number</td><td>64,507,790</td><td>18,416,487</td><td>206,003,021</td></tr><tr><td><strong>Uniquely mapped reads<sup>a</sup> %</strong></td><td><strong>36%</strong></td><td><strong>35%</strong></td><td><strong>60%</strong></td></tr><tr><td>Average mapped length</td><td>175</td><td>176</td><td>196</td></tr><tr><td>Number of splices*: Total</td><td>496,170</td><td>267,268</td><td>1,700,235</td></tr><tr><td>Number of splices*: Annotated (sjdb)</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Number of splices*: GT/AG</td><td>258,208</td><td>105,551</td><td>1,162,885</td></tr><tr><td>Number of splices*: GC/AG</td><td>10,749</td><td>5,693</td><td>79,495</td></tr><tr><td>Number of splices*: AT/AC</td><td>1,486</td><td>2,192</td><td>1,840</td></tr><tr><td>Number of splices*: Non-canonical</td><td>225,727</td><td>153,832</td><td>456,015</td></tr><tr><td>Mismatch rate per base %</td><td>0.50%</td><td>0.53%</td><td>0.59%</td></tr><tr><td>Deletion rate per base</td><td>0.03%</td><td>0.03%</td><td>0.03%</td></tr><tr><td>Deletion average length</td><td>2.72</td><td>2.53</td><td>3.02</td></tr><tr><td>Insertion rate per base</td><td>0.02%</td><td>0.02%</td><td>0.03%</td></tr><tr><td>Insertion average length</td><td>1.93</td><td>1.86</td><td>1.91</td></tr><tr><td colspan="4"><strong>MULTI-MAPPING READS:</strong></td></tr><tr><td>Number of reads mapped to multiple loci</td><td>98,694,222</td><td>29,366,122</td><td>108,669,657</td></tr><tr><td><strong>% of reads mapped to multiple loci<sup>a</sup></strong></td><td><strong>56%</strong></td><td><strong>56%</strong></td><td><strong>32%</strong></td></tr><tr><td>Number of reads mapped to too many loci</td><td>4,652,918</td><td>1,555,704</td><td>1,541,238</td></tr><tr><td>% of reads mapped to too many loci</td><td>2.63%</td><td>2.98%</td><td>0.45%</td></tr><tr><td colspan="4"><strong>UNMAPPED READS:</strong></td></tr><tr><td>% of reads unmapped: too many mismatches</td><td>0%</td><td>0%</td><td>0%</td></tr><tr><td>% of reads unmapped: too short</td><td>5.25%</td><td>5.43%</td><td>7.75%</td></tr><tr><td>% of reads unmapped: other</td><td>0%</td><td>0%</td><td>0%</td></tr><tr><td colspan="4"><strong>CHIMERIC READS:</strong></td></tr><tr><td>Number of chimeric reads</td><td>0</td><td>0</td><td>0</td></tr><tr><td>% of chimeric reads</td><td>0%</td><td>0%</td><td>0%</td></tr></tbody></table>

and helped us optimise the pipeline. The alignments within groups showed differences that can be attributed to biological diversity, e.g. SNPs and INDELS as well as alternative splicing<sup>61</sup>.

To estimate the proportion of transcripts originating from organisms other than potato, we performed a taxonomic classification of all cultivar-specific transcriptome sequences using Centrifuge v1.0.4<sup>62</sup> and the NCBI nt database. We used pavian<sup>63</sup> to generate classification summary reports (Auxiliary File 3<sup>64</sup>). The transcriptomes include altogether less than 1% bacterial, viral, fungal and protozoan transcripts. Other non-plant sequences are from common plant pests such the arachnid *Tetranychus urticae* in PW363 transcriptome (cca. 4%) and insects *Leptinotarsa decemlineata* (cca. 0.1%) and *Myzus persicae* (cca. 0.1%) in the Désirée transcriptome<sup>64</sup>.

**Sanger sequencing confirmation of assembled transcripts.** To further validate the quality of the constructed Rywal reference transcriptome, the cultivar-specific assembled transcript coding for NPR1-1 protein in cv. Rywal was compared to sequences amplified from isolated cDNA and sequenced by Sanger method. Total RNA was isolated from 4-week-old cv. Rywal plants using RNeasy Plant Mini kit (Qiagen). Residual genomic DNA was digested with DNase I (Qiagen). Treated total RNA was reversed transcribed using High Capacity cDNA Reverse Transcription Kit (Applied Biosystems). The full-length CDS of NPR1-1 (Sotub07g016890.1.1, ITAG genome annotation) was amplified from the cv. Rywal cDNA with forward (5'-ATGGAGAGCGGCCACGAGA-3') and reverse (5'-CTACTTTTTTCTAAACTTGTGACTGACATT-3') primers. The PCR product was inserted into pJET1.2/blunt vector with CloneJET PCR Cloning Kit (Thermo Scientific) and introduced into OneShot TOP10 Chemically Competent *Escherichia coli* cells. The plasmids were isolated from 8 transformed colonies, grown on selection media, using the GenElute Plasmid Miniprep Kit (Sigma-Aldrich). Inserts were sequenced (GATC Services, Eurofins Genomics) using pJET1.2 sequencing primers (Thermo Scientific) as well as forward (5'-CTCCAAGGTTGTGAACGAGGTACTT-3') and reverse (5'-AAGTACCTCGTTCACAACCTTGGAG-3') insert-specific primers, designed to ensure full sequence coverage in both directions.

A multiple sequence alignment comparing Sanger sequences with NPR1-1 coding sequence from the assembled cv. Rywal transcriptome (paralogue cluster stCuSTr-R_29366) and Phureja DM gene model (Sotub07g016890.1.1) was constructed using Geneious Prime 2020.1.1 (https://www.geneious.com)<sup>65</sup>. The Sanger
<table><caption>**Table 4.** Prior and post-filtering transcriptome summary statistics for potato cultivar-specific coding sequences generated by TransRate. The largest contig size at which at least 90% or 50% of bases are contained in contigs at least this length. *Reference-based summary statistics (merged Phureja DM coding sequences were used as reference). *Proportion of reference proteins with at least N% of their bases covered by a Conditional Reciprocal Best Blast (CRBB) hit.</caption><thead><tr><th rowspan="2">TransRate metrics</th><th colspan="2">Désirée</th><th colspan="2">PW363</th><th colspan="2">Rywal</th></tr><tr><th>Pre-filter (initial)</th><th>Post-filter</th><th>Pre-filter (initial)</th><th>Post-filter</th><th>Pre-filter (initial)</th><th>Post-filter</th></tr></thead><tbody><tr><td colspan="7"><strong>CONTIG METRICS</strong></td></tr><tr><td>No. sequences</td><td>350,271</td><td>197,839</td><td>273,216</td><td>159,278</td><td>134,755</td><td>79,095</td></tr><tr><td>Sequence mean length</td><td>504</td><td>792</td><td>516</td><td>775</td><td>459</td><td>707</td></tr><tr><td>No. sequences under 200 nt</td><td>125,465</td><td>25,330</td><td>88,230</td><td>17,370</td><td>52,653</td><td>13,198</td></tr><tr><td>No. sequences over 1000 nt</td><td>57,679</td><td>55,837</td><td>44,508</td><td>42,571</td><td>19,175</td><td>18,748</td></tr><tr><td>No. sequences over 10000 nt</td><td>23</td><td>23</td><td>3</td><td>3</td><td>1</td><td>1</td></tr><tr><td>'n90</td><td>369</td><td>444</td><td>366</td><td>429</td><td>351</td><td>390</td></tr><tr><td>'n50</td><td>1,194</td><td>1,209</td><td>1,110</td><td>1,131</td><td>1,227</td><td>1,218</td></tr><tr><td>GC %</td><td>41%</td><td>42%</td><td>42%</td><td>42%</td><td>42%</td><td>42%</td></tr><tr><td>Ambiguous nucleotide (N) %</td><td>0%</td><td>0%</td><td>0%</td><td>0%</td><td>0%</td><td>0%</td></tr><tr><td colspan="7"><strong>COMPARATIVE METRICS</strong></td></tr><tr><td>No. seq. with CRBB hits*</td><td>160,295</td><td>138,131</td><td>138,443</td><td>116,834</td><td>66,258</td><td>55,239</td></tr><tr><td>No. reference seq. with CRBB hits*</td><td>29,858</td><td>27,642</td><td>25,739</td><td>23,839</td><td>23,549</td><td>22,163</td></tr><tr><td>coverage50*%</td><td>25,991</td><td>24,586</td><td>21,875</td><td>20,620</td><td>20,258</td><td>19,538</td></tr><tr><td>coverage95*%</td><td>19,329</td><td>18,246</td><td>15,664</td><td>14,727</td><td>14,967</td><td>14,470</td></tr><tr><td>Reference coverage*</td><td>65%</td><td>63%</td><td>56%</td><td>54%</td><td>53%</td><td>52%</td></tr></tbody></table>

<table><caption>**Table 5.** Summary statistics for potato cultivar-specific representative transcript sequences generated by TransRate. The largest contig size at which at least 90% or 50% of bases are contained in contigs at least this length. *Reference-based summary statistics (merged Phureja DM coding sequences were used as reference). *Proportion of reference proteins with at least N% of their bases covered by a Conditional Reciprocal Best Blast (CRBB) hit. PGSC\_DM\_v3.4\_transcript-update\_representative\_basta.zip file from Spud DB was used for Phureja-specific representative transcript sequences (PGSC).</caption><thead><tr><th>TransRate metrics</th><th>Désirée</th><th>PW363</th><th>Rywal</th><th>PGSC<sup>+</sup></th></tr></thead><tbody><tr><td colspan="5"><strong>CONTIG METRICS</strong></td></tr><tr><td>No. sequences</td><td>57,943</td><td>43,883</td><td>36,336</td><td>39,031</td></tr><tr><td>Sequence mean length</td><td>922</td><td>926</td><td>1,028</td><td>1,283</td></tr><tr><td>No. sequences under 200 nt</td><td>875</td><td>1,377</td><td>1,310</td><td>87</td></tr><tr><td>No. sequences over 1000 nt</td><td>18,500</td><td>14,545</td><td>14,307</td><td>20,226</td></tr><tr><td>No. sequences over 10000 nt</td><td>13</td><td>6</td><td>2</td><td>0</td></tr><tr><td>'n90</td><td>369</td><td>387</td><td>440</td><td>645</td></tr><tr><td>'n50</td><td>1,566</td><td>1,535</td><td>1,673</td><td>1,726</td></tr><tr><td>GC %</td><td>40%</td><td>41%</td><td>41%</td><td>40%</td></tr><tr><td>Ambiguous nucleotide (N) %</td><td>0%</td><td>0%</td><td>0%</td><td>0%</td></tr><tr><td colspan="5"><strong>COMPARATIVE METRICS</strong></td></tr><tr><td>No. seq. with CRBB hits*</td><td>38,034</td><td>30,826</td><td>28,389</td><td>38,600</td></tr><tr><td>No. reference seq. with CRBB hits*</td><td>25,094</td><td>21,751</td><td>21,299</td><td>37,534</td></tr><tr><td>coverage50*%</td><td>12,799</td><td>10,693</td><td>7,909</td><td>36,379</td></tr><tr><td>coverage95*%</td><td>8,053</td><td>6,430</td><td>5,053</td><td>30,187</td></tr><tr><td>Reference coverage*</td><td>33%</td><td>28%</td><td>20%</td><td>75%</td></tr></tbody></table>

sequencing revealed the presence of several distinct gene variants in the analysed colonies, differing in the SNP pattern. The PBdnRY1_9437 sequence is validated by complete sequence identity with the two colonies (seq. 7 and 8), while the CLCdnRY_10265 shares two SNPs with seq. 4–6, matching the SNP pattern partially. Although the Phureja DM gene model also shares two SNPs with some of the colonies, its overall SNP pattern differs significantly from cv. Rywal, distinguishing the cultivar-specific transcripts from that of Phureja DM (Fig. 4).

## Usage Notes

**Insights into variability of potato transcriptomes.** Based on the comparison of cultivar-specific transcriptomes we identified cca. 23,000, 13,000, and 7,500 paralogue groups of transcripts in cv. Désirée, breeding clone PW363 and cv. Rywal, respectively, not present in the merged Phureja DM gene model. The addition of
<table><caption>**Table 6.** Assessment of completeness of constructed transcriptomes. Percentage of BUSCOs identified in each transcriptome assembly step. (S): Complete and single-copy BUSCOs %; (D): Complete and duplicated BUSCOs % (C): Complete BUSCOs (S + D) % (F): Fragmented BUSCOs % (M): Missing BUSCOs % representative alt: alternative *Database size: 1440.</caption><thead><tr><th>cv. Désirée</th><th>initial rep+alt</th><th>post 1st filtering rep+alt</th><th>final rep+alt</th></tr></thead><tbody><tr><td>(S)</td><td>37.8</td><td>37.8</td><td>37.4</td></tr><tr><td>(D)</td><td>59.4</td><td>59.2</td><td>58.4</td></tr><tr><td>(C)</td><td><strong>97.2</strong></td><td><strong>97.0</strong></td><td><strong>95.8</strong></td></tr><tr><td>(F)</td><td>1.1</td><td>1.2</td><td>1.4</td></tr><tr><td>(M)</td><td>1.7</td><td>1.8</td><td>2.8</td></tr><tr><th>breeding clone PW363</th><th>initial rep+alt</th><th>post 1st filtering rep+alt</th><th>final rep+alt</th></tr><tr><td>(S)</td><td>39.9</td><td>39.2</td><td>38.4</td></tr><tr><td>(D)</td><td>51.7</td><td>51.2</td><td>50.9</td></tr><tr><td>(C)</td><td><strong>91.6</strong></td><td><strong>90.4</strong></td><td><strong>89.3</strong></td></tr><tr><td>(F)</td><td>2.9</td><td>3.4</td><td>3.5</td></tr><tr><td>(M)</td><td>5.6</td><td>6.2</td><td>7.2</td></tr><tr><th>cv. Rywal</th><th>initial rep+alt</th><th>post 1st filtering rep+alt</th><th>final rep+alt</th></tr><tr><td>(S)</td><td>55.8</td><td>55.8</td><td>55.1</td></tr><tr><td>(D)</td><td>35.2</td><td>34.8</td><td>34.7</td></tr><tr><td>(C)</td><td><strong>91.0</strong></td><td><strong>90.6</strong></td><td><strong>89.8</strong></td></tr><tr><td>(F)</td><td>2.4</td><td>2.6</td><td>2.7</td></tr><tr><td>(M)</td><td>6.5</td><td>6.9</td><td>7.5</td></tr><tr><th>pan-transcriptome</th><th>rep</th><th>alt</th><th>rep+alt</th></tr><tr><td>(S)</td><td>92.2</td><td>11.0</td><td>3.9</td></tr><tr><td>(D)</td><td>6.1</td><td>85.9</td><td>95.6</td></tr><tr><td>(C)</td><td><strong>98.3</strong></td><td><strong>96.9</strong></td><td><strong>99.4</strong></td></tr><tr><td>(F)</td><td>1.4</td><td>1.3</td><td>0.3</td></tr><tr><td>(M)</td><td>0.3</td><td>1.7</td><td>0.3</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Fig. 4** Sanger sequencing validates the constructed cultivar specific transcriptome. Multiple sequence alignment of NPR1-1 coding sequence obtained from eight *E. coli* colonies (NPR1-1 seq. 1–8) by the Sanger method, assembled short or long-read cv. Rywal transcripts and Phureja DM gene model (Sotub07g016890.1.1). Grey - sequence identity, black - SNPs. The alignment was prepared and visualised with Geneious Prime 2020.1.1<sup>65</sup>.</figcaption></figure>

Iso-Seq dataset in the case of cv. Rywal confirms that long reads contribute to less fragmentation of the *de novo* transcriptome. It is therefore recommended to generate at least a subset of data with one of the long-read technologies to complement the short-read RNA-seq. As it can be seen by the reduction rate in PW363 (24-fold), producing additional short-read assemblies does not contribute to the transcriptome quality as much as having several tissues or a combination of 2<sup>nd</sup> and 3<sup>rd</sup> generation sequencing (38-fold Rywal).

From all four genotypes, cv. Désirée has the highest number of cultivar-specific representative transcripts, which can be attributed to having the most diverse input dataset used for the *de novo* assemblies in terms of tissues sequenced (stem, seedlings and roots) and experimental conditions covered. Cv. Désirée also benefited from the inclusion of a DSN Illumina library to capture low level expressed transcripts. However, even the leaf-specific reference transcriptomes of cv. Rywal and breeding clone PW363 include thousands of specific genes, indicating that cultivar specific gene content is common. Remarkably, we identified several interesting features when inspecting paralogue groups of transcripts, demonstrating the variability of sequences in potato haplotypes and the presence of the alternative splicing variants that contribute to the pan-transcriptome (Fig. 5, Auxiliary File 2<sup>61</sup>).
<figure><img src="image_7.png" /><figcaption>**Fig. 5** Transcript variants present in pan-transcriptome paralogue gene groups. a) Alignment part of stPanTr_010101 with two PW363-specific SNPs marked by red dots. Such SNPs can be used to design cultivar- or allele-specific qPCR assays. b) Alignment part of stPanTr_074336 showing an alternative splice variant in Désirée, (VdnDe4_33782). Both multiple sequence alignments were made using ClustalOmega v 1.2.1⁶⁹ and visualized with MView v 1.66⁶⁰. The remaining alignments can be found in Auxiliary file 2⁶¹.</figcaption></figure>
<table><thead><tr><th>Reference</th><th>Désirée</th><th>ITAG/PGSC</th><th>PGSC</th><th>Désirée</th><th>ITAG/PGSC</th><th>PGSC</th></tr></thead><tbody><tr><td>Mapping statistics/Sample</td><td>SRR10416847</td><td></td><td></td><td>SRR10416848</td><td></td><td></td></tr><tr><td>Number of input reads</td><td>14,953,659</td><td></td><td></td><td>14,610,172</td><td></td><td></td></tr><tr><td>Average input read length</td><td>252</td><td></td><td></td><td>252</td><td></td><td></td></tr><tr><td colspan="7">UNIQUE READS:</td></tr><tr><td>Uniquely mapped reads<sup>#</sup> %</td><td>74%</td><td>69%</td><td>71%</td><td>73%</td><td>67%</td><td>69%</td></tr><tr><td>Average mapped length</td><td>246</td><td>244</td><td>246</td><td>246</td><td>244</td><td>246</td></tr><tr><td>Number of splices: Total</td><td>123,490</td><td>104,658</td><td>159,842</td><td>127,414</td><td>110,423</td><td>165,595</td></tr><tr><td>Number of splices: Non-canonical</td><td>40,960</td><td>43,764</td><td>45,750</td><td>40,898</td><td>45,459</td><td>42,978</td></tr><tr><td>Mismatch rate per base %</td><td>0.75%</td><td>0.94%</td><td>1.00%</td><td>0.74%</td><td>0.93%</td><td>1.00%</td></tr><tr><td>Deletion rate per base</td><td>0.05%</td><td>0.04%</td><td>0.06%</td><td>0.05%</td><td>0.04%</td><td>0.06%</td></tr><tr><td>Deletion average length</td><td>3.42</td><td>3.46</td><td>2.88</td><td>3.41</td><td>3.47</td><td>2.89</td></tr><tr><td>Insertion rate per base</td><td>0.03%</td><td>0.02%</td><td>0.04%</td><td>0.03%</td><td>0.02%</td><td>0.04%</td></tr><tr><td>Insertion average length</td><td>2.21</td><td>2.99</td><td>2.52</td><td>2.24</td><td>3.00</td><td>2.53</td></tr><tr><td colspan="7">MULTI-MAPPING READS:</td></tr><tr><td>% of reads mapped to multiple loci<sup>#</sup></td><td>4.94%</td><td>1.95%</td><td>2.74%</td><td>4.86%</td><td>1.68%</td><td>2.52%</td></tr><tr><td>% of reads mapped to too many loci</td><td>0%</td><td>0%</td><td>0%</td><td>0%</td><td>0%</td><td>0%</td></tr><tr><td colspan="7">UNMAPPED READS:</td></tr><tr><td>% of reads unmapped: too many mismatches</td><td>0%</td><td>0%</td><td>0%</td><td>0%</td><td>0%</td><td>0%</td></tr><tr><td>% of reads unmapped: too short</td><td>21%</td><td>29%</td><td>27%</td><td>22%</td><td>31%</td><td>29%</td></tr><tr><td>% of reads unmapped: other</td><td>0%</td><td>0%</td><td>0%</td><td>0%</td><td>0%</td><td>0%</td></tr><tr><td>Mapping statistics/Sample</td><td>SRR10416849</td><td></td><td></td><td>all samples</td><td></td><td></td></tr><tr><td>Number of input reads</td><td>14,755,430</td><td></td><td></td><td>44,319,261</td><td></td><td></td></tr><tr><td>Average input read length</td><td>252</td><td></td><td></td><td>252</td><td></td><td></td></tr><tr><td colspan="7">UNIQUE READS:</td></tr><tr><td>Uniquely mapped reads<sup>#</sup> %</td><td>49%</td><td>44%</td><td>46%</td><td>66%</td><td>60%</td><td>62%</td></tr><tr><td>Average mapped length</td><td>245</td><td>243</td><td>246</td><td>246</td><td>244</td><td>246</td></tr><tr><td>Number of splices: Total</td><td>95,409</td><td>77,103</td><td>115,083</td><td>346,313</td><td>292,184</td><td>440,520</td></tr><tr><td>Number of splices: Non-canonical</td><td>33,269</td><td>31,224</td><td>31,065</td><td>115,127</td><td>120,447</td><td>119,793</td></tr><tr><td>Mismatch rate per base %</td><td>0.75%</td><td>0.94%</td><td>1.01%</td><td>0.75%</td><td>0.94%</td><td>1.00%</td></tr><tr><td>Deletion rate per base</td><td>0.06%</td><td>0.05%</td><td>0.07%</td><td>0.05%</td><td>0.05%</td><td>0.06%</td></tr><tr><td>Deletion average length</td><td>3.23</td><td>3.25</td><td>2.83</td><td>3.36</td><td>3.41</td><td>2.87</td></tr><tr><td>Insertion rate per base</td><td>0.03%</td><td>0.03%</td><td>0.04%</td><td>0.03%</td><td>0.02%</td><td>0.04%</td></tr><tr><td>Insertion average length</td><td>2.27</td><td>3.15</td><td>2.56</td><td>2.24</td><td>3.04</td><td>2.54</td></tr><tr><td colspan="7">MULTI-MAPPING READS:</td></tr><tr><td>% of reads mapped to multiple loci<sup>#</sup></td><td>3.60%</td><td>1.19%</td><td>1.94%</td><td>4.47%</td><td>1.61%</td><td>2.40%</td></tr><tr><td>% of reads mapped to too many loci</td><td>0%</td><td>0%</td><td>0%</td><td>0%</td><td>0%</td><td>0%</td></tr><tr><td colspan="7">UNMAPPED READS:</td></tr><tr><td>% of reads unmapped: too many mismatches</td><td>0%</td><td>0%</td><td>0%</td><td>0%</td><td>0%</td><td>0%</td></tr><tr><td>% of reads unmapped: too short</td><td>47%</td><td>55%</td><td>52%</td><td>30%</td><td>38%</td><td>36%</td></tr><tr><td>% of reads unmapped: other</td><td>0%</td><td>0%</td><td>0%</td><td>0%</td><td>0%</td><td>0%</td></tr></tbody></table>

**Table 7.** Mapping of independent dataset to newly assembled cultivar specific reference transcriptome. Mapping statistics for Désirée leaf samples under drought stress to Désirée, ITAG/PGSC merged and PGSC representative transcriptome sequences is shown. RNA-seq data from Désirée leaf samples under drought stress retrieved from the GEO Series GSE140083 – “Transcriptome profiles of contrasting potato (Solanum tuberosum L.) genotypes under water stress”. No chimeric reads detected. <sup>#</sup>Relevant % of mapped reads: % of uniquely mapped reads + % of reads mapped to multiple loci.

It should be noted, that the reconstructed transcriptomes include also the meta-transcriptome stemming from microbial communities present in sampled potato tissues (Auxiliary File 3<sup>64</sup>). We decided not to apply any filter on these transcripts. Inclusion of meta-transcripts makes it possible to also investigate the diversity of plant-associated endo- and epiphytes. The majority of these microbial transcripts will have microbial annotations, facilitating their future removal when necessary for other experiments.

**Cultivar-specific transcriptomes can improve high-throughput sequencing analyses.** Most gene expression studies have been based on either potato UniGenes, assembled from a variety of potato expressed sequence tags (e.g. StGI, POCI), or the reference DM genome transcript models. Studies based on any of these resources have provided useful information on potato gene expression, but also have major drawbacks.
When using the DM genome as a reference for mapping RNA-Seq reads, the potato research community faces the existence of two overlapping, but not identical, gene model predictions. When using either of available GFFs, we are missing some of the genes known to be encoded in the assembled scaffold. The newly generated merged DM-based GTF helps to circumvent this problem, but even in the merged GTF the cultivar-specific genes and variations are not considered. Differences in expression and important marker transcripts can therefore be missed. In addition, the computational prediction of DM transcript isoforms is incomplete and, in some cases, gene models are incorrectly predicted. On the other hand, the inherent heterogeneity and redundancy of UniGenes or similarly combined transcript sets causes the mapping of short reads to multiple transcripts, thus making the result interpretation more difficult. The cultivar-specific transcriptomes presented here are an improvement, as they include several transcripts missing in the Phureja DM transcriptome. The transcriptomes are also a valuable asset for other high-throughput sequencing applications, such as sRNA-Seq, Degradome-Seq or ATAC-Seq, as we now have more detailed information also on transcript variability within one locus.

The benefit of using cultivar-specific transcriptomes was demonstrated through mapping statistics for Désirée leaf samples under drought stress⁶⁶ to Désirée, ITAG/PGSC merged and PGSC representative transcriptome sequences. Taking all three samples together, 5.5% more reads were uniquely mapped to Désirée than the ITAG/ PGSC Phureja DM gene models (Table 7). From the reads mapped to Désirée transcriptome, 5.3% mapped to Désirée-specific transcripts and 8.2% to transcripts specific to Group Tuberosum genotypes (Auxiliary Table 8⁶⁷).

Cultivar-specific transcriptomes may also help improve mass-spectroscopy based proteomics. Comprehensive protein databases, obtained from transcriptomic data, offer a higher chance of finding significant targets with peptide spectrum match algorithms, thus enhancing the detection and sensitivity of protein abundance measurements⁶⁸.

**Using transcriptomes to inform qPCR amplicon design.** Aligning transcript coding sequences from a StPanTr paralogue cluster can be used to inform qPCR primer design in order to study gene expression in different cultivars or specific isoforms by selecting variable transcript regions (Fig. 5). On the other hand, when qPCR assays need to cover multiple cultivars, the nucleotide alignments can be inspected for conservative regions for design.

**Limitations.** Although the here presented transcriptomes are of higher quality than assemblies produced by a single assembler they are neither comprehensive nor perfectly accurate. Firstly, for cv. Rywal and breeding clone PW363 reference transcriptomes are leaf-specific. Additionally, some transcripts are represented by partial sequence(s) and some might be misassembled. Due to very low expression, some transcripts are missed. One also has to note that some true transcripts with very high similarity to the representative transcript, potentially expressed from distinct alleles differing in the SNP pattern, were removed by filtering. Furthermore, most long non-coding mRNAs are discarded in the EvidentialGene step.

When mapping reads back to the initial transcriptomes a fraction of reads was spliced (Table 3), indicating that the transcripts to which these reads map might originate from splice variants that are not present in the transcriptome. We would however like to point out that contiguity of this reads was not confirmed thus they might represent misassemblies. In Désirée 14,094 transcripts had uniquely mapped spliced reads (7% of the filtered transcriptome), and similarly 7,493 (9%) transcripts in Rywal and 10,220 (6%) in PW363. Low number of non-canonical sites was however confirmed by STARlong mapping: 118 for Désirée, 59 for Rywal and 101 for PW363.

An indication of potential problems with some of the assemblies is that the number of transcripts with length less than 200 nt is 10x higher in constructed cultivar-specific transcriptomes than in DM Phureja PGSC representative transcripts (875, 1377 and 1310 vs. 87 for Désirée, PW363 and Rywal vs. DM Phureja) while the transcript length is about 300 nt shorter in average (Table 5). To check for any misassemblies for genes of interest the users are advised to check the MSA of the paralogue cluster (examples for core transcriptome subset can be found in Auxiliary File 2⁶¹). Some of the assembled sequences are also not of potato origin. Vecscreen implied vector segment contamination in 46 Désirée, 9 PW363 and 31 Rywal sequences, while Centrifuge classified some sequences as originating from common plant pests, among which are the meta-transcriptome sequences of potato microbial communities and the complete genome of the Potato Virus Y (PVY), present due to experimental treatments.

**Code availability**

All used Bash, Perl, Python, and R/Markdown custom code and scripts complemented with intermediate and processed data (input and output files), and all other supporting information that enable reproduction and re-use are available at FAIRDOMHub under project name _p_stRT (fairdomhub.org/projects/161) under CC BY 4.0 licenses. Data were locally stored in a ISA-tab compliant project directory tree generated using pISA-tree (github.com/NIB-SI/pISA) and uploaded to FAIRDOMHub repository using FAIRDOMHub API and R package pisar (github.com/NIB-SI/pisar).

Received: 27 November 2019; Accepted: 19 June 2020;

Published online: 24 July 2020

**References**

1. Hardigan, M. A. et al. Genome diversity of tuber-bearing Solanum uncovers complex evolutionary history and targets of domestication in the cultivated potato. *Proc. Natl. Acad. Sci. USA* **114**, E9999–E10008 (2017).

2. Gan, X. et al. Multiple reference genomes and transcriptomes for Arabidopsis thaliana. *Nature* **477**, 419–423 (2011).

3. Hirsch, C. N. et al. Insights into the maize pan-genome and pan-transcriptome. *Plant Cell* **26**, 121–135 (2014).

4. Jin, M. et al. Maize pan-transcriptome provides novel insights into genome complexity and quantitative trait variation. *Sci. Rep.* **6**, 18936 (2016).
5. Zhao, Q. et al. Pan-genome analysis highlights the extent of genomic variation in cultivated and wild rice. *Nat. Genet.* **50**, 278–284 (2018).

6. Montenegro, J. D. et al. The pangenome of hexaploid bread wheat. *Plant J.* **90**, 1007–1013 (2017).

7. Li, Y. H. et al. De novo assembly of soybean wild relatives for pan-genome analysis of diversity and agronomic traits. *Nat. Biotechnol.* **32**, 1045–1052 (2014).

8. Xu, X. et al. Genome sequence and analysis of the tuber crop potato. *Nature* **475**, 189–195 (2011).

9. Liu, Y. et al. Comparative transcriptome analysis of white and purple potato to identify genes involved in anthocyanin biosynthesis. *PLoS One* **10**, e0129148 (2015).

10. Sato, S. et al. The tomato genome sequence provides insights into fleshy fruit evolution. *Nature* **485**, 635–641 (2012).

11. Hölzer, M. & Marz, M. De novo transcriptome assembly: A comprehensive cross-species comparison of short-read RNA-Seq assemblers. *GigaScience* **8**, giz039 (2019).

12. Gilbert, D. G. Genes of the pig, Sus scrofa, reconstructed with EvidentialGene. *PeerJ* **7**, e6374 (2019).

13. Hirsch, C. D. et al. Spud DB: A resource for mining sequences, genotypes, and phenotypes to accelerate potato breeding. *Plant Genome* **7**, 1–12 (2014).

14. Zagorščak, M. Merged ITAG and PGSC gene models for *S. tuberosum* Group Phureja DM genome v4.04. *FAIRDOMHub* https://doi.org/10.15490/FAIRDOMHUB.1.DATAFILE.3408.1 (2020).

15. Zerbino, D. R. Using the Velvet de novo assembler for short-read sequencing technologies. *Curr. Protoc. in Bioinformatics* **31**, 11.5.1–11.5.12 (2010).

16. Crusoe, M. R. et al. The khmer software package: enabling efficient nucleotide sequence analysis. *F1000Res.* **4**, 900 (2015).

17. Tseng, E. cdna_cupcake v9.0.1 (2019).

18. Grabherr, M. G. et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. *Nat. Biotechnol.* **29**, 644–652 (2011).

19. Schulz, M. H., Zerbino, D. R., Yingron, M. & Birney, E. Oases: Robust de novo RNA-seq assembly across the dynamic range of expression levels. *Bioinformatics* **28**, 1086–1092 (2012).

20. Bushmanova, E., Antipov, D., Lapidus, A. & Prjibelski, A. D. rnaSPAdes: a de novo transcriptome assembler and its application to RNA-Seq data. *GigaScience* **8**, giz100 (2019).

21. Zagorščak, M. & Petek, M. Detailed de novo assemblies information table listing parameters used for short-read de novo assembly generation. *FAIRDOMHub* https://doi.org/10.15490/FAIRDOMHUB.1.DATAFILE.3091.1 (2020).

22. He, B. et al. Optimal assembly strategies of transcriptome related to ploidies of eukaryotic organisms. *BMC Genomics* **16**, 65 (2015).

23. Slater, G. S. C. & Birney, E. Automated generation of heuristics for biological sequence comparison. *BMC Bioinformatics* **6**, 31 (2005).

24. Fu, L., Niu, B., Zhu, Z., Wu, S. & Li, W. CD-HIT: accelerated for clustering the next-generation sequencing data. *Bioinformatics* **28**, 3150–3152 (2012).

25. Dobin, A. et al. STAR: Ultrafast universal RNA-seq aligner. *Bioinformatics* **29**, 15–21 (2013).

26. Wolstencroft, K. et al. FAIRDOMHub: a repository and collaboration environment for sharing systems biology research. *Nucleic Acids Res.* **45**, D404–D407 (2017).

27. Jones, P. et al. InterProScan 5: Genome-scale protein function classification. *Bioinformatics* **30**, 1236–1240 (2014).

28. Buchfink, B., Xie, C. & Huson, D. H. Fast and sensitive protein alignment using DIAMOND. *Nat. Methods* **12**, 59–60 (2014).

29. Schäffer, A. A. et al. VecScreen_plus_taxonomy: imposing a tax(onomy) increase on vector contamination screening. *Bioinformatics* **34**, 755–759 (2017).

30. Zagorščak, M. Désirée biological evidence filtering results. *FAIRDOMHub* https://doi.org/10.15490/FAIRDOMHUB.1.DATAFILE.3407.1 (2020).

31. Zagorščak, M. PW363 biological evidence filtering results. *FAIRDOMHub* https://doi.org/10.15490/FAIRDOMHUB.1.DATAFILE.3406.1 (2020).

32. Zagorščak, M. Rywal biological evidence filtering results. *FAIRDOMHub* https://doi.org/10.15490/FAIRDOMHUB.1.DATAFILE.3405.1 (2020).

33. De Nooy, W., Mrvar, A. & Batagelj, V. *Exploratory Social Network Analysis With Pajek* 3rd edn. (Cambridge University Press, 2018)

34. Zagorščak, M. EvidentialGene Summary Statistics for PacBio sequences. *FAIRDOMHub* https://doi.org/10.15490/FAIRDOMHUB.1.DATAFILE.3363.1 (2020).

35. Zagorščak, M. Paralogue cluster information for cultivar-specific and pan-transcriptome sequences extended with annotations and quality classification. *FAIRDOMHub* https://doi.org/10.15490/FAIRDOMHUB.1.DATAFILE.3721.2 (2020).

36. Khan, A. W. et al. Super-Pangenome by Integrating the Wild Side of a Species for Accelerated Crop Improvement. *Trends Plant Sci.* **25**, 148–158 (2020).

37. NCBI Sequence Read Archive https://identifiers.org/ncbi/insdc.sra:SRP220411 (2019).

38. NCBI Sequence Read Archive https://identifiers.org/ncbi/insdc.sra:SRP040682 (2015).

39. NCBI Sequence Read Archive https://identifiers.org/ncbi/insdc.sra:SRP220250 (2019).

40. NCBI Sequence Read Archive https://identifiers.org/ncbi/insdc.sra:SRP220356 (2019).

41. NCBI Sequence Read Archive https://identifiers.org/ncbi/insdc.sra:SRP172523 (2019).

42. NCBI Sequence Read Archive https://identifiers.org/ncbi/insdc.sra:SRP069961 (2016).

43. NCBI Sequence Read Archive https://identifiers.org/ncbi/insdc.sra:SRP083083 (2017).

44. European Nucleotide Archive https://identifiers.org/ena.embl:ERP003480 (2014).

45. Zagorščak, M. & Petek, M. Detailed sample information table used to generate the de novo transcriptome assemblies. *FAIRDOMHub* https://doi.org/10.15490/FAIRDOMHUB.1.DATAFILE.3090.1 (2020).

46. Tomaž, Š., Gruden, K. & Coll, A. *GenBank* https://identifiers.org/ncbi/insdc:MT210578 (2020).

47. Tomaž, Š., Gruden, K. & Coll, A. *GenBank* https://identifiers.org/ncbi/insdc:MT210579 (2020).

48. Tomaž, Š., Gruden, K. & Coll, A. *GenBank* https://identifiers.org/ncbi/insdc:MT210580 (2020).

49. Tomaž, Š., Gruden, K. & Coll, A. *GenBank* https://identifiers.org/ncbi/insdc:MT210581 (2020).

50. Tomaž, Š., Gruden, K. & Coll, A. *GenBank* https://identifiers.org/ncbi/insdc:MT210582 (2020).

51. Tomaž, Š., Gruden, K. & Coll, A. *GenBank* https://identifiers.org/ncbi/insdc:MT210583 (2020).

52. Tomaž, Š., Gruden, K. & Coll, A. *GenBank* https://identifiers.org/ncbi/insdc:MT210584 (2020).

53. Tomaž, Š., Gruden, K. & Coll, A. *GenBank* https://identifiers.org/ncbi/insdc:MT210585 (2020).

54. Blejec, A. & Zagorščak, M. Annotated fasta files. *FAIRDOMHub* https://doi.org/10.15490/FAIRDOMHUB.1.ASSAY.1268.2 (2020).

55. Smith-Unna, R., Boursnell, C., Patro, R., Hibberd, J. M. & Kelly, S. TransRate: Reference-free quality assessment of de novo transcriptome assemblies. *Genome Res.* **26**, 1134–1144 (2016).

56. Aubry, S., Kelly, S., Kümpers, B. M. C., Smith-Unna, R. D. & Hibberd, J. M. Deep evolutionary comparison of gene expression identifies parallel recruitment of trans-factors in two independent origins of C4 photosynthesis. *PLoS Genet.* **10**, e1004365 (2014).

57. Simão, F. A., Waterhouse, R. M., Ioannidis, P., Kriventseva, E. V. & Zdobnov, E. M. BUSCO: Assessing genome assembly and annotation completeness with single-copy orthologs. *Bioinformatics* **31**, 3210–3212 (2015).

58. Waterhouse, R. M. et al. BUSCO Applications from quality assessments to gene prediction and phylogenomics. *Mol. Biol. Evol.* **35**, 543–548 (2017).
59. Nakamura, T., Yamada, K. D., Tomii, K. & Katoh, K. Parallelization of MAFFT for large-scale multiple sequence alignments. *Bioinformatics* **34**, 2490–2492 (2018).

60. Brown, N. P., Leroy, C. & Sander, C. MView: A web-compatible database search or multiple alignment viewer. *Bioinformatics* **14**, 380–381 (1998).

61. Zagorščak, M. Paralogue clusters selection multiple sequence alignments using ClustalOmega v1.2.1 or MAFFT v7.271 and MView v1.66. *FAIRDOMHub* https://doi.org/10.15490/FAIRDOMHUB.1.DATAFILE.3691.1 (2020).

62. Kim, D., Song, L., Breitwieser, F. P. & Salzberg, S. L. Centrifuge: rapid and sensitive classification of metagenomic sequences. *Genome Res.* **26**, 1721–1729 (2016).

63. Breitwieser, F. P. & Salzberg, S. L. Pavian: interactive analysis of metagenomics data for microbiome studies and pathogen identification. *Bioinformatics* **36**, 1303–1304 (2019).

64. Zagorščak, M. Evaluation of constructed reference transcriptomes for presence of the metagenome using Pavian visualisation of the Centrifuge taxonomic classification program output. *FAIRDOMHub* https://doi.org/10.15490/FAIRDOMHUB.1.DATAFILE.3509.1 (2020).

65. Kearse, M. et al. Geneious Basic: An integrated and extendable desktop software platform for the organization and analysis of sequence data. *Bioinformatics* **28**, 1647–1649 (2012).

66. NCBI Sequence Read Archive (SRA), https://identifiers.org/ncbi/insdc.sra:SRP229087 (2019).

67. Zagorščak, M. Read count summary for Désirée drought samples mapped to the representative Phureja DM and Désirée reference transcriptomes. *FAIRDOMHub* https://doi.org/10.15490/FAIRDOMHUB.1.DATAFILE.3722.2 (2020).

68. Luge, T., Fischer, C. & Sauer, S. Efficient application of de novo RNA assemblers for proteomics informed by transcriptomics. *J. Proteome Res.* **15**, 3938–3943 (2016).

69. Sievers, F. & Higgins, D. G. Clustal Omega for making accurate alignments of many protein sequences. *Protein Sci.* **27**, 135–145 (2018).

## Acknowledgements

We thank Robin Buell for potato PGSC-ITAG gene model reference table, Thomas Doak and Don Gilbert for advice on using EvidentialGene scripts, Henrik Krnec for BLAST output parser, Andrej Blejec for FAIRDOMhub API usage in R, Špela Baebler for providing DSN Illumina-sequenced samples and Karmen Pogačar for providing Rywal cDNA samples for technical validation. This project was supported by the Slovenian Research Agency (grants P4-0165, J4-4165, J4-7636, J4-8228 and J4-9302), COST actions CA15110 (CHARME) and CA15109 (COSTNET).

## Author contributions

M.P., K.G., Ž.R. and M. Zagorščak participated in study design and evaluation of transcriptomes. A.C. provided Rywal Illumina-sequenced samples. M.P. collected and pre-processed RNA-Seq datasets and produced CLC, Velvet/Oases and rnaSPAdes de novo assemblies. M. Zouine produced Trinity assemblies. E.T. processed Iso-Seq data. M.P. and M. Zagorščak run tr2aacds scripts and transcriptome annotation, filtering and evaluation software. Ž.R. merged PGSC and ITAG gene models of reference potato genome. M.P., Ž.R. and M. Zagorščak generated the pan-transcriptome. Š.T. conducted bacterial transformation and Sanger sequencing analysis. K.G. secured funding, and managed the project. M.P., M. Zagorščak, Ž.R., Š.T. and K.G. wrote and edited the manuscript. S.S. helped with interpretation of EvidentialGene results, provided advice on filtering of transcriptomes and language editing. All authors have read and commented the manuscript and approved the final submission. Conceptualization: K.G., M.P., Ž.R., M. Zagorščak; Supervision: K.G.; Project Administration: K.G.; Investigation: M.P., M. Zagorščak; Formal Analysis: M.P., Ž.R., M. Zagorščak; Software: M.P., Ž.R., E.T., M. Zagorščak; Methodology: K.G., M.P., Ž.R., M. Zagorščak; Validation: M.P., Ž.R., S.S., Š.T., M. Zagorščak; Data Curation: M.P., M. Zagorščak; Resources: A.C., K.G., M.P., E.T., Š.T., M. Zouine; Funding Acquisition: K.G.; Writing - Original Draft Preparation: M.P., M. Zagorščak; Writing - Review & Editing: all authors; Visualization: M.P., M. Zagorščak.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** is available for this paper at https://doi.org/10.1038/s41597-020-00581-4.

**Correspondence** and requests for materials should be addressed to M.P. or M.Z.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_8.png" /></figure>

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

The Creative Commons Public Domain Dedication waiver http://creativecommons.org/publicdomain/zero/1.0/ applies to the metadata files associated with this article.

© The Author(s) 2020