Genomic tree scans identify loci underlying adaptive peaks in Antirrhinum

Genomic tree scans identify loci underlying adaptive peaks in Antirrhinum | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (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];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Genomic tree scans identify loci underlying adaptive peaks in Antirrhinum Daniel M. Richardson , Desmond Bradley , Lucy Copsey , Annabel Whibley , Monique Burrus , Christophe Andalo , Sihui Zhu , David L. Field , Yongbiao Xue , View ORCID Profile Enrico Coen doi: https://doi.org/10.1101/2025.02.12.637406 Daniel M. Richardson 1 Department of Cell and Developmental Biology, John Innes Centre , Colney Lane, Norwich, NR4 7UH Find this author on Google Scholar Find this author on PubMed Search for this author on this site Desmond Bradley 1 Department of Cell and Developmental Biology, John Innes Centre , Colney Lane, Norwich, NR4 7UH Find this author on Google Scholar Find this author on PubMed Search for this author on this site Lucy Copsey 1 Department of Cell and Developmental Biology, John Innes Centre , Colney Lane, Norwich, NR4 7UH Find this author on Google Scholar Find this author on PubMed Search for this author on this site Annabel Whibley 2 School of Biological Sciences, University of Auckland , Auckland, New Zealand Find this author on Google Scholar Find this author on PubMed Search for this author on this site Monique Burrus 3 Centre de recherches sur la biodiversité et l’environnement CRBE, CNRS, UMR 5300, Université Toulouse III Paul Sabatier , F-31062, Toulouse, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site Christophe Andalo 3 Centre de recherches sur la biodiversité et l’environnement CRBE, CNRS, UMR 5300, Université Toulouse III Paul Sabatier , F-31062, Toulouse, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site Sihui Zhu 4 Center for Genomics and Biotechnology, Fujian Provincial Key Laboratory of Haixia Applied Plant Systems Biology, Key Laboratory of Genetics, Breeding and Multiple Utilization of Corps, Ministry of Education, Fujian Agriculture and Forestry University , Fuzhou, 350002, China Find this author on Google Scholar Find this author on PubMed Search for this author on this site David L. Field 5 Applied BioSciences, Macquarie University , NSW 2109, Australia Find this author on Google Scholar Find this author on PubMed Search for this author on this site For correspondence: enrico.coen{at}jic.ac.uk ybxue{at}genetics.ac.cn david.field{at}mq.edu.au Yongbiao Xue 6 Institute of Genetics and Developmental Biology; Chinese Academy of Sciences , Beijing, China Find this author on Google Scholar Find this author on PubMed Search for this author on this site For correspondence: enrico.coen{at}jic.ac.uk ybxue{at}genetics.ac.cn david.field{at}mq.edu.au Enrico Coen 1 Department of Cell and Developmental Biology, John Innes Centre , Colney Lane, Norwich, NR4 7UH Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Enrico Coen For correspondence: enrico.coen{at}jic.ac.uk ybxue{at}genetics.ac.cn david.field{at}mq.edu.au Abstract Full Text Info/History Metrics Supplementary material Preview PDF Abstract Species can be considered to occupy different peaks in a fitness landscape, visualised as an undulating surface, with height corresponding to fitness and horizontal coordinates to genotype 1 . A key evolutionary problem is to understand how populations traverse fitness valleys to reach different peaks. Hypothetical solutions include shifting balance theory 2 , peak movements due to environmental changes 1 , 3 , or peak connectivity in high-dimensional fitness spaces 4 , 5 . However, testing hypotheses is not straightforward because identifying loci and epistatic interactions underlying fitness peaks remains challenging 6 . Here we analyse the genetic architecture of a fitness landscape by studying barrier loci that reduce gene flow between populations. We identify barrier loci in an unbiased way, by scanning the genome for regions with deep-rooted trees with a common topology. Applying this method to populations of Antirrhinum , we identify seven barrier loci, all of which control flower colour. We propose these loci underlie two peaks in a fitness landscape, the heights and positions of which vary depending on genotype-environment interactions, allowing them to have been accessed without traversing a fitness valley. Our findings provide a model for how reproductive barriers between species can arise through populations taking diverse paths in fitness space. Main One approach to identifying loci underlying peaks in a fitness landscape is to study subspecies or races in which reproductive isolation is incomplete. Loci contributing to fitness peaks are expected to exhibit reduced gene flow between populations, as hybrid genotypes fall into fitness valleys. Candidate genes conferring reduced gene flow, termed barrier loci, have been identified 7 , 8 , 9 , 10 , 11 . However, barrier loci may reflect different processes, such as adaptation to different environments or genetic incompatibilities 12 , so the relationship between barrier loci and fitness landscapes is unclear. One limitation of previous studies is that they typically involve comparisons between populations with clear phenotypic differences, so they may bias analysis towards traits that are readily observed. If spatial heterogeneity in the environment underlies barrier loci, populations may be subdivided in multiple ways, by traits not readily visible to the human eye. A further limitation is that by using relative divergence (F ST ) as the primary metric for pairwise genome scans, barrier loci may be confounded with loci undergoing differential selective sweeps among populations that are not in reproductive contact 13 . The limitations of pairwise genome scans have been addressed by using phylogenetic methods to compare multiple taxa simultaneously 14 , 15 , 16 , 17 , 18 . However, these methods can only be applied to a small number of taxa, typically subsets of four, selected based on phenotype, ecology, or geography. Here we develop an aphenotypic approach that allows trees of 20 or more taxa or populations to be compared, using absolute divergence (d XY ) 19 as the primary metric, which is insensitive to selective sweeps when calculated over large genomic regions. By utilising the linear correlation between trees, it is possible to determine their similarity based on the underlying genetic distance matrices. We apply this method to identify putative barrier loci between populations of Antirrhinum and consider the implications of our findings for genotype, phenotype and fitness spaces. The Antirrhinum genus comprises about 25 Mediterranean species and subspecies, thought to have arisen within the last 5 million years 20 . The species exhibit diverse flower colour patterns, each pattern providing a distinct floral guide by highlighting the position of pollinator entry with contrasting colours. This diversity correlates with ecological niche (growth habit and environment) 21 . Alpine species growing on steep cliffs have a trailing habit and typically have white flowers with a magenta and/or yellow highlight. Ruderal species growing on slopes have an erect habit, and typically have magenta flowers with a yellow highlight, or yellow flowers with a magenta highlight 22 , 23 . Analysis of hybrid zones between two ruderal varieties – yellow-flowered Antirrhinum majus majus var. striatum and magenta-flowered A. m. m. var. pseudomajus – has identified three major barrier loci, ROSEA ( ROS ), ELUTA ( EL ) and SULFUREA ( SULF ), that confer reduced gene flow between populations. ROS and EL encode MYB transcription factors that control the distribution of magenta pigment 24 ; whereas SULF encodes sRNAs that destabilise transcripts of a yellow pigment biosynthetic gene 25 . More recently, two further barrier loci have been identified: FLAVIA ( FLA ), which affects yellow and is the target of SULF 26 , and RUBIA ( RUB ), which affects magenta and likely encodes flavonol synthase 27 . Whole genome trees are polyphyletic To explore the genomic relationship between A. m. m. var. pseudomajus (magenta flowers with yellow highlight) and A. m. m. var. striatum (yellow flowers with magenta highlight), nine populations of each variety were sampled across their range ( Fig. 1a ). A tree based on geographic distance between populations, gave a polyphyletic grouping for the varieties, showing that population phenotype was not strongly associated with geographic proximity ( Fig. 1b ). Download figure Open in new tab Figure 1: Population locations and whole-genome trees. a Topographic map showing the locations of populations of A. m. m. var. pseudomajus (magenta spots) and A. m. m. var. striatum (yellow spots) populations. GPS coordinates and population sizes are shown in Supplementary Table 1. b Tree constructed using pairwise geographic distance between populations, based on GPS coordinates of population sampling sites. Magenta labels are A. m. var. pseudomajus populations, and yellow labels are A. m. m. var. striatum populations. c Whole-genome d XY tree based on pairwise comparison of all populations over all genomic biallelic SNPs. The shortest root branch (SRB) is highlighted in red. d Whole-genome d XY tree after simulating a selective sweep across all polymorphic sites in the MP11 population. E Whole genome D tree based on pairwise comparison of all populations over all genomic biallelic SNPs. f Whole genome D tree after simulating a selective sweep across all polymorphic sites in the MP11 population. g Maximum likelihood whole-genome tree, with inference carried out on the consensus sequences of genomic SNPs. h Maximum likelihood tree generated after simulating a selective sweep in the MP11 population. To construct genomic trees, leaves were sampled from 20-60 individuals from each population and pooled (Supplementary Table 1). DNA was extracted from each pool, sequenced, and mapped to the Antirrhinum majus reference genome 28 . SNPs were filtered to ensure all were biallelic. Whole-genome trees based on d XY ( Fig. 1c ) or Nei’s D ( Fig. 1e ) gave polyphyletic grouping for the varieties, but with different topologies. The Maximum Likelihood (ML) method yielded three clades, two comprising only A. m. m. var. striatum populations, and one comprising A. m. m. var. pseudomajus and an A. m. m. var. striatum population (YP1), located adjacent to a hybrid zone ( Fig. 1g ). To illustrate the sensitivity of these trees to inbreeding or sweeps, one of the populations, MP11, was subjected to a genome-wide sweep (equivalent to a severe bottleneck) by randomly sampling SNPs at each position and setting their frequency to 1 ( i.e. fixing them). The resulting d XY tree showed no change in topology ( Fig. 1d ), whereas the D tree ( Fig. 1f ) and ML tree ( Fig. 1h ) were modified such that MP11 was shifted to become the outgroup. Thus, d XY trees offer the possibility of identifying candidate barrier loci without confounding effects of sweeps and inbreeding. Genome scan identifies an outlier forest Genomic regions harbouring loci that act as barriers between two population groups, A and B, would be expected to give d XY trees with deep divisions between these groups. To screen for such regions without phenotypic bias, we divided the genome into 50 kb windows, with a 25 kb overlap between adjacent windows, yielding 18,556 windows in total. The corresponding 18,556 d XY trees were classified according to their cophenetic correlation coefficient, c , a measure of tree similarity 29 . Calculating c for all pairwise combinations of the 18,556 trees would be computationally intensive. We therefore used a progressive elimination approach in which a random tree, termed the seed tree, was selected and compared to all other trees. All trees with c > 0.5 were assigned to the same group, or forest, as the seed tree, and removed from further comparisons. Another seed tree was then randomly selected from the remaining unclassified trees. This process was iterated until no trees remained, yielding 499 forests. To identify forests containing d XY trees with deep divisions between two groups of populations, the length of the shortest root branch (SRB, highlighted in red in Fig. 1b-e ) was calculated for each tree. A high value of SRB indicates a deeply divided tree (the longest root branch was not used to avoid highlighting trees with a single outlier population). Plotting mean SRB length against forest size revealed one well-populated outlier forest with a mean SRB about 6 times greater than the mean (labelled (i) in Fig. 2a ). This outlier forest contained 45 trees (0.24 % of all d XY trees). The mean tree of this forest bifurcated at its root into two groups, A (black circles) and B (white circles), each comprising nine populations ( Fig. 2b ). By contrast, the mean tree for the largest forest (labelled (ii) in Fig. 2a ) was polyphyletic for A and B ( Fig. 2c ), as was the whole genome tree (Supplementary Fig. 1a). Download figure Open in new tab Figure 2: Identification and analysis of an outlier forest. a Forests from a genomic d XY tree scan, with the number of trees within each forest (forest size) plotted against the mean shortest root branch (SRB). The outlier forest, with the greatest mean SRB of all forests with >3 members, is labelled in red (i), and the largest forest is labelled in blue (ii). b Mean of all trees appearing within the outlier forest (i). Root division of the tree defines two population groups: A (black) and B (white). c Mean of all trees appearing within the largest forest (ii). d Frequency with which different genomic regions occurred within the most outlying forest (based on mean SRB) of size >3, across 100 tree-scan bootstrap replicates. Each labelled peak corresponds to a region present in over 50 % of replicates (grey dotted line). Chromosome indicated by alternating grey-black. Trees summarising the results of the bootstrap analysis are given in Supplementary Fig. 1. Running the above tree-scanning algorithm 100 times, using independent random seed trees, consistently yielded one or two outlier forests showing the highest SRB value and a forest size > 3. In > 50 % of the replicates, 67 genomic windows (0.36 % of all genomic windows) appeared in the most outlying of these forests. These 67 windows yielded a mean tree that was monophyletic for A and B (Supplementary Fig. 1b). Analysis of the 67 individual trees from these windows showed that 29 gave the same A and B root division as above ( i.e. were monophyletic for A and B) (Supplementary Fig. 1c). Other genomic windows gave near-monophyletic root divisions in which one of the A populations was grouped with B, or vice versa . Thus, the 67 genomic windows are candidates for harbouring barrier loci between the A and B population groups, with perhaps recombination or introgression events accounting for modulation in tree structure. Identification of monophyletic islands The 67 genomic 50kb windows giving monophyletic or near-monophyletic trees for the A and B grouping mapped to six chromosome regions, numbered 1-6 ( Fig. 2d , Supplementary Table 2). To obtain a higher resolution picture, we scanned the genome in 5 kb windows for trees that gave a monophyletic or near-monophyletic A/B grouping. Monophyletic windows were clustered into regions ranging from approximately 0.06 Mb to more than 1 Mb ( Fig. 3a, c, e, g ; Supplementary Fig. 3). We refer to these regions as monophyletic islands. Variation in island size may reflect recombination rate and/or the number of barrier loci harboured by an island. Download figure Open in new tab Figure 3: Genome scans for monophyletic density and F ST and allele frequency cline at island 1. a, c, e, g, i Frequencies of monophyletic (purple) and near-monophyletic (blue) SNPs on chromosomes 1, 2, 4, 5, and 6 summed over 50 kb windows with 25 kb overlaps, mapped to the A. majus reference genome. Numbers 1-6 indicate the positions of the respective monophyletic islands. Close-ups of monophyletic islands are shown in Supplementary Fig. 3. b, d, f, h, j Mean F ST for pairwise comparisons of all group A and group B populations, averaged in 10 kb windows with 9 kb overlaps across chromosomes 1, 2, 4, 5, and 6. Chromosomes 3, 7, and 8 are shown in Supplementary Fig. 2. k Cline in allele frequency for SNPs corresponding to the pseudomajus allele of monophyletic island 1 ( CRE ), across a hybrid zone. To determine whether monophyletic islands correlated with F ST peaks, we carried out F ST genomic scans for all 36 A/B population comparisons mapped to an A. m. m. var. pseudomajus genome assembly. The mean value of F ST for all 36 pairwise comparisons was then calculated for each genomic window. Four of the monophyletic islands exhibited F ST peaks, while two (islands 2 and 5) did not ( Fig. 3b,d,f,h,j ). Thus, our tree-scanning approach identified 2 regions harbouring candidate barrier loci that would have been missed through F ST analysis. Monophyletic islands harbour colour loci The above results show that through aphenotypic analysis we identified a single major subdivision of the 18 populations into two groups, A and B, with candidate barrier loci between them mapping to six chromosome regions. Incorporation of phenotypic information showed that the A group corresponded to populations classified as A. m. m. var. pseudomajus , and the B group to A. m. m. var. striatum . Thus, the taxonomic division into varieties 22 reflects the operation of candidate barrier loci located within monophyletic islands comprising 0.39 % of the 508 Mb genome. Because our approach was aphenotypic, the results imply that there are no other major subdivisions causing detectable differences in d XY between populations. Genetic, F ST and cline analysis has shown that regions corresponding to monophyletic islands 3-6 harbour flower colour loci: island 3 harbours FLA 26 , island 4 harbours SULF 25 , island 5 harbours RUB 27 and island 6 harbours ROS and EL 24 . Analysis of allele frequencies in population pools sampled across a hybrid zone indicate that islands 1 and 2 also show geographic clines, confirming they harbour barrier loci 27 . To confirm this finding, and to measure cline width more accurately, we determined allele frequencies at the hybrid zone using SNPs from island 1 26 . These SNPs showed a sharp increase in the frequency of A. m. m. var. pseudomajus alleles from the western to the eastern flank of the hybrid zone ( Fig. 3k ), supporting the hypothesis that this monophyletic island harbours a barrier locus. The barrier loci harboured by islands 1 and 2 may correspond to flower colour loci, or they could affect other traits. To determine which of these hypotheses was correct, we analysed flower colour variation in an F 2 population from a cross between A. m. m. var. striatum and A. m. m. var. pseudomajus . To remove colour variation caused by the ROS , EL , and SULF loci, we genotyped the F 2 for these loci, giving six main genotypes (Supplementary Figs. 4-5). For each genotype, we ranked flowers according to either spread of yellow or magenta intensity and assigned them to high or low bins. We also genotyped for SNPs from islands 1 and 2. SNP allele frequencies were significantly different between the high and low yellow bins ( Fig. 4a-b ) but showed no significant differences between high and low magenta bins ( Fig. 4c-d ). Thus, both islands 1 and 2 likely harbour yellow flower colour loci. Download figure Open in new tab Figure 4: Flower colour ranking for islands 1 and 2. a Frequency of Island 1 genotypes for high and low yellow bins from 52 F 2 plants from A. m. m. var. pseudomajus crossed to A. m. m. var. striatum . Plants were genotyped for ROS , EL , and SULF , and shared genetic backgrounds were grouped together. Photographs within each group were ranked according to yellow intensity, yielding a low and high bin. Frequency of plants homozygous for A. m. m. var. pseudomajus (p) or A. m. m. var. striatum (s) SNPs at monophyletic island 1 were determined for each bin. The p -value reflects the results of a contingency χ² test between the low and high yellow bins. b As (a) but for 59 F 2 plants genotyped for island 2 SNPs. c As (a) but with ranking based on magenta colour intensity rather than yellow. d As (b) with ranking based on magenta colour intensity rather than yellow. e Frequency of island 1 genotypes for yellow quartile bins from an F 4 population with a FLA s /fla p ros s /ros s EL s /EL s colour gene background. Photographs were ranked according to yellow intensity, yielding four quartiles of increasing yellow. The median photograph from each quartile is shown as representative. f As (e) but for island 2 genotypes. Ranked images are given in Supplementary Figs. 4-11. As a further test, we ranked and genotyped F 4 populations that were homozygous for ros s EL s (superscript ‘s’ refers to striatum allele) and segregating for the SNPs at islands 1 or 2. For both islands, the striatum SNP allele frequency was significantly elevated in the highest compared to the lowest yellow quartile ( Fig. 4e and f ). For island 1, striatum homozygotes were found only in the highest two quartiles, indicating that the pseudomajus allele was largely dominant. For island 2, homozygotes for the pseudomajus allele were only found in the lowest two quartiles, indicating the striatum allele was largely dominant. These results were confirmed by replicate rankings, carried out by different observers (Supplementary Figs. 4-11). We propose the names CREMOSA ( CRE ) and AURINA ( AUN ) for the yellow colour loci harboured by islands 1 and 2, respectively. Genes underlying previously identified flower colour loci are differentially expressed between A. m. m. var. pseudomajus and A. m. m. var. striatum 24 , 25 . To identify candidate genes underlying CRE and AUN , we therefore performed differential expression analysis using DESeq2 30 on RNA extracted from two pools of A. m. m. var. pseudomajus plants and three pools of A. m. m. var. striatum (three plants per pool), mapped to a pseudomajus assembly. Out of 45,067 predicted coding sequences, 3,109 (6.9 %) were differentially expressed between the pseudomajus and striatum pools (p-adjusted < 0.01) ( Fig. 5 ). Many of these differentially expressed loci may reflect the limited number of individuals pooled, rather than representing consistent differences across the geographic range of the varieties. BLASTN searches against the reference genome suggested that 75 differentially expressed genes mapped to monophyletic islands. 63 of these mapped to islands harbouring previously identified flower colour loci: 30 to Island 3 ( FLA ), 18 to Island 4 (including SULF ), two to Island 5 (including RUB ) and 13 to Island 6 (including ROS and EL ). Download figure Open in new tab Figure 5: Transcript comparisons between A. m. m. var. pseudomajus and A. m. m. var. striatum . DESeq2 differential expression results between two A. m. m. var. pseudomajus and three A. m. m. var. striatum pools, each containing three individuals, mapped to an A. m. m. var. pseudomajus genome assembly. Positive fold changes reflect elevated expression in A. m. m. var. striatum , and negative values in A. m. m. var. pseudomajus . All differentially expressed transcripts are shown in grey ( p < 0.01) with a subset of points coloured according to monophyletic island membership. Island 1 candidate genes from BLAST searches are annotated, along with the island 2 aureusidin synthase transcripts. AS1a / b = Aureusidin synthase, CGT = Chalcone glucosyltransferase, EL = ELUTA, FLA = FLAVIA, HIOMTa / b = Hydroxyindole-O-methyltransferase, ROMT = Trans-resveratrol di-O-methyltransferase, ROS1 = ROSEA 1, RUB = RUBIA, HCT = Shikimate O-hydroxycinnamoyltransferase, SCL8 -like = SCARECROW-like-8-like, QPT = Nicotinate-nucleotide pyrophosphorylase, UGT = UDP-glucuronosyltransferase. Of the remaining 12, 7 mapped to island 1 ( CRE ). NCBI BLASTX searches suggested that these included two hydroxyindole-O-methyltransferase genes ( HIOMTa and HIOMTb ), a trans-resveratrol di-O-methyltransferase gene ( ROMT ), a shikimate O-hydroxycinnamoyltransferase ( HCT ), a SCARECROW-like gene ( SCL8-like ), a nicotinate-nucleotide pyrophosphorylase (QPT), and a gene of unknown function. Activity of O-methyltransferases has been implicated in flower colour, through interactions with the anthocyanin biosynthesis pathway 31 , 32 , 33 . CRE may therefore encode an O-methyltransferase which affects yellow by altering flux of upstream aurone precursors or by methylation of aurones 34 . The SCARECROW-like gene is another possible candidate for CRE as SCARECROW-like transcription factors have been implicated in a range of developmental processes 35 . The remaining 5 differentially-expressed genes mapped to island 2 ( AUN ), two of which were the nearest homologues of the aureusidin synthase gene AmAS1 in the pseudomajus assembly . AmAS1 is responsible for synthesising yellow aurone pigments from chalcones 36 and is therefore a strong candidate for AUN . Origin of flower colour differences Our aphenotypic analysis shows that < 1 % of genomic regions give d XY trees that divide Antirrhinum populations into two distinct groups. The genomic regions map to 6 monophyletic islands, all of which harbour flower colour loci that interact to confer the phenotypic differences between A. m. m. var. striatum and A. m. m. var. pseudomajus. The loci exhibit clines in allele frequency at hybrid zones indicating they are barrier loci that reduce gene flow of linked regions between A. m. m. var. pseudomajus and A. m. m. var. striatum . Thus, the distinction between A. m. m. var striatum and A. m. m. var. pseudomajus is maintained by selection acting on a small group of barrier loci. These findings raise two key questions: why are there are so few barrier loci, and why do they all correspond to flower colour genes? To help address these questions, we first represent the phenotypes conferred by the flower colour loci with a genotype-phenotype space ( Fig. 6a ). The pseudomajus genotype is positioned at the top left, with alleles conferring spread magenta and restricted yellow, whereas the striatum genotype is positioned at the bottom right, with alleles conferring spread of yellow, and restricted magenta. Hybrid genotypes, such as those conferring white or orange flower colours, occupy other positions. Download figure Open in new tab Figure 6: Evolutionary scenarios for the origin of monophyletic islands in Antirrhinum . a Genotype-phenotype space for flower colour. Magenta with a restricted yellow highlight (top left) corresponds to A. m. m. var. pseudomajus ; whereas yellow with restricted magenta (bottom right) corresponds to A. m. m. var. striatum . b,c Genotype-fitness spaces for local adaptation hypothesis where one condition favours the A. m. m. var. pseudomajus genotype (b) and another favours the A. m. m. var. striatum genotype (c). Dark blue indicates lowest fitness, with fitness increasing up the colour gradient to a peak of red. d,e Genotype-fitness space viewed from the top (d) or side (e) where the phenotypes of A. m. m. var. pseudomajus and A. m. m. var. striatum occupy distinct adaptive peaks in a common environment. f Scenario for the origin of magenta and yellow flowered Antirrhinum species from a white-flowered ancestor, with subsequent secondary contact and gene flow (dark grey). The merging path following secondary contact reflects extensive gene flow, with the outer paths representing the fraction of the genome that resists gene flow. Percentages reflect approximate amount of the genome following each path. g 3D genotype-environment-fitness (GEF) space illustrating scenario in which A. m. m. var. pseudomajus and A. m. m. var. striatum descended from a white-flowered alpine ancestor. Vertical axis represents niche (growth habit/environment), and horizontal axes represent flower colour genotype. High-fitness regions are coloured red. One hypothesis for how flower colour genes act as barrier loci is local adaptation. In environmental conditions where A. m. m. var. pseudomajus grows, the top left position in genotype-phenotype space may have the highest fitness ( Fig. 6b ); whereas in conditions where A. m. m. var. striatum grows, the bottom right position may correspond to a fitness peak ( Fig. 6c ). Divergent selection on the flower colour loci would therefore prevent flow of alleles between the varieties. However, this hypothesis does not readily account for barrier loci operating at the hybrid zone, where there are no significant differences in Antirrhinum pollinators between the two sides 37 , or obvious abiotic environmental differences 38 . An alternative hypothesis is that the fitness landscape for a common environment has two adaptive peaks ( Fig. 6d and e ), and hybrid genotypes fall into a fitness valley. This hypothesis has two requirements. First, the phenotypes at each peak should represent alternative adaptive solutions to a common environmental challenge. Second, paths should exist that allowed populations to access one or the other peak on a microevolutionary timescale. Flower colour pattern likely satisfies the first requirement of providing alternative adaptive solutions because several colour patterns may act as effective pollinator guides. Thus, a yellow highlight on magenta, or a magenta highlight on yellow, could both be effective solutions for signalling pollinator entry. Each signal requires coadaptation of several components (the different genes controlling flower colour), allowing more than one adaptive peak to be produced through reciprocal sign epistasis 39 . Signalling in hybrid genotypes may be less effective because overlapping yellow and magenta reduces colour contrast (e.g. orange flower at top right of Fig. 6a ), or a white background may be less salient to pollinators (bottom left, Fig. 6a ). The second requirement – evolutionary accessibility of the peaks – may be satisfied because optimality of floral guides can vary with environment. Although flowers with a white background (bottom left, Fig. 6a ) are not prevalent in ruderal species (e.g. A. m. m. var. pseudomajus and A. m. m. var. striatum ), this phenotype is prevalent in alpine species, such as Antirrhinum molle . These alpine species have a trailing (procumbent) habit and grow in crevices of cliffs, in contrast to ruderal A. m. m. var. pseudomajus and A. m. m. var. striatum which have an erect habit and grow on slopes. A white background might be more salient to pollinators when contrasted against the solid rock of a cliff-face, but less salient than magenta or yellow in a more open or ruderal setting. A plausible scenario is that a common ancestor of A. m. m. var. pseudomajus and A. m. m. var. striatum was a cliff-dwelling alpine, with a white flower background ( Fig. 6f ). In colonising slope environments, descendants evolved an erect habit, leading to occupancy of a new niche (ruderal) and a corresponding change in the fitness landscape from a single peak (white background) to two peaks (magenta or yellow background). Selection could have driven some geographically isolated populations along a fitness path leading to a magenta background ( Fig. 6f , left), and others along a path leading to a yellow background ( Fig. 6f , right). In the absence of gene flow, genomes in different populations would undergo divergence (increase in d XY ). Following secondary contact (dark grey, Fig. 6f ), gene flow would erode genome divergence except at or near flower colour loci, leading to deep-rooted d XY trees for <1 % of the genome. Some of the barrier regions would exhibit F ST peaks, particularly those that underwent strong selective sweeps 24 . Connectivity of adaptive peaks in this scenario can be summarised with a three-dimensional genotype-environment-fitness (GEF) space, in which the vertical axis represents growth habit and environment (collectively defining an ecological niche), and the horizontal axes genotype at flower colour loci ( Fig. 6g ). Regions of high fitness are colour-coded red. The two adaptive peaks in A. m. m. var. pseudomajus and A. m. m. var. striatum correspond to an upper slice through this GEF space. They are connected through a higher dimensional path of high fitness, or evolutionary wormhole 40 , that depends on genes controlling growth habit (vertical axis) and genes controlling flower colour (horizontal axes). Such bifurcating fitness paths leading to different positions in a common environment may be rarely followed on a microevolutionary timescale. Flower colour patterns, which can provide alternative coadapted solutions that vary in efficacy according to environment, may provide one such rare example. Thus, the paucity of barrier loci and their correspondence to flower colour genes may be accounted for by the rarity of bifurcating fitness paths. The situation illustrated by A. m. m. var. pseudomajus and A. m. m. var. striatum , with a few barrier loci all relating to one trait, may represent one end of a speciation continuum 41 . In other cases, which involve adaptation to different environments, flower colour is only one of several traits that distinguish morphs and many more genomic regions of reduced gene flow are likely involved 42 , 43 , 44 , 9 , 45 , 46 . A small number of barrier loci has been described for wing colour patterns in butterflies, which can act as warning (aposematic) signals for predators such as birds 47 , 48 . As with floral guides, this represents a situation in which alternative coadapted solutions to a common problem (avoidance of predation in this case). The method we describe for identifying candidate barrier loci through shared deep-rooted trees complements approaches based on cline and F ST analysis. Employing d XY , rather than other measures, has the advantage of identifying candidate loci without confounding effects of inbreeding or selective sweeps. Our approach also has the advantage of being aphenotypic, with the potential to reveal barrier loci controlling hidden traits. The results indicate that frequent involvement of colour patterns in barrier loci is not solely a consequence of their visual saliency to geneticists and ecologists. Rather, it may arise because colour patterns lend themselves to the formation of multiple peaks in a fitness landscape, the heights, and positions of which vary depending on genotype-environment interactions, allowing different peaks to be accessed over a microevolutionary timescale. The reproductive isolation conferred by these peaks forms a speciation continuum, providing a model for how reproductive barriers between species can arise and be maintained through populations taking diverse paths in gene-environment-fitness spaces. Methods DNA extraction from pooled leaf tissue for whole genome sequencing Genomic DNA was isolated using a cetyltrimethylammonium bromide (CTAB) method on 2-5g of leaves harvested as final weight for pooled samples as described 49 . Field samples were either placed in bags of silica or water-moist paper towel and kept cool at 4°C until they were courier posted by overnight delivery to the John Innes Centre and, upon arrival, frozen at - 80°C. Greenhouse samples were stored in foil at −80°C. A GPS reading was taken at each wild population sampling location (recorded in Supplementary Table 1). Short read sequencing (2 x 150bp) of DNA extracted from leaf tissue was carried out as described 50 . Long read PacBio sequencing of A. m. m. var. pseudomajus sample Ac1099 was also carried out as described 50 . RNA extraction from petals Total RNA was isolated from petal tissues using the Qiagen RNeasy Plant Mini Kit, including DNaseI treatment. For a comparison of A. m. m. var. pseudomajus (accession Ac1266) and A. m. m. var. striatum (accession Ac1125) we harvested petal lobes from dissected flowers just before opening. Three independent samples were used for reproducibility, and each sample was a pool of 3 individuals, each contributing 1-2 flowers. Samples were sequenced by Novogene using an Illumina HiSeq 2500. Growth conditions and crosses Greenhouse cultivated accessions were grown in JI compost soil mixes as described 51 , with supplemental lights in winter giving 12 to 16 hour days. Plants were grown outside in the summer using the same soil mixes, in pots or plugs in trays placed upon raised benches. Plants stocks were maintained through inter-sibling crosses. To carry out crosses, pre-anthesis floral buds were opened with forceps and young anthers were removed. Once the flower was open (within two to four days) pollen from another plant was applied to the stigma, using the forceps to carry pollen or the whole stamen. Seed capsules were harvested after four to six weeks. Flower photography Flowers were placed on a black velvet background with a scale bar and a Small Grey & Colour Separation Chart (Danes Picta BST13) for monitoring light level, colour and white balance. We used an Olympus XZ-1 (10 Megapixels) Camera, with side-on lighting from table lamps fitted with halogen 42W 630 lumen (2800k) warm white light bulbs. Camera settings were set to the closest White Balance of 3000K, with no flash, an Exposure Time of 1/20-1/40 sec, F stop 8.0, ISO 200, RAW images, Macro On, aspect 4:3, high definition. KASP genotyping KASP Genotyping was performed as described 25 . Fluorescent signals were detected using a BioRad CFX96 light cycler, and data processed using BioRad CFX Manager v3.1. AFLP methods used standard PCR and agarose gel analysis as described 26 . KASP / AFLP oligos are as described 26 . Assembly of A. m. m. var. pseudomajus genome A. m. m. var. pseudomajus accession Ac1099 was grown from seed, and leaf tissue was harvested and sequenced. Raw Illumina reads were trimmed using Trim Galore! with parameters -q 25 and --stringency 3 to remove low-quality sequence, and sequence corresponding to adapters. To infer genome properties prior to assembly, GenomeScope2 52 was used to evaluate raw reads based on the k-mer spectrum, heterozygosity rate and haploid genome length. Contig-level assemblies of PacBio reads were carried out using the Canu package (version 1.9) 53 . Duplicate contigs were removed using purge_dups 54 . RepeatMasker 55 was used to identify and mask repetitive sequence. High-quality plant proteins were retrieved from SwissProt 56 and aligned to the genome using ProtHint 57 . RNAseq data were also aligned to the assembly, and gene prediction was carried out using Braker2 58 . The overall quality of the assembly, and gene annotations, was assessed using LAI 59 and BUSCO 60 . Mapping reads to reference genome and generating SYNC files FASTQ reads for each of the sequenced pools were mapped to the A. majus reference genome V3.0 (Genome Warehouse accession number GWHBJVT00000000) and to the A. m. m. var. pseudomajus genome using BWA-MEM, with the -M flag set for downstream compatibility with Picard ( http://broadinstitute.github.io/picard/ ). Output SAM files were sorted using SAMtools, and Picard MarkDuplicates was used to remove duplicate reads. Local realignment of reads around indels was carried out using GATK version 3.5.0 61 . GATK RealignerTargetCreator was used to generate a list of intervals for realignment using GATK IndelRealigner. Processed BAM files were combined to generate MPILEUP files using SAMtools mpileup 62 . The minimum mapping quality threshold ( -q ) was set to 40, with a minimum base quality threshold ( -Q ) of 30. Orphaned reads were included in variant calling by using the -A flag. The -B flag was used to disable probabilistic realignment in base alignment quality calculations, as this can increase misalignments and false SNP calls. One MPILEUP file was generate for each of the eight A. majus chromosomes. MPILEUP files were converted to the Popoolation2 SYNC format using mpileup2sync.jar, which offers a more concise representation of allelic depth across populations 63 . Population genomic analyses Within-population diversity, π w , was estimated as π w = p 1 q 1 , where p 1 and q 1 refer to the frequencies of alleles p and q within a single population, population 1. d XY , between-population diversity, was estimated using , where 1 and 2 refer to a pair of populations, 1 and 2. F ST , relative divergence, is estimated from d XY and πw as , where π t (total diversity) is the sum of d XY and π w . Nei’s standard genetic distance, D , was calculated as d XY − π w . treeXY population analyses Population genomic statistics were calculated using the treeXY software version 1.1. Statistics were summarised in windows – unless otherwise stated, analyses detailed here used a 50 kb window size, with 25 kb overlaps between adjacent windows. treeXY was run on all SYNC files, yielding windows for all eight chromosomes. Window size was set to 50 kb using -w , and window overlap to 25 kb using -o . To output treeXY-filtered SYNC files, -- write_sync was set, and --compute_trees was used to enable computation of SNP trees. --ignore_multiallelic was set such that multiallelic sites would be skipped. SNPs were required to be present in all 18 samples ( -A 18 ). Otherwise, default parameters were used (minimum depth = 15, maximum depth = 200, minimum allele depth = 2). Constructing mean d XY and D trees Distance matrices were populated with between-population statistics calculated for each genomic window. UPGMA hierarchical clustering trees were generated for each distance matrix using the base R hclust function, with clustering method set to average. To generate mean trees across forests, the mean was taken across all windows using the Reduce function. Whole genome sweep analysis A whole genome selective sweep was simulated within the one population, MP11, by fixing all polymorphic sites for one allele, randomly sampled based on the existing allele frequencies. To do this, treeXY was run with default settings, with the --write_sync flag enabled to output a depth-filtered SYNC file. A custom Python script, artificial_sweep.py , was run with the filtered SYNC file as input. At all genomic sites where more than one allele showed depth >= 2 (the default allelic depth threshold,) one allele was randomly sampled. Frequencies of ambiguous (N) alleles were ignored. The probability of sampling an allele was weighted according to its frequency prior to sweeping, with common alleles more likely to be sampled than rare ones. After sweeping, the depth of the chosen allele was changed to equal the sum of the total site depth, and all other allele frequencies were set to 0. The modified SYNC file, with the swept MP11 population, was written to a new file for downstream analyses. Whole genome trees with d XY , D, and maximum likelihood For d XY and D trees, a modified version of the treeXY script was run, which calculated the mean for both genetic distance measures across all biallelic sites and between all pairs of populations. Results from all eight SYNC files were averaged, and used to populate distance matrices. These matrices were then summarised as UPGMA trees using the base R hclust function. For maximum likelihood (ML) trees, filtered SYNC files obtained using the treeXY --write_sync feature were concatenated into a whole genome biallelic SYNC file. At each site, and in each population, a consensus base call was made by selecting the allele with the highest frequency. This yielded a set of aligned consensus sequences, which was converted into FASTA format. To generate an ML tree from this alignment, RAxML-NG 64 version 0.9.0 was used, with the arguments --all --model GTR+G --tree pars{10} --bs-trees 100 . ML trees were drawn in RStudio, using the phangorn library. d XY , D, and ML trees were generated before and after applying the whole genome selective sweep within the MP11 population Calculating shortest root branch and cophenetic correlation coefficient The shortest root branch (SRB) of a given tree was calculated from its cophenetic matrix, obtained using the base R cophenetic function. SRB is equal to the maximum value within the cophenetic matrix, subtract the second highest value. The similarity of hierarchical clustering trees was estimated using the cophenetic correlation coefficient, c 29 . To compare two trees, the Pearson correlation coefficient was calculated between their cophenetic matrices, where: x and y correspond to the distances in each matrix. When carrying out the grouping tree scan to populate forests, a c threshold of 0.5 was chosen to capture trees showing moderate or high topological similarity. To test the replicability of grouping tree scan results, a bootstrapping approach was implemented. Each bootstrap replicate consisted of one grouping tree scan. For each scan, all forests were summarised based on the mean SRB of their trees. The forest showing the highest mean SRB was deemed the outlier forest, and recorded. Once all bootstrap replicates had been completed, all outlier forest regions were recorded, along with the frequency with which they had been detected. Classifying trees based on root division Root division was determined using the base R cutree function, with k = 2, which bisects trees at their topmost branch to yield the two outermost clades. To classify trees based on root division, populations within each clade were recorded, and compared to a user specified signature. Geographic cline analysis Genotyping of plants across the Planoles hybrid zone, and clinal analysis, was carried out as described 27 . Colour ranking of flower photographs Flower photographs processed in Photoshop. Flower photographs included a Danes Picta BST13 colour chart, which was selected with the eyedropper tool to standardise white balance. Adjusted images were saved in the JPEG format. Images were cropped to retain only the flowers, before being loaded into a Photoshop canvas. Each image was, in turn, visually appraised for the intensity of its yellow colour, and positioned within the rank accordingly. Positions were adjusted as more photographs were incorporated, until a final ranking was reached. Upon completion, the ranking was split into either two halves, or four quartiles. Within each group, the CRE and AUN genotypes were checked, and the number of A. m. m. var. pseudomajus and A. m. m. var. striatum alleles was inferred. RNAseq differential expression analysis using DESeq2 Ribo-depleted total RNA libraries from A. m. m. var. pseudomajus and A. m. m. var. striatum were mapped to the A. m. m. var. pseudomajus assembly using HISAT2 65 . Output SAM files were sorted using Samtools, and transcripts were assembled using StringTie 66 , operating in expression estimation mode ( -e ) and reporting gene abundances ( -A ). To prepare data tables for analysis, read counts were extracted across all samples using the StringTie prepDE.py Python script. Differential expression analysis was carried out using DESeq2 version 1.32.0 30 using default parameters. Data availability Raw DNA and RNA data have been uploaded to SRA under accession number SUB15081578. The A. m. m. var. pseudomajus assembly and GFF annotations have been uploaded to NCBI WGS under accession number SUB15081867. The A. majus reference genome V3.0 is available at the NGDC Genome Warehouse under accession number GWHBJVT00000000. Code availability The treeXY software and documentation are available at https://github.com/DR-Antirrhinum/treeXY . Other scripts used to carry out the analyses presented here are available at https://github.com/DR-Antirrhinum/phylogenetic_forests . Author contributions DR designed and carried out bioinformatic analyses of DNA and RNA, developed research software, managed data, and redrafted the manuscript with EC. DB carried out genotyping, DNA / RNA extraction for sequencing, and flower photography, and assisted with flower colour ranking and manuscript discussion. LC designed and carried out crosses to provide experimental plant populations. AW designed RNA and DNA mapping pipelines, and contributed to manuscript discussions. MB and CA together undertook fieldwork to identify and collect tissue from wild plant populations. SZ carried out the assembly of the A. m. m. var. pseudomajus genome, helped with data transfer and organisation. DF carried out cline analysis of the CRE locus, and contributed to project supervision and manuscript discussions. YX administered and supervised genome sequencing, assembly and data transfer. EC conceived and administered the project, wrote and redrafted the manuscript with DR, designed experiments and contributed to methodology. Competing interest declaration The authors declare no competing interests. Additional information Supplementary Information is available for this paper. Correspondence and requests for materials should be addressed to Enrico Coen, Yongbiao Xue, or David Field Acknowledgements We thank the John Innes Centre (JIC) Horticultural Services for providing growth facilities and maintenance of plant material, and JIC Research Computing for use of High Performance Computing facilities. We also thank Tingting Li for carrying out a replicate flower colour ranking analysis and assisting with photography. This work was supported by the Biotechnology and Biological Sciences Research Council (grants BB/S009256/1, BB/G009325/1, BBS/E/JI/230002C, and BBS/E/J/000PR9773 to EC, and Norwich Research Park Biosciences Doctoral Training Partnership grant BB/M011216/1 to DR), the Natural Science Foundation of China (grant 32030007, to YX), and an ANR funded French Laboratory of Excellence project (‘LabEx TULIP’, to MB). This research was also funded in whole or in part by the Austrian Science Fund (FWF) [P 32166] (to DF). References 1. ↵ Wright , S . The roles of mutation, inbreeding, crossbreeding and selection in evolution . Proc. 6th int. Cong. Genet . 1 , 356 – 366 ( 1932 ). OpenUrl 2. ↵ Coyne , J. A. , Barton , N. H. & Turelli , M . Perspective: a critique of Sewall Wright’s shifting balance theory of evolution . Evolution 51 , 643 – 671 ( 1997 ). OpenUrl CrossRef PubMed Web of Science 3. ↵ Steinberg B. , Ostermeier M . Environmental changes bridge evolutionary valleys . Sci. Adv . 2 , e1500921 ( 2016 ). OpenUrl FREE Full Text 4. ↵ Gavrilets , S . Fitness Landscapes and the Origin of Species ( Princeton University Press , Princeton , 2004 ). 5. ↵ Gerald , N. , Dutta , D. , Brajesh , R. , G. & Saini , S. Mathematical modeling of movement on fitness landscapes . BMC Syst. Biol . 13 , 25 ( 2019 ). OpenUrl CrossRef PubMed 6. ↵ Vigué , L. et al. Deciphering polymorphism in 61,157 Escherichia coli genomes via epistatic sequence landscapes . Nat. Commun . 13 , 4030 ( 2022 ). OpenUrl CrossRef PubMed 7. ↵ Christmas , M. J. et al. Genetic Barriers to Historical Gene Flow between Cryptic Species of Alpine Bumblebees Revealed by Comparative Population Genomics . Mol. Biol. Evol . 38 , 3126 – 3143 ( 2021 ). OpenUrl CrossRef PubMed 8. ↵ Enbody , E. D. et al. Community-wide genome sequencing reveals 30 years of Darwin’s finch evolution . Science 381 , eadf6218 ( 2023 ). OpenUrl CrossRef PubMed 9. ↵ Nelson , T. C. et al. Ancient and recent introgression shape the evolutionary history of pollinator adaptation and speciation in a model monkeyflower radiation ( Mimulus section Erythranthe ) . PLOS Genet . 17 , e1009095 ( 2021 ). OpenUrl CrossRef PubMed 10. ↵ Carneiro , M. et al. Steep clines within a highly permeable genome across a hybrid zone between two subspecies of the European rabbit . Mol. Ecol . 22 , 2511 – 2525 ( 2013 ). OpenUrl CrossRef Web of Science 11. ↵ Haenel , Q. et al. Clinal genomic analysis reveals strong reproductive isolation across a steep habitat transition in stickleback fish . Nat. Commun . 12 , 4850 ( 2021 ). OpenUrl CrossRef PubMed 12. ↵ Ravinet , M. et al. Interpreting the genomic landscape of speciation: a road map for finding barriers to gene flow . J. Evol. Biol . 30 , 1450 – 1477 ( 2017 ). OpenUrl CrossRef PubMed 13. ↵ Cruickshank , T. E. & Hahn , M. W . Reanalysis suggests that genomic islands of speciation are due to reduced diversity, not reduced gene flow . Mol. Ecol . 23 , 3133 – 3157 ( 2014 ). OpenUrl CrossRef PubMed Web of Science 14. ↵ Martin , S. H. & Van Belleghem , S. M . Exploring Evolutionary Relationships Across the Genome Using Topology Weighting . Genetics 206 , 429 – 438 ( 2017 ). OpenUrl Abstract / FREE Full Text 15. ↵ Fontaine , M. C. et al. Extensive introgression in a malaria vector species complex revealed by phylogenomics . Science 347 , 1258524 ( 2015 ). OpenUrl Abstract / FREE Full Text 16. ↵ Pease , J. B. , Brown , J. W. , Walker , J. F. , Hinchliff , C. E. & Smith , C. A . Quartet Sampling distinguishes lack of support from conflicting support in the green plant tree of life . Am J. Bot . 105 , 385 – 403 ( 2018 ). OpenUrl CrossRef PubMed 17. ↵ Rose , J. P. , Kriebel , R. , Sytsma , K. J. & Drew , B. T . Phylogenomic perspectives on speciation and reproductive isolation in a North American biodiversity hotspot: an example using California sages ( Salvia subgenus Audibertia : Lamiaceae) . Ann. Bot . 134 , 295 – 310 ( 2024 ). OpenUrl CrossRef PubMed 18. ↵ Rosser N. et al. Hybrid speciation driven by multilocus introgression of ecological traits . Nature 628 , 811 – 817 ( 2024 ). OpenUrl CrossRef PubMed 19. ↵ Nei , M. & Li , W. H . Mathematical model for studying genetic variation in terms of restriction endonucleases . Proc. Natl. Acad. Sci. U. S. A 76 , 5269 – 5273 ( 1979 ). OpenUrl Abstract / FREE Full Text 20. ↵ Vargas , P. , Carrió , E. , Guzmán , B. , Amat , E. & Güemes , J . A geographical pattern of Antirrhinum (Scrophulariaceae) speciation since the Pliocene based on plastid and nuclear DNA polymorphisms . J. Biogeogr . 36 , 1297 – 1312 ( 2009 ). OpenUrl CrossRef 21. ↵ Durán-Castillo , M. , Hudson , A. , Wilson , Y. , Field , D. L. & Twyford , A. D . A phylogeny of Antirrhinum reveals parallel evolution of alpine morphology . New Phytol . 233 , 1426 – 1439 ( 2022 ). OpenUrl CrossRef PubMed 22. ↵ Rothmaler , W . Taxonomische Monographie der Gattung Antirrhinum ( Akademie-Verlag , Berlin , 1956 ). 23. ↵ Sutton , D. A . A Revision of the Tribe Antirrhineae ( Oxford Univ. Press , Oxford , 1988 ). 24. ↵ Tavares , H. et al. Selection and gene flow shape genomic islands that control floral guides . Proc. Natl. Acad. Sci. U. S. A . 115 , 11006 – 11011 ( 2018 ). OpenUrl Abstract / FREE Full Text 25. ↵ Bradley , D. et al. Evolution of flower color pattern through selection on regulatory small RNAs . Science 358 , 925 – 928 ( 2017 ). OpenUrl Abstract / FREE Full Text 26. ↵ Bradley , D. et al. Shaping of developmental gradients through selection on gene interactions in Antirrhinum. 27. ↵ Field , D. L. et al. Genome-wide cline analysis identifies new locus contributing to a barrier to gene flow across an Antirrhinum hybrid zone . 28. ↵ Li , M. et al. Genome structure and evolution of Antirrhinum majus L . Nat. Plants 5 , 174 – 183 ( 2019 ). OpenUrl CrossRef PubMed 29. ↵ Sokal , R. R. & Rohlf , F. J . The Comparison of Dendrograms By Objective Methods . Taxon 11 , 33 – 40 ( 1962 ). OpenUrl CrossRef 30. ↵ Love , M. I. , Huber , W. & Anders , S . Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2 . Genome Biol . 15 , 550 ( 2014 ). OpenUrl CrossRef PubMed 31. ↵ Akita , Y. et al. Isolation and characterization of the fragrant cyclamen O-methyltransferase involved in flower coloration . Planta 234 , 1127 – 1136 ( 2011 ). OpenUrl CrossRef PubMed Web of Science 32. ↵ Du , H. et al. Methylation mediated by an anthocyanin, O-methyltransferase, is involved in purple flower coloration in Paeonia . J. Exp. Bot . 66 , 6563 – 6577 ( 2015 ). OpenUrl CrossRef PubMed 33. ↵ Okitsu , N. , Mizuno , T. , Matsui , K. , Choi , S. H. & Tanaka , Y . Molecular cloning of flavonoid biosynthetic genes and biochemical characterization of anthocyanin o-methyltransferase of Nemophila menziesii Hook. and Arn . Plant Biotechnol . 35 , 9 – 16 ( 2018 ). OpenUrl CrossRef 34. ↵ Mazziotti , I. , Petrarolo , G. & La Motta , C . Aurones: Golden resource for active compounds . Molecules 27 , 2 ( 2022 ). 35. ↵ Cenci , A. & Rouard , M . Evolutionary Analyses of GRAS Transcription Factors in Angiosperms . Front. Plant. Sci . 8 , 273 ( 2017 ). 36. ↵ Nakayama , T. et al. Aureusidin Synthase: A Polyphenol Oxidase Homolog Responsible for Flower Coloration . Science 290 , 1163 – 1166 ( 2000 ). OpenUrl Abstract / FREE Full Text 37. ↵ Andalo , C. , Burrus , M. , Sandrine , B. , Lauzeral , C. & Field , D. L . Prevalence of legitimate pollinators and nectar robbers and the consequences for fruit set in an Antirrhinum majus hybrid zone . Bot. Lett . 166 , 1 – 13 ( 2018 ). OpenUrl 38. ↵ Khimoun , A. et al. Ecology predicts parapatric distributions in two closely related Antirrhinum majus subspecies . Evol. Ecol . 27 , 51 – 64 ( 2012 ). OpenUrl 39. ↵ Poelwijk , F. J. , Tǎnase-Nicola , S. , Kiviet , D. J. & Tans , S. J. Reciprocal sign epistasis is a necessary condition for multi-peaked fitness landscapes . J. Theor. Biol . 272 , 141 – 144 ( 2011 ). OpenUrl CrossRef PubMed Web of Science 40. ↵ Prusinkiewicz , P. , Erasmus , Y. , Lane , B. , Harder , L. D. & Coen , E . Evolution and development of inflorescence architectures . Science 316 , 1452 – 1455 ( 2007 ). OpenUrl Abstract / FREE Full Text 41. ↵ Mallet , J. & Dasmahapatra , K. K . Hybrid zones and the speciation continuum in Heliconius butterflies . Mol. Ecol . 21 , 5643 – 5645 ( 2012 ). OpenUrl CrossRef Web of Science 42. ↵ Stankowski , S. , Sobel , J. M. , Streisfeld , M. A . The geography of divergence with gene flow facilitates multitrait adaptation and the evolution of pollinator isolation in Mimulus aurantiacus . Evolution 69 , 3054 – 3068 ( 2015 ). OpenUrl CrossRef PubMed 43. ↵ Bradshaw , H. D. & Schemske , D. W . Allele substitution at a flower colour locus produces a pollinator shift in monkeyflowers . Nature 426 , 176 – 178 ( 2003 ). OpenUrl CrossRef PubMed Web of Science 44. ↵ Liang , M. et al. Taxon-specific, phased siRNAs underlie a speciation locus in monkeyflowers . Science 379 , 576 – 582 ( 2023 ). OpenUrl CrossRef PubMed 45. ↵ Wessinger , C. A. et al. A few essential genetic loci distinguish Penstemon species with flowers adapted to pollination by bees or hummingbirds . PLOS Biol . 21 , e3002294 ( 2023 ). OpenUrl CrossRef PubMed 46. ↵ Stankowski , S. , Chase , M. A. , McIntosh , H. & Streisfeld , M. A . Integrating top-down and bottom-up approaches to understand the genetic architecture of speciation across a monkeyflower hybrid zone . Mol. Ecol . 32 , 2041 – 2054 ( 2023 ). OpenUrl CrossRef 47. ↵ Seymoure , B. M. , Raymundo , A. , McGraw , K. J. , McMillan , W. O. & Rutowski , R. L . Environment-dependent attack rates of cryptic and aposematic butterflies . Curr. Zool . 64 , 663 – 669 ( 2018 ). OpenUrl CrossRef PubMed 48. ↵ Nadeau , N. J. et al. Population genomics of parallel hybrid zones in the mimetic butterflies, H. melpomene and H. erato . Genome Res . 24 , 1316 – 1333 ( 2014 ). OpenUrl Abstract / FREE Full Text Methods references 49. ↵ Coen , E. S. , Carpenter , R. & Martin , C . Transposable elements generate novel spatial patterns of gene expression in Antirrhinum majus . Cell 47 , 285 – 296 ( 1986 ). OpenUrl CrossRef PubMed Web of Science 50. ↵ Zhu , S. et al. The Snapdragon Genomes Reveal the Evolutionary Dynamics of the S-Locus Supergene . Mol. Biol. Evol . 40 , msad080 ( 2023 ). OpenUrl CrossRef PubMed 51. ↵ Martin , C. et al. The control of floral pigmentation in Antirrhinum . Biochem. Soc. Trans . 15 , 14 – 17 ( 1987 ). OpenUrl FREE Full Text 52. ↵ Ranallo-Benavidez , T. R. , Jaron , K. S. & Schatz , M. C . GenomeScope 2.0 and Smudgeplot for reference-free profiling of polyploid genomes . Nat. Commun . 11 , 1432 ( 2020 ). OpenUrl CrossRef PubMed 53. ↵ 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 ). OpenUrl Abstract / FREE Full Text 54. ↵ Guan , D. et al. Identifying and removing haplotypic duplication in primary genome assemblies . Bioinformatics 36 , 2896 – 2898 ( 2020 ). OpenUrl CrossRef PubMed 55. ↵ Smit , A. F. A. , Hubley , R. & Green , P. RepeatMasker Open-4.0 . ( 2013 ). 56. ↵ The UniProt Consortium UniProt: the Universal Protein Knowledgebase in 2023 , Nucleic Acids Res . 51 , D523 – D531 ( 2023 ). OpenUrl CrossRef PubMed 57. ↵ Brůna , T. , Lomsadze , A. & Borodovsky , M . GeneMark-EP+: Eukaryotic gene prediction with self-training in the space of genes and proteins . NAR Genom. Bioinform . 2 , 1 – 14 ( 2020 ). OpenUrl 58. ↵ Brůna , T. , Hoff , K. J. , Lomsadze , A. , Stanke , M. & Borodovsky , M . BRAKER2: Automatic eukaryotic genome annotation with GeneMark-EP+ and AUGUSTUS supported by a protein database . NAR Genom. Bioinform . 3 , 1 – 11 ( 2021 ). OpenUrl 59. ↵ Ou , S. , Chen , J. & Jiang , N . Assessing genome assembly quality using the LTR Assembly Index (LAI) . Nucleic Acids Res . 46 , e126 ( 2018 ). OpenUrl PubMed 60. ↵ Simão , F. A. , Waterhouse , R. M. , Ioannidis , P. , Kriventseva , E. V. & Zdobnov , E. M . BUSCO: Assessing genome assembly and annotation completeness with single-copy orthologs . Bioinformatics 31 , 3210 – 3212 ( 2015 ). OpenUrl CrossRef PubMed 61. ↵ McKenna , A. et al. The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data . Genome Res . 20 , 1297 – 303 ( 2010 ). OpenUrl Abstract / FREE Full Text 62. ↵ Danecek , P. et al. Twelve years of SAMtools and BCFtools . Gigascience 10 , giab008 ( 2021 ). OpenUrl CrossRef PubMed 63. ↵ Kofler , R. , Pandey , R. V. & Schlötterer , C . PoPoolation2: Identifying differentiation between populations using sequencing of pooled DNA samples (Pool-Seq) . Bioinformatics 27 , 3435 – 3436 ( 2011 ). OpenUrl CrossRef PubMed Web of Science 64. ↵ Kozlov , A. M. , Darriba , D. , Flouri , T. , Morel , B. & Stamatakis , A . RAxML-NG: A fast, scalable and user-friendly tool for maximum likelihood phylogenetic inference . Bioinformatics 35 , 4453 – 4455 ( 2019 ). OpenUrl CrossRef PubMed 65. ↵ Kim , D. , Paggi , J. M. , Park , C. , Bennett , C. & Salzberg , S. L . Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype . Nat. Biotechnol . 37 , 907 – 915 ( 2019 ). OpenUrl CrossRef PubMed 66. ↵ Shumate , A. , Wong , B. , Pertea , G. & Pertea , M . Improved transcriptome assembly using a hybrid of long and short reads with StringTie . PLoS Comput. Biol . 18 , 1 – 18 ( 2022 ). OpenUrl CrossRef PubMed View the discussion thread. Back to top Previous Next Posted February 16, 2025. Download PDF Supplementary Material Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Genomic tree scans identify loci underlying adaptive peaks in Antirrhinum Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Genomic tree scans identify loci underlying adaptive peaks in Antirrhinum Daniel M. Richardson , Desmond Bradley , Lucy Copsey , Annabel Whibley , Monique Burrus , Christophe Andalo , Sihui Zhu , David L. Field , Yongbiao Xue , Enrico Coen bioRxiv 2025.02.12.637406; doi: https://doi.org/10.1101/2025.02.12.637406 Share This Article: Copy Citation Tools Genomic tree scans identify loci underlying adaptive peaks in Antirrhinum Daniel M. Richardson , Desmond Bradley , Lucy Copsey , Annabel Whibley , Monique Burrus , Christophe Andalo , Sihui Zhu , David L. Field , Yongbiao Xue , Enrico Coen bioRxiv 2025.02.12.637406; doi: https://doi.org/10.1101/2025.02.12.637406 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Evolutionary Biology Subject Areas All Articles Animal Behavior and Cognition (7635) Biochemistry (17691) Bioengineering (13892) Bioinformatics (41937) Biophysics (21452) Cancer Biology (18588) Cell Biology (25504) Clinical Trials (138) Developmental Biology (13378) Ecology (19899) Epidemiology (2067) Evolutionary Biology (24320) Genetics (15609) Genomics (22506) Immunology (17736) Microbiology (40394) Molecular Biology (17181) Neuroscience (88605) Paleontology (666) Pathology (2832) Pharmacology and Toxicology (4824) Physiology (7641) Plant Biology (15156) Scientific Communication and Education (2045) Synthetic Biology (4294) Systems Biology (9825) Zoology (2271)