Three segregating F2 populations were planted and analyzed at the Newe Ya’ar Research Center in northern Israel: (1) ‘P.G.160-61 × P.G.100-1’ (n = 399) was constructed from a cross between ‘P.G.100-1’ (‘Wonderful’) a deciduous cultivar and (P.G.160-61) an evergreen cultivar originated from India. F1 plants, which were the progenitors of the F2 populations, were planted in 2006 and first flowers were self-pollinated at 2007. Five F2 populations were planted in 2008. (2) ‘Nana × Black’ (n = 295) was constructed from a cross between ‘Black’ (P.G.127–128) a deciduous cultivar and ‘Nana’ seedling selection (P.G.232-243, Punica granatum var. Nana) an evergreen non-deciduous cultivar. F1 plants, which were the progenitors of the F2 populations, were planted in 2008 and were self-pollinated at 2009. Four F2 population were planted in 2010. (3) ‘Nana × White’ (n = 267) was constructed from a cross between ‘White’ (P.G.164-65) a deciduous cultivar and evergreen ‘Nana’ seedling selection (P.G.232-243, P. granatum var. Nana). F1 plants, which were the progenitors of the F2 populations, were planted in 2010 and were self-pollinated at 2011. Three F2 populations were planted in 2012. The parents are part of the ARO collection that is located at the Newe Ya’ar Research Center¹ (Holland et al., 2009) and were called according to the nomenclature as reported in Ophir et al. (2014).

Advanced F3 (n = 355) and F4 (n = 176) populations were developed for ‘P.G.160-61 × P.G.100-1.’ Five F2 and two F3 heterozygous plants were self-hybridized to develop F3 and F4 generations, respectively.

Additional evergreen varieties from the ARO’s pomegranate germplasm that were studied: P.G.145-46, P.G.149-50 (additional ‘Nana’ seedling selection) and P.G.227-238 (cultivar with sour fruits with “Black” peel).

A plant was described as having the evergreen habit when leaves did not turn yellow and shed during autumn-winter months (November–January). In addition, when flowering and budding started early, toward the end of winter (February).

The evergreen habit of the ‘P.G.160-61 × P.G.100-1’ F2 population was monitored over 3 years (2009–2012). The ‘Nana × Black’ F2 population was monitored over 5 years (2011–2016) and the ‘Nana × White’ F2 populations was monitored over 2 years (2014–2015). F3 and F4 populations were described through autumn and winter seasons after 2–5 years from planting.

In order to analyze the evergreen trait as a quantitative trait, dates when almost all leaves of the tree (estimated 60%) became yellow (Figure 1) and when budding started were recorded. The dates for ‘P.G.160-61 × P.G.100-1’ F2 population were indicated at 2011–2012 season and for ‘Nana × Black’ F2 population at 2011–2012 and 2012–2013 seasons. The dates were inverted to quantitative values by calculating days from November 1st to the recorded date. In addition, subtraction of the recording date when leaves became yellow from the date of budding gave the number of days that trees were without leaves.

Chlorophyll content in the leaves was measured by chlorophyll content meter (CCM-200 plus, apogee, Logan, UT) along 21 months. Chlorophyll content is presented as chlorophyll concentration index (CCI) units. At least once a month, 35 representative leaves of two trees of the varieties P.G.100-1, P.G.160-61, and two F1 hybrids were measured.

Cutting of P.G.145-46, P.G.160-61, P.G.100-1 were taken from trees in the orchard after the accumulation of 47 chilling hours (less than 7°C), to a warmed room (26°C) with 16 and 8 h of light and dark, respectively. Each cultivar had three replicates with ten stems of 30 cm length. Date of 50% budding was recorded. Budding capability was assessed as the number of days to 50% budding.

The DNA-extraction protocol was based on Porebski et al. (1997) with few modifications as published in Ophir et al. (2014), young leaves (0.5 g) were resuspended in a 6 ml extraction buffer and the chloroform–octanol solution was replaced with chloroform–isoamyl alcohol. DNA was precipitated with sodium acetate instead of sodium chloride.

Crude DNA was extracted from young leaf (10 mg) in 96-well PCR plate according to a protocol based on Xin et al. (2003). DNA was diluted 1:100 for further analysis.

Means, standard deviations, trait distribution, pairwise correlation, ANOVA and Welch ANOVA, Shapiro–Wilk W Test, and Wilcoxon/Kruskal–Wallis Tests (Rank Sums) analyses were conducted with the JMP, v. 7.0 software (SAS Institute Inc., Cary, NC, United States).

