Skip to main content
BMC is moving to Springer Nature Link. Visit this journal in its new home.
Download PDF

High resolution diel transcriptomes of autotetraploid potato reveal expression and sequence conservation among rhythmic genes

BMC Genomics volume 26, Article number: 925 (2025) Cite this article

Abstract

Background

See More

Photoperiodic changes in diel cycles of gene expression are pervasive in plants. The timing of circadian regulators, together with light signals, regulate multiple photoperiod-dependent responses such as growth, flowering or tuber formation. However, for most genes, the importance of cyclic mRNA levels is less clear. We analyzed the diel transcriptome of modern cultivated potato, a highly heterozygous autotetraploid. Clonal propagation and limited meiosis have led to the accumulation of deleterious alleles, making tetraploid potato an ideal model system to investigate the conservation of cyclic expression and cyclic genes during artificial selection and clonal propagation.

Results

Our results indicate that rhythmic alleles of cultivated potato are more highly expressed than non-rhythmic genes and are highly co-expressed not only under diel cycles but also across tissues, developmental stages, and stress conditions. Moreover, the smaller ratio of non-synonymous to synonymous differences within rhythmic versus non-rhythmic allelic groups indicates that cyclic genes, in general, have more conserved core functions than non-cyclic genes. In accordance with this observation, fully rhythmic allelic groups are highly enriched in photosynthesis and ribosome biogenesis genes, which have core functions in plants. Furthermore, we investigated differences in cyclic expression patterns between photoperiods identifying potential regulators for the strong changes in phase of expression of ribosome biogenesis and pathogen response genes. Finally, analyses of genes involved in tuber formation suggests that the regulation of CO gene transcription is not the only factor enabling tuberization under long days in modern cultivated potato.

Conclusions

This study not only provides high quality diel transcriptomic datasets of cultivated potato but also provides important insight on the role of allelic diversity in rhythmic expression in plants.

Peer Review reports

Background

Cultivated potato, Solanum tuberosum L. Group Tuberosum (2n = 4x = 48), is a heterozygous autotetraploid that is vegetatively propagated. Polyploidy has been associated with domestication and there is evidence that polyploidy preceded domestication for many crops [1]. In addition, variation in circadian clock-associated genes has been linked to crop domestication [2,3,4,5] indicating that cyclic gene expression is an important factor targeted in domestication.

Most of the previous work on cyclic transcriptomics in polyploid crops has relied on haploid genome assemblies or on allopolyploids where variability can be attributed to ancestry such as Brassica rapa and hexaploid wheat. B. rapa is a mesohexaploid, generated through an ancient three-species hybridization event after its split from Arabidopsis 13–43 million years ago, followed by subsequent diploidization [6, 7]. Hexaploid wheat (AABBDD) formed through interspecific hybridization between a domesticated tetraploid and a diploid species 8,500–9,000 years ago [8]. In B. rapa, 42% of paralog pairs have differential expression patterns under light/dark conditions, and most display differences in their median level of expression [9]. In hexaploid wheat, 64.15% of triads also display differences in rhythmicity patterns [10]. These results indicate functional divergence of these cyclic homeologs, which might be partly caused by their different origins. There is no information on how rhythmic expression correlates across alleles in autopolyploids such as potato.

Potato species have undergone more recent polyploidization events. Sexual polyploidization via 2n gamete production has been recurrent in wild as well as cultivated potato species [11]. Recent studies indicate that modern potato cultivation was initiated with tetraploid Andean landraces and that the breeding process has included admixture with tetraploid Chilean genotypes and wild species [5, 12]. Due to its poor fertility caused by different factors including polyploidy [11], cultivated potato is primarily maintained vegetatively, leading to an accumulation of deleterious alleles and inbreeding depression, resulting in a high genetic load [13]. Thus, tetraploid potato is a good model system to investigate the conservation of cyclic expression and cyclic genes throughout artificial selection and extended clonal propagation.

Modern cultivated potatoes are able to grow in a wide range of latitudes [14] indicating flexibility to photoperiod adjustment for both physiological and developmental processes. In photoperiod sensitive potatoes, tuberization is inhibited in days longer than a critical day length due to elevated CONSTANS (CO) gene expression. COs act as activators of SP5G (SELF PRUNING 5G), which is in turn represses SP6A, encoding a phloem mobile protein that induces tuber formation. Under long day conditions, CDF1, an inhibitor of CO, is actively degraded by the F-box E3 ubiquitin ligase FLAVIN-BINDING KELCH REPEAT F-BOX 1 (FKF1). Thus, only under short-day conditions,when CDF1 proteins are allowed to accumulate, is repression on SP6A released, and thus tuberization is able to occur. In modern cultivars, truncated CDF1 alleles, which encode proteins that no longer interact with FKF1, accumulate and repress CO expression under both short- and long-day conditions, thus allowing for tuberization even under long photoperiods [5, 14]. It is believed that the circadian clock is also involved in this process by generating oscillations of expression of FKF1, CDF1 and CO in an analogous manner to the regulation of flowering time in Arabidopsis. This current model of photoperiod control of tuberization was developed in S. tuberosum Group Andigena, which requires short day conditions to tuberize. It is not clear how photoperiod affects the expression of these genes in cultivated potato adapted to higher latitudes that are able to tuberize under long days yet still retain some photoperiod sensitivity.

The goal of this study was to perform a comprehensive analysis of allele specific rhythmic expression in an autotretraploid plant using cultivated potato as a model. Our previous studies on cultivated potato showed that both leaves and tubers are able to maintain robust diel and free running rhythms [15]. However, these experiments were performed using unphased assemblies of the highly heterozygous tetraploid potato using 3'Tag-Seq, and thus lack the ability to distinguish the expression of the individual alleles. Here, we generated high time resolution datasets of diel expression under long and short photoperiods on the cultivated potato S. tuberosum cv. “Atlantic” using a haplotype resolved genome assembly and annotation. These diel datasets together with a comprehensive potato gene expression atlas were used to investigate allele specific expression of rhythmic genes under different photoperiods, tissues, and stress conditions. Taking advantage of our high resolution diel transcriptomes and circadian reporter lines, we also studied tissue specific phase changes between leaves and tubers. Finally, we characterized the expression of an extensive list of circadian, photoperiod and tuberization related genes. We believe this work is a valuable resource for the understanding of potato growth and development.

Methods

Plant growth and sample collection conditions for diel datasets

Solanum tuberosum cv. "Atlantic"(referred to herein at "Atlantic") tissue culture plantlets were transplanted to soil (Suremix, Michigan Grower Products, Galesburg, MI) in 14 cm deep × 9 cm square pots and fertilized with Peter’s 20–20-20 weekly. Plants were grown in a BioChambers High-Light Hi/Lo CO2 GRC-40 growth chamber under either short days (SD) [12 h light (22°C, 350 mol s−1 m−2), 12 h dark (18°C)], or long days (LD) [16 h light (22°C, 350 μmol s−1 m−2), 8 h dark (18°C)]. Sampling occurred every 2 h for 24 h. Dawn samples (ZT 0, ZT 24) were collected in the dark and dusk samples (ZT 12) were collected in the light. At each time point, tissue was collected from three plants (replicates). For leaf tissue, three terminal leaflets were collected from the third and fourth newest fully expanded leaves per replicate. Under short days, tissue was also collected from the three largest tubers per plant, using a 5 mm punch and discarding epidermal tissue. Tubers were dig out of the soil just before tissue collection. For samples collected during the day tubers were harvested in the light, for night samples tubers were collected in the dark. All tissue was snap-frozen in liquid nitrogen and stored at –80°C prior to RNA isolation. At the time of harvest, short day plants were 9 weeks old and long day plants were 4.5 weeks old.

Identification of syntenic allelic groups

Syntelogs were identified between the genomes of S. tuberosum Group Phureja DM 1–3 516 R44 (DM) v6.1 [16], Atlantic v3 [17], Solanum chacoense v5 [17], Solanum candolleanum v1 [17], and Solanum lycopersicum M82 v1 using GENESPACE (v1.2.3) [18] with the representative working model proteins for each genome downloaded from SpudDB [19] and setting the ploidy to 1, 4, 1, 1, 1 respectively. The syntelog sets were extracted from the pangenome database file from each GENESPACE run. Further processing was performed using a custom python script to select allelic groups with four or fewer alleles, requiring that alleles were found on the same chromosome but without duplication on a single haplotype. For genes with fewer than four alleles in phased chromosomes, we selected alleles from unphased chromosomes when available. Allelic groups used in this study can be found in Dataset S1.

Total RNA isolation, library preparation, and sequencing of diel expression datasets

Total RNA was isolated using a modified hot borate method [20] followed by treatment with TURBO DNase (Invitrogen, Waltham, MA) to remove any residual DNA. RNA integrity was verified via gel electrophoresis and Fragment Analyzer (Agilent Technologies, Santa Clara, CA). Short day RNA-Seq libraries were prepared using the Illumina Stranded mRNA Prep, Ligation kit (San Diego, CA) with IDT for Illumina Unique Dual indexes (Coralville, IA), and sequenced on the NovaSeq 6000 in paired end mode 100 nt by the Research Technology Support Facility Genomics Core at Michigan State University. RNA-Seq libraries for long day samples were prepared using the PerkinElmer NEXTFLEX Rapid Directional RNA-Seq Kit 2.0 with NEXTFLEX Poly(A) Beads 2.0 and NEXTFLEX® RNA-Seq 2.0 Unique Dual Index Barcodes (Waltham, MA) and sequenced on the NovaSeq 6000 in paired end mode 100 nt or 150 nt by the Texas A&M AgriLife Research: Genomics and Bioinformatics Service.

RNA-Seq data processing of diel expression datasets

To check for contamination in sequencing data, Kraken2 v1.2 [21] was run using the pre-built database k2_pluspfp_20220908 [22]. RNA-Seq reads were then cleaned using Cutadapt v4.1 [23] to discard poor quality bases and adapter sequences by providing 3’ adapter sequences for trimming from both read one and two (-a, -A) along with the following parameters: -m 70 -q 30 –trim-n –times 2 -l 100. FastQC [24] and MultiQC [25] were both used to visualize quality and additional sequencing metrics for all libraries before and after cleaning. Kraken2 results, percent of reads discarded by Cutadapt cleaning, percent GC content, and over-represented sequence(s) outliers were flagged and individual libraries discarded based on results. After Cutadapt cleaning, if a library had greater than 40 million read pairs then reads were subsampled using seqtk v1.3-r106 [26] with the sample function and the parameters -s100 40000000. Transcript abundances were quantified using Kallisto quant v0.48.0 [27] with the Atlantic v3 high confidence representative gene models [28], the –rf-stranded option, and a k-mer size of 21. Pearson correlation coefficients were calculated between biological replicates using R “pairwise.complete.obs”, method “pearson”.

