Thousands of potato cultivars exists and most of them are tetraploid (2n = 4x = 48; FAO, 2008). However, chromosome level genome assemblies for tetraploid potato clones became only recently available (Hoopes et al., 2022), after the analyses for this study were finalized. Instead, we used four tuber-bearing diploid potato clones belonging to cultivated, non-cultivated, and wild potato species for which chromosome-level genome assemblies are available. The cultivated potato Solanum tuberosum ssp. tuberosum L. is represented in our study by a diploid clone (hereafter referred to as dAg) which was derived from the tetraploid elite potato cultivar Agria (tAg; Freire et al., 2021). The non-cultivated potato clones include Solanum tuberosum L. DM1-3516 R44 (hereafter referred to as DM), which is a doubled monoploid clone derived from the group Phureja (Pham et al., 2020) and Solanum tuberosum L. RH89-039-16 (hereafter referred to as RH) which is a diploid clone derived from a cross between a dihaploid and a diploid potato (Zhou et al., 2020). The wild clone in our study is Solanum chacoense M6 (hereafter referred to as M6; Leisner et al., 2018). Both the sequence (genome and gene) and annotation files for DM (version 6.1), RH, and M6 were downloaded from http://potato.plantbiology.msu.edu and for dAg was obtained from Freire et al. (2021). In addition, we obtained transposable elements (TEs) annotation for dAg and DM from the above-mentioned sources. Due to the lack or absence of chromosome-level TE annotation for M6 and RH, we excluded TE annotation for these two genotypes from the analyses.

In this study, we improved the existing gene annotation for dAg using the PASA pipeline v2.5.0 (Haas et al., 2003) followed by classifying the resulting annotation into full-length and partial gene models using AGAT v0.8.0 (Dainat et al., 2022). In that procedure, we used transcriptome datasets generated as part of dAg genome sequencing (Freire et al., 2021) to improve the existing gene annotations.

For reasons of consistency, we have performed functional annotation for the longest isoform of high-confidence genes of all four potato genomes using the AHRD pipeline.¹ Orthologs among the four potato genomes were predicted by feeding protein sequences of the longest isoform of high-confidence genes to OrthoFinder v2.5.4 (Emms and Kelly, 2019). As, high-confidence genes, we considered those genes for which expression/functional evidence was available and that had full-length without internal stop-codons. Hence, annotated partial/pseudogenes were ignored in our study. Furthermore, we annotated transposon (TE) or TE-related genes in all four potato genomes using the approach described by Jayakodi et al. (2020). Briefly, this approach involves two stages to annotate TEs in all annotated genes. The first stage involves searching for keywords and PFAM IDs related to TEs in the functional annotation obtained from AHRD and classify each gene as either TE or non-TE gene. The second stage involves combining the orthogroup (OG) information for each gene obtained from OrthoFinder with the curated genes of the first stage. In the last step, we evaluated whether each OG is classified as non-TE OG based on the criteria that the OG contains less than 30% of TE genes and the mean AHRD score of OG is ≥2, otherwise the OG is classified as TE OG.

Tandemly duplicated genes clusters were identified among the non-TE genes of each potato genome separately using the methodology described by Jayakodi et al. (2020). In our study, we restricted our analyses to non-TE genes with a known chromosomal location. Briefly, this approach involves the identification of homologous gene pairs present on the same chromosome using all vs. all BlastN (Altschul et al., 1990) of coding sequences (CDS) of the longest iso-form of high-confidence gene models. This is followed by finding TDG clusters based on the following thresholds: e-value cut off of 1e⁻¹⁰, bit score ratio of ≥30%, and coverage of both query and subjects of ≥50%. In this study, we defined TDGs as the homologous genes present on the same chromosome which are intervened by up to 10 genes. TDGs (or TDG pair) correspond to single genes (or gene pairs) that belong to a TDG cluster. A TDG cluster corresponds to a group of TDG pairs.