Simple-sequence repeats markers were developed from genomic pomegranate sequences (Supplementary Table 3). SSR amplification reaction was carried out in 20 μl volume containing template ∼30 ng gDNA, 2.0X Apex Taq RED master mix kit (Genesee Scientific, El Cajon, CA, United States) and 0.4 μM each primer (Supplementary Table 3). Thermal cycling conditions were: 5 min at 94°C, followed by 30 amplification cycles at 94°C for 15 s, 58°C for 30 s, and 72°C for 30 s and then a final extension at 72°C for 45 min. Genotyping was performed using GENESCAN-500-LIZ size standard (Applied Biosystems, Foster City, CA, United States) on a 3130 xl Genetic Analyzer (Applied Biosystems, Foster City, CA, United States) and analyzed using the GeneMapper software (Applied Biosystems, Foster City, CA, United States).

Mapping was performed using JoinMap 3.0 software (Van Ooijen and Voorrips, 2001). Markers were grouped at a minimum LOD score of 4.0 and a recombination frequency value of 0.4. JoinMap 3.0 software uses Kosambi mapping functions to translate recombination frequency into map distance.

Using the MapQTL^® 5 software (Van Ooijen, 2004), QTLs and their significance were calculated using interval mapping (IM), multiple QTL model (MQM), and permutation analysis. A QTL was determined significant when its LOD score was higher than the calculated threshold (1000 permutation at p = 0.05).

Two varieties, P.G.100-1 and P.G.160-60, the parents of the segregating populations, as well as eighteen representatives of the F3 population (nine evergreen and nine deciduous) were selected for RNAseq analysis. Selected F3 plants had a homozygous genotype at the c26365, highly associated marker to the evergreen trait. In addition, they were described as evergreen or deciduous plants during 2 years (2013–2015). Chlorophyll content was examined for P.G.100-1 and P.G.160-60 and the selected F3 plants at three dates: 24/8/2016, 25/9/2016, and 9/11/2016.

A sample consisted of leaves collected at 25/9/2016 from different branches on the same tree. Tissues immediately frozen in liquid nitrogen were crushed and kept at -80°C until further analysis. Total RNA was extracted as described in Trainin et al. (2021) for each plant individually. Samples for RNAseq analysis consisted of three and two biological replicates (taken from different trees) of P.G.100-1 and P.G.160-61, respectively, and three mixed RNA bulks of evergreen and deciduous seedlings, each hold three different plants. Construction and sequencing of mRNA-seq libraries were performed at the Technion Genome Center.³ Briefly, mRNAseq libraries were generated from total RNA, using the TruSeq RNA protocol. The indexed libraries were pooled and subjected to sequencing on an Illumina HiSeq2500 instrument (Illumina, San Diego, CA, United States), 50 bp single read.

Raw-reads were subjected to a filtering and cleaning procedure. First, raw reads were filtered and using Trimmomatic (Bolger et al., 2014) to remove adapters. Next, the FASTX Toolkit⁴ was used to trim read-end nucleotides with quality scores < 30, using the FASTQ Quality Trimmer, and to remove reads with less than 70% base pairs with a quality score ≤ 30 using the FASTQ Quality Filter. Data was submitted to SRA bioproject PRJNA802080. Clean reads were mapped to the reference genome of P. granatum⁵ using Cufflinks RNA-Seq workflow⁶ using Tophat2 software (Kim et al., 2013) (v. 2.1) with an average mapping rate of 96.3%. Gene abundance estimation was performed using Cufflinks (Trapnell et al., 2010), (v. 2.2) combined with gene annotations (GCA_002201585.1_ASM220158v1). Gene expression values were computed as FPKM (Fragments per kilo base per million mapped reads). Differential expression analysis was completed using the DESeq2 R package (Love et al., 2014). Genes that were more than twofold differentially expressed with false discovery-corrected statistical significance of at most 0.05 were considered differentially expressed (Benjamini and Hochberg, 1995).

The gene sequences were used as a query term for a search of the NCBI non-redundant (nr) protein database that was carried out with the DIAMOND program (Buchfink et al., 2014). The search results were imported into Blast2GO version 4.0 for gene ontology (GO) assignments.

Transcripts were analyzed by the Integrative Genomics Viewer (IGV) interactive tool to visual explore their structure.