Atlantic developmental gene expression atlas

A replicated developmental and stress gene expression atlas of the S. tuberosum Atlantic spanning 16 tissues and 8 treatments was used to examine expression breadth of rhythmic genes [19, 29] (NCBI BioProject PRJNA753086). Expression abundances of cleaned subsampled reads were determined using Kallisto quant (v0.4.8) [27] using the Atlantic high confidence representative gene models as described in [29].

Normalization of RNA-seq experiments

RNA-seq read counts were normalized using the rlog function from the R package DEseq2 [30]. We defined expressed transcripts as those with at least one sample with an rlog > 0. All three photoperiodic datasets (short day leaf, short day tuber, and long day leaf) were co-normalized to allow for comparisons in expression level across samples. Tissue and stress samples from the Developmental Gene Expression Atlas were normalized separately. The samples, Closed Flower, Cold Leaf Control, Hooked Stolon, Immature Fruit, Mature Fruit, Open Flower, Root Control, Sprout, Stem Control, Swollen Stolon, Tuber S1, Tuber S2, Tuber S3, Tuber S4, Tuber S5 and Young Leaf 10am, were included in the "Tissue" dataset and samples Methyl Jasmonate Control, Methyl Jasmonate Treatment, BTH Control, BTH Treatment, Salt Leaf Control, Salt Leaf Treatment, Cold Leaf Control, Cold Leaf Treatment, Drought Leaf Control, Drought Leaf Treatment, Heat Control, Heat Treatment, Root Control, Salt Root Treatment and Drought Root Treatment included in the "Stress" dataset. The diurnal time course samples from the Developmental Gene Expression Atlas were not included in this study. Normalized gene expression can be accessed via Dataset S2-4.

Rhythmicity analysis

For the analysis of rhythmicity, rlog-normalized expression values were analyzed using the JTK-CYCLE algorithm in the MetaCycle r-package [31] using the parameters of minper = 20, max_per = 26, and adjustPhase ="predictedPer". Transcripts were considered rhythmic if the JTK p-value was less than 0.001. Phase of rhythmic gene expression were used as provided by JTK when phase < 24, and circularized (phasecircularized = phaseJTK – 24) when phase ≥ 24. Rhythmicity data can be found in Dataset S5-7.

Determination of correlation of gene expression

Normalized gene expression profiles generated as described in “Rhythmicity Analysis” were further normalized using z-scoring using the python package scipy (1.11.3) [32] to compare for the daily patterns of expression rather than the abundances. z-score normalization was performed within each individual time course or experiment, rather than across experiments. The Pearson correlation of the z-scored expression patterns was calculated for each pair of alleles using the series.corr() method in the python package pandas (2.1.1) [33]. When classifying allelic pairs by rhythmicity for the correlations based on the tissue and stress datasets from the Atlantic Developmental Gene Expression Atlas, the rhythmicity values from the LD datasets were used since these experiments were performed under 15 h light/9 h dark conditions. Correlation data be found in Dataset S8-10.

Differential gene expression analysis

Differential gene expression between photoperiods was performed using the R-package DEseq2 [30]. For this analysis, the time component of RNAseq samples was disregarded, and samples were assigned to the photoperiod from which they originated. Transcripts were considered differentially expressed if the Benjamini–Hochberg false discovery rate adjusted p-value was < 0.05 and the abs(log2Fold Change) > 1. Differential expression data be found in Dataset S11-12.

Tissue specificity analysis

We calculated the Tau index as described previously [34]. We defined expressed transcripts as those with an average rlog across replicates greater than zero in at least one tissue. We set tissues with rlog equal or smaller than zero as not expressed, i.e. equal to zero.

Determination of synonymous and nonsynonymous substitution rates

We used SynMap within the CoGe online software suit [35] using default settings to calculate synonymous (Ks) and nonsynonymous (Kn) substitution rates between syntenic alleles in Atlantic [36]. We performed syntenic analysis using SynMap using the annotation file ATL_v3.working_models.gff3 available in [17]. We filtered syntenic allelic pairs using our previous syntenic analysis. Only Ks values larger than 0.001 and smaller than 3 were used for the Kn/Ks ratio calculations.

Functional enrichment

We first determined the functional protein annotation using Mercator4 v7.0 [37, 38] and the Atlantic high confidence representative gene models (Dataset S13). Functional enrichment was done using Mercator4 BIN enrichment analysis online tool [38] that uses the MapMan4 categories. The selection of background genes for the enrichment analyses depended on specific research question and is provided in the text and/or the respective figure legend. Over-representation was calculated using a one-sided Fischer's exact test. The FDR-adjusted p-value cutoff was set to 0.05.

Promoter cis-element enrichment analysis

Identification of binding sites within Atlantic promoters was performed using the annotatePeaks.pl program from the HOMER suite [39]. In brief, candidate binding sites from equivalent transcription factors from A. thaliana were downloaded from the Plant Transcription Factor Database [40] and converted into the HOMER motif format, prior to being compared to the promoter regions of the Atlantic genome. Promoter regions were defined as 1,500 bases upstream to 100 bases downstream of the annotated transcription start sites. One-sided Fisher’s exact tests were performed using the python package scipy [32] in order to determine whether promoters of query genes were enriched for these binding sites.

Identification of circadian clock, photoperiod and tuberization associated syntenic groups

Using our computational pipeline, we associated 460 circadian clock, photoperiod and tuberization (CPT) alleles to 138 allelic groups. We determined that the main reason for the lack of allelic group identification for 39 of our CPT transcripts was copy number variation and therefore, manually associated them to'putative syntenic allelic groups'(Supplemental Table S1). For several of these genes (CONSTANS, EARLY FLOWERING 4, SELF PRUNING 5G LIKE, SUGARS WILL EVENTUALLY BE EXPORTED TRANSPORTER), syntenic allelic groups were likely not identified using our automatic pipeline as these genes had local tandem duplications in the Atlantic genome. In addition, the allelic groups of CO, ELF4 and SPA1 (SUPPRESOR OF PHYA-105) contain currently annotated non-syntenic duplicates (Supplemental Table S1). For example, in addition to the syntenic triplet of CO genes on haplotype 2 of chromosome 2 there is an additional "triplet" of CO genes, which we named CO1b, CO2b, CO3b. To investigate whether these additional haplotype specific copies were annotated in other cultivated potato genomes, we searched for non-syntenic orthologs in the genomes of Otava [41] and Cooperation-88 [42] using BLASTP. We did not find any non-syntenic copies of CPT genes in these genomes. These results indicate that these additional copies in Atlantic may be mislocalized in the current Atlantic haplotype assembly or represent cultivar-specific structural variation. After this manual curation we linked the 496 CPT gene models to 153 allelic groups in Atlantic, of which, 53% had four alleles, 25% had three and only 14% and 8% had two and one allele respectively (Supplemental Table S1). It is important to notice that gene structural annotation can differ between alleles, which can lead to apparent gene expression differences. For example, significant differences in gene models between the alleles of LNK3/4 are most likely the explanation for the strong apparent allele specific expression patterns (Supplemental Fig. 1).

Graphs and statistical analysis

Graphs were generated in R version 4.3.1, with the exception of cis-element motifs. If not otherwise indicated, statistical tests were implemented in R version 4.3.1 and described in the respective figure legends. Code to generate figures can be accessed here: https://github.com/efarre/Autotetraploid_potato_diel_transcriptome

Results

Fully rhythmic allelic groups display stronger rhythms and higher gene expression in both leaves and tubers

In order to evaluate allele specific expression, we first needed to identify syntenic allelic groups or genes across the four S. tuberosum cv. “Atlantic” (referred to herein as “Atlantic”) haplotypes. To identify these high confidence syntenic allelic groups we used synteny between Atlantic haplotypes, the doubled monoploid S. tuberosum Group Phureja DM 1–3 516 R44 (DM), the two diploid wild potato species S. chacoense and S. candeollanum, and S. lycopersicum (Dataset S1). Of these genes, 2,827 (11.9%) had a single allele in the Atlantic genome, 3,871 (16.3%) had two alleles, 5,436 (22.9%) had three alleles, and 11,599 (48.9%) had four alleles. This distribution is similar to the previous reported distribution of genes in allelic groups [28]. We defined "expressed" genes as those with at least one expressed allele.

To address the nature of differences in diel expression across haplotypes in cultivated tetraploid potato, we generated short day (SD, 12-h light/12-h dark) and long day (LD, 16-h light/8-h dark) 24-h gene expression time courses in Atlantic during the tuber bulking phase. We chose a 12-h photoperiod for our short day condition, as it is a photoperiod encountered during the growing season of cultivated potatoes and allows tuber induction in most potato species [43]. Leaf tissue was collected from soil-grown plants every 2 h and used to construct RNA-seq libraries in which the resulting reads were aligned to the phased assembly of the Atlantic genome [28]. Tuber tissue was collected only under the short-day photoperiod.

Differences in the phase of expression or overall expression level between alleles could have masked cycling expression patterns in the non-haploid resolved genome of Atlantic [28]. Using the haploid resolved Atlantic genome assembly, we first identified rhythmic transcripts using the JTK algorithm [44] (Datasets S2-4). We selected one representative transcript model for each allele of each gene using the Atlantic v3 high confidence representative gene models, and transcripts were considered rhythmic if they resulted in an adjusted p-value of less than 0.001 (Fig. 1A). We selected this stringent p-value to ensure high accuracy in the estimation of cyclic parameters, such as phase. Using these criteria, we found 30,986 rhythmic transcripts in the leaf in either one or both of our conditions, corresponding to 33.6% of all expressed transcripts. Of these strongly rhythmic transcripts, 8,591 (27.7%) were identified in both SD and LD experiments. In contrast, in tubers, only 996 (1.1%) of the expressed transcripts cycled. Even among cyclic transcripts, those cycling in the tuber had larger adjusted p-values than those cycling in leaf tissue (Supplemental Fig. 2A). Most (889) of the tuber transcripts were also expressed in leaves and 52.3% of tuber rhythmic transcripts also cycled in leaf tissue in SD (Fig. 1A).