Further, the variation in density of TDGs across the respective genomes were explained by fitting a general linear model against the density of various genomic features such as all non-TE genes, DNA transposable elements (TEs), and RNA TEs using R v3.6.1.² In the next step, the residuals of these models were tested against a uniform distribution using a Kolmogorov–Smirnov (KS) test in order to evaluate whether the TDGs were distributed uniformly across the genome after adjusting for the distribution effects of the above-mentioned genomic features. In addition, Pearson’s correlation coefficients between density of TDGs and the above-mentioned genomic features were computed. All density calculations were performed in 1 Mb windows across the genome. For TDG also, 1.5 Mb windows were considered. Density was defined as the proportion of bases in each window that corresponded to the respective genomic feature. TDG clusters were categorized based on their level of sharing across four potato genomes into core, i.e., present in all four potato genomes, shared, i.e., present in more than one potato genome but absent in at least one potato genome, and private (or lineage-specific) clusters, i.e., present in a single potato genome.

Both Pfam domain and Gene Ontology (GO) term information was extracted from the output of AHRD for each potato genome. For each detected Pfam domain, we calculated the number of proteins present among the proteins of TDGs followed by performing the protein domain enrichment analysis using Fisher exact test (Fisher, 1992). FDR corrected values of p < 0.05 were considered as significant. Enrichment of GO terms was performed using GOATOOLS (Klopfenstein et al., 2018).

For this analysis, publicly available RNA-Seq datasets from NCBI SRA³ were used (Supplementary Tables S5–S8). The raw-reads were filtered for low-quality and trimmed adapter sequences using Trimmomatic v0.39 (Bolger et al., 2014) with the following parameters: (1) adapters were removed using for pair-end: ILLUMINACLIP:TruSeq3-PE.fa:2:30:10:8:true, for single-end: ILLUMINACLIP:TruSeq3-SE.fa:2:30:10:8; (2) removing leading and trailing low-quality or N bases using LEADING:3 TRAILING:3; (3) scanning the read with a four-base wide sliding window, cutting when the average quality per base drops below 20 (SLIDINGWINDOW:4:20); and (4) selecting reads with at least 50 (for pair-end: MINLEN:50) or 36 bases long (for single-end: MINLEN:36). Only, RNA-Seq datasets with at least 50% of high-quality reads after trimming were chosen for further analysis. The high-quality reads were used to estimate transcripts abundances on a gene level (i.e., averaging across the alleles present at the respective gene) using Kalisto v0.46.1 (Bray et al., 2016) with default parameters (for single-ends: -l 200 -s 20) and obtained transcript per million (TPM) values. For all examined four genomes, we classified a gene as expressed, if TPM was >0.5 in at least one RNA-Seq dataset, otherwise classified as unexpressed. In the next step, we evaluated expression breadth, i.e., the number of RNA-Seq datasets in which the gene was expressed. Further, the expressed genes were used to calculate the expression level and expression specificity. The expression level was defined as the mean value of TPM across RNA-Seq datasets for each gene. The expression specificity was measured as described by Yang and Gaut (2011) and ranged from 0 to 1, with a higher value indicating higher specificity, i.e., higher variation in expression across RNA-Seq datasets. If a gene is expressed in a single library only, the expression specificity is 1. In contrast, if a gene is expressed equally in all RNA-Seq datasets, the expression specificity is 0.

For each potato genome and each gene, we considered the non-overlapping 1 Kb sequence upstream of the transcription start site (TSS) as putative promoter sequence and retrieved it from the respective potato genomes using BEDTools v2.27.1 (Quinlan and Hall, 2010). The promoter sequences for each tandem duplicate gene pair were aligned using the matcher program of EMBOSS v6.6.0.0 (Rice et al., 2000) to compute promoter similarity (Ps). We excluded promoter sequences that contained unknown nucleotides (N).

For each TDG pair, the amino acid sequences were aligned using MAFFT v7.453 (Katoh and Standley, 2013) followed by reverse translation into nucleotide alignment using PAL2NAL v14 (Suyama et al., 2006). Finally, the nucleotide alignments were used to compute K_a, K_s, and K_a/K_s values using the Gamma-MYN method of K_aK_s_Calculator v2.0 (Wang et al., 2010). As the high number of reversions or multiple substitutions at synonymous sites reduces accuracy and reliability for rate estimation (Vanneste et al., 2013), we excluded K_s values >2 from the analysis.

