CHARACTERIZING MORPHOLOGICAL AND GENETIC VARIATION IN PRUNUS AND MALUS TO UNDERSTAND FLOWERING TIME REGULATION IN ROSACEOUS TREES By Charity Z Goeckeritz A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Plant Breeding, Genetics, and Biotechnology – Horticulture – Doctor of Philosophy 2023 ABSTRACT The Rosaceae plant family is one of the most economically important families for fruit crop production. Members include domesticated apple (Malus × domestica Borkh.) and sweet and sour cherry (Prunus avium and P. cerasus, respectively). Bloom time is an agronomically significant trait in these species since flowering quality directly translates to yield. If bloom time is too early or too late, the economic losses may be severe. A possible strategy to maximize flowering and fruit set in either scenario is to understand the genetic basis of bloom time and use this information to create cultivars with bloom times customized to certain climates. The genes and processes of flower development they affect are well-studied in model species, which has aided our understanding in species with similar floral phenology. However, the control of flowering time in fruit trees, which exhibit perennial growth and dormancy phases, is poorly understood. In the first chapter of this dissertation, the entire flowering process is reviewed for rosaceous fruit trees from molecular, physiological, and genetic perspectives. Gaps in our knowledge are identified and future directions of the field are described. In chapter two, I characterize the morphology of diverse Malus accessions showing extreme bloom time variation and speculate on the potential regulatory mechanisms considering the results. Chapter three describes the construction of a valuable resource for genomics-assisted breeding: a chromosome- scale genome of P. cerasus cv. Montmorency, which is the predominant sour cherry cultivar grown in the United States. In chapter four, I holistically compare floral biology in a bloom time segregating population of sour cherry for several years to identify where siblings deviate in their developmental journeys. Using this information, the ‘Montmorency’ genome is used to identify molecular changes underlying the bloom time differences with RNA-sequencing, differential expression comparisons, and network analysis. Moreover, I identify genomic variation associated with a late-blooming QTL haplotype, k, on Prunus chromosome 4. These variants are then intersected with the differential expression comparisons to obtain a set of high-quality candidate genes for the basis of late- bloom in this segregating population of sour cherry. This work allows a better understanding of the molecular and genetic control of flower development and bloom time in trees and lays a strong foundation for future work. This knowledge will inform researchers and fruit tree breeders on ways to modify, select for, or engineer desired bloom times to maximize crop yields. They will also provide insight into how temperate fruit trees may respond to a changing climate. Copyright by CHARITY Z GOECKERITZ 2023 ACKNOWLEDGEMENTS When I started college at the University of Utah in 2012, I went to the biology advisor, Dr. Dave Gard, to set up my curriculum. He asked me what I wanted to do, and I simply responded ‘plants.’ From then on ‘Plants’ was the header of every one of my advising summaries. I had grown up loving those little things – admiring them for their form, function, and resilience. So, the first to receive an acknowledgement will be plants – yes all plants – every single one that’s helped me get to where I am today. From the tiny orange tree that my baby brother ripped all the leaves off, to the shrubs, flowers, and trees I recommended for customers while working at Glover Nursery. I hope they are someplace flourishing and making folks happy. And of course, I am eternally grateful for the majestic trees from the ‘Balaton’ x ‘Surefire’ cross created by Dr. Amy Iezzoni, the Tart Cherry Queen Mother. Sorry pop4 for killing thousands of your flowers during my PhD – to be fair, most of you wouldn’t have become fruits anyway on account of your terrible fertility. I’d like to thank my family – especially my dear sisters. Shannon, Kim, Stacey, Katie, Hope, Faith, and Grace have always been in my corner. We’ve all become strong, accomplished women and I’m so proud of us. I’d like to thank my mothers as well, Laurie and Laura. I have always admired Laurie for her kind and selfless spirit, and Laura for her hard-working nature and determination. Laurie introduced me to gardening, and Laura’s childhood nickname for me was ‘cherry pie.’ I suppose my fate was written in the stars. Thank you to Jeremy, Scott, and Eric for getting me into videogames. They sparked my love for puzzles, and (perhaps detrimentally) cultivated my competitive nature. Mikey, Micah, and Timmy – thanks for being one of the best things about visiting home. When I got accepted into Michigan State in February 2017, I remember my dad asking me ‘Are you really going to move to Michigan!?’ ‘Yes Dad, I really am!’ A month later, he unexpectedly passed away. I so very wish he could see me now – I hope he would be proud. I remember my dad supported me financially for my first semester and he always encouraged me academically. I also owe him for nurturing my love for the natural world. I have very fond memories of camping with my family, as it was one of Dad’s favorite things to do. Wherever you are, Dad, I miss you, I love you – and thanks for helping me with my high school math homework. iv I also have many friends to thank. First, my high school mates Shelley and Tina: thanks for always making me laugh. During my time at Glover Nursery, Monica Moser, Cassandra Becker, and Stacie Gerdes were my rock – I think of you all often. Thank you to Katie Sanbonmatsu and Christina Kilbane for pulling through that Biology degree with me. Thank you to my fellow teaching assistants, particularly the ones I met while teaching Biochemistry; grading papers until 2:00am really trauma bonds a group, and I’m still not sure the free food was worth it. To all the folks I met volunteering and touring in Thailand, I hope you are doing well. I may not remember all your names, but I remember your faces and I’ll always remember Nightmare Island. A mentor of mine told me once that you’ll make some of the best friends of your life in graduate school – he was right. A special shoutout to my dear friends Alan Yocca, Beth Alger, Kevin Bird, Jeremy Pardo, Phil Engelgau, and Davis Mathieu – you all held me together during my first couple years of graduate school. I miss our karaoke nights. Thank you to my friend, colleague, confidant, and doppelganger Kathleen Rhoades. I am eternally grateful for my Hollender lab mate and sister in science, Andrea Kohler, and for her reassurance and feedback all these years. Thank you to my fellow Utahn, Sarah Lee. And finally, thank you to my lifelong friend Jonathan Quincy. You are truly one of the kindest people I know, and I would never be where I am without you. Although there are more folks to thank, it is required that I have some scientific writing in this dissertation – so if you feel I have wrongfully excluded you, please let me know and I will promptly deliver apology cookies to you. I have had exceptional mentors in my life thus far. My 4th grade teacher, Kaye Flanery, was a true cheerleader to me. Dr. Dave Gard encouraged me to go to graduate school and would frequently say ‘remember that science is done by people.’ To believe science exists outside of the human condition is senseless – we must constantly recognize our own biases, acknowledge our mistakes, and be kind and inclusive. Drs. Chris Gottschalk, Ben Mansfeld, and Joseph Hill helped me find my footing when I first started graduate school. And of course, thank you to my outstanding committee members for the guidance over the years. Drs. Courtney Hollender, Rebecca Grumet, Robert VanBuren, Amy Iezzoni, and Chad Niederhuth – I admire each and every one of you for your kindness, intelligence, and dedication. I have truly loved getting to know you during my graduate career, and I feel honored to have been your student. And another special thank you to my advisor, Dr. Courtney Hollender. Your enthusiasm is infectious, and your support has been unparalleled. v Last but not least, thank you to my little family: my three cats Luna, Jiji, and Emerson, and my life partner Indiana (Indy) Uribe. Partners of graduate students truly deserve more credit – Indy was the one out there lighting propane torches all night to save my data one year, and he’s brought me back from a stress implosion or two. Thank you all for shaping the person and scientist I am today. vi TABLE OF CONTENTS CHAPTER 1…………………….………………………………………………………………...1 CHAPTER 2………………………………………………………………………………………3 CHAPTER 3………………………………………………………………………………………5 CHAPTER 4 ……………………………………………………………………………………...7 REFERENCES…………………………………………………………………………………...67 APPENDIX………………………………………………………………………………………78 vii CHAPTER 1 The work presented in this chapter is part of the final publication: Goeckeritz C & Hollender CA. There is more to flowering than those DAM genes: the biology behind bloom in rosaceous fruit trees. Curr Opin Plant Biol 2021;59:1–10. doi.org/10.1016/j.pbi.2020.101995 Author contributions: Both CG and CH conceptualized and prepared the manuscript. CG generated the figures. 1 Abstract The regulation of bloom time in deciduous fruit trees is an area of increasing interest due to the negative impact of climate change on fruit production. Although flower development has been well-studied in model species, there are many knowledge gaps about this process in perennial fruit trees, whose floral development spans the four seasons and includes many temperature-driven transitions. To develop solutions for minimizing crop loss, a comprehensive research strategy is needed to understand flower development and bloom time in deciduous fruit trees. This approach must incorporate genetic, physiological, and phenological strategies which include morphological and molecular analyses. Here, we describe key floral development events for rosaceae family fruit trees, highlight recent molecular and genetic discoveries, and discuss future directions for this field. Summary For this chapter and published article, I conducted a thorough literature review to identify gaps in our knowledge of bloom time in Rosaceous fruit trees. Courtney and I worked together to decide physiological and genetic discussion topics from floral initiation until anthesis, a broad perspective not often taken in the field. Here we argue such a holistic view is needed to fully understand and design strategies for manipulating bloom time to ensure a reliable, abundant crop each year. I generated all figures. 2 CHAPTER 2 The work presented in this chapter is part of the final publication: Goeckeritz CZ, Gottschalk C, van Nocker S and Hollender CA. Malus species with diverse bloom times exhibit variable rates of floral development. J Amer Soc Hort Sci 2023;148:64–73. doi.org/10.21273/JASHS05236-22 Author contributions: CAH and SVN conceptualized the experiments. CZG and CG collected materials and CZG performed the experiments. All authors participated in writing the manuscript. 3 Abstract In response to challenges caused by climate change, apple (Malus ×domestica) breeding programs must quickly develop more resilient cultivars. One strategy is to breed for various bloom times. Members of the genus Malus, including domesticated apple, wild species, and hybrids, exhibit striking variations in the bloom date. Although bloom time is strongly influenced by chilling requirements, other aspects of floral development in Malus and their contributions to bloom time are less known. The purpose of this study was to investigate potential connections between predormancy flower development and final bloom time in Malus species. We performed a phenological analysis of flower development in wild and domesticated apple with extreme differences in bloom time over the course of one developmental season. We tracked histological changes in the floral apex of representatives of three early-blooming Malus genotypes (M. ×domestica ‘Anna’ PI 280400, M. orthocarpa PI 589392, M. sylvestris PI 633824) and three late-blooming genotypes (M. angustifolia PI 589763, M. angustifolia PI 613880, M.×domestica ‘Koningszuur’ PI 188517). Our study documented their floral meristem progression and organ development and expanded on current staging systems for apple flower development to describe the changes observed. The developmental trajectories of each genotype did not group according to bloom category, and we observed variations in the floral development stage at the time of dormancy onset. Summary In this chapter and published article, my fellow co-authors and I investigate the question of whether Malus species showing contrasting final bloom times show differences in development before dormancy. Such a question has been mostly overlooked in the field, and researchers have historically focused much more on dormancy’s contributions to bloom. Our results suggest the developmental phases prior to winter should also be accounted for when considering the biology of bloom time. SVN and CH conceived the study, I collected materials with CG, and I performed all sectioning experiments and interpretations. I expanded upon a current staging system for apple to properly compare the different species, and CH and I created the figures. All authors participated in writing the manuscript. 4 CHAPTER 3 The work presented in this chapter is part of the final publication: Goeckeritz CZ, Rhoades KE, Childs KL, Iezzoni AF, VanBuren R and Hollender CA. Genome of tetraploid sour cherry (Prunus cerasus L.) ‘Montmorency’ identifies three distinct ancestral Prunus genomes. Hortic Res 2023;10:uhad097. doi.org/10.1093/hr/uhad097 Author contributions: CAH, AFI, and RV conceptualized the experiments. CZG performed genome assembly, annotation, subgenome assignment using orthologs, and divergence time estimate analyses. KBR performed the k-mer hierarchal clustering, Ks analysis, synteny, DAM gene phylogenetic analyses, and S-allele alignments. KLC provided expertise and assistance with annotation and contributed code. CZG, KBR, and AFI wrote the manuscript. All authors assisted with editing the manuscript. 5 Abstract Sour cherry (Prunus cerasus L.) is a valuable fruit crop in the Rosaceae family and a hybrid between progenitors closely related to extant Prunus fruticosa (ground cherry) and Prunus avium (sweet cherry). Here we report a chromosome-scale genome assembly for sour cherry cultivar Montmorency, the predominant cultivar grown in the USA. We also generated a draft assembly of P. fruticosa to use alongside a published P. avium sequence for syntelog-based subgenome assignments for ‘Montmorency’ and provide compelling evidence P. fruticosa is also an allotetraploid. Using hierarchal k-mer clustering and phylogenomics, we show ‘Montmorency’ is trigenomic, containing two distinct subgenomes inherited from a P. fruticosa- like ancestor (A and A’) and two copies of the same subgenome inherited from a P. avium-like ancestor (BB). The genome composition of ‘Montmorency’ is AA’BB and little-to-no recombination has occurred between progenitor subgenomes (A/A’ and B). In Prunus, two known classes of genes are important to breeding strategies: the self-incompatibility loci (S- alleles), which determine compatible crosses, successful fertilization, and fruit set, and the Dormancy Associated MADS-box genes (DAMs), which strongly affect dormancy transitions and flowering time. The S-alleles and DAMs in ‘Montmorency’ and P. fruticosa were manually annotated and support subgenome assignments. Lastly, the hybridization event ‘Montmorency’ descended from was estimated to have occurred less than 1.61 million years ago, making sour cherry a relatively recent allotetraploid. The ‘Montmorency’ genome highlights the evolutionary complexity of the genus Prunus and will inform future breeding strategies for sour cherry, comparative genomics in the Rosaceae, and questions regarding neopolyploidy. Summary Michigan is the number one producer of sour cherry in the United States, with the cultivar ‘Montmorency’ being most widely grown. It is also home to the only breeding program for the crop, so a group of researchers, me included, worked together to create a genome for sour cherry to facilitate genomics-assisted improvement. I performed genome assembly, annotation, quality assessments, subgenome assignments, and evolutionary divergence dating. The genome for sour cherry turned out to be more complicated than originally thought, having three subgenomes instead of the expected two – therefore I also conducted quite a bit of troubleshooting. AF, KBR, and I wrote the manuscript, and all authors edited it. KBR and I generated the figures. AF, CAH, and RV conceptualized the work. 6 CHAPTER 4 The work presented in this chapter is being prepared for the publication: Goeckeritz CZ, Grabb C, Grumet R, Iezzoni A, Hollender CA. A thorough assessment of flower development in sour cherry (Prunus cerasus) siblings with variable bloom times and associated genetic variation. In prep Author contributions: CZG, CH, AI, and RG designed the experiments. CZG completed all experiments and data analysis, and CG contributed to histological work. CZG wrote the manuscript. 7 Abstract Bloom time is central to reproductive fitness, critical for agriculture, and a quantitatively- inherited trait resulting from numerous processes. It is acutely important for the sour cherry (Prunus cerasus) industry, which periodically suffers huge crop losses due to late spring freezes. Prunus species, like other rosaceous fruit trees, begin floral development in the summer nearly eight months before blooming the following spring. Closely studying this floral biology, especially in segregating populations, will inform future cultural practices and genetic strategies to adjust bloom times for certain climates. Here we build from a previous study where a major QTL on linkage group 4 was identified in a population of sour cherries segregating for bloom time. An examination of the full development of flowers of some of these trees indicated early- blooming siblings show developmental differences as soon as floral commitment; they maintain this lead until blooming the following spring. Using this information, we sequenced RNA from three trees in each bloom group to capture gene expression changes as siblings transitioned from a vegetative to reproductive apex. A gDNA coverage analysis revealed the late-blooming haplotype of this QTL, k, is located on sour cherry subgenome A, a subgenome donated by the late-blooming P. fruticosa-like progenitor. Furthermore, a weighted gene co-expression network analysis and gene ontology enrichment tests were performed to provide an overview of the biological processes modified between bloom groups. Four modules were recognized as majorly or moderately associated with the bloom time phenotype since genes within these groups showed distinctive expression patterns between early- vs late-blooming siblings. The strongest candidate genes for the late bloom phenotype (k) include those that show unique genetic variation in individuals with the k haplotype, show differential expression between bloom groups at the time of developmental divergence, and are part of a network module showing patterns specific to each bloom group. This includes several terpene synthases (Pcer_097523, Pcer_024813, Pcer_097525) and transcription factors with homology to the REproductive Meristem (REM) B3 domain-containing proteins (Pcer_097451, Pcer_028683, Pcer_024686). We discuss these results in the context of the complex phenotype observed in early- versus late-blooming siblings, relating our findings to a recently fine-mapped chromosome 4 QTL in sweet cherry. We provide evidence that the basis of k is separate from the candidate gene(s) found in the sweet cherry mapping study – suggesting several major regulators of bloom time are located on Prunus chromosome 4. 8 Introduction Prunus cerasus L. (sour cherry) is an agronomically important temperate fruit tree in the Rosaceae family, which also contains several other crops including sweet cherry, peach, apple, and pear [1]. Sour cherries are prized for their uniquely sweet and acidic flavor, and their processing qualities make them ideal for juice, jam, compote and pie. A critical trait for this industry is bloom time. In addition to its significance for successful cross pollination, blooming too early or too late may result in devastating crop loss. For instance, in 2012 a series of unexpectedly warm days in late March pushed flowers far past the typical stage of development in fruit trees across the Midwest. Michigan, the number one producer of sour cherries in the nation, was no exception. Shortly after this hot spell, a late spring freeze killed vulnerable flowers and annihilated over 90% of the expected crop [2]. Unfortunately, this was not an isolated incident, as spring freezes caused significant losses for 2002, 2007, 2020 and 2021. Moreover, the unpredictability of climate change is anticipated to increase the frequency of these events [2–4]. Therefore, strategies ranging from cultural practices to genetic improvement are sorely needed for a healthy and stable market. Sour cherry is a recent allotetraploid with progenitors closely related to Prunus avium (sweet cherry) and Prunus fruticosa (ground cherry) [5–8]. The recent publication of a reference genome for the cultivar Montmorency revealed three distinct subgenomes: A, A’, and B, with subgenome B at twice the dosage of A and A’. The two copies of subgenome B were inherited from the P. avium-like ancestor while subgenomes A and A’ originated from the P. fruticosa-like ancestor, which was also shown to be an allotetraploid of two unknown Prunus species [9]. Prunus avium originates from the Mediterranean region south of the Black Sea, while P. fruticosa is typically found at more northern latitudes including Ukraine, Russia, and Kazakhstan [10]. The latter shows some of the latest bloom times of the entire Prunus genus; in a common orchard, some P. avium accessions will bloom weeks earlier than P. fruticosa. Given their native habitats, this is consistent with observations in other species that show environmental signals such as photoperiod, temperature, and latitude have major implications for bloom time [11–13]. Comparatively, sour cherry genotypes exhibit diverse bloom times spanning the ranges of both progenitors. Therefore, fine-tuning bloom time in sour cherry using genetic approaches is a promising approach, especially given its newly-available genome sequence. 9 Bloom time is a highly heritable, quantitative trait that is a collection of many sequential developmental processes [14,15]. A clear, comprehensive view of floral development from initiation to anthesis is needed to attribute phenotypic variation to genetic variation. This tactic also provides information useful for designing new cultural practices to manipulate bloom time. As for other Rosaceous tree fruits, sour cherry begins its flower development in the summer [16,17]. The apices in lateral and spur positions of the new year’s growth are initially vegetative but have the potential to eventually become flowers [16]. Although the molecular triggers for the floral switch are not fully understood, they have become clearer in the past decade thanks to transcriptomic data and hormonal profiling. In transcriptome comparisons of biennial versus regular-bearing apple cultivars (Malus × domestica), researchers implicated a low gibberellic acid (GA) to cytokinin (CK) ratio and sugar-based signaling for the floral transition [18,19]. RNA-sequencing and hormone profiling in sweet cherry also support the positive role of CK for floral initiation in Rosaceous fruit trees, as well as positive associations of increased abscisic acid (ABA) and indole-acetic acid (IAA; auxin) [20,21]. Other transcriptional studies and applications of plant growth regulators further support GA’s inhibitory role on floral initiation [18,22–24]. However, such studies should be interpreted with caution since hormones’ effects are spatially, temporally, and structurally dependent [25–27]. Finally, a few studies have shown variation between years for floral development of apple and peach, suggesting photoperiod may not be a strong environmental cue for initiation – at least for these genotypes [28–30]. However, physiological evidence indicates that there are strong temperature effects on this developmental process, with different species exhibiting distinct optimum temperatures for floral initiation [31,32]. The question remains whether genotypes within a species show variation for this optimum temperature. The summer is spent producing floral organs, which involves the activity of MADS-box transcription factors including members of the SHORT VEGETATIVE PHASE (SVP), SEPALLATA (SEP), PISTILLATA (PI), and AGAMOUS (AG) subfamilies [20]. Additionally, APETALA 2/ETHYLENE RESPONSE FATOR (AP2/ERF), MYB, basic HELIX-LOOP-HELIX (bHLH), and NAC transcriptions factors are also expressed at this time of development [20,21,33]. ABA and CK are at their highest concentrations around the time of floral induction but then decline as organogenesis proceeds – IAA mirrors these trends in sweet cherry but shows a similar one in apple [18,20,21]. As the days grow shorter and temperatures decline, 10 development gradually slows, and trees enter a state of dormancy. For Malus and Pyrus, low temperatures appear to be the key environmental inducer of dormancy, while Prunus seems to monitor both temperature and photoperiod [34,35]. Pioneering work in poplar exploring tree physiology during dormancy transitions has greatly furthered our understanding of this biological process. Decreasing daylength causes the build-up of callose in phloem plasmodesmata, obstructing transport of growth-promoting signals to the apical meristem. In poplar, this process was shown to be positively regulated by ABA through a homolog of SVP, which moves on to negatively affect the expression of GA20 oxidases (gibberellin biosynthetic enzymes) and FT1 which is required for proper budbreak [36–39]. Traditionally winter dormancy has been categorized as two discrete phases: endodormancy and ecodormancy [40]. To relieve the growth inhibition wrought by endodormancy, chilling must be sensed by the quiescent meristem; for ecodormancy, warmth must return for growth resumption and budbreak. There is evidence to suggest both phases are controlled by genetics, although chilling, or length of endodormancy, has a much higher and more reliable broad sense heritability [41–45]. ‘Chill requirements’ refers to the amount of cold a genotype must experience to reach 50% budbreak after some specified amount of time in forcing conditions, and it is usually calculated using one of three models: the Chilling Hours model, the Utah Model, or the Dynamic Model [46]. The Chilling Hours model is the simplest: hours within a certain range (usually 0°C – 7.2°C) equally count towards a meristem’s chilling requirements [47]. The Utah Model also tallies chilling hours within a certain range, but each temperature is weighted differently; moreover, those above the range may contribute negatively towards chilling requirements [48]. The most complex yet most accurate model across environments is the Dynamic Model, which measures chill in units called ‘chill portions’ [49]. A chill portion is a theoretical, biological product that is created during optimum temperatures for chilling. Once a chilling portion is created, it cannot be destroyed; however, as the product is actively made, warm temperatures have the potential to negate the process [50,51]. Irrespective of the model chosen to calculate chill, researchers typically switch to tallying ‘heat requirements’ once chilling needs are deemed fulfilled, usually for practical purposes and to simplify comparisons between genotypes. However, abundant evidence demonstrates further chilling beyond the defined requirements increases heat’s effectiveness on development, and that historical estimations of chilling are insufficient for full, synchronized bloom and commercial productivity 11 [52,53]. The gradual opening of plasmodesmata with chilling provides an intriguing explanation for controlling the transport of growth-promoting signals and subsequent responsiveness to heat. This process may act as a natural throttle for starch and hexose mobilization, which fluctuates during different dormancy phases [54–57]. As alluded to previously, ABA is considered a dormancy maintenance hormone while GA appears to be involved in dormancy breakage through activation of genes essential to basic metabolism and growth [58–60]. Plentiful transcriptomic data in a variety of Rosaceous fruit trees support this declaration, and the antagonism between these two hormones in the context of dormancy is a key discussion point of many reviews [21,33,58,59,61–64]. Less studied but gaining traction is the proposed influence reactive oxygen species (ROS) have on the ABA:GA interplay, dormancy, and the cell cycle. For years, the nut and fruit tree industry have used a commercial product called Dormex® to increase synchronized budbreak in the absence of sufficient chilling. The main active ingredient in Dormex® is hydrogen cyanimide, and molecular investigations into its mechanism of action links a ROS burst to the upregulation of enzymes involved in ROS scavenging and the synthesis of antioxidant compounds [65,66]. Indeed, at non-lethal levels, ROS such as H2O2 and O2 • - have been likened to hormones since they have potent effects on genes related to cell wall-loosening, expansion, and cycling [65–68]. Several studies in perennials have shown ROS accumulation associates with deeper dormancy and steadily declines as the meristem becomes more responsive to growth-promoting temperatures, a trend quite like ABA and opposite to GA [57,58,69–71]. Of course, no discussion of dormancy is complete without mentioning the six tandemly- arrayed Dormancy Associated MADs-box genes (DAMs); these were discovered in a peach genotype that failed to maintain dormancy nearly 30 years ago and were later mapped to a large deletion on chromosome 1 [72,73]. Since then, these genes’ roles in dormancy transitions have been intensely studied through transcriptomic, epigenomic, and transgenic means [74–80]. These genes are members of the SVP clade of MADs-box transcription factors, and orthologs have been identified and functionally characterized in several other species aside from peach [38,81,82]. Although tandemly-arrayed, these genes appear to have distinct expression patterns, with several connected to higher chilling needs and later bloom times [76,78,80,83,84]. Several findings have shown positive feedback loops between ABA levels and DAM gene expression [38,58,80]. In contrast, cytokinin signaling was linked with the downregulation of an apple DAM1 ortholog and 12 dormancy release [85]. Given these complex mechanisms controlling dormancy, it is certainly more appropriate to think of it as not a stall in growth, but an active process to prevent the waste of resources and certain death in the face of unfavorable conditions. After the meristem accumulates sufficient chill and warmth returns the following spring, flowers grow exponentially until final bloom. GA, sugar transport, and genes involved in basic metabolic processes like cell division and differentiation are all crucial to growth resumption in spring [21,60,70]. In Prunus, the DAM genes have been independently attributed to chilling perception, dormancy transitions, and contrasting bloom times in Quantitative Trait Locus (QTL) mapping and other genetic association studies [41,80,84]. However, they are certainly not the full story. QTLs have been found on all Prunus chromosomes, with large effect QTLs frequently found on chromosome 4 [41,42,44,45,86,87]. Prunus is well-known for its genome synteny – so like the DAM genes, it is hoped that future gene discoveries will expand our understanding of bloom time in sister species as well [45,84]. Indeed, some of the QTLs on chromosome 4 overlap in syntenic regions, perhaps hinting to common genes and mechanisms controlling bloom [42,45,86–89]. Attempts have been made to determine candidate genes behind the basis of bloom time on chromosome 4, but long generational times and few genomic resources for Prunus have made this endeavor challenging. Without genomes of rosaceous fruit trees, these earlier studies were limited in their scope and narrowly-focused on homologs known to affect flowering time in annual species. For instance, a SQUAMOSA PROMOTER-BINDING LIKE gene (SPL5), LFY, SUPPRESSOR OF OVEREXPRESSION OF CONSTANS (SOC1), and PHOTOPERIOD INSENSITIVE EARLY-FLOWERING (PIE1) have all been cherry-picked as candidates over the years [88,90–94]. Fortunately, resources are rapidly being generated as genome sequencing becomes routine [9,95–99]. Recently, researchers successfully fine-mapped a chromosome 4 QTL in Prunus avium (sweet cherry) to just 12 genes within a 68 kb region in the ‘Regina’ genome [45]. Alas, this locus does not overlap with any of the candidate genes proposed in previous studies – a cautionary tale for choosing candidate genes based solely on homology [100]. In the following work, we take a holistic approach to identify candidate genes in a population of sour cherry individuals segregating for bloom time [87]. This F1 population (n = 76) is the result of a cross between ‘Balaton’ and ‘Surefire,’ and a large (c. 8 Mb) but robust QTL on linkage group 4 named qP-BD4.1m explained 27.9% of the phenotypic variation in 13 bloom time. Marker phasing led to the identification of seven haplotypes of this QTL. One haplotype in particular, k, was shown to delay average bloom time in an additive manner by approximately three days. Each tetraploid parent contained one copy of k – thus, individuals inheriting 2k were considered late-blooming (LB) while those inheriting 0k were considered early-blooming (EB). In contrast with past studies, we thoroughly assessed the floral biology in 2k and 0k full-siblings for three years to understand when and how these individuals differed in their developmental trajectories. The comprehensive picture of flowering phenology in these individuals led to the creation of a staging system that will be easily adaptable for other Prunus species. We also sequenced the DNA of some of these individuals to identify genetic variants unique to the k haplotype. Lastly, based on our dissection of the bloom time phenotype, we sequenced the RNA of LB and EB siblings and conducted differential gene expression and co- expression network analyses to select high-quality candidate genes underlying the late-blooming trait. We believe our results significantly enhance our understanding of this complex trait in Prunus since it was guided by a complete picture of the flowering process. Candidate genes are anticipated to be steppingstones for further research. In short, our work suggests there is indeed more to flowering than those DAM genes. 14 Figure 4.1 Floral bud photos of three late-blooming (inherited 2 k haplotypes) and three early- blooming individuals (inherited 0 copies of k) on the same day for two years. The top three individuals for each date are Early 1, Early 2, and Early 4, respectively. The bottom three individuals for each date are Late 3, Late 1, Late 2, respectively. 15 Depiction Stage Description and Number Depiction Stage Description and Number Depiction Stage Description and Number 1 Meristem is vegetative, continually initiating bud scales. There are distinct L1 (purple), L2 (blue), and L3 (pink) layers. The central zone (yellow) is highlighted. 2 Increased cellular divisions in the central zone widen and raise the meristem. There are distinct L1 (purple), L2 (blue), and L3 (pink) layers. The central zone (yellow) is highlighted. 3 The inflorescence meristem initiates lateral floral meristems (green). 4 Several floral meristems (green) make up the inflorescence meristem. Bracts are associated with individual floral meristems (yellow) and the whole inflorescence (blue). 5 Sepals (green) are initiated on the outermost regions of each floral meristem. Floral and inflorescence bracts are shown in yellow and blue, respectively. 6 The meristematic region (pink) widens and flattens while sepals elongate (green). Floral bracts are shown in yellow. 7 Sepals continue to elongate (green) while petal primordia (pink) are initiated. 8 As sepals (green) and petals (pink) continue development, stamen primordia (orange) initiate. Typically there will be 2 whorls of 10 stamens. 9 As sepals (green), petals (pink), and stamens (orange) continue to develop; carpel primordia (blue) initiate. 10 Sepals (green), petals (pink), and stamens (orange) continue development as the plicate carpel (blue) begins to fuse near the apical region. 11 The carpel (blue) closes, placing the carpel margin meristem (CMM; yellow arrow) in the adaxial region of the fusion area. Anther (orange) and filament (brown) primordia are distinct. 12 Inside the ovary (blue), the CMM (lavender) is distinguishable from the placental wall. Elongated cells may emerge from the CMM. In the anthers, sporogenous tissue (yellow) is distinct from other layers. 13 Pollen mother cells (PMCs; yellow) enlarge. Inside the ovary, the ovule primordia (lavender) emerge. Papillae cells (pale blue) may be distinct from neighboring epidermal cells. Elongated cells are usually seen by this stage. 14 Ovule primordia (lavender) are large, round protuberances. PMCs (yellow) degrade their cell walls in preparation for microsporogenesis, indicated by their rounded shape. The tapetum (pink) cytoplasm becomes more dense. 15 The integument (INT) primordia (lavender) initiate; the megaspore mother cell (MMC; dark green) is in the central position of the nucellus. Microsporogenesis begins; PMCs (yellow) have a thick wall of callose. 16 INTs (lavender) begin surrounding the nucellus; the MMC enlarges (dark green). Microsporogenesis is nearly complete; tetrads (yellow) are present and the tapetum (pink) begins to degrade. 17 INTs (lavender) continue to engulf nucellus; the MMC (dark green) initiates megasporogenesis. Developing pollen grains (yellow) assume a tricolpate shape and the tapetum (pink) further degrades. 18 INTs (lavender) nearly engulf nucellus; The MMC has completed meiosis I when two thick-walled cells appear (dark green). Tapetum (pink) is negligible, pollen grains (yellow) continue development. 19 The INTs are distinguishable at micropyle (inner in lavender, outer in dark lavender); the MMC completes meiosis II and one functional megaspore (magenta) arises from the tetrad (dark green). Tapetum is fully degraded; pollen (yellow) development continues. 20 The functional megaspore undergoes mitotic divisions to form the embryo sac (magenta); 1 - 6 nuclei may be seen in addition to remnants of the tetrad (dark green). Pollen grains (yellow) are mature. 21 The fully developed embryo sac (magenta) is formed, with 7 - 8 nuclei; the synergids and egg cell at the micropylar end have their own membranes. Pollen grains (yellow) are mature. Figure 4.2 A floral development staging system for sour cherry describing development pre- floral initiation to just before anthesis. A stage number and description are given to the right of each respective depiction. Colors are used to highlight features important for stage classification. 16 Figure 4.3 Photos of histological sections and hand dissections illustrating each floral stage described in Figure 4.2. An * indicates these section photos are identical to those in Figure 4.8E. A ** means the anther / pollen photograph is identical to that in Figure 4.8A, and the embryo sac is a different section of the same sample also in Figure 4.8A. Scale bars are 200 µm unless stated otherwise. 17 Results The average range in bloom times between individuals of the segregating population with 0k and 2k was 72 GDDs and 10 days for the three years of this study. Three representative early- bloomers (0k) and three late-bloomers (2k) photographed on the same day in two of these years are shown in Figure 4.1. Unique identifiers for each tree along with bloom times (50% open flowers) in calendar days and Growing Degree Days (GDDs) in the years they were scored are shown in Table S4.1. Materials collected and corresponding dates are fully documented in Tables S4.2 – S4.5. Our preliminary observations of segregating individuals in the late winter of 2018 revealed slight differences in internal development between early-blooming (EB) and late- blooming (LB) siblings even as the flowers were still dormant (data not shown). For example, EB anthers were clearly further advanced based on a more defined shape and differentiation of cell layers. Therefore, flowers/apices in the floral position of several early-blooming and late- blooming sour cherry trees from the same population were collected from June to May of the 2018-2019 and 2019-2020 seasons to determine how and when these developmental changes occurred. For effective comparisons between early- and late-blooming siblings, we developed a novel floral development staging system for sour cherry totaling 21 stages that span the entire life of the flower (Figure 4.2). This system was created since the literature for Rosaceous fruit trees’ floral morphology is somewhat fractionated. Previous works usually describe only a portion of floral development and are often lacking in resolution since they are based on external phenology or microdissections [16,20,21,101–108]. More detailed studies have used Scanning Electron Microscopy (SEM), but even this technique is limited when it comes to observing internal cellular processes such as micro- and megasporogenesis [16,17,104,105,109]. Therefore, it was prudent to design a system for flower development from start to finish at high resolution. We anticipate this system will be applicable to other Prunus species since basic developmental patterns of individual flowers are often similar, although the structure of the inflorescence may vary [17,101,104,109]. Stage designations were based on previous literature and noticeable developmental landmarks/patterns that were visible from histological sections and hand dissections (Figure 4.3) [16,17,104,109–116]. Prior to dormancy, developmental stages for each tree were based on internal anatomy, as determined by histological sections. For buds harvested during dormancy 18 until bloom, staging was based on both histology and pistil measurements. Pistil measurements were recorded for every year of the study as a quantitative means to track floral bud development. Early- and late-blooming individuals differed in development as early as floral initiation Histological comparisons between bloom groups revealed the first detectable developmental differences in EB and LB siblings of this segregating population appear as early as floral initiation and persist throughout the year. In June 2018 and 2019, the apices from both bloom groups were vegetative, continually initiating bud scale primordia (Stage 1) (Figure 4.4: 6/12/19, 6/29/19; Figure S4.1). In August 2018, apices of the different bloom groups were well into organogenesis and EB trees were more advanced than their late siblings by an average of three stages (Figure S4.2). This necessitated a finer series of collections for the summer of 2019, which revealed that on average, EB individuals committed to floral development (Stage 2+) before their LB counterparts (Figure 4.4: 7/18/19). As organogenesis progressed into dormancy these slight disparities were maintained, and the LB individuals did not surpass their EB siblings in either year (Figure 4.5, Figures S4.2 – S4.5). Still, even as relative development was constant between the bloom groups, the same individuals observed on 8 Aug 2018, were more advanced compared to 8 Aug 2019 (Stage 5 – 9 vs Stage 4 – 7). 19 Figure 4.4 Early-bloomers are the first to make the vegetative to floral transition. Histological sections illustrate the stage of the apex from 6/12/19 – 7/18/19 based on the system in Figure 4.2. Sampling date is indicated on the y-axis. Individuals A – C are Early 1, Early 2, and Early 3. Individuals D – F are Late 1, Late 2, and Late 3. All apices appeared morphologically vegetative until 18 Jul, when more early-blooming apices on average had made the floral transition. These three dates and individuals of each bloom group were used for the RNA-sequencing analyses. Scale bars are 200 µm unless stated otherwise. As the flowers entered dormancy, EB trees were at similar stages in both years (Stages 12 – 13; Figure 4.5: 11/16/19; Figures S4.5 – S4.6). Sporogenous tissue was easily discernible from the other anther layers, and the protuberance from the placental wall marked the carpel margin meristem (CMM; Stage 12) and the start of ovule formation (Stage 13). Late-bloomers were also at similar stages between the two years (Stages 11 – 12), but with sporogenous tissue barely distinguishable from other anther layers and a rudimentary CMM at best (Figure 4.5: 11/16/19; Figures S4.5 – S4.6). Winter proceeded with little development until flowers prepared to exit dormancy (Figure 4.6; Figures S4.7 – S4.8). An upsurge in metabolic activity was marked by increasing amounts of elongated cells emerging from the CMM sometime between carpel closing and the appearance of stylar tracts (Figure 4.6B, E, H; Figure S4.7B, C; Supplementary Figure S4.8A – C). Given the timing of their appearance, it is speculated they may be involved in the development of stylar tissue. 20 Figure 4.5 Floral buds from early-blooming trees were more advanced as organogenesis proceeded to dormancy. Histological sections illustrate the stage of floral development from 8/8/19 – 11/16/19 based on the system in Figure 2. Sampling date is indicated on the y-axis. Individuals A – C are Early 1, Early 2, and Early 3, respectively. Individuals D – F are Late 1, Late 2, and Late 3, respectively. Note that only four individuals were sectioned for 11/16/19. Early 1 (A) on 8/22/19 is represented by two photos to show the variability in this individual. Black arrows in 11/16/19 indicate sporogenous tissue, which is more defined in the early- bloomers. On all these dates, early-bloomers keep their lead in development. Scale bars are 200 µm unless stated otherwise. Growth resumed exponentially in the spring and flowers underwent microsporogenesis before megasporogenesis (Figures 4.7 – 4.8; Figures S4.9 – S4.10). As seen at other time points, the EB trees persisted ahead of their LB siblings and bloomed first. On 6 Apr 2020, EB individuals continued pollen development and the megaspore mother cell became distinct from the surrounding tissues in the nucellus (Stages 15 – 17; Figure 4.7A – D). Meanwhile, LBs were just initiating microsporogenesis and integument primordia were rarely seen (Figure 4.7E – H). Occasionally, the stage of microsporocyte development was slightly out of sync compared to other features of the flower (Figure 4.7A, D; Figure 4.8G). For EBs on 2 May 2020, pollen was mature, the tapetum was fully degraded, and the megagametophyte was in the process of forming the embryo sac (Stages 20 – 21; Figure 4.8A – D). However, remnants of the tapetum endured in 21 the majority of LBs; and in all LBs, the megasporophyte ranged from just beginning megasporogenesis to just completing it (Stages 17 – 19; Figure 4.8E – H). Although megasporogenesis was not captured in the sectioning collection dates of the 2018-2019 season, assuming the EBs completed it prior to the LBs is a reasonable interpolation given the EBs bloomed first (Table S4.1). A graphical summary of the histological results is presented in Figure 4.9 for the seasons of 2018-2019 and 2019-2020. Figure 4.6 Floral buds from early-blooming trees remained ahead in their development when both bloom groups exited dormancy. Histological sections illustrate developmental stages on 2/22/2020, just before warm temperatures return. A – D represent Early 1, Early 2, Early 4, and Early 5, respectively, while E – H represent Late 1, Late 2, Late 3, and Late 4, respectively. Dotted lines indicate separate sections of the same flower to show the internal state of the ovary. The cell layers within the early-blooming anthers are more differentiated compared to those in late-blooming flowers, although sporogenous tissue is at least visible for all genotypes. The ovule is a noticeable protrusion in early-bloomers while most late-bloomer flowers have rudimentary carpel margin meristems. Scale bars are 200 µm unless stated otherwise. 22 Figure 4.7 Floral buds from early-blooming and late-blooming trees were at different stages of ovule formation and microsporogenesis in early spring. Histological sections illustrate developmental stages on 4/6/2020. A – D represent Early 1, Early 2, Early 4, and Early 5, respectively, while E – H represent Late 1, Late 2, Late 3, and Late 4, respectively. Arrows indicate the emerging megaspore mother cell at the center of the nucellus for several early- bloomers; this cell is not yet visible in late-bloomer flowers. Most late-bloomer flowers have not undergone microsporogenesis while most early-bloomer flowers have. Scale bars are 200 µm unless stated otherwise. 23 Figure 4.8 Floral buds from early-blooming and late-blooming individuals are at different stages of megasporogenesis in mid-spring. Histological sections illustrating developmental stages on 5/2/2020, the day before full bloom for Early 1. A – D represent Early 1, Early 2, Early 4, and Early 5, respectively, while E – H represent Late 1, Late 2, Late 3, and Late 4, respectively. Arrows indicate the megagametophyte, which ranges in development from just initiating meiosis (late-bloomer) to a fully-developed embryo sac (early-bloomer). Scale bars are 200 µm unless stated otherwise. 24 Figure 4.9 A summary of the histological results for early- and late-blooming siblings for two years show differences as early as floral initiation that persist until bloom. Qualitative developmental stage (see Figure 2) is plotted against calendar date. The lull in each sigmoidal curve represents dormancy, when very little development occurred and very few GDDs accumulated. Lastly, dissection photos in the bottom right illustrate floral development of a late- bloomer (top row) and early-bloomer (bottom row) on select dates; some of these dates are indicated on the graph above. Early-blooming siblings, those that inherited no k haplotypes, begin developmental divergence around the time of floral initiation and maintain this lead for the rest of the year. Scale bars are 200 µm unless stated otherwise. Pistil progression supports consistent differences between EB and LB individuals To quantify and compare development between EB and LB siblings, pistil growth was measured for the 2018-2019, 2019-2020, and 2021-2022 seasons (Supplementary Files 4.1 – 4.3). Since the pistil primordia does not emerge until late summer or fall, this limited possible measurements to this developmental window. Despite the variability in sample size between years, the pattern of development was consistent for each season (Figure 4.10). In all three years, small differences between bloom group prior to dormancy were detected, although they were statistically insignificant. Once growth resumed in the spring with the accumulation of GDDs, the rate of development was significantly faster for the EBs compared to the LBs. Boxed dates in Figure 4.10 emphasize when statistical significance was detected at the p < 0.05 level between the bloom groups according to a Tukey’s HSD test. Individuals reaching 0 on the log scale (representing an open 25 flower) with fewer GDDs bloomed earlier. Significance between the bloom groups was detected between 2711 – 2739 GDDs in 2019, 2564 – 2637 GDDs in 2020, and only at 3184 GDDs in 2022. The corresponding calendar dates for the three years were 4 May – 8 May, 24 Apr – 5 May, and 8 May, respectively. The reason for the shrunken range in bloom times and the spike in development across fewer collections in late spring of 2022 (~3100 – 3200 GDDs) was likely due to a string of unseasonably warm days that accelerated development. The quantitative pistil measurements are consistent with histological observations since they imply a small amount of variation between bloom groups was established prior to dormancy onset for all years. Although we did not have the power to identify statistical significance at these dates, it is possible these minor differences translated to the major (significant) disparities later in development. Indeed, since the growth rate in the spring is exponential, raising trivial differences to some power would result in greater absolute differences. 26 Figure 4.10 Pistil measurements of early- and late-blooming siblings for three years support small pre-dormancy differences and variation in spring growth rate. Pistil measurements are indicated on the left, with growing degree days (GDDs; base temperature 4.4 ºC) since January 1st of the year of floral initiation plotted on the x-axis and pistil completion (log scale) on the y- axis. Trees are colored by bloom group (blues = EB; pinks = LB) with means plotted as large symbols. Smaller symbols indicate distinct measurements for the respective individual. Slopes of growth trends for each bloom group are approximated from dormancy to the date the first tree had reached bloom. The inset in the top left demonstrates how pistils were measured. Measurements were done from 9/20/18 – 5/16/19, 10/13/19 – 5/15/20, 12/29/21 – 5/13/22. Gold boxes with asterisks indicate collections where a statistically significant difference between the bloom groups was detected (Type III ANOVA; p < 0.05). In all years, GDDs, bloom group, and their interaction significantly impacted developmental rate. Notice the different x-axis scale for 2021 – 2022. The photos on the right are representatives of Early 4 and Late 4’s floral buds on select dates of the 2019 – 2020 season. 27 Figure 4.10 (cont’d) Budbreak forcing experiments suggest contrasting temperature perception between early- and late-blooming siblings Classical forcing experiments were performed to assess the depth of dormancy and compare temperature perception of early- and late-blooming siblings. Dormant branch cuttings with floral buds were collected from individuals in each bloom group periodically for two winters. Each collection represented the amount of chill floral buds had sensed up to the date of collection. Budbreak on cuttings was recorded as the percent of flower buds reaching a certain stage on the BBCH scale, which is a decimal coding system for the phenology of woody perennials [101]. Stages 51 – 54 were modeled as a function of bloom group, heat (GDDs), chill (chill portions, CP, according to the Dynamic Model), and their interactions; inset photos in Figure 4.11 provide an example of these stages in sour cherry. The BBCH system was used since it is simple, non-destructive, and does not require flower dissections. Only results for stages 51 – 54 are reported as this was sufficient to describe flowers’ initial response to heat and progression to later stages of budbreak. Moreover, the health of the branches steeply declined during prolonged forcing (>12 – 14 days). Observations between 26.6 – 91.5 and 20.7 – 67.7 chill portions are reported for several EB and LB individuals in 2019 – 2020 and 2021 – 2022, respectively. EB and LB siblings were always harvested together on the same dates (Table S4.4). Chill was coded as an ordered factor while GDDs was considered a quantitative variable. Therefore, trends presented in Figure 4.11 demonstrate the bloom groups’ responses given a 28 certain level of chill (x-axis) and a constant, average amount of GDDs across time points for each year (121 GDDs for 2019 – 2020, 96 GDDs for 2021 – 2022). BBCH Stage 51 is characterized by swelling of the dormant bud and the first sign of green tissues beneath the bud scales. In 2019 – 2020, the main effects of bloom group, chill, and GDDs had a significant effect on whether flowers progressed to stage 51, with EB individuals having an average of 11% more flowers progressing to this stage compared to the LB trees (p < 0.05, Type III ANOVA and Tukey’s HSD; Figure 4.11). The interaction of GDDs and chill was also statistically significant. In 2021 – 2022, the main effects of bloom group, GDDs, their interaction, and the interaction between GDDs and chill were also significant. However, the main effect of chill was not (p = 0.30). After a given amount of heat, EBs averaged 21% more flowers reaching stage 51 compared to LBs. The average differences seen in both years are reflected in the growth trends: when flowers begin responding to forcing temperatures, the mean budbreak percentage for the early-bloomers is always slightly higher (Figure 4.11). This observation suggests EBs tend to respond to forcing before their LB siblings irrespective of how much chilling has accumulated. BBCH Stages 52 and 53 are defined by further bud swelling and exposure of green tissues. BBCH 52 was not a formal stage prior to the current study, so it was created here as an intermediate between 51 and 53 to add resolution between stages. In both years, the main effects of GDDs and its interaction with chill were major developmental drivers to these intermediate stages. However, the main effect of bloom group only made the significance cutoff in 2019 – 2020 for BBCH 52. For stage 53 of this year and both stages in 2021 – 2022, bloom group’s interaction with GDDs was statistically significant in driving budbreak. Chill was significant at the p < 0.05 level only in 2019 – 2020 as well. The reason for the discrepancies among the two years for the main effects of chill and bloom group may be due to increased variability as more individuals were added in the 2021 – 2022 season. Moreover, the design was less balanced since a few individuals were not sampled at every date (chill level; see Table S4.4). However, one can conclude the response trends agree for both stages and years (Figure 4.11). Lastly, BBCH Stage 54 is ‘budbreak’ in the most literal sense as it marks the first opening of the flower bud. In both seasons, the same effectors were statistically significant in influencing development to this stage: the main effects of bloom group and GDDs, the interaction of bloom group and GDDs, the interaction of bloom group and chill, and the 29 interaction of chill and GDDs. Effectively what this means is that given a certain amount of chilling, forcing translates to a faster rate of development for early-bloomers. Interestingly, these forcing results echo those of the pistil growth analyses. However, the analyses differ in their sensitivities: the more quantitative assessment of pistil growth suggested minor differences prior to dormancy, but development according to the BBCH scale placed the EB and LB groups at the same stage. This stresses the consequences of scale when it comes to developmental comparisons and their subsequent conclusions. Still, in traditional terms, one may conclude the EBs have lower chilling requirements than their LB siblings, and they perceive temperature differently. 30 Figure 4.11 Early-bloomers respond first to warm temperatures and additional chill equates to a faster rate of growth compared to the late-bloomers. Forcing experiments were conducted for two winters. The GDDs for each response are a constant average within years (121 GDDs for 2019 – 2020; 96 GDDs for 2021 – 2022) Each column represents the winter of one year, and rows represent flower percentages reaching stages 51, 52, 53, and 54. Stages were modeled separately to show the gradual budbreak responses of the bloom groups given the same amount of chill and forcing conditions. Gold boxes with asterisks indicate statistical differences between the bloom groups when the means were separated by chill level. Axes are consistent for visualization of the trends between years. Means for each bloom group are shown as solid blue squares or pink circles and the error bars are the 95% confidence intervals of the statistical model. Smaller symbols represent the means for each tree. Inset photos in the upper right corners of the 2019-2020 data are a representative photo for the respective BBCH stage. Note that since 2019 – 2020’s average GDD accumulation was slightly higher, this likely explains the slight increase in budbreak for all chill levels and stages compared to 2021 – 2022 (i.e., an average of 121 GDDs vs 96 GDDs). 31 Figure 4.11 (cont’d) Higher quantities of ROS associate with depth of dormancy irrespective of bloom group Reactive oxygen species (ROS) have gained traction as a potential growth regulator for dormancy transitions, with higher ROS content being associated with deeper dormancy [117]. Therefore, we were interested in whether the segregating individuals differed in their ROS content during dormancy. Hydrogen peroxide (H2O2), a relatively stable ROS, was quantified in floral buds collected during the late fall and winter of 2019 – 2020 using an adaptation of the xylenol orange assay (see Methods) [118]. Flowers were harvested from 13 Oct 2019 (6.6 CP since 1 Oct 2019) and 23 Mar 2020 (93.4 CP) for five dates (Table S4.5, Supplementary File 4.6). According to a type III ANOVA, bloom group did not significantly influence H2O2 concentrations (p = 0.15); however, chilling (synonymous with collection date) did (p = 4.9 x 10- 5; Figure 4.12). µmoles of H2O2 per gram of fresh weight were similarly lower on 13 Oct, 22 Feb, and 23 Mar compared to 16 Nov and 20 Dec. The higher concentrations of H2O2 coincided with deeper dormancy, as determined by the response in forcing conditions. The lower concentrations appeared to be linked with active metabolism and visible growth and / or with the speed at which flower buds respond to forcing temperatures. It should be noted that H2O2 showed 32 a weak negative correlation with temperature at the time of collection (Figure S4.11). However, this association should be interpreted with caution since this is to be expected if chill accumulation is equivalent to the time of year. Moreover, the range of temperatures during our collections was narrow; therefore, temperature was excluded from the statistical model as a more thorough assessment is needed to detect any sort of relationship. Nonetheless, these observations support the hypothesis that higher ROS content reinforces dormancy and is inhibitory to growth. Figure 4.12 Concentration of hydrogen peroxide (H2O2), a reactive oxygen species, associates with depth of dormancy. H2O2 content (y-axis) was measured five times from fall (13 Oct) to spring (23 Mar) to determine if there were differences between the bloom groups as they accumulated chill (x-axis) and navigated dormancy. While bloom group did not significantly affect H2O2 content, chill levels (collection date) did. The significantly higher concentrations (signified with ‘b’) represent early winter collections when flowers hardly responded to forcing conditions, which was determined by forcing experiments in 2019 – 2020. Means for each genotype are shown as blue squares or pink circles and the error bars are the 95% confidence intervals of the statistical model. 33 A genomic coverage analysis of EB and LB individuals suggests k is derived from Prunus cerasus subgenome A Recently there was evidence found that sour cherry as a species contains various dosages of three subgenomes (A, A’, and B) due to regular homoeologous exchange (personal communication with Iezzoni and Rhoades). Thus, we reasoned that the parents and individuals of the segregating population could have different dosages at the QTL locus, and exploring the subgenome content of this region may yield insight into which subgenome(s) the late-bloom k haplotype is derived from. Therefore, we sequenced 12 individuals, seven with 0k (EB) and five with 2k (LB), to a depth of approximately 5x per haplotype and aligned their DNA sequences to the three subgenomes of ‘Montmorency.’ All seven haplotypes identified at the chr4 QTL (k, a, h, n, o, i, and g) from the previous work [87] were included amongst these sequenced individuals. Therefore the haplotype composition for each EB sibling is (x, x, x, x), where x is any haplotype other than k. The composition for each LB sibling is (x, x, k, k). Knowing sour cherry has four alleles per site (i.e., it is tetraploid), we examined the relative coverage of each subgenome within the QTL region for these 12 individuals alongside their parents, ‘Balaton’ and ‘Surefire’ (Figure 4.13; Supplementary File 4.7). Through pairwise comparisons and deductive reasoning, we were able to determine that k originated from subgenome A, a subgenome donated by the late-blooming Prunus fruticosa-like progenitor. For example, by relative coverage analysis, Early 1 is determined to have dosage A’BBB at the QTL (Figure 4.13). The previous mapping study showed this individual to have a haplotype composition of (a, h, o, o). Because B is the only subgenome with at least 2X dosage, haplotype o must come from subgenome B. Furthermore, since there is no subgenome A dosage in the QTL region for this genotype, haplotypes a and h do not originate from subgenome A (Supplementary File 4.7). Notably, the sequenced late-blooming siblings showed no dosage of A’, meaning they did not inherit haplotypes a or g. This might partially be explained by the preferential pairing of subgenomes during meiosis, a well-documented occurrence in sour cherry [119,120]. For example, during meiosis in ‘Balaton,’ k may preferentially pair with a, and progeny would most likely inherit one or the other – although subgenome composition outside of the QTL should be accounted for when considering meiotic pairing. Similarly to k and a, n and h are not often inherited with one another (save for Early 7). Therefore, since progeny inheriting 2k were specifically selected for 34 comparisons in the present study, it would be expected these particular individuals would be unlikely to inherit an A’ haplotype. Figure 4.13 A comparison of the coverage of each subgenome at the bloom time QTL suggests the late-blooming haplotype k originates from subgenome A. Subgenome dosage at the bloom time QTL was determined by mapping gDNA Illumina reads of segregating individuals and their parents to the ‘Montmorency’ reference genome. Using haplotype information from Cai et al. 2018, pairwise comparisons of dosage and haplotype composition revealed the subgenome each haplotype most likely originates from (see Supplementary File 4.7). The bloom time QTL for the 12 siblings’ plots is left unshaded and separated by the rest of the chromosome with dotted lines. Subgenome dosages are given at the top right corner or top center (Early 7) of each plot. 35 Genetic variation relevant to the k haplotype After aligning gDNA reads of the individuals and parents to the three Montmorency subgenomes, single nucleotide polymorphism (SNPs) and small insertions and deletions (INDELs < 30 bp) were called in an ~8 Mb region encompassing the QTL – only calling variants in individuals with some dosage of subgenome A. Ultimately this meant comparisons were between haplotypes k, i, and the latter part of g. After hard-filtering for quality of the variant calls, the dataset was further filtered based on k’s uniqueness from the other subgenome A haplotypes. This analysis resulted in 8101 variants. According to snpEff, 7726 of these variants were classified as ‘modifiers,’ which contained upstream, downstream, 5’ and 3’ UTR, intragenic, intergenic, and intronic variants. 128 variants were predicted to have low effects, such as start codons gained and synonymous mutations. 222 were considered moderate effect variants, which includes nonsynonymous mutations and codon insertions and deletions. Finally, there were 25 high effect variants, which includes those that are predicted to cause a new splice site, stop codon, or frameshift in the protein product. Altogether, these variants were within or associated with 755 genes. The complete list of these smaller variants is provided in Supplementary File 4.8, which includes a separate tab for curated high effect variants (confirmed in reads viewed in a genome browser). We also attempted to identify larger structural variants (SVs; > 30bp) associated with the k haplotype. Reads mapping to subgenome A from individuals only containing haplotype i or haplotype k were pooled to increase accuracy of the calling. ‘Surefire’ and Early 7, both containing haplotype g, were used to rule out SVs whenever possible, but relatively low coverage limited their utility (see Methods). Variants were still kept if the k allele was homozygous and the other allele information was missing. Post-filtering, 169 variants were found, ranging in size from 31 to nearly 81000 bp (Supplementary File 4.8). Due to their length, some of these variants were predicted to have effects from those mentioned above on several different genes, which included the majority of genes within the QTL (n = 590). This included intergenic effects, which was predicted to affect 231 genes alone. To focus time and resources, only structural variants with effects predicted to be high effect or upstream of protein- coding genes were further considered. 98 of these variants were predicted to have at least one upstream effect and 18 were predicted to have high effects. 124 genes were anticipated to have 36 upstream SV effects, 28 were predicted to have high effect SVs, and 25 genes were predicted to be affected by both. Differentially expressed genes between EB and LB siblings To narrow the list of candidate genes for the late-blooming phenotype of haplotype k, gene expression differences between early- and late-blooming groups around the time of floral meristem initiation were investigated. RNA from three EB and three LB individuals’ apices was sequenced on 12 Jun 2019, 29 Jun 2019, and 18 Jul 2019. Early 1, Early 2, and Early 3 constituted the EB group while Late 1, Late 2, and Late 3 constituted the LB group. Based on the haplotype key shown in Figure 4.13, subgenome dosages for the EB group trees were A’BBB, AA’BB, and AA’BB, while dosages of the LB group trees were AABB, AAAB, and AAAB, respectively, within the QTL. Prunus cerasus subgenome A was used as a reference for alignment given the likelihood that it contains k. Since trees varied in their subgenome A dosages at the chr4 QTL, we relaxed parameters so that overall alignment rates to subgenome A were similar among individuals and time points (see Methods). We were also reassured gene expression quantification biases would be minimized since Prunus genomes are extremely syntenic [121]. For the first two time points, the apex in the floral position appeared vegetative; however, on 18 Jul 2019, more apices from the EB individuals began to proceed to stage 2 and 3 (Figure 4.4). A principal component analysis (PCA) revealed a clear separation between bloom groups, demonstrating these individuals were already distinct at the molecular level on 12 June and 29 June despite their similar morphology (Figure 4.14; Figure 4.4). Markedly, Early 1 was substantially removed from its fellow EBs – suggesting unique aspects of subgenome A, especially at the chr4 bloom time locus, drive developmental changes within these tissues. As bloom time is a quantitative trait, even subtle changes in gene expression could be responsible for the phenotypic differences; therefore, genes with an absolute log2 fold change of 0.5 or higher and a p-value less than 0.05 were considered differentially expressed genes (DEGs; Supplementary File 4.9). While the results here concentrate on intercomparisons between the bloom groups, DEGs and associated GO terms for intra-comparisons (e.g., 12 June early vs 29 June early) were also recorded (Supplementary Files 4.9 – 4.10). 37 Figure 4.14 Principal Component Analysis (PCA) for three early-blooming and three late- blooming individuals suggests fundamental molecular differences between bloom groups. The PCA was produced with RNAseq data from 12 June 2019, 29 June 2019, and 18 July 2019 of EB and LB siblings’ apices in the floral position. PC1 appears to contain most of the variation seen from the different sampling dates while PC2 is at least partially driven by bloom group. RNAseq reads were mapped to subgenome A and two outliers were removed prior to producing the PCA. On 12 Jun 2019, 1371 genes were differentially expressed (Figure 4.15A). Approximately 13% (n = 183) of the DE genes for this date were expressed within the chromosome 4 QTL, despite this region only being approximately 3% of the genome (Figure 4.15B). The top 20 gene ontology (GO) terms enriched in DEGs more highly expressed in EBs (blue) and LBs (pink) are given in Figure 4.16. These top 20 terms support the notion that EBs were undergoing developmental processes mediated by a suite of hormones such as jasmonic acid (JA), salicylic acid (SA), abscisic acid (ABA) and gibberellic acid (GA). The GO terms of genes more highly expressed in the LBs also included JA-related terms, as well as brassinosteroid (BR) and ethylene associated terms (Supplementary File 4.10). On 29 Jun 2019, 879 genes were found to be differentially expressed, and nearly 15% of these DE genes (n = 129) were located within the bloom time QTL on chromosome 4 (Figure 4.15A – B). The top 20 GO terms for genes expressed at higher levels in the EBs implicate 38 auxin- (IAA) and GA-mediated pathways at this stage of development (Figure 4.16), and a venture to the top 50 connects SA and ABA regulated processes. For genes more highly expressed in LBs, GO terms appear to be more related to basic metabolic and growth-related activities: for example, ‘lignin biosynthetic process,’ ‘regulation of cell death,’ and ‘protoderm histogenesis,’ might describe high rates of cell turnover. On 18 Jul 2019, EB individuals showed a higher proportion of apices making the transition from a vegetative to reproductive meristem. 1379 DEGs were identified on this date, with 185 of these housed within the chromosome 4 QTL (Figure 4.15A – B). The genes more highly expressed in the EBs support the morphological observations, with ‘inflorescence morphogenesis’ and ‘radial axis specification’ among the top 20 GO terms (Figure 4.16). For genes more highly expressed in the LBs, several of the top 20 terms are repeated from 29 June (e.g., ‘regulation of cell death,’ ‘protein autophosphorylation,’ and ‘callose localization’), perhaps suggesting fewer major developmental changes since 29 June. Markedly, 366 DEGs were common to all three dates, with 93 (25%) of these located in the chromosome 4 QTL (Figure 4.15A – C). 39 Figure 4.15 Venn diagrams and a heatmap illustrating the number of differentially expressed genes between early-blooming and late-blooming trees by date in the full genome and chromosome 4 QTL. A) Number of genes expressed by date throughout the genome, using subgenome A as a reference. B) Number of genes expressed by date only within the chromosome 4 QTL. C) A heatmap showing relative fold changes of DEGs at every date within the QTL. Blue shading represents genes more highly expressed in late-blooming individuals, while red shading represents genes more highly expressed in the early-blooming individuals. Gene numbers beginning in 097 were manually edited with Apollo (see Methods) and thus do not reflect relative order on the chromosome. 40 Figure 4.16 Results of the Gene Ontology enrichment analysis for each date’s comparison. The top 20 most enriched GO terms of genes more highly expressed in early-bloomers (blue) or late- bloomers (pink) are shown. The bolded number above each bar is the number of genes in the DEGs list with the respective GO annotation, while the number below is the amount expected based on the number annotated in the full database (DESeq2 object). The log value of the inverse p-value is plotted on the y-axis – thus, the taller the bar, the more significant the p-value according to a Fisher’s Exact Test. All p-values are below 0.05. 41 Uniting genetic variation with differentially expressed genes To obtain quality candidate genes for the basis of k, DEGs information was intersected with genetic variation. For smaller sources of variation (SNPs, INDELs < 30bp), we chose to focus on variation categorized as high effect (stop codon, frame shift, splice site changes), moderate effect (nonsynonymous amino acid changes, codon insertions/deletions), and certain modifier effects (upstream, 5’ UTR, and 3’ UTR) since these types would presumably affect protein function or gene regulation. For larger variants, only upstream and high effect structural variants were considered. In total, 186 DEGs (for any date) in the QTL region were predicted to have at least one of the above listed kinds of variation (Supplementary File 4.11). 130 of these were expressed on 12 June 2019, 89 on 29 June 2019, 133 on 18 July 2019, and 62 were differentially expressed at all three dates. Weighted Gene Correlation Network Analysis Next, we explored co-expression networks to identify groups of genes with expression profiles that distinguish the two bloom groups from one another. After filtering for genes with a coefficient of variation of at least 0.2 among the samples and an FPKM (Fragments Per Kilobase of transcript per Million mapped reads) greater than five, 2053 genes were grouped into 19 modules ranging in size from 33 to 400 genes (excludes the grey module as these genes by definition do not fit well into any module; Figure S4.12; Supplementary File 4.12). To visualize the expression profiles for each of these modules, the module eigengene values were plotted for each sample across the three time points (Figure 4.17). Several modules exhibited interesting trends relevant to the bloom time phenotypes, and the top 20 GO terms associated with these are provided in Figure 4.18. According to the GO enrichment analysis, the hotpink module (n = 85) was suggestive of terpene-based hormone synthesis, with genes in this module showing higher average expression in EBs compared to LBs on every date. The gold module (n = 92) showed similar biological terms but exhibited distinctions regarding abiotic stress (e.g., cellular response to water deprivation) and possibly hormone-regulated processes as it was enriched for strigolactone biosynthetic processes. It also displayed a similar expression pattern to that of hotpink with the exception of Late 1, whose expression profile was far less removed from the EBs compared to its fellow LBs. The darkgreen module (n = 102), which showed lower expression on 12 Jun compared to 29 Jun and 18 Jul for all samples, was enriched for terms related to auxin signaling and specialized tissue development. For example, suberin and lignin 42 are likely components of the bud scales surrounding the apex (Figure 4.18). At any given date, EBs presented greater average expression of the genes in this module compared to LBs. Finally, LBs showed slightly higher expression of genes within the lavender module (n = 216), although Early 1 is somewhat of an outlier. This module was enriched for terms related to cell turnover, most notably ‘homeostasis of number of meristem cells.’ In agreement with the enriched GO terms, the hotpink module includes seven terpene synthases, a cytochrome P450 protein, and two AP2/ERF transcription factors. 15 genes in this module have some form of genetic variation, are differentially expressed, and are within the chr4 QTL. The gold module includes two MYB transcription factors, an ortholog of the transcription factor DEMETER, two cytochrome P450 proteins, and four terpene synthases. The darkgreen contained genes including two MADS-box containing genes, three auxin-responsive proteins, and four cytochrome P450 proteins. For the lavender module genes, seven INSENSITIVE 1-ASSOCIATED RECEPTOR KINASES and a FANTASTIC FOUR-like protein were identified. The latter has been shown to regulate meristem size in arabidopsis [122]. 43 Figure 4.17 Module Eigengene (ME) expression trends for each module identified in the weighted gene correlation network analysis. MEs represent groups of genes with similar expression patterns and can help describe the differences in individuals. For example, the gene expression of the seashell module sets Late 1 apart from all individuals, while lightcoral sets Late 2 apart. The dotted gold circles draw attention to the hotpink, gold, darkgreen, and lavender modules, as all four separate the bloom groups to some degree. The number of genes housed in each module is given in parentheses next to the heading of each plot. A total of 2053 genes throughout the full genome were grouped into MEs. The hotpink module, by far, contains genes that most effectively separate the EB and LB siblings’ expression profiles. 44 Figure 4.18 Results of the Gene Ontology enrichment analysis for select modules from the network analysis. The top 20 most enriched GO terms of genes found in the hotpink, darkgreen, gold, and lavender modules are shown. The bolded number above each bar is the number of genes in this module with the respective GO annotation, while the number below is the amount expected based on the number annotated in the full database (DESeq2 object). The log value of the inverse p-value is plotted on the y-axis – thus, the taller the bar, the more significant the p- value according to a Fisher’s Exact Test. All p-values are below 0.05. Candidate genes for the basis of k Gene candidates for the late bloom phenotype were required to be located within the chr4 QTL. They were then prioritized based on the results of the differential expression analysis, their amount and type of genetic variation, their co-expression with genes of certain modules (which distinguished early- and late-blooming individuals transcriptionally; Figure 4.17), and their predicted gene function (Table 4.1). According to the PCA, bloom groups were separable transcriptionally by 12 June despite all apices being vegetative, suggesting the k haplotype may impact phenology even before floral initiation. Thus we reasoned the greatest chance of capturing the activity of gene(s) behind the basis of k would be on the earliest date, 12 June. Genes that were differentially expressed on at least 12 June (i.e., expressed only on 12 June, 12 45 June + [29 June or 18 July] or at all time points) were favored. Since the hotpink module of the network most cleanly delineated the bloom groups, 8 of the 13 highest ranked candidates were co-expressed with genes in this module. The fine-mapping results for a chromosome 4 bloom time QTL discovered in a population developed from P. avium cultivars ‘Regina’ and ‘Garnet’ were also taken into account when considering candidates [45]. This sweet cherry QTL was associated with a delay in bloom time of nearly 4 days [88]. The peak position of the chr4 sour cherry QTL was barely outside the 95% confidence interval of the one identified the sweet cherry, but since their confidence intervals overlapped, they were previously considered the same QTL [87]. We therefore identified the syntenic region of the fine-mapped sweet cherry QTL on sour cherry chromosome 4A (approx. Pcer_chr4A: 9717998 - 9802438) and surveyed genes within this region for differential expression for any dates sampled in the current study. However, there were none – suggesting these QTL are distinct from one another, and the effect of k may be due to separate gene(s). Among the topmost ranked candidates are several terpene synthases and B3 domain- containing transcription factors. The terpene synthases are of considerable interest since hormones known to be involved in several phenological and developmental transitions are terpene-based (e.g., GA, ABA) [20–22,123]. One terpene synthase (Pcer_097523-RA) has a single base pair deletion anticipated to induce a frame-shift in the k allele – consistent with this prediction, the expression has been reduced to nearly nothing in the late-blooming individuals. Several of the B3 domain-containing transcription factors share homology with arabidopsis AT5G66980, which has a strikingly floral-specific expression pattern and belongs to a B3 domain-containing subfamily of proteins called REPRODUCTIVE MERISTEM (REMs) [124]. Like other members of this subfamily, these three candidates exist in a tandem array of six genes within an 18 kb region (Pcer_chr4A: 13481626 - 13500091). However, the predicted protein sequence of the fourth reads through one stop codon and was initially missed during annotation (Pcer_024686-RA). 46 Table 4.1 High-quality candidate genes for the basis of the k late-bloom phenotype. These contain genetic variation and are in most cases expressed in modules of interest. Genes beginning with 097 were manually edited with Apollo and thus are not indicative of relative chromosome position. 47 Discussion Various studies have tried to identify connections between early floral development and final bloom time in fruit trees. In 1984 and 1985, Raseira and Moore grouped nine peach cultivars into low (< 500 chilling hours below 7.2 ºC), medium (550 – 750), and high (800 – 1050) chill requirements and assessed floral bud initiation times. Between the two years, it was extremely variable: in 1984, the low-chill group was the first to initiate floral buds, whereas the high-chill group initiated theirs first in 1985 [29]. In 1989, Hoover et al. tracked floral development in four Malus × domestica cultivars and found cultivar-specific rates of organogenesis, but no clear association between initiation and bloom time [30]. Fadón et al. investigated early floral development in five sweet cherry (Prunus avium) cultivars. They noted there was asynchronous progression among the cultivars for all three years, but that a common stage was reached by dormancy – albeit at different times [102]. More recently, Goeckeritz et al. described the pre-dormancy development of six apple accessions – four species total – and showed a wide range of stages at approximate dormancy [108]. Aside from other limitations of each study – including the scale at which flowers were examined (e.g., cellular level vs macrodissections) – each work fails to control for (Goeckeritz et al.) or even address (Raseira and Moore, Hoover et al., Fadon et al.) the wide genetic basis behind comparisons. This genetic variation likely clouds our ability to identify subtle changes in phenology between genotypes, which could explain why an association for floral initiation and final bloom time has yet to be found. Many genes affect bloom time, so it is quite probable that several mechanisms are confounded when contrasting diverse materials. But confronting this issue is easier said than done. Research orchards are expensive to maintain, so access to segregating populations is rare. To make matters more complicated, ‘domesticated’ fruit trees are sometimes not so far genetically removed from the wild populations in which they originate – meaning the pedigrees of some cultivars are nigh impossible to disentangle. Still, we cannot ignore that breakthroughs in gene discovery, including those involved in bloom time regulation in rosaceous fruit trees, often come from classical mapping approaches [45,73,84]. Therefore, it is in our best interest to devote the time and resources necessary to control for genetic background in future studies and acknowledge the constraints of our interpretations. The work detailed here takes advantage of the close relations in a population of sour cherry created at the Michigan State University Clarksville Research Center. We built on 48 previous work where a major QTL was located on linkage group 4, a hotspot for bloom time regulation in Prunus [87]. In other studies that searched for candidate genes within Prunus chromosome 4 QTLs, researchers automatically geared their focus toward dormancy, which is indisputably critical to bloom time in Rosaceous fruit trees [45,91,125]. However, doing so fails to address the developmental events that have occurred up until that point. An analogy for an annual species that grows and reproduces in a single season might be to attribute bloom time variation only to the developmental window from the first visible flower to anthesis. Our morphological and quantitative data show EBs surpass their LB siblings as soon as floral development begins, and that their lead is maintained even up until bloom the following spring. This cannot simply be the byproduct of these individuals blooming earlier the previous flowering cycle, and getting a head start on growth in the new season. If this were the case, the separation between these bloom groups would not be as consistent across years. The purpose of extensively dissecting this population’s flowering trajectory was to increase our resolution of the events altered and provide clues for candidate gene function and selection. Furthermore, assessing molecular distinctions as close to the time of divergence as possible increases the chances of capturing expression of the gene(s) driving the phenotypic differences while controlling for confounding effects of changes that may have already occurred. RNA-sequencing data indicated that, despite their morphological similarities, gene expression profiles were fundamentally different between the bloom groups at all dates sampled – suggesting the effect of k may act well before the apex has committed to the floral transition. Based on the gene expression data and subsequent GO term analyses, the EB and LB individuals fine tune their hormone balances during the vegetative to reproductive phase transition, which is in accordance with observations in sweet cherry and apple [18,20]. Examination of the GO terms found in upregulated DEGs as EBs proceeded from 12 June to 29 June included ‘response to ABA’ and ‘response to IAA’ paired with ‘vegetative to reproductive phase transition,’ lending support to the positive influences ABA and IAA may have during this transition in fruit trees [20]. The LBs show a strikingly diverse profile – in addition to ABA and IAA, they show enrichment terms related to BR, JA, SA, CK, and GA, the latter of which is robustly connected to maintenance of the vegetative phase in fruit trees [22–24]. Of course, without directly quantifying these hormones, we could speculate endlessly on their ratios and crosstalk. Still, from a broad perspective of these data, especially when including GO 49 information for co-expression modules, a consistent feature seems to separate the bloom groups: terpene synthesis. Terpene compounds are greatly diversified in plants and take part in everything from biotic defense responses to phytohormone synthesis [126]. GA, ABA, and BR all have isoprene units common to their synthesis, and isoprenoids are the building blocks of terpenoids [127]. Gibberellins constitute a group of tetracyclic diterpenes (C20), ABA is synthesized from a carotenoid precursor and classified as a tetraterpenoid (C40), while BR is made from the triterpene (C30) squalene [123]. Terpene synthases, the enzymes responsible for these chemical reactions, may form species-specific paralagous gene clusters with tissue- specialized gene expression [124,126,128]. The hotpink module identified in the network analysis (n = 85) most effectively separated the bloom groups and contained seven terpene synthases as well as several hormone-related GO terms (Figure 4.18). Variable GA responses and/or synthesis as the basis of k is an interesting theory to entertain in light of these individuals’ phenotypic differences. Not only is GA linked to the floral transition, but its relationship with ABA is also a key indicator of dormancy depth [33]. In the two seasons of forcing experiments, it was shown that EBs responded to warm temperatures first; this implies a shallower state of dormancy and the possibility that k may affect temperature responses as well. In 2019 – 2020, we attempted to saturate the chilling response for both bloom groups by placing branches in a 5 ºC chamber before forcing (CP = 91.5). However, as shown in Figure 4.11, this only served to further separate these individuals as additional chill caused faster growth rates for the EBs compared to LBs. This trend is most striking at later stages of budbreak (BBCH 54), which is also observed in the gradual divergence of the pistil’s growth (Figure 4.10). The apparent differences in temperature-sensing mechanisms among the bloom groups fits well with GA’s known involvement in these pathways [129]. It would be intriguing to examine gene expression at these stages of development to ascertain whether certain genes in the chr4 QTL are differentially expressed between the bloom groups at these times as well. GA-modulated flowering phenology is also intriguing given sour cherry’s evolutionary history, since k shows dosage effects and originates from a subgenome donated by the cold-adapted, late-blooming P. fruticosa-like ancestor [5,7,9]. The apparent paradox that higher GA levels negatively regulate the floral transition in rosaceous fruit trees but positively regulate dormancy transitions and spring growth could be explained by their spatiotemporal context. For example, in arabidopsis, GA signaling first encourages the floral transition but then 50 later prevents flower formation in the lateral primordia of the inflorescence [27]. Taken together, these results advise a closer examination of terpene-based hormone synthesis as a potential means of regulating bloom time in sour cherry. Also among the promising candidate genes in Table 4.1 are four differentially expressed B3 domain-containing transcription factors, two of which are co-expressed with genes in the hotpink module. AUXIN RESPONSE FACTORS (ARFs), ABSCISIC ACID INSENSITIVE 3 (ABI3), HIGH LEVEL EXPRESSION OF SUGAR INDUCIBLE (HSI), RELATED TO ABI3/VP1 (RAV), and REPRODUCTIVE MERISTEM (REM) genes make up five subfamilies within the B3 domain-containing superfamily [130]. The candidates noted here all share closest homology with the REM subfamily, which coincidentally has the most sequence divergence among this family and the highest rate of expansion [130,131]. Usually, they exist in a tandem array, with REM gene clusters showing more sequence similarity within themselves than between clusters [130]. As their name implies, a survey of the expression patterns in arabidopsis shows they are tissue- specific – at times, only expressed in a few cells – with various roles in flower development [124,132]. However, due to their duplicative nature, functional redundancy and tight genetic linkage has made identification of specific phenotypes difficult [124]. The best functionally- characterized member of the REM subfamily is undoubtedly VERNALIZATION 1 (VRN1), a gene required for the silencing of FLC and the accelerated flowering time seen after exposing arabidopsis to low temperatures for several weeks [93,133]. Certain residues are crucial to VRN1’s ability to bind DNA – but like the rest of its subfamily members, there seems to be no clear motif it binds to [134]. More recently, REM17 was shown to be a positive effector of flowering time that is negatively targeted by FLC and SVP1 and positively targeted by SOC1 after vernalization [135]. Furthermore, REM16 was identified as a flowering time regulator through its interaction with ACTIN DEPOLYMERIZING FACTOR 1 (ADF1); REM16 was shown to directly bind to both the SOC1 and FT promoters to enhance their expression [136]. Evidently, REM proteins have essential functions in regulating flowering time. The B3 domain- containing/REM proteins identified as candidates for the late-blooming phenotype in sour cherry are therefore quite promising; but certain strategies such as overexpression lines, RNA interference, and/or CRISPR knockouts of genes in a tandem array may be necessary to conclude their functions. 51 As discussed earlier, there is evidence genes regulating bloom time in one Prunus species may regulate bloom time in another [83,137]. Therefore we compared the expression of candidate genes for an overlapping bloom time QTL in sweet cherry to the sour cherry chr4 QTL DEGs identified in the present work – however, the syntelogs showed no clear expression patterns between our bloom groups. It is possible we were unable to capture the expression of gene(s) conferring the late bloom time seen in this sour cherry population; but if we take the number DEGs expressed in the QTL as a proxy for flowering time regulators, then this region is especially enriched. Therefore, it is reasonable to suspect more than one major regulator of bloom time is present on Prunus chromosome 4. Originally it was thought these regulators may be the same since extant sweet cherry and sour cherry share a very recent common ancestor [5,9]. However, it was determined here that k is derived from subgenome A, a P. fruticosa-like subgenome. Given the contrasting selective pressures in the native habitats of P. fruticosa and P. avium, we now have incentive to believe the bloom time variation in these species is conferred by separate genes. In addition to k, other alleles in sour cherry have unique potential to additively delay bloom in other Prunus species since some accessions of sour cherry are the latest bloomers of the entire genus (Iezzoni personal communication) [87]. It will be further thought-provoking to observe what effects these alleles may have on the morphological progression of the flower in Prunus species with solitary (e.g., peach, plum) versus corymbose (e.g., sour cherry, sweet cherry) inflorescence structures. In conclusion, exciting new prospects for curating bloom time in Prunus lay ahead. These findings are projected to assist the industry both culturally and genetically. For example, sour cherry may benefit from GA synthesis inhibitor applications during deep dormancy to slow down future budbreak. And, since the late bloom gene(s) controlling k is probably unique from that found in sweet cherry, there may be more raw genetic material with large effects on Prunus bloom time than originally thought. 52 Future directions for the present chapter The next logical direction after a study like this is functional validation of the most attractive candidate genes. Functional validation via transformation into arabidopsis is suitable for the short-term – but alongside this goal it would be worthwhile to make plans to functionally validate these genes in a fruit tree as well. Transformations have been successful in apple; and for reassurance this will be fruitful in Malus, look to the DAM genes. DAM1 orthologs have been transformed in apple and show expected dormancy-related phenotypes [138]. Prunus would be ideal, but it likely would involve more troubleshooting and polyploid genetics. I believe the genes within the hotpink module should be prioritized first. From the discussion, it is apparent which candidates I find especially intriguing: the terpene synthases and B3 domain-containing proteins. For the former, I think it would be best to directly transform the cherry sequence(s) (including their promoters) into the model species of choice. Terpene synthases are diverse enough that obtaining SALK lines or other insertion mutants of homologs may not be a fair functional assessment of the protein in question. The specific expression patterns of candidates should also be carefully considered. For example, Pcer_097523-RA is essentially not expressed at all in the LB trees. This may be unusual for what we would expect for k. The EB siblings show expression for this gene, presumably from other alleles/subgenomes, most of which they share. Thus it would be expected to show some expression if the knockout was truly unique to k – unless this knockout was somehow affecting the protein product from other alleles. However, that would imply a dominant genetic effect of k, and we already know it to be additive [87]. In addition to making overexpression (OE) and RNAi / CRISPR constructs with the genes of interest, it would also be telling to know the spatial expression patterns of candidates. This may be relatively simple using a GUS reporter system driven by the cherry promoters. For subcellular location, a fluorescent tag attached to the protein would be required, such as GFP. If the hypothesis is centered on hormone synthesis, it would be useful to quantify hormones in mutants, especially terpene-based ones (ABA, BR, GA); and, depending on the construct, focus efforts on the tissue locations identified using the reporter systems. Furthermore, treatments with bioactive terpene-based hormones while monitoring the expression of candidate gene(s) (with qPCR) in wild-type might reveal regulatory mechanisms. The B3 domain-containing proteins are worth investigation as well, perhaps even more so than the terpene synthases. And it is worth noting that these two groups of candidates are co- 53 expressed in the hotpink module. This suggests they both may be regulating and/or be regulated by various hormone pathways as well. Although circumstantial, I also find it intriguing that the B3 domain-containing proteins are in close linkage with a marker that co-localizes with the k phenotype: RosBREED_snp_tart_cherry_f_Pp4_10832168. The best mapping for this marker (90.2% sequence identity) on chromosome 4A is at 13840124 bp. The tandem array of B3 domain-containing genes is located approximately from 13481626 to 13500091. To investigate the functions of the B3 domain-containing transcription factors, many of the above-mentioned experiments would apply here as well. On the applied side of things, it would be worth experimenting with certain plant growth regulators that may stall development at certain stages and cause later bloom. Start with dormancy, since the industry has a successful history of this already with minimal deleterious effects (e.g., Dormex®). Things to carefully consider would be the mode of action, concentration, duration of activity, health and environmental effects, and expense. Moreover, timing and the number of applications will be critical. The most biologically accurate understanding we have of chilling is through the Dynamic Model, which describes a theoretical product that builds to eventually become a chilling portion. Once accumulated, chilling portions are irreversible [50,51]. Therefore, to effectively increase the length of dormancy, repeated applications of a longer-acting substance will probably be needed, and it should be applied regularly during earlier phases of dormancy. Regulators increasing active ABA concentrations and/or decreasing active GAs or their biosynthesis may be promising candidates. Finally, I think the field could benefit immensely through basic experiments aimed at tracking hormone localization and transportation as the flower develops. Doing so will inform future decisions on when PGRs are applied and perhaps make them more targeted. Fluorescent biosensors for certain hormones may be worth exploring [139]. An alternative might be RNA in situ hybridization studies targeting rate-limiting enzymes of hormonal pathways; however, these would provide limited information on hormone transportation. Other important directions for the field I would suggest could be interpreted more politically than scientifically. Mainly, breeding programs ought to be given more funding and resources not only to develop elite germplasm, but to create segregating populations for gene discovery. This is particularly significant for fruit trees since the long juvenility phases pose an extra challenge. 54 Methods Sampling procedures for morphological, physiological, gDNA-sequencing, and RNA-sequencing analyses All materials were collected from trees at Michigan State University’s (MSU) Clarksville Research Center (42.87361004512462, -85.25873447306927). The number of individuals included in our observational experiments from year to year largely depended on factors such as tree health, size, and available time and resources. Tier, row, and tree numbers are given in Table S4.1, alongside their identity in the present work. Whole apices/flower buds for tissue fixation and sectioning, pistil growth measurements, and ROS extraction were sampled randomly on the east and west sides of the trees (the orchards run north to south). Branches for forcing experiments were also collected randomly but with the additional goal of having at least 25 flower buds total per genotype per collection date from at least two separate branches. We had high confidence collected apices/buds were or would become flowers since Prunus species tend to flower and fruit in consistent positions from year to year [16]. Replicates intended for ROS quantification were a pool of four adjacent flower buds. For RNA-sequencing, three biological replicates for six trees (3 early-blooming and 3 late-blooming) were sampled randomly as a pool of 6 neighboring apices at three dates (12 June 2019, 29 June 2019, 18 July 2019). Sampling occurred at approximately the same time of day for each collection (12 – 2pm) to mitigate gene expression noise between time points. Whole apices were excised with a razor blade, immediately flash-frozen in liquid nitrogen in the field, then stored at -80°C until RNA extraction. Dates for pistil and histological analyses, branches for chill assessments, flower buds for ROS (H2O2) quantification, and tissue collection for RNA-sequencing are given in Table S4.2 – S4.5, with corresponding growing degree days (GDDs) and Dynamic Model chill portions (where applicable). For the present study, growing degree hours / days (GDHs / GDDs) were measured using the method developed by Richardson et al. (1975): GDDscumulative = sum((hourlytemp—basetemp)/24) with a base temperature of 4.4 °C and no max temperature limit [48]. Chill portions were calculated using chillR [49]. All hourly temperature data was downloaded from the Clarksville, MI, MSU weather station: https://legacy.enviroweather.msu.edu/weather.php?stn=clr [140]. 55 Tissue fixation, embedding, histological sectioning, and staging Fixation, embedding and sectioning were done using the protocol described by Goeckeritz et al 2023 [108]. Briefly, flowers were collected in the field and either placed in a 2.5% glutaraldehyde / 2.5% formaldehyde / 0.1 M sodium cacodylate fixative solution (pH 7.2) after excising the tip of the bud to maximize infiltration, or placed in fixative after bud scales were removed in the lab. They were then subject to at least two 30 min vacuum infiltrations before being stored at 4 °C until dehydration and embedding. Sectioning was done at a 1 µm thickness using a PMC XL Ultramicrotome (Boeckeler Instruments, Tucson, AZ, USA), and sections were heat-mounted to glass slides and stained with a diluted toluidine blue solution (catalog no. 14950, Electron Microscopy Sciences). The staging system for flower development in Prunus cerasus was created by studying histological patterns in development among individuals throughout two seasons (2018-19; 2019- 20) and a review of the literature [115,116]. Features such as the emergence of organ primordia (sepals, petals, stamens, pistil) and the emergence of specific cell layers were used to demarcate staging. Particular attention was paid to the male (anthers) and female (ovule) gametophytes since their development was more easily defined by discrete stages; in contrast, sepal and petal advancement mainly manifested as relatively subtle cell divisions and elongation. SEM photos of related species’ development (P. laurocerasus and P. serotina) assisted with 3-dimensional visualization of the pistil/carpel and was critical to determining the approximate timing of closure [109]. Illustrations in Esau 1965, Bierhorst 1971, Foster and Gifford 1974, Mauseth 1988, and Fahn 1990 were fundamental in accurately recognizing features of micro- and megagametogenesis; for example, anther cell layers and the general progression of the embryo sac. Stages were hand drawn using an iPad Air (5th generation) and the ProCreate version 5.3.5 software. Pistil measurements and analysis Pistil measurements for trees were done for years 2018-2019 (three trees per bloom group), 2019-2020 (six EB, four LB), and 2021-2022 (seven EB, five LB) as growing degree days accumulated. 5 – 10 flowers from separate buds per tree were dissected using a Nikon SMZ800N stereo microscope and imaged with a Nikon DS-Fi3 color camera and Nikon NIS- Elements BR 4.60.00 software (Nikon, Tokyo, Japan). Pistils were measured from the tip of the stigma to the ovary base (see Figure 4.9 inset). On the day of bloom (approximately 50% open 56 flowers on the tree), 30 open flowers were collected for each genotype to normalize measurements according to the final pistil length of each individual; thus, the data is expressed as fraction of final pistil length. Due to the different numbers of individuals and collection times per year, each year was modeled separately as: Log(pistil completion fraction) ~ Group+GDDs+Group*GDDs+(1│Genotype:Group)+ (1│Genotype:Group:GDDs)+error A logarithmic transformation was necessary to approach statistical normality. ‘Group’ was a fixed effect factor with early- and late-blooming levels and ‘GDDs’ was a fixed effect, ordered factor (to detect significant differences at specific dates) with as many levels as collection time points that year. ‘Group * GDDs’ represented the developmental effect of the interaction of bloom group with accumulating heat. The ‘(1 | Genotype:Group)’ represented the random effects that each individual tree, which was nested within bloom group, contributed to development. Lastly, the ‘(1 | Genotype:Group:GDDs)’ term represented the random effect on development that each tree had at each collection time point. The raw data (pistil measurements / pistil completion fraction) is provided as Supplementary Files 4.1 – 4.3 and R scripts used to analyze the data and create the figures are available at https://github.com/goeckeritz. Forcing experiments and analysis Forcing experiments took place in the dormant seasons of 2019-2020 and 2021-2022. Branches for forcing were collected from six EB and three LB individuals in 2019-2020 at eight dates (26.6 – 91.5 CP), and from 6 – 7 EB and 3 – 4 LB individuals in 2021-2022 at nine dates (20.7 – 67.7 CP). Budbreak scoring was done using a slightly modified BBCH scale, a decimal coding system describing the phenology of woody perennials [101]. The dependent variable for each statistical model was the proportion of flowers reaching BBCH stages 51, 52, 53, and 54 (i.e., four statistical models were created per year); data after stage 54 is not presented since the health of the branches seemingly declined after this developmental stage. The raw formatted data is provided as Supplementary Files 4.4 – 4.5. The ends of each collected branch were given fresh cuts, placed in distilled (DI) water, and forced out of dormancy in a chamber held at a constant a 22°C (+/- less than 1 °C) with fluorescent lighting set to a 12-hour photoperiod. Branches received fresh cuts and DI water thereafter every 1 – 2 days, when scoring of each flower bud occurred. Briefly, stage 50 was described as dormant with no visible bud swelling. Stage 51 was characterized by first swelling, with some evident green tissue emerging from the 57 breaking bud. Stage 52 was added to the scoring system as an intermediate between stages 51 and 53. A flower bud was considered stage 53 when green tissue could be seen at least halfway down opposite sides of the bud. Stage 54 was reached when the bud cracked open due to the enlargement of the flowers inside. If any field GDDs had accumulated after January 1st in the year of study, these GDDs were added at the start of the forcing experiment for that collection date (e.g., for the collection on 10 Jan 2020, 88.1 GDHs, or 3.7 GDDs, had accumulated since 1 Jan 2020). However, this value was negligible as field heat accumulation was no higher than 178.2 GDHs / 7.4 GDDs for any collection after 1 Jan for both years. Chill was calculated using the Dynamic Model in the chillR package starting 1 October [49–51,141]. The two seasons were analyzed independently since there was not complete overlap in the accessions and chill portions at each collection date between years. For 2019-2020, 5 early- blooming genotypes (with 0k) and 3 late-blooming genotypes (with 2k) were scored and analyzed. For 2021-2020, 7 early-blooming genotypes (5 overlapping between years) and 4 late- blooming genotypes (3 overlapping between years) were scored and analyzed. For each year, separate statistical models, including both fixed and random effects, were constructed for BBCH stages 51, 52, 53, and 54. Models were not created for later stages due to concerns about branch health and thus, reliability of the data. Proportion of flowers at or above a particular stage (X) were modeled as such: Proportion BBCH stage X+ ~ Group+chill portions+GDDs+chill portions*GDDs+(1 | dummy/rep) +error where ‘Group’ was either early- or late-blooming (0k or 2k, respectively), and ‘chill portions’ was the amount of chill accumulated between Oct 1st and the collection date according to the Dynamic Model. ‘Chill portions’ was coded as an ordered factor as it was equivalent with collection date. GDDs was a continuous variable. To account for the repeated measurements within each collection, a dummy variable was first coded as a combination of an individual tree within each timepoint. Then, using the ‘corSymm’ function of R package ‘nlme’, the correlation structure was modeled by nesting the repeated measurement within this dummy variable [142]. Raw data for the two seasons of forcing experiments are given as Supplementary File 4.4 – 4.5, and R scripts for data analysis and figure generation can be found at https://github.com/goeckeritz. 58 Measuring H2O2 content in early- and late-blooming individuals during dormancy transitions We tried several fluorometric / spectrometric methods to measure hydrogen peroxide; in the event this saves researchers time and effort in the future, they are briefly described here. Two of these three methods did not pass rigorous controls. As the first control, samples were measured at 5 - 30 min intervals for at least a half hour to ensure stable signal and therefore, stable concentration of H2O2. For the second control, known concentrations of H2O2 were added to the sample extract and the correlation between the expected and observed quantities of H2O2 were noted. The first method used Amplex Red™ as a fluorophore and is described in Chakraborty et al. 2017 [143]. We tried this method with and without a Dowex Ion Exchange Resin 1X8-200 (cat no. 217425, Sigma-Aldrich) to reduce contaminants; however, in either case, repeated fluorometric measurements of the samples at an excitation set to a wavelength of 550 nm and emission at 590 nm showed a linear increase in signal over time, suggesting the extracted sample H2O2 was not stable. This auto generation of signal was also reflected in the second control method, as spike-ins of known H2O2 concentrations with sample extract did not yield expected values (the relationship between expected and observed was not linear). Next, we tried a spectrometric method using 3-methyl-2-benzothiazolinone hydrazone and 3-(dimethylamino) benzoic acid (MBTH-DMAB) [144,145]. The signal stabilized over time, but the assay failed the second control: adding known volumes of H2O2 did not yield expected concentrations at a 590 nm absorption. Therefore, for both the Amplex Red™ and MTBH- DMAB assays, we posited that unknown compounds interfered with true signal caused by H2O2. Lastly, we adapted the xylenol-orange method of Cheeseman et al. 2006 to measure samples in a 96-well plate [118]. This method passed both controls but was sensitive to temperature and humidity. This may have been due to the solubility of KCN, a component of the extraction buffer. Therefore, the plate (Greiner Bio-One, Ref 655903) was pre-cooled prior to sample loading and kept on ice at all other times, being careful not to smudge the underside as this would affect signal readings. Additionally, all reagents and sample extracts were constantly kept on ice and all centrifugation steps were performed at 4º C. The humidity was monitored with a dehumidifier and kept under 30% during extraction to prevent condensation buildup on liquid N2 chilled mortars and pestles. Periodically during the dormant season (October – March; Table S4.5), at least three biological replications of four adjacent flower buds were collected per tree and immediately frozen in liquid nitrogen. Flower buds of four late-blooming and six early- 59 blooming individuals were collected at all five time points. Given reports that ROS degrade over time even when materials are stored at -80º C, all samples were kept in liquid nitrogen in a large Dewar until extraction [118]. On the day of extraction, the fresh weight of each sample (four flower buds) was determined by first measuring each tube containing 5 mL extraction buffer without sample, weighing the tube again after adding ground sample material (thawed), and calculating the difference. Samples were ground to a fine powder using a mortar and pestle before quickly adding each to 5 mL of cold NaH2PO4 buffer at pH 6.4 + 5 mM KCN, which was then vortexed immediately. Eight fresh H2O2 standards between concentrations of 0 – 4 µM were prepared alongside samples by serial diluting a 30% H2O2 solution (Thermo Scientific, H325-500) in extraction buffer. Vortexed samples were then centrifuged at 4º C at max speed for 10 minutes. Meanwhile, a working solution was prepared containing 500 µM ferrous (II) ammonium sulphate, 200 µM sorbitol, 50 mM H2SO4, 2% ethanol and 200 µM of xylenol orange. A ‘pre- working solution,’ without xylenol orange, was stored at room temperature; xylenol orange was not added until the day of extraction. After centrifuging samples, the supernatants were carefully transferred to new tubes already on ice. The chilled 96-well plate was then loaded by first adding 100 µL of the working solution to the empty wells. Then, 100 µL of the standards and samples were added to the wells with working solution. Each standard and sample had four technical replicates. After loading, the plate was placed back in a 4º C fridge and allowed to incubate for 60 – 90 minutes (the reaction is slower in the cold). Multiple readings of the plate were sometimes warranted if the standard curve had not yet linearized, indicating the reaction was incomplete. The plate was then removed from the fridge and allowed to sit at room temperature for 30 seconds before measuring the absorption at 550 and 800 nM using a Synergy H1 microplate reader (BioTek Instruments, Inc). The data was exported to an Excel spreadsheet with Gen5 software version 2.01.14 (BioTek Instruments, Inc) and µmoles of H2O2 per gram of fresh weight for every well were determined using the standard curve. In 4 of 14 cases, the 0 µM standard was dropped since it significantly skewed the curve. Standard curve r2 values ranged from 0.86 – 0.99. After all extractions, there were 2-3 biological replications for each tree per time point. The detailed protocol can be found at https://github.com/goeckeritz. 60 First, the raw data was filtered by dropping negative values (representing readings where background noise was higher than the signal from the sample). Then the data was modeled as: µmoles of H2O2 ~ Group + chill + Group*chill + (1 | Genotype:Group) + (1 | dummy) + error Where the fixed factor ‘Group’ has levels early- and late-blooming, and the fixed factor ‘chill’ is how much field chill has accumulated since October 1st, 2019 (synonymous with collection date). ‘(1 | Genotype:Group)’ represents the random variation introduced by individual genotypes within a bloom group, and ‘(1 | dummy)’ represents the technical replicate variation, which is a combination of repeated measurements nested within a plate, genotype, and collection date (chill level). Manual curation of gene models within the chr4 QTL region Prior to predicting genetic variation effects on gene models within the chr4 QTL region on subgenome A (pipeline explained below), the genes within this region were manually examined against the RNAseq data from the current study using Apollo v 2.6.5 [146]. The RNAseq reads were aligned to ‘Montmorency’ subgenome A as described below. Manual curation for these genes was done to refine their quality for more accurate candidate gene selection. Genome sequencing, coverage assessments, SNP-calling, and structural variant analysis DNA of young, unexpanded leaves was extracted from 14 individuals (7 early- blooming, 5 late-blooming, and the two parents) with a DNeasy Mini Plant Kit (Qiagen, Hilden, Germany) and quantified with a Quant-iTTM PicoGreenTM dsDNA assay kit (Thermo Fisher Scientific, Waltham, MA) or using the Qubit DNA Broad Range kit and Qubit 4 Fluorometer (Thermo Fisher Scientific Inc, Waltham, MA). Quality was assessed with an Agilent 4200 TapeStation HS DNA1000 and Kapa Biosystems Illumina Library Quantification qPCR assays. Illumina TruSeq Nano DNA 150 bp paired-end (PE) libraries were prepared at MSU’s Research Technology Support Facility (RTSF) and sequenced at 20 – 40X coverage per sample on one lane of a HiSeq4000 instrument. Between 22.6 and 35.8 million paired end reads were produced for each sample and checked for quality using FASTQC [147]. After quality checking, reads were trimmed using Trimmomatic v 0.39 and aligned to the ‘Montmorency’ reference genome (including subgenomes A, A’, and B) using BWA v 0.7.17 [148]. The resulting .bam files were sorted, and duplicates were removed with picard tools v 2.25.0 [149]. 61 Previously, seven haplotypes were discovered at a major bloom time locus on linkage group 4 [87]. ‘Montmorency’ contains phased subgenomes A, A’ and B, and recent findings suggest diverse sour cherries contain variable dosages of these subgenomes (Iezzoni and Rhoades, personal communication). This merited exploration of subgenome dosage at the locus associated with bloom time in individuals of the segregating population. Knowing the haplotypes of the 12 sequenced genotypes as well as the parents ‘Balaton’ and ‘Surefire,’ we used coverage plots to deduce which subgenome the k haplotype most likely originated from. To produce the plots, first read depth was calculated for every position using SAMtools v 1.11 [150]. Then, a custom bash loop was written to average the read depth over every 500 kb window for all 14 samples (see coverage_loops.sh, https://github.com/goeckeritz). These files were plotted in R using ggplot2 to visualize and determine the subgenome dosage at the bloom time locus [141,151]. From there, logical statements comparing haplotypes and subgenome composition amongst individuals were developed to identify which haplotype originated from which subgenome (described fully in Supplementary File 4.7). Since k was determined to be located on subgenome A, we excluded Early 1 (A’BBB) from genetic variation comparisons. SNP calls were made only in the QTL region (Pcer_chr4A: 8324991- 16386209), and ploidy was adjusted for every individual according to their subgenome A dosage. This required two separate regions to be called for ‘Surefire’ and Early 7, as their dosage increased (due to the g haplotype) at approximately 12583308 bp along the length of chromosome 4 (determined by examining depth files). We used HaplotypeCaller within GATK v 4.2.0.0 and a pedigree file to call SNPs for all samples [152]. Next, CombineGVCFs was used to merge the two parts of ‘Surefire’ and Early 7; then, these complete files were combined with all others individuals’ gvcf files. The variants were then hard-filtered based on quality metrics (see variant_filtration.sh, https://github.com/goeckeritz) and assigned a predicted effect on the gene model with snpEff [153]. VariantsToTable was used to convert this file to a table, which was read into R for further filtering based on known inheritance patterns. For example, variants of interest include those that are unique to the individuals with 2k (late-blooming); and, if k is an individual’s only source of subgenome A (i.e., it is 2X subgenome A and does not have haplotype i or g in addition to 2k), these late-bloomers’ genotype calls should be identical. In this situation, each late-blooming individual will have two alleles (e.g., T/T, A/A, C/T, etc.), one of which should come from ‘Balaton’ and the other from ‘Surefire.’ Using such logic, we further 62 filtered the dataset (see further_filtering_and_plotting.R). Regular expressions were then used to extract lists of genes with different classifications of variants (e.g., high effect mutations, 5’UTR, upstream, etc.) from the filtered file. To identify structural variants upstream and within protein-coding genes, reads mapped to the chr4A QTL region were pooled so that higher coverage would represent both the k and i haplotypes. These pools were genotyped alongside ‘Surefire’, which was sequenced to approximately 60X (15X per allele), and Early 7 (20X, 5X per allele). Surefire represented haplotypes i and k in the first half of the QTL and haplotypes i, k, and g in the latter half. Early 7 was genotyped separately for another representative of the g haplotype. Variants were called with Delly v 0.7.8 [154]. Limited coverage representing the g haplotypes meant fewer variants were called compared to i and k, and the reads representing g in Surefire were confounded with i and k. Nonetheless, several criteria were imposed to filter variants. First, k’s unique genotype calls were anticipated to be homozygous and different than i calls, which were also expected to be homozygous. ‘Surefire’ was useful when i calls were missing in the first half of the QTL, because homozygous calls in Surefire meant i and k were identical and such a variant could be ruled out as causative for late bloom. The resulting 169 structural variants is a liberal estimate since in many instances, calls for haplotypes other than k were missing. To compensate for this, called variants predicted to affect DEGs of interest were manually verified in IGV v 2.13.0 by examining coverage of individuals with i, k, and g [155]. RNA sequencing, differential gene expression, and weighted gene correlation network analyses RNA from summer apices was extracted using a previously-developed protocol treated with a DNA-freeTM Kit (Thermo Fisher Scientific, Waltham, MA) to remove DNA, and assessed for initial quality using a Nanodrop OneC (Thermo Fisher Scientific, Waltham, MA) [156]. If purity was insufficient (260/230 < 1.8), samples were further cleaned and concentrated using the RNA Clean & ConcentratorTM-25 kit (Zymo Research, Seattle, WA). Concentration of each sample was measured using the Qubit RNA Broad Range kit and Qubit 4 Fluorometer (Thermo Fisher Scientific Inc, Waltham, MA), and quality was assessed with an Agilent High Sensitivity RNA ScreenTape System (Agilent Technologies, Santa Clara, CA). Illumina stranded mRNA 150 bp PE libraries were prepared at MSU’s (RTSF) before they were pooled and sequenced on approximately 0.8 lanes of a NovaSeq 6000, producing at least 39.8 million reads per library. A total of 48 libraries were sequenced: three replicates per tree for 12 June, two replicates per tree 63 for 29 June, and three replicates per tree for 18 July. The same three trees were used as representatives of each bloom group (EB and LB) for each collection date (i.e., six trees were used). Each resulting fastq file was confirmed to be of high quality with FASTQC. Reads were first processed with Trimmomatic v 0.39, then aligned with STAR v 2.7.3a to P. cerasus ‘Montmorency’ subgenome A since the coverage assessment results suggested the gene(s) responsible for the k late bloom phenotype were unique alleles of subgenome A [157,158]. We were reassured that reads from subgenome B and A’ would align to orthologous genes on subgenome A since Prunus genomes are highly syntenic and collinear, and because alignment parameters were relaxed to allow for expected mismatches while still retaining sufficient mapping quality. Briefly, parameters –outFilterScoreMinOverLread (allowable ratio of mismatches to mapped length) --outFilterMatchNminOverLread, (minimum requirement of matched bases, normalized by read length) and --outFilterMatchNmin (minimum requirement of matched bases) were all set to zero to ensure maximum mapping. No more than 30% of each read’s mapped length was allowed to mismatch (--outFilterMismatchNoverLmax). These parameters, however, had the potential of keeping very short (and potentially inaccurate) mappings; therefore, we checked the average input read length versus average mapped read length. The average input read length was 294 and above for each library (2 x 150 paired-end reads), and the lowest average mapped length was 251; thus, on average, at least 85% of each read’s full length was successfully mapped. Uniquely-mapped reads ranged from 87 – 92% per library, and there was no obvious mapping bias among certain trees or time points. These alignments were used with Montmorency subgenome A gene models (.gff3) to assess the expression of each gene. These values were determined using StringTie v 2.1.3 and expressed as counts for the required input for the DESeq2 R package, and Fragments Per Kilobase per transcript per Million mapped reads (FPKM) as the required input for the WGCNA (Weighted Gene Correlation Network Analysis) R package [141,159–161] For the differential expression (DE) analysis using DESeq2, transcripts with no expression across all samples were dropped from the analysis. Using a regularized log transformation, normalized expression counts for each library were visualized by a Principal Component Analysis (PCA). Two of the 48 samples were identified as outliers (27_04_34_6_12_19_R2 and 27_03_46_6_29_19_R3) and removed from further analysis. Counts were normalized and modeled using the DESeqDataSetFromMatrix function with the 64 following independent variables: bloom group (early or late / 0k or 2k), the individual nested within the bloom group, date, and the interaction of bloom and date. Seven comparisons were performed to assess gene expression changes: 1) early_v_late_Jun12, 2) early_v_late_Jun29, 3) early_v_late_Jul18, 4) early_Jun12_Jun29, 5) early_Jun29_Jul18, 6) late_Jun12_Jun29, and 7) late_Jun29_Jul18. These comparisons correspond to gene expression changes between early- and late-bloomers on 12 June, 29 June, and 18 July, then early-bloomers from 12 June to 29 June and 29 June to 18 July, and late-bloomers from 12 June to 29 June and 29 June to 18 July, respectively. Differentially expressed genes with adjusted p-values < 0.05 and a log2FoldChange > 0.5 were kept for further analysis. Genome-wide results and the subset of results within the QTL region on chromosome 4 can be found in Supplementary File 4.9. Venn diagrams were created using R package VennDiagram [162]. For network analyses, genes with an FPKM value < 5 across all samples were excluded. Log-transformed FPKM values were further filtered to reduce contributions of genes that varied little between samples in the analysis (coefficient of variation < 0.2). Using hierarchal clustering, the two samples identified as outliers during the differential expression analysis using DESeq2 were again removed from the dataset prior to weighted gene correlation network analysis. A signed adjacency matrix with a power threshold of 14 was created and used to calculate the topological overlap dissimilarity matrix. This matrix was used to create a gene tree, which was used as input for cutreeDynamic. DeepSplit was set to 2, the minClusterSize (minimum number of genes in each module) was set to 30, and the cut height was 0.996. Next, the eigengenes were clustered to identify modules with similar expression patterns and modules below cutHeight < 0.20 were merged. The eigengene expression for each of the final modules was then plotted by date to identify those with patterns associating with bloom group (Figure 4.17). GO database creation and enrichment tests GO (gene ontology) IDs for each Montmorency subgenome A gene were determined with InterProScan v 5.33-72.0 as previously described and BLAST+ v 2.9 [163,164]. We formatted a BLAST+ database for the arabidopsis proteins, and these were used as query while sour cherry protein sequences were used as the target [165]. Only hits with p-value < 1.0 x 10-10 were kept and all other parameters were set to default. Arabidopsis GO information was downloaded from arabidopsis.org and a simple R script was used to merge the InterProScan, 65 arabidopsis to sour cherry blast hits, and arabidopsis GO terms to create the final database (see https://github.com/goeckeritz). The database is available as Supplementary File 4.14. topGO was used to determine which biological processes (BP) were enriched in the differentially expressed genes and WGCNA modules [166]. topGOfunctions.R is a function wrapper created previously by Mansfeld et al. [167]. Since the dds object from the differential expression analysis with DESeq2 was used as input for all GO enrichment tests, R code performing both the DE and WGCNA module tests can be found at the end of the DESeq2 R script. Using the ‘weight01’ algorithm, terms were considered significant if the Fisher test had a p-value < 0.05 [168]. Minimum node size was set to 50. The results of these tests can be found in Supplementary File 4.11 (DE analysis) and Supplementary File 4.12 (WGCNA modules). Comparing Regina QTL syntelogs A syntenic analysis between ‘Regina’ chromosome 4 and ‘Montmorency’ chromosome 4A was conducted using MCScan, and the sour cherry syntelogs of genes within the sweet cherry QTL region were used as anchors to identify this region in sour cherry [45]. All genes within this region were searched for differential expression between early- and late-blooming individuals for any sampling date in the current study. 66 REFERENCES 1. Xiang Y, Huang CH, Hu Y et al. Evolution of Rosaceae fruit types based on nuclear phylogeny in the context of geological times and genome duplication. Mol Biol Evol 2017;34:262–81. 2. Kistner E, Kellner O, Andresen J et al. Vulnerability of specialty crops to short-term climatic variability and adaptation strategies in the Midwestern USA. Climatic Change 2018;146:145–58. 3. Unterberger C, Brunner L, Nabernegg S et al. Spring frost risk for regional apple production under a warmer climate. PLoS ONE 2018;13:1–18. 4. Marino GP, Kaiser DP, Gu L et al. Reconstruction of false spring occurrences over the southeastern United States, 1901–2007: an increasing risk of spring freeze damage? Environ Res Lett 2011;6:024015. 5. Bird KA, Jacobs M, Sebolt A et al. Parental origins of the cultivated tetraploid sour cherry (Prunus cerasus L.). Plants People Planet 2022;4:444–50. 6. Iezzoni AF, Hancock AM. A comparison of pollen size in sweet and sour cherry. HortScience 1984;19:560–2. 7. Brettin TS, Karle R, Crowe EL et al. Chloroplast inheritance and DNA variation in sweet, sour, and ground cherry. J Hered 2000;91:75–9. 8. Olden EJ, Nybom N. On the origin of Prunus cerasus L. Hereditas 1968;59:327–59. 9. Goeckeritz CZ, Rhoades KE, Childs KL et al. Genome of tetraploid sour cherry (Prunus cerasus L.) ‘Montmorency’ identifies three distinct ancestral Prunus genomes. Hortic Res 2023;10:uhad097. 10. Quero-García J, Iezzoni AF, Pulawska J et al. eds. Cherries: Botany, Production, and Uses. Boston, MA: CABI International, 2017. 11. Hung H-Y, Shannon LM, Tian F et al. ZmCCT and the genetic basis of day-length adaptation underlying the postdomestication spread of maize. Proc Natl Acad Sci 2012;109:E1913–21. 12. Rohde A, Bastien C, Boerjan W. Temperature signals contribute to the timing of photoperiodic growth cessation and bud set in poplar. Tree Physiol 2011;31:472–82. 13. Hardigan MA, Laimbeer FPE, Newton L et al. Genome diversity of tuber-bearing Solanum uncovers complex evolutionary history and targets of domestication in the cultivated potato. Proc Natl Acad Sci 2017:201714380. 14. Salomé PA, Bomblies K, Laitinen RAE et al. Genetic architecture of flowering-time variation in Arabidopsis thaliana. Genet 2011;188:421–33. 15. Chardon F, Virlon B, Moreau L et al. Genetic architecture of flowering time in maize as inferred from quantitative trait loci meta-analysis and synteny conservation with the rice genome. Genet 2004;168:2169–85. 67 16. Diaz D, Rasmussen H, Dennis F. Scanning electron microscope examination of flower bud differentiation in sour cherry. J Am Soc Hortic 1981;106:513–5. 17. Guimond C, Andrews P, Lang G. Scanning electron microscopy of floral initiation in sweet cherry. J Amer Soc Hort Sci 1998;123:509–12. 18. Li Y, Zhang D, Zhang X et al. A transcriptome analysis of two apple (Malus × domestica) cultivars with different flowering abilities reveals a gene network module associated with floral transitions. Sci Hortic 2018;239:269–81. 19. Li Y, Zhang D, An N et al. Transcriptomic analysis reveals the regulatory module of apple (Malus × domestica) floral transition in response to 6-BA. BMC Plant Biol 2019;19, DOI: 10.1186/s12870-019-1695-0. 20. Villar L, Lienqueo I, Llanes A et al. Comparative transcriptomic analysis reveals novel roles of transcription factors and hormones during the flowering induction and floral bud differentiation in sweet cherry trees (Prunus avium L. cv. Bing). PLOS ONE 2020;15:e0230110. 21. Vimont N, Fouché M, Campoy JA et al. From bud formation to flowering: transcriptomic state defines the cherry developmental phases of sweet cherry bud dormancy. BMC genomics 2019;20:974. 22. Bradley MV, Crane JC. Gibberellin-induced inhibition of bud development in some species of Prunus. Science Magazine 1960;131:825–6. 23. Gottschalk C, Zhang S, Schwallier P et al. Genetic mechanisms associated with floral initiation and the repressive effect of fruit on flowering in apple (Malus × domestica Borkh). PLoS ONE 2021;16, DOI: 10.1371/journal.pone.0245487. 24. Zhang S, Gottschalk C, Van Nocker S. Genetic mechanisms in the repression of flowering by gibberellins in apple (Malus × domestica Borkh.). BMC Genom 2019;20, DOI: 10.1186/s12864- 019-6090-6. 25. Hu J, Mitchum MG, Barnaby N et al. Potential sites of bioactive gibberellin production during reproductive growth in Arabidopsis. Plant Cell 2008;20:320–36. 26. Rieu I, Eriksson S, Powers SJ et al. Genetic analysis reveals that C19-GA 2-Oxidation is a major gibberellin inactivation pathway in Arabidopsis. Plant Cell 2008;20:2420–36. 27. Yamaguchi N, Winter CM, Wu M et al. Gibberellin acts positively then negatively to control onset of flower formation in Arabidopsis. Science 2014;344:638–41. 28. Kofler J, Milyaev A, Capezzone F et al. High crop load and low temperature delay the onset of bud initiation in apple. Sci Rep 2019;9, DOI: 10.1038/s41598-019-54381-x. 29. Raseira MCB, Moore JN. Time of flower bud initiation in peach cultivars differing in chilling requirement. HortSci 1987;22:216–8. 30. McArtney S, Hoover E, Booking IR. Seasonal variation in the onset and duration of flower development in ‘Royal Gala’ apple buds. J Hortic Sci Biotechnol 2001;75:536–40. 68 31. Heide OM, Rivero R, Sønsteby A. Temperature control of shoot growth and floral initiation in apple (Malus × domestica Borkh.). CABI 2020;1, DOI: 10.1186/s43170-020-00007-6. 32. Sønsteby A, Heide OM. Temperature effects on growth and floral initiation in sweet cherry (Prunus avium L.). Sci Hortic 2019;257:108762. 33. Goeckeritz C, Hollender CA. There is more to flowering than those DAM genes: the biology behind bloom in rosaceous fruit trees. Current Opinion in Plant Biology 2021;59:101995. 34. Heide OM. Interaction of photoperiod and temperature in the control of growth and dormancy of Prunus species. Sci Hortic 2008;115:309–14. 35. Heide OM, Prestrud AK. Low temperature, but not photoperiod, controls growth cessation and dormancy induction and release in apple and pear. Tree Physiol 2005;25:109–14. 36. Tylewicz S, Petterle A, Marttila S et al. Photoperiodic control of seasonal growth is mediated by ABA acting on cell-cell communication. Science 2018;360:212–5. 37. Singh RK, Maurya JP, Azeez A et al. A genetic network mediating the control of bud break in hybrid aspen. Nat Commun 2018;9, DOI: 10.1038/s41467-018-06696-y. 38. Singh RK, Miskolczi P, Maurya JP et al. A tree ortholog of SHORT VEGETATIVE PHASE floral repressor mediates photoperiodic control of bud dormancy. Curr Biol 2019;29:128-133.e2. 39. André D, Marcon A, Lee KC et al. FLOWERING LOCUS T paralogs control the annual growth cycle in Populus trees. Curr Biol 2022;32:2988–96. 40. Lang G, Early J, Martin G et al. Endo-, para-, and ecodormancy: physiological terminology and classification for dormancy research. HortScience 1987;22:371–7. 41. Fan S, Bielenberg DG, Zhebentyayeva TN et al. Mapping quantitative trait loci associated with chilling requirement, heat requirement and bloom date in peach (Prunus persica). New Phytol 2010;185:917–30. 42. Sánchez-Pérez R, Dicenta F, Martínez-Gómez P. Inheritance of chilling and heat requirements for flowering in almond and QTL analysis. Tree Genet Genomes 2012;8:379–89. 43. Allard A, Bink MCAM, Martinez S et al. Detecting QTLs and putative candidate genes involved in budbreak and flowering time in an apple multiparental population. Journal of Experimental Botany 2016;67:2875–88. 44. Kitamura Y, Habu T, Yamane H et al. Identification of QTLs controlling chilling and heat requirements for dormancy release and bud break in Japanese apricot (Prunus mume). Tree Genet Genomes 2018;14, DOI: doi.org/10.1007/s11295-018-1243-3. 45. Branchereau C, Quero-García J, Zaracho-Echagüe NH et al. New insights into flowering date in Prunus: fine mapping of a major QTL in sweet cherry. Hortic Res 2022;9:uhac042. 46. Okie WR, Blackburn B. Increasing chilling reduces heat requirement for floral budbreak in peach. HortSci 2011;46:245–52. 69 47. Weinberger J. Chilling Requirements of Peach Varieties. American Society of Horticultural Science 1950;56:122–8. 48. Richardson E, Seeley S, Walker R. A model estimating the completion of rest for “Red Haven” and “Elberta” peach. HortSci 1974;9:331–2. 49. Luedeling E, Kunz A, Blanke MM. Identification of chilling and heat requirements of cherry trees - a statistical approach. Int J Biometeorol 2013;57:679–89. 50. Fishman S, Erez A, Couvillon G. The temperature dependence of dormancy breaking in plants: Computer simulation of processes studied under controlled temperatures. J Theor Biol 1987;126:309–21. 51. Fishman S, Erez A, Couvillon G. The temperature dependence of dormancy breaking in plants: Mathematical analysis of a two-step model involving a cooperative transition. J Theor Biol 1987;124:473–83. 52. Campoy JA, Darbyshire R, Dirlewanger E et al. Yield potential definition of the chilling requirement reveals likely underestimation of the risk of climate change on winter chill accumulation. Int J Biometeorol 2019;63:183–92. 53. Luedeling E. Climate change impacts on winter chill for temperate fruit and nut production: A review. Sci Hortic 2012;144:218–29. 54. Hillmann L, Elsysy M, Goeckeritz C et al. Preanthesis changes in freeze resistance, relative water content, and ovary growth preempt bud phenology and signify dormancy release of sour cherry floral buds. Planta 2021;254:1–12. 55. Fadón E, Herrero M, Rodrigo J et al. Dormant flower buds actively accumulate starch over winter in sweet cherry. Front Plant Sci 2018;9, DOI: 10.3389/fpls.2018.00171. 56. Kaufmann H, Blanke M. Changes in carbohydrate levels and relative water content (RWC) to distinguish dormancy phases in sweet cherry. Journal of Plant Physiology 2017, DOI: 10.1016/j.jplph.2017.07.004. 57. Sapkota S, Liu J, Islam MT et al. Changes in reactive oxygen species, antioxidants and carbohydrate metabolism in relation to dormancy transition and bud break in apple (Malus × domestica Borkh.) cultivars. Antioxidants 2021;10, DOI: 10.3390/antiox10101549. 58. Yang Q, Yang B, Li J et al. ABA-responsive ABRE-BINDING FACTOR3 activates DAM3 expression to promote bud dormancy in Asian pear. Plant Cell Environ 2020;43:1360–75. 59. Zhuang W, Gao Z, Wang L et al. Comparative proteomic and transcriptomic approaches to address the active role of GA4 in Japanese apricot flower bud dormancy release. J Exp Bot 2013;64:4953–66. 60. Zhuang W, Gao Z, Wen L et al. Metabolic changes upon flower bud break in Japanese apricot are enhanced by exogenous GA4. Hortic Res 2015;2, DOI: 10.1038/hortres.2015.46. 70 61. Falavigna V da S, Guitton B, Costes E et al. I want to (bud) break free: The potential role of DAM and SVP-like genes in regulating dormancy cycle in temperate fruit trees. Front Plant Sci 2019;10:1–17. 62. Beauvieux R, Wenden B, Dirlewanger E. Bud dormancy in perennial fruit tree species: a pivotal role for oxidative cues. Front Plant Sci 2018;9:1–13. 63. Fadón E, Fernandez E, Behn H et al. A conceptual framework for winter dormancy in deciduous trees. Agron 2020;10:241. 64. Cooke JEK, Eriksson ME, Junttila O. The dynamic nature of bud dormancy in trees: environmental control and molecular mechanisms. Plant Cell Environ 2012;35:1707–28. 65. Ionescu IA, López-Ortega G, Burow M et al. Transcriptome and metabolite changes during hydrogen cyanamide-induced floral bud break in sweet cherry. Front Plant Sci 2017;8, DOI: 10.3389/fpls.2017.01233. 66. Sudawan B, Chang C-S, Chao H et al. Hydrogen cyanamide breaks grapevine bud dormancy in the summer through transient activation of gene expression and accumulation of reactive oxygen and nitrogen species. BMC Plant Biol 2016;16:202. 67. Roussos PA. Adventitious root formation in plants: The implication of hydrogen peroxide and nitric oxide. Antioxidants 2023;12:862. 68. Velappan Y, Signorelli S, Considine MJ. Cell cycle arrest in plants: What distinguishes quiescence, dormancy and differentiated G1? Ann Bot 2017;120:495–509. 69. Wang D, Gao Z, Du P et al. Expression of ABA metabolism-related genes suggests similarities and differences between seed dormancy and bud dormancy of peach (Prunus persica). Front Plant Sci 2016;6:1–17. 70. Lv L, Huo X, Wen L et al. Isolation and role of PmRGL2 in GA-mediated floral bud dormancy release in Japanese apricot (Prunus mume Siebold et Zucc.). Front Plant Sci 2018;9, DOI: 10.3389/fpls.2018.00027. 71. Wen Z, Cao X, Hou Q et al. Expression profiling and function analysis highlight the positive involvement of sweet cherry PavTCP17 in regulating flower bud dormancy. Sci Hortic 2023;318, DOI: 10.1016/j.scienta.2023.112138. 72. Rodriguez A, Sherman W, Scorza R et al. ‘Evergreen’ peach, its inheritance and dormant behavior. J Am Soc Hortic 1994;119:789–92. 73. Bielenberg DG, Wang Y, Fan S et al. A deletion affecting several gene candidates is present in the evergrowing peach mutant. J Hered 2004;95:436–44. 74. Prudencio AS, Díaz-Vivancos P, Dicenta F et al. Monitoring the transition from endodormancy to ecodormancy in almond through the analysis and expression of a specific class III peroxidase gene. Tree Genet Genomes 2019;15, DOI: 10.1007/s11295-019-1351-8. 71 75. de la Fuente L, Conesa A, Lloret A et al. Genome-wide changes in histone H3 lysine 27 trimethylation associated with bud dormancy release in peach. Tree Genet Genomes 2015;11, DOI: 10.1007/s11295-015-0869-7. 76. Leida C, Conejero A, Arbona V et al. Chilling-dependent release of seed and bud dormancy in peach associates to common changes in gene expression. PLOS ONE 2012;7, DOI: 10.1371/journal.pone.0035777. 77. Bai S, Saito T, Sakamoto D et al. Transcriptome analysis of Japanese pear (Pyrus pyrifolia ’Nakai’) flower buds transitioning through endodormancy. Plant Cell Physiol 2013;54:1132–51. 78. Rothkegel K, Sánchez E, Montes C et al. DNA methylation and small interference RNAs participate in the regulation of MADS-box genes involved in dormancy in sweet cherry (Prunus avium L.). Tree Physiol 2017;37:1739–51. 79. Zhu H, Chen PY, Zhong S et al. Thermal-responsive genetic and epigenetic regulation of DAM cluster controlling dormancy and chilling requirement in peach floral buds. Hortic Res 2020;7, DOI: 10.1038/s41438-020-0336-y. 80. Zhao Y-L, Li Y, Cao K et al. MADS-box protein PpDAM6 regulates chilling requirement- mediated dormancy and bud break in peach. Plant Physiol 2023:1–18. 81. Horvath DP, Sung S, Kim D et al. Characterization, expression and function of DORMANCY ASSOCIATED MADS-BOX genes from leafy spurge. Plant Molecular Biology 2010;73:169–79. 82. Wu R, Tomes S, Karunairetnam S et al. SVP-like MADS box genes control dormancy and budbreak in apple. Front Plant Sci 2017;8, DOI: 10.3389/fpls.2017.00477. 83. Jimenez S, Reighard GL, Bielenberg DG. Gene expression of DAM5 and DAM6 is suppressed by chilling temperatures and inversely correlated with bud break rate. Plant Mol Biol 2010;73:157–67. 84. Calle A, Cai L, Iezzoni A et al. Genetic dissection of bloom time in low chilling sweet cherry (Prunus avium L.) using a multi-family QTL approach. Front Plant Sci 2020;10, DOI: 10.3389/fpls.2019.01647. 85. Cattani AM, Falavigna V da S, Silveira CP et al. Type-B cytokinin response regulators link hormonal stimuli and molecular responses during the transition from endo- to ecodormancy in apple buds. Plant Cell Rep 2020;39:1687–703. 86. Ballester J, Socias I Company R, Arus P et al. Genetic mapping of a major gene delaying blooming time in almond. Plant Breeding 2001;120:268–70. 87. Cai L, Stegmeir T, Sebolt A et al. Identification of bloom date QTLs and haplotype analysis in tetraploid sour cherry (Prunus cerasus). Tree Genet Genomes 2018;14, DOI: 10.1007/s11295- 018-1236-2. 88. Dirlewanger E, Quero-García J, Dantec LL et al. Comparison of the genetic determinism of two key phenological traits, flowering and maturity dates, in three Prunus species: peach, apricot and sweet cherry. J Hered 2012;109:280–92.; 72 89. Quilot B, Wu BH, Kervella J et al. QTL analysis of quality traits in an advanced backcross between Prunus persica cultivars and the wild relative species P. davidiana. Theor Appl Genet 2004;109:884–97. 90. Silva C, Garcia-Mas J, Sánchez AM et al. Looking into flowering time in almond (Prunus dulcis (Mill) D. A. Webb): The candidate gene approach. Theor Appl Genet 2005;110:959–68. 91. Castede S, Campoy JA, Dantec LL et al. Mapping of candidate genes involved in bud dormancy and flowering time in sweet cherry (Prunus avium). PLOS ONE 2015;10:1–18. 92. Moon J, Suh S, Lee H et al. The SOC1 MADS-box gene integrates vernalization and gibberellin signals for flowering in Arabidopsis. TPJ 2003;35:613–23. 93. Noh Y, Amasino RM. PIE1, an ISWI family gene, is required for FLC activation and floral repression in Arabidopsis. Plant Cell 2003;15:1671–82. 94. Wu G, Poethig RS. Temporal regulation of shoot development in Arabidopsis thaliana by miR156 and its target SPL3. Development 2006;133:3539–47. 95. Verde I, Jenkins J, Dondini L et al. The Peach v2.0 release: high-resolution linkage mapping and deep resequencing improve chromosome-scale assembly and contiguity. BMC Genom 2017;18, DOI: 10.1186/s12864-017-3606-9. 96. Jiang F, Zhang J, Wang S et al. The apricot (Prunus armeniaca L.) genome elucidates Rosaceae evolution and beta-carotenoid synthesis. Hortic Res 2019;6, DOI: 10.1038/s41438- 019-0215-6. 97. Yi XG, Yu XQ, Chen J et al. The genome of Chinese flowering cherry (Cerasus serrulata) provides new insights into Cerasus species. Hortic Res 2020;7, DOI: 10.1038/s41438-020- 00382-1. 98. Cao K, Peng Z, Zhao X et al. Chromosome-level genome assemblies of four wild peach species provide insights into genome evolution and genetic basis of stress resistance. BMC Biol 2022;20, DOI: 10.1186/s12915-022-01342-y. 99. Wang J, Liu W, Zhu D et al. Chromosome-scale genome assembly of sweet cherry (Prunus avium L.) cv. Tieton obtained using long-read and Hi-C sequencing. Hortic Res 2020;7:1–11. 100. Baxter I. We aren’t good at picking candidate genes, and it’s slowing us down. Curr Opin Plant Biol 2020;54:57–60. 101. Fadón E, Herrero M, Rodrigo J. Flower development in sweet cherry framed in the BBCH scale. Sci Hortic 2015;192:141–7. 102. Fadón E, Rodrigo J, Herrero M. Is there a specific stage to rest? Morphological changes in flower primordia in relation to endodormancy in sweet cherry (Prunus avium L.). Trees 2018;32:1583–94. 103. Fernandez E, Luedeling E, Behrend D et al. Mild water stress makes apple buds more likely to flower and more responsive to artificial forcing — Impacts of an unusually warm and dry summer in Germany. Agronomy 2020;10:1–19. 73 104. Engin H, Ünal A. Examination of flower bud initiation and differentiation in sweet cherry and peach by scanning electron microscope. Turk J Agric For 2007;31:373–9. 105. Foster T, Johnston R, Seleznyova A. A morphological and quantitative characterization of early floral development in apple (Malus × domestica Borkh .). Annals of Botany 2003;92:199– 206. 106. Yarur A, Soto E, Leo G et al. The sweet cherry (Prunus avium) FLOWERING LOCUS T gene is expressed during floral bud determination and can promote flowering in a winter-annual Arabidopsis accession. Plant Reproduction 2016;29:311–22. 107. Mimida N, Ureshino A, Tanaka N et al. Expression patterns of several floral genes during flower initiation in the apical buds of apple (Malus x domestica Borkh.) revealed by in situ hybridization. Plant Cell Rep 2011;30:1485–92. 108. Goeckeritz CZ, Gottschalk C, van Nocker S et al. Malus species with diverse bloom times exhibit variable rates of floral development. J Amer Soc Hort Sci 2023;148:64–73. 109. Wang X, Gong JZ, Li QJ et al. Floral organogenesis of Prunus laurocerasus and P. serotina and its significance for the systematics of the genus and androecium diversity in rosaceae. Botany 2019;97:71–84. 110. Esau K. Plant Anatomy. Second. USA: John Wiley & Sons, Inc., 1965. 111. Mauseth JD. Plant Anatomy. First. Crowley A, Williams R, Faust L, et al. (eds.). Menlo Park, CA: Benjamin/Cummings Publishing Company, Inc, 1988. 112. Bierhorst DW. Morphology of Vascular Plants. First. Giles NH, Torrey JG (eds.). New York, USA: The Macmillan Company, 1971. 113. Foster A, Gifford EMJr. Comparative Morphology of Vascular Plants. Second. USA: W. H. Freeman and Company, 1974. 114. Fahn A. Plant Anatomy. Fourth. Great Britain: Pergamon Press, 1990. 115. Endress PK. Patterns of angiospermy development before carpel sealing across living angiosperms: diversity, and morphological and systematic aspects. Bot J Linn Soc 2015;178:556–91. 116. Becker A. A molecular update on the origin of the carpel. Curr Opin Plant Biol 2020;53:15–22. 117. Javier F, Noriega X, Rubio S. Hydrogen peroxide increases during endodormancy and decreases during budbreak in grapevine (Vitis vinifera L.) Buds. 2021. 118. Cheeseman JM. Hydrogen peroxide concentrations in leaves under natural conditions. J Exp Bot 2006;57:2435–44. 119. Tsukamoto T, Hauck NR, Tao R et al. Molecular and genetic analyses of four nonfunctional S haplotype variants derived from a common ancestral S haplotype identified in sour cherry (Prunus cerasus L.). Genetics 2010;184:411–27. 74 120. Beaver JA, Iezzoni AF. Allozyme inheritance in tetraploid sour cherry (Prunus cerasus L.). J Am Soc Hortic 1993;118:873–7. 121. Dirlewanger E, Graziano E, Joobeur T et al. Comparative mapping and marker-assisted selection in Rosaceae fruit crops. PNAS 2004;101:9891–6. 122. Wahl V, Brand LH, Guo Y-L et al. The FANTASTIC FOUR proteins influence shoot meristem size in Arabidopsis thaliana. BMC Plant Biol 2010;10:285. 123. Taiz L, Zeiger E. Plant Physiology. 5th edition. Massachusetts, USA: Sinauer Associates, Inc, 2010. 124. Mantegazza O, Gregis V, Mendes MA et al. Analysis of the arabidopsis REM gene family predicts functions during flower development. Ann Bot 2014;114:1507–15. 125. Romeu JF, Monforte AJ, Sánchez G et al. Quantitative trait loci affecting reproductive phenology in peach. BMC Plant Biol 2014;14, DOI: 10.1186/1471-2229-14-52. 126. Tholl D. Terpene synthases and the regulation, diversity and biological roles of terpene metabolism. Curr Opin Plant Biol 2006;9:297–304. 127. Pichersky E, Raguso RA. Why do plants produce so many terpenoid compounds? New Phytol 2018;220:692–702. 128. Jia Q, Brown R, Köllner TG et al. Origin and early evolution of the plant terpene synthase family. Proc Natl Acad Sci USA 2022;119:e2100361119. 129. Stavang JA, Gallego-Bartolomé J, Gómez MD et al. Hormonal regulation of temperature- induced growth in Arabidopsis. TPJ 2009;60:589–601. 130. Romanel EAC, Schrago CG, Couñago RM et al. Evolution of the B3 DNA binding superfamily: New insights into REM family gene diversification. Ingvarsson PK (ed.). PLOS ONE 2009;4:e5791. 131. Franco-Zorrilla JM, Cubas P, Jarillo JA et al. AtREM1, a member of a new family of B3 domain-containing genes, is preferentially expressed in reproductive meristems. Plant Physiol 2002;128:418–27. 132. Caselli F, Beretta VM, Mantegazza O et al. REM34 and REM35 control female and male gametophyte development in Arabidopsis thaliana. Front Plant Sci 2019;10:1351. 133. Chandler J, Wilson A, Dean C. Arabidopsis mutants showing an altered response to vernalization. TPJ 1996;10:637–44. 134. King GJ, Chanson AH, McCallum EJ et al. The Arabidopsis B3 domain protein VERNALIZATION1 (VRN1) Is involved in processes essential for development, with structural and mutational studies revealing its DNA-binding surface. J Biol Chem 2013;288:3198–207. 135. Richter R, Kinoshita A, Vincent C et al. Floral regulators FLC and SOC1 directly regulate expression of the B3-type transcription factor TARGET OF FLC AND SVP 1 at the Arabidopsis shoot apex via antagonistic chromatin modifications. Yu H (ed.). PLOS Genet 2019;15:e1008065. 75 136. Yu Y, Qiao L, Chen J et al. Arabidopsis REM16 acts as a B3 domain transcription factor to promote flowering time via directly binding to the promoters of SOC1 and FT. TPJ 2020;103:1386–98. 137. Calle A, Grimplet J, Dantec LL et al. Identification and characterization of DAMs mutations associated with early blooming in sweet cherry, and validation of DNA-based markers for selection. Front Plant Sci 2021;12, DOI: 10.3389/fpls.2021.621491. 138. Moser M, Asquini E, Miolli GV et al. The MADS-Box Gene MdDAM1 Controls Growth Cessation and Bud Dormancy in Apple. Frontiers in Plant Science 2020;11:1–13. 139. Balcerowicz M, Shetty KN, Jones AM. Fluorescent biosensors illuminating plant hormone research. Plant Physiol 2021;187:590–602. 140. Michigan State University. Enviroweather: Weather-based pest, natural resources, and production management tools. MSU Enviroweather 2023. 141. R Core Team. R: A language and environment for statistical computing. 2022. 142. Pinheiro J, Bates D. nlme: Linear and nonlinear mixed effects models. 2023. 143. Chakraborty S, Hill AL, Shirsekar G et al. Quantification of hydrogen peroxide in plant tissues using Amplex Red. Methods 2016;109:105–13. 144. Noctor G, Mhamdi A, Foyer CH. Oxidative stress and antioxidative systems: recipes for successful data collection and interpretation. Methods in Oxidative Stress 2016;39:1140–60. 145. Ngo T, Lenhoff H. A sensitive and versatile chromogenic assay for peroxidase and peroxidase-coupled reactions. Anal Chem 1980;105:389–97. 146. Dunn N, Unni D, Diesh C et al. Apollo: Democratizing genome annotation. PLOS Comput Biol 2019;15:e1006790. 147. Andrews S. FastQC: A quality control tool for high throughput sequence data. [Online] 2010, DOI: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/. 148. Li H, Durbin R. Fast and accurate short read alignment with Burrows – Wheeler transform. J Bioinform 2009;25:1754–60. 149. PicardTools. Broad Institute, GitHub Repository 2018;v 2.18.1. 150. Li H, Handsaker B, Wysoker A et al. The Sequence Alignment/Map format and SAMtools. J Bioinform 2009;25:2078–9. 151. Wickham H. Ggplot2: Elegant Graphics for Data Analysis. New York: Springer-Verlag 152. Mckenna A, Hanna M, Banks E et al. The Genome Analysis Toolkit: A MapReduce framework for analyzing next-generation DNA sequencing data. Genome Research 2010;20:1297–303. 76 153. Cingolani P, Platts A, Wang LL et al. A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w 1118; iso-2; iso-3. Fly 2012;6:80–92. 154. Rausch T, Zichner T, Schlattl A et al. DELLY: structural variant discovery by integrated paired-end and split-read analysis. Bioinform 2012;28:i333–9. 155. Thorvaldsdóttir H, Robinson JT, Mesirov JP. Integrative Genomics Viewer (IGV): High- performance genomics data visualization and exploration. Brief Bioinform 2013;14:178–92. 156. Gasic K, Hernandez A, Korban SS. RNA extraction from different apple tissues rich in polyphenols and polysaccharides for cDNA library construction. Plant Mol Biol Rep 2004;22:437–8. 157. Bolger AM, Lohse M, Usadel B. Trimmomatic: A flexible trimmer for Illumina sequence data. J Bioinform 2014;30:2114–20. 158. Dobin A, Davis CA, Schlesinger F et al. STAR: Ultrafast universal RNA-seq aligner. J Bioinform 2013;29:15–21. 159. Kovaka S, Zimin AV, Pertea GM et al. Transcriptome assembly from long-read RNA-seq alignments with StringTie2. Genome Biol 2019;20:278. 160. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 2014;15, DOI: 10.1186/s13059-014-0550-8. 161. Langfelder P, Horvath S. WGCNA: An R package for weighted correlation network analysis. BMC Bioinform 2008;9, DOI: 10.1186/1471-2105-9-559. 162. Boutros P. Generate High-Resolution Venn and Euler Plots. 2022. 163. Camacho C, Coulouris G, Avagyan V et al. BLAST+: Architecture and applications. BMC Bioinform 2009;10, DOI: 10.1186/1471-2105-10-421. 164. Quevillon E, Silventoinen V, Pillai S et al. InterProScan: Protein domains identifier. Nucleic Acids Res 2005;33:W116–20. 165. Carbon S, Douglass E, Good BM et al. The Gene Ontology resource: Enriching a GOld mine. Nucleic Acids Res 2021;49:D325–34. 166. Alexa A, Rahnenfuhrer J. topGO: Enrichment Analysis for Gene Ontology. 2023. 167. Mansfeld BN, Colle M, Zhang C et al. Developmentally regulated activation of defense allows for rapid inhibition of infection in age-related resistance to Phytophthora capsici in cucumber fruit. BMC Genomics 2020;21, DOI: 10.1186/s12864-020-07040-9. 168. Fisher RA. Statistical Methods for Research Workers. 12th ed. Edinburgh: Oliver and Boyd, 1954. 77 Supplementary Figures APPENDIX Figure S4.1 Sections of early-blooming and late-blooming siblings on 14 June 2018. A) Early 2, B) Early 3, C) Late 2, D) Late 3. All apices were vegetative, showing no obvious morphological differences between bloom groups. Scale bars are 200 μm unless stated otherwise. 78 Figure S4.2 Sections of early-blooming and late-blooming siblings on 8 August 2018. A) Early 2, B) Early 3, C) Late 2, D) Late 3. By August, early-bloomers were further ahead than their two late siblings – the former showed development of 3 – 4 floral whorls while sepal primordia was the only whorl visible for the latter. Scale bars are 200 μm unless stated otherwise. 79 Figure S4.3 Sections of early-blooming and late-blooming siblings on 20 September 2018. A) Early 1, B) Early 2, C) Early 3, D) Late 1, E) Late 2, F) Late 3. As summer fades into fall, the early-blooming individuals remain ahead in their floral development compared to their late- blooming siblings. The carpel is closed in the former but just beginning to close (or initiate, in the case of Late 2) in the latter. Scale bars are 200 μm unless stated otherwise. 80 Figure S4.4 Sections of early-blooming and late-blooming siblings on 30 October 2018. A) Early 1, B) Early 2, C) Early 3, D) Late 2, E) Late 3. As flowers enter dormancy, some early- bloomers show signs of ovule emergence while late-bloomers’ carpels have a rudimentary carpel margin meristem (CMM) at best. Scale bars are 200 μm unless stated otherwise. 81 Figure S4.5 Sections of early-blooming and late-blooming siblings on 27 November 2018. A) Early 2, B) Early 3, C) Late 1, D) Late 2, E) Late 3. Like the previous date, early-bloomers remain further advanced during the dormant season. Scale bars are 200 μm unless stated otherwise. 82 Figure S4.6 Sections of early-blooming and late-blooming siblings on 27 December 2018. A) Early 1, B) Early 2, C) Early 3, D) Late 1, E) Late 2, F) Late 3. A month later, development has progressed slightly for both bloom groups; more early-bloomers show ovule emergence and more late-bloomers’ carpel margin meristem is visible. Dashed lines indicate the rest of the pistil is present, but not visible in the photo due to the angle of sectioning. Scale bars are 200 μm unless stated otherwise. 83 Figure S4.7 Sections of early-blooming and late-blooming siblings on 21 February 2019. A) Early 1, B) Early 2, C) Early 3, D) Late 1, E) Late 2, F) Late 3. Development is similar to that seen in December, and the return of warm temperatures will facilitate progression. Scale bars are 200 μm unless stated otherwise. 84 Figure S4.8 Sections of early-blooming and late-blooming siblings on 18 March 2019. A) Early 1, B) Early 2, C) Early 3, D) Late 1, E) Late 2, F) Late 3. A cold spring resulted in little advancement since February for both bloom groups. Scale bars are 200 μm unless stated otherwise. 85 Figure S4.9 Sections of early-blooming and late-blooming siblings on 5 April 2019. A) Early 1, B) Early 2, C) Early 3, D) Late 1, E) Late 2, F) Late 3. By April, growing degree days (GDDs) finally begin to accumulate. Early-bloomers’ pollen mother cells initiate microsporogenesis; the ovule is a large mass of cells emerging from the placental wall. The ovule primordia of late- bloomers are just emerging. Scale bars are 200 μm unless stated otherwise. 86 Figure S4.10 Sections of early-blooming and late-blooming siblings on 15 April 2019. A) Early 1, B) Early 2, C) Early 3, D) Late 1, E) Late 2, F) Late 3. The early-bloomers are still further ahead. The flowers from this group are in the midst of or just finishing microsporogenesis. The integuments on the ovule are clearly visible, and in some cases – so is the megaspore mother cell. Meanwhile, the late-bloomers are still preparing for microsporogenesis, and the ovule’s integuments are not present or they are just emerging. Scale bars are 200 μm unless stated otherwise. 87 Figure S4.11 Association between H2O2 content, a reactive oxygen species (ROS), and the temperature at the time of collection. The equation of the trendline and r2 value is given at the top right of the plot. A weak association was identified; however, since ROS levels declined with chilling, by spring (when warmer temperatures return), all trees have sensed sufficient chill for a high rate of growth. Therefore, this trend may be a result of seasonal progression. More research is required to understand these data. 88 Figure S4.12 The hierarchal clustering of gene expression and module groupings indicate strong clusters were identified during the WGCNA analysis. 19 modules were identified; genes in the hotpink module most effectively separated the bloom time groups. Portions of the hotpink module are noted with an *. 89 Supplementary Tables b l o o m t i m e s , i d e n t i f i c a t i o n n u m b e r f o r t h e i r b l o o m g r o u p , a n d b l o o m d a t e s f o r s e v e r a l y e a r s . T a b l e S 4 . 1 T h e i d e n t i t i e s o f e a c h t r e e i n t h e r e s e a r c h o r c h a r d u s e d f o r t h e p r e s e n t s t u d y , t h e i r k d o s a g e , r e l a t i v e 90 Table S4.2 Dates flowers were collected for pistil measurements. The number of trees representing each bloom group is given in the first two columns and corresponding growing degree days (GDDs) for each respective collection date is given in the last column. If a number is within ( ), this indicates how many trees in that bloom group had bloomed by that date. At least 5 – 10 pistils were measured for every tree and collection. 91 Fixation/Sectioning/Histology # trees w/ no k (Early) # trees w/ 2k (Late) replicates per tree 2 2 3 3 2 3 3 3 3 3 3 3 3 1 1 3 3 3 3 2 4 4 4 2 2 3 2 3 3 3 3 3 3 3 3 3 1 1 3 3 3 3 2 4 4 4 2 2 2 2 2 2 2 (one late tree had only 1) 2 2 2 3 2 3 3 3 5 5 5 (one early tree had 7) 5 2 2 2 2 Date 6/14 8/8 9/20 10/30* 11/28* 12/27* 2/21* 3/18* 4/5 4/15* 6/7 6/12** 6/29** 7/7 7/11 7/18** 8/8 8/22 9/7 11/16* 2/22* 4/6* 5/2* Year GDDs since Jan 1st of starting year 2018 - 19 2018 - 19 2018 - 19 2018 - 19 2018 - 19 2018 - 19 2018 - 19 2018 - 19 2018 - 19 2018 - 19 2019 - 20 2019 - 20 2019 - 20 2019 - 20 2019 - 20 2019 - 20 2019 - 20 2019 - 20 2019 - 20 2019 - 20 2019 - 20 2019 - 20 2019 - 20 684 1639 2293 2538 2558 2559 2566 2572 2593 2628 472 539 780 937 1007 1135 1493 1717 1926 2423 2456 2516 2616 Table S4.3 Dates flowers were collected for histological assessments. The number of trees representing each bloom group is given in the first two columns and corresponding growing degree days (GDDs) for each respective collection date is given in the last column. Bold date* = has corresponding pistil measurements; bold date** = RNA-sequencing was done for these dates. 92 Table S4.4 Dates branches with floral buds were collected for the forcing experiments. The number of trees representing each bloom group is given in the first two columns and corresponding growing degree days (GDDs) for each respective collection date is given in the last column. Chill portions according to the dynamic model is also provided for both years. Bold date* = has corresponding pistil measurements. Table S4.5 Dates flowers were collected for the H2O2 measurements. The number of trees representing each bloom group is given in the first two columns and corresponding growing degree days (GDDs) for each respective collection date is given in the last column. Chill portions according to the dynamic model is also provided for both years. Bold date* = has corresponding pistil measurements. 93