Development of homozygous tetraploid potato and whole genome doubling-induced the enrichment of H3K27ac and potentially enhanced resistance to cold-induced sweetening in tubers

Abstract Polyploid plants typically display advantages on some agronomically important traits over their diploid counterparts. Extensive studies have shown genetic, transcriptomic, and epigenetic dynamics upon polyploidization in multiple plant species. However, few studies have unveiled those alternations imposed only by ploidy level, without any interference from heterozygosity. Cultivated potato is highly heterozygous. Thus, in this study, we developed two homozygous autotetraploid lines and one homozygous diploid line in parallel from a homozygous diploid potato. We confirmed their ploidy levels using chloroplast counting and karyotyping. Oligo-FISH and genome re-sequencing validated that these potato lines are nearly homozygous. We investigated variations in phenotypes, transcription, and histone modifications between two ploidies. Both autotetraploid lines produced larger but fewer tubers than the diploid line. Interestingly, each autotetraploid line displayed ploidy-related differential expression for various genes. We also discovered a genome-wide enrichment of H3K27ac in genic regions upon whole-genome doubling (WGD). However, such enrichment was not associated with the differential gene expression between two ploidies. The tetraploid lines may exhibit better resistance to cold-induced sweetening (CIS) than the diploid line in tubers, potentially regulated through the expression of CIS-related key genes, which seems to be associated with the levels of H3K4me3 in cold-stored tubers. These findings will help to understand the impacts of autotetraploidization on dynamics of phenotypes, transcription, and histone modifications, as well as on CIS-related genes in response to cold storage.


Introduction
Polyploid plants typically display enlarged organ, altered architecture, enhanced growth vigor, and improved tolerance to abiotic stresses [1][2][3], as well as novel phenotypes compared to their diploid progenitors [4], which thus may not only provide them a fitness advantage for short-term adaptation [3], but also offer a potential breeding approach to increase yields in crops without gene manipulation [5][6][7]. Polyploids typically include autopolyploids and allopolyploids [3,8]. Autopolyploids are considered to arise within a single species by the doubling of homologous genome, whereas allopolyploids arise through the merging of genomes from different species (hybridization) and subsequent doubling [3,8].
A number of allopolyploid plant species, including important crops, have been well studied, such as bread wheat [9], Brassica juncea [10], strawberry [11], and cotton [12]. Investigations of allopolyploidization in such plant species mainly focus on genetic variation, changes in gene expression and DNA methylation [9,10,12,13]. Since the merge of distinct subgenomes affecting allopolyploid offspring on both ploidy level and heterozygosity, it is challenging to discuss the effect of ploidy level on such changes alone without any interferences from heterozygosity. The example of maize suggests that gene expression is altered by both hybridization and ploidy change [14]. Therefore, whole-genome doubling (WGD) of homozygous diploid plants may provide an ideal system to assess the ploidy level on genetic, gene expression, and epigenetic variations. Compared to allopolyploidization, few studies have been dedicated to investigating the impacts of autopolyploidization on genetic and epigenetic changes [15][16][17][18][19]. In addition, most of these studies only focus on alternation of DNA methylation as the epigenetic impact caused by autopolyploidization, whereas evidence for changes in histone modifications are limited, such as active histone marks H3K4me3 [20] and H3K27ac [21], as well as repressive histone mark H3K27me3 [22]. Thus, the dynamics of histone modifications upon autopolyploidization remain elusive.
Potato is the most important non-grain food crop. After harvesting, potato tubers are usually stored at cold temperature to control diseases and prevent sprouting. However, along with the cold storage, starch is broken down into reducing sugars (glucose and fructose) in potato tubers. This process is called cold-induced sweetening (CIS) [23]. The accumulation of reducing sugars further results in an unsatisfactory dark color and accumulation of acrylamide in processing potato products with high temperature [24][25][26]. CIS is a long-lasting challenge to the potato industry. Successful attempts to improve the resistance to CIS are achieved mainly by silencing/mutating a key gene called vacuolar invertase (VINV) using RNAi [26], TALENs [27], and CRISPR-Cas9 [28]. However, the regulation of a number of other key genes involved in CIS is unclear [29]. Polyploid plants are often associated with better tolerance to abiotic stresses over their diploid progenitors [3], such as cold [30]. Thus, understanding the expression and epigenetic regulation of CIS-related genes at different ploidy levels may provide clues for improving resistance to CIS in potato.
Cultivated potato (Solanum tuberosum, 2n = 4x = 48) is highly heterozygous [31]. However, a doubled monoploid, DM1-3516 R44 (S. tuberosum group Phureja, 2n = 2x = 24, referred to as DM 1-3) with its homozygous and well-sequenced genome [32] has been a valuable source for genetic, cytogenetic, and epigenetic studies. DM 1-3 and its derived tetraploids provide an ideal system and allow us to study the impacts of ploidy levels on genetic and epigenetic variations on a genome-wide scale, with minimal interference from heterozygosity. In this study, we developed two tetraploid potato lines from DM 1-3 through WGD using callus culture. Both tetraploid derivatives displayed phenotypic polymorphisms, including larger tuber size, over the parallel diploid line, which failed in WGD. Chloroplast counting and karyotyping confirmed their ploidy levels. No structural variations were observed between the tetraploid and diploid lines at chromosome level using Oligo-FISH. Genome re-sequencing revealed limited DNA sequence variations among these lines, indicating these potato lines are nearly homozygous. Interestingly, each homozygous tetraploid line displayed ploidy-related differential expression for various genes. We discovered a genome-wide enrichment of H3K27ac in genic regions upon WGD. However, that enrichment was not associated with differential gene expression between two ploidies. The tetraploid lines may exhibit better resistance to CIS than the diploid line, indicated by lower levels of reducing sugars in cold tubers. The differences in reducing sugar content among three potato lines are coincident with the expression of key genes involving in CIS, which seems to be associated with the levels of H3K4me3. Our study will provide useful germplasm and help to understand the impacts of WGD on multiple aspects from phenotypes, chromosomes, to genetic and epigenetic variations, as well as on CIS-related genes in response to cold storage.