Pearson’s correlation coefficient (r) was calculated across all TDG pairs between (i) expression and promoter divergence, (ii) expression divergence and age of duplicate pairs, (iii) expression divergence and K_a/K_s ratios, and (iv) promoter divergence and age of duplicate pairs using SciPy library (Virtanen et al., 2020) in Python v3.8.5.⁴

The recently published gene annotation of the diploid clone dAg derived from the elite variety Agria does not contain information about isoforms and partial gene models. Hence, we first improved the existing gene annotation of dAg using the PASA pipeline (Haas et al., 2003) by utilizing the available Iso-Seq data of Agria (tAg; Freire et al., 2021) and obtained 44,464 gene models with 58,734 isoforms. We then filtered out partial gene models using AGAT tools (Dainat et al., 2022) and obtained 39,088 full-length gene models with 53,352 isoforms (referred as full-length set in Table 1A). This new annotation contains a significantly higher number of gene models than the high-confidence annotation of DM v6.1 (Pham et al., 2020), and a slightly higher number of gene models than the annotation of M6 (Leisner et al., 2018) and RH (Zhou et al., 2020; Table 1B).

To ensure that the results of our analyses can be compared across all four potato genomes, the functional annotation for all four potato genomes was performed using the AHRD pipeline.⁵ This pipeline assigns a quality score in the form of a three-character string, where each character is either “*” if respective criteria is met or “-” otherwise, for each annotated gene to indicate how confident the assigned annotation is. The “*” in first position indicates bit score, and e-value of the blast result are >50 and 1e⁻¹⁰, respectively. The “*” in second position indicates overlap of the blast result is >60%, and the “*” in third position indicates top token score of assigned Human-Readable-Description is >0.5. In our study, we selected annotation with a quality score of at least two stars as best annotation, i.e., at least two out of three criteria should meet by the annotated gene. Overall, AHRD assigned functions to 84.31% of the genes of all four potato genomes with at least two stars, whereas, AHRD was unable to annotate 13.13% of the genes of all four potato genomes (Supplementary Table S1). We also estimated orthologs and orthogroups among the four potato genomes using OrthoFinder (Emms and Kelly, 2019) and obtained 28,647 orthogroups representing 93.5% of the genes of all potato genomes. A total of 16,107 orthogroups contained genes from all four potato genomes, of which 9,563 were single-copy orthogroups (Supplementary Table S2). Finally, we annotated 10.11% of the genes (Table 1B) of all four potato genomes as TE or TE-related genes using the approach of Jayakodi et al. (2020).

Tandemly duplicated genes clusters were identified among the non-TE genes of each potato genome. In total, 2,090, 1,867, 1,661, and 1,832 TDG clusters were identified in dAg, DM, M6, and RH, respectively. The percentage of annotated genes in TDG clusters of the total number of genes were 18.67, 18.52, 16.83, and 16.06% in dAg, DM, M6, and RH, respectively (Table 2). The availability of multiple high-quality de novo potato genome assemblies allowed us to determine the consistency of various characteristics of TDG clusters across potato genomes. The identified TDGs were dispersed throughout the genome and shown to have a similar distribution across the four potato genomes (Figure 1). Moreover, the density distribution of TDGs with both 1 and 1.5 Mb sliding-windows resulted in same density distribution pattern across the four potato genomes (Figure 1; Supplementary Figure S5). Overall, the density of TDGs across the genome was significantly (value of p < 2.2e⁻¹⁶) associated with the density of non-TE genes across the respective potato genomes. The correlation coefficients were 0.6, 0.68, 0.61, and 0.63 for dAg, DM, M6, and RH, respectively. The same observation was made for the densities of both DNA and RNA TEs. After correcting for these densities using a general linear model, the density of TDG showed a significant (value of p < 2.2e⁻¹⁶) deviation from a uniform distribution in each potato genome. The KS test statistics (D) are 0.38, 0.38, 0.41, and 0.57 for dAg, DM, M6, and RH, respectively. Chromosome 1 of all potato genomes harbored the highest number of TDG clusters, while Chromosomes 6, 5, 8, and 7 harbored the lowest number of TDG clusters in dAg, DM, M6, and RH, respectively (Supplementary Table S3). A similar distribution in terms of the size of TDG clusters was observed across the four potato genomes, and the majority of the TDG clusters comprised two genes (Figure 2A). Further, the majority of TDG pairs within TDG clusters did not contain intervening genes, and moreover, DM contained the highest number of TDG pairs without intervening genes among the four potato genomes (Figure 2B). For dAg, a higher proportion of TDGs with two exons was observed compared to the other three genomes. For the latter, the proportion of TDG with single exons was higher compared to that of non-tandemly duplicated genes (Supplementary Figure S1).