Fig. 1
figure 1

Diel cycling rates in cultivated potato. A Number of cycling transcripts per condition. In parenthesis, the percent cycling rate. B Number of cycling genes per condition. In parenthesis, the percent cycling rate. C Distribution of cycling genes with either all expressed alleles cycling (Full rhythmicity) or at least one not cycling allele (Mixed rhythmicity). Genes with no cycling alleles are not indicated. (*) Indicate p-value < 0.01 for Chi-square tests for pairwise comparisons of gene ratios between 2 and 3, and 3 and 4 alleles in each condition. D Likelihood of one allele cycling within cyclic genes with 2–4 expressed alleles. Dark purple dots represent minimum likelihood of cycling in the respective category, e.g. for a cyclic gene with two expressed alleles the minimum likelihood of cycling is 0.5. (*) Indicate that the measured ratios are larger than the minimum likelihood based on a Chi-square test (p-value < 0.05). Chi-square tests were implements with the R function chisq.test

Full size image

We defined cyclic genes as those with at least one cyclic allele. Overall, the rate of rhythmicity in genes across different conditions was similar to the rate of rhythmicity among transcripts (Fig. 1A, B). Structural variation among the four haplotypes occurs [28] and we tested whether cyclic genes had a different number of alleles than non-cyclic ones. There was only a very minor, but statistically significant, difference in the number of alleles of cyclic and non-cyclic genes, with an average 3.5 expressed alleles for non-cyclic and 3.3 for cyclic genes (p-value < 2.2e-16, Kruskal–Wallis rank sum test).

We next asked how many genes exhibited "mixed rhythmicity", i.e. they have at least one non-cyclic and one cyclic allele. We observed that the majority of rhythmic genes with at least two expressed alleles existed in a mixed rhythmicity state, rather than being fully rhythmic (all expressed alleles were rhythmic) (Fig. 1C). Genes with more alleles were more likely to have mixed rhythmicity (Fig. 1C). A Chi-square test showed that this tendency was statistically significant in leaves but not in tubers. In tubers, there was a slightly higher percentage of genes with mixed rhythmicity, which might be due to the lower robustness of tuber cyclic expression (Supplemental Fig. 2A). These observations agree with previous studies in tetraploid potato cultivars showing that genes with more alleles have higher variation in expression [28]. However, the likelihood of a transcript being rhythmic was still higher if it belonged to a cyclic gene, such that cyclic genes contained more cyclic alleles than the expected minimum (Fig. 1D). These results indicate that rhythmicity might be relatively well conserved between alleles.

Alleles from fully rhythmic genes had more robust rhythmicity than rhythmic alleles of genes with mixed rhythmicity, as quantified using the adjusted p-value of our JTK analyses (Fig. 2A). It is possible that low expression is associated with higher noise, reducing the sensitivity of rhythmicity detection algorithms specifically for lower expressed transcripts [45]. We also observed that within allelic pairs, the non-rhythmic allele had a higher average coefficient of variation between replicates than the rhythmic allele (Supplemental Fig. 2B), which might make it difficult to detect rhythmicity. We therefore performed the comparison using only genes in which all alleles were well expressed and had small coefficients of variation between replicates (Supplemental Fig. 2C). In this selection of high expressed genes we still observed more robust rhythmicity in fully rhythmic genes. The transcripts from fully rhythmic genes, when compared to all rhythmic transcripts in multiallelic genes, were enriched in functions related to the light reactions of photosynthesis and ribosome biogenesis (Fig. 2B, Supplemental Table S2), indicating that rhythmic mRNA levels might be especially important for these functions.

Fig. 2
figure 2

Allele specific cycling in potato tissues. A Cycling strength, as determined by JTK calculated FDR adjusted mean p-value of cyclic transcripts in genes with either full or mixed rhythmicity. SD, short days; LD, long days. B Functional enrichment of transcripts within fully cyclic genes of multiallelic genes in leaves in both SD and LD compared to all transcripts that were rhythmic in leaves in both conditions belonging to multiallelic genes. C Average allele expression within genes with no cyclic alleles (Non), at least one cyclic allele (Mixed) or all alleles cycling (Full). D Degree of tissue specific expression as determined by Tau within genes with mixed rhythmicity. “Yes” indicates the allele which cycles, while “No” indicates the allele that does not. Tau ranges between 0 (constitutively expressed) to 1 (tissue specific). For C and D Kruskal–Wallis multiple groups test (p-value < 0.0001) was performed and the horizontal bar indicates the mean. For A, C and D: Wilcoxon signed-rank test with Bonferroni correction was used for pairwise comparisons, such that adjusted p-value < 0.0001 (****), < 0.001 (***), < 0.01(**), < 0.05 (*) and ns (not significant)

Full size image

Finally, diel rhythmicity was also associated with higher expression. Transcripts from cyclic genes were more highly expressed than those from non-cyclic genes (Fig. 2C). We still observed a significantly higher expression of cyclic than non-cyclic transcripts within low noise well-expressed genes (Supplemental Fig. 2D). It has been hypothesized in tetraploid wheat, that arrhythmic alleles may be silenced and therefore not rhythmic [10]. However, within mixed-rhythmicity allelic pairs selected for low noise and high expressed transcripts, the cyclic allele was, on average, more highly expressed than the non-cyclic one (Supplemental Fig. 2E). Our results indicate that there are other mechanisms causing variation in rhythmic gene expression other than silencing.

We then tested whether cyclic transcripts might also be more widely expressed across tissues and developmental stages non-cyclic transcripts. We used a publicly available allele-specific dataset of Atlantic, composed of samples from 16 different tissues and developmental stages, including stolon, tuber root, stem, leaf, flower, fruit and sprout tissues to calculate the Tau index [34]. A higher Tau index indicates stronger tissue specific expression. The non-rhythmic alleles of mixed rhythmicity allelic pairs (Fig. 2D) and non-rhythmic transcripts had slightly more tissue specific expression than rhythmic ones (Supplemental Fig. 2F). For example, in leaves ~ 8.7% of non-rhythmic transcripts had a Tau larger than 0.8, which indicates a high tissue specificity, but only ~ 5.0% transcripts rhythmic in either SD or LD did. In tubers, ~ 6.1% of non-rhythmic transcripts had a Tau larger than 0.8, but only ~ 1.2% rhythmic transcripts did.

Rhythmic alleles display high correlation of expression across multiple tissues and stress conditions

We investigated whether there were differences in the timing of expression between rhythmic alleles. Most fully rhythmic allelic pairs (> 81%) displayed phase differences of less than 2 h, which was our sampling time resolution (Fig. 3A). Only 6–7% of those pairs had phase differences greater than two hours. To further assess the similarity of timing of expression rhythms within genes we calculated the correlation of expression between allelic pairs using z-score normalized expression. We observed strong correlation of expression among pairs with two rhythmic alleles (Fig. 3B), such that more than 80% of those pairs had a correlation of expression greater than 0.85. In contrast, allelic pairs with either none or only one rhythmic allele had significantly weaker correlations. For example, only 28–34% of pairs with only one rhythmic allele had correlations larger than 0.85. Taken together, we observed that although the majority of rhythmic allelic groups have mixed rhythmicity, the alleles that cycle display a high degree of rhythmic coherence. Thus, our data suggests that, despite the significant sequence and structural changes that occur within potato homologous chromosomes and alleles [28, 41, 42, 46], the timing of expression of rhythmic genes is highly conserved across alleles.

Fig. 3
figure 3

Similarities of expression among cyclic alleles. A Phase differences between cyclic alleles. Data contains all three datasets, leaf (short day and long day) and tuber (short day). B Correlation of z-score normalized gene expression between allelic pairs in leaf tissue. SD, short day; LD, long day. Rhythmicity was determined in the respective (SD or LD dataset). C Correlation of z-score normalized gene expression between allelic pairs across different tissues (Tissue), or stress/hormone treatments (Stress). Rhythmicity was determined in the LD dataset, since the Tissue and Stress experiments were performed under a 15 h light photoperiod. D The ratio of non-synonymous (Kn) to synonymous (Ks) substitution rates between allelic pairs. Rhythmicity was determined in the SD dataset. For B-D: The horizontal bar indicates the mean; Kruskal–Wallis multiple groups test (p-value < 0.0001) was followed by a post-hoc Wilcoxon signed-rank test with Bonferroni correction and letters indicate significant differences adjusted p-value < 0.0001

Full size image

To investigate whether the tight co-expression of rhythmic allelic pairs within single tissues corresponds to a wider degree of co-expression across tissues and stresses, we used the tissue and development expression dataset ("Tissue" dataset), and another publicly available dataset composed of eight hormone and abiotic stress treatments and their respective controls ("Stress" dataset). Interestingly, in all datasets co-expression was higher for fully rhythmic allelic pairs than for pairs in which one or fewer alleles cycled (Fig. 3C). For example, in these tissue and stress datasets ~ 70% of fully rhythmic pairs had correlations larger than 0.85, but only ~ 45% of mixed rhythmicity pairs did.

To estimate whether these highly co-expressed cyclic allelic pairs are under purifying selection, we examined the ratio of nonsynonymous (Kn) to synonymous (Ks) substitution rate within our syntenic allelic pairs. The distribution of Ks between these pairs in Atlantic was similar to that of other tetraploid potatoes independent of their rhythmicity [46] (Supplemental Fig. 3). We observed a smaller Kn/Ks ratio within fully rhythmic syntenic allelic pairs than within pairs with either none or only one rhythmic allele (Fig. 3D). This observation indicates that fully cyclic allelic pairs might maintain higher functional similarity in addition to higher expression similarity than non-cyclic genes.

Tuber phase is delayed with respect to the leaf under temperature cycles