Development of tetraploid potato lines from the homozygous diploid DM 1-3
The homozygous diploid potato DM 1-3 and its homozygous tetraploid derivatives would provide a great opportunity to investigate the impact of ploidy level on phenotypes, gene expression and histone modifications, with minimal interferences from heterozygosity. Therefore, we developed two tetraploid potato lines (DM4X-13 and DM4X-17) with doubled chromosome number from DM 1-3 using callus culture [33]. In addition, we selected one diploid line (DM2X-new) in parallel, which failed in chromosome doubling, to eliminate the inf luence from the hormone treatment during callus culture.
Typical variations in phenotypes from different ploidy levels were observed, including increased plant height, stem diameter, leaf size and upright plant architecture with less branch number in both DM4X-13 and DM4X-17, compared to those in DM2X-new ( Fig. 1a and b; Table S1, see online supplementary material). Both tetraploid lines displayed larger anthers and styles than DM2Xnew (Fig. S1, see online supplementary material). These morphological differences between the diploid and tetraploid potato lines are similar to other plant species, such as tomato [34], Arabidopsis [35], and watermelon [36].
Tetraploid potato tends to produce fewer but larger tubers than diploid ones [37,38]. A similar trend was also observed between the tetraploid lines and DM2X-new ( Fig. 1c; Table S1, see online supplementary material). Both DM4X-13 and DM4X-17 produced significantly fewer tubers (Fig. S2a, see online supplementary material) and lower total tuber yield than DM2X-new (t test, P < 0.01) (Fig. S2b, see online supplementary material). However, the average tuber weight of both tetraploid lines, particularly for DM4X-17, tends to be higher compared to that of DM2X-new ( Fig. S2c, Table S1, see online supplementary material). Thus, the elevation of the ploidy level from diploid to tetraploid may be associated with an increase in the single tuber weight, rather than the total yield.

Cytogenetic verification using chloroplast counting and oligo-FISH
Because the average ploidy level can be quickly estimated by using the chloroplast number in the guard cells of a stoma [39], we evaluated these newly developed potato lines and found that the chloroplast number of DM4X-13 and DM4X-17 ranged from 15-21 and from 16-26, respectively (Fig. S3, see online supplementary material), while the chloroplast number of DM2X-new varied from 8 to 14, which is similar to the diploid progenitor DM 1-3 with chloroplast number 7-14 ( Fig. S3, see online supplementary material). This result indicates that the number of the chromosomes is likely doubled in both DM4X-13 and DM4X-17, while it remains at diploid level in DM2X-new. We further confirmed the ploidy level for each line by karyotyping. Both DM4X-13 and DM4X-17 contain 48 chromosomes ( Fig. 2a and b), while DM2X-new has 24 chromosomes (Fig. 2c). In addition, we employed Oligo-FISH with potato chromosome-specific probes [40] to evaluate structural variations at chromosome level between the diploid and tetraploid lines. The FISH results displayed no visible changes in signal localization among DM2X-new, DM4X-13, and DM4X-17 (Fig. 2d), suggesting both tetraploid lines contain the full set of the doubled genome and likely there is no large chromosomal structural variation among these lines.

