Genomic evidence for supergene control of Darwin’s “complex marriage arrangement” – the tristylous floral polymorphism

Genomic evidence for supergene control of Darwin’s “complex marriage arrangement” – the tristylous floral polymorphism | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Article Genomic evidence for supergene control of Darwin’s “complex marriage arrangement” – the tristylous floral polymorphism Haoran Xue, Yunchen Gong, Stephen Wright, Spencer Barrett This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-4232430/v1 This work is licensed under a CC BY 4.0 License Status: Posted Version 1 posted You are reading this latest preprint version Abstract Tristyly is a polymorphism characterized by three flower morphs with reciprocal stigma and anther heights controlled by two epistatically interacting diallelic loci ( S and M ), hypothesized to be supergenes. Chromosome-level genome assemblies of Eichhornia paniculata identified the S- and M- loci. The S -locus is a supergene consisting of two divergent alleles: The S -allele (2.51Mb) with three S -allele specific genes hemizygous in most S-morph plants and the s -allele (596kb) with five s -allele specific genes. Two of the S -allele specific genes, LAZY1-S and HRGP-S, were specifically expressed in styles and stamens, respectively, making them tristyly candidate genes. The M -locus contained one gene (LAZY1-M), homologous to LAZY1-S, that was present in the M -allele but absent from the m -allele. Estimates of allele ages are consistent with the prediction that the S- locus evolved before the M -locus. Re-use of the same gene family highlights the potential role of gene duplication in the evolution of epistatic multilocus polymorphisms. Biological sciences/Evolution Biological sciences/Plant sciences/Plant evolution Biological sciences/Plant sciences/Plant genetics tristyly plants genomics supergene Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Introduction Flowering plants possess a remarkable diversity of floral strategies promoting cross-pollination and outcrossed mating. This variation attracted the attention of Charles Darwin who devoted most of his classic book “ The Different Forms of Flowers on Plants of the Same Species ” 1 to the floral polymorphism heterostyly in which populations are characterized by two (distyly) or three (tristyly) floral morphs, which possess a reciprocal arrangement of stigma and anther heights 2 . Darwin described heterostyly as a ‘complex marriage arrangement’ promoting cross-pollination between the floral morphs by animal pollinators owing to segregated pollen placement on their bodies (see Fig. 1 in Ref. 3). Contemporary experimental evidence generally supports this hypothesis also finding that heterostyly limits sexual interference between female and male function thus reducing pollen wastage 3-4 . Heterostyly is reported from 28 animal-pollinated angiosperm families 5-6 and represents a remarkable example of adaptive convergence in form and function. Controlled crosses on over a dozen unrelated species have established that distyly is governed by a single diallelic Mendelian locus ( S ) with dominance, whereas tristyly was controlled by two diallelic epistatically interacting loci ( S , M ), both with dominance 7 (Fig. 1). To resolve the discrepancy between the complexity of the heterostylous syndrome and its relatively simple inheritance, it has been proposed that heterostyly loci are in fact supergenes comprised of clusters of functionally related genes 8 . Evidence indicates that these genes are tightly linked to each other due to suppression of recombination and are inherited as a single unit 7 . Supergenes are also known to control several other complex polymorphic traits, including Batesian mimicry in Heliconius butterflies, social organizations of fire ants, and plumage variation and behaviour in ruff 9-11 . Of particular importance for understanding the evolution of supergenes is to determine their genetic architecture and the mechanisms causing suppression of recombination. Identification of the causal loci can also provide important insights into the age of the alleles, to assess the evidence for long-term balanced polymorphisms maintaining variation. Recent efforts have been made to elucidate the genetic architecture of heterostylous supergenes in a handful of distylous systems (reviewed in 12-13 ). Significantly, genomic studies on Fagopyrum, Gelsemium , Linum , Nymphoides , Primula, and Turnera , each belonging to a different angiosperm lineage, have demonstrated that distyly is controlled by several hemizygous genes present in the short-styled morph (S-morph) but absent from the long-styled morph (L-morph) 14-20 . For example, in Primula , the S -locus is a hemizygous region containing five genes controlling the expression of different components of the distylous polymorphism 14-15, 21 . All molecular investigations of heterostyly have focused on distylous species and the supergene hypothesis has yet to be investigated in any of the six angiosperm families in which the more complex tristyly occurs 22-23 . Charlesworth 24 suggested that alternatives to the supergene hypothesis should be considered for tristyly because of the absence of empirical evidence for supergenes. Thus, the genetic basis and evolutionary assembly of this two-locus epistatically interacting genetic system remains unclear. Here, we generate chromosome-level de novo genome assemblies with long-read sequences and integrate them with data from population genomic and transcriptomic analysis to investigate the genetic architecture of tristyly in Eichhornia paniculata , an annual aquatic native to the Neotropics. Earlier molecular studies demonstrated that the S- and M -loci in this species are closely linked (separated by a map distance of 2.7 cM), and that the M -locus is located on linkage group 5 in a region containing 204 genes 25 , but the genomic location of the S -locus was not investigated, and they did not resolve whether the M -locus is a supergene or merely contains one pleiotropic gene 12 . To investigate this issue and provide novel insight into the molecular basis of tristyly we addressed the following questions: 1) What is the genomic architecture of the S- and M- loci and are they indeed supergenes? 2) What are the likely candidate genes at the two loci controlling different components of tristyly? 3) What is the evolutionary timing of the origins of the S- and M- loci? We note that in the only theoretical analysis of the evolution of tristyly, the origin of the S -locus was inferred to have preceded the origin of the M -locus 24 . Our results are consistent with this stepwise origin of the two loci, with recurrent gene duplication from the same ancestral gene playing a key role in the formation of the epistatic interaction between the two loci. Results Comparisons of chromosome-level genome assemblies We generated chromosome-level genome assemblies from a heterozygous S-morph individual ( SsMm ) and an L-morph individual ( ssmm ) combining long-read PacBio sequencing with proximity ligation (Omni-C and Hi-C) scaffolding (see methods and Fig. S1). In bothassemblies, most of the two genomes were assembled into eight major chromosomal scaffolds (Table 1). Comparisons of the two assemblies demonstrated large-scale concordance, with small-scale differences likely reflecting inversions and/or small assembly errors (Fig. S2a). There was clear evidence for a recent whole genome duplication event in which each chromosome showed strong synteny with one or more additional chromosomes. Based on the distribution of synonymous substitution rate ( Ks ) values between paralogs, the most recent whole genome duplication event has a mode of Ks =0.38 (Fig. S2b). Assuming a substitution rate of 7 x 10 -9 per year 26 and a whole genome duplication event, this would correspond to two sets of homeologous chromosomes that diverged less than roughly 30 million years ago. Since the base chromosome number in Pontederiaceae was inferred to be n =8 27 , the duplication event was likely shared across the family, with additional lineage-specific duplications due to additional rounds of polyploidy occurring in other lineages of the family 28 . Unless otherwise noted, our analysis below uses our SsMm assembly as the reference genome. The Marey maps (plots of the linkage map against the physical map) indicated that there was a large (likely pericentromeric) region of low recombination in the middle of each chromosome (Fig. 2a). In chromosomes 3 and 7, the region extended to one end of the chromosome, likely corresponding to submetacentric or acrocentric chromosomes, whereas the remaining chromosomes appeared to be more metacentric. As expected, the regions of reduced recombination were associated with a high repeat density (particularly LTR retrotransposons), low gene density, and low nucleotide diversity, whereas the more highly recombining chromosome tips were gene-dense, comprised of fewer repeats, and has higher nucleotide diversity (Fig. 2b-f, Fig. S3). Both tristyly loci occur in a large low-recombining region Using a genome-wide association study (GWAS) we identified 4,640 SNPs that were significantly associated with the S-morph ( S ---) versus non-S-morphs ( ss --, i.e. L- and M-morphs) ( p < 2.52397779 x 10 -8 , Fig. 3a). 96.8% of these SNPs occurred in a very large 63.1 Mb region on chromosome 2, embedded in an extensive region of reduced recombination on this chromosome (Fig. 3a). As expected from earlier genetic work 25 , the M- locus region is partially linked to the S -locus. We found 4,541 SNPs significantly associated with the M- ( ssM- ) versus L- ( ssmm ) morphs ( p < 2.23826411 x 10 -8 , Fig. 3b). 98.0% of these SNPs were in a 35.1 Mb region on chromosome 2, at the other end of the region of reduced recombination from the S -locus region. Thus, both the S- and M -loci governing tristyly occur in a large, putatively pericentromeric region of low recombination on chromosome 2 causing extensive and significant SNP associations across many Mb of sequence. The S -locus is a supergene consisting of divergent alleles with hemizygous genes The S- and M- locus associated regions we identified were both extremely large with very low rates of recombination. To further narrow down possible candidate genes and to investigate the potential hemizygous regions associated with tristyly, we identified regions with significant differences in read coverage between the floral morphs. We identified 50 consecutive 50-kb windows that were specific to the S -allele (total length: 2.5 Mb; genome coordinates: Chr2:114,300,000-116,800,000), with read coverage significantly higher in 20 S-morph ( S ---) plants than 40 non-S-morph ( ss --) plants (criteria: the ratio between the average read coverage of non-S-morph and S-morph plants was lower than 0.4, with p -value of the t- test < 9.73217066 × 10 -7 , Fig. 4a). Downstream adjacent to this 2.5 Mb region were 13 consecutive 50-kb windows (total length: 650 kb; genome coordinates: Chr2:116,800,000-117,450,000) that were specific to the s -allele and had read coverage significantly lower in 20 S-morph ( S ---) plants than 40 non-S-morph ( ss --) plants (criteria: the ratio between the average read coverage of non-S-morph and S-morph plants was higher than 1.6, with the p -value of the t -test < 9.73217066 × 10 -7 , Fig. 4a). We hypothesized that given the high divergence between alleles, our collapsed genome assembly may have created a chimeric assembly in this region, leaving the s -allele sequence adjacent to the S -allele region in our assembled scaffold and contributing to the lack of a signal of high diversity when mapping to this reference genome. To investigate this hypothesis further, we compared our L-morph assembly ( ssmm ) to this assembly in this region using whole genome alignments. As expected, Anchorwave 29 alignments revealed the presence of a large 2.57 Mb deletion in the ssmm assembly relative to the SsMm assembly (Fig. S4), whereas the adjacent region showed homology. Thus, we conclude that the 2.5 Mb region roughly corresponds to the S -allele and the 650 kb region roughly corresponds to the s -allele. We refined the boundaries of the two alleles using the coordinates of the closest genes in their flanking regions (Fig. 4a). Two genes (TIE/SPL and MPC2) were shared between the S- and s -alleles. The order of the two genes are opposite in the two alleles. In the S -allele, both genes are on the negative strand, whereas in the s –allele they are on the positive strand. This finding suggests the existence of an inversion containing the two genes (Fig. 4a). The rest of the genes were allele specific: three genes (CYCD6;1-S, HRGP-S, and LAZY1-S) were specific to the S -allele, while five genes were specific to the s -allele (Fig. 4a). LAZY1 is a gene family that regulates the distribution of auxin 30 , HRGP is a superfamily of genes that forms covalently cross-linked networks in primary cell walls 31 , and CYCD6;1 is a gene that regulates the formative division of cells 32 . The M -locus contains one hemizygous gene We identified five consecutive 50-kb windows (total length: 250 kb; genome coordinates: Chr2:13,000,000-13,250,000) that had read coverage significantly higher in the 20 M-morph ( ssM- ) plants than in the 20 L-morph ( ssmm ) plants (criteria: the ratio between the average read coverage of M-morph and L-morph plants was higher than 1.6, with the p -value of the t -test < 9.73217066 × 10 -7 , Fig. 4b). There were nine coding genes in these five windows. However, among them only one gene was missing from the ssmm genome assembly suggesting that it is absent from the m -allele (Fig. 4b). Strikingly, this absent gene, LAZY1-M, also shows strong similarity to genes in the LAZY1 family, raising the possibility that this locus acts competitively in the same pathway as LAZY1-S. Consistent with the SNP association mapping, the hemizygous M -locus region was located on the other side of the low-recombing region on chromosome 2 (Fig. 3e). Genetic degeneration at the S -allele In species with supergenes, signatures of genetic degeneration, including gene loss and accumulation of TEs, are expected for the dominant haplotypes of the supergenes, because dominant haplotypes are only present in heterozygotes and thus do not exhibit recombination 33 . The S -allele specific gene loss (the five s -allele specific genes) was in accord with this prediction. In addition, we found that the density of repetitive elements in the S - and M -alleles was significantly higher than the flanking regions ( p = 0.00001985 and p = 0.04731, respectively, Wilcoxon signed-rank test) (Fig. 4f and h), whereas the differences between the s -allele and the flanking regions were not significantly different ( p = 0.847, Wilcoxon signed-rank test) (Fig. 4g). The expression pattern of hemizygous genes Among the three S -allele specific genes, LAZY1-S and HRGP-S were expressed in floral organs of the S-morph but not non-S-morph plants. LAZY1-S was expressed in styles but not stamens, whereas in contrast HRGP-S was expressed in stamens but not styles (Fig. 5a). CYCD6;1-S has low expression levels in both the styles and stamens of the S-morph (Fig. 5a). The M -allele specific gene LAZY1-M was expressed in both styles and stamens (Fig. 5b). Genes that were differentially expressed in different morphs For genes that were differentially expressed genome-wide between floral morphs (see Methods and Table S1 for details of comparisons), three sets of genes showed functional enrichment: ( 1) genes that were upregulated in S-morph styles (compared with non-S-morph styles) in early buds (enriched in three non-overlapping functions), (2) genes that were upregulated in S-morph styles (compared with non-S-morph styles) in late buds (enriched in one function), and (3) genes that were upregulated in M-morph styles (compared with L-morph styles) in late buds (enriched in one function) (Table S2). However, none of the five functions inferred had a direct connection to growth or flower organs, suggesting that only a small number of genes downstream to the S - and M -locus candidate genes are involved in the phenotypic differences between the floral morphs. Other differentially expressed genes may simply be incidental associations. The origins of the S and M loci To further understand the origins of the S- and M- loci, we identified the closest paralogs of the S - and M -specific genes. We identified one gene that was paralogous to LAZY1-S and LAZY1-M. Hereafter, we refer to it as “LAZY1-P1”. LAZY1-P1 was located 223 kb downstream to LAZY1-M (Fig. 4b). The closest paralog of HRGP-S was found on the same chromosome but not in the same region (Fig. 4c). Two paralogs of CYCD6;1-S were annotated in both genome assemblies, and they are on chromosomes 3 and 5, respectively (hereafter CYCD6;1-Chr3 and CYCD6;1-Chr5). CYCD6;1-Chr3 shared higher sequence similarity with CYCD6;1-S (Fig. S5) and was located in a homologous region between chromosomes 2 and 3 (Fig. S2a), which was probably a result of the whole genome duplication estimated to have happened < 30 million years ago, and CYCD6;1-Chr5 is likely the ancestral copy of both CYCD6;1-S and CYCD6;1-Chr3. The distance between the three ancestral genes therefore suggests that the S -locus was formed via stepwise duplication rather than segmental duplication (Fig. 4c). Based on our gene trees and assuming that the hemizygous loci represent the derived gene copies, LAZY1-S duplicated from LAZY1-P1 first, followed by subsequent gene duplication from LAZY1-P1 to form LAZY1-M. Phylogenetic dating with BEAST v2.7.6 34 suggested that the divergence between LAZY1-S and LAZY-P1 occurred approximately 52.4 million years ago (95% confidence interval: 44.9-59.7) and that the subsequent divergence between LAZY1-M and LAZY-P1 was around 43.0 million years ago (95% confidence interval: 36.4-50.0) (Fig. 6a). Discussion We show that in tristylous Eichhornia paniculata , both the S- and M -loci are located in a low recombining region of Chromosome 2. The S -locus is comprised of two divergent alleles: the 2.51 Mb S -allele with three S -allele specific genes that are hemizygous in most S-morph plants, and the 596 kb s -allele with five s -allele specific genes. In the S -allele we identified two strong tristyly candidate genes: LAZY1-S and HRGP-S. LAZY1-S was expressed in styles but not stamens suggesting that it may inhibit style elongation in the S-morph by decreasing auxin concentration in styles. HRGP-S was expressed in stamens but not styles suggesting that it may be involved in filament elongation or possibly pollen function. Our evidence for the S -locus therefore strongly supports its status as a supergene, and its origin in a large region of low recombination may have been important during the initial origins of these genes via gene duplication from more physically distant regions (Fig. 4 ). We found that the M -locus contained only one hemizygous gene – LAZY1-M – and this gene was hemizygous in -- Mm plants. Strikingly, this gene is in the same gene family as LAZY1-S suggesting the possible role of gene duplication and neofunctionalization in the origin of the epistatic interaction between the two loci. Unlike LAZY1-S, which was only expressed in styles, LAZY1-M was expressed in both styles and stamens and with a higher expression level in stamens. This suggests that while LAZY1-S may only be involved in the shortening of styles in the S-morph, LAZY1-M may also regulate auxin distribution and organ growth in both the styles and stamens of the M-morph. Given that LAZY-M is the only M -allele-specific hemizygous gene, the M-locus may not represent a supergene at this stage in its evolutionary history but simply contains LAZY1-M, a pleiotropic modifier. LAZY1-M was also expressed in the stamens but not the styles of the S-morph (Fig. 5 b), suggesting that LAZY1-S and LAZY1-M may compete for the same transcriptional activator and that the style expression of LAZY1-M is shut down if LAZY1-S is expressed. The expression level of LAZY1-S in S-morph styles is also much higher than LAZY1-M in M-morph styles. For example, at early development stages the average normalized count of LAZY1-M in samples of M-morph styles is 0.02% of that of LAZY1-S in S-morph styles. At late development stages, the ratio is 1.9%, indicating that the promoter region of LAZY1-S may have stronger affinity to the transcriptional activator that they compete for, potentially driving differences in style length in the S- and M-morphs through different degrees of auxin suppression. In a theoretical treatment of the evolutionary assembly of tristyly 24 , Charlesworth predicted that the M-morph originated last spreading in dimorphic populations of the L- and S-morph and eventually giving rise to the tristylous polymorphism. Our results on the relative ages of LAZY1-S and LAZY1-M are consistent with the hypothesis that the M -locus evolved after the S -locus providing empirical support for this prediction. Our results also provide evidence suggesting a major role for ‘re-use’ through recurrent duplication of the same gene (LAZY1). Recruitment of this gene explains the origin of the M -locus and birth of the third floral morph in the tristylous system (Fig. 6 b) resulting in the morphological and mating complexity that characterizes Darwin’s ‘complex marriage arrangement’. Online Methods Study system . Plants used in our study originated from seed collected in wild populations of E. paniculata in N.E. Brazil, Cuba and Mexico. The localities of plant material and details of the five populations used in our study are presented in Table S4. Information on the reproductive ecology, mating systems and population genetics of E. paniculata are available in Ref. 35 and references cited therein. We germinated seeds using PRO-MIX Premium Potting Mix (Premier Tech, Rivière-du-Loup, QC, Canada) and seedlings were transferred to individual 7 cm plastic pots filled with the same soil mix and placed in flooded trays in a glasshouse at the University of Toronto (Toronto, Canada) with temperatures of 25–28°C (day) and 18–21°C (night) and a 16-hour day length. Plants were then grown and flowered regularly for several years facilitating crosses and collection of leaves and floral buds for molecular analyses. Genome assemblies. We conducted two chromosome-level genome assemblies of E. paniculata. The first assembly was obtained from an individual of genotype SsMm (inferred from selfed progeny, Table S3) obtained from an inter-population cross between an S-morph plant ( SsM- ) from N.E. Brazil and an M-morph plant ( ssMm ) from Cuba (see Table S4 for information on the two parents). From this individual, we obtained 79 Gb (105X coverage) PacBio CLR sequences (sequenced by Arizona Genomics Institute, Tucson, AZ, USA), 55 Gb (45X coverage) Hi-C (sequenced by Phase Genomics, Seattle, WA, USA) sequences, and 117 Gb (156X coverage) 10X Chromium sequences (sequenced by The Hospital for Sick Children, Toronto, ON, Canada). The PacBio sequences were assembled with Canu 36 and the assembly was curated with Purge Haplotigs 37 , both by the Canadian Center for Computational Genomics (C3G) (Montreal, QC, Canada). Then the genome was phased with Falcon-Phase 38 , then scaffolded into chromosome-level scaffolds with the Hi-C data, both by Phase Genomics (Seattle, WA, USA). To identify misjoins and characterize the relationship between the physical distance and the genetic distance on each chromosome, we constructed a linkage map with genotyping-by-sequencing data obtained from a backcross pedigree from an L-morph plant and an M-morph plant, detailed in 25 , accessible under Accession no. SRP069136 at the Sequence Read Archive (SRA), ( http://www.ncbi.nlm.nih.gov/sra , accessed 17 January 2023), and the associated BioProject alias PRJNA310331). The reads were trimmed with Scythe 39 and Seqtk ( https://github.com/lh3/seqtk ). We aligned the trimmed reads to the SsMm genome assembly with NextGenMap V0.5.5 40 , and the BAM files were processed with Picard ( http://broadinstitute.github.io/picard ) using the four functions FixMateInformation, CleanSam, SortSam, and AddOrReplaceReadGroups. SAMtools 41 was used to exclude the secondary alignments. We constructed the linkage map with the alignment files and the software program Lep-Map3 42 . Markers were split into eight large linkage groups using a logarithm of the odds (LOD) score limit of 70 and a Theta value of 0.25 for the initial split (SeparateChromosomes2). Most of the remaining single markers were put into these groups using a LOD score limit of 30 and a LOD difference of 5 (JoinSingles2All). We calculated marker order for each linkage group with OrderMarkers2 (with default settings). The first ten linkage groups contained 2,124 markers. We then used Juicebox 43 to correct the misjoins with the linkage map. The genome was polished with the 10X Chromium sequences, using the software Pilon 44 . The repetitive elements were annotated with RepeatModeller V1.0.11 45 and RepeatMasker V4.0.7 46 . The second assembly of E. paniculata used an L-morph ( ssmm ) plant from Mexico (see Table S4) using sequencing and assembly services from Dovetail Genomics (Scotts Valley, CA, USA). 59.6 Gb (~ 50X coverage) PacBio HiFi reads were generated and the reads were assembled with hifiasm 47 . The assembly was then scaffolded with Omni-C data using the HiRise software 48 . We used the linkage map previously generated to concatenate scaffolds that belong to the same linkage group. We also annotated TEs with RepeatModeller V1.0.11 45 and RepeatMasker V4.0.7 46 . For both of our genome assemblies, the distribution of scaffold length was assessed with QUAST V5.2.0 49 , and the completeness of assemblies was assessed with BUSCO (Benchmarking Universal Single-Copy Orthologs) V5.4.4 50 with MetaEuk 51 and the embryophyta_odb10 (2024-01-08) dataset (Table S5). Genome-wide associations of single nucleotide polymorphisms (SNPs) and read coverage (RCs). We extracted genomic DNA from 60 E. paniculata plants (20 of each floral morph) from a population in N.E. Brazil (hereafter “Pop 1”, see Table S4), using the QIAGEN DNeasy Plant Mini Kit (QIAGEN, Germantown, MD, USA). Library preparation and sequencing were conducted by the McGill University and Génome Québec Innovation Centre (Montréal, QC, Canada) using the 150-bp paired-end protocol on Illumina HiSeq 4000. The samples had an average sequencing depth of ~ 8X. The Illumina sequencing reads were trimmed with Scythe 39 and Seqtk ( https://github.com/lh3/seqtk ). We aligned the trimmed reads to the SsMm genome assembly with NextGenMap V0.5.5 52 and the BAM files were processed with Picard ( http://broadinstitute.github.io/picard/ ) using the four functions FixMateInformation, CleanSam, SortSam, and AddOrReplaceReadGroups. We used SAMtools 41 to exclude the secondary alignments. We used freebayes-parallel 53,54 to call SNPs jointly for all samples. We then filtered the VCF file with VCFtools 55 and vcflib2 56 following a published guideline 57 . We removed sites if more than 50% of the individuals had missing data, had a minor allele count of < 6, or had a quality score lower than 30. Genotypes with fewer than 4 reads were treated as "missing" and individuals with more than 50% missing data were not considered further. Then, sites with strong within-individual allelic bias at heterozygous sites (0.01 < AB 0.75), a discrepancy in the properly paired status for reads supporting reference or alternate alleles, a low quality relative to depth, a mean depth that was > 13, a quality score < 30, or that occurred within 5 bp of an indel, were excluded from our analysis. The resulting VCF file had 7,358,165 filtered variant sites. We conducted a Genome-Wide Association Study (GWAS) with the variant sites using GEMMA 58 , where a relatedness matrix and a linear mixed model were used and p -values of variant sites were calculated with the Wald test. To detect hemizygous regions associated with the S- and M -loci, we conducted another GWAS with read-coverage, following the pipeline in Ref. 59. We split the reference genome into 25,688 windows, each of which was 50-kb, and calculated the average depth of each window in individuals using the software mosdepth 60 . For each individual, we normalized the coverage by dividing the coverage of each window with the average coverage across all windows, then multiplied it by 2. Therefore, the expected normalized coverage is 1 for hemizygous regions and 2 for the rest of the genome. We conducted two t -tests for normalized coverage for each window: one for the comparison between S-morph versus non-S-morph plants, and the other for the comparison between the M-morph versus L-morph plants. We identified tristyly allele-specific windows with the following criteria: 1) p value of the t -test was smaller than 9.73217066 × 10 − 7 , i.e. Bonferroni correction with α = 0.05. 2) We set thresholds on the ratio between the average read coverage of plants of different morphs: For S -allele specific windows: The read coverage of a S -allele is at least 4 times of a s -allele, i.e. the ratio between the average read coverage of non-S-morph plants and S-morph plants was lower than 0.4: \(\frac{2RC\left(s\right)}{RC\left(S\right)+RC\left(s\right)}=\frac{2\times 1}{4+1}=0.4\) . For s -allele specific windows: The read coverage of a s -allele is at least 4 times of a S -allele, i.e. the ratio between the average read coverage of non-S-morph plants and S-morph plants was higher than 1.6: \(\frac{2RC\left(s\right)}{RC\left(S\right)+RC\left(s\right)}=\frac{2\times 4}{4+1}=1.6\) For M -allele specific windows: The read coverage of an M -allele is at least 4 times of an m -allele, i.e. the ratio between the average read coverage of non-M-morph plants and M-morph plants was lower than 0.25, assuming all M-morph plants were homozygous: \(\frac{2RC\left(m\right)}{RC\left(M\right)+RC\left(m\right)}=\frac{2\times 1}{4+1}=0.4\) . For m -allele specific windows: The read coverage of an s -allele is at least 4 times of a M -allele, i.e. the ratio between the average read coverage of non-M-morph plants and M-morph plants was higher than 1.6: \(\frac{2RC\left(m\right)}{RC\left(M\right)+RC\left(m\right)}=\frac{2\times 4}{4+1}=1.6\) Only regions with at least two consecutive windows that were allele-specific were considered. As a result, a 2.5 Mb S -allele, a 650 kb s -allele, and a 250 kb M -allele were identified. We split each allele into 10 equal-sized windows (window sizes: S -allele: 250 kb; s -allele: 65 kb; M -allele: 25 kb), calculated the repeat content of each window, then compared it with the 10 windows of the same size flanking each allele (5 windows upstream and 5 windows downstream of the locus) with Wilcoxon rank-sum tests. Differential gene expression analysis and genome annotation with transcriptome data. To infer genes controlling different components of the tristyly syndrome, we conducted a floral organ-specific transcriptomic analysis. We dissected floral buds of the three morphs at two developmental stages: “early”, with buds between 3–5 mm in length, and “late”, with buds between 5–8 mm. We used sterilized forceps under a stereo dissecting microscope to open the floral buds and placed the styles and stamens into 2 mL centrifuge tubes which were then frozen by putting them into liquid nitrogen. We extracted RNA from the frozen tissues with Spectrum™ Plant Total RNA-Kit (Sigma-Aldrich, St. Louis, MO, USA) following the manufacturer’s instructions. For each type of sample, we collected four biological replicates, each containing styles or stamens from 5–27 buds sampled from 5–17 individuals from each of two populations (Pop. 1, 2) originating from N.E. Brazil (Table S4, S6). We prepared sequencing libraries and sequenced on the Illumina HiSeq4000 platform (100-bp paired-end reads) by the McGill University and Génome Québec Innovation Centre (Montréal, QC, Canada). Then we aligned the Illumina sequences to the SsMm genome assembly with HISAT2 61 and assembled the transcripts with StringTie 62 . We conducted functional annotations of the floral organ transcriptome with Mercator4 63 . We used StringTie 62 to estimate the abundance of each transcript in each sample. Then we used the R package DESeq2 64 to conduct differential gene expression analysis and plot the normalized count in each sample for the genes at the S- and M- loci. We conducted eight differential gene expression analyses for each type of morph difference ( S vs. non-S and M vs. L), each tissue (style vs. stamen), and each development stage (early bud vs. late bud). Only genes with FPKM (Fragments Per Kilobase of transcript per Million mapped reads) > 1 and TPM (Transcript Per Million) > 1 in at least one sample were included in the differential gene expression analyses. In each analysis, we identified the genes that are upregulated in each morph with false discovery rate 2 or < -2. We then did a Mercator BIN enrichment analysis 61 with the list of upregulated genes for each comparison. We also annotated the repeat-masked versions of the two genome assemblies with the floral organ transcriptomes and also leaf and floral bud transcriptomes from Ref. 65, available under accession number SRP049636 on the Sequence Read Archive (SRA) ( http://www.ncbi.nlm.nih.gov/sra , accessed 20 January 2023) and the associated BioProject alias is PRJNA266681 ( http://www.ncbi.nlm.nih.gov/bioproject/ , accessed 20 January 2023), and also Ref. 66, available under accession number SRP060405 at the Sequence Read Archive and associated BioProject alias (PRJNA288861) using MAKER v3.01.03 67 and InterProScan 5 68 . This resulted in a total of 35,224 annotated genes for the SsMm genome and 29,592 genes for the ssmm genome. Origin of the S and M loci. According to the annotation of the SsMm and ssmm genome assemblies, there is only one LAZY1 homolog (LAZY1-P1) that exists in both genome assemblies. We then downloaded the coding sequences of LAZY1 in Oryza sativa and Zea mays from NCBI ( https://www.ncbi.nlm.nih.gov/ , accession numbers: NM_001138862.1 and DQ855268.1), and two paralogs of LAZY1 in Musa acuminata 'DH-Pahang' from the website Banana Genome Hub ( https://banana-genome-hub.southgreen.fr/ ). We conducted a multiple sequence alignment of these sequences and the coding sequences of three LAZY1 genes in E. paniculata , using MAFFT v7 69 with L-INS-i method. We then estimated the time of divergence between LAZY1-S and LAZY1-P1 and between LAZY1-M and LAZY1-P1 with BEAST v2.7.6 34 . First, we did model selection of jModelTest 70 using Akaike information criterion (AIC) as the criteria. The GTR + I + Gamma model with nCat = 4 was selected. Then we set the time to the most recent ancestors of Oryza sativa and Zea mays (81.43 ± 0.3 MYA) 71 , Zingiberales and Commelinales (79.8 ± 9.87 MYA) 72 , and Poales and Commelinales (106.7 ± 8.3 MYA) 72 as calibration points. Ten independent Markov Chain Monte Carlo runs, each with 1 x 10 8 generations (first 10% was burn-in) were conducted. We combined the output of the eight runs using LogCombiner v2.7.6 34 with a sample frequency of 5,000. Then we used TreeAnnotator v2.7.6 34 to generate the maximum clade credibility tree, which was visualized in FigTree v1.4.4 ( http://tree.bio.ed.ac.uk/software/figtree/ ) (Fig. 6 ). There are two CYCD6;1 paralogs that exists in both genome assemblies, and they are on chromosomes 3 and 5, respectively (hereafter CYCD6;1-Chr3 and CYCD6;1-Chr5). We then downloaded the coding sequence of CYCD6;1 in O. sativa from NCBI ( https://www.ncbi.nlm.nih.gov/ , accession numbers: XM_015789871.2). We conducted a multiple sequence alignment of the coding sequences of the three CYCD6;1 genes using MAFFT v7 69 with the L-INS-i method. Then we constructed a Neighbor Joining tree with all the conserved sites, Jukes-Cantor substitution model, and bootstrap resampling with 100 samples. (Fig. S). HRGP-S was not annotated by MAKER or InterProScan. We aligned the transcript of HRGP-S against the SsMm and ssmm genome assemblies with the software program blastn 73 , to identify its closest paralog. Synteny analyses, whole genome alignments and structural variants . We used COGE’s SynMap2 algorithm 74 with the MAKER and InterProScan annotations to examine the synteny between the two genome assemblies and identify structural variation. To obtain finer-scale whole genome alignment of syntenic regions to localize putative insertion/deletion differences between the S- and s -alleles and between the M- and m -alleles, we used Anchorwave 29 to align the two genome assemblies, with the “minimap2” option and the “proali” function. Population genomic statistics . We calculated the nucleotide diversity ( π ) of 1-Mb sliding windows for the 60 individuals (20 of each morph) used in the GWAS analyses. We used SAMtools 41 to remove the reads with mapping quality < 10 in the filtered BAM files generated in the section “Genome-wide associations of single nucleotide polymorphisms (SNPs) and read coverage (RCs)”. Then we used ANGSD 75 to calculate the nucleotide diversity ( π ) of the sliding windows. Declarations Acknowledgements This work was supported by Discovery Grants from the Natural Sciences and Engineering Research Council of Canada (NSERC) to S.C.H.B. and S.I.W. and student fellowships from the University of Toronto to H.X. The Canu assembly and Purge Haplotigs pipeline of the SsMm genome was conducted by the Canadian Center for Computational Genomics (C3G), which is a Genomics Technology Platform (GTP) supported by the Canadian Government through Genome Canada. We thank William Cole and Tom Gludovacz for glasshouse maintenance, Joanna Tran, SungA Choe, Ryan Ruan, and Jianuo Feng for help with growing plants, Claire Ellis and Po-Yuan Huang for help with tissue collection, Baharul Choudhury, Iuliia Pelikh, and Yihe Zheng for help with DNA extraction, and Joanna Rifkin, Christian Kappel, Mathias Scharmann, Julia Kreiner, Jingshan Yang, Meng Yuan, Liangjiao Xue, Zoë Humphries, Tyler Kent, and Ting Liu for advice on data analysis. Author contributions H.X., S.I.W., and S.C.H.B. conceived and designed the study and wrote the manuscript. S.C.H.B. obtained the seeds used for the study. H.X. grew the plants, collected the tissues, and extracted DNA and RNA. Y.G. performed the genome annotation. S.I.W. performed the Anchorwave alignment. H.X. performed all remaining analyses. H.X., S.I.W., and S.C.H.B. wrote the manuscript. All authors read and approved the final manuscript. Competing interests None declared. Data Availability Statement The PacBio CLR sequences used to assemble the Eichhornia paniculata SsMm genome are available under Accession no. SRR21038607 at the Sequence Read Archive (SRA) (http://www.ncbi.nlm.nih.gov/sra, accessed 09 February 2024), and the associated BioProject alias is PRJNA863567. The Hi-C data used to assemble the E. paniculata SsMm genome are available under Accession no. SRR21038606 at the Sequence Read Archive (SRA) (http://www.ncbi.nlm.nih.gov/sra, accessed 09 February 2024), and the associated BioProject alias is PRJNA863567. The 10X Chromium data used to polish the E. paniculata SsMm genome are available under Accession no. SRR21038605 at the Sequence Read Archive (SRA) (http://www. ncbi.nlm.nih.gov/sra, accessed 09 February 2024), and the associated BioProject alias is PRJNA863567. The PacBio CCS sequences used to assemble the E. paniculata ssmm genome are available under Accession no. SRRXXXXXXXX at the Sequence Read Archive (SRA) (http://www.ncbi.nlm.nih.gov/sra, accessed XX XX XXXX), and the associated BioProject alias is PRJNAXXXXXX. The Omni-C data used to assemble the E. paniculata ssmm genome are available under Accession no. SRRXXXXXX at the Sequence Read Archive (SRA) (http://www.ncbi.nlm.nih.gov/sra, accessed XX XX XXXX), and the associated BioProject alias is PRJNAXXXXXX. The whole genome Illumina sequences of the 60 E. paniculata individuals in this study are available under Accession no. SRRXXXXXXXX - SRRXXXXXXXX at the Sequence Read Archive (SRA) (http://www.ncbi.nlm.nih.gov/sra, accessed XX XX XXXX), and the associated BioProject alias is PRJNA862356. The RNA-Seq data of styles and stamens are available under Accession no. SRRXXXXXXXX - SRRXXXXXXXX at the Sequence Read Archive (SRA) (http://www.ncbi.nlm.nih.gov/sra, accessed XX XX XXXX), and the associated BioProject alias is PRJNAXXXXXX. The E. paniculata SsMm genome and its MAKER and InterProScan annotation are available under Genome ID XXXXX on CoGe (https://genomevolution.org/coge/, accessed XX XX XXXX). The E. paniculata ssmm genome and its MAKER and InterProScan annotation are available under Genome ID XXXXX on CoGe (https://genomevolution.org/coge/, accessed XX XX XXXX). References Darwin, C. R. The different forms of flowers on plants of the same species. (John Murray, 1877). Barrett, S. C. H. Heterostylous genetic polymorphisms: model systems for evolutionary analysis. In Evolution and function of heterostyly (ed. Barrett, S. C. H.) 1–29 (Springer, 1992). Barrett, S. C. H. The evolution of plant sexual diversity. Nat. Rev. Genet. 3 , 274-284 (2002). Lloyd, D. G. & Webb, C. The selection of heterostyly. In Evolution and function of heterostyly (ed. Barrett, S. C. H.) 179–207 (Springer, 1992). Lloyd, D. G. & Webb, C. J. The evolution of heterostyly. In Evolution and function of heterostyly (ed. Barrett, S. C. H.) 151–178 (Springer, 1992). Barrett, S. C. H., Jesson, L. K. & Baker, A. M. The evolution and function of stylar polymorphisms in flowering plants. Ann. Bot. 85 , 253-265 (2000). Lewis, D. & Jones, D. The genetics of heterostyly. In Evolution and function of heterostyly (Ed. Barrett, S. C. H.) 129–150 (Springer, 1992). Mather, K. The genetical architecture of heterostyly in Primula sinensis . Evolution 4 , 340–352 (1950). Schwander, T., Libbrecht, R. & Keller, L. Supergenes and complex phenotypes. Current Biology 24 , 288–294 (2014). Thompson, M. J. & Jiggins, C. Supergenes and their role in evolution. Heredity 113 , 1-8 (2014). Maney, D. L. & Küpper, C. Supergenes on steroids. Philos. Trans. R. Soc. B 377 , 20200507 (2022). Kappel, C., Huu, C. N. & Lenhard, M. A short story gets longer: recent insights into the molecular basis of heterostyly. J. Exp. Bot. 68 , 5719-5730 (2017). Barrett, S. C. H. ‘A most complex marriage arrangement’: recent advances on heterostyly and unresolved questions. New Phytol. 224 , 1051-1067 (2019). Huu, C. N. et al. Presence versus absence of CYP734A50 underlies the style-length dimorphism in primroses. eLife 5 , e17956 (2016). Li, J. et al. Genetic architecture and evolution of the S locus supergene in Primula vulgaris . Nat. Plants 2 , 1-7 (2016). Shore, J. S. et al. The long and short of the S -locus in Turnera (Passifloraceae). New Phytol. 224 , 1316–1329 (2019). Gutiérrez-Valencia, J. et al. Genomic analyses of the Linum distyly supergene reveal convergent evolution at the molecular level. Curr. Biol. 32 , 4360–4371 (2022). Zhao, Z. et al. Genomic evidence supports the genetic convergence of a supergene controlling the distylous floral syndrome. New Phytol. 237 , 601–614 (2023). Yang, J. et al. Haplotype-resolved genome assembly provides insights into the evolution of S -locus supergene in distylous Nymphoides indica . New Phytol. 240 , 2058-2071 (2023). Fawcett, J. A. et al. Genome sequencing reveals the genetic architecture of heterostyly and domestication history of common buckwheat. Nature Plants 9 , 1236-1251 (2023). Potente, G. et al. Comparative genomics elucidates the origin of a supergene controlling floral heteromorphism. Mol. Biol. Evol. 39 , msac035 (2022). Barrett, S. C. H. The evolutionary biology of tristyly. In Oxford Surveys in Evolutionary Biology (Eds. D. Futuyma & J. Antonovics) vol. 9 , 283–326 (Oxford University Press, 1993). Thompson, J. D., Pailler, T., Strasberg, D. & Manicacci, D. Tristyly in the endangered Mascarene Island endemic Hugonia serrata (Linaceae). Am. J. Bot. 83 , 1160–1167 (1996). Charlesworth, D. The evolution and breakdown of tristyly. Evolution 33 , 489-498. (1979) Arunkumar, R., Wang, W., Wright, S. I. & Barrett, S. C. H. The genetic architecture of tristyly and its breakdown to self-fertilization. Mol. Ecol. 26 , 752–765 (2017). Ossowski, S. et al. The rate and molecular spectrum of spontaneous mutations in Arabidopsis thaliana . Science 327 , 92–94 (2010). Eckenwalder, J. E. & Barrett, S. C. H. Phylogenetic systematics of Pontederiaceae. Syst. Bot. 11 , 373–391 (1986). Ness, R. W., Siol, M. & Barrett, S. C. H. De novo sequence assembly and characterization of the floral transcriptome in cross-and self-fertilizing plants. BMC Genomics 12 , 298 (2011). Song, B. et al. AnchorWave: Sensitive alignment of genomes with high sequence diversity, extensive structural polymorphism, and whole-genome duplication. Proc. Natl. Acad. Sci. 119 , e2113075119 (2022). Jiao, Z., Du, H., Chen, S., Huang, W. & Ge, L. LAZY gene family in plant gravitropism. Front. Plant Sci. 11 , 606241 (2021). Chen, Y. et al. Identification of the abundant hydroxyproline-rich glycoproteins in the root walls of wild-type Arabidopsis , an ext3 mutant line, and its phenotypic revertant. Plants 4 , 85–111 (2015). Sozzani, R. et al. Spatiotemporal regulation of cell-cycle genes by SHORTROOT links patterning and growth. Nature , 466 , 128-132 (2021). Gutiérrez-Valencia, J. et al. The genomic architecture and evolutionary fates of supergenes. Genome Biology and Evolution 13 , evab057 (2021). Bouckaert, R. et al. BEAST 2.5: An advanced software platform for Bayesian evolutionary analysis. PLoS Comput. Biol . 15 , e1006650 (2019). Ness, R. W., Wright, S. I. & Barrett, S. C. H. Mating-system variation, demographic history and patterns of nucleotide diversity in the tristylous plant Eichhornia paniculata . Genetics 184 , 381–392 (2010). Koren, S. et al. Canu: scalable and accurate long-read assembly via adaptive k-mer weighting and repeat separation. Genome Res . 27 , 722–736 (2017). Roach, M. J., Schmidt, S. A., & Borneman, A. R. Purge Haplotigs: allelic contig reassignment for third-gen diploid genome assemblies. BMC bioinformatics 19 , 460 (2018). Kronenberg, Z. N . et al. Extended haplotype-phasing of long-read de novo genome assemblies using Hi-C. Nature Communications 12 , 1935 (2021). Thress, K., Henzel, W., Shillinglaw, W. & Kornbluth, S. Scythe: a novel reaper-binding apoptotic regulator. EMBO J. 17 , 6135–6143 (1998). Sedlazeck, F. J., Rescheneder, P. & Von Haeseler, A. NextGenMap: fast and accurate read mapping in highly polymorphic genomes. Bioinformatics 29 , 2790–2791 (2013). Danecek, P. et al. Twelve years of SAMtools and BCFtools. GigaScience 10 , (2021). Rastas, P. Lep-MAP3: robust linkage mapping even for low-coverage whole genome sequencing data. Bioinformatics 33 , 3726–3732 (2017). Durand, N. C. et al. Juicebox provides a visualization system for Hi-C contact maps with unlimited zoom. Cell Syst. 3 , 99–101 (2016). Walker, B. J. et al. Pilon: an integrated tool for comprehensive microbial variant detection and genome assembly improvement. PLoS One 9 , e112963 (2014). Smit, A. F. A. & Hubley, R. RepeatModeler Open-1.0. 2008–2015. (2008). Smit, A. F. A., Hubley, R. & Green, P. RepeatMasker Open-4.0. (2013). Cheng, H., Concepcion, G. T., Feng, X., Zhang, H. & Li, H. Haplotype-resolved de novo assembly using phased assembly graphs with hifiasm. Nat. Methods 18 , 170–175 (2021). Putnam, N. H. et al. Chromosome-scale shotgun assembly using an in vitro method for long-range linkage. Genome Res. 26 , 342–350 (2016). Gurevich, A., Saveliev, V., Vyahhi, N. & Tesler, G. QUAST: quality assessment tool for genome assemblies. Bioinformatics 29 , 1072–1075 (2013). Manni, M. et al. BUSCO update: novel and streamlined workflows along with broader and deeper phylogenetic coverage for scoring of eukaryotic, prokaryotic, and viral genomes. Mol. Biol. Evol. 38 , 4647-4654 (2021). Levy Karin, E., Mirdita, M., & Söding, J. MetaEuk—sensitive, high-throughput gene discovery, and annotation for large-scale eukaryotic metagenomics. Microbiome 8 , 48 (2020). Sedlazeck, F. J., Rescheneder, P. & Von Haeseler, A. NextGenMap: fast and accurate read mapping in highly polymorphic genomes. Bioinformatics 29 , 2790-2791 (2013). Garrison, E. & Marth, G. Haplotype-based variant detection from short-read sequencing. ArXiv Prepr. , ArXiv12073907 (2012). Tange, O. Gnu parallel: the command-line power tool. Usenix Mag . 36 , 42 (2011). Danecek, P. et al. The variant call format and VCFtools. Bioinformatics 27 , 2156–2158 (2011). Garrison, E., Kronenberg, Z. N., Dawson, E. T., Pedersen, B. S. & Prins, P. A spectrum of free software tools for processing the VCF variant call format: vcflib, bio-vcf, cyvcf2, hts-nim and slivar. PLoS Comput. Biol. 18 , e1009123 (2022). Puritz, J. B., Hollenbeck, C. M. & Gold, J. R. dDocent: a RADseq, variant-calling pipeline designed for population genomics of non-model organisms. PeerJ 2 , e431 (2014). Zhou, X. & Stephens, M. Genome-wide efficient mixed-model analysis for association studies. Nat. Genet . 44 , 821–824 (2012). Xue, L. et al. Evidences for a role of two Y-specific genes in sex determination in Populus deltoides . Nature Communications 11 , 5893 (2020). Pedersen, B. S. & Quinlan, A. R. Mosdepth: quick coverage calculation for genomes and exomes. Bioinformatics 34 , 867–868 (2018). Kim, D., Langmead, B. & Salzberg, S. L. HISAT: a fast spliced aligner with low memory requirements. Nat. Methods 12 , 357–360 (2015). Pertea, M. et al. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat. Biotechnol. 33 , 290–295 (2015). Lohse, M . et al. Mercator: a fast and simple web server for genome scale functional annotation of plant sequence data. Plant, Cell & Environment 37 , 1250-1258 (2014). Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15 , 1–21 (2014). Arunkumar, R., Ness, R. W., Wright, S. I. & Barrett, S. C. H. The evolution of selfing is accompanied by reduced efficacy of selection and purging of deleterious mutations. Genetics 199 , 817–829 (2015). Arunkumar, R., Maddison, T. I., Barrett, S. C. H. & Wright, S. I. Recent mating-system evolution in Eichhornia is accompanied by cis-regulatory divergence. New Phytol. 211 , 697–707 (2016). Cantarel, B. L. et al. MAKER: an easy-to-use annotation pipeline designed for emerging model organism genomes. Genome Res. 18 , 188–196 (2008). Jones, P. et al. InterProScan 5: genome-scale protein function classification. Bioinformatics 30 , 1236–1240 (2014). Katoh, K. & Standley, D. M. MAFFT Multiple Sequence Alignment Software Version 7: Improvements in Performance and Usability. Mol. Biol. Evol. 30 , 772-780 (2013). Posada, D. jModelTest: Phylogenetic Model Averaging. Mol. Biol. Evol. 25 , 1253-1256 (2008). Huang, W. et al. A well-supported nuclear phylogeny of Poaceae and implications for the evolution of C4 photosynthesis. Mol. Plant 15 , 755-777 (2022). Magallón, S., Gómez-Acevedo, S., Sánchez-Reyes, L. L. & Hernández-Hernández, T. A metacalibrated time-tree documents the early rise of flowering plant phylogenetic diversity. New Phytol. 207 , 437-453 (2015). Camacho, C. et al. BLAST+: architecture and applications. BMC Bioinformatics 10 , 421 (2009). Haug-Baltzell, A., Stephens, S. A., Davey, S., Scheidegger, C. E. & Lyons, E. SynMap2 and SynMap3D: web-based whole-genome synteny browsers. Bioinformatics 33 , 2197–2198 (2017). Korneliussen, T. S., Albrechtsen, A., & Nielsen, R. ANGSD: analysis of next generation sequencing data. BMC bioinformatics 15 , 356 (2014). Additional Declarations There is NO Competing Interest. Supplementary Files GenomicevidenceforsupergenecontrolofthetristylousfloralpolymorphismSupplementaryMaterials.docx Cite Share Download PDF Status: Posted Version 1 posted You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-4232430","acceptedTermsAndConditions":true,"allowDirectSubmit":true,"archivedVersions":[],"articleType":"Article","associatedPublications":[],"authors":[{"id":288771696,"identity":"89373c53-1ca2-4b74-bddd-88ef276d6a82","order_by":0,"name":"Haoran Xue","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAABBklEQVRIiWNgGAWjYBACAwST+QCQa0GSFrYEIFcCSKNL4NbCA2ITocWc/fCxDz93MCRuON7zdcOPAgkGg/vNxz583PFHjoH98ANsWix70pJn9p5hMDY4c3bbzR6gwwyOsSXPnHnGwJiBJw2rVQYHcowZeNsY5Axu5G67wQPWwmPMzNtmkNgggd11BuffGDP+bQN65P6bZzf/oGph/4BVy40ckAKQLTxst9Fs4cFuy41nycyybRLGkmfSzG7LGEjwSB5LS2ac2WZszMaTU4DdYcmHGd+22ST2HT/87OabPzZyfIcPH2b42CYnx89+fAPWUIYACTiLB85iw6N+FIyCUTAKRgF+AAAH7Vhx6WNnXwAAAABJRU5ErkJggg==","orcid":"https://orcid.org/0000-0003-1429-0996","institution":"University of Toronto","correspondingAuthor":true,"prefix":"","firstName":"Haoran","middleName":"","lastName":"Xue","suffix":""},{"id":288771697,"identity":"ee5edf18-cca4-4913-baf0-84cc626f5e71","order_by":1,"name":"Yunchen Gong","email":"","orcid":"","institution":"University of Toronto","correspondingAuthor":false,"prefix":"","firstName":"Yunchen","middleName":"","lastName":"Gong","suffix":""},{"id":288771698,"identity":"a9b4e833-e20e-4129-9b73-5c09058f1740","order_by":2,"name":"Stephen Wright","email":"","orcid":"https://orcid.org/0000-0001-9973-9697","institution":"University of Toronto","correspondingAuthor":false,"prefix":"","firstName":"Stephen","middleName":"","lastName":"Wright","suffix":""},{"id":288771699,"identity":"a04cb406-72df-446d-8810-25f8d87bc5d5","order_by":3,"name":"Spencer Barrett","email":"","orcid":"","institution":"University of Toronto","correspondingAuthor":false,"prefix":"","firstName":"Spencer","middleName":"","lastName":"Barrett","suffix":""}],"badges":[],"createdAt":"2024-04-07 17:45:11","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-4232430/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-4232430/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":54471963,"identity":"cc136fbd-9572-4318-8b16-c9eff465cb8a","added_by":"auto","created_at":"2024-04-11 05:10:49","extension":"png","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":135916,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eTristyly in \u003c/strong\u003e\u003cem\u003e\u003cstrong\u003eEichhornia paniculata\u003c/strong\u003e\u003c/em\u003e\u003cstrong\u003e. a, \u003c/strong\u003eFlower of the S-morph\u003cem\u003e\u003cstrong\u003e \u003c/strong\u003e\u003c/em\u003eillustrating exerted long- and mid-level stamens. The short style is concealed in the perianth tube.\u0026nbsp; \u003cstrong\u003eb,\u003c/strong\u003e The arrangement\u003cstrong\u003e \u003c/strong\u003eof styles and stamens in the three floral morphs of a tristylous\u003cem\u003e \u003c/em\u003especies. L, M and S refer to the long-, mid- and short-styled morphs, respectively. The arrows indicate “legitimate” cross-pollinations between anthers and stigmas of equivalent height. In most tristylous species these are the only compatible matings because of trimorphic incompatibility. However, in \u003cem\u003eE.\u003c/em\u003e \u003cem\u003epaniculata\u003c/em\u003e “illegitimate” pollinations are permitted as the species is self- and intramorph compatible\u003csup\u003e35\u003c/sup\u003e. Genotypes for the floral morphs of \u003cem\u003eE. paniculata\u003c/em\u003e are indicated.\u003c/p\u003e","description":"","filename":"1.png","url":"https://assets-eu.researchsquare.com/files/rs-4232430/v1/b686d85050a8d44dedcf069a.png"},{"id":54471960,"identity":"22f1f561-8cf9-45c2-a21e-9de7cbcc503f","added_by":"auto","created_at":"2024-04-11 05:10:49","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":272576,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eThe recombination landscape, genome content, and nucleotide diversity (π) across the chromosomes of \u003c/strong\u003e\u003cem\u003e\u003cstrong\u003eEichhornia paniculata\u003c/strong\u003e\u003c/em\u003e\u003cstrong\u003e.\u003c/strong\u003e \u003cstrong\u003ea,\u003c/strong\u003eMarey maps plotting the position on the linkage map (cM) against the position on the physical map (Mb). \u003cstrong\u003eb-e,\u003c/strong\u003eDensities of genes (\u003cstrong\u003eb\u003c/strong\u003e), all repetitive elements (\u003cstrong\u003ec\u003c/strong\u003e), LTR retrotransposons (\u003cstrong\u003ed\u003c/strong\u003e), and non-LTR repeats (\u003cstrong\u003ee\u003c/strong\u003e), all measured as proportions of 1-Mb windows. \u003cstrong\u003eh,\u003c/strong\u003e Nucleotide diversity (π) of 1-Mb windows.\u003c/p\u003e","description":"","filename":"2.png","url":"https://assets-eu.researchsquare.com/files/rs-4232430/v1/074259853283ab738ccb04f1.png"},{"id":54471973,"identity":"0d8c1a2f-56d1-498b-961d-a2ecba995815","added_by":"auto","created_at":"2024-04-11 05:10:50","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":136910,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eSNPs associated with tristyly in \u003c/strong\u003e\u003cem\u003e\u003cstrong\u003eEichhornia paniculata\u003c/strong\u003e\u003c/em\u003e\u003cstrong\u003e and morph-specific read coverage, recombination landscape, and repeat content of the \u003c/strong\u003e\u003cem\u003e\u003cstrong\u003eS-\u003c/strong\u003e\u003c/em\u003e\u003cstrong\u003e and \u003c/strong\u003e\u003cem\u003e\u003cstrong\u003eM-\u003c/strong\u003e\u003c/em\u003e\u003cstrong\u003eloci. a, \u003c/strong\u003eSNPs associated with S- vs. non-S-morphs. The blue lines indicate \u003cem\u003ep\u003c/em\u003e = 2.52397779 x 10\u003csup\u003e-8\u003c/sup\u003e.\u003cstrong\u003e \u003c/strong\u003eInset panels zoom in on chromosome 2. \u003cstrong\u003eb,\u003c/strong\u003e SNPs associated with M- vs. L-morphs. The blue lines indicate \u003cem\u003ep\u003c/em\u003e = 2.23826411 x 10\u003csup\u003e-8\u003c/sup\u003e. Inset panels zoom in on chromosome 2. \u003cstrong\u003ec\u003c/strong\u003e, Average normalized coverages of the L- (red), M- (blue), and S- (cyan) morph plants at the \u003cem\u003eS\u003c/em\u003e-locus (\u003cem\u003eS\u003c/em\u003e- and \u003cem\u003es\u003c/em\u003e-alleles assembled next to each other)\u003cem\u003e.\u003c/em\u003e The data were based on 20 plants per morph. \u003cstrong\u003ed, \u003c/strong\u003eThe average normalized coverages of the L- (red), M- (blue), and S- (cyan) morph plants at the \u003cem\u003eM\u003c/em\u003e-locus\u003cem\u003e.\u003c/em\u003e The data were based on 20 plants per morph. \u003cstrong\u003ee,\u003c/strong\u003e The recombination landscape of chromosome 2. \u003cstrong\u003ef-h, \u003c/strong\u003eComparisons of proportion of repetitive elements between the \u003cem\u003eS\u003c/em\u003e-allele (ten 250-kb windows), \u003cem\u003es\u003c/em\u003e-allele (ten 65-kb windows), and \u003cem\u003eM\u003c/em\u003e-allele (ten 25-kb windows), and their flanking regions (five windows upstream and five windows downstream to the locus, with window size equal to the allele being compared with).\u003c/p\u003e","description":"","filename":"3.png","url":"https://assets-eu.researchsquare.com/files/rs-4232430/v1/16cf76edb9b33dee2f99da58.png"},{"id":54471962,"identity":"16949327-04f3-4a13-adc4-5c7004adbd74","added_by":"auto","created_at":"2024-04-11 05:10:49","extension":"png","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":93666,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eThe composition and proposed model of the evolution of the \u003c/strong\u003e\u003cem\u003e\u003cstrong\u003eS-\u003c/strong\u003e\u003c/em\u003e\u003cstrong\u003e and \u003c/strong\u003e\u003cem\u003e\u003cstrong\u003eM-\u003c/strong\u003e\u003c/em\u003e\u003cstrong\u003eloci in \u003c/strong\u003e\u003cem\u003e\u003cstrong\u003eEichhornia paniculata\u003c/strong\u003e\u003c/em\u003e\u003cstrong\u003e.\u003c/strong\u003e \u003cem\u003e\u003cstrong\u003ea,\u003c/strong\u003e\u003c/em\u003e Genomic organization of the \u003cem\u003eS\u003c/em\u003e-locus genes. Black dashed lines represent the boundaries of the tristyly alleles inferred from the read coverage analysis. Blue dashed lines represent inferred boundaries of an inversion.\u003cstrong\u003e b,\u003c/strong\u003e Genomic organization of the genes around LAZY1-M. Green dashed lines represent the boundaries of LAZY1-M and its paralog LAZY1-P1. \u003cstrong\u003ec,\u003c/strong\u003eHypothesized evolutionary pathway for the organization of the \u003cem\u003eS-\u003c/em\u003e and \u003cem\u003eM\u003c/em\u003e-loci.\u003c/p\u003e","description":"","filename":"4.png","url":"https://assets-eu.researchsquare.com/files/rs-4232430/v1/b78ae9d910ad6cc64a129100.png"},{"id":54471971,"identity":"541d9036-cebe-4f18-a25a-1b0761c9b593","added_by":"auto","created_at":"2024-04-11 05:10:50","extension":"png","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":72607,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eExpression patterns of the three genes that are specific to the \u003c/strong\u003e\u003cem\u003e\u003cstrong\u003eS\u003c/strong\u003e\u003c/em\u003e\u003cstrong\u003e- or \u003c/strong\u003e\u003cem\u003e\u003cstrong\u003eM\u003c/strong\u003e\u003c/em\u003e\u003cstrong\u003e-allele and expressed in the styles or stamens of \u003c/strong\u003e\u003cem\u003e\u003cstrong\u003eEichhornia paniculata\u003c/strong\u003e\u003c/em\u003e\u003cstrong\u003e. a, Genes specific to the \u003c/strong\u003e\u003cem\u003e\u003cstrong\u003eS\u003c/strong\u003e\u003c/em\u003e\u003cstrong\u003e-allele. b, LAZY1-M, which is specific to the \u003c/strong\u003e\u003cem\u003e\u003cstrong\u003eM\u003c/strong\u003e\u003c/em\u003e\u003cstrong\u003e-allele.\u003c/strong\u003e Samples of different morphs are separated by solid lines. For each morph, samples of different floral organs (styles or stamens) are separated by dashed lines. For each combination of morph and floral organ, the samples at the early and late development stages are on the left and right, respectively. Each type of sample had four biological replicates.\u003c/p\u003e","description":"","filename":"5.png","url":"https://assets-eu.researchsquare.com/files/rs-4232430/v1/f4fd3e3a78a4ebb8b79f95f8.png"},{"id":54471974,"identity":"d08d8a33-562f-4ec6-82b0-3581a8914728","added_by":"auto","created_at":"2024-04-11 05:10:50","extension":"png","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":69772,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003ea, Phylogenetic inference of the dates of divergence between LAZY1-P1 and LAZY1-S and between LAZY1 and LAZY1-M in the ancestors of \u003c/strong\u003e\u003cem\u003e\u003cstrong\u003eEichhornia paniculata\u003c/strong\u003e\u003c/em\u003e\u003cstrong\u003e. \u0026nbsp;\u003c/strong\u003eBayesian tree of LAZY1 genes in \u003cem\u003eOryza sativa\u003c/em\u003e, \u003cem\u003eZea mays\u003c/em\u003e, \u003cem\u003eMusa acuminata\u003c/em\u003e, and \u003cem\u003eE. paniculata\u003c/em\u003e. Numbers outside of parentheses represent the estimated time of divergence. Blue rectangles and numbers in parentheses are the 95% confidence intervals of the estimated time of divergence. \u003cstrong\u003eb, The postulated sequence of the evolution to trimorphism from monomorphism.\u003c/strong\u003e\u003c/p\u003e","description":"","filename":"6.png","url":"https://assets-eu.researchsquare.com/files/rs-4232430/v1/b3e44621fefd24412fe53c83.png"},{"id":73329143,"identity":"0f770685-fdde-4506-9b9a-7fc5b5483ebc","added_by":"auto","created_at":"2025-01-09 02:28:57","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":2157590,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-4232430/v1/a69c4305-21e0-4bc1-aae6-bede3c087b47.pdf"},{"id":54471964,"identity":"adc98580-b428-4919-b6bb-5ba05101b952","added_by":"auto","created_at":"2024-04-11 05:10:49","extension":"docx","order_by":1,"title":"","display":"","copyAsset":false,"role":"supplement","size":1385974,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cbr\u003e\u003c/p\u003e","description":"","filename":"GenomicevidenceforsupergenecontrolofthetristylousfloralpolymorphismSupplementaryMaterials.docx","url":"https://assets-eu.researchsquare.com/files/rs-4232430/v1/909437d3a08486d61e0c3f77.docx"}],"financialInterests":"There is \u003cb\u003eNO\u003c/b\u003e Competing Interest.","formattedTitle":"Genomic evidence for supergene control of Darwin’s “complex marriage arrangement” – the tristylous floral polymorphism","fulltext":[{"header":"Introduction","content":"\u003cp\u003eFlowering plants possess a remarkable diversity of floral strategies promoting cross-pollination and outcrossed mating. This variation attracted the attention of Charles Darwin who devoted most of his classic book \u0026ldquo;\u003cem\u003eThe Different Forms of Flowers on Plants of the Same Species\u003c/em\u003e\u0026rdquo;\u003csup\u003e1\u003c/sup\u003e to the floral polymorphism heterostyly in which populations are characterized\u0026nbsp;by two (distyly) or three (tristyly) floral morphs, which possess a reciprocal arrangement of stigma and anther heights\u003csup\u003e2\u003c/sup\u003e. \u0026nbsp;Darwin described heterostyly as a \u0026lsquo;complex marriage arrangement\u0026rsquo; promoting cross-pollination between the floral morphs by animal pollinators owing to segregated pollen placement on their bodies (see Fig. 1 in Ref. 3). Contemporary experimental evidence generally supports this hypothesis also finding that heterostyly limits sexual interference between female and male function thus reducing pollen wastage\u003csup\u003e3-4\u003c/sup\u003e. Heterostyly is reported from 28 animal-pollinated angiosperm families\u003csup\u003e5-6\u003c/sup\u003e and represents a remarkable example of adaptive convergence in form and function. Controlled crosses on over a dozen unrelated species have established that distyly is governed by a single diallelic Mendelian locus (\u003cem\u003eS\u003c/em\u003e) with dominance, whereas tristyly was controlled by two diallelic epistatically interacting loci (\u003cem\u003eS\u003c/em\u003e, \u003cem\u003eM\u003c/em\u003e), both with dominance\u003csup\u003e7\u003c/sup\u003e (Fig. 1). To resolve the discrepancy between the complexity of the heterostylous syndrome and its relatively simple inheritance, it has been proposed that heterostyly loci are in fact supergenes comprised of clusters of functionally related genes\u003csup\u003e8\u003c/sup\u003e. Evidence indicates that these genes are tightly linked to each other due to suppression of recombination and are inherited as a single unit\u003csup\u003e7\u003c/sup\u003e.\u0026nbsp;\u003c/p\u003e\n\u003cp\u003eSupergenes are also known to control several other complex polymorphic traits, including Batesian mimicry in \u003cem\u003eHeliconius\u003c/em\u003e butterflies, social organizations of fire ants, and plumage variation and behaviour in ruff\u003csup\u003e9-11\u003c/sup\u003e. Of particular importance for understanding the evolution of supergenes is to determine their genetic architecture and the mechanisms causing suppression of recombination. Identification of the causal loci can also provide important insights into the age of the alleles, to assess the evidence for long-term balanced polymorphisms maintaining variation. Recent efforts have been made to elucidate the genetic architecture of heterostylous supergenes in a handful of distylous systems (reviewed in\u003csup\u003e12-13\u003c/sup\u003e).\u0026nbsp;Significantly, genomic studies on \u003cem\u003eFagopyrum,\u003c/em\u003e \u003cem\u003eGelsemium\u003c/em\u003e, \u003cem\u003eLinum\u003c/em\u003e, \u003cem\u003eNymphoides\u003c/em\u003e, \u003cem\u003ePrimula,\u003c/em\u003e and \u003cem\u003eTurnera\u003c/em\u003e, each belonging to a different angiosperm lineage, have demonstrated that distyly is controlled by several hemizygous genes present in the short-styled morph (S-morph) but absent from the long-styled morph (L-morph)\u003csup\u003e14-20\u003c/sup\u003e. For example, in \u003cem\u003ePrimula\u003c/em\u003e, the \u003cem\u003eS\u003c/em\u003e-locus is a hemizygous region containing five genes controlling the expression of different components of the distylous polymorphism\u003csup\u003e14-15, 21\u003c/sup\u003e.\u003c/p\u003e\n\u003cp\u003eAll molecular investigations of heterostyly have focused on distylous species and the supergene hypothesis has yet to be investigated in any of the six angiosperm families in which the more complex tristyly occurs\u003csup\u003e22-23\u003c/sup\u003e. Charlesworth\u003csup\u003e24\u003c/sup\u003e suggested that alternatives to the supergene hypothesis should be considered for tristyly because of the absence of empirical evidence for supergenes. Thus, the genetic basis and evolutionary assembly of this two-locus epistatically interacting genetic system remains unclear.\u003c/p\u003e\n\u003cp\u003eHere, we generate chromosome-level \u003cem\u003ede novo\u003c/em\u003e genome assemblies with long-read sequences and integrate them with data from population genomic and transcriptomic analysis to investigate the genetic architecture of tristyly in \u003cem\u003eEichhornia paniculata\u003c/em\u003e, an annual aquatic native to the Neotropics. Earlier molecular studies demonstrated that the \u003cem\u003eS-\u003c/em\u003e and \u003cem\u003eM\u003c/em\u003e-loci in this species are closely linked (separated by a map distance of 2.7 cM), and that the \u003cem\u003eM\u003c/em\u003e-locus is located on linkage group 5 in a region containing 204 genes\u003csup\u003e25\u003c/sup\u003e, but the genomic location of the \u003cem\u003eS\u003c/em\u003e-locus was not investigated, and they did not resolve whether the \u003cem\u003eM\u003c/em\u003e-locus is a supergene or merely contains one pleiotropic gene\u003csup\u003e12\u003c/sup\u003e.\u0026nbsp;\u003c/p\u003e\n\u003cp\u003eTo investigate this issue and provide novel insight into the molecular basis of tristyly we addressed the following questions: 1) What is the genomic architecture of the \u003cem\u003eS-\u003c/em\u003e and \u003cem\u003eM-\u003c/em\u003eloci and are they indeed supergenes? 2) What are the likely candidate genes at the two loci controlling different components of tristyly? 3) What is the evolutionary timing of the origins of the \u003cem\u003eS-\u003c/em\u003e and \u003cem\u003eM-\u003c/em\u003eloci? We note that in the only theoretical analysis of the evolution of tristyly, the origin of the \u003cem\u003eS\u003c/em\u003e-locus was inferred to have preceded the origin of the \u003cem\u003eM\u003c/em\u003e-locus\u003csup\u003e24\u003c/sup\u003e. Our results are consistent with this stepwise origin of the two loci, with recurrent gene duplication from the same ancestral gene playing a key role in the formation of the epistatic interaction between the two loci. \u0026nbsp;\u003c/p\u003e"},{"header":"Results","content":"\u003cp\u003e\u003cstrong\u003eComparisons of chromosome-level genome assemblies\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eWe generated chromosome-level genome assemblies from a heterozygous S-morph individual (\u003cem\u003eSsMm\u003c/em\u003e) and an L-morph individual (\u003cem\u003essmm\u003c/em\u003e) combining long-read PacBio sequencing with proximity ligation (Omni-C and Hi-C) scaffolding (see methods and Fig. S1). In bothassemblies, most of the two genomes were assembled into eight major chromosomal scaffolds (Table 1). Comparisons of the two assemblies demonstrated large-scale concordance, with small-scale differences likely reflecting inversions and/or small assembly errors (Fig. S2a). There was clear evidence for a recent whole genome duplication event in which each chromosome showed strong synteny with one or more additional chromosomes. Based on the distribution of synonymous substitution rate (\u003cem\u003eKs\u003c/em\u003e) values between paralogs, the most recent whole genome duplication event has a mode of \u003cem\u003eKs\u003c/em\u003e=0.38 (Fig. S2b). Assuming a substitution rate of 7 x 10\u003csup\u003e-9\u003c/sup\u003e per year\u003csup\u003e26\u003c/sup\u003e and a whole genome duplication event, this would correspond to two sets of homeologous chromosomes that diverged less than roughly 30 million years ago. Since the base chromosome number in Pontederiaceae was inferred to be \u003cem\u003en\u003c/em\u003e=8\u003csup\u003e27\u003c/sup\u003e, the duplication event was likely shared across the family, with additional lineage-specific duplications due to additional rounds of polyploidy occurring in other lineages of the family\u003csup\u003e28\u003c/sup\u003e. Unless otherwise noted, our analysis below uses our \u003cem\u003eSsMm\u003c/em\u003e assembly as the reference genome.\u0026nbsp;\u003c/p\u003e\n\u003cp\u003eThe Marey maps (plots of the linkage map against the physical map) indicated that there was a large (likely pericentromeric) region of low recombination in the middle of each chromosome (Fig. 2a). In chromosomes 3 and 7, the region extended to one end of the chromosome, likely corresponding to submetacentric or acrocentric chromosomes, whereas the remaining chromosomes appeared to be more metacentric. As expected, the regions of reduced recombination were associated with a high repeat density (particularly LTR retrotransposons), low gene density, and low nucleotide diversity, whereas the more highly recombining chromosome tips were gene-dense, comprised of fewer repeats, and has higher nucleotide diversity (Fig. 2b-f, Fig. S3). \u0026nbsp;\u0026nbsp;\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eBoth tristyly loci occur in a large low-recombining region\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eUsing a genome-wide association study (GWAS) we identified 4,640 SNPs that were significantly associated with the S-morph (\u003cem\u003eS\u003c/em\u003e---) versus non-S-morphs (\u003cem\u003ess\u003c/em\u003e--, i.e. L- and M-morphs) (\u003cem\u003ep\u003c/em\u003e \u0026lt; 2.52397779 x 10\u003csup\u003e-8\u003c/sup\u003e, Fig. 3a). 96.8% of these SNPs occurred in a very large 63.1 Mb region on chromosome 2, embedded in an extensive region of reduced recombination on this chromosome (Fig. 3a). As expected from earlier genetic work\u003csup\u003e25\u003c/sup\u003e, the \u003cem\u003eM-\u003c/em\u003elocus region is partially linked to the \u003cem\u003eS\u003c/em\u003e-locus. We found 4,541 SNPs significantly associated with the M- (\u003cem\u003essM-\u003c/em\u003e) versus L- (\u003cem\u003essmm\u003c/em\u003e) morphs (\u003cem\u003ep\u003c/em\u003e \u0026lt; 2.23826411 x 10\u003csup\u003e-8\u003c/sup\u003e, Fig. 3b). 98.0% of these SNPs were in a 35.1 Mb region on chromosome 2, at the other end of the region of reduced recombination from the \u003cem\u003eS\u003c/em\u003e-locus region. Thus, both the \u003cem\u003eS-\u003c/em\u003e and \u003cem\u003eM\u003c/em\u003e-loci governing tristyly occur in a large, putatively pericentromeric region of low recombination on chromosome 2 causing extensive and significant SNP associations across many Mb of sequence.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eThe \u003cem\u003eS\u003c/em\u003e-locus is a supergene consisting of divergent alleles with hemizygous genes\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe \u003cem\u003eS-\u0026nbsp;\u003c/em\u003eand\u003cem\u003e\u0026nbsp;M-\u003c/em\u003elocus associated regions we identified were both extremely large with very low rates of recombination. To further narrow down possible candidate genes and to investigate the potential hemizygous regions associated with tristyly, we identified regions with significant differences in read coverage between the floral morphs.\u0026nbsp;\u003c/p\u003e\n\u003cp\u003eWe identified 50 consecutive 50-kb windows that were specific to the \u003cem\u003eS\u003c/em\u003e-allele (total length: 2.5 Mb; genome coordinates: Chr2:114,300,000-116,800,000), with read coverage significantly higher in 20 S-morph (\u003cem\u003eS\u003c/em\u003e---) plants than 40 non-S-morph (\u003cem\u003ess\u003c/em\u003e--) plants (criteria: the ratio between the average read coverage of non-S-morph and S-morph plants was lower than 0.4, with \u003cem\u003ep\u003c/em\u003e-value of the \u003cem\u003et-\u003c/em\u003etest \u0026lt; 9.73217066\u0026nbsp;×\u0026nbsp;10\u003csup\u003e-7\u003c/sup\u003e, Fig. 4a). Downstream adjacent to this 2.5 Mb region were 13 consecutive 50-kb windows (total length: 650 kb; genome coordinates: Chr2:116,800,000-117,450,000) that were specific to the \u003cem\u003es\u003c/em\u003e-allele and had read coverage significantly lower in 20 S-morph (\u003cem\u003eS\u003c/em\u003e---) plants than 40 non-S-morph (\u003cem\u003ess\u003c/em\u003e--) plants (criteria: the ratio between the average read coverage of non-S-morph and S-morph plants was higher than 1.6, with the \u003cem\u003ep\u003c/em\u003e-value of the \u003cem\u003et\u003c/em\u003e-test \u0026lt; 9.73217066\u0026nbsp;×\u0026nbsp;10\u003csup\u003e-7\u003c/sup\u003e, Fig. 4a).\u0026nbsp;\u003c/p\u003e\n\u003cp\u003eWe hypothesized that given the high divergence between alleles, our collapsed genome assembly may have created a chimeric assembly in this region, leaving the \u003cem\u003es\u003c/em\u003e-allele sequence adjacent to the \u003cem\u003eS\u003c/em\u003e-allele region in our assembled scaffold and contributing to the lack of a signal of high diversity when mapping to this reference genome. To investigate this hypothesis further, we compared our L-morph assembly (\u003cem\u003essmm\u003c/em\u003e) to this assembly in this region using whole genome alignments. As expected, Anchorwave\u003csup\u003e29\u003c/sup\u003e alignments revealed the presence of a large 2.57 Mb deletion in the \u003cem\u003essmm\u003c/em\u003e assembly relative to the \u003cem\u003eSsMm\u003c/em\u003e assembly (Fig. S4), whereas the adjacent region showed homology. Thus, we conclude that the 2.5 Mb region roughly corresponds to the \u003cem\u003eS\u003c/em\u003e-allele and the 650 kb region roughly corresponds to the \u003cem\u003es\u003c/em\u003e-allele. We refined the boundaries of the two alleles using the coordinates of the closest genes in their flanking regions (Fig. 4a).\u0026nbsp;\u003c/p\u003e\n\u003cp\u003eTwo genes (TIE/SPL and MPC2) were shared between the \u003cem\u003eS-\u003c/em\u003e and \u003cem\u003es\u003c/em\u003e-alleles. The order of the two genes are opposite in the two alleles. In the \u003cem\u003eS\u003c/em\u003e-allele, both genes are on the negative strand, whereas in the \u003cem\u003es\u003c/em\u003e–allele they are on the positive strand. This finding suggests the existence of an inversion containing the two genes (Fig. 4a). The rest of the genes were allele specific: three genes (CYCD6;1-S, HRGP-S, and LAZY1-S) were specific to the \u003cem\u003eS\u003c/em\u003e-allele, while five genes were specific to the \u003cem\u003es\u003c/em\u003e-allele (Fig. 4a). LAZY1 is a gene family that regulates the distribution of auxin\u003csup\u003e30\u003c/sup\u003e, HRGP is a superfamily of genes that forms covalently cross-linked networks in primary cell walls\u003csup\u003e31\u003c/sup\u003e, and CYCD6;1 is a gene that regulates the formative division of cells\u003csup\u003e32\u003c/sup\u003e.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eThe \u003cem\u003eM\u003c/em\u003e-locus contains one hemizygous gene\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eWe identified five consecutive 50-kb windows (total length: 250 kb; genome coordinates: Chr2:13,000,000-13,250,000) that had read coverage significantly higher in the 20 M-morph (\u003cem\u003essM-\u003c/em\u003e) plants than in the 20 L-morph (\u003cem\u003essmm\u003c/em\u003e) plants (criteria: the ratio between the average read coverage of M-morph and L-morph plants was higher than 1.6, with the \u003cem\u003ep\u003c/em\u003e-value of the \u003cem\u003et\u003c/em\u003e-test \u0026lt; 9.73217066\u0026nbsp;×\u0026nbsp;10\u003csup\u003e-7\u003c/sup\u003e, Fig. 4b). There were nine coding genes in these five windows. However, among them only one gene was missing from the \u003cem\u003essmm\u003c/em\u003e genome assembly suggesting that it is absent from the \u003cem\u003em\u003c/em\u003e-allele (Fig. 4b). Strikingly, this absent gene, LAZY1-M, also shows strong similarity to genes in the LAZY1 family, raising the possibility that this locus acts competitively in the same pathway as LAZY1-S. Consistent with the SNP association mapping, the hemizygous \u003cem\u003eM\u003c/em\u003e-locus region was located on the other side of the low-recombing region on chromosome 2 (Fig. 3e).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eGenetic degeneration at the \u003cem\u003eS\u003c/em\u003e-allele\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eIn species with supergenes, signatures of genetic degeneration, including gene loss and accumulation of TEs, are expected for the dominant haplotypes of the supergenes, because dominant haplotypes are only present in heterozygotes and thus do not exhibit recombination\u003csup\u003e33\u003c/sup\u003e. The \u003cem\u003eS\u003c/em\u003e-allele specific gene loss (the five \u003cem\u003es\u003c/em\u003e-allele specific genes) was in accord with this prediction. In addition, we found that the density of repetitive elements in the \u003cem\u003eS\u003c/em\u003e- and \u003cem\u003eM\u003c/em\u003e-alleles was significantly higher than the flanking regions (\u003cem\u003ep\u003c/em\u003e = 0.00001985 and \u003cem\u003ep\u003c/em\u003e = 0.04731, respectively, Wilcoxon signed-rank test) (Fig. 4f and h), whereas the differences between the \u003cem\u003es\u003c/em\u003e-allele and the flanking regions were not significantly different (\u003cem\u003ep\u003c/em\u003e = 0.847, Wilcoxon signed-rank test) (Fig. 4g).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eThe expression pattern of hemizygous genes\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eAmong the three \u003cem\u003eS\u003c/em\u003e-allele specific genes, LAZY1-S and HRGP-S were expressed in floral organs of the S-morph but not non-S-morph plants. LAZY1-S was expressed in styles but not stamens, whereas in contrast HRGP-S was expressed in stamens but not styles (Fig. 5a). CYCD6;1-S has low expression levels in both the styles and stamens of the S-morph (Fig. 5a). The \u003cem\u003eM\u003c/em\u003e-allele specific gene LAZY1-M was expressed in both styles and stamens (Fig. 5b).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eGenes that were differentially expressed in different morphs\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eFor genes that were differentially expressed genome-wide between floral morphs (see Methods and Table S1 for details of comparisons), three sets of genes showed functional enrichment:\u003cstrong\u003e\u0026nbsp;(\u003c/strong\u003e1) genes that were\u0026nbsp;upregulated in S-morph styles (compared with non-S-morph styles) in early buds (enriched in three non-overlapping functions), (2)\u0026nbsp;genes that were\u0026nbsp;upregulated in S-morph styles (compared with non-S-morph styles) in late buds (enriched in one function), and (3)\u0026nbsp;genes that were\u0026nbsp;upregulated in M-morph styles (compared with L-morph styles) in late buds (enriched in one function)\u0026nbsp;(Table S2). However, none of the five functions inferred had a direct connection to growth or flower organs, suggesting that only a small number of genes downstream to the \u003cem\u003eS\u003c/em\u003e- and \u003cem\u003eM\u003c/em\u003e-locus candidate genes are involved in the phenotypic differences between the floral morphs. Other differentially expressed genes may simply be incidental associations.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eThe origins of the \u003cem\u003eS\u003c/em\u003e and \u003cem\u003eM\u003c/em\u003e loci\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eTo further understand the origins of the \u003cem\u003eS-\u003c/em\u003e and \u003cem\u003eM-\u003c/em\u003eloci, we identified the closest paralogs of the \u003cem\u003eS\u003c/em\u003e- and \u003cem\u003eM\u003c/em\u003e-specific genes. We identified one gene that was paralogous to LAZY1-S and LAZY1-M. Hereafter, we refer to it as “LAZY1-P1”. LAZY1-P1 was located 223 kb downstream to LAZY1-M (Fig. 4b). The closest paralog of HRGP-S was found on the same chromosome but not in the same region (Fig. 4c). Two paralogs of CYCD6;1-S were annotated in both genome assemblies, and they are on chromosomes 3 and 5, respectively (hereafter CYCD6;1-Chr3 and CYCD6;1-Chr5). CYCD6;1-Chr3 shared higher sequence similarity with CYCD6;1-S (Fig. S5) and was located in a homologous region between chromosomes 2 and 3 (Fig. S2a), which was probably a result of the whole genome duplication estimated to have happened \u0026lt; 30 million years ago, and CYCD6;1-Chr5 is likely the ancestral copy of both CYCD6;1-S and CYCD6;1-Chr3. The distance between the three ancestral genes therefore suggests that the \u003cem\u003eS\u003c/em\u003e-locus was formed via stepwise duplication rather than segmental duplication (Fig. 4c).\u003c/p\u003e\n\u003cp\u003e\u0026nbsp; \u0026nbsp; \u0026nbsp; \u0026nbsp; \u0026nbsp; \u0026nbsp;\u0026nbsp;Based on our gene trees and assuming that the hemizygous loci represent the derived gene copies, LAZY1-S duplicated from LAZY1-P1 first, followed by subsequent gene duplication from LAZY1-P1 to form LAZY1-M. Phylogenetic dating with BEAST v2.7.6\u003csup\u003e34\u003c/sup\u003e suggested that the divergence between LAZY1-S and LAZY-P1 occurred approximately 52.4 million years ago (95% confidence interval: 44.9-59.7) and that the subsequent divergence between LAZY1-M and LAZY-P1 was around 43.0 million years ago (95% confidence interval: 36.4-50.0) (Fig. 6a).\u003c/p\u003e"},{"header":"Discussion","content":"\u003cp\u003eWe show that in tristylous \u003cem\u003eEichhornia paniculata\u003c/em\u003e, both the \u003cem\u003eS-\u003c/em\u003e and \u003cem\u003eM\u003c/em\u003e-loci are located in a low recombining region of Chromosome 2. The \u003cem\u003eS\u003c/em\u003e-locus is comprised of two divergent alleles: the 2.51 Mb \u003cem\u003eS\u003c/em\u003e-allele with three \u003cem\u003eS\u003c/em\u003e-allele specific genes that are hemizygous in most S-morph plants, and the 596 kb \u003cem\u003es\u003c/em\u003e-allele with five \u003cem\u003es\u003c/em\u003e-allele specific genes. In the \u003cem\u003eS\u003c/em\u003e-allele we identified two strong tristyly candidate genes: LAZY1-S and HRGP-S. LAZY1-S was expressed in styles but not stamens suggesting that it may inhibit style elongation in the S-morph by decreasing auxin concentration in styles. HRGP-S was expressed in stamens but not styles suggesting that it may be involved in filament elongation or possibly pollen function. Our evidence for the \u003cem\u003eS\u003c/em\u003e-locus therefore strongly supports its status as a supergene, and its origin in a large region of low recombination may have been important during the initial origins of these genes via gene duplication from more physically distant regions (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003e).\u003c/p\u003e \u003cp\u003eWe found that the \u003cem\u003eM\u003c/em\u003e-locus contained only one hemizygous gene \u0026ndash; LAZY1-M \u0026ndash; and this gene was hemizygous in --\u003cem\u003eMm\u003c/em\u003e plants. Strikingly, this gene is in the same gene family as LAZY1-S suggesting the possible role of gene duplication and neofunctionalization in the origin of the epistatic interaction between the two loci. Unlike LAZY1-S, which was only expressed in styles, LAZY1-M was expressed in both styles and stamens and with a higher expression level in stamens. This suggests that while LAZY1-S may only be involved in the shortening of styles in the S-morph, LAZY1-M may also regulate auxin distribution and organ growth in both the styles and stamens of the M-morph. Given that LAZY-M is the only \u003cem\u003eM\u003c/em\u003e-allele-specific hemizygous gene, the M-locus may not represent a supergene at this stage in its evolutionary history but simply contains LAZY1-M, a pleiotropic modifier.\u003c/p\u003e \u003cp\u003eLAZY1-M was also expressed in the stamens but not the styles of the S-morph (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eb), suggesting that LAZY1-S and LAZY1-M may compete for the same transcriptional activator and that the style expression of LAZY1-M is shut down if LAZY1-S is expressed. The expression level of LAZY1-S in S-morph styles is also much higher than LAZY1-M in M-morph styles. For example, at early development stages the average normalized count of LAZY1-M in samples of M-morph styles is 0.02% of that of LAZY1-S in S-morph styles. At late development stages, the ratio is 1.9%, indicating that the promoter region of LAZY1-S may have stronger affinity to the transcriptional activator that they compete for, potentially driving differences in style length in the S- and M-morphs through different degrees of auxin suppression.\u003c/p\u003e \u003cp\u003eIn a theoretical treatment of the evolutionary assembly of tristyly\u003csup\u003e24\u003c/sup\u003e, Charlesworth predicted that the M-morph originated last spreading in dimorphic populations of the L- and S-morph and eventually giving rise to the tristylous polymorphism. Our results on the relative ages of LAZY1-S and LAZY1-M are consistent with the hypothesis that the \u003cem\u003eM\u003c/em\u003e-locus evolved after the \u003cem\u003eS\u003c/em\u003e-locus providing empirical support for this prediction. Our results also provide evidence suggesting a major role for \u0026lsquo;re-use\u0026rsquo; through recurrent duplication of the same gene (LAZY1). Recruitment of this gene explains the origin of the \u003cem\u003eM\u003c/em\u003e-locus and birth of the third floral morph in the tristylous system (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eb) resulting in the morphological and mating complexity that characterizes Darwin\u0026rsquo;s \u0026lsquo;complex marriage arrangement\u0026rsquo;.\u003c/p\u003e"},{"header":"Online Methods","content":"\u003cp\u003e \u003cb\u003eStudy system\u003c/b\u003e. Plants used in our study originated from seed collected in wild populations of \u003cem\u003eE. paniculata\u003c/em\u003e in N.E. Brazil, Cuba and Mexico. The localities of plant material and details of the five populations used in our study are presented in Table S4. Information on the reproductive ecology, mating systems and population genetics of \u003cem\u003eE. paniculata\u003c/em\u003e are available in Ref. 35 and references cited therein. We germinated seeds using PRO-MIX Premium Potting Mix (Premier Tech, Rivi\u0026egrave;re-du-Loup, QC, Canada) and seedlings were transferred to individual 7 cm plastic pots filled with the same soil mix and placed in flooded trays in a glasshouse at the University of Toronto (Toronto, Canada) with temperatures of 25\u0026ndash;28\u0026deg;C (day) and 18\u0026ndash;21\u0026deg;C (night) and a 16-hour day length. Plants were then grown and flowered regularly for several years facilitating crosses and collection of leaves and floral buds for molecular analyses.\u003c/p\u003e \u003cp\u003e \u003cb\u003eGenome assemblies.\u003c/b\u003e We conducted two chromosome-level genome assemblies of \u003cem\u003eE. paniculata.\u003c/em\u003e The first assembly was obtained from an individual of genotype \u003cem\u003eSsMm\u003c/em\u003e (inferred from selfed progeny, Table S3) obtained from an inter-population cross between an S-morph plant (\u003cem\u003eSsM-\u003c/em\u003e) from N.E. Brazil and an M-morph plant (\u003cem\u003essMm\u003c/em\u003e) from Cuba (see Table S4 for information on the two parents). From this individual, we obtained 79 Gb (105X coverage) PacBio CLR sequences (sequenced by Arizona Genomics Institute, Tucson, AZ, USA), 55 Gb (45X coverage) Hi-C (sequenced by Phase Genomics, Seattle, WA, USA) sequences, and 117 Gb (156X coverage) 10X Chromium sequences (sequenced by The Hospital for Sick Children, Toronto, ON, Canada). The PacBio sequences were assembled with Canu\u003csup\u003e36\u003c/sup\u003e and the assembly was curated with Purge Haplotigs\u003csup\u003e37\u003c/sup\u003e, both by the Canadian Center for Computational Genomics (C3G) (Montreal, QC, Canada). Then the genome was phased with Falcon-Phase\u003csup\u003e38\u003c/sup\u003e, then scaffolded into chromosome-level scaffolds with the Hi-C data, both by Phase Genomics (Seattle, WA, USA).\u003c/p\u003e \u003cp\u003eTo identify misjoins and characterize the relationship between the physical distance and the genetic distance on each chromosome, we constructed a linkage map with genotyping-by-sequencing data obtained from a backcross pedigree from an L-morph plant and an M-morph plant, detailed in\u003csup\u003e25\u003c/sup\u003e, accessible under Accession no. SRP069136 at the Sequence Read Archive (SRA), (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttp://www.ncbi.nlm.nih.gov/sra\u003c/span\u003e\u003cspan address=\"http://www.ncbi.nlm.nih.gov/sra\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e, accessed 17 January 2023), and the associated BioProject alias PRJNA310331). The reads were trimmed with Scythe\u003csup\u003e39\u003c/sup\u003e and Seqtk (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://github.com/lh3/seqtk\u003c/span\u003e\u003cspan address=\"https://github.com/lh3/seqtk\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e). We aligned the trimmed reads to the \u003cem\u003eSsMm\u003c/em\u003e genome assembly with NextGenMap V0.5.5\u003csup\u003e40\u003c/sup\u003e, and the BAM files were processed with Picard (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttp://broadinstitute.github.io/picard\u003c/span\u003e\u003cspan address=\"http://broadinstitute.github.io/picard\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e) using the four functions FixMateInformation, CleanSam, SortSam, and AddOrReplaceReadGroups. SAMtools\u003csup\u003e41\u003c/sup\u003e was used to exclude the secondary alignments. We constructed the linkage map with the alignment files and the software program Lep-Map3\u003csup\u003e42\u003c/sup\u003e. Markers were split into eight large linkage groups using a logarithm of the odds (LOD) score limit of 70 and a Theta value of 0.25 for the initial split (SeparateChromosomes2). Most of the remaining single markers were put into these groups using a LOD score limit of 30 and a LOD difference of 5 (JoinSingles2All). We calculated marker order for each linkage group with OrderMarkers2 (with default settings). The first ten linkage groups contained 2,124 markers. We then used Juicebox\u003csup\u003e43\u003c/sup\u003e to correct the misjoins with the linkage map. The genome was polished with the 10X Chromium sequences, using the software Pilon\u003csup\u003e44\u003c/sup\u003e. The repetitive elements were annotated with RepeatModeller V1.0.11\u003csup\u003e45\u003c/sup\u003e and RepeatMasker V4.0.7\u003csup\u003e46\u003c/sup\u003e.\u003c/p\u003e \u003cp\u003eThe second assembly of \u003cem\u003eE. paniculata\u003c/em\u003e used an L-morph (\u003cem\u003essmm\u003c/em\u003e) plant from Mexico (see Table S4) using sequencing and assembly services from Dovetail Genomics (Scotts Valley, CA, USA). 59.6 Gb (~\u0026thinsp;50X coverage) PacBio HiFi reads were generated and the reads were assembled with hifiasm\u003csup\u003e47\u003c/sup\u003e. The assembly was then scaffolded with Omni-C data using the HiRise software\u003csup\u003e48\u003c/sup\u003e. We used the linkage map previously generated to concatenate scaffolds that belong to the same linkage group. We also annotated TEs with RepeatModeller V1.0.11\u003csup\u003e45\u003c/sup\u003e and RepeatMasker V4.0.7\u003csup\u003e46\u003c/sup\u003e. For both of our genome assemblies, the distribution of scaffold length was assessed with QUAST V5.2.0\u003csup\u003e49\u003c/sup\u003e, and the completeness of assemblies was assessed with BUSCO (Benchmarking Universal Single-Copy Orthologs) V5.4.4\u003csup\u003e50\u003c/sup\u003e with MetaEuk\u003csup\u003e51\u003c/sup\u003e and the embryophyta_odb10 (2024-01-08) dataset (Table S5).\u003c/p\u003e \u003cp\u003e \u003cb\u003eGenome-wide associations of single nucleotide polymorphisms (SNPs) and read coverage (RCs).\u003c/b\u003e We extracted genomic DNA from 60 \u003cem\u003eE. paniculata\u003c/em\u003e plants (20 of each floral morph) from a population in N.E. Brazil (hereafter \u0026ldquo;Pop 1\u0026rdquo;, see Table S4), using the QIAGEN DNeasy Plant Mini Kit (QIAGEN, Germantown, MD, USA). Library preparation and sequencing were conducted by the McGill University and G\u0026eacute;nome Qu\u0026eacute;bec Innovation Centre (Montr\u0026eacute;al, QC, Canada) using the 150-bp paired-end protocol on Illumina HiSeq 4000. The samples had an average sequencing depth of ~\u0026thinsp;8X.\u003c/p\u003e \u003cp\u003eThe Illumina sequencing reads were trimmed with Scythe\u003csup\u003e39\u003c/sup\u003e and Seqtk (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://github.com/lh3/seqtk\u003c/span\u003e\u003cspan address=\"https://github.com/lh3/seqtk\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e). We aligned the trimmed reads to the \u003cem\u003eSsMm\u003c/em\u003e genome assembly with NextGenMap V0.5.5\u003csup\u003e52\u003c/sup\u003e and the BAM files were processed with Picard (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttp://broadinstitute.github.io/picard/\u003c/span\u003e\u003cspan address=\"http://broadinstitute.github.io/picard/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e) using the four functions FixMateInformation, CleanSam, SortSam, and AddOrReplaceReadGroups. We used SAMtools\u003csup\u003e41\u003c/sup\u003e to exclude the secondary alignments. We used freebayes-parallel\u003csup\u003e53,54\u003c/sup\u003e to call SNPs jointly for all samples. We then filtered the VCF file with VCFtools\u003csup\u003e55\u003c/sup\u003e and vcflib2\u003csup\u003e56\u003c/sup\u003e following a published guideline\u003csup\u003e57\u003c/sup\u003e. We removed sites if more than 50% of the individuals had missing data, had a minor allele count of \u0026lt;\u0026thinsp;6, or had a quality score lower than 30. Genotypes with fewer than 4 reads were treated as \"missing\" and individuals with more than 50% missing data were not considered further. Then, sites with strong within-individual allelic bias at heterozygous sites (0.01\u0026thinsp;\u0026lt;\u0026thinsp;AB\u0026thinsp;\u0026lt;\u0026thinsp;0.25 or AB\u0026thinsp;\u0026gt;\u0026thinsp;0.75), a discrepancy in the properly paired status for reads supporting reference or alternate alleles, a low quality relative to depth, a mean depth that was \u0026gt;\u0026thinsp;13, a quality score\u0026thinsp;\u0026lt;\u0026thinsp;30, or that occurred within 5 bp of an indel, were excluded from our analysis. The resulting VCF file had 7,358,165 filtered variant sites.\u003c/p\u003e \u003cp\u003eWe conducted a Genome-Wide Association Study (GWAS) with the variant sites using GEMMA\u003csup\u003e58\u003c/sup\u003e, where a relatedness matrix and a linear mixed model were used and \u003cem\u003ep\u003c/em\u003e-values of variant sites were calculated with the Wald test.\u003c/p\u003e \u003cp\u003eTo detect hemizygous regions associated with the \u003cem\u003eS-\u003c/em\u003e and \u003cem\u003eM\u003c/em\u003e-loci, we conducted another GWAS with read-coverage, following the pipeline in Ref. 59. We split the reference genome into 25,688 windows, each of which was 50-kb, and calculated the average depth of each window in individuals using the software mosdepth\u003csup\u003e60\u003c/sup\u003e. For each individual, we normalized the coverage by dividing the coverage of each window with the average coverage across all windows, then multiplied it by 2. Therefore, the expected normalized coverage is 1 for hemizygous regions and 2 for the rest of the genome. We conducted two \u003cem\u003et\u003c/em\u003e-tests for normalized coverage for each window: one for the comparison between S-morph versus non-S-morph plants, and the other for the comparison between the M-morph versus L-morph plants. We identified tristyly allele-specific windows with the following criteria:\u003c/p\u003e \u003cp\u003e1) \u003cem\u003ep\u003c/em\u003e value of the \u003cem\u003et\u003c/em\u003e-test was smaller than 9.73217066 \u0026times; 10\u003csup\u003e\u0026minus;\u0026thinsp;7\u003c/sup\u003e, i.e. Bonferroni correction with α\u0026thinsp;=\u0026thinsp;0.05.\u003c/p\u003e \u003cp\u003e2) We set thresholds on the ratio between the average read coverage of plants of different morphs:\u003c/p\u003e \u003cp\u003eFor \u003cem\u003eS\u003c/em\u003e-allele specific windows: The read coverage of a \u003cem\u003eS\u003c/em\u003e-allele is at least 4 times of a \u003cem\u003es\u003c/em\u003e-allele, i.e. the ratio between the average read coverage of non-S-morph plants and S-morph plants was lower than 0.4: \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\frac{2RC\\left(s\\right)}{RC\\left(S\\right)+RC\\left(s\\right)}=\\frac{2\\times 1}{4+1}=0.4\\)\u003c/span\u003e\u003c/span\u003e.\u003c/p\u003e \u003cp\u003eFor \u003cem\u003es\u003c/em\u003e-allele specific windows: The read coverage of a \u003cem\u003es\u003c/em\u003e-allele is at least 4 times of a \u003cem\u003eS\u003c/em\u003e-allele, i.e. the ratio between the average read coverage of non-S-morph plants and S-morph plants was higher than 1.6: \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\frac{2RC\\left(s\\right)}{RC\\left(S\\right)+RC\\left(s\\right)}=\\frac{2\\times 4}{4+1}=1.6\\)\u003c/span\u003e\u003c/span\u003e\u003c/p\u003e \u003cp\u003eFor \u003cem\u003eM\u003c/em\u003e-allele specific windows: The read coverage of an \u003cem\u003eM\u003c/em\u003e-allele is at least 4 times of an \u003cem\u003em\u003c/em\u003e-allele, i.e. the ratio between the average read coverage of non-M-morph plants and M-morph plants was lower than 0.25, assuming all M-morph plants were homozygous: \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\frac{2RC\\left(m\\right)}{RC\\left(M\\right)+RC\\left(m\\right)}=\\frac{2\\times 1}{4+1}=0.4\\)\u003c/span\u003e\u003c/span\u003e.\u003c/p\u003e \u003cp\u003eFor \u003cem\u003em\u003c/em\u003e-allele specific windows: The read coverage of an \u003cem\u003es\u003c/em\u003e-allele is at least 4 times of a \u003cem\u003eM\u003c/em\u003e-allele, i.e. the ratio between the average read coverage of non-M-morph plants and M-morph plants was higher than 1.6: \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\frac{2RC\\left(m\\right)}{RC\\left(M\\right)+RC\\left(m\\right)}=\\frac{2\\times 4}{4+1}=1.6\\)\u003c/span\u003e\u003c/span\u003e\u003c/p\u003e \u003cp\u003eOnly regions with at least two consecutive windows that were allele-specific were considered. As a result, a 2.5 Mb \u003cem\u003eS\u003c/em\u003e-allele, a 650 kb \u003cem\u003es\u003c/em\u003e-allele, and a 250 kb \u003cem\u003eM\u003c/em\u003e-allele were identified. We split each allele into 10 equal-sized windows (window sizes: \u003cem\u003eS\u003c/em\u003e-allele: 250 kb; \u003cem\u003es\u003c/em\u003e-allele: 65 kb; \u003cem\u003eM\u003c/em\u003e-allele: 25 kb), calculated the repeat content of each window, then compared it with the 10 windows of the same size flanking each allele (5 windows upstream and 5 windows downstream of the locus) with Wilcoxon rank-sum tests.\u003c/p\u003e \u003cp\u003e \u003cb\u003eDifferential gene expression analysis and genome annotation with transcriptome data.\u003c/b\u003e To infer genes controlling different components of the tristyly syndrome, we conducted a floral organ-specific transcriptomic analysis. We dissected floral buds of the three morphs at two developmental stages: \u0026ldquo;early\u0026rdquo;, with buds between 3\u0026ndash;5 mm in length, and \u0026ldquo;late\u0026rdquo;, with buds between 5\u0026ndash;8 mm. We used sterilized forceps under a stereo dissecting microscope to open the floral buds and placed the styles and stamens into 2 mL centrifuge tubes which were then frozen by putting them into liquid nitrogen. We extracted RNA from the frozen tissues with Spectrum\u0026trade; Plant Total RNA-Kit (Sigma-Aldrich, St. Louis, MO, USA) following the manufacturer\u0026rsquo;s instructions. For each type of sample, we collected four biological replicates, each containing styles or stamens from 5\u0026ndash;27 buds sampled from 5\u0026ndash;17 individuals from each of two populations (Pop. 1, 2) originating from N.E. Brazil (Table S4, S6).\u003c/p\u003e \u003cp\u003eWe prepared sequencing libraries and sequenced on the Illumina HiSeq4000 platform (100-bp paired-end reads) by the McGill University and G\u0026eacute;nome Qu\u0026eacute;bec Innovation Centre (Montr\u0026eacute;al, QC, Canada). Then we aligned the Illumina sequences to the \u003cem\u003eSsMm\u003c/em\u003e genome assembly with HISAT2\u003csup\u003e61\u003c/sup\u003e and assembled the transcripts with StringTie\u003csup\u003e62\u003c/sup\u003e. We conducted functional annotations of the floral organ transcriptome with Mercator4\u003csup\u003e63\u003c/sup\u003e. We used StringTie\u003csup\u003e62\u003c/sup\u003e to estimate the abundance of each transcript in each sample. Then we used the R package DESeq2\u003csup\u003e64\u003c/sup\u003e to conduct differential gene expression analysis and plot the normalized count in each sample for the genes at the \u003cem\u003eS-\u003c/em\u003e and \u003cem\u003eM-\u003c/em\u003eloci. We conducted eight differential gene expression analyses for each type of morph difference (\u003cem\u003eS\u003c/em\u003e vs. non-S and M vs. L), each tissue (style vs. stamen), and each development stage (early bud vs. late bud). Only genes with FPKM (Fragments Per Kilobase of transcript per Million mapped reads)\u0026thinsp;\u0026gt;\u0026thinsp;1 and TPM (Transcript Per Million)\u0026thinsp;\u0026gt;\u0026thinsp;1 in at least one sample were included in the differential gene expression analyses. In each analysis, we identified the genes that are upregulated in each morph with false discovery rate\u0026thinsp;\u0026lt;\u0026thinsp;0.05 and LFC (log fold change)\u0026thinsp;\u0026gt;\u0026thinsp;2 or \u0026lt; -2. We then did a Mercator BIN enrichment analysis\u003csup\u003e61\u003c/sup\u003e with the list of upregulated genes for each comparison.\u003c/p\u003e \u003cp\u003eWe also annotated the repeat-masked versions of the two genome assemblies with the floral organ transcriptomes and also leaf and floral bud transcriptomes from Ref. 65, available under accession number SRP049636 on the Sequence Read Archive (SRA) (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttp://www.ncbi.nlm.nih.gov/sra\u003c/span\u003e\u003cspan address=\"http://www.ncbi.nlm.nih.gov/sra\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e, accessed 20 January 2023) and the associated BioProject alias is PRJNA266681 (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttp://www.ncbi.nlm.nih.gov/bioproject/\u003c/span\u003e\u003cspan address=\"http://www.ncbi.nlm.nih.gov/bioproject/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e, accessed 20 January 2023), and also Ref. 66, available under accession number SRP060405 at the Sequence Read Archive and associated BioProject alias (PRJNA288861) using MAKER v3.01.03\u003csup\u003e67\u003c/sup\u003e and InterProScan 5\u003csup\u003e68\u003c/sup\u003e. This resulted in a total of 35,224 annotated genes for the \u003cem\u003eSsMm\u003c/em\u003e genome and 29,592 genes for the \u003cem\u003essmm\u003c/em\u003e genome.\u003c/p\u003e \u003cp\u003e \u003cb\u003eOrigin of the\u003c/b\u003e \u003cb\u003eS\u003c/b\u003e \u003cb\u003eand\u003c/b\u003e \u003cb\u003eM\u003c/b\u003e \u003cb\u003eloci.\u003c/b\u003e According to the annotation of the \u003cem\u003eSsMm\u003c/em\u003e and \u003cem\u003essmm\u003c/em\u003e genome assemblies, there is only one LAZY1 homolog (LAZY1-P1) that exists in both genome assemblies. We then downloaded the coding sequences of LAZY1 in \u003cem\u003eOryza sativa\u003c/em\u003e and \u003cem\u003eZea mays\u003c/em\u003e from NCBI (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://www.ncbi.nlm.nih.gov/\u003c/span\u003e\u003cspan address=\"https://www.ncbi.nlm.nih.gov/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e, accession numbers: NM_001138862.1 and DQ855268.1), and two paralogs of LAZY1 in \u003cem\u003eMusa acuminata\u003c/em\u003e 'DH-Pahang' from the website Banana Genome Hub (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://banana-genome-hub.southgreen.fr/\u003c/span\u003e\u003cspan address=\"https://banana-genome-hub.southgreen.fr/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e). We conducted a multiple sequence alignment of these sequences and the coding sequences of three LAZY1 genes in \u003cem\u003eE. paniculata\u003c/em\u003e, using MAFFT v7\u003csup\u003e69\u003c/sup\u003e with L-INS-i method.\u003c/p\u003e \u003cp\u003eWe then estimated the time of divergence between LAZY1-S and LAZY1-P1 and between LAZY1-M and LAZY1-P1 with BEAST v2.7.6\u003csup\u003e34\u003c/sup\u003e. First, we did model selection of jModelTest\u003csup\u003e70\u003c/sup\u003e using Akaike information criterion (AIC) as the criteria. The GTR\u0026thinsp;+\u0026thinsp;I\u0026thinsp;+\u0026thinsp;Gamma model with nCat\u0026thinsp;=\u0026thinsp;4 was selected. Then we set the time to the most recent ancestors of \u003cem\u003eOryza sativa\u003c/em\u003e and \u003cem\u003eZea mays\u003c/em\u003e (81.43\u0026thinsp;\u0026plusmn;\u0026thinsp;0.3 MYA)\u003csup\u003e71\u003c/sup\u003e, Zingiberales and Commelinales (79.8\u0026thinsp;\u0026plusmn;\u0026thinsp;9.87 MYA)\u003csup\u003e72\u003c/sup\u003e, and Poales and Commelinales (106.7\u0026thinsp;\u0026plusmn;\u0026thinsp;8.3 MYA)\u003csup\u003e72\u003c/sup\u003e as calibration points. Ten independent Markov Chain Monte Carlo runs, each with 1 x 10\u003csup\u003e8\u003c/sup\u003e generations (first 10% was burn-in) were conducted. We combined the output of the eight runs using LogCombiner v2.7.6\u003csup\u003e34\u003c/sup\u003e with a sample frequency of 5,000. Then we used TreeAnnotator v2.7.6\u003csup\u003e34\u003c/sup\u003e to generate the maximum clade credibility tree, which was visualized in FigTree v1.4.4 (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttp://tree.bio.ed.ac.uk/software/figtree/\u003c/span\u003e\u003cspan address=\"http://tree.bio.ed.ac.uk/software/figtree/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e) (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003e).\u003c/p\u003e \u003cp\u003eThere are two CYCD6;1 paralogs that exists in both genome assemblies, and they are on chromosomes 3 and 5, respectively (hereafter CYCD6;1-Chr3 and CYCD6;1-Chr5). We then downloaded the coding sequence of CYCD6;1 in \u003cem\u003eO. sativa\u003c/em\u003e from NCBI (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://www.ncbi.nlm.nih.gov/\u003c/span\u003e\u003cspan address=\"https://www.ncbi.nlm.nih.gov/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e, accession numbers: XM_015789871.2). We conducted a multiple sequence alignment of the coding sequences of the three CYCD6;1 genes using MAFFT v7\u003csup\u003e69\u003c/sup\u003e with the L-INS-i method. Then we constructed a Neighbor Joining tree with all the conserved sites, Jukes-Cantor substitution model, and bootstrap resampling with 100 samples. (Fig. S).\u003c/p\u003e \u003cp\u003eHRGP-S was not annotated by MAKER or InterProScan. We aligned the transcript of HRGP-S against the \u003cem\u003eSsMm\u003c/em\u003e and \u003cem\u003essmm\u003c/em\u003e genome assemblies with the software program blastn\u003csup\u003e73\u003c/sup\u003e, to identify its closest paralog.\u003c/p\u003e \u003cp\u003e \u003cb\u003eSynteny analyses, whole genome alignments and structural variants\u003c/b\u003e. We used COGE\u0026rsquo;s SynMap2 algorithm\u003csup\u003e74\u003c/sup\u003e with the MAKER and InterProScan annotations to examine the synteny between the two genome assemblies and identify structural variation. To obtain finer-scale whole genome alignment of syntenic regions to localize putative insertion/deletion differences between the \u003cem\u003eS-\u003c/em\u003e and \u003cem\u003es\u003c/em\u003e-alleles and between the \u003cem\u003eM-\u003c/em\u003e and \u003cem\u003em\u003c/em\u003e-alleles, we used Anchorwave\u003csup\u003e29\u003c/sup\u003e to align the two genome assemblies, with the \u0026ldquo;minimap2\u0026rdquo; option and the \u0026ldquo;proali\u0026rdquo; function.\u003c/p\u003e \u003cp\u003e \u003cb\u003ePopulation genomic statistics\u003c/b\u003e. We calculated the nucleotide diversity (\u003cem\u003eπ\u003c/em\u003e) of 1-Mb sliding windows for the 60 individuals (20 of each morph) used in the GWAS analyses. We used SAMtools\u003csup\u003e41\u003c/sup\u003e to remove the reads with mapping quality\u0026thinsp;\u0026lt;\u0026thinsp;10 in the filtered BAM files generated in the section \u0026ldquo;Genome-wide associations of single nucleotide polymorphisms (SNPs) and read coverage (RCs)\u0026rdquo;. Then we used ANGSD\u003csup\u003e75\u003c/sup\u003e to calculate the nucleotide diversity (\u003cem\u003eπ\u003c/em\u003e) of the sliding windows.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003e \u003c/p\u003e"},{"header":"Declarations","content":"\u003cp\u003e\u003cstrong\u003eAcknowledgements\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThis work was supported by Discovery Grants from the Natural Sciences and Engineering Research Council of Canada (NSERC) to S.C.H.B. and S.I.W. and student fellowships from the University of Toronto to H.X. The Canu assembly and Purge Haplotigs pipeline of the \u003cem\u003eSsMm\u003c/em\u003e genome was conducted by the Canadian Center for Computational Genomics (C3G), which is a Genomics Technology Platform (GTP) supported by the Canadian Government through Genome Canada. We thank William Cole and Tom Gludovacz for glasshouse maintenance, Joanna Tran, SungA Choe, Ryan Ruan, and Jianuo Feng for help with growing plants, Claire Ellis and Po-Yuan Huang for help with tissue collection, Baharul Choudhury, Iuliia Pelikh, and Yihe Zheng for help with DNA extraction, and Joanna Rifkin, Christian Kappel, Mathias Scharmann, Julia Kreiner, Jingshan Yang, Meng Yuan, Liangjiao Xue, Zoë Humphries, Tyler Kent, and Ting Liu for advice on data analysis.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAuthor contributions\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eH.X., S.I.W., and S.C.H.B. conceived and designed the study and wrote the manuscript. S.C.H.B. obtained the seeds used for the study. H.X. grew the plants, collected the tissues, and extracted DNA and RNA. Y.G. performed the genome annotation. S.I.W. performed the Anchorwave alignment. H.X. performed all remaining analyses. H.X., S.I.W., and S.C.H.B. wrote the manuscript. All authors read and approved the final manuscript.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eCompeting interests\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eNone declared.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eData Availability Statement\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe PacBio CLR sequences used to assemble the \u003cem\u003eEichhornia paniculata\u003c/em\u003e \u003cem\u003eSsMm\u003c/em\u003e genome are available under Accession no. SRR21038607 at the Sequence Read Archive (SRA) (http://www.ncbi.nlm.nih.gov/sra, accessed 09 February 2024), and the associated BioProject alias is PRJNA863567. The Hi-C data used to assemble the \u003cem\u003eE. paniculata\u003c/em\u003e \u003cem\u003eSsMm\u003c/em\u003e genome are available under Accession no. SRR21038606 at the Sequence Read Archive (SRA) (http://www.ncbi.nlm.nih.gov/sra, accessed 09 February 2024), and the associated BioProject alias is PRJNA863567. The 10X Chromium data used to polish the \u003cem\u003eE. paniculata SsMm\u003c/em\u003e genome are available under Accession no. SRR21038605 at the Sequence Read Archive (SRA) (http://www. ncbi.nlm.nih.gov/sra, accessed 09 February 2024), and the associated BioProject alias is PRJNA863567. The PacBio CCS sequences used to assemble the \u003cem\u003eE. paniculata ssmm\u003c/em\u003e genome are available under Accession no. SRRXXXXXXXX at the Sequence Read Archive (SRA) (http://www.ncbi.nlm.nih.gov/sra, accessed XX XX XXXX), and the associated BioProject alias is PRJNAXXXXXX. The Omni-C data used to assemble the \u003cem\u003eE. paniculata ssmm\u003c/em\u003e genome are available under Accession no. SRRXXXXXX at the Sequence Read Archive (SRA) (http://www.ncbi.nlm.nih.gov/sra, accessed XX XX XXXX), and the associated BioProject alias is PRJNAXXXXXX. The whole genome Illumina sequences of the 60 \u003cem\u003eE. paniculata\u003c/em\u003e individuals in this study are available under Accession no. SRRXXXXXXXX - SRRXXXXXXXX at the Sequence Read Archive (SRA) (http://www.ncbi.nlm.nih.gov/sra, accessed XX XX XXXX), and the associated BioProject alias is PRJNA862356. The RNA-Seq data of styles and stamens are available under Accession no. SRRXXXXXXXX - SRRXXXXXXXX at the Sequence Read Archive (SRA) (http://www.ncbi.nlm.nih.gov/sra, accessed XX XX XXXX), and the associated BioProject alias is PRJNAXXXXXX. The \u003cem\u003eE. paniculata SsMm\u003c/em\u003e genome and its MAKER and InterProScan annotation are available under Genome ID XXXXX on CoGe (https://genomevolution.org/coge/, accessed XX XX XXXX). The \u003cem\u003eE. paniculata ssmm\u003c/em\u003e genome and its MAKER and InterProScan annotation are available under Genome ID XXXXX on CoGe (https://genomevolution.org/coge/, accessed XX XX XXXX).\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\n\u003cli\u003eDarwin, C. R. \u003cem\u003eThe different forms of flowers on plants of the same species.\u003c/em\u003e (John Murray, 1877).\u003c/li\u003e\n\u003cli\u003eBarrett, S. C. H. Heterostylous genetic polymorphisms: model systems for evolutionary analysis. In \u003cem\u003eEvolution and function of heterostyly\u003c/em\u003e (ed. Barrett, S. C. H.) 1\u0026ndash;29 (Springer, 1992).\u003c/li\u003e\n\u003cli\u003eBarrett, S. C. H. The evolution of plant sexual diversity. \u003cem\u003eNat. Rev. Genet.\u003c/em\u003e \u003cstrong\u003e3\u003c/strong\u003e, 274-284 (2002).\u003c/li\u003e\n\u003cli\u003eLloyd, D. G. \u0026amp; Webb, C. The selection of heterostyly. In \u003cem\u003eEvolution and function of heterostyly\u003c/em\u003e (ed. Barrett, S. C. H.) 179\u0026ndash;207 (Springer, 1992).\u003c/li\u003e\n\u003cli\u003eLloyd, D. G. \u0026amp; Webb, C. J. The evolution of heterostyly. In \u003cem\u003eEvolution and function of heterostyly\u003c/em\u003e (ed. Barrett, S. C. H.) 151\u0026ndash;178 (Springer, 1992).\u003c/li\u003e\n\u003cli\u003eBarrett, S. C. H., Jesson, L. K. \u0026amp; Baker, A. M. The evolution and function of stylar polymorphisms in flowering plants. \u003cem\u003eAnn. Bot.\u003c/em\u003e \u003cstrong\u003e85\u003c/strong\u003e, 253-265 (2000).\u003c/li\u003e\n\u003cli\u003eLewis, D. \u0026amp; Jones, D. The genetics of heterostyly. In \u003cem\u003eEvolution and function of heterostyly\u003c/em\u003e (Ed. Barrett, S. C. H.) 129\u0026ndash;150 (Springer, 1992).\u003c/li\u003e\n\u003cli\u003eMather, K. The genetical architecture of heterostyly in \u003cem\u003ePrimula sinensis\u003c/em\u003e. \u003cem\u003eEvolution\u003c/em\u003e \u003cstrong\u003e4\u003c/strong\u003e, 340\u0026ndash;352 (1950).\u003c/li\u003e\n\u003cli\u003eSchwander, T., Libbrecht, R. \u0026amp; Keller, L. Supergenes and complex phenotypes. \u003cem\u003eCurrent Biology\u003c/em\u003e \u003cstrong\u003e24\u003c/strong\u003e, 288\u0026ndash;294 (2014).\u003c/li\u003e\n\u003cli\u003eThompson, M. J. \u0026amp; Jiggins, C. Supergenes and their role in evolution. \u003cem\u003eHeredity\u003c/em\u003e \u003cstrong\u003e113\u003c/strong\u003e, 1-8 (2014).\u003c/li\u003e\n\u003cli\u003eManey, D. L. \u0026amp; K\u0026uuml;pper, C. Supergenes on steroids. \u003cem\u003ePhilos. Trans. R. Soc. B\u003c/em\u003e \u003cstrong\u003e377\u003c/strong\u003e, 20200507 (2022).\u003c/li\u003e\n\u003cli\u003eKappel, C., Huu, C. N. \u0026amp; Lenhard, M. A short story gets longer: recent insights into the molecular basis of heterostyly. \u003cem\u003eJ. Exp. Bot.\u003c/em\u003e \u003cstrong\u003e68\u003c/strong\u003e, 5719-5730 (2017).\u003c/li\u003e\n\u003cli\u003eBarrett, S. C. H. \u0026lsquo;A most complex marriage arrangement\u0026rsquo;: recent advances on heterostyly and unresolved questions. \u003cem\u003eNew Phytol.\u003c/em\u003e \u003cstrong\u003e224\u003c/strong\u003e, 1051-1067 (2019).\u003c/li\u003e\n\u003cli\u003eHuu, C. N. \u003cem\u003eet al.\u003c/em\u003e Presence versus absence of CYP734A50 underlies the style-length dimorphism in primroses. \u003cem\u003eeLife\u003c/em\u003e \u003cstrong\u003e5\u003c/strong\u003e, e17956 (2016).\u003c/li\u003e\n\u003cli\u003eLi, J. \u003cem\u003eet al.\u003c/em\u003e Genetic architecture and evolution of the \u003cem\u003eS\u003c/em\u003e locus supergene in \u003cem\u003ePrimula vulgaris\u003c/em\u003e. \u003cem\u003eNat. Plants\u003c/em\u003e \u003cstrong\u003e2\u003c/strong\u003e, 1-7 (2016).\u003c/li\u003e\n\u003cli\u003eShore, J. S. \u003cem\u003eet al.\u003c/em\u003e The long and short of the \u003cem\u003eS\u003c/em\u003e-locus in \u003cem\u003eTurnera\u003c/em\u003e (Passifloraceae). \u003cem\u003eNew Phytol.\u003c/em\u003e \u003cstrong\u003e224\u003c/strong\u003e, 1316\u0026ndash;1329 (2019).\u003c/li\u003e\n\u003cli\u003eGuti\u0026eacute;rrez-Valencia, J. \u003cem\u003eet al.\u003c/em\u003e Genomic analyses of the \u003cem\u003eLinum\u003c/em\u003e distyly supergene reveal convergent evolution at the molecular level. \u003cem\u003eCurr. Biol.\u003c/em\u003e \u003cstrong\u003e32\u003c/strong\u003e, 4360\u0026ndash;4371 (2022).\u003c/li\u003e\n\u003cli\u003eZhao, Z. \u003cem\u003eet al.\u003c/em\u003e Genomic evidence supports the genetic convergence of a supergene controlling the distylous floral syndrome. \u003cem\u003eNew Phytol.\u003c/em\u003e \u003cstrong\u003e237\u003c/strong\u003e, 601\u0026ndash;614 (2023).\u003c/li\u003e\n\u003cli\u003eYang, J. \u003cem\u003eet al.\u003c/em\u003e Haplotype-resolved genome assembly provides insights into the evolution of \u003cem\u003eS\u003c/em\u003e-locus supergene in distylous \u003cem\u003eNymphoides indica\u003c/em\u003e. \u003cem\u003eNew Phytol.\u003c/em\u003e \u003cstrong\u003e240\u003c/strong\u003e, 2058-2071 (2023).\u003c/li\u003e\n\u003cli\u003eFawcett, J. A. \u003cem\u003eet al.\u003c/em\u003e Genome sequencing reveals the genetic architecture of heterostyly and domestication history of common buckwheat. \u003cem\u003eNature Plants\u003c/em\u003e \u003cstrong\u003e9\u003c/strong\u003e, 1236-1251 (2023).\u003c/li\u003e\n\u003cli\u003ePotente, G. \u003cem\u003eet al.\u003c/em\u003e Comparative genomics elucidates the origin of a supergene controlling floral heteromorphism. \u003cem\u003eMol. Biol. Evol.\u003c/em\u003e \u003cstrong\u003e39\u003c/strong\u003e, msac035 (2022).\u003c/li\u003e\n\u003cli\u003eBarrett, S. C. H. The evolutionary biology of tristyly. In \u003cem\u003eOxford Surveys in Evolutionary Biology\u003c/em\u003e (Eds. D. Futuyma \u0026amp; J. Antonovics) vol. \u003cstrong\u003e9\u003c/strong\u003e, 283\u0026ndash;326 (Oxford University Press, 1993).\u003c/li\u003e\n\u003cli\u003eThompson, J. D., Pailler, T., Strasberg, D. \u0026amp; Manicacci, D. Tristyly in the endangered Mascarene Island endemic \u003cem\u003eHugonia serrata\u003c/em\u003e (Linaceae). \u003cem\u003eAm. J. Bot.\u003c/em\u003e \u003cstrong\u003e83\u003c/strong\u003e, 1160\u0026ndash;1167 (1996).\u003c/li\u003e\n\u003cli\u003eCharlesworth, D. The evolution and breakdown of tristyly. \u003cem\u003eEvolution\u003c/em\u003e \u003cstrong\u003e33\u003c/strong\u003e, 489-498. (1979)\u003c/li\u003e\n\u003cli\u003eArunkumar, R., Wang, W., Wright, S. I. \u0026amp; Barrett, S. C. H. The genetic architecture of tristyly and its breakdown to self-fertilization. \u003cem\u003eMol. Ecol.\u003c/em\u003e \u003cstrong\u003e26\u003c/strong\u003e, 752\u0026ndash;765 (2017).\u003c/li\u003e\n\u003cli\u003eOssowski, S. \u003cem\u003eet al.\u003c/em\u003e The rate and molecular spectrum of spontaneous mutations in \u003cem\u003eArabidopsis thaliana\u003c/em\u003e. \u003cem\u003eScience\u003c/em\u003e \u003cstrong\u003e327\u003c/strong\u003e, 92\u0026ndash;94 (2010).\u003c/li\u003e\n\u003cli\u003eEckenwalder, J. E. \u0026amp; Barrett, S. C. H. Phylogenetic systematics of Pontederiaceae. \u003cem\u003eSyst. Bot.\u003c/em\u003e \u003cstrong\u003e11\u003c/strong\u003e, 373\u0026ndash;391 (1986).\u003c/li\u003e\n\u003cli\u003eNess, R. W., Siol, M. \u0026amp; Barrett, S. C. H. \u003cem\u003eDe novo\u003c/em\u003e sequence assembly and characterization of the floral transcriptome in cross-and self-fertilizing plants. \u003cem\u003eBMC Genomics\u003c/em\u003e \u003cstrong\u003e12\u003c/strong\u003e, 298 (2011).\u003c/li\u003e\n\u003cli\u003eSong, B. \u003cem\u003eet al.\u003c/em\u003e AnchorWave: Sensitive alignment of genomes with high sequence diversity, extensive structural polymorphism, and whole-genome duplication. \u003cem\u003eProc. Natl. Acad. Sci.\u003c/em\u003e \u003cstrong\u003e119\u003c/strong\u003e, e2113075119 (2022).\u003c/li\u003e\n\u003cli\u003eJiao, Z., Du, H., Chen, S., Huang, W. \u0026amp; Ge, L. \u003cem\u003eLAZY\u003c/em\u003e gene family in plant gravitropism. \u003cem\u003eFront. Plant Sci.\u003c/em\u003e \u003cstrong\u003e11\u003c/strong\u003e, 606241 (2021).\u003c/li\u003e\n\u003cli\u003eChen, Y. \u003cem\u003eet al.\u003c/em\u003e Identification of the abundant hydroxyproline-rich glycoproteins in the root walls of wild-type \u003cem\u003eArabidopsis\u003c/em\u003e, an ext3 mutant line, and its phenotypic revertant. \u003cem\u003ePlants\u003c/em\u003e \u003cstrong\u003e4\u003c/strong\u003e, 85\u0026ndash;111 (2015).\u003c/li\u003e\n\u003cli\u003eSozzani, R. \u003cem\u003eet al.\u003c/em\u003e Spatiotemporal regulation of cell-cycle genes by SHORTROOT links patterning and growth. \u003cem\u003eNature\u003c/em\u003e, \u003cstrong\u003e466\u003c/strong\u003e, 128-132 (2021).\u003c/li\u003e\n\u003cli\u003eGuti\u0026eacute;rrez-Valencia, J. \u003cem\u003eet al.\u003c/em\u003e The genomic architecture and evolutionary fates of supergenes. \u003cem\u003eGenome Biology and Evolution\u003c/em\u003e \u003cstrong\u003e13\u003c/strong\u003e, evab057 (2021).\u003c/li\u003e\n\u003cli\u003eBouckaert, R. \u003cem\u003eet al.\u003c/em\u003e BEAST 2.5: An advanced software platform for Bayesian evolutionary analysis. \u003cem\u003ePLoS Comput. Biol\u003c/em\u003e\u003cstrong\u003e.\u003c/strong\u003e \u003cstrong\u003e15\u003c/strong\u003e, e1006650 (2019).\u003c/li\u003e\n\u003cli\u003eNess, R. W., Wright, S. I. \u0026amp; Barrett, S. C. H. Mating-system variation, demographic history and patterns of nucleotide diversity in the tristylous plant \u003cem\u003eEichhornia paniculata\u003c/em\u003e. \u003cem\u003eGenetics\u003c/em\u003e \u003cstrong\u003e184\u003c/strong\u003e, 381\u0026ndash;392 (2010).\u003c/li\u003e\n\u003cli\u003eKoren, S. \u003cem\u003eet al.\u003c/em\u003e Canu: scalable and accurate long-read assembly via adaptive k-mer weighting and repeat separation. \u003cem\u003eGenome Res\u003c/em\u003e. \u003cstrong\u003e27\u003c/strong\u003e, 722\u0026ndash;736 (2017).\u003c/li\u003e\n\u003cli\u003eRoach, M. J., Schmidt, S. A., \u0026amp; Borneman, A. R. Purge Haplotigs: allelic contig reassignment for third-gen diploid genome assemblies. \u003cem\u003eBMC bioinformatics\u003c/em\u003e \u003cstrong\u003e19\u003c/strong\u003e, 460 (2018).\u003c/li\u003e\n\u003cli\u003eKronenberg, Z. N\u003cem\u003e. et al.\u003c/em\u003e Extended haplotype-phasing of long-read de novo genome assemblies using Hi-C. \u003cem\u003eNature Communications\u003c/em\u003e \u003cstrong\u003e12\u003c/strong\u003e, 1935 (2021).\u003c/li\u003e\n\u003cli\u003eThress, K., Henzel, W., Shillinglaw, W. \u0026amp; Kornbluth, S. Scythe: a novel reaper-binding apoptotic regulator. \u003cem\u003eEMBO J.\u003c/em\u003e \u003cstrong\u003e17\u003c/strong\u003e, 6135\u0026ndash;6143 (1998).\u003c/li\u003e\n\u003cli\u003eSedlazeck, F. J., Rescheneder, P. \u0026amp; Von Haeseler, A. NextGenMap: fast and accurate read mapping in highly polymorphic genomes. \u003cem\u003eBioinformatics \u003c/em\u003e\u003cstrong\u003e29\u003c/strong\u003e, 2790\u0026ndash;2791 (2013).\u003c/li\u003e\n\u003cli\u003eDanecek, P. \u003cem\u003eet al.\u003c/em\u003e Twelve years of SAMtools and BCFtools. \u003cem\u003eGigaScience\u003c/em\u003e \u003cstrong\u003e10\u003c/strong\u003e, (2021).\u003c/li\u003e\n\u003cli\u003eRastas, P. Lep-MAP3: robust linkage mapping even for low-coverage whole genome sequencing data. \u003cem\u003eBioinformatics\u003c/em\u003e \u003cstrong\u003e33\u003c/strong\u003e, 3726\u0026ndash;3732 (2017).\u003c/li\u003e\n\u003cli\u003eDurand, N. C. \u003cem\u003eet al.\u003c/em\u003e Juicebox provides a visualization system for Hi-C contact maps with unlimited zoom. \u003cem\u003eCell Syst.\u003c/em\u003e \u003cstrong\u003e3\u003c/strong\u003e, 99\u0026ndash;101 (2016).\u003c/li\u003e\n\u003cli\u003eWalker, B. J. \u003cem\u003eet al.\u003c/em\u003e Pilon: an integrated tool for comprehensive microbial variant detection and genome assembly improvement. \u003cem\u003ePLoS One\u003c/em\u003e \u003cstrong\u003e9\u003c/strong\u003e, e112963 (2014).\u003c/li\u003e\n\u003cli\u003eSmit, A. F. A. \u0026amp; Hubley, R. RepeatModeler Open-1.0. 2008\u0026ndash;2015. (2008).\u003c/li\u003e\n\u003cli\u003eSmit, A. F. A., Hubley, R. \u0026amp; Green, P. RepeatMasker Open-4.0. (2013).\u003c/li\u003e\n\u003cli\u003eCheng, H., Concepcion, G. T., Feng, X., Zhang, H. \u0026amp; Li, H. Haplotype-resolved de novo assembly using phased assembly graphs with hifiasm. \u003cem\u003eNat. Methods\u003c/em\u003e \u003cstrong\u003e18\u003c/strong\u003e, 170\u0026ndash;175 (2021).\u003c/li\u003e\n\u003cli\u003ePutnam, N. H. \u003cem\u003eet al.\u003c/em\u003e Chromosome-scale shotgun assembly using an \u003cem\u003ein vitro\u003c/em\u003e method for long-range linkage. \u003cem\u003eGenome Res.\u003c/em\u003e \u003cstrong\u003e26\u003c/strong\u003e, 342\u0026ndash;350 (2016).\u003c/li\u003e\n\u003cli\u003eGurevich, A., Saveliev, V., Vyahhi, N. \u0026amp; Tesler, G. QUAST: quality assessment tool for genome assemblies. \u003cem\u003eBioinformatics\u003c/em\u003e \u003cstrong\u003e29\u003c/strong\u003e, 1072\u0026ndash;1075 (2013).\u003c/li\u003e\n\u003cli\u003eManni, M. \u003cem\u003eet al.\u003c/em\u003e BUSCO update: novel and streamlined workflows along with broader and deeper phylogenetic coverage for scoring of eukaryotic, prokaryotic, and viral genomes. \u003cem\u003eMol. Biol. Evol.\u003c/em\u003e \u003cstrong\u003e38\u003c/strong\u003e, 4647-4654 (2021).\u003c/li\u003e\n\u003cli\u003eLevy Karin, E., Mirdita, M., \u0026amp; S\u0026ouml;ding, J. MetaEuk\u0026mdash;sensitive, high-throughput gene discovery, and annotation for large-scale eukaryotic metagenomics. \u003cem\u003eMicrobiome\u003c/em\u003e \u003cstrong\u003e8\u003c/strong\u003e, 48 (2020).\u003c/li\u003e\n\u003cli\u003eSedlazeck, F. J., Rescheneder, P. \u0026amp; Von Haeseler, A. NextGenMap: fast and accurate read mapping in highly polymorphic genomes. \u003cem\u003eBioinformatics\u003c/em\u003e \u003cstrong\u003e29\u003c/strong\u003e, 2790-2791 (2013).\u003c/li\u003e\n\u003cli\u003eGarrison, E. \u0026amp; Marth, G. Haplotype-based variant detection from short-read sequencing. \u003cem\u003eArXiv Prepr.\u003c/em\u003e, ArXiv12073907 (2012).\u003c/li\u003e\n\u003cli\u003eTange, O. Gnu parallel: the command-line power tool. \u003cem\u003eUsenix Mag\u003c/em\u003e.\u003cstrong\u003e 36\u003c/strong\u003e, 42 (2011).\u003c/li\u003e\n\u003cli\u003eDanecek, P. \u003cem\u003eet al.\u003c/em\u003e The variant call format and VCFtools. \u003cem\u003eBioinformatics\u003c/em\u003e \u003cstrong\u003e27\u003c/strong\u003e, 2156\u0026ndash;2158 (2011).\u003c/li\u003e\n\u003cli\u003eGarrison, E., Kronenberg, Z. N., Dawson, E. T., Pedersen, B. S. \u0026amp; Prins, P. A spectrum of free software tools for processing the VCF variant call format: vcflib, bio-vcf, cyvcf2, hts-nim and slivar. \u003cem\u003ePLoS Comput. Biol.\u003c/em\u003e \u003cstrong\u003e18\u003c/strong\u003e, e1009123 (2022).\u003c/li\u003e\n\u003cli\u003ePuritz, J. B., Hollenbeck, C. M. \u0026amp; Gold, J. R. dDocent: a RADseq, variant-calling pipeline designed for population genomics of non-model organisms. \u003cem\u003ePeerJ\u003c/em\u003e \u003cstrong\u003e2\u003c/strong\u003e, e431 (2014).\u003c/li\u003e\n\u003cli\u003eZhou, X. \u0026amp; Stephens, M. Genome-wide efficient mixed-model analysis for association studies. \u003cem\u003eNat. Genet\u003c/em\u003e. \u003cstrong\u003e44\u003c/strong\u003e, 821\u0026ndash;824 (2012).\u003c/li\u003e\n\u003cli\u003eXue, L. \u003cem\u003eet al.\u003c/em\u003e Evidences for a role of two Y-specific genes in sex determination in \u003cem\u003ePopulus deltoides\u003c/em\u003e. \u003cem\u003eNature Communications\u003c/em\u003e \u003cstrong\u003e11\u003c/strong\u003e, 5893 (2020).\u003c/li\u003e\n\u003cli\u003ePedersen, B. S. \u0026amp; Quinlan, A. R. Mosdepth: quick coverage calculation for genomes and exomes. \u003cem\u003eBioinformatics\u003c/em\u003e \u003cstrong\u003e34\u003c/strong\u003e, 867\u0026ndash;868 (2018).\u003c/li\u003e\n\u003cli\u003eKim, D., Langmead, B. \u0026amp; Salzberg, S. L. HISAT: a fast spliced aligner with low memory requirements. \u003cem\u003eNat. Methods\u003c/em\u003e \u003cstrong\u003e12\u003c/strong\u003e, 357\u0026ndash;360 (2015).\u003c/li\u003e\n\u003cli\u003ePertea, M. \u003cem\u003eet al.\u003c/em\u003e StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. \u003cem\u003eNat. Biotechnol.\u003c/em\u003e \u003cstrong\u003e33\u003c/strong\u003e, 290\u0026ndash;295 (2015).\u003c/li\u003e\n\u003cli\u003eLohse, M\u003cem\u003e. et al.\u003c/em\u003e Mercator: a fast and simple web server for genome scale functional annotation of plant sequence data. \u003cem\u003ePlant, Cell \u0026amp; Environment\u003c/em\u003e \u003cstrong\u003e37\u003c/strong\u003e, 1250-1258 (2014).\u003c/li\u003e\n\u003cli\u003eLove, M. I., Huber, W. \u0026amp; Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. \u003cem\u003eGenome Biol.\u003c/em\u003e \u003cstrong\u003e15\u003c/strong\u003e, 1\u0026ndash;21 (2014).\u003c/li\u003e\n\u003cli\u003eArunkumar, R., Ness, R. W., Wright, S. I. \u0026amp; Barrett, S. C. H. The evolution of selfing is accompanied by reduced efficacy of selection and purging of deleterious mutations. \u003cem\u003eGenetics \u003c/em\u003e\u003cstrong\u003e199\u003c/strong\u003e, 817\u0026ndash;829 (2015).\u003c/li\u003e\n\u003cli\u003eArunkumar, R., Maddison, T. I., Barrett, S. C. H. \u0026amp; Wright, S. I. Recent mating-system evolution in \u003cem\u003eEichhornia\u003c/em\u003e is accompanied by cis-regulatory divergence. \u003cem\u003eNew Phytol.\u003c/em\u003e \u003cstrong\u003e211\u003c/strong\u003e, 697\u0026ndash;707 (2016).\u003c/li\u003e\n\u003cli\u003eCantarel, B. L. \u003cem\u003eet al.\u003c/em\u003e MAKER: an easy-to-use annotation pipeline designed for emerging model organism genomes. \u003cem\u003eGenome Res.\u003c/em\u003e \u003cstrong\u003e18\u003c/strong\u003e, 188\u0026ndash;196 (2008).\u003c/li\u003e\n\u003cli\u003eJones, P. \u003cem\u003eet al.\u003c/em\u003e InterProScan 5: genome-scale protein function classification. \u003cem\u003eBioinformatics\u003c/em\u003e \u003cstrong\u003e30\u003c/strong\u003e, 1236\u0026ndash;1240 (2014).\u003c/li\u003e\n\u003cli\u003eKatoh, K. \u0026amp; Standley, D. M. MAFFT Multiple Sequence Alignment Software Version 7: Improvements in Performance and Usability. \u003cem\u003eMol. Biol. Evol.\u003c/em\u003e \u003cstrong\u003e30\u003c/strong\u003e, 772-780 (2013).\u003c/li\u003e\n\u003cli\u003ePosada, D. jModelTest: Phylogenetic Model Averaging. \u003cem\u003eMol. Biol. Evol. \u003c/em\u003e\u003cstrong\u003e25\u003c/strong\u003e, 1253-1256 (2008).\u003c/li\u003e\n\u003cli\u003eHuang, W. \u003cem\u003eet al.\u003c/em\u003e A well-supported nuclear phylogeny of Poaceae and implications for the evolution of C4 photosynthesis. \u003cem\u003eMol. Plant\u003c/em\u003e \u003cstrong\u003e15\u003c/strong\u003e, 755-777 (2022).\u003c/li\u003e\n\u003cli\u003eMagall\u0026oacute;n, S., G\u0026oacute;mez-Acevedo, S., S\u0026aacute;nchez-Reyes, L. L. \u0026amp; Hern\u0026aacute;ndez-Hern\u0026aacute;ndez, T. A metacalibrated time-tree documents the early rise of flowering plant phylogenetic diversity. \u003cem\u003eNew Phytol.\u003c/em\u003e \u003cstrong\u003e207\u003c/strong\u003e, 437-453 (2015).\u003c/li\u003e\n\u003cli\u003eCamacho, C. \u003cem\u003eet al.\u003c/em\u003e BLAST+: architecture and applications. \u003cem\u003eBMC Bioinformatics\u003c/em\u003e \u003cstrong\u003e10\u003c/strong\u003e, 421 (2009).\u003c/li\u003e\n\u003cli\u003eHaug-Baltzell, A., Stephens, S. A., Davey, S., Scheidegger, C. E. \u0026amp; Lyons, E. SynMap2 and SynMap3D: web-based whole-genome synteny browsers. \u003cem\u003eBioinformatics\u003c/em\u003e \u003cstrong\u003e33\u003c/strong\u003e, 2197\u0026ndash;2198 (2017).\u003c/li\u003e\n\u003cli\u003eKorneliussen, T. S., Albrechtsen, A., \u0026amp; Nielsen, R. ANGSD: analysis of next generation sequencing data. \u003cem\u003eBMC bioinformatics\u003c/em\u003e \u003cstrong\u003e15\u003c/strong\u003e, 356 (2014).\u003c/li\u003e\n\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":true,"hideJournal":true,"highlight":"","institution":"","isAcceptedByJournal":false,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true},"keywords":"tristyly, plants, genomics, supergene","lastPublishedDoi":"10.21203/rs.3.rs-4232430/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-4232430/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eTristyly is a polymorphism characterized by three flower morphs with reciprocal stigma and anther heights controlled by two epistatically interacting diallelic loci (\u003cem\u003eS\u003c/em\u003e and \u003cem\u003eM\u003c/em\u003e), hypothesized to be supergenes. Chromosome-level genome assemblies of \u003cem\u003eEichhornia paniculata \u003c/em\u003eidentified the \u003cem\u003eS-\u003c/em\u003e and \u003cem\u003eM-\u003c/em\u003eloci. The \u003cem\u003eS\u003c/em\u003e-locus is a supergene consisting of two divergent alleles: The \u003cem\u003eS\u003c/em\u003e-allele (2.51Mb) with three \u003cem\u003eS\u003c/em\u003e-allele specific genes hemizygous in most S-morph plants and the \u003cem\u003es\u003c/em\u003e-allele (596kb) with five \u003cem\u003es\u003c/em\u003e-allele specific genes. Two of the \u003cem\u003eS\u003c/em\u003e-allele specific genes, LAZY1-S and HRGP-S, were specifically expressed in styles and stamens, respectively, making them tristyly candidate genes. The \u003cem\u003eM\u003c/em\u003e-locus contained one gene (LAZY1-M), homologous to LAZY1-S, that was present in the \u003cem\u003eM\u003c/em\u003e-allele but absent from the \u003cem\u003em\u003c/em\u003e-allele. Estimates of allele ages are consistent with the prediction that the \u003cem\u003eS-\u003c/em\u003elocus evolved before the \u003cem\u003eM\u003c/em\u003e-locus. Re-use of the same gene family highlights the potential role of gene duplication in the evolution of epistatic multilocus polymorphisms.\u003c/p\u003e","manuscriptTitle":"Genomic evidence for supergene control of Darwin’s “complex marriage arrangement” – the tristylous floral polymorphism","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2024-04-11 05:10:40","doi":"10.21203/rs.3.rs-4232430/v1","editorialEvents":[{"type":"communityComments","content":0}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"ea1654e9-2c9c-4315-b593-70fb7cfde4a3","owner":[],"postedDate":"April 11th, 2024","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"posted","subjectAreas":[{"id":30401336,"name":"Biological sciences/Evolution"},{"id":30401337,"name":"Biological sciences/Plant sciences/Plant evolution"},{"id":30401338,"name":"Biological sciences/Plant sciences/Plant genetics"}],"tags":[],"updatedAt":"2025-01-09T02:20:49+00:00","versionOfRecord":[],"versionCreatedAt":"2024-04-11 05:10:40","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-4232430","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-4232430","identity":"rs-4232430","version":["v1"]},"buildId":"8U1c8b4HqxoKbykW_rLl7","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}