Our high temporal resolution data sets enabled us to perform accurate phase determinations between leaves and tubers under short day conditions. Transcript expression peaked at all times throughout the day (Fig. 4A-C), although in leaves most rhythmic transcripts peaked at dawn or during the second half of the light period (after ZT 8), with fewer transcripts peaking in the late night (ZT18-22). In the tuber, the majority of cyclic genes peaked during the light period (ZT 3–9). Interestingly, of the 521 transcripts rhythmic in both leaves and tubers, 47% displayed at least a 2 h delay in tuber (Fig. 4D), with an average delay of 4.4 h. In contrast, only 19% of those rhythmic transcripts had an advanced phase equal or larger than 2 h in tubers with respect to leaves. Functional categories enriched among transcripts with delayed phase in tuber included circadian, photoperiod and light signaling, and proteasomal degradation (Supplemental Fig. 4). Differences in phase between above ground and below ground tissues have been reported in A. thaliana [47,48,49,50].

Fig. 4
figure 4

Changes in timing of expression across tissues and photoperiods. A Distribution of phases of gene expression under short day conditions in leaves. B Distribution of phases of gene expression under long day conditions in leaves. C Distribution of phases of gene expression under short day conditions in tubers. D Comparison of phase of expression of genes cycling in both leaves and tubers under short day conditions. E Expression of the GH3::LUC bioluminescence reporter (average ± SE, n = 6–16). Grey areas indicate dark/subjective dark periods. LD, light/dark conditions. Imaging was performed under temperature cycles, 18°C for dark/subjective dark and 22°C for light/subjective light periods

Full size image

Since our plants were grown under both light and temperature cycles, phase differences between underground tubers with respect to leaves could be due to delays of the entrainment signals coming from the shoot, delays in the perception of environmental temperature oscillations due to the temperature buffering capacity of the soil, or differences in entrainment phase in the soil due to the absence of direct light signals. We used the previously characterized expression reporter, GH3::LUC, to determine the phase of expression in leaves and tubers [15]. We observed that there were no significant differences in phase between leaves and attached or detached tubers when both tissues were exposed to the same light and temperature cycles (Fig. 4E). However, when tubers were kept in the dark, the phase of both attached and detached tubers was delayed with respect to leaves, even if both tissues were exposed to temperature cycles. These results indicate that light signals are able to entrain the tuber clock directly, and in their absence, there is a delay in phase in that tissue with respect to leaves. Under natural field conditions, temperature cycles are offset from light/dark cycles such that the minimum temperature generally coincides with dawn and the maximum temperature with the end of the day. This shift, together with temperature buffering by the soil, would suggest that tuber rhythms of field grown plants might be several hours phase shifted from leaf rhythms.

Photoperiod-dependent changes in gene expression in potato leaves

We next evaluated photoperiod-dependent changes in rhythmic gene expression in potato leaves. Among transcripts rhythmic in both SD and LD, we observed that a higher percentage had significantly elevated average expression under short days (~ 11%) than under long days (~ 5%) (Fig. 5A). We also observed that the average amplitude of gene expression was larger under SD than under LD (Fig. 5B). However, only ~ 23% of transcripts with larger amplitudes in SD had significantly elevated expression under short photoperiods (Supplemental Fig. 5 A). The larger amplitude could explain why we detected a higher number of rhythmic transcripts in SD (Fig. 1A,B). Transcripts with larger amplitudes under short days were expressed almost exclusively during the light in both photoperiods (Supplemental Fig. 5B) and, in comparison to other cyclic transcripts, were strongly enriched in photosynthesis-related functions (FDR-adjusted p-value < 0.001) (Supplemental Table S2). We identified several transcription regulators involved in photoperiod and stress responses among genes with larger SD amplitudes such as REVEILLE 2 (RVE2), CDF2 and HEAT SHOCK TRANSCTRIPTION FACTOR A6B (HSFA6B) (Supplemental Fig. 5 C).

Fig. 5
figure 5

Expression changes between short (SD) and long photoperiods (LD). A Volcano plot of differentially expressed transcripts between photoperiods. Blue indicates upregulated expression under SD and pink indicates upregulated expression under LD. B Pairwise amplitude (determined by JTK) of transcripts rhythmic in SD and LD. (****) Indicates adjusted p-value < 0.0001, Wilcoxon signed-rank test with Bonferroni correction. C Phase of expression in short and long days among genes cycling in both conditions in potato leaves. D Most informative functional enrichment terms of genes with changes in expression of six hours or greater under long days with respect to short days in potato leaves, when compared to all genes rhythmic in both SD and LD. The full set of enriched terms are in Supplemental Table S2. E Expression pattern of genes in the MapMan4 17.1 Ribosome biogenesis category with large changes in phase. Note that two TRB4/5a transcripts (Soltu.Atl_v3.03_0G008090 and Soltu.Atl_v3.03_0G008150) have the same sequence and therefore the same estimated expression. TRB4/5a-like transcript, Soltu.Atl_v3.03_0G008210 was not plotted, because it was very lowly expressed. Expression pattern of genes in the MapMan4 26.10 Pathogen and the 15.5.7.5 WRKY transcription factor categories with large changes in phase. Expression values were z-score normalized before plotting; the average of three biological replicates is shown. All times represent hours after dawn (ZT). F Summary motif for TRB-binding sites enriched in the Ribosome biogenesis category of panel C. G Summary motif for WRKY-binding sites enriched in genes in the Pathogen category of panel C

Full size image

As in other studies [51] we observed that many rhythmic genes had delayed phases under long days relative to short days (Fig. 5C). About 51% of genes that cycled in both photoperiods were delayed in phase by at least 2 h under long days, with an average delay of 3.7 h. We did not find any significant differences in functional enrichment in those transcripts when compared to all transcripts rhythmic in both conditions. In contrast, only about 12% of cycling genes had an advanced phase under long days.

Strong differences in phase of expression between photoperiods has been hypothesized to be a mechanism to optimize translation under lower energetic resources under short days in photosynthetic organisms [52]. We found 716 genes with phase delays under long days of at least 6 h. These genes were strongly enriched in functions related to ribosome biogenesis, pathogen responses, pyrimidine catabolism and WRKY transcription factor activity (Fig. 5D, Supplemental Table S2).

Protein biosynthesis genes counted for 28% of the 716 transcripts with at least a six hour phase delay under long days, the majority of which were related to ribosome biogenesis (Supplemental Table S2). The mRNA level of these genes peaked around dawn under short days and the second half of the light period under long days (Fig. 5E). Moreover, although not all were detected as exhibiting phase shifts, in our experiments the majority of transcripts in the functional categories of ribosome biogenesis-large ribosomal subunit (17.1.2) and ribosome biogenesis-small ribosomal subunit (17.1.3) followed similar expression patterns (Supplemental Fig. 6). Similar photoperiod dependent changes in expression have been observed in some ribosomal genes in A. thaliana [53]. In A. thaliana, TELOMERE BINDING PROTEINS (TRB) have been shown to bind to the promoters of ribosome biogenesis genes via telobox motifs [54]. We identified several TRB genes in potato (Supplemental Fig. 7), one of which, TRB4/5a, displayed maximum mRNA levels at the end of the day in long days and around dawn under short days (Fig. 5E). We also observed that TRB-related cis-elements were enriched in upstream regions of the ribosomal genes with strong changes in phase (Fig. 5F, Supplemental Table S3). Taken together, these photoperiod-controlled TRB proteins could be involved in the observed photoperiod dependent shift in expression.

Of the 76 genes related to pathogen responses that were rhythmic under both photoperiods, 41 displayed at least six hour delays in phase under long days. These genes included 20 genes potentially involved in effector-triggered immunity and 10 genes related to systemic acquired resistance (Supplemental Table S2). These pathogen response genes had peaks of expression at the beginning of the night period under short days and end of the night under long days (Fig. 5E). Interestingly, we found that WRKY transcription factors were also enriched within the strong phase shifted genes and had a very similar expression pattern as the pathogen related genes. Moreover, cis-element motifs recognized by WRKY transcription factors in A. thaliana were enriched in the upstream regions of the phase shifted pathogen related genes (Fig. 5G, Supplemental Table S3), indicating that the WRKY genes displaying changes in phase are involved in pathogen responses in potato.

Expression characteristics of clock, photoperiod and tuberization associated genes

Diel and circadian regulation of gene expression plays a key role in the regulation of growth and photoperiod responses in plants and are hypothesized to be involved in tuber formation in potato [55,56,57]. We identified 496 potential circadian, photoperiod or tuberization (CPT) transcripts in the Atlantic genome based on similarity to A. thaliana genes (Supplemental Table S1). Using both our computational pipeline and manual curation we linked the 496 CPT genes to 153 allelic groups in Atlantic, of which, 53% had four alleles, 25% had three alleles, and only 14% and 8% had two and one allele, respectively.

We classified our CPT genes into five categories based on their characterized function in A. thaliana, potato or tomato. The "Clock" set of genes included genes associated with core clock function such as the pseudo-response regulators (PRRs), or EARLY FLOWERING TIME 3. "Clock_Aux", included genes with pleiotropic functions that can affect circadian control, such as casein kinases or phytochrome interacting factors. The "Clock_Photoperiod" category includes genes involved in both clock and photoperiod control such as GIGANTEA (GI). The "Photoperiod" category includes genes mainly involved in photoperiod signaling such CO or the CDFs. Finally, the "Tuberization" category includes genes involved in tuber formation and growth, such as BEL-like transcription factors, gibberellin oxidases and genes involved in starch accumulation.

We observed that CPT genes in all categories displayed a higher percentage of rhythmicity and had higher rates of full rhythmicity than the genome-wide averages (Fig. 6A,B). Genes associated with core circadian clock function or photoperiodic regulation had a higher percentage of rhythmic transcripts than groups merely associated with the tuber formation process (Fig. 6A). Interestingly, the group of genes we named clock auxiliary had a lower number of rhythmic transcripts than the group of core clock genes in all conditions. Genes associated with core clock function were also more likely to have all their alleles rhythmic than genes in the other categories (Fig. 6B). Among the core clock associated allelic groups that were fully rhythmic in both leaf and tubers we found REVEILLE 3/5 (RVE3/5), RVE4/8, RVE1a, PRR5, PRR3a, PRR3b, LUX ARRHYTHMO/BROTHER OF LUX (LUX-BOAb) and LATE ELONGATED HYPOCOTYL (LHY) (Fig. 6C) [58].

Fig. 6
figure 6

Rhythmic expression of clock, photoperiod and tuberization related genes. Rhythmicity of transcripts (A) and genes (B) Rhythmicity was quantified using JTK. Fully rhythmic genes were genes with all their alleles cycling. Horizontal lines indicate the respective genome wide percentages in each condition. C Expression of fully rhythmic circadian clock related genes in short days. Values were z-score normalized before plotting; the averages of three biological replicates for each allele are shown. All times represent hours after dawn (ZT)