Sequence variation evaluated by genome re-sequencing
To further identify variation at DNA sequence level, genomes from the three newly developed lines were re-sequenced (Table S2, see online supplementary material). A total of 277 million (M) raw pair-end (PE) reads were obtained for these lines. Approximately, 90% of the raw reads were mapped to the genome assembly DM 1-3 (v4.04) [41], covering 20X for DM2X-new, 31X for DM4X-13 and 32X for DM4X-17, respectively. DM2X-new displayed 94.1% homozygosity, while DM4X-13 and DM4X-17 showed 89.8% and 89.3% homozygosity, respectively (Fig. 2e), suggesting these three potato lines are nearly entirely homozygous. Compared to the reference genome of DM 1-3, we identified a limited number of sequence variants in DM2X-new, including 11 217 structural variations (SVs), 195 121 single nucleotide polymorphisms (SNPs) and 2128 insertion/deletions (Indels) (Fig. S4, centromere regions indicated [42], see online supplementary material). In the tetraploid lines, we detected a slightly larger number of sequence variants (Fig. S4, see online supplementary material). Noticeably, in DM2X-new, 6.0% (2358/39400) and 4.7% (1860/39400) of the potato protein-coding genes were associated with sequence variants in the putative promoters In addition, comparison between DM2X-new and the tetraploid lines revealed that only a small number of genes (2.2% of 39 400) harboring sequence variants in exons were unique to DM2X-new, Together, these results indicate that there are limited differences in DNA sequence between the newly developed diploid and tetraploid lines, particularly in the putative promoters and the coding regions.