Tandemly duplicated genes account for about 18% of the total non-TE genes. Thus, it would be interesting to gain insights into the evolutionary forces that affect the TDGs. Therefore, we examined the sequence divergence in TDGs of each potato genome by estimating K_a (number of substitutions per nonsynonymous site), K_s (number of substitutions per synonymous site), and K_a/K_s ratios, and compared their distributions across the four potato genomes. We observed pronounced peaks at 0.1, between 0.1 and 0.15, as well as between 0.3 and 0.4 for K_a, K_s, and K_a/K_s distributions, respectively, for all potato genomes (Figures 3A–C). Further, all four potato genomes showed a higher K_a/K_s values compared to K_a values (Figures 3D,F). An average of 92.97% TDG pairs showed K_a/K_s < 1.0 (i.e., negative or purifying selection), while an average of 6.47% showed a K_a/K_s value >1.0 (i.e., positive selection). In addition, we observed that the cultivated potato genotype dAg contained the highest number of TDG pairs (916) that were under positive selection (i.e., K_a/K_s > 1), while the wild potato genotype M6 contained the least number (345; Supplementary Table S4).

Furthermore, we investigated the functional specificities of the TDGs by an enrichment analysis to answer whether evolutionary forces drive TDGs in potato toward a specific biological function. First, we used the Pfam protein domain information for an enrichment analysis (DEA). A total of 83.42, 91.47, 86.73, and 78.68% of the identified TDGs contain Pfam domains in dAg, DM, M6, and RH, respectively. Further, a total of 917, 805, 706, and 756 unique Pfam domains were identified in TDGs of dAg, DM, M6, and RH, respectively. Across the four potato genomes, TDGs showed a similar distribution of the number of protein domains they harbor. The majority of TDGs contained a single Pfam domain only (Supplementary Figure S2). The DEA identified that 46, 69, 60, and 61 protein domains in TDGs of dAg, DM, M6, and RH, respectively, were significantly (FDR < 0.05 and a minimum number of 10 TDGs/protein domain) over-represented. The most important protein domains that were enriched included NB-ARC, leucine-rich repeat, pathogenesis-related proteins, UDP-glucosyl transferase, glutathione S-transferase, and auxin-responsive protein IPP transferase. Interestingly, all the enriched protein domains were also differentially enriched (z-score) across the four potato genomes (Table 2 and Supplementary Figure S3). Further, a total of 65.21, 66.66, 50, and 39.34% of enriched protein domains were present in positively selected (K_a/K_s > 1) TDGs of dAg, DM, M6, and RH, respectively (Figure 4A). As alternative approach to identify functional specificities of the TDGs, we performed a GO-term enrichment analysis (GOEA) to identify over-represented gene ontology (GO) terms in TDGs. GOEA identified 190, 230, 199, and 154 statistically significant (FDR < 0.01) biological processes (BP) in dAg, DM, M6, and RH, respectively (Table 2). Further, a total of 56.84, 70.43, 63.81, and 48.05% of enriched BP were present in positively selected (K_a/K_s > 1) TDGs of dAg, DM, M6, and RH, respectively. The majority of the top 30 BP were associated with defense responses against various biotic conditions (bacteria, fungus, and virus), and stress responses against various abiotic (hypoxia, cadmium, heat, light, and UV-B) stress conditions (Figure 4B). In line with the domain enrichment, the enriched BP were also differentially enriched (fold enrichment) across four potato genomes (Figure 4B).