Full size image

Expression of photoperiod regulating genes in cultivated potato suggest the presence of a CO-independent mechanism of SP6A induction

Like other modern potato cultivars, Atlantic tuberizes under long photoperiods. In photoperiod sensitive potato (S. tuberosum Group Andigena), CO is a repressor of tuberization and its expression is repressed by CDF1 [55, 56]. The Atlantic genome contains four CDF1 alleles. Two alleles (haplotype 1 and haplotype 0) are identical and encode for truncated proteins missing the C-terminus required for interaction with the F-box protein FKF1 in the A. thaliana ortholog (Supplemental Fig. 8). This interaction has been hypothesized to mediate CDF1 degradation under long days. Therefore, truncated CDF1 proteins are likely to have enhanced protein stability and are hypothesized to mediate the repression of CO transcription even under long days [56]. We observed that CO1 and CO2, the highest expressed CO genes in Atlantic, displayed a phase delay and slightly lower expression under long days in spite of the presence of deregulated CDF1 alleles (Fig. 7). This phase difference has also been observed in other potatoes, both photoperiod sensitive and less sensitive species [55, 59]. In contrast, the lower-expressed CO3 had slightly higher expression under long-day conditions (Fig. 7).

Fig. 7
figure 7

Expression patterns under short and long photoperiods. Each panel represents the mean expression of one allele (average ± standard deviation, n = 3). SD, short days; LD, long days. For the CO genes, we also included the non-syntenic alleles CO1b (Soltu.Atl_v3.02_2G030560), CO2b (Soltu.Atl_v3.02_2G030580) and CO3b (Soltu.Atl_v3.02_2G030600)

Full size image

The average expression of Atlantic SP5G and SP6A genes was significantly sensitive to photoperiod (Fig. 7). COs are believed to activate SP5G expression, which in turn represses SP6A, the mobile tuber inducing factor [55, 60]. There are four alleles of SP5G annotated in the Atlantic genome, three of which were expressed in the leaf, at a low level. All four alleles were predominantly expressed in tubers in our experiments (Supplemental Fig. 9). In leaves, SP5G alleles were more than four-fold higher expressed under long days. The two SP6A alleles found in Atlantic were about 30-fold higher expressed under short days than under long days. However, SP6A expression was still detectable in long days in Atlantic leaves, which might be sufficient to drive tuberization in this cultivar under those conditions.

Discussion

Conserved highly expressed genes are more likely to display cyclic expression at the whole organ level

Polyploidy in plants represents a good model to investigate relative short-term selection on gene neofunctionalization and the importance of rhythmic expression for gene function. Here, we observed high similarity of expression between alleles under diel cycles. Ferrari et al. [61] demonstrated that the timing of orthologous gene expression is strongly conserved across taxa. We had previously shown cyclic gene expression to be predominantly retained after gene duplication [62]. We observed high conservation of cyclic expression within potato alleles, such that most cyclic allelic pairs displayed negligible differences in phase and high correlations gene expression under diel cycles (Fig. 3). These results indicate that the timing of gene expression is a key part of the function of cycling genes and is retained across different evolutionary scales.

Moreover, Ferrari et al. [61] had shown that genes conserved throughout the Archaeplastida were more likely to cycle than expected by chance. We observe that rhythmic allelic pairs have smaller Kn/Ks ratios and therefore are more likely to have the same function than non-cyclic allelic pairs, indicating higher functional conservation. Potato allelic pairs display highly similar patterns of expression across tissues and multiple stress conditions than non-cyclic allelic pairs (Fig. 3). Taken together these results support the hypothesis that cyclic gene expression occurs more often in conserved genes playing core roles and that there is less neofunctionalization of cyclic genes than non-cyclic genes.

We observed that rhythmicity was associated with higher expression across genes and also within alleles (Fig. 1C, Supplemental Fig. 2B). Similarly, cyclic paralogs in B. rapa [9] and cyclic homeologs in hexaploid wheat [10] are higher expressed than non-cyclic ones. These observations can be partly explained by the bias of rhythmic detection algorithms toward higher expressed genes that display low noise to signal ratio [45] and we observed that non-rhythmic alleles had higher noise than rhythmic ones (Supplemental Fig. 2B). Differences in phase and/or transcript content across cell types [63] can also generate low and noisy expression patterns at the organ level and will bias the types of transcripts considered rhythmic. However, our analyses comparing only well-expressed and low-noise genes (Supplemental Fig. 2D) indicated that cyclic genes might be on average higher expressed. This is a similar conclusion reached by a metanalysis of several animal species evaluated with different rhythmic signal detection algorithms [45]. In other organisms cyclic mRNA levels have been associated with energetically "expensive" genes, i.e. genes with high costs of transcription and translation, mainly because cyclic genes have higher expression [63]. The authors of these studies, hypothesized that cyclic expression enables the generation of high gene expression at time of need while minimizing the overall cost of protein production [64]. Interestingly, in plants there might also be a correlation between stronger rhythms of higher amplitude with higher growth and metabolic rate. For example, in both gymnosperms and angiosperm trees diel oscillations only occur in the summer and not in the winter [65, 66]; young A. thaliana plants have a greater number of cyclic genes and higher amplitude rhythms than in older plants [67], and pathogen responses that restrict growth are also associated with a decrease in amplitude [68]. Unicellular eukaryotic Archaeplastida also have higher amplitude rhythms of gene expression than multicellular ones, which may be linked to a higher growth rate [61]. In contrast, we observed a higher number of cyclic transcripts and an increase in amplitude under short photoperiods in cultivated potato (Figs. 1 and 5). This increase in amplitude was mainly due to changes in photosynthesis related genes. In Camelina sativa, although biomass accumulation is decreased under short days, photosynthesis rates were increased, which has been hypothesized to compensate for the shortened light period [69]. Higher amplitude rhythms could enable this rate increase without the additional cost of producing more RNA. Further studies are needed to test the role of metabolic activity on rhythmic gene expression in photosynthetic organisms.

Some genes display especially strong changes in the phase of expression between photoperiods

It has been proposed that genes expressed at dawn display a decrease in their respective protein concentration as the photoperiod becomes longer due to the increases in vegetative growth rates in long days [52]. In contrast, evening expressed genes maintain or increase protein abundance with longer photoperiods. We observed that most ribosome biogenesis genes in cultivated potato displayed a strong difference in phase between photoperiods, such that mRNA content peaked around dawn under short days and at the end of the day under long days. In A. thaliana, ribosome associated genes have similar photoperiod-dependent changes in expression patterns [53] and appear to be preferentially translated at the end of the day under long photoperiods [70], indicating that, at least in long days, the phase of mRNA correlates with translation rate. The shift in phase of expression of ribosomal genes could represent a mechanism to maintain protein content and therefore plant growth under long days in cultivated potato and other plants. The mechanism by which these changes in gene expression are regulated is yet unclear. As in A. thaliana, we observe several TRB related genes that display similar changes in phase. TRB genes have been associated with ribosomal gene transcription in Arabidopsis [54, 71], but it remains to be tested whether they are also involved in the photoperiod specific differences observed in potato.

In our experiments, more than half of the transcripts in the functional category "external stimuli response-pathogen responses" that displayed cyclic expression had very strong differences in phase, such that their mRNA peaked early in the night under short days but at the end of the night under long days. Interestingly these genes were enriched in potential effector-triggered-immunity and systemic acquired resistance function. Light signals and the circadian clock have been shown to modulate time-dependent changes in sensitivity to pathogens in A. thaliana [72,73,74,75]. In cultivated potato, the upstream sequences of these photoperiod regulated pathogen-response genes contained cis-elements that potentially mediate regulation via WRKY transcription factors. We identified three WRKY factors tightly co-expressed with these pathogen-related genes that are strong candidates for controlling these responses. They belong to three genes that previous studies in S. tuberosum Group Phureja named StWRKY30, StWRKY38, StWRKY63 [76]. Potato WRKY transcription factors have been associated with stress responses including pathogen responses [76, 77]. However, there is currently no direct evidence for the involvement of these three WRKYs in pathogen responses in potato.

Photoperiod control of tuberization in cultivated potato

Studies using different potato genotypes indicate that there are still unanswered questions on the CDF1-CO-SP5G-SP6A model of the photoperiod control of tuberization. Truncated alleles of CDF1 are the principal mechanism for strongly reduced photoperiod sensitivity of tuberization in cultivated potatoes [5, 12, 56]. These truncated alleles lead to elevated levels of CDF1 protein and have been proposed to constitutively repress CO genes, leading to photoperiod independent tuberization [56]. Potato CO proteins, like A. thaliana CO, are stabilized by light, and photoperiod-dependent differences in mRNA peak levels have been linked to large changes in CO protein between short and long days in photoperiod sensitive potato [55, 78, 79]. However, the characteristic phase difference of CO genes between 8-h light (maximum during the dark period) and 16-h light (maximum early in the day) photoperiods is maintained in both photoperiod sensitive potatoes such as S. tuberosum spp Andigena and less photoperiod sensitive potatoes such as S. neotuberosum [59]. We observed similar shifts in the peak of the expression in Atlantic, which has weak photoperiod sensitivity, when comparing12-hour light and 16-h light conditions (Fig. 7). S. neotuberosum and Atlantic contain endogenous CDF1 truncated alleles [59], and therefore, more detailed analyses are needed to understand how truncated CDF1 alleles influence the expression patterns of CO genes in different potatoes.

There are differences in expression of SP5G and SP6A between different potato genotypes. In Atlantic, we detected photoperiod dependent changes in both SP5G and SP6A expression, although only SP5G was rhythmic (Fig. 7). In contrast, S. neotuberosum displays constitutive low levels of SP5G but photoperiod dependent changes in SP6A mRNA content [59]. A lack of apparent correlation between SP5G and SP6A expression has also been observed in transgenic S. tuberosum spp. Andigena plants overexpressing CO1. These transgenic lines have similar overall levels of SP5G under both photoperiods early in the day but only show an increase in SP6A under short days and accordingly still display photoperiod sensitivity of tuberization [78]. These results suggests that there other levels of SP6G regulation in potato leaves. BEL5, encoded by a phloem mobile mRNA, is also able to induce SP6A expression in potato leaves [80], and in Atlantic leaves BEL5 expression is induced about two-fold under short days (Supplemental Fig. 10). It has been proposed [56] that CDF1 could directly regulate SP6A in potato, however, SP6A has not been identified as a direct CDF1 target in recent work on cultivated potato [81]. Therefore, further studies are needed to clarify how photoperiod signals regulate SP6A expression to trigger tuberization.