Impact of WGD on gene expression
To explore any variation in gene expression pattern between the diploid and tetraploid lines, we conducted RNA-seq for each line using room-temperature stored tubers (named RT tubers) ( Table S3, see online supplementary material). Gene expression level between two biological replicates from each line was highly correlated (Pearson's correlation R 2 = 0.96-0.98, P < 2.2e-16) (   supplementary material). It is interesting that the differentially expressed genes (DEGs) from the two independent WGD events (DM2X-new vs. DM4X-13; DM2X-new vs. DM4X-17) overlapped only in a small proportion (six up-regulated genes and nine down-regulated genes), resulting in a number of DEGs that were unique to each WGD event (Fig. 3), leading us to suggest that ploidy elevation from diploid to tetraploid does not affect gene expression in tubers in consistent ways across the tetraploid lines. Comparison between the original line DM 1-3 and the tetraploid lines also revealed a considerable number of DEGs (14%-30%) were unique to each WGD event. Similar patterns are evidenced by WGD-wild potato species [2].
Interestingly, investigation of the six up-regulated genes in both WGD events (Fig. 3) revealed that four genes have been reported to be involved in abiotic stresses [43][44][45]. In addition, gene ontology (GO) analysis for 110 genes, which were specifically up-regulated in DM4X-17 in the two WGD events (Fig. 3), showed that most of them were enriched in functions related to responses to stresses and stimulus (Fig. S7, see online supplementary material). Similar up-regulation of gene expression related to enhanced stress tolerance was found in autotetraploid Mulberry [19] and autotetraploid Birch [46]. Polyploid plants are thought to be associated with improved ability to tolerate harsh environments [47]. Thus, these results may indicate that elevation of ploidy from diploid to tetraploid enhances the ability for potato to tolerate multiple stresses, especially for DM4X-17.

WGD induced the enrichment of active histone modification H3K27ac
We conducted chromatin immunoprecipitation followed by sequencing (ChIP-seq) for each potato line, with antibodies against active histone marks of H3K4me3 and H3K27ac, as well as a repressive histone mark H3K27me3 (Table S4, see online supplementary material). Comparison of the histone modification levels in genic regions of all potato genes revealed that DM4X-13 displayed generally higher levels of H3K4me3 than DM2X-new (Fig. 4a). Both tetraploid lines showed higher levels of H3K27ac than DM2X-new (Fig. 4b), while less difference in the levels of H3K27me3 was observed (Fig. 4c), indicating the tetraploid lines seem to be enriched with the active histone modifications in genic regions. In addition, the diploid progenitor DM 1-3 was associated with the lowest levels of H3K4me3 among all tested lines (Fig. S8a, see online supplementary material), while the levels of H3K27ac in DM 1-3 fell between DM2X-new and the tetraploid . It is interesting that a marked difference in genomic distribution of H3K27ac peaks was found in exons and intergenic regions between two ploidies (z test, P < 0.01) (Fig. 4d). Approximately, two-thirds more H3K27ac peaks from the tetraploid lines (29% for DM4X-13 and 32% for DM4X-17) were observed in exons comparing to that from DM2X-new (18%) (Fig. 4d). In contrast, the distribution of both H3K4me3 and H3K27me3 peaks was generally similar between DM2X-new and the tetraploid lines ( Fig. 4e and f).
We further observed a considerable amount of H3K27ac peaks specific to DM4X-13 (53.7%, 27 759) and DM4X-17 (35.3%, 8097) in two WGD events (Fig. S10a, see online supplementary material). These specific peaks as well as their cognate genes overlapped only in a small proportion between the two tetraploid lines (21.3% of peaks and 26.5% of genes relative to DM4X-13). A similar trend was also found for H3K27ac peaks located in exons (Fig. S10b, see online supplementary material). In contrast, only a small proportion of H3K4me3 and H3K27me3 peaks were found to be specific to the tetraploid lines ( Fig. S10c and d, see online supplementary material). Thus, it is indicated that WGD triggers the emergence of a considerable amount of H3K27ac in genic regions, but a limited number of H3K4me3 and H3K27me3 enriched regions. Unexpectedly, in comparison to DM2X-new, 99.9% (12 406/12420) and 99.5% (4189/4209) of genes associated with tetraploid-specific H3K27ac peaks were either silenced or constitutively expressed in DM4X-13 and DM4X-17, respectively. Collectively, these results suggest that WGD, at least for potato, may induce the enrichment of the active histone modification H3K27ac, which seems to have limited effect on alternation of transcription.

WGD-induced differential expression tends to be independent to sequence variations
To interrogate whether differential gene expression between ploidies is potentially associated with variation in DNA sequence, we identified sequence variants (SVs, indels and SNPs) in exon regions as well as putative promoters of the DEGs from each WGD event. In comparison with DM2X-new, we found sequence variants in only a single (3.7% of 27) DEG in DM4X-13, while the remaining 26 DEGs (96.3%) were not associated with any sequence variants in the corresponding regions (Fig. 5a). Similarly, we found six (3.8% of 158) DEGs were associated with sequence variants in DM4X-17 compared to DM2X-new, while the sequence from both regions of 152 (96.2% of 158) DEGs in DM4X-17 were identical to that in DM2X-new (Fig. 5b). This result indicates that variation in DNA sequence does not tend play a major role in the regulation of the differential gene expression between the homozygous diploid and tetraploid lines. Further observation on histone modifications for these DEGs without any sequence variants revealed that the alternation of the expression between two ploidies tends to be independent to the enrichment of the histone modifications (Fig. S11, see online supplementary material), indicating these DEGs might be regulated by other epigenetic modifications rather than the tested ones.

Dynamics of gene expression in response to cold storage
To investigate any differences in gene expression pattern in response to cold storage of tubers between the diploid and tetraploid lines, we developed additional transcriptome data from tubers stored under cold temperature (4 • C) for each line (named cold tubers) (Table S3, see online supplementary material). The overall gene expression in cold tubers between two biological replicates of each line was highly correlated (Pearson's correlation r = 0.81-0.93, P < 2.2e-16) (Fig. S6a, see online supplementary material). In DM2X-new, we found that a total of 2201 and 2302 genes were statistically significantly up-and down-regulated upon cold storage, respectively (Fig. S12, see online supplementary material). A slightly larger number of DEGs (2750 up-regulated genes and 2736 down-regulated genes) was observed in DM4X-13 upon cold storage, while only 1292 upregulated genes and 1481 down-regulated genes were identified in DM4X-17 (Fig. S12, see online supplementary material). We categorized these DEGs into four groups based on their expression behavior upon cold storage, including genes: G1, up-regulated in three lines (n = 900) (Fig. 6a); G2, down-regulated in three lines (n = 1075) (Fig. 6b); G3, up-regulated in the tetraploid lines, but not up-regulated (down-regulated, constitutively expressed and silenced) in the diploid line (n = 138) (Fig. 6c); and G4, up-regulated in the diploid line, but not up-regulated in the tetraploid lines (n = 476) (Fig. 6d).
We found a large proportion of up-(G1) and down-regulated genes (G2) behaved consistently in the three lines upon cold storage ( Fig. 6a and b), which was substantially more than DEGs from the other two groups (Fig. 6c and d), suggesting the direction of the differential gene expression in response to cold storage tends to be similar among the three lines regardless of the ploidy level. The GO terms for both G1 and G2 were mainly related to organelles and ribonucleoprotein complex ( Fig. S13a and b, see online supplementary material), indicating that cold storage appears to affect the function of the basic cell components. Noticeably, DM4X-13 showed the largest changes in gene expression in both G1 and G2 upon cold storage than the other two lines ( Fig. S14a and b, see online supplementary material), whereas DM4X-17 displayed the least variation in gene expression ( Fig. S14a and b, see online supplementary material). The difference in the degree of the differential expression among the three lines was also supported by individual coldresponsive genes [48][49][50][51][52][53][54] (Fig. S14d, see online supplementary material).
For G3, these genes displayed larger changes in expression in DM4X-13 than in DM4X-17 in response to cold storage (Fig. S14c, see online supplementary material). They were found to be enriched in functions related to response to stress and metabolic processes (Fig. S13c, see online supplementary material). Two heat-shock protein genes (PGSC0003DMG400003122 and PGSC0003DMG400028744) from G3, whose closest homologs participate in responding heat stress in peppers [50], tomatoes [51], potato [51,52], and cold stress in soybeans [53], were also detected to be significantly up-regulated in the tetraploid lines (Fig. S14e, see online supplementary material). Therefore, it is indicated that the tetraploid lines might be associated with better tolerance to cold stress to some degree.
In G4, 476 DEGs were found to be only up-regulated in DM2Xnew compared to the tetraploid lines upon cold storage (Fig. 7d). Over 98% of them were either constitutively expressed or silenced in both tetraploid lines. The GO analysis revealed that these DEGs were mainly related to thylakoids (Fig. S13d, see online supplementary material). Reduction in efficiency of photosynthesis can be frequently observed in multiple plant species as a mean of protection for the photosynthetic apparatus in response to environmental stress, such as cold stress for Vitis vinifera [55], Cistus albidus and Quercus ilex [56]. However, the potato tuber is not a primarily photosynthetic organ [57]. The relatively stable gene expression in the tetraploid lines might be associated with less greening in tuber skin.

Dynamics of histone modifications in response to cold storage
We conducted additional ChIP-seq with antibodies against H3K4me3, H3K27ac, and H3K27me3, separately, using the coldstored tubers from each line (Table S4, see online supplementary material). Upon cold storage, both DM2X-new and DM4X-17 displayed overall increased levels of active histone modifications (H3K4me3 and H3K27ac) around TSSs of the DEGs and constitutively expressed genes (Fig. 7a-h). In contrast, DM4X-13 showed generally reduced levels of both active histone modifications around TSSs of these genes (Fig. 7i-l). For the repressive mark  19.56% (10 457/53462) peaks were specific to cold tubers from DM2X-new, DM4X-13, and DM4X-17, respectively (Fig. S16a, see online supplementary material). In contrast to H3K4me3, a large number of cold-specific peaks were discovered for H3K27ac (39%-72%) (Fig. S16b, see online supplementary material) as well as for H3K27me3 (46%-59%) (Fig. S16c, see online supplementary material) in each line. Thus, cold stress induces a considerable number of H3K27ac and H3K27me3 enriched regions, but a limited number of H3K4me3 enriched regions, indicating H3K27ac and H3K27me3 may play a regulatory role in response to cold storage. Furthermore, comparing to cold-specific H3K4me3 peaks from DM2X-new, a number of cold-specific H3K4me3 peaks were only detected in tetraploid lines (n = 3990 and n = 7181), and a large proportion of these cold-specific peaks were unique to each tetraploid line (n = 2579 and n = 6407) (Fig. S16a, see online supplementary material). Similar patterns were observed for H3K27ac (Fig. S16b, see online supplementary material) and H3K27me3 (Fig. S16c, see online supplementary material), suggesting that the tetraploid lines were associated with diverse histone modifications in response to cold storage.
We further analysed the changes of histone modifications for genes from G1-G4 upon cold storage and found that each group displayed similar patterns with the least variations of H3K4me3 ( Fig. S17a- The cold-induced variation in histone modifications was not associated with the changes of the gene expression, suggesting that other epigenetic modifications may potentially play a major role in the regulation of differential expression in response to cold.

Resistance of cold-induced sweetening in the tetraploid lines
To evaluate the degree of CIS in the newly developed potato lines, we measured the contents of starch, sucrose, and the reducing sugar (glucose and fructose) in both RT and cold tubers using high-performance liquid chromatography (HPLC). Upon cold storage, the starch content was reduced in both diploid and tetraploid lines (Fig. 8a), while the sucrose and the reducing sugar were increased (Fig. 8b and c). Coincident with cold-induced degradation of starch, most of the starch synthesisrelated genes were down-regulated upon cold storage (Fig. 8d),  Figure 7. Profiles of histone modifications associated with potato genes in response to cold storage. The levels of H3K4me3 (a, b) and H3K27ac (c, d) associated with potato genes in RT tubers (a, c) and cold tubers (b, d) of DM2X-new. The levels of H3K4me3 (e, f) and H3K27ac (g, h) associated with potato genes in RT tubers (e, g) and cold tubers (f, h) of DM4X-17. The levels of H3K4me3 (i, j) and H3K27ac (k, l) associated with potato genes in RT tubers (i, k) and cold tubers (j, l) of DM4X-13. Genes were divided into 100 bins and aligned together from TSSs to TTSs. Genes f lanking regions (±1 kb) were analysed in 100 bins. Histone modification signal was measured by quantifying the mid-point of ChIP-seq PE reads within each bin and normalizing by ChIP-seq read number per bp genome region per million mapped reads.
while the starch degradation-related genes were generally up-regulated (Fig. 8d). In addition, the genes related to reducing sugar synthesis were overall up-regulated (Fig. 8d). For reducing sugar degradation, Glucose 6-Phosphate/phosphate Transporter (GPT2.1, PGSC0003DMG400005269), a highly tuber-specific expressed gene involved in transporting Glucose 6-Phosphate (G6P) [29], was substantially down-regulated in DM2X-new and DM4X-13, but remained similarly expressed in DM4X-17, detected by both RNA-seq (Fig. 8d) and quantitative real-time PCR (qRT-PCR) (Fig. S18, see online supplementary material). Collectively, these results indicate that cold storage triggers CIS in both diploid and tetraploid lines possibly through the differential expression of the key genes involved in the carbohydrate metabolism pathway.
Although both diploid and tetraploid lines contained similar level of starch in both RT and cold tubers (Fig. 8a), the accumulation of the sucrose and reducing sugar was obviously different upon cold storage (Fig. 8b and c). In cold tubers, DM4X-13 contained the highest amount of sucrose, while DM 4X-17 had the lowest amount of sucrose (Fig. 8b). The Sucrose Synthesis 4 (SuSy4) (PGSC0003DMG400002895), a major tuber-specific SuSy gene [58,59], encode proteins that catalyze the breakdown of the sucrose in tubers [60]. The sucrose contents of these lines are coincident with the expression levels of the SuSy4 in cold tubers, with transcripts per kilobase million (TPM) = 345, 201, and 3472 in DM2X-new, DM4X-13, and DM4X-17, respectively ( Fig. 8d; Fig. S18, see online supplementary material). Both tetraploid lines accumulated less reducing sugar than DM2Xnew in cold tubers (t test P < 0.05) (Fig. 8c). In addition to SuSy4, the degradation of sucrose into glucose and fructose is also mediated by vacuolar acid invertase (VINV), a key enzyme involved in CIS [61]. Coincident with the accumulation of the reducing sugar in cold tubers of each line, the expression of VINV (PGSC0003DMG400013856) showed the highest level (TPM = 266) in DM2X-new compared to that in the tetraploid lines (TPM = 178 in DM4X-13 and TPM = 73 in DM4X-17) (Fig. 8d; Fig. S18, see online supplementary material). Collectively, the reduction in accumulation of the reducing sugar in tetraploid lines implies that these homozygous tetraploid lines might be associated with better resistance to CIS than the diploid line through the differential expression of the key genes in the starch metabolism pathway.
We further examined if the differential expression of the tested carbohydrate metabolism-related genes is potentially associated with sequence variation. None of these genes were associated with any types of sequence variants in their putative promoters and genic regions, except for the putative promoters of the Starch Synthesis I (SS1) and the Debranching Enzyme 9 (DBE9). Therefore, the difference in the cold-induced differential expression of these genes among the three potato lines is likely independent to the genetic variation observed here. Investigation of histone modifications for the key genes SuSy4 and VINV revealed their expression levels were associated with the levels of H3K4me3. In cold tubers, the H3K4me3 levels of SuSy4 in DM2X-new (3.97), DM4X-13 (2.8), and DM4X-17 (4.45) were coincident with the expression levels of the gene in each line (Fig. S19, see online supplementary material). Similarly, VINV, which displayed the highest levels of H3K4me3 in cold tubers of DM2X-new, was expressed at the highest level than the tetraploid lines (Fig. S19, see online supplementary material). Collectively, we speculate that the homozygous tetraploid lines may display better resistance to CIS via epigenetic-mediated ectopic expression of some key CIS genes upon WGD.  Table S5.

Phenotypic variations upon WGD
WGD in both autopolyploids and allopolyploids can lead to immediate phenotypic, genetic, and epigenetic changes and such polyploids might have adaptive advantage compared to their diploid progenitors [3]. In this study, both autotetraploid potato lines were associated with various phenotypic changes compared to their parallel diploid line ( Fig. 1; Table S1, see online supplementary material), which are typical differences commonly found between diploid and polyploid from other plant species [34][35][36]. General observations show that tetraploid potato cultivars are associated with higher yield than diploid ones [7]. However, the total tuber yield of two autotetraploid lines in this study was significantly lower than that of the diploid line (Fig. S2b, Table S1, see online supplementary material), although both autotetraploid lines produced larger but fewer tubers ( Fig. 1c; Table S1, see online supplementary material). Thus, the elevation of the ploidy level from homozygous diploid to tetraploid might contribute to the average tuber size rather than the total yield. A recent pioneer study in hybrid potato showed that F1 hybrids produced significant higher single plant tuber yield than their homozygous diploid parents [62], indicating the heterosis may impact the tuber yield. Thus, heterozygous diploid potato may produce comparable tuber yield to the tetraploid one. Deciphering the regulation of key genes between ploidies as well as between diverged subgenomes may provide clues for breeding targets.

WGD impacts epigenetic modifications in the autotetraploid lines
The autotetraploid lines did exhibit sequence differences compared to the diploid line upon WGD (Fig. S4, see online supplementary material). Thus, the autotetraploid lines harbor more genetic diversity over the diploid line. Interestingly, 96% of the DEGs between two ploidies lacked any sequence polymorphisms in their putative promoters and coding regions (Fig. 5), suggesting epigenetic modifications may play a role in regulation of differential expression. However, the variations of all tested histone modifications appear to be not associated with the direction of the differential expression (Fig. S11, see online supplementary material), indicating the DEGs are likely to be regulated by other epigenetic modifications. DNA methylation has been discovered to play a regulatory role in differential expression between ploidies in multiple plant species, such as autotetraploid apple [63], Aegilops tauschii [64], Artemisia annua [65], Arabidopsis thaliana [66], and rice [17], implying that it is possible that the differential expression between the autotetraploid and diploid potato lines might be potentially mediated by differential DNA methylation.
Both autotetraploid lines displayed similar enrichment of H3K27ac in genic regions upon WGD (Fig. 4b). Histone modification H3K27ac is associated with transcription activation in plants [21]. However, both up-and down-regulated genes showed increased levels of H3K27ac upon WGD (Fig. S11a, see online supplementary material), suggesting the independence of H3K27ac to the alternation of transcription. Histone acetylation generally enhances chromatin accessibility by neutralization of the basic charge in histones [67,68]. In plant, H3K27ac is enriched and peaked around TSSs, where chromatin is more accessible [69,70]. We speculate that WGD-induced enrichment of H3K27ac might provide more accessible chromatin environment facilitating the access of transcription factors (TFs).

WGD impacts gene expression unpredictably between each tetraploid line
An interesting question raised by this study is why WGD-induced phenotypically similar autotetraploid individuals were associated with divergent gene expression. Noticeably, the ploidy-related DEGs in DM4X-17 (n = 158) are close to six times more than the DEGs in DM4X-13 (n = 27) (Fig. S6b, see online supplementary material), and these DEGs overlapped in a very small proportion between the tetraploid lines (Fig. 3). Thus, WGD may affect expression of different genes in each autotetraploid lineage. Alteration of gene expression upon polyploidization has been well studied in both animal and plant systems and the alternation tends to be species-specific [3]. For instance, in Arabidopsis, two accessions (C24 and Ws) exhibited opposite direction of the expression change upon autotetraploidization [15]. Intriguingly, four autotetraploid individuals derived from two wild diploid potato species, respectively, displayed line-specific expression patterns with a large proportion of DEGs specific to each autotetraploid line [2], which is similar to our findings in the homozygous potato lines. Coincidently, the GO terms of the ploidy-related DEGs from the wild potato study as well as from our current study were mainly enriched in functions related to responses to stresses and cell components. None of these terms were associated with plant architecture, development, or growth. One possible explanation for the discrepancy between similar phenotypes and diverse gene expression between autotetraploid individuals would be that key genes related to core functions for development remain genetically and epigenetically identical between the derived autotetraploid individuals.

Autotetraploid lines may be associated with better resistance to CIS
WGD-induced polyploids were associated with enhanced tolerance to abiotic stresses, such as cold [30], drought [16], and salt [17]. Similarly, the autotetraploid potato lines in this study may be associated with better tolerance to cold stress. We found that a number of genes involved in the carbohydrate metabolism pathway displayed diverse pattern of differential expression among the three potato lines in response to cold stress (Fig. 8d). The differences in differential expression pattern among the three lines are coincident with their amount of sucrose and reducing sugar contents (Fig. 8b and c). It is noted that cold tubers from both autotetraploid lines contained lower levels of reducing sugars than cold tubers from the diploid line, while the elevation of reducing sugars upon cold storage was less in the autotetraploid lines (Fig. 8c), implying that the autotetraploid lines may be associated with better tolerance to CIS. The differences in sucrose and reducing sugars contents among the three potato lines seem to be coincident with their differences in differential expression of key genes for sucrose degradation and reducing sugar synthesis, such as SuSy4 and VINV ( Fig. 8d; Fig. S18, see online supplementary material). In addition, their expression tended to be associated with the level of histone modification H3K4me3. Thus, understanding the epigenetic regulation of such key genes involved in CIS between ploidies may provide clues for the improvement of resistance to CIS.

Plant materials
Tetraploid potato lines were developed from the homozygous diploid potato DM 1-3516 R44 (S. tuberosum Group phureja, 2n = 2x = 24) using callus culture [33] with hormones of NAA, BA, and GA 3 . The newly developed tetraploid lines and their diploid counterpart were subject to phenotyping, cytogenetic validation, genome re-sequencing, RNA-seq, and ChIP-seq. All potato lines were grown in a greenhouse under the photoperiod of 16 h/24 • C daylight and 8 h/18 • C darkness. Tubers were harvested 80 days after planting in the soil and were grouped for treatments at room temperature (22 • C) and cold temperature (4 • C) for 14 days, respectively. Upon temperature treatment, tubers were ground into fine powder under liquid nitrogen and stored at −80 • C for experiments.

Phenotyping
Plant architecture, plant height, stem diameter, branch number, leaf area, and f lower were evaluated for all potato lines 60 days after planting. Tuber yield, number, and size were measured at harvest (80 days after planting).

Chloroplast counting
The number of chloroplasts was scored in two guard cells of a stomata. Eight stomata were examined for each plant. Brief ly, epidermal peels were taken from the abaxial side of leaves from 60-day old plants and immediately put on a slide with 1% iodine solution for staining. Cover slip was mounted after two minutes staining. Chloroplasts were then counted using a light microscope at 100 and 400× magnification.

Chromosome counting and oligo-FISH
Root tips were cut from young potato plants and immediately treated with nitrous oxide under a pressure of 160 psi for 1 h for chromosome condensation. Fixation for root tips was performed immediately in the ethanol/acetic acid (3:1) solution, followed by a thorough wash three times using autoclaved distilled water. Root tips were then digested with 2% cellulase (Yakult Pharmaceutical, Tokyo, Japan) and 1% pectinase (Sigma, St. Louis, MO, USA) at 37 • C for 90 minutes. The digested root tips were applied for chromosome preparation following the dropping method [71]. The number of chromosomes was counted in 30 individual cells from at least five plants for each line. The potato specific Oligo-FISH probes [40] were labeled using digoxigenin [72] and the hybridization was performed following the standard protocol [73]. Signal detection was conducted using rhodamine anti-digoxigenin with the DAPI antifade solution (Vector Laboratories, Newark, CA, USA). FISH images were captured using a f luorescence microscope (Zeiss Axio Scope A1) with a CCD camera (ORCA Flash4.0).

Starch and sugar contents analysis
Both RT and cold tubers from each potato line were analysed for the starch, sucrose, glucose, and fructose contents. Starch was quantified based on the GB5009.9-2016 method. Sugar (sucrose, glucose, and fructose) contents were measured using HPLC [74], with the refractive index detector (RID) at 40 • C, an amino column at 35 • C, 10 μL injection volume, and 1 mL/min f low rate.

RNA-seq and ChIP-seq
Tuber samples, which were used for starch and sugar contents analysis, were used for both RNA-seq and ChIP-seq. Tubers harvested from six individual plants were pooled as a biological replicate for each line. Two biological replicates of RNA-seq libraries were constructed for each temperature treatment and sequenced using an Illumina NovaSeq 6000 platform with 150 bp PE mode. All Illumina sequencing in this study was conducted using the same platform and mode in Novagene.

Reverse transcription and qRT-PCR
cDNA was synthesized using All-in-One First-Strand cDNA Synthesis SuperMix (EasyScript, #AE341-02) with oligo(dT) 20 . qRT-PCR for each gene was performed with SYBR Green I T5 Fast qPCR Mix (TsingKe, #TSE202), using a RT-PCR cycler (CFX Connect Bio-Rad). The qRT-PCR parameters were set with initial denaturation of 94 • C for 30 s, and 45 cycles of 94 • C for 5 s, 54 • C for 15 s, and final extension of 72 • C for 10 s. Three technical replicates and two biological replicates from each treatment were used for quantifying the expression. The reference gene EF1α was employed to normalize the expression of each gene, using the calculation of 2 -Ct . The gene-specific primers used for each gene are listed in Table S6 (see online supplementary material).
Similarly, cleaned ChIP-seq reads were mapped to the same genome assembly [41] using BWA 'mem' algorithm [81]. Reads were retained for further analyses if they were mapped to a unique position in the genome. The histone modification signal was defined as the mid-point of the uniquely mapped PE reads. The level of a histone modification within an interval was quantified by calculating histone modification signals and normalizing to length of the interval, PE read number per million uniquely mapped reads, as well as IgG. Histone modification enriched regions were identified using MACS2 sofware [82].

Gene ontology enrichment
Potato genes were used to find homologous sequences with the highest similarity in Arabidopsis, using Blastp program (BLAST 2.9.0). The enriched GO terms were screened with the program of agriGO v2.0 [83], using Arabidopsis homologous protein sequences. Fisher's exact test was applied for the GO enrichment test with Benjamini-Hochberg FDR P-value normalization. GO terms from all annotated Arabidopsis genes were set as background for each enrichment test.

Variant calling
Genomic DNA isolated from young leaves was subject to library construction and sequencing. The raw reads were first cleaned using fastp program with default parameters, followed by the alignment to the PGSC v4.04 [41] using the BWA 'mem' [81] with default parameters. The mapped reads were further filtered using Samtools [84] with Q20 followed by 'markdup' [84] with default parameter for PCR duplication removal. The GATK program (https://gatk.broadinstitute.org) was used to detect SNPs and InDels. SNPs were hard-filtered with QD less than 2.0, MQ less than 40.0, FS greater than 60.0, SOR greater than 3.0, MQRankSum less than −12.5, ReadPosRankSum less than −8.0. InDels were filtered with QD less than 2.0, FS greater than 200.0, SOR greater than 10.0, MQRankSum less than −12.5, ReadPosRankSum less than −20.0. InDels included insertion or deletion with length of 1 to 10 bp. Structural variations were detected using Delly software [85] with default parameters. Structural variations included insertion, deletion, inter/intra chromosomal translocation and inversion.