The PgPolyQ-MADS full length transcript sequence of pomegranate was cloned and sequenced from ripe arils of P.G.200-211, a deciduous cultivar with a light pink aril color originated from Spain (Harel-Beja et al., 2019). Primers used for the full-length transcribed gene were 5610_9 and 5610_6 (Supplementary Table 4). PCR products were sequenced with 3130 Genetic Analyzer (Applied Biosystems, Foster City, CA, United States). NCBI Blastp was used to find homologous proteins and conserved domains.⁷

Pomegranate PgPolyQ-MADS was cloned into the expression vector PBI121 (GenBank entry AF485783) under the control of the 35S promoter by Bio Basic Inc. (El Cajon, ON, Canada) (Supplementary Figure 1), the plasmid was designated PBI121-5610.

Agrobacterium tumefaciens strain GV3101 containing PBI121-5610 and PBI121 (control) were used to transform Arabidopsis (Col-0) by the floral dip in filtration method (Clough and Bent, 1998). Selection of transformants was conducted on 0.8% agar containing Murashige and Skoog (MS) salts (2.2 g/l) and kanamycin (50 μg/ml). Kanamycin-resistant seedlings (five of PBI121-5610 transformant and two of PBI121) were then transferred into soil to set seed. Progeny from self-fertilized primary transformants were grown in pots for phenotypic analyses at 20°C, 12 h light. PBI121-5610 transformant plants were verified by PCR with the primers 5610-5 and 35S_out (Supplementary Table 4). Control plants, transformants with PBI121, were verified by PCR with primers uidA forward and uidA reverse (Supplementary Table 4). Thirty-seven transgenic with PBI121-5610 Arabidopsis lines were obtained and characterized. Seeds were acquired from twelve positive lines for further analysis. Experiments were continued and repeated with two transgenic lines (#22 and #31) that showed strong characteristic phenotypes. Flowering time was recorded when first bolt had open flowers.

Structure model was constructed using the intensive mode of Phyre2 (Kelley et al., 2015). Visualization and analyses were done using PyMol⁸ (Schrödinger and DeLano, 2020).

The leaves of deciduous pomegranate change their color in the autumn, while the evergreen varieties retain their green leaves and continue to generate new buds even during the winter (Figure 1). The latter phenomenon was titled evergreen. In order to quantify the differences between the different pomegranate varieties through time, chlorophyll content in the leaves was measured in two varieties along 21 months (June 2014–February 2016). Measurements are represented as chlorophyll concentration index (CCI). While the leaves of P.G.100-1 and those of the F1 hybrids lose their chlorophyll from September through February, the leaves of the evergreen variety, P.G.160-61, remain green during all year (Figure 2). However, when temperatures dropped below 10°C for few days or more, the leaves of the evergreen defined cultivar P.G.160-61 lost their chlorophyll due to cold injury and dropped. Four additional accessions (P.G.145-45, P.G.227-238, P.G.149-50, and P.G.160-61) in the ARO pomegranate collection showed the same behavior. This phenomenon led to determine pomegranate as a conditional deciduous plant.

The differential effect of chilling on the budding of evergreen and deciduous pomegranates was studied on stem cuttings. Cutting of evergreen varieties P.G.145-46 and P.G.160-61 and deciduous variety P.G.100-1 were cut from trees in the orchard after accumulating 47 chilling hours (less than 7°C), and then warmed to 26°C. Date of 50% budding was recorded. Budding capability was assessed as the number of days to 50% budding. The evergreen varieties showed significantly earlier budding capability (Tukey–Kramer test) than the deciduous variety. 50% of the branches of P.G.145-46 and P.G.160-61 had buds after 23 and 39 days, respectively, while 50% branches of P.G.100-1 had buds only after 58 days (Supplementary Table 4).

In order to study the genetics of the evergreen trait, three different F2 populations, originated from five different parents (’P.G.160-61 × P.G.100-1,’ ‘Nana × Black,’ and ‘Nana × White’) were established. The distribution and segregation of the evergreen phenotype was analyzed in these populations through autumn and winter over 2 or 3 years (‘P.G.160-61 × P.G.100-1’ 2009–2012, ‘Nana × Black’ 2010–2012, and ‘Nana × White’ 2013–2015). Plants were described as evergreen when leaves remained green and did not drop during autumn-winter months (November–January), while flowering and budding started early, before spring (February). The evergreen phenotype was segregated in all the different F2 populations as a one major gene with significant (Chi square p > 0.05) 1:3 ratio (Table 1), where evergreen phenotype was recessive to the deciduous phenotype.

In order to find strongly linked markers to the evergreen trait, two divergent groups of plants were selected from the ‘P.G.160-61 × P.G.100-1’ F2 population. One group contained plants that expressed the evergreen phenotype and the second group expressed the deciduous phenotype. Comparative hybridization between pooled DNA from evergreen group and the deciduous group was performed on the SNP array, which was developed for pomegranate (Harel Beja et al., 2015). Differential signals were sorted from the strongest associated SNP to the least by sorting the differential signal ascending, yielding 133 statistically significant differential signals, which depicted SNPs (adjusted-p < 0.05), whereas 40 of them had greater confidence at (adjusted-p < 0.001) (Supplementary Table 6).

In addition, sequential SNP genotyping was performed using the Fluidigm technology with 201 SNPs for the ‘P.G.160-61 × P.G.100-1’ F2 population (n = 92). The markers used in this analysis included 45 SNPs with the highest scores (p < 0.0025) from the bulk array hybridization described above, three additional new SNPs and 153 SNPs from the pomegranate SNP subset (Harel Beja et al., 2015). Altogether, a total number of 181 SNPs were found appropriate for further analysis (Supplementary Table 2).

A linkage map was built for the ‘P.G.160-61 × P.G.100-1’ F2 population (Supplementary Table 7) in addition to the published ‘Nana × Black’ map (Harel Beja et al., 2015). The new map contains 167 markers and is divided to 12 linkage groups (LGs). 130 markers were common with the published map and their order was identical except for a few mismatches: one in LG3, three in LG5 and three in LG8. LG6 of ‘P.G.160-61 × P.G.100-1’ map contained three markers that were mapped at LG10 and one at LG11 on the ‘Nana × Black’ map.

The evergreen trait was mapped as a qualitative trait to LG3 at both linkage maps, at 98 and 163 cM (centimorgan) on ‘P.G.160-61 × P.G.100-1’ and ‘Nana × Black,’ respectively (Figure 3 and Supplementary Table 7). Four DNA markers were mapped at ‘P.G.160-61 × P.G.100-1’ map and two at the ‘Nana × Black’ map in the vicinity of the evergreen locus within approximately one cM.

A total of 70 and 38 of the bulk segregated associated SNP markers were mapped on ‘Nana × Black’ and ‘P.G.160-61 × P.G.100-1’ pomegranate linkage maps, respectively. They were mapped on 10 LGs at the ‘Nana × Black’ map, and at five LGs on the ‘P.G.160-61 × P.G.100-1’ map (Supplementary Table 7). The highly linked markers (p < 5^E–07) were mapped to LG3 at the region where the evergreen trait was mapped (153–172 cM at the ‘Nana × Black’ and 79–107 at the ‘P.G.160-61 × P.G.100-1’ map).

The evergreen phenotype was divided to three traits: leaves shedding, budding and the time between leaves shedding and budding. The dates of each trait were indicated for ‘P.G.160-61 × P.G.100-1’ F2 population at 2011–2012 season and for ‘Nana × Black’ F2 population at 2011–2012 and 2012–2013 seasons. Each trait was quantified by calculating the number of days between November, 1st and the date the trait was recorded. All three traits were mapped, by MQM analysis, as one major QTL at LG3 for both F2 segregating populations (Figure 4) with 56–89% explanation of the trait by IM analysis, depending on the character measured, year and population (Supplementary Table 8). The highest LOD values (LOD > 7.4) were related to the EverGreen loci and to adjacent markers. Additional minor QTLs (LOD > 3.8) for the three traits were analyzed for ‘P.G.160-61 × P.G.100-1’ F2 population, by IM analysis, at LG7 (Supplementary Table 8).

Markers that were mapped close to the evergreen loci, 22 and 12 cM, above and below, respectively (Figure 5A), were blasted (BLASTN) to two pomegranate genome sequences: GCA_002201585.1 assembly of ‘Dabenzi’ pomegranate and ASM765513v2 of ‘Tunisia’ (Supplementary Table 9). The markers were ordered on MTKT01004399 of GCA_002201585.1 assembly and to NC_045128.1 scaffold (Chromosome 2) of ASM765513v2 according to their order on the genetic map (Figure 5C).

The 34 cM interval that included the “EverGreen” loci covered a region of 4,099,398 bp at NC_045128.1 (905,484–5,004,872 bp) scaffold, including 577 genes. Identification of the relevant genomic sequences facilitated the development of additional DNA markers that were used to genotype plants within advanced segregating populations (Supplementary Table 3).

Eucalyptus and pomegranate are members of the Myrtales order. It was of interest to study the synteny between the genomes at the evergreen region. Sequences of genes of the evergreen trait on the genetic map (‘P.G.160-61 × P.G.100-1’ 73 –99 cM), that comprised markers of around 25 cM were blasted (BLASTN) to the Eucalyptus grandis genome V1.0. Interestingly, the order of these gene markers on the pomegranate genetic map and on the pomegranate ‘Tunisia’ genome was in synteny to that of the genes at chromosome 9 of E. grandis (scaffold9:24.4-27.2 Mb) (Figure 5B and Supplementary Table 9).

The matching genomic pomegranate and Eucalyptus sequences enabled the development of additional markers, including six SSR and eight SNPs (Supplementary Table 3). Genotyping with those markers enriched the marker collection within the close region for additional accuracy of the genetic map by additional recombinations.

Mapping the two different F2 populations, ‘P.G.160-61 × P.G.100-1’ and ‘Nana × Black,’ suggested that the evergreen trait was located between two close markers at LG3: c26365 and c8339. In order to reduce the area by additional recombinations, advanced segregating populations were developed. Five F2 and two F3 plants from ‘P.G.160-61 × P.G.100-1,’ that were found as heterozygous at the closed markers of the trait, were selected for self-hybridization to develop additional F3 and F4 populations with 355 and 176 plants, respectively. All populations segregated for the evergreen trait.

The highly associated markers, c26365 and c8339, were selected for genotyping the advanced populations. Interestingly, there were no recombinants between c26365 and the trait, within the F2 ‘P.G.160-61 × P.G.100-1’ population, although they were mapped with a gap of 5 cM. However, c8339 was mapped only one cM apart from the trait to the other direction with two recombinations. Moreover, within the F2 ‘Nana × Black’ population there were two recombinants between c26365 and one between c8339 and the trait, respectively (Table 2). The physical distance between those markers, according to of pomegranate genome assembly ASM765513v2 is 176,366 bp.

Within the advanced populations, two plants out of 243 and three out of 176 were recombinants between c26365 marker and the evergreen trait among the F3 and F4, respectively (Table 2). Additional markers (Supplementary Table 3) were developed at the associated region according to corresponding Eucalyptus and pomegranate genomic sequences. Genotyping with close developed markers was performed for selected recombinant plants consistent with the recombination position (Table 2). Altogether, eighteen DNA markers were genotyped within the evergreen region (Supplementary Table 3). All the markers that were genotyped, exposed recombination with the trait, except for marker 5610. Therefore, it is suggested that 5610 is highly genetically associated to evergreen trait.

Closely linked markers’ sequences were blasted (BLASTN) to P. granatum isolate Tunisia-2019, ASM765513v2 whole genome shotgun sequence. Based on the recombinants at the advanced segregating populations the genomic region was reduced to 19,000 bp on chromosome 2 (NC_045128.1), between marker 195_155 and 195_174 (Supplementary Table 3). The marker 5610 was associated to the evergreen trait over all recombinant plants. The genomic associated region included only one gene and it is LOC116195279 (XM_031524351.1; protein XP_031380211), which is annotated as agamous-like MADS-box protein AGL27. The gene is composed of seven exons that include two consensus sequences, a MADS_MEF2_like domain and a K-box consensus sequence (Figure 6). Interestingly, the gene contains an intron of almost 8,000 bp, which caused the separation of the gene into two genes (CDL15_Pgr005610 and CDL15_Pgr005611) at other pomegranate genome assembly (GCA_002201585.1 of ‘Dabenzi’ pomegranate). A repeat of eight successive glutamines (polyQ) is located at the N-terminal site of the deduced amino acid sequence, twenty-two amino acids subsequent the first methionine (Figure 6C). The gene was named PgPolyQ-MADS. The position of this small polyQ at the beginning of the protein could not be found in any other MADS-box proteins, by searching with BLASTP the available data bases.

Two polymorphic genetic changes were found between the deciduous P.G.100-1 and the evergreen P.G.160-61 at the PgPolyQ-MADS sequence. These changes are located within the CAG repeat located at the 5’ end of the gene (Figures 6D,E). The first change was that the evergreen P.G.160-61 had thirteen repeats, while P.G.100-1 had eight. Interestingly, the second change was at the sixth repeat of P.G.160-61, where a C was replaced by T, resulting in the formation of a stop codon (TAG) in the evergreen phenotype. P.G.160-61 and all other evergreen varieties present in the ARO pomegranate collection (P.G.145-46, P.G.149-50, P.G.227-238, and P.G.232-243) had the same genetic variation that caused a stop codon. Screening pomegranate accessions in the collection revealed that among the deciduous cultivars that were genotyped, twenty had identical allele to that of the PgPolyQ-MADS gene as P.G.100-1, one accession (P.G 221–232) had seven repeats and two seedlings from ‘Daru’ (P.G.263-274 and P.G.264-275) had nine repeats. However, none of the deciduous accessions showed the nucleotides replacements, which caused a stop codon (Supplementary Table 10).

In order to understand how the PgPolyQ-MADS box protein is influenced by the addition of the polyQ domain at the N terminal of the protein, the Phyere 2 program was used to model the protein (Figure 7). The monomeric protein is composed of two domains–DNA-binding domain and keratin like domain. The two domains were modeled with high confidence score, presenting up to 46% amino acid sequence identity to the templates used by Phyre2. The N-terminal polyQ tail and the C-terminal tail were predicted as disordered due to lack of available homologous structures for these areas. The model predicts that the polyQ domain is in close proximity to the DNA binding site and that its presence could influence its DNA binding capability. The polyQ domain might also play a role in targeting of other factors that could influence expression of the genes it binds to.

Arabidopsis transformants (Col-0) containing the PgPolyQ-MADS gene originated from a deciduous cultivar (P.G.200-211) had a distinctive phenotype (Figure 8A). Effects included smaller plants, fewer, smaller and rolled leaves, and early flowering. Twenty-five transformed plants showed at least one of the effects and twelve had no effect. Days to flowering were reduced significantly (p < 0.01) from 52 ± 30.3 days at the control plants (T1, n = 8) to 37 ± 1.43 for the PgPolyQ-MADS transformants (T1, n = 36, one plant did not flower) (Supplementary Table 11). Seedlings of two transformed lines with the PgPolyQ-MADS gene that had a distinguished phenotype were grown for two additional experiments. Interestingly, those lines were originated from two different T0 plants. The extreme phenotype repeated (Figure 8A) and the number of days to flowering was reduced significantly (Figure 8B).

Tracking chlorophyll content in the leaves along the year (Figure 2) showed that at the end of September deciduous varieties started to lose chlorophyll, while evergreen variety did not. That period was selected to study the different transcriptomes of deciduous and evergreen plants. P.G.100-1 and deciduous F3 progenies’ CCI values at 25.9.2016 (53.64 ± 12.52) were in-between the values at 24.8.2016 and 9.11.2016 (Figure 9). However, evergreen progenies and P.G.160-61 had mean CCI values of 67.76 ± 1.60 and were constant for all dates.

Transcriptome analysis was conducted for RNA extractions of leaves collected on 25.9.2016 from three P.G.100-1 trees and two P.G.160-61 trees. In addition, six bulks of F3 progenies from the ‘P.G.160-61 × Wonderful’ population were analyzed. Each bulk combined RNA from three different plants. Three bulks were from deciduous plants and three from evergreen plants. An average of 17,730,933 reads was obtained for each sample (Supplementary Table 12), with no considerable differences. Reads were mapped to pomegranate genome (GCA_002201585.1) where 32,582 contigs were analyzed.

Differential expression analysis showed forty-one transcripts that were differentially (p < 0.01) expressed between P.G.100-1 and P.G.160-61. These transcripts were also differentially expressed within the bulks of deciduous and bulks of the evergreen plants (Supplementary Table 13). Interestingly, these different genes included three genes from scaffold MTKT01004399.1 (GCA_002201585.1), which is associated to the evergreen trait. Those genes were annotated as acyltransferase and kinesin light chain. However, their location on the scaffold was not at the same loci as the associated markers were. Differential genes contain five transcription factors, including MOTHER of FT and TF 1 (CDL15_Pgr019923, LOC116212806) which is located at MTKT01006319.1 (GCA_002201585.1).

Interestingly, transcript CDL15_Pgr005611 deduced to encode for the PgPolyQ-MADS was not differentially expressed between deciduous and evergreen plants, although the transcript in the evergreen plants contained the additional nucleotide assumed to cause for a stop codon. Thus, the evergreen alteration does not seem to change the transcription level of the gene. Transcription levels for transcript CDL15_Pgr005611 were with an average of 51.6 and 66.7 FPKM (Fragments per kilo base per million mapped reads) in P.G.100-1 and the deciduous bulks, respectively. The expression of the evergreen plants was 44.6 and 35.3 for P.G.160-61 and evergreen bulks, respectively.