Conclusion

Our analysis of diel expression across haplotypes in cultivated tetraploid potato showed that rhythmic alleles not only retain high degree of co-expression under diel cycles but also across tissues, developmental stages and stress conditions, when compared to non-cyclic alleles. Moreover, smaller Kn/Ks ratios also suggest that rhythmic allelic pairs are likely to retain higher functional similarity. In general, rhythmic transcripts appeared higher and more constitutively expressed than non-rhythmic transcripts. Taken together these results suggest that rhythmic expression plays a key role in core plant cell functions. Finally, our results also showed that genes related to core circadian and photoperiod sensing functions are more likely to be rhythmic than genes involved in tuber formation. Moreover, photoperiod-dependent changes in expression of genes associated with photoperiod sensing reveal open questions about the regulation of tuberization in cultivated potato.

Data availability

The diel expression data can be accessed under the Sequence Read Archive BioProjects PRJNA957457 (short day) and PRJNA1093480 (long day). Atlantic Developmental Gene Expression Atlas expression data were obtained from NCBI under BioProject PRJNA753086. The code used in this study is available at GitHub [82]. The large datasets (Dataset S1-S14) are available here [83].

References

  1. Salman-Minkov A, Sabath N, Mayrose I. Whole-genome duplication as a key factor in crop domestication. Nat Plants. 2016;2: 16115.

    Article  CAS  PubMed  Google Scholar 

  2. Muller NA, Zhang L, Koornneef M, Jimenez-Gomez JM. Mutations in EID1 and LNK2 caused light-conditional clock deceleration during tomato domestication. Proc Natl Acad Sci U S A. 2018;115:7135–40.

    Article  PubMed  PubMed Central  Google Scholar 

  3. Steed G, Ramirez DC, Hannah MA, Webb AAR. Chronoculture, harnessing the circadian clock to improve crop yield and sustainability. Science. 2021;372:eabc9141.

    Article  CAS  PubMed  Google Scholar 

  4. McClung CR. Circadian clock components offer targets for crop domestication and improvement. Genes. 2021;12(3): 12:374.

    Article  PubMed  Google Scholar 

  5. Hardigan MA, Laimbeer FPE, Newton L, Crisovan E, Hamilton JP, Vaillancourt B, et al. Genome diversity of tuber-bearing Solanum uncovers complex evolutionary history and targets of domestication in the cultivated potato. Proc Natl Acad Sci U S A. 2017;114: E9999–10008.

    Article  Google Scholar 

  6. Wang X, Wang H, Wang J, Sun R, Wu J, Liu S, et al. The genome of the mesopolyploid crop species Brassica rapa. Nat Genet. 2011;43:1035–9.

    Article  CAS  PubMed  Google Scholar 

  7. Cheng F, Mandáková T, Wu J, Xie Q, Lysak MA, Wang X. Deciphering the diploid ancestral genome of the mesohexaploid Brassica rapa. Plant Cell. 2013;25:1541–54.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Levy AA, Feldman M. Evolution and origin of bread wheat. Plant Cell. 2022;34:2549–67.

    Article  PubMed  PubMed Central  Google Scholar 

  9. Greenham K, Sartor RC, Zorich S, Lou P, Mockler TC, McClung CR. Expansion of the circadian transcriptome in Brassica rapa and genome-wide diversification of paralog expression patterns. Elife. 2020;9: e58993.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Rees H, Rusholme-Pilcher R, Bailey P, Colmer J, White B, Reynolds C, et al. Circadian regulation of the transcriptome in a complex polyploid crop. PLoS Biol. 2022;20: e3001802.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Jansky SH, Spooner DM. The Evolution of Potato Breeding. In: Goldman I, editor. Plant Breeding Reviews. 1st ed. Wiley; 2018. p. 169–214.

    Chapter  Google Scholar 

  12. Gutaker RM, Weiß CL, Ellis D, Anglin NL, Knapp S, Luis Fernández-Alonso J, et al. The origins and adaptation of European potatoes reconstructed from historical genomes. Nat Ecol Evol. 2019;3:1093–101.

    Article  PubMed  Google Scholar 

  13. Zhang C, Wang P, Tang D, Yang Z, Lu F, Qi J, et al. The genetic basis of inbreeding depression in potato. Nat Genet. 2019;51:374–8.

    Article  CAS  PubMed  Google Scholar 

  14. Hijmans RJ. Global distribution of the potato crop. Am J Potato Res. 2001;78:403–12.

    Article  Google Scholar 

  15. Hoopes GM, Zarka D, Feke A, Acheson K, Hamilton JP, Douches D, et al. Keeping time in the dark: potato diel and circadian rhythmic gene expression reveals tissue-specific circadian clocks. Plant Direct. 2022;6: e425.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Pham GM, Hamilton JP, Wood JC, Burke JT, Zhao H, Vaillancourt B, et al. Construction of achromosome-scale long-read reference genome assembly for potato. Gigascience. 2020;9(9):giaa100. https://doi.org/10.1093/gigascience/giaa100.

  17. Spud DB Potato Genomics Resource. https://spuddb.uga.edu/. Accessed 1 May 2025.

  18. Lovell JT, Sreedasyam A, Schranz ME, Wilson M, Carlson JW, Harkess A, et al. Genespace tracks regions of interest and gene copy number variation across multiple genomes. Elife. 2022;11:e78526.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Hamilton JP, Brose J, Buell CR. SpudDB: a database for accessing potato genomic data. Genetics. 2025;229:iyae205.

    Article  CAS  PubMed  Google Scholar 

  20. Wan CY, Wilkins TA. A modified hot borate method significantly enhances the yield of high-quality RNA from cotton (Gossypium hirsutum L.). Anal Biochem. 1994;223:7–12.

    Article  CAS  PubMed  Google Scholar 

  21. Wood DE, Lu J, Langmead B. Improved metagenomic analysis with kraken 2. Genome Biol. 2019;20:257.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Kraken 2. https://benlangmead.github.io/aws-indexes/k2. Accessed 1 May 2025.

  23. Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J. 2011;17:10.

    Article  Google Scholar 

  24. Babraham Bioinformatics. https://www.bioinformatics.babraham.ac.uk/projects/fastqc/. Accessed 1 May 2025.

  25. Ewels P, Magnusson M, Lundin S, Käller M. Multiqc: summarize analysis results for multiple tools and samples in a single report. Bioinformatics. 2016;32:3047–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Seqtk. https://github.com/lh3/seqtk. Accessed 1 May 2025.

  27. Bray NL, Pimentel H, Melsted P, Pachter L. Near-optimal probabilistic RNA-seq quantification. Nat Biotechnol. 2016;34:525–7.

    Article  CAS  PubMed  Google Scholar 

  28. Hoopes G, Meng X, Hamilton JP, Achakkagari SR, de Alves Freitas Guesdes F, Bolger ME, et al. Phased, chromosome-scale genome assemblies of tetraploid potato reveal a complex genome, transcriptome, and predicted proteome landscape underpinning genetic diversity. Mol Plant. 2022;15:520–36.

    Article  CAS  PubMed  Google Scholar 

  29. Brose J,Martin D, Wang Y-W, Wood JC, Vaillancourt B, Hamilton JP, et al. An allelic resolution gene atlasfor tetraploid potato provides insights into tuberization and stress resilience. bioRxiv [Preprint].2025. https://doi.org/10.1101/2025.06.26.661617.

  30. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.

    Article  PubMed  PubMed Central  Google Scholar 

  31. Wu G, Anafi RC, Hughes ME, Kornacker K, Hogenesch JB. MetaCycle: an integrated R package to evaluate periodicity in large scale data. Bioinformatics. 2016;32:3351–3.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Virtanen P, Gommers R, Oliphant TE, Haberland M, Reddy T, Cournapeau D, et al. Scipy 1.0: fundamental algorithms for scientific computing in Python. Nat Methods. 2020;17:261–72.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. McKinney W. Data Structures for Statistical Computing in Python. In: Walt S van der, Millman J, editors. Proceedings of the 9th Python in Science Conference. 2010. p. 56–61.

  34. Kryuchkova-Mostacci N, Robinson-Rechavi M. A benchmark of gene expression tissue-specificity metrics. Brief Bioinform. 2017;18:205–14.

    CAS  PubMed  Google Scholar 

  35. Genome Evolution - CoGe. https://genomevolution.org/coge/. Accessed 1 May 2025.

  36. Albert VA, Krabbenhoft TJ. Navigating the CoGe Online Software Suite for Polyploidy Research. In: Van de Peer Y, editor. Polyploidy: Methods and Protocols. New York, NY: Springer US; 2023. p. 19–45.

  37. Bolger M, Schwacke R, Usadel B. Mapman visualization of RNA-Seq data using mercator4 functional annotations. Methods Mol Biol. 2021;2354:195–212.

    Article  PubMed  Google Scholar 

  38. Mercator 4. www.plabipd.de/mercator_main.html. Accessed 1 May 2025.

  39. Heinz S, Benner C, Spann N, Bertolino E, Lin YC, Laslo P, et al. Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities. Mol Cell. 2010;38:576–89.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Jin J, Tian F, Yang DC, Meng YQ, Kong L, Luo J, et al. PlantTFDB 4.0: toward a central hub for transcription factors and regulatory interactions in plants. Nucleic Acids Res. 2017;45:D1040-5.

    Article  CAS  PubMed  Google Scholar 

  41. Sun H, Jiao WB, Krause K, Campoy JA, Goel M, Folz-Donahue K, et al. Chromosome-scale and haplotype-resolved genome assembly of a tetraploid potato cultivar. Nat Genet. 2022;54:342–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Bao Z, Li C, Li G, Wang P, Peng Z, Cheng L, et al. Genome architecture and tetrasomic inheritance of autotetraploid potato. Mol Plant. 2022;15:1211–26.

    Article  CAS  PubMed  Google Scholar 

  43. Ewing EE, Struik PC. Tuber Formation in Potato: Induction, Initiation, and Growth. In: Janick J, editor. Horticultural Reviews. 1st ed. Wiley; 1992. p. 89–198.

    Chapter  Google Scholar 

  44. Hughes ME, Hogenesch JB, Kornacker K. JTK_cycle: an efficient nonparametric algorithm for detecting rhythmic components in genome-scale data sets. J Biol Rhythms. 2010;25:372–80.

    Article  PubMed  PubMed Central  Google Scholar 

  45. Laloum D, Robinson-Rechavi M. Methods detecting rhythmic gene expression are biologically relevant only for strong signal. PLoS Comput Biol. 2020;16:e1007666.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Wang F, Xia Z, Zou M, Zhao L, Jiang S, Zhou Y, et al. The autotetraploid potato genome provides insights into highly heterozygous species. Plant Biotechnol J. 2022;20:1996–2005.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Endo M. Tissue-specific circadian clocks in plants. Curr Opin Plant Biol. 2016;29:44–9.

    Article  PubMed  Google Scholar 

  48. James AB, Monreal JA, Nimmo GA, Kelly CL, Herzyk P, Jenkins GI, et al. The circadian clock in Arabidopsis roots is a simplified slave version of the clock in shoots. Science. 2008;322:1832–5.

    Article  CAS  PubMed  Google Scholar 

  49. Takahashi N, Hirata Y, Aihara K, Mas P. A hierarchical multi-oscillator network orchestrates the Arabidopsis circadian system. Cell. 2015;163:148–59.

    Article  CAS  PubMed  Google Scholar 

  50. Bordage S, Sullivan S, Laird J, Millar AJ, Nimmo HG. Organ specificity in the plant circadian system is explained by different light inputs to the shoot and root clocks. New Phytol. 2016;212:136–49.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  51. Michael TP, Mockler TC, Breton G, McEntee C, Byer A, Trout JD, et al. Network discovery pipeline elucidates conserved time-of-day-specific cis-regulatory modules. PLoS Genet. 2008;4:e14.

    Article  PubMed  PubMed Central  Google Scholar 

  52. Seaton DD, Graf A, Baerenfaller K, Stitt M, Millar AJ, Gruissem W. Photoperiodic control of the Arabidopsis proteome reveals a translational coincidence mechanism. Mol Syst Biol. 2018;14:e7962.

    Article  PubMed  PubMed Central  Google Scholar 

  53. Wang Q, Liu W, Leung CC, Tarté DA, Gendron JM. Plants distinguish different photoperiods to independently control seasonal flowering and growth. Science. 2024;383:eadg9196.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Schrumpfová PP, Vychodilová I, Hapala J, Schořová Š, Dvořáček V, Fajkus J. Telomere binding protein TRB1 is associated with promoters of translation machinery genes in vivo. Plant Mol Biol. 2016;90:189–206.

    Article  PubMed  Google Scholar 

  55. Abelenda JA, Cruz-Oro E, Franco-Zorrilla JM, Prat S. Potato StCONSTANS-like1 suppresses storage organ formation by directly activating the FT-like StSP5G repressor. Curr Biol. 2016;26:872–81.

    Article  CAS  PubMed  Google Scholar 

  56. Kloosterman B, Abelenda JA, Gomez Mdel M, Oortwijn M, de Boer JM, Kowitwanich K, et al. Naturally occurring allele diversity allows potato cultivation in northern latitudes. Nature. 2013;495:246–50.

    Article  CAS  PubMed  Google Scholar 

  57. Abelenda JA, Navarro C, Prat S. Flowering and tuberization: a tale of two nightshades. Trends Plant Sci. 2014;19:115–22.

    Article  CAS  PubMed  Google Scholar 

  58. Harmer SL. The time machine: feedback loops, post-transcriptional regulation, and environmental integration in the plant circadian oscillator. Plant J. 2025;122: e70275.

    Article  CAS  PubMed  Google Scholar 

  59. Morris WL, Hancock RD, Ducreux LJ, Morris JA, Usman M, Verrall SR, et al. Day length dependent restructuring of the leaf transcriptome and metabolome in potato genotypes with contrasting tuberization phenotypes. Plant Cell Environ. 2014;37:1351–63.

    Article  CAS  PubMed  Google Scholar 

  60. Navarro C, Abelenda JA, Cruz-Oro E, Cuellar CA, Tamaki S, Silva J, et al. Control of flowering and storage organ formation in potato by flowering locus T. Nature. 2011;478:119–22.

    Article  CAS  PubMed  Google Scholar 

  61. Ferrari C, Proost S, Janowski M, Becker J, Nikoloski Z, Bhattacharya D, et al. Kingdom-wide comparison reveals the evolution of diurnal gene expression in Archaeplastida. Nat Commun. 2019. https://doi.org/10.1038/s41467-019-08703-2.

    Article  PubMed  PubMed Central  Google Scholar 

  62. Panchy N, Wu G, Newton L, Tsai CH, Chen J, Benning C, et al. Prevalence, evolution, and cis-regulation of diel transcription in Chlamydomonas reinhardtii. G3 Genes|Genomes|Genetics. 2014;4:2461–71.

    Article  PubMed  PubMed Central  Google Scholar 

  63. Qin Y, Liu Z, Gao S, Martínez-Vasallo C, Long Y, Zhu X, et al. 48-hour and 24-hour time-lapse single-nucleus transcriptomics reveal cell-type specific circadian rhythms in Arabidopsis. Nat Commun. 2025;16:4171.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  64. Wang G-Z, Hickey SL, Shi L, Huang H-C, Nakashe P, Koike N, et al. Cycling transcriptional networks optimize energy utilization on a genome scale. Cell Rep. 2015;13:1868–80.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  65. Ramos A, Perez-Solis E, Ibanez C, Casado R, Collada C, Gomez L, et al. Winter disruption of the circadian clock in chestnut. Proc Natl Acad Sci U S A. 2005;102:7037–42.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  66. Nose M, Watanabe A. Clock genes and diurnal transcriptome dynamics in summer and winter in the gymnosperm Japanese cedar (Cryptomeria japonica (L.f.) D.Don). BMC Plant Biol. 2014;14: 308.

    PubMed  PubMed Central  Google Scholar 

  67. Jung S, Kim H, Lee J, Kang MH, Kim J, Kim JK, et al. The genetically programmed rhythmic alteration of diurnal gene expression in the aged Arabidopsis leaves. Front Plant Sci. 2024;15:1481682.

    Article  PubMed  PubMed Central  Google Scholar 

  68. Li Z, Bonaldi K, Uribe F, Pruneda-Paz JL. A localized Pseudomonas syringae infection triggers systemic clock responses in Arabidopsis. Curr Biol. 2018;28:630-639.e4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  69. Xu Y, Koroma AA, Weise SE, Fu X, Sharkey TD, Shachar-Hill Y. Daylength variation affects growth, photosynthesis, leaf metabolism, partitioning, and metabolic fluxes. Plant Physiol. 2023;194:475–90.

    Article  PubMed  PubMed Central  Google Scholar 

  70. Missra A, Ernest B, Lohoff T, Jia Q, Satterlee J, Ke K, et al. The circadian clock modulates global daily cycles of mRNA ribosome loading. Plant Cell. 2015;27:2582–99.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  71. Amiard S, Feit L, Vanrobays E, Simon L, Le Goff S, Loizeau L, et al. The telomere repeat binding proteins TRB4 and TRB5 function as transcriptional activators of PRC2-controlled genes to regulate plant development. Plant Commun. 2024;5: 100890.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  72. Griebel T, Zeier J. Light regulation and daytime dependency of inducible plant defenses in Arabidopsis: phytochrome signaling controls systemic acquired resistance rather than local defense. Plant Physiol. 2008;147:790–801.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  73. Bhardwaj V, Meier S, Petersen LN, Ingle RA, Roden LC. Defence responses of Arabidopsis thaliana to infection by Pseudomonas syringae are regulated by the circadian clock. PLoS One. 2011;6:e26968.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  74. Cortleven A, Roeber VM, Frank M, Bertels J, Lortzing V, Beemster GTS, et al. Photoperiod stress in Arabidopsis thaliana induces a transcriptional response resembling that of pathogen infection. Front Plant Sci. 2022;13:838284.

    Article  PubMed  PubMed Central  Google Scholar 

  75. Ingle RA, Stoker C, Stone W, Adams N, Smith R, Grant M, et al. Jasmonate signalling drives time-of-day differences in susceptibility of Arabidopsis to the fungal pathogen Botrytis cinerea. Plant J. 2015;84:937–48.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  76. Zhang C, Wang D, Yang C, Kong N, Shi Z, Zhao P, et al. Genome-wide identification of the potato WRKY transcription factor family. PLoS One. 2017;12: e0181573.

    Article  PubMed  PubMed Central  Google Scholar 

  77. Villano C, Esposito S, D’Amelia V, Garramone R, Alioto D, Zoina A, et al. WRKY genes family study reveals tissue-specific and stress-responsive TFs in wild potato species. Sci Rep. 2020;10:7196.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  78. Plantenga FDM, Heuvelink E, Rienstra JA, Visser RGF, Bachem CWB, Marcelis LFM. Coincidence of potato CONSTANS (StCOL1) expression and light cannot explain night-break repression of tuberization. Physiol Plant. 2019;167:250–63.

    Article  CAS  PubMed  Google Scholar 

  79. González-Schain ND, Díaz-Mendoza M, Zurczak M, Suárez-López P. Potato CONSTANS is involved in photoperiodic tuberization in a graft-transmissible manner. Plant J. 2012;70:678–90.

    Article  PubMed  Google Scholar 

  80. Sharma P, Lin T, Hannapel DJ. Targets of the StBEL5 transcription factor include the FT ortholog StSP6A. Plant Physiol. 2016;170:310–24.

    Article  CAS  PubMed  Google Scholar 

  81. Shaikh MA, Ramírez-Gonzales L, Franco-Zorrilla JM, Steiner E, Oortwijn M, Bachem CWB, et al. StCDF1: A “jack of all trades” clock output with a central role in regulating potato nitrate reduction activity. New Phytol. 2025;245:282–98.

    Article  CAS  PubMed  Google Scholar 

  82. Autotetraploid potato diel transcriptome code in GitHub. https://github.com/efarre/Autotetraploid_potato_diel_transcriptome. Accessed 1 Aug 2025.

  83. Datasets for this manuscript. https://doi.org/10.5061/dryad.x69p8czwp. Accessed 1 May 2025.

  84. Boxshade. 2025. https://junli.netlify.app/apps/boxshade/.