In this study, we examined patterns of expression divergence between TDGs in four potato genomes using publicly available RNA-Seq datasets of the respective potato genomes except for dAg. In the public domain, only one RNA-Seq dataset was available for tAg, the tetraploid ancestor of dAg, but the dataset did not pass our selection criteria after trimming out low-quality reads to include in the analysis. Consequently, the RNA-Seq datasets generated under various stress conditions (drought, salt, heat, and cold) belonging to different potato cultivars were used to estimate expression of dAg (Supplementary Tables S5–S8). We used log10 transformed transcripts per million (TPM) values obtained from Kallisto (Bray et al., 2016) as a proxy for expression levels. In the next step, we classified each TDG as expressed if TPM > 0.5 in at least one RNA-Seq dataset, otherwise classified as unexpressed. Based on this criterion, we found that 68.82, 81.14, 76.98, and 74.35% of TDGs were expressed in dAg, DM, M6, and RH, respectively (Table 2). Across all potato genomes, TDGs showed higher expression specificities than non-TDGs (Figure 5A), while both expression breadth and levels were lower than non-TDGs (Figures 5B,C).

For each potato genome, we classified each TDG pair as expressed TDG pair if both gene copies were expressed, otherwise classified as unexpressed TDG pair. Based on this criterion, we found 46.57, 67.24, 64.29, and 69.08% of TDG pairs were classified as “expressed TDG pairs” in dAg, DM, M6, and RH, respectively (Supplementary Table S9). Further, for each potato genome, we selected expressed TDG pairs and calculated Pearson’s correlation coefficient (r) between expression profiles of both gene copies. As comparison, we calculated r for the same number of randomly selected non-TDG pairs. The TDG pairs of DM, RH, and dAg revealed the same approximate normal distribution of correlation coefficients as the control sample (Figures 6A,C,D), while the TDG pairs of M6 showed for both a distribution of the correlation coefficients that deviated from a normal distribution (Figure 6B). Further, the 95% quantile of the correlation coefficient r of randomly selected gene pairs was used as threshold for determining if the two gene copies of a TDG pair have diverged expression, i.e., if r > = 95% quantile, then the TDG pair have a similar expression, while r < 95% quantile, then the TDG pair have diverged expression. Based on this criterion, 79.56, 69.95, 74.55, and 73.66% of expressed TDG pairs showed diverged expression in dAg, DM, M6, and RH, respectively (Supplementary Table S9).

As shown above, the TDG pairs exhibited significant transcriptional divergences that prompted us to undertake a systematic investigation of variation present in their promoters. In order to do that we retrieved non-overlapping 1 Kb sequence upstream of the transcription start site for each gene of a TDG pair as a putative promoter sequence and measured promoter similarity (Ps). We measured Ps for the same number of randomly selected non-TE gene pairs of respective potato genomes to represent the background level of Ps that is expected to be observed by chance. We found a similar distribution in Ps of tandemly duplicated gene pairs across four potato genomes (Figure 7). On average, Ps for randomly selected gene pairs was 0.159, 0.152, 0.146, and 0.139%, while Ps for TDG pairs was 0.34, 0.29, 0.25, and 0.31% in dAg, DM, M6, and RH, respectively (Supplementary Table S10). As mentioned above, 95% quantile in the Ps of randomly selected gene pairs was used to classify the promoter sequences of TDG pairs as either conserved or diverged. Based on this criterion, 66.06, 70.83, 75.61, and 66.59% of TDG pairs showed diverged promoters in dAg, DM, M6, and RH, respectively (Supplementary Table S10).