Download references

Acknowledgements

We thank Marisa Laitinen for help harvesting plants and tissue culture. We also thank Patrick Edger for feedback on the manuscript and advise on the determination of synonymous and nonsynonymous substitution rates.

Large datasets

(https://doi.org/10.5061/dryad.x69p8czwp).

Dataset S1. Allelic groups.

Dataset S2. Diel expression rlog.

Dataset S3. Expression of Tissue samples from the Developmental Gene Expression Atlas rlog.

Dataset S4. Expression of Stress samples from the Developmental Gene Expression Atlas rlog.

Dataset S5. Leaf short day tissue cycling parameters as determined per JTK.

Dataset S6. Leaf long day cycling parameters as determined per JTK.

Dataset S7. Tuber short day cycling parameters as determined per JTK.

Dataset S8. Leaf pairwise allelic expression correlations in short days.

Dataset S9. Leaf pairwise allelic expression correlations in long days.

Dataset S10. Pairwise allelic expression correlations in Tissue samples from the Developmental Gene Expression Atlas.

Dataset S11. Pairwise allelic expression correlations in Stress samples from the Developmental Gene Expression Atlas.

Dataset S12. Differential expression short vs. long day determined by DEseq.

Dataset S13. Differential expression leaf vs. tuber under short days determined by DEseq.

Dataset S14. Functional annotation of Atlantic using Mercator4.

Funding

This project was funded by a National Science Foundation Plant Genome Research Project award to C.R.B and E.M.F. (IOS-1950376) and awards to C.R.B. (IOS-1929982, IOS-2140176).

Author information

Authors and Affiliations

  1. Department of Plant Biology, Michigan State University, East Lansing, MI, 48824, USA

    Ann Feke, Brieanne Vaillancourt, Kaitlyn Acheson, Julia Brose, John P. Hamilton, Joshua C. Wood, C. Robin Buell & Eva M. Farré

  2. Center for Applied Genetic Technologies, University of Georgia, Athens, GA, 30602, USA

    Brieanne Vaillancourt, Julia Brose, John P. Hamilton, Yi-Wen Wang, Joshua C. Wood & C. Robin Buell

  3. Department of Crop & Soil Sciences, University of Georgia, Athens, GA, 30602, USA

    John P. Hamilton & C. Robin Buell

  4. Department Genetics, University of Georgia, Athens, GA, 30602, USA

    Dionne Martin

  5. Institute of Plant Breeding, Genetics & Genomics, University of Georgia, Athens, GA, 30602, USA

    C. Robin Buell

  6. The Plant Center, University of Georgia, Athens, GA, 30602, USA

    C. Robin Buell

Authors
  1. Ann Feke
    View author publications

    Search author on:PubMed Google Scholar

  2. Brieanne Vaillancourt
    View author publications

    Search author on:PubMed Google Scholar

  3. Kaitlyn Acheson
    View author publications

    Search author on:PubMed Google Scholar

  4. Julia Brose
    View author publications

    Search author on:PubMed Google Scholar

  5. John P. Hamilton
    View author publications

    Search author on:PubMed Google Scholar

  6. Dionne Martin
    View author publications

    Search author on:PubMed Google Scholar

  7. Yi-Wen Wang
    View author publications

    Search author on:PubMed Google Scholar

  8. Joshua C. Wood
    View author publications

    Search author on:PubMed Google Scholar

  9. C. Robin Buell
    View author publications

    Search author on:PubMed Google Scholar

  10. Eva M. Farré
    View author publications

    Search author on:PubMed Google Scholar

Contributions

C.R.B and E.M.F conceived the research. A.F, B.V., K.A., J.W, C.R.B, E.M.F, designed and carried out the experiments. A.F., B.V, J.H., J.C.W, D.M, Y.W.W. and E.M.F performed the data analyses. A.F and E.M.F wrote the manuscript. All authors revised and edited the manuscript.

Corresponding author

Correspondence to Eva M. Farré.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

12864_2025_11945_MOESM1_ESM.pdf

Supplementary Material 1: Supplemental Figure 1. Expression of circadian clock related genes in short days. Values were z-score normalized before plotting; the averages of three biological replicates for each allele are shown. All times represent hours after dawn (ZT). Supplemental Figure 2. Allele specific cycling in potato tissues. A. Rhythmic strength, as determined by JTK-calculated FDR-adjusted p-value of cycling transcripts. The horizontal line indicates the mean. B. Comparison of the average coefficient of variation (CV) in each time point between rhythmic and non-rhythmic allelic pairs. C. Cycling strength, as determined by JTK calculated FDR adjusted mean p-value of cyclic transcripts in genes with either full or mixed rhythmicity. Only transcripts with average expression rlog >5 and average coefficient of variation between 0 and 0.1 were used. D. Average allele expression of only highly expressed transcripts (rlog >5) with average coefficient of variation between 0 and 0.1, within genes with no rhythmic alleles (Non), at least one rhythmic allele (Mixed), all alleles rhythmic (Full). The horizontal line indicates the mean. E. Comparison of the average expression within mixed-rhythmicity allelic pairs.“Yes” indicates the allele which cycles, while “No” indicates the allele that does not. Only transcripts with average expression rlog >5 and average coefficient of variation between 0 and 0.1 were used. F. Tau index determining tissue-specific expression of Atlantic transcripts either rhythmic (“Yes”) or non-rhythmic (No) in the different conditions; Tau ranges between 0 (constitutively expressed) to 1 (tissue specific). The horizontal line indicates the mean. For A and D, Kruskal-Wallis multiple groups test was performed before pairwise analyses. All pairwise comparisons used a Wilcoxon signed-rank tests with Bonferroni correction and in all instances, p-value < 0.0001 (****),<0.001 (***), < 0.01(**), <0.05 (*) and ns (not significant). SD, short days; LD, long days. Supplemental Figure 3. Distribution of Ks of syntenic allelic pairs of Solanum tuberosum cv. Atlantic, with either 0, 1 or 2 rhythmic alleles as determined by JTK using the short-day diel dataset. Supplemental Figure 4. Functional enrichment of rhythmic transcripts that had a delayed phase in tubers with respect to leaves. Transcripts with at least 2 h delayed phase were compared to all tuber expressed transcripts. Supplemental Figure 5. Expression characteristics of transcripts with larger amplitudes under short days. A Volcano plot of transcripts with larger amplitudes under SD than LD (Δ Amp>0.5). Blue indicates upregulated expression under SD and pink upregulated expression under LD. B. Phase of expression of genes with amplitudes higher in SD than LD (Δ Amp>0.5). C. Average leaf expression of transcriptional regulators with larger amplitudes under SD than LD (Δ Amp>1) (average n = 3). Supplemental Figure 6. Expression patterns of all genes in the Mercator4 bin categories related to Protein biosynthesis.ribosome biogenesis.large ribosomal subunit (LSU) 17.1.2 and Protein biosynthesis.ribosome biogenesis.small ribosomal subunit (SSU) 17.1.3. Supplemental Figure 7. TRB proteins in A. thaliana and cultivated potato. Peptide sequences were aligned using clustalW. Shading was implemented using BoxShade [84]. Supplemental Figure 8. CDF1 alleles in Atlantic. A. Expression patterns CDF1 alleles in short and long day conditions. Values represent the average ± standard deviation (n = 3). B. Alignment of the C-terminus end of CDF1 alleles. Supplemental Figure 9. Expression patterns in tubers under short photoperiods. Each panel represents the mean expression of one allele (average ± standard deviation, n = 3). CO genes include the non-syntenic alleles CO1b (Soltu.Atl_v3.02_2G030560),CO2b (Soltu.Atl_v3.02_2G030580), CO3b (Soltu.Atl_v3.02_2G030600). Supplemental Figure 10. Expression patterns BEL5 alleles in short and long day conditions in leaf tissue. Values represent the average ± standard deviation (n = 3).

12864_2025_11945_MOESM2_ESM.xlsx

Supplementary Material 2: Table S1. Clock, tuberization, photoperiod genes (CPT). A. CPT_ATL_DM-syntenic associations between Atlantic and DM CPT genes."Gene"indicates allelic group identified using the automatic pipeline (see Materials and Methods)."Gene_manual"contains additional associations after manual evaluation. H1-H4, indicate the haplotype, and H0 indicates gene models on the unphased assembly. B. CPT_C88_blast, five best blast hits of Atlantic CPT protein models against the Cooperation-88 (C88) protein models. C. CPT_Otava_blast, five best blast hits of Atlantic CPT protein models against the Otava protein models

12864_2025_11945_MOESM3_ESM.xlsx

Supplementary Material 3: Table S2. Functional enrichment analysis determined using Mercator4. A. FullyRhymicGenes: Transcripts from multiallelic genes fully rhythmic in either SD or LD compared to all rhythmic transcripts from multiallelic genes in leaf tissue. B. TuberDelayed: Tuber cycling genes with delays of at least two hours with respect to all tuber expressed genes. C. LD_delayed_6h: Transcripts with phase delays of at least six hours with respect to short days compared to all transcripts rhythmic in both photoperiods. D. SDLD_Amplitude_Difference: Transcripts displaying increase amplitude under short than under long day conditions (>0.5) with respect to all transcripts rhythmic under both photoperiods in leaves

12864_2025_11945_MOESM4_ESM.xlsx

Supplementary Material 4: Table S3. Identification of putative binding sites in potato genes. Binding sites from equivalent A. thaliana transcription factors were compared to the Atlantic promoter regions. One-sided Fisher’s exact tests were performed for testing enrichment in specific gene groups. A. TRB_statistics_summary: Enrichment of putative TRB binding sites in genes among genes in the MapMan 17.1 Ribosome biogenesis functional category with large changes in phase (Figure 3E). B. wrky_statistics_summary: Enrichment of putative WRKY binding sites in genes among genes MapMan4 26.10 Pathogen functional category with large changes in phase (Figure 3E)

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, 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 licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence 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 licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Feke, A., Vaillancourt, B., Acheson, K. et al. High resolution diel transcriptomes of autotetraploid potato reveal expression and sequence conservation among rhythmic genes. BMC Genomics 26, 925 (2025). https://doi.org/10.1186/s12864-025-11945-8

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-025-11945-8

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us