As shown above, expression, age of TDG pairs measured as K_s, and their associated promoters exhibited significant similarities as well as divergences in TDGs. It would be interesting to know the correlations among them. Therefore, to test whether the divergence of promoter similarities correlates with the expression divergence, we computed Pearson correlation coefficients (r) between the expressed TDG pairs against their promoter similarities. The correlation coefficients were low but significantly positive and ranged for the four genomes from 0.09 (value of p = 6.73 × 10⁻¹⁰) for dAg to 0.19 (value of p = 2.54 × 10⁻³⁹) for RH (Figure 8A). Further, to test whether promoter divergence correlates with coding sequence divergence measured as K_s, we computed Pearson correlation coefficients between K_s of TDG pairs and their respective promoter similarities. The correlation coefficients were low but significantly negative and ranged from −0.28 (value of p = 1.93 × 10⁻¹⁹⁸) for DM to −0.23 (value of p = 7.47 × 10⁻⁷⁸) for M6 (Figure 8B). However, the distribution of promoter divergence across K_s (Figure 8B) suggests a continuous expression divergence has occurred over the evolutionary time (K_s) between TDG pairs. To verify this, we computed Pearson correlation coefficients between K_s of TDG pairs and their expression. Similarly, a continuous expression divergence over the time between TDG pairs occurred with a significantly negative correlation ranging from −0.17 (value of p = 4.97 × 10⁻³⁰) for M6 to −0.07 (value of p = 6.72 × 10⁻⁰⁷) for RH (Figure 8C). Further, we computed Pearson correlation coefficients between K_a/K_s ratios of TDG pairs and their expression to reveal the type of selection that caused the divergence in expression between the duplicated genes of a TDG pair. The results indicated that the divergence in gene expression between duplicated genes across TDG pairs was due to purifying selection (Figure 8D). Only for M6 a positive correlation of 0.08 (value of p = 5.5 × 10⁻⁸) was observed.

A total of 7,450 TDG clusters were identified across four potato genomes. To determine if TDG clusters were shared across the four potato genomes, we used the orthology information that linked the non-TE gene models of all four potato genomes and categorized them into core, shared, and private (or lineage-specific). Based on this categorization, on average, 25.02, 29.94, and 45.03% of all TDG clusters were private, core, and shared clusters, respectively, across the four potato genotypes (Figure 9A; Table 3). In addition, the non-cultivated potato genotype DM contained the highest proportion of shared clusters (51.96%), while the cultivated potato genotype dAg contained the lowest proportion of shared clusters (40.57%). Conversely, the cultivated genotype dAg contained the highest proportion of private clusters (32.92%), while the non-cultivated potato genotype DM contained the lowest proportion of private clusters (18.37%; Table 3). An average of 52.04% of Pfam protein domains enriched in all TDG clusters was present in private clusters. In addition, the private clusters of the cultivated potato genotype dAg showed with the highest proportion of enriched Pfam protein domains (about 70%), while the private clusters of the non-cultivated potato genotype RH revealed the lowest proportion of enriched Pfam protein domains (about 41%). Furthermore, an average of 62.24% of Pfam protein domains which were enriched within the private clusters were present in positively selected TDGs (K_a/K_s > 1). In addition, the cultivated potato genotype dAg contained a high proportion of Pfam protein domains (about 66%) which were enriched within the private clusters and that were present in positively selected TDGs (K_a/K_s > 1). In contrast, the wild potato genotype contained the lowest proportion of Pfam protein domains (about 59%) which were enriched within the private clusters and that were present in positively selected TDGs (K_a/K_s > 1; Supplementary Table S11). In general, private clusters showed low-expression specificities but higher expression breadth compared to shared and core clusters. The private clusters were shown to have lower K_a and K_s values but higher K_a/K_s values compared to shared and core clusters (Figure 9C). Further, we observed that the majority of private clusters localized toward the centromere, while the majority of core and shared clusters tended to localize toward the ends of chromosome-arms (Supplementary Figure S4).

Tandem duplications are widespread in plant genomes and contribute significantly to the evolution of genomes. By using the available well-annotated multiple de novo genome assemblies of potatoes, we observed that the TDGs in potato were dispersed throughout the genome with similar distribution across four potato genomes. TDGs accounted for about 18% of all non-TE genes in potatoes. This number is considerably higher than in rice (about 15.1%; Jiang et al., 2013), maize (average of 10.3%; Kono et al., 2018), and pear (about 11.1%; Qiao et al., 2018). The differences in proportions of TDGs among potatoes, rice, maize, and pear might be generated by species-specific gene duplications as observed across 141 plant species (Qiao et al., 2019). The distribution and chromosomal localization of TDGs of potato observed in our study are similar to the TDGs of rice (Jiang et al., 2013) and maize (Kono et al., 2018). Further, our results indicate that there is variation among potato genotypes for the content of TDGs in the genome (Table 2). The higher number of TDGs in dAg compared to the other three genotypes is likely due to the availability of a more accurate and larger genome assembly. For example, more than 20% of annotated non-TE genes of M6 were present on the unknown chromosome, and hence these genes were excluded from the analysis which might in turn lead to the identification of the lowest number of TDGs among the four potato genomes. However, further studies with larger number of genotypes are required in order to link the observed variation in the content of TDGs in the genome to the history of the examined genetic material.

Gene duplication is a mechanism that creates functional innovation and novelty in the genome. Here, we explored the relationships between TDGs and functional specificities across cultivated and wild genotypes of potatoes. Protein domain enrichment revealed enrichment of several important protein domains related to genes involved in disease resistance [NB-ARC (PF00931) and leucine-rich repeat (PF00560; Jupe et al., 2012; Prakash et al., 2020); pathogenesis-related proteins (PF00407; Lakhssassi et al., 2020)], stress-responsive [UDP-glucosyl transferase (PF00201; Rehman et al., 2018), glutathione S-transferase (PF02798; Islam et al., 2018)], auxin-responsive protein (PF02519; Jain et al., 2006), and various biosynthetic pathways [IPP transferase (PF01715; Lindner et al., 2014)] in TDGs across four potato genomes. Our results highlighted that the potato inhibitor 1 family protein domain (PF00280) containing genes were enriched only in the wild potato genotype (i.e., M6) and were shown to be under positive selection, i.e., K_a/K_s > 1 (Figure 4A). The potato inhibitor 1 family protein domain containing genes are naturally occurring plant serine proteinase inhibitors. They act as both endogenous and defense-related plant regulators in potato under wounding and nematode infection (Turrà et al., 2009). In line with protein domain enrichment, GO enrichment also revealed that these TDGs were involved in biological processes related to defense responses against various pathogens, such as bacteria, fungi, and virus, stress responses against various abiotic stress conditions (such as light, auxin, cadmium, heat, and UV), and biosynthetic pathways (lignin, saponin, phenylpropanoid, di-terpenoid, flavonoid, glucosinolate, and indole). Moreover, both protein domains and GO processes were differentially enriched across four potato genomes and TDGs encoding these functional specificities were under positive selection, i.e., K_a/K_s > 1 (Figures 4A,B). These results suggested that the bias in functional specificities coupled with positive selection might play an important role in the retention of TDGs in potatoes (Shiu et al., 2006; Ren et al., 2014). This finding is consistent with a previous study that found that the retention of TDGs favors genes involved in certain important functions to maintain the fitness of the organism (Blanc and Wolfe, 2004; Edger and Pires, 2009).

Potato underwent at least two rounds of genome duplication, 185 and 67 million years ago (Potato Genome Sequencing Consortium et al., 2011), and retained 6,078, 5,817, 4,748, and 5,018 TDGs for dAg, DM, M6, and RH, respectively. Our study highlighted a number of striking patterns in sequence, expression, and regulatory divergences between gene copies of TDG pairs across four potato genomes (Figures 3, 5–8). Based on these patterns, we propose and distinguish multiple models such as sub-functionalization (Force et al., 1999), genetic redundancy (Panchy et al., 2016), and neo-functionalization (Ohno, 1970) that may contribute to the retention of TDGs in potato genomes.

Based on the orthology information, we found a significant proportion (an average of 25.02% of all TDG clusters) of private or lineage-specific TDG clusters across four potato genotypes (Figure 9A; Table 3). The majority of them localized in pericentromeric regions which was not observed for core and shared clusters (Supplementary Figure S4). The reason for this observation might be the same that is responsible for an over-representation of presence absence variation (PAV) genes in Arabidopsis thaliana (Tan et al., 2012) in pericentromeric regions. The low extent of recombination in those regions of the genome might prevent the spread of present TDGs in a population, and thus the TDGs remains private. Our results highlighted that the tandem duplication generates a significantly varying proportion of private clusters across four potato genomes (Table 3). In addition, we found that the cultivated genotype dAg contained a high proportion of enriched Pfam protein domains which are present in positively selected TDGs of private clusters compared to non-cultivated as well wild potato genotypes (Supplementary Table S11). This observation might be due to breeder’s selection to combine positive alleles for many traits. These results may also indicate that the tandem duplication generates lineage-specific TDGs with functional bias between evolutionarily closed species, such as the four potato genotypes, which is similar to that of generation of lineage-specific TDGs with functional bias between evolutionarily distant plant species (Hanada et al., 2008).

In general, our results highlight that the private TDG clusters showed a lower expression specificity and higher expression breadth compare to the shared and core clusters (Figure 9B). Our observation indicates that the private clusters were involved in tissue-specific functional specificities. This result is in contrast to results of legume species (Xu et al., 2018) where private TDGs showed higher expression specificity and lower expression breadth. The reason for that remains elusive. Furthermore, an average of 30.99% of private TDG pairs showed divergence in expression and have K_a/K_s < 1.0 (Supplementary Table S14) indicating purifying selection at the nucleotide level which in turn suggests that sub-functionalization of duplicated genes across tissues has been occurred to retain the duplicated genes across four potato genotypes (Force et al., 1999; Lynch and Force, 2000). In addition, a vast majority of these retained private TDG pairs (an average of 84.29%) had an identical annotation in terms of protein domains between duplicated gene copies of respective TDG pairs (Supplementary Table S14). These results together reinforced that the retention of a majority of TDGs of private clusters were occurred through sub-functionalization.

In addition, an average of 19.05% of private TDG pairs showed similarity in expression profiles, of which a majority of them (an average of 61.6%) are under purifying selective pressure (i.e., K_a/K_s < 1.0), and a vast majority of them (an average of 87.62%) contained identical annotation in terms of protein domains between duplicated genes of respective private TDG pairs, across four potato genotypes (Supplementary Table S14). These results indicated that these private TDG pairs might have been retained through selection for genetic redundancy that may be beneficial in a way that is similar to fail-safe in an engineered system (Hanada et al., 2009; Zhang, 2012; Panchy et al., 2016). Alternatively, these private TDG pairs might have been retained simply because there has been insufficient time for one copy to be removed, because they are evolving relatively neutrally (Panchy et al., 2016).

We also found that an average of 0.97% only of private TDG pairs have different annotation in terms of protein domain composition with K_a/K_s > 1.0 between gene copies of respective TDG pairs across four potato genotypes (Supplementary Table S14). These results indicate that these 0.97% of private TDG pairs might have been retained through neo-functionalization (Ohno, 1970). In addition, a total of eight enriched Pfam protein domains are present in these retained 0.96% of private TDG pairs in all potato genotypes and these enriched protein domains were mainly involved in important biological processes such as disease resistance (NB-ARC domain: PF00931; Leucine rich repeat: PF13855; Jupe et al., 2012; Prakash et al., 2020) and self-incompatibility (S-locus glycoprotein domain: PF00954; Xing et al., 2013; Supplementary Table S15).

Our results indicated that the tandem duplication contributed to lineage-specific expansion of several gene families across potato genotypes. For example, NBS-LRR, Cytochrome P450, UDP-glucosyl transferase, and 2OG-Fe (II) oxygenase gene families were differentially expanded by tandem duplication across potato genotypes (Figure 4A). Furthermore, the GO enrichment revealed a functional bias of TDGs across the four potato genotypes (Figure 4B). This is supported by recent studies of specific gene families in potato (Herath and Verchot, 2020; Liu et al., 2020; Yang et al., 2020; Xuanyuan et al., 2022) and provided an important source for genetic diversity in plants for adaptive evolution against various environmental stimuli. These results are similar to a previous study conducted on two maize genotypes (such as B73 and PH207) where more than 49% of B73’s and 40% of PH207’s TDGs were lineage-specific (Kono et al., 2018). Furthermore, the importance of lineage-specific expansion of TDGs was also studied in A. thaliana against various abiotic stress stimuli and found a strong correlation between tandem duplication and abiotic stress conditions (Hanada et al., 2008). Thus, the lineage-specific expansion of gene families by tandem duplication coupled with functional bias might significantly contribute to potato’s genotypic diversity. However, to understand their effect on phenotypic characters requires further research.