Three chromosome-scale Papaver genomes reveal punctuated patchwork evolution of the morphinan and noscapine biosynthesis pathway

For millions of years, plants evolve plenty of structurally diverse secondary metabolites (SM) to support their sessile lifestyles through continuous biochemical pathway innovation. While new genes commonly drive the evolution of plant SM pathway, how a full biosynthetic pathway evolves remains poorly understood. The evolution of pathway involves recruiting new genes along the reaction cascade forwardly, backwardly, or in a patchwork manner. With three chromosome-scale Papaver genome assemblies, we here reveal whole-genome duplications (WGDs) apparently accelerate chromosomal rearrangements with a nonrandom distribution towards SM optimization. A burst of structural variants involving fusions, translocations and duplications within 7.7 million years have assembled nine genes into the benzylisoquinoline alkaloids gene cluster, following a punctuated patchwork model. Biosynthetic gene copies and their total expression matter to morphinan production. Our results demonstrate how new genes have been recruited from a WGD-induced repertoire of unregulated enzymes with promiscuous reactivities to innovate efficient metabolic pathways with spatiotemporal constraint.

The authors have presented a comparative genomic study of three species of poppy using new genome assemblies to study how to genes involved in alkaloid metabolism have evolved over time. This is quite an interesting area of research and it's good to see the new genome assemblies. Figure 4 is the most important contribution of this paper and on this piece of evidence alone, this paper merits eventual publication --but this figure has some problems and could be substantially clarified (see comment #5). I think greater care should be taken to not overstate the significance of any results and to take a balanced view of what they mean. It is more interesting to just show us the data and allow the data to speak for themselves, rather than trying to build elaborate stories about how something might have happened when the evidence isn't really conclusive supporting a given sequence of rearrangement/duplication/fusion events (e.g. comment #2). I find the authors are often seeming to push a little too hard on the importance or strength of their conclusions, perhaps as they are worried about being perceived as sufficiently novel/impactful for Nature Communications? This is detrimental to the paper. Also, the question of cluster evolution has been considered extensively by other papers and if this is the aim of this paper, more connection should be made to this literature (see comments 7 & 8). A major question in the evolution of gene clusters is how commonly they have evolved from tandem duplication vs. rearrangement. It is critical that this aspect be very rigourously assessed. I think with some substantial re-working this manuscript could be suitable for publication, but it would really need to push less on promoting one particular story and take a much more balanced view of the data. I would appreciate if the authors approached their data with a bit more hesitancy and acknowledgement of the areas of uncertainty --this is fine! Having two new genomes provides lots of scope for assessing how rearrangements happened and is sufficiently novelty to merit publication in Nature Communications, but only if this is discussed in a more balanced and careful way.
Major comments: 1. I found it very difficult to understand the results of the paper without going to the supplementary materials. I think this paper is much more suitable for a longer and more in-depth format in a specialty journal. Most of the analyses are descriptive with little in the way of hypothesis testing --the reconstructions of orthology relationships don't conclusively support one model of gene evolution over another, as this would require sequencing of more species.
2. The previous paper by Guo et al (Science 2018) focused extensively on evolution of STORR. It is unclear what is contributed by this further analysis --can the authors more clearly describe what this shows that was not known previously? I found the paragraphs and figures on STORR particularly confusing and unclear. The model of "Fusion-Translocation" seems like one of many possible sequences of events and this seems more like evolutionary story-telling than hypothesis testing. In fact, given the figures as they are presented, it seems less likely than an alternative hypothesis. As I understand it, here is the presented FT-hypothesis: -two parent genes leading to STORR were present on same chromosome prior to WGD -WGD occurred in ancestor of setigerum and somniferum (not rhoeas) -some "Fusion, translocation (FT)" event happened (which is not described but says it may involve some translocation event) The authors claim the evidence in Figure 3 and Figure S17 supports this claim. But I fail to see how this is more parsimonious than simply having some deletion between the intergenic regions of the "pre-fusion" modules? This would be SO much clearer if the paralogs to the pre-fusion modules within somniferum were shown on Figure 4. Where are they? I can only assume they must be on some other chromosome, which would actually be much clearer support for the "translocation" part of the story than the convoluted Figure 3 and accompanying explanation. Alternatively, expanding Figure 3 to show the chromosomal context more clearly might be useful. It's clear from Figure 4 that there are translocations involved becuase the noscapine cluster got fused with the morphinan cluster. I generally feel like this manuscript spends a lot of words trying to prove hypotheses that are not particularly complicated, and that are not really "provable" anyway with the given data. We won't be able to understand conclusively how rearrangements happened, but showing Figure 4 makes things as clear as they can be. By contrast, Figure 3, which only shows STORR with no surrounding loci, is opaque. I think the manuscript could be considerably condensed with a lot less discussion about STORR --it was already previously studied quite extensively by Guo et al. 2018. Currently the entire section on "recruitment of new genes to BIA gene cluster locus" is only discussing STORR, which is really less interesting than the other stuff shown on Figure 4. With a good figure, this entire sentence could be condensed to a single sentence.
3. This paper uses a new and untested method to reconstruct the fissions and fusions in evolution. There is not a lot analysis presented about how the new "MO Solving" framework should work or testing of this on simulated/known data. It is good to see innovation within this paper, so I don't mean to knock this, but I would like to know whether the results are robust to possible errors there. What if you just counted breakpoints? This seems much simpler and less prone to inference errors --and it is known that such reconstruction is exceedingly error prone/lacks information for strong inference.
4. Line 189: While I agree that transposable elements may have facilitated rearrangements as those involved in STORR evolution, Figure S19 doesn't actually show anything interesting or provocative there so it should not be listed as evidence somehow in support of this idea. Unless it is showing some pattern that is not immediately obvious, in which case this should be explained more clearly. Figure 4: This is an amazing figure! The clear comparison among species is really the novelty that is introduced by this paper. It's very interesting to show that parts of this cluster exists in the other outgroups from somniferum. But, I'm concerned the authors are attempting to over-explain the results shown here and that this figure doesn't go far enough. How can you be so confident that there is only one tandem duplication involved here, and the rest are all dispersed duplications? This would seem to require some contrasting of the %identity between all copies present on the same chromosome vs. present on different chromsomes. This could be clarified by extending Figure 4 to show where the paralog copies are residing for the "dispersed duplicates". Given that some of these putative dispersed duplicates have gene names that are pretty close to each other, this needs to be clearly established. For example: -PSMT3 is "dispersed" but it's only 1 gene away from PSMT1. This could be tandem duplication of PSMT1 + CYP719A21 followed by deletion. -THS is classified as "dispersed" but it's syntenic with the paralogs of SALSYN and SALR on chr8 in somniferum, so where did the duplication happen?

5.
7. The discussion about cluster evolution needs to be improved. The study has not cited another paper published in Nature Communications on this subject that covers similar questions (evolution of clustering) and identified more extensive gene clustering, as well as patterns of differential expression (Li et al. 2020, Nat Comm; https://doi.org/10.1038/s41467-020-15040-2) and also focused on testing hypotheses about the evolution of these clusters. Some comparison of results is warranted, at the bare minimum. What did they conclude? How do their conclusions contrast with the conclusions here?
The current study has more extensive data with new genome assemblies, so it should be clarified whether the new results contradict or are concordant with previous results. Also, the paper should cite some of the other extensive work on cluster evolution in other species --for example by Osbourn, Nutzmann, and colleagues (many others as well). The paper should discuss the various hypotheses and evaluate more clearly what kind of data are needed to test them. If this is not done, then I don't see this paper as being suitable for publication --connections to existing literature are critical.
8. Figure 5A. The text on lines 221-228 seems important but is unclear --it is stated that the genes "arrived at the recipient locus perhaps pre-equipped with stem-specific promoter or regulatory elements" but this really doesn't make sense --if that were the case, wouldn't the donor loci also have the same pattern? How would they be pre-equipped if there wasn't a similar pattern in the donor locus? It seems like the "recipient prior" loci also have the same pattern in somniferum, so wouldn't it seem more parsimonious that the new loci took on the expression patterns of their local genomic environment? In any case, this is very murky. Li et al. (2020) found that some genes within a cluster had low patterns of co-expression. Do you find evidence that contradicts their claims? Why are the other downstream morphine production genes not clustered? 9. The Hi-C data does not seem conclusive. Hi-C is very noisy and trying to identify loops based on slight enrichment seems vague at best. I can see that it uses the tool HICCUPS but how confident can we be that this is correct? Is there any statement possible about the likelihood of fit to this model over others? Looking visually at the figures, I would have guessed other areas to be enriched and plotted other domains as being TADs/loops, etc. Unless there is some kind of rigouorus analysis of this (more replication, ideally), it should be acknowledged as highly speculative. I would like to see how chromatin loops have evolved in these different species but this would require more data. This is another area where I feel the authors are pushing too hard on establishing novelty and not enough on careful interpretation. 11. It would be useful to have a final summary somewhere of the genes involved in the morphinan/noscapine clusters and whether they are syntenic, tandem duplicates, or rearranged.
10. The manuscript needs extensive language editing throughout.
Minor comments: -Fig S14-S15: I don't see the value in synteny painting. It is somewhat interesting that it can be done, but basically all it shows is that there are many rearrangements. This does not lend itself to quantification or hypothesis testing.
-What tests are being used for the enrichment one line 130? Is this chisq test of all simultaneous? Does it correct for chromosome length? - Fig 1C: What do the colours on the dot plots mean? - Figure 4b: the two colours of blue for "new gene from dispersed duplciation" and "fusion translocation" are easy to confuse. Change the colour scheme.
- Figure 5A. It should be stated in the figure that this is for STORR, as it's only clear in the text.
Reviewer #2 (Remarks to the Author): In this study, Yang et al. generated chromosome-level genomes of three different poppy species that display extremes in noscapine and morphinan productions to study the evolution of these biosynthetic pathways. Cutting edge technologies were used to sequence and assemble high-quality sequences. These allow for the investigations of different whole-genome duplication events from around 7 million years ago. The authors have used different pieces of evidence to support their claims on the ancestral genomes and their evolution, rearrangements and recruitments of new genes to BIA gene cluster locus. Intriguingly, 15 genes involved in the two distinct pathways were assembled into such a compact gene cluster on one single chromosome 4 in opium poppy. Of particular interest is the investigation of the formation of the fusion STORR enzyme through fusion translocation events, possibly thanks to transposable elements. This article could be of great interest to audiences who are interested in natural product pathway evolution. Some minor comments: -This manuscript lacks an abstract and introduction.
-Line 54, page 4: what are the levels of noscapine in P. setigerum and morphine levels in P. rhoeas? Can the authors provide a figure or a table in the supplementary regarding the noscapine and morphine profiles of these three species? -Please define patchwork model more clearly.
-It is unclear whether the formation of the fusion STORR enzyme happens prior to or after the formation of the morphinan gene cluster. The author might want to add a statement regarding this event in the main text. - Fig. 4: the label of new genes from dispersed and new genes from "fusion, translocation" are quite similar. It is suggested that the authors choose different colours.
-It is interesting to see that new genes that arrived at the recipient locus might be pre-equipped with the stem-specific promoter or regulatory elements. With the high sequence quality that the authors have, it would be interesting to know what these promoters and regulatory elements would be. Can the authors provide more information regarding this? -Can the authors elaborate more on the potential epigenetic regulatory mechanism for simultaneously regulating genes from the two branches? It would significantly increase the novelty of the work and would be of great interest for the audience who read this article to have deeper information regarding these mechanisms.
-What are the levels of identity of the genes from these different species? -Is this patchwork model common in other biosynthetic pathways? Can the authors comment on this fact on the reported clusters so far?
Reviewer #3 (Remarks to the Author): Reviewers comment for Nature communication article-In the present study under review, titled, "Three chromosome-scale Papaver genomes reveal punctuated patchwork pathway evolution leading to morphinan and noscapine biosynthesis", authors reported whole-genome assembly for two Papaver species, namely P. setigerum and P. rhoeas, and worked to improve previously published genome for P. somniferum. These three species represent extreme in the sense of morphine and noscapine synthesis. The authors then used these three genome assemblies to reconstruct the ancestral genome for Papaver species and derived a hypothesis for the evolution of BIAs in plants. The authors proposed the role of whole-genome duplication, followed by chromosomal rearrangement, fusion, and translocation as key events towards the innovation of secondary metabolites. The study is well written, and all figures are of good quality. I also need to mention here that overall, the manuscript is super-condensed and sometimes hard to follow, especially with the transition from ancestral genome construction towards the proposed hypothesis of the role of fusion and translocation towards the evolution of specialized metabolites. I highly recommend to consider expanding the explanation, simplification of text, and improving the flow if the argument with more detail in order to attract general and specialized readers towards this study. One of the weaknesses of this study is the generality, which is entirely missing. The offered method for ancestral genome construction, the whole hypothesis of evolution of specialized metabolites, and all discussions are Papaver specific. I feel that the authors have not attempted to draw comparisons with other studies or to draw a general point-of-view based on their results and other high-quality genomes that have also offered different approaches towards the evolution of specialized metabolites. This is one of the significant weaknesses of this study that I feel must be taken care of. While going through this study, I find several questions unanswered. I am listing these in the same flow as these are described in the manuscript. Major comments-1. I am wondering as why authors choose to use P. somniferum genome that they improved in this study while a recent publication has already shown a much-improved genome assembly of P. somniferum (https://www.nature.com/articles/s41467-020-15040-2#Sec25). In this study, authors showed significant improvement in overall genome assembly compared to what were published in the Science article (Guo et al., 2018). Furthermore, the contig N50, which I personally believe is one of the best criteria to judge a complete and well-ordered genome, is significantly better in Li. et al, 2020 (7.6 Mb compared to 1.74 Mb of this P. somniferum genome reported in this study). To me, not including that genome for comparative analysis is not a good idea, and I wonder what authors have the reason to not include it in this study. I am not sure if using this genome changes any of the hypothesis that authors have proposed in this study. 2. Authors described different levels of morphine and noscapine across the three species (Supplementary Figure 2). As the entire manuscript premise is about evolution of the specialized metabolites biosynthesis using these three plants, I feel that looking at overall metabolites and metabolite intermediates to morphine, noscapine across these species will be important. I am particularly curious as how about the intermediates from the earlier steps of biosynthesis pathways of key metabolites and their levels across these three species. Probably such quantification would allow to associate evolution of genes with presence/absence of the metabolites. If these metabolites/intermediates have already been reported and described for these plant species in the past, would be nicer to include as part of discussion and authors views on association of presence/absence of genes with the metabolites. If metabolites are not identified yet, authors should consider identifying and quantifying these intermediates. 3. I am curious about the difference between the different level of contiguity among the genome assemblies reported in this study. Based on my perception, I expected P. sentigerum genome to have more difficulties in terms of getting assigned to the scaffolds while P. rhoeas to be the easiest based on the genome size, number of chromosomes and so on. But this is just opposite, and 97.6%, 92.4% and 87.9% of contigs for P. setigerum, P. somniferum and P. rhoeas, respectively. I wanted to ask authors their opinion on this. I am assuming that this has to do with raw reads quality (ONT raw read N50), and some other factors. If this could be discussed in the supplementary method section, that would be helpful for people who may want to read this study to get inspiration for a genome assembly project. 4. This is a minor point, but in Table 1, authors reported number of protein coding genes across these three species. The numbers are quite comparable between P. somniferum and P. rhoeas, while the former have undergone whole genome duplication while the latter has not. While we know that whole genome duplication is further followed by massive genome readjustment, including deletion and rejection of gene families and so on, I am curious to know as what genes specifically were gained and lost across these species. Are they somehow being favored for BIA biosynthesis? What kind of processes were lost in P. somniferum? Authors have briefly described it in section 10 of supplementary information, but I feel that a direct comparison between these three species in terms of understanding gene gain and loss would be interesting as they do have a contrasting levels of BIAs despite being closely related. I also feel that including or expanding this aspect in the discussion will make the discussion section better to understand evolution of BIAs as emerging through this study. 5. Authors mentioned in the method and supplementary method section that they used purge_dup to exclude any redundancies from the genome assembly. I think that its important to provide parameters used. As purge_dup may also exclude some of the genes that may be real and not artifacts, the used cutoff threshold needs to be mentioned. It may be a good idea if authors could provide all scripts and parameters for different stages of analysis/assembly/ comparative genomics as a github repository as they have provided for reconstruction of ancestral genome. 6. This is just a minor comment. For phylogenetic tree, authors have used eight angiosperms. I feel that including more plant species for the phylogenetic tree by including plants from different lineages and from different time of evolution would be nicer to provide a relative view of the evolution of Papaver species with others. This is simply my preference, and I leave this up to the authors to agree or disagree. 7. I am not sure if I have understood Fig1d. I know that none of the Papaver species have undergone whole genome triplication, hence no peak around 1.5~2 Ks. Normally, we expect a peak that corresponds to a whole genome duplication as probability of sharp increase in Ks value increased with it, hence, I was expecting two peaks for P. setigerum, while single peak for P. somniferum (this is inline with what we have observed in case of Arabidopsis, which shows three peaks, one that represents whole genome triplication while rest of the two peaks corresponds to the two whole genome duplication events). If authors do not split P. setigerum genome as WGD1 and WGD2 set, how does the peak looks like? 8. Much of this article, including hypothesis and interpretations for comparative genome analysis replies on the constructed ancestral genome of Papaver species. At the beginning of section, "Ancestral genomes…", author says, "We developed a novel bottom-up framework to reconstruct preand post-WGD-1 ancestors of Papaver genomes". I have followed work mainly from Jerome Salse, and I wonder how this method is different from what his group have been using and have proposed before? A comparison in terms of what difference between this approach and what previously have been used and providing further discussion on advantages and novelty would be very important. 9. Authors provided github repository for the codes used in this study, and they say in the text as well as in the description of github repository that this code is suitable for reconstruction of ancestral genome of Papaver species. I wonder if the codes could be used for other plant species. Does this code be applicable for including and comparing different plant species and to derive ancestral genomes for a certain lineage? Authors need to clarify this in the text, and if this is specific to Papaver species, why? This is a topic that would be of wider audience interest, and if addressed would be able to make this paper even stronger and would add general appeal. 10. Authors reported six and eleven protochromosomes for the pre-and post-WGD-1 event. A widely believed concept using reconstructed ancestral genomes using species known to have undergone minimal rearrangements proposed seven protochromosomes and described as ancestral eudicot genome having seven chromosomes while ancestral monocot genomes having five chromosomes prewhole genome duplication (Murat et al., 2017). I am curious as how these reconstructed genomes are different from what authors have reported here. Authors have compared AEK with these plants, but how about the reconstructed genome. A comparison is essential as a lot of hypothesis and interpretations depends on this aspect as well. Also, when authors would predict the origin of ancestral genomes of Papaver species (I mean around what Millions of years ago on the time scale). 11. Authors have described chromosomal fissions and chromosomal fusions associated with the transformation of pre-WGD-1 protochromosomes to modern chromosomes. Given that the genome continuity is vastly different between these species, on what basis we can assign a confidence score for the detection of fusion event sites, which could very well be due to misassembly. Another point, in a recent study on Chromosome assembly of Ophiorrhiza pumila, authors experimentally validated genome (assembly gap as 21). In this study, authors identified orientational error across the assembly gap, which were not identified using Bionano and Hi-C datasets. My question is this: How authors could order or talk about the fusion event sites and the shuffling breakpoints when one may very well question the correctness of contigs orientation within assigned scaffolds, and when the percentage of contigs assigned to chromosomes ranges between 87.9% to 97.6%? 12. Authors have very elegantly described the origin of STORR in the section "Recruitment of new genes to BIA gene cluster locus". However, as a reader, I would be interested if similar events and similar processes have ever been reported in other plant species. In terms of situation in other plants, is there any report which have observed such scenarios, or this is the first result that have observed means of evolution? 13. Role of transposons in the evolution of BIAs in these plant species are not explored. This is another week point that I wish authors should work. I agree that the synteny analysis and translocation of fusion protein does make sense, but for me, its hard to believe that the transposons have not played any major role when almost 3/4th of the entire genome of all these species are constituted of transposons. A old study published in PNAS on wild tobacco and comparative genome analysis showed role of transposons towards evolution of nicotine, and that study showed that incorporation of specific transposons to the promoter regions promoted or disrupted expression of key genes. I believe that an entire new section on how transposons distribution across promoter and genes of these species, and its interpretation needs to be explored. I disagree that WGD alone could derive such a fascinating metabolic pathway, and including other players are important to make this study complete. Authors have one sentence in page 15, line number 188, and a supplementary Fig.  19, but I find it not described enough and missing interpretations. A more detailed discussion about authors interpretation on role of TEs in these species would be very helpful. 14. I wanted to ask about the Hi-C experiment that authors used to describe epigenetic factors in the last section. I did my best to find out number of replicates that were used for this interpretation but could not. I hope that authors used replicates in the Hi-C experiment, and only then, the interpretations were derived. I will not like to believe this data if replicates are not included for the interpretation. 15. Since expression of genes involved in the BIAs biosynthesis are known to be tissue species, and this is what authors also reported by saying that the epigenetic factor played a role in getting these genes highly expressed in a tissue specific manner. Do authors have any reference to support their point other than the comparative HiC data for these three species? What I mean is that it would be clearer if authors could show Hi-C data for tissues know to have no expression of BIA biosynthesis associated genes and metabolites, and tissues that show highest accumulation of BIAs and expression of genes. Those comparisons will be better to identify interactions at chromatin level (even if this is done for P. somniferum), which then could be explored across these three species to derive the conclusion that authors have made here. 16. I feel that discussion need a comprehensive overview of current study standing with what has been done previously. For example, previously published study on P. somniferum (https://www.nature.com/articles/s41467-020-15040-2#Sec25) reported copy-number variations for key enzymes towards biosynthesis of BIAs across producing species. Authors should discuss as what they think about those key enzymes described in that study, and if possible, evolutionary scenarios for them. Further, standing of this study with other specialized metabolic pathway is needed to provide a more general point of view on the evolution of secondary metabolism. I could find several interesting aspects being presented in the supplementary information, which should be included in the discussion section.
These are few of the minor comments and are not necessary in the order they appear. Minor comments-1. Scale bar is missing in the Supplementary Fig 1. 2. All Hi-C figure panels that authors have shown across all supplementary figures are not clear and difficult to get a sense of quality of the genome. Supplementary Fig. 4b, 6b, and 7b are not clear at all, and ideally should have, axis labels. The authors of the article "Three chromosome-scale Papaver genomes reveal punctuated patchwork pathway evolution leading to morphinan and noscapine biosynthesis" provide the genomes of Papaver rhoeas and Papaver setigerum and by focusing on the evolutionary fate of genes involved in benzylisoquinoline alkaloid (BIA) biosynthesis they also provide very interesting insights into how a cascaded pathway evolves. They conclude, that BIA biosynthesis follows the patchwork model: metabolic pathways assemble via the recruitment of primitive, promiscuous enzymes that can react with a wide range of chemically related substrates. To completely follow this conclusion, it is necessary to elaborate in more detail on the phytochemical space of P. rhoeas, P. setigerum and P. somniferum. P. rhoeas resembles the ancestral state of the three species most, as it did not experience a whole genome duplication events, which often accelerates evolution. Which type of BIAs are produced in P. rhoeas? (see specific comments). Does P. rhoeas possess a primitive metabolic pathway that enables the production of noscapine and/or morphinan derivatives? I also suggest to the authors to elaborate more on the biosynthetic steps. Possibly provide a schematic representation of the noscapine and morphinan pathway and as a phytochemist, I miss the chemical structure of at least the major BIAs, e.g. noscapine and morphine. These structures illustrate best the differences/similarities of the major BIAs, and the metabolic steps that "separate" them. The genome sequencing and data analysis is profound. I only can't completely follow the reconstruction of the ancestral genomes, which might be due to my limited expertise in this field. Still, I would recommend to revise this part, as the authors want to address a broad audience. Also -I am not a native speaker -but I would recommend language/linguistic editing of the manuscript. I highlighted some of my language/grammar concerns in the specific comments. Summarizing, I believe the data generated in this work is very valuable and interesting to the scientific community, but some points should be revised -see my specific comments.
Specific comments: Page 3: Line 22: "While single variant could create a new gene" -this is very unvague, what is meant exactly? Line 25: "various numbers of whole genome duplications" -this phrase is a bit exaggerated, please be more precise -one WGD shared by P. somniferum and P. setigerum, and a second only in P. setigerum. Line 27: "nonrandom distribution towards innovation of secondary metabolism" -in P. somniferum! See results part, page 10, line 119. In P. setigerum -genes in breakpoint vicinity were enriched in plant-pathogen interactions. Page 4: Line 40: please elaborate a bit more on the morphinan and noscapine alkaloids -as mentioned above. Biosynthesis, chemical structure, occurrence in the three Papaver species. Line 47: "and compared them with P. somniferum genome". Line 52: Chapter "Genome assembly and annotation" -describing the pattern of morphine and noscapine accumulation in the model species P. somniferum, P. rhoeas and P. setigerum better fits to the introduction (see above). Also, check the sentence and cite references! And the cross reference to Fig. 1a is a bit misleading at this point (line 55), as the Fig. 1. does not illustrate alkaloid accumulation. Which, by the way, could be helpful. For example, the "alkaloid-type" could be included in Fig. 1e Fig.1c: Please include color code/legend that was used to color the syntenic blocks. I assume its Ks? Fig.1d: What should the close up illustrate? And its very hard to differentiate between the colored lines representing syntenic orthologs of P.som. -P. set. and P. set. -P.rho. in the figure. Please adjust the colors.   Fig. 3b, the grey boxes in the recipient loci -are they necessary? Graphic legend Fig. 3a: the blue color of STORR is hardly to distinguish from the oxidoreductase pre-fusion module, anyway -would it not be more correct to give the same symbol for STORR in the graphic legend as in the figure itself (lilac rectangle together with blue arrow)? Fig. 4: The color code is confusing in in panel b -new gene from "fusion, translocation" same color as new gene from dispersed duplication. "dispersed duplication" -are paralogs present somewhere else in the genome?? Please explain what you mean and please elaborate the origin of this genes in more detail in the text! Extended Data Fig. 3: Please explain the abbreviation GMP, and pPG. In general, the figure is not very intuitive. The figure should be able to stand alone. I would advice to revise it. Extended Data Fig. 4: Please include color code/legend that was used to color the syntenic block. And could you elaborate it bit on the findings, that we can get from the dotplots in the figure legend? E.g. in 4b -P. somniferum plotted against the post WGD-1 ancestor -chromosome 1 shows a syntenic block ration of 1:2, or in other words: a syntenic blocks in chromosome 1 P. somniferum "matches" chromosome 1 and chromosome 6 in the post-WGD-1 ancestor. What is the explanation? Small-scale duplication?

With kind regards! Elisabeth Kaltenegger
Thank you for all your valuable comments. We have provided a response letter addressing all the issues raised by the reviewers. For clarity, all reviewer comments or quoted contents are in italicized fonts. A point-to-point response to each comment is provided in normal fonts. References to revised manuscript contents are also provided where needed. Please notice that the

Reviewer #1 (Remarks to the Author):
The authors have presented a comparative genomic study of three species of poppy using new genome assemblies to study how to genes involved in alkaloid metabolism have evolved over time. This is quite an interesting area of research and it's good to see the new genome assemblies. Figure 4 is the most important contribution of this paper and on this piece of evidence alone, this paper merits eventual publication --but this figure has some problems and could be substantially clarified (see comment #5). I think greater care should be taken to not overstate the significance of any results and to take a balanced view of what they mean. It is more interesting to just show us the data and allow the data to speak for themselves, rather than trying to build elaborate stories about how something might have happened when the evidence isn't really conclusive supporting a given sequence of rearrangement/duplication/fusion events (e.g. comment #2). I find the authors are often seeming to push a little too hard on the importance or strength of their conclusions, perhaps as they are worried about being perceived as sufficiently novel/impactful for Nature Communications? This is detrimental to the paper. Also, the question of cluster evolution has been considered extensively by other papers and if this is the aim of this paper, more connection should be made to this literature (see comments 7 & 8). A major question in the evolution of gene clusters is how commonly they have evolved from tandem duplication vs. rearrangement. It is critical that this aspect be very rigourously assessed. I think with some substantial re-working this manuscript could be suitable for publication, but it would really need to push less on promoting one particular story and take a much more balanced view of the data. I would appreciate if the authors approached their data with a bit more hesitancy and acknowledgement of the areas of uncertainty --this is fine! Having two new genomes provides lots of scope for assessing how rearrangements happened and is sufficiently novelty to merit publication in Nature Communications, but only if this is discussed in a more balanced and careful way. Figure 4 shows the BIA gene cluster evolutionary history and show the detail evidences of each gene in Supplementary Figures  22-30 in the revision manuscript. In addition, we improved the previous draft genome of P. somniferum and de novo assembled two additional Papaver genomes, P. setigerum and P. rhoeas. We evaluated the genomes and found they were high quality and in chromosomal-level. We did extensive synteny analysis on these three genomes, and decoded the evolutionary history of these three Papaver species. We found two rounds of whole genome duplication (WGD) events with one being shared by P. somniferum and P. setigerum at around 7.2 Ma and another being P. setigerum specific at around 4.0 Ma. Moreover, thanks for the three chromosomal-level genomes, we decoded the evolutionary history of BIA gene cluster and found it can be explained by the patchwork model.

RESPONSE: Thanks for your interests on our work and indeed
We summarized the BIA gene cluster evolutionary history in Figure 4 and show the detail evidences of each genes in Supplementary Figures 22-30. Thanks for your comments, and we have revised the Figure as your suggestion. In addition, we constructed the identity heatmap of all genes in noscapine and morphinan clusters, as well as their close paralogs ( Figure R4) and added it as a new Supplementary Figure (Supplementary Figure 31). The identity patterns father supports our conclusion of the evolution of genes in BIA genes cluster.
We have revised our manuscript to tone down our statement. Based on our data, we did extensive analysis to reach the conclusions, e.g. the evolutionary history of STORR and the BIA gene cluster. In our opinion, these results are the reasonable inference and stories based on the data.
Compare to previous paper, e.g. Guo, 2018, Science; Li et al. 2020, Nat Comm, we rescaffold P. somniferum draft assembly using Hi-C data from the same cultivar (HN1), and de novo assembled two additional Papaver species. Based on the three high-quality genomes, we investigated the evolutionary history of BIA gene cluster and reveals novel insights into how BIA gene cluster was formed and evolved among different species, which is unexplored by previous paper. We have revised the discussion to link the available researches on gene cluster evolution.
We investigated the origin of genes in BIA gene cluster by integrating evidence of multiple sources including synteny, phylogeny and WGD. We have provided Supplementary Figures 22-30  and Supplementary Tables 13-15 in the revision manuscript to show the original results. We summarized the results in Figure 4. From our results, we found 6 genes in the BIA gene cluster were already presented in the MRCA of three Papaver species. During the evolution, nine genes were assembled to this gene cluster. Except STORR was resulted from "fusion, translocation" events, we have inferred that PSCXE1, PSAT1, PSMT3, SALAT, and THS were resulted from "dispersed duplication", and CYP82X2 was resulted from "tandem duplication". Moreover, two genes, PSMT2 and CYP82Y1, was not obtained a clear origin. We defined the "tandem duplication" as the one genes' origin is adjacency with it, and others defined as the "dispersed duplication". We also tried our best to find the TEs associated with BIA related genes, and marked them on Supplementary  We have revised our manuscript as suggested, and provided a point-by-point response to your questions.

Major comments:
1. I found it very difficult to understand the results of the paper without going to the supplementary materials. I think this paper is much more suitable for a longer and more in-depth format in a specialty journal. Most of the analyses are descriptive with little in the way of hypothesis testing --the reconstructions of orthology relationships don't conclusively support one model of gene evolution over another, as this would require sequencing of more species.
RESPONSE: Thanks for your comments. As you suggested, we have revised our manuscript by adding more details. We apologize for any unclarity in our manuscript, which we intended to write in a way that keeps the main manuscript concise and provides abundant supporting evidence in the supplementary materials. In this way, the readers may grasp the main conclusions of the main text, while those who are interested in more technical details are directed to read the supplementary materials.
We agree with you that sequencing more species could bring more intermediate states into the analysis, probably allowing us to capture a finer picture of evolutionary events for each gene in the cluster. However, based on the cutting-edge sequencing data and our extensive analysis of three high-quality Papaver genomes, we believe the comparison among them already yields sufficient evidence that the so called "patchwork model" is a decent working model to explain the evolutionary history of benzylisoquinoline alkaloid (BIA) gene cluster. That being said, we expect that the evolutionary history of BIA gene cluster would be updated after more Papaver species being sequenced, which will be our future research direction. We have revised the manuscript in a more open-minded manner, and tuned down our conclusion about the BIA gene cluster evolutionary model.

The previous paper by Guo et al (Science 2018) focused extensively on evolution of STORR. It is unclear what is contributed by this further analysis --can the authors more clearly describe what this shows that was not known previously?
RESPONSE: In our previous Science paper (Guo et al. 2018), we, for the first time, reported a chromosome-level genome assembly of the P. somniferum, from which a large gene cluster harboring 15 genes involved in the BIA biosynthetic pathway was discovered. In addition, we also showed opium poppy underwent a whole genome duplication at around 7.8 million years ago that could contribute to the clustering of BIA genes. However, due to a lack of genome sequences for additional Papaver species, in the 2018 paper we were not able to get a clear picture of the evolutionary steps leading to the assembly of the entire gene cluster formed.
In this current manuscript (Yang et al.), we significantly improved this draft genome assembly by incorporating Hi-C data and increased the chromosome anchor rate from 81.6% to 92.4%. This has effectively allowed us to correct assembly errors in the original assembly and enabled the investigation of the complex events involved in evolution of BIA gene cluster. For example, regarding the evolution of STORR, we identified the pre-fusion module of STORR in 2018 Science paper ( Figure 2D in Guo et al. 2018). Using P. somniferum genome only, we inferred the "STORR and its closest paralogs show amino acid sequence identity of 75 and 82% for the P450 and oxidoreductase modules, respectively, which suggests that the duplication leading to the STORR gene fusion occurred earlier than the WGD event" (Guo et al. 2018). Given only one genome at that time, we could only rely on the sequence identity information to draw the conclusion, which was the best explanation at that time. As you correctly pointed out in Comment #1 that sequencing more genomes will improve the evolutionary model, we assembled another two Papaver genomes, and did an extensive syntenic analysis among these three species. We found that the "duplication" in the Science 2018 STORR evolutionary model "duplication leading to the STORR gene fusion occurred earlier than the WGD event" is actually the WGD shared by P. setigerum and P. somniferum (Figure 3a, Supplementary Figure 20 in the revision manuscript). We come up with an updated hypothesis of STORR evolution, e.g., a "fusion, translocation" event after WGD leading to STORR formation at BIA gene cluster ( I found the paragraphs and figures on STORR particularly confusing and unclear. The model of "Fusion-Translocation" seems like one of many possible sequences of events and this seems more like evolutionary story-telling than hypothesis testing. In fact, given the figures as they are presented, it seems less likely than an alternative hypothesis. As I understand it, here is the presented FT-hypothesis: -two parent genes leading to STORR were present on same chromosome prior to WGD -WGD occurred in ancestor of setigerum and somniferum (not rhoeas) -some "Fusion, translocation (FT)" event happened (which is not described but says it may involve some translocation event) The authors claim the evidence in Figure 3 and Figure S17 supports this claim. But I fail to see how this is more parsimonious than simply having some deletion between the intergenic regions of the "pre-fusion" modules? This would be SO much clearer if the paralogs to the pre-fusion modules within somniferum were shown on Figure 4. Where are they? I can only assume they must be on some other chromosome, which would actually be much clearer support for the "translocation" part of the story than the convoluted Figure 3 and accompanying explanation. Alternatively, expanding Figure 3 to show the chromosomal context more clearly might be useful. It's clear from Figure 4 that there are translocations involved becuase the noscapine cluster got fused with the morphinan cluster. I generally feel like this manuscript spends a lot of words trying to prove hypotheses that are not particularly complicated, and that are not really "provable" anyway with the given data. We won't be able to understand conclusively how rearrangements happened, but showing Figure 4 makes things as clear as they can be. By contrast, Figure 3, which only shows STORR with no surrounding loci, is opaque. I think the manuscript could be considerably condensed with a lot less discussion about STORR --it was already previously studied quite extensively by Guo et al. 2018. Currently the entire section on "recruitment of new genes to BIA gene cluster locus" is only discussing STORR, which is really less interesting than the other stuff shown on Figure 4. With a good figure, this entire sentence could be condensed to a single sentence.

RESPONSE:
We appreciate your comments and suggestions. The evolutionary model of BIA genes we proposed in this manuscript is derived from a thoughtful explanation of our observation from the multi-genome comparison. Based on thorough analysis of three high-quality genomes, we rejected the evolution model of STORR presented in Guo et al, Science, 2018 and come up a "fusion, translocation" model. Here, we will present a brief explanation of why our proposed model is the more likely hypothesis and why the hypothesis of "simply having some deletions" conflicted with what the data shows.
First, if STORR formation is caused by a "simply some deletions between the intergenic regions", the STORR flanking genes must have synteny relationships with the flanking genes of pre-fusion modules ( Figure R1). However, we did not detect any such syntenic relations between the pre-fusion model loci and STORR loci (the donor loci and recipient loci in our manuscript) (  Second, since the "deletion" hypothesis does not stand, we proposed the "fusion, translocation" model based on the three chromosomal-level genome comparison, the supporting evidences as following:  Lack of synteny between the STORR loci and pre-fusion module loci in P. somniferum and P.   As for the section "Recruitment of new genes to BIA gene cluster locus ", we believe "FT" event leading to STORR at BIA gene cluster is one of the most significant findings in this work. To our knowledge, this phenomenon is the rare in gene cluster evolution. Moreover, the translocation of STORR from loci of largely low and non-tissue-specific expression to the BIA gene cluster displaying a high and coordinated expression in stem is striking. Thus, we mainly focus on STORR in this section.

This paper uses a new and untested method to reconstruct the fissions and fusions in evolution.
There We simulated the evolutionary scenario from top to bottom with two WGDs consistent with three Papaver species ( Figure R2a) under infinite sites assumption (IS), which means a mutation does not occur at the same locus more than once during evolution and is commonly used in evolutionary studies (Aganezov S. et al, Genome Res. 2020). We simulated block sequences with some random block adjacencies change between each species. Here, we required that the endpoints involved in changes do not overlap to satisfy IS assumption. And then, we applied our model to reconstruct each middle species, e.g. Species 2, Species 3, and Species 5 in Figure R2a. We repeated 200 times and found that Species 2 (simulated pre-WGD-1 ancestor) and Species 3 (simulated post-WGD-1 ancestors) can be reconstructed with 100% block adjacency accuracy, and Species 5 (simulated pre-WGD-2 ancestor) can be correctly reconstructed with average 99.68% block adjacency accuracy ( Figure R2b). This result indicated the accuracy and robustness of our framework under IS assumption.
Next, we evaluated the block adjacency reliability for three ancestors in real Papaver evolutionary scenarios. We found all block endpoints in the reconstructed pre-WGD-2 and post-WGD-1 ancestors satisfied IS assumption. We inferred that the block adjacency reliability of both pre-WGD-2 and post-WGD-1 ancestors were 99.68% (pre-WGD-2 ancestor is 99.68% and post-WGD-1 ancestor is 99.68%×100%) based on the simulated results under IS assumption. We adjusted the block adjacency reliability by accumulated multiplication bottom-to-up. However, the pre-WDG-1 ancestor has 11.67% non-IS block endpoint. So, we simulated the pre-WGD-1 ancestor reconstruction under non-IS assumption 1000 times with non-IS block ratio from 0 to 100% ( Figure R2c). We used quadratic polynomial to fit the correlation between non-IS endpoint rate and endpoint adjacency inconsistence rate, and obtain the fitting curve with R 2 of about 0.99. Finally, we estimated the reconstructed endpoint adjacency inconsistence rate of pre-WGD-1 ancestor being 5.70%. Therefore, the adjacency reliability for this ancestor is 94.0% (99.68%×100%×(1-5.7%)). So, pre-WGD-1 ancestor may have two endpoint adjacencies inconsistence ((1-94%)*36=2.16).
We have revised Supplementary Note 9 in Supplementary Materials. We added Figure

What if you just counted breakpoints? This seems much simpler and less prone to inference errors --and it is known that such reconstruction is exceedingly error prone/lacks information for strong inference.
In some simple cases, counting breakpoint seems to be an easy way to reconstruct ancestor. However, WGDs in Papaver species created multi-copies of syntenic blocks. The reconstructed ancestor must include all block endpoints from three Papaver species ( Figure R3a). In addition, the ancestral genome reconstruction is a global optimization process rather than local optimization, while simple counting breakpoints is a local optimization strategy. Here, we proposed two examples to show the difference between breakpoints counting and our method ( Figure R3). The first one shows that if we only count and select the breakpoints with the most occurrences by a greedy strategy, some endpoints will be isolated leading to errors in reconstruction ( Figure R3a). Therefore, the reconstruction process should have constraints to ensure every endpoint be considered. Secondly, if we just count breakpoints, the final reconstruction may include errors. For example, there are four blue-red breakpoints and three blue-green breakpoints ( Figure R3b). If we just count breakpoints, the final ancestral connection will be blue-red due to the higher count ( Figure R3b). However, in our bottom-to-up process, we are able to get the accurate ancestral connection with blue-green and the intermediate stats (e.g. Species 5 and Species 3). If the final state is blue-red, there are at least three shuffling events to reach the final states. However, if the final state is blue-green, only two shuffling events are required to reach the final states. Therefore, we think the ancestral genome as blue-green is the correct one and the breakpoint counting seems not the best method.
In summary, we added the endpoint constrains in integer programming model to avoid isolated endpoints and make all endpoints in final reconstructions. And our ancestor reconstruction followed a bottom-to-up process based on MO framework which can avoid fall into local optimum. Although we are confident about our ancestor reconstruction pipeline and its associated results, if reviewer #1 and the editor insist our ancestral genome reconstruction is not robust, we are willing to remove the content related to it.

Line 189: While I agree that transposable elements may have facilitated rearrangements as
those involved in STORR evolution, Figure S19 doesn't actually show anything interesting or provocative there so it should not be listed as evidence somehow in support of this idea. Unless it is showing some pattern that is not immediately obvious, in which case this should be explained more clearly.

RESPONSE:
Thanks for your comments. We detected the DNA transposable elements around STORR and pre-fusion module loci in three Papaver species ( Figure S19, which is updated as Supplementary Figure 21 in the revision manuscript). Papaver species are rich in transposable elements, and we detected 17 types of transposable elements distributed around the donor and recipient loci. We speculated that these TEs may have been involved in the genome rearrangement around these loci. However, we do not have solid evidence about which TE contributed to the STORR formation. Therefore, to be cautious, we did not make overstatement about the role of TEs in the STORR formation. For CYP82X2, we found the reciprocal best BLASTp hit is CYP82X1 with sequence identity of 58% (Supplementary Figure 29, Supplementary Table 14 in the revision manuscript). Since the two genes are adjacent to each other in the genome, we reason that CYP82X1 or CYP82X2 could arise from a tandem duplication event.

Figure 4: This is an amazing figure! The clear comparison among species is really the novelty
For PSMT3, we found the best BLASTp hit is Pso05G43960.0 with sequence identity of 80%, located on chr5 (Supplementary Figure 28, Supplementary Table 14 in the revision manuscript). The sequence identity between PSMT3 and PSMT1 is as low as 30.32%. In addition, we detected the synteny relations of Pso05G43960.0, and found one syntenic pair in P. somniferum, four syntenic pairs in P. setigerum, and one syntenic pair in P. rhoeas, indicating Pso05G43960.0 was present in the MRCA of three Papaver species. Taken together, it indicates that PSMT3 was duplicated from Pso05G43960.0 by a dispersed duplication event.
THS's non-syntenic best BLASTp hit is Pso04G09740.0 located 34Mb downstream of BIA gene cluster with percentage identity of 67% (Supplementary Figure 23, Supplementary Table 14 in the revision manuscript). In addition, we detected six synteny pairs of Pso04G09740.0 in three species (one in P. somniferum, four in P. setigerum, and one in P. rhoeas), indicating Pso04G09740.0 is present in MRCA of three species. These results suggest THS was likely duplicated from Pso04G09740.0 by a dispersed duplication event. And for SALSYN, and SALAT, we have shown that they were already present in MRCA of three species. In summary, six genes in BIA gene cluster including PSSDR1, CYP82X1, CYP719A21, PSMT1, SALAYN, and SALR already existed in the MRCA of three Papaver species. Five genes, including PSCXE1, PSAT1, PSMT3, SALAT, ant THS, were assembled into BIA gene cluster by "dispersed duplication". CYP82X2 was tandem duplicated from CYP82X1, and STORR was created by a "fusion, translocation" event.

RESPONSE:
Thanks for your comments. We totally agree that a comparison between our work and the work by Li et al. 2020 Nature communications you mentioned should be conducted. In fact, the analysis conducted by Li et al. only focused on P. somniferum cultivars, while our work comprehensively compared genomes of three Papaver species. Therefore, many of our findings in genome evolution of Papaver genus in this manuscript (Yang et al.) are simply out of reach for Li et al. due to the differences of plant materials and computational analysis involved. That said, we would like to point out several common and different findings in the two papers as you rightly suggested.
Here are the main findings or conclusions from Li et al paper: 1. A rescaffolded genome assembly of P. somniferum based on our draft genome (Guo, et al. Science, 2018) was achieved using Hi-C data, improving the proportion of contigs anchored to chromosomes; 2. Gene clustering for BIA biosynthesis genes and co-expressed in stems and root tissues; 3. Co-expression of BIA genes is correlated with alkaloid contents in P. somniferum tissues; 4. Copy number variation is found in BIA biosynthesis genes of P. somniferum cultivars, well correlated with BIA accumulation in these cultivars. Now we would like to make a pairwise comparison of these conclusions between Li et al paper and our current manuscript.
Regarding the first, second and third conclusions, we have published these conclusions in our Science 2018 paper, and Li et al. improved our draft genome assembly using Hi-C (generated from a completely different cultivar PS7) and echoed our finding of the BIA gene cluster, but how BIA gene cluster was formed and evolved in Papaver genus is unexplored, simply due to a lack of genome sequences for additional Papaver species. In this manuscript, our Hi-C data were generated from the same cultivar that was used for the initial genome assembly published on Science 2018, therefore, we chose not to use the Li et al. assembly, and trust our own improved genome assembly of P. somniferum, because we believe that it is best to rescaffold the genome using Hi-C data of the same cultivar.
Regarding the fourth conclusions, what Li et al did nicely is to include resequencing data for nine P. somniferum cultivars and detected copy number variants. They showed that copy numbers of morphinan biosynthetic genes are positively correlated with morphinan accumulation in these cultivars. In our analysis of three Papaver genomes, we noticed there are different copy numbers of morphine biosynthesis genes in P. rhoeas, P. somniferum and P. setigerum ( Figure R5). The least copy numbers (zero or one) found in P. rhoeas account for the low accumulation of morphinans, whereas P. somniferum that has more copy numbers than P. rhoeas does accumulates morphine abundantly ( Figure R5, Supplementary Figure 2). We also found that P. setigerum has more copy numbers than P. somniferum for many morphine biosynthesis genes ( Figure R5). However, the production of morphinans is lower than P. somniferum ( Figure R5, Supplementary  Figure 2), given that the accumulated expression levels of the P. setigerum copies are lower than the corresponding genes in P. somniferum ( Figure R5). Therefore, our results show that the accumulated expression levels, in additional to the copy numbers, are well correlated with morphinan production in different Papaver species, updated the conclusions in various P. somniferum cultivars by Li et al. paper. We added the Figure R5 as Figure 6 in our revised manuscript and added a section "Accumulated gene expression contributes to morphinan production" in the revised manuscript.
In summary, given the genome data from two extra Papaver species, our study reveals novel insights into how BIA gene cluster was formed and evolved among different species, which could not be accessed by using a single species analysis in Li et al. paper.
We have made revision to our discussion in the manuscript on Page 22-23. Also, the paper should cite some of the other extensive work on cluster evolution in other species --for example by Osbourn, Nutzmann, and colleagues (many others as well). The paper should discuss the various hypotheses and evaluate more clearly what kind of data are needed to test them. If this is not done, then I don't see this paper as being suitable for publication --connections to existing literature are critical.

RESPONSE:
As you suggested, we have revised our manuscript, and added discussion with the previous work, including Li et al 2020, and other cluster evolution related work. The discussion has been revised as the following on Page 22-23 line 355-364. Figure 5A. The text on lines 221-228 seems important but is unclear --it is stated that the genes "arrived at the recipient locus perhaps pre-equipped with stem-specific promoter or regulatory elements" but this really doesn't make sense --if that were the case, wouldn't the donor loci also have the same pattern? How would they be pre-equipped if there wasn't a similar pattern in the donor locus? It seems like the "recipient prior" loci also have the same pattern in somniferum, so wouldn't it seem more parsimonious that the new loci took on the expression patterns of their local genomic environment?

RESPONSE:
We think this comment is probably due to some misunderstandings of our statement "arrived at the recipient locus perhaps pre-equipped with stem-specific promoter or regulatory elements". The statement means the recipient locus carrying some cis-regulatory elements (CREs) that will enable transcription of genes (donor loci) jumping into or around the recipient locus. Therefore, under this premise, the recipient loci prior and post have a similar high and stem-specific expression (Figure 5a, Extended Data Figure 6). We wouldn't expect the donor loci (prior or post) have a similar expression pattern with the recipient loci (Figure 5a, Extended Data Figure 6). This statement is supported by what we observed from the data. We found genes at donor loci (either prior or post) showed low and non-tissues-specific expression pattern (Figure 5a and Extended Data Figure 6), suggesting the lack of stem-specific CREs in donor loci. By contrast, we found genes at recipient loci no matter prior or post showed high and stem-specific expression pattern (Figure 5a and Extended Data Figure 6), suggesting stem-specific CREs were equipped at recipient loci no matter prior or post statues. These results suggested that BIA relate genes were assembled from loci without stem-specific CREs (donor loci) to loci pre-equipped stem-specific CREs (recipient loci). And the pre-equipped stem-specific CREs give the BIA related genes a high and stem-specific gene expression in post recipient loci.
In any case, this is very murky. Li et al. (2020) found that some genes within a cluster had low patterns of co-expression. Do you find evidence that contradicts their claims?

RESPONSE:
In the paper we published on Science (Guo et al. 2018), we showed the 15 morphinan biosynthesis genes within the BIA gene cluster are co-expressed in stem and root. However, BIA gene cluster also contains genes that have no known roles in biosynthetic pathways of BIA. We showed that these non-morphinan genes have low expression patterns (Supplementary Table S15  Why are the other downstream morphine production genes not clustered?

RESPONSE:
We have the same question. Unlike the BIA gene cluster, these downstream genes (CODM, COR, and T6ODM) are located on different chromosomes, where they also went through several tandem duplications. In our Science paper in 2018, we found these downstream genes have high expression in capsule and stem, whereas the BIA gene cluster has strong co-expression in stem and root. Considering these observations, and the fact that the two groups of genes are involved in different stages of the biosynthetic pathway, we speculate that opium poppy has evolved a modular structure of the morphinan biosynthetic pathway with different tissue-specificity of local regulation. This has allowed majority of the pathway genes (leading to thebaine) in a gene cluster, and the rest of the pathway genes (converting thebaine to codeine and morphine) in a different genomic location. By separating the two modules genomically, it may help achieving a highly efficient and finely regulated accumulation of morphinans in the designated tissue types. Another possibility is that the opium poppy is still dynamically evolving as Li et al. (2020) and others pointed out, so that these downstream genes could be recruited into the BIA gene cluster in the future through either natural or artificial selection. , suggesting that the Hi-C loop between morphine branch genes and noscapine branch genes lacks robustness. Therefore, we removed the Figure 5b into Supplementary Figure 32 in the revision manuscript and tune down our conclusion accordingly.
We also did the Hi-C interaction comparison between morphinan gene copies on chr15 and that on chr8 of each P. setigerum replicate (replicate 1 coverage is about 58x, and replicate 2 coverage is about 85x). We found that the Hi-C interactions on chr15 copy was significantly larger than that on chr8 copy on both replicates ( Figure R7), indicating the conclusion about different Hi-C contacts between two copies of morphine branch genes in P. setigerum is robust. Therefore, we added Figure R7 as Supplementary Figure 33 in the revision manuscript.
We agree that more data will be helpful to investigate the intricate chromatin interactions within different Papaver species. This will be our future work.  setigerum. The p-value is calculated by Wilcoxon rank-sum test. e. Hi-C interaction heatmap of regions including two copies of morphinan gene cluster on chr15 (left) and chr8 (right) of P. setigerum merged replicate. The morphinan gene cluster regions are marked as orange boxes. f. The comparison of the interactions from merged replicate between two morphinan gene cluster copies in P. setigerum. The p-value is calculated by Wilcoxon rank-sum test.

It would be useful to have a final summary somewhere of the genes involved in the
morphinan/noscapine clusters and whether they are syntenic, tandem duplicates, or rearranged.

RESPONSE:
Thanks for your comments. Based on the syntenic and BLASTp analysis, we found six genes, including PSSDR1, CYP82X1, CYP719A21, PSMT1, SALSYN, and SALR, were already presented in the MRCA of three Papaver species. During the evolution, PSCXE1, PSAT1, PSMT3, SALAT, and THS were assembled into the gene cluster by dispersed duplications; CYP82X2 was assembled into the gene cluster by tandem duplication; and STORR was assembled into the gene cluster by "fusion, translocation" event. For each gene, we provided the evidence of synteny, phylogeny and WGD from Supplementary Figures 22 -30, and Supplementary Tables 13 -15 in the revision manuscript. The evolutionary model of the BIA gene cluster is summarized in Figure  4. We summarized the evolutionary history of gene in BIA gene cluster in our manuscript from line 221 to line 233 on Page 13-14.
11. The manuscript needs extensive language editing throughout.

RESPONSE:
We have carefully revised our manuscript, and corrected some grammar errors.

Minor comments:
-Fig S14-S15: I don't see the value in synteny painting. It is somewhat interesting that it can be done, but basically all it shows is that there are many rearrangements. This does not lend itself to quantification or hypothesis testing.

RESPONSE:
We agree with your opinion on these two Supplementary figures. We have removed the synteny painting plot in Fig S14-S15, which are updated as Supplementary Figures 15 and 16 in the revision manuscript.
-What tests are being used for the enrichment one line 130? Is this chisq test of all simultaneous? Does it correct for chromosome length?

RESPONSE:
The enrichment test is z-test, which has been corrected by chromosome length. Thanks for pointing this. We have revised this in our manuscript in the legend of Figure 2.
- Fig 1C: What do the colours on the dot plots mean?

RESPONSE:
The color in the dotplot of Figure 1c is generated by MCScanX automatically. Different colors indicate different synteny blocks. We have revised this in our manuscript in the legend of Figure 1.
- Figure 4b: the two colours of blue for "new gene from dispersed duplciation" and "fusion translocation" are easy to confuse. Change the colour scheme.

RESPONSE:
We have revised the color scheme of Figure 4b.
- Figure 5A. It should be stated in the figure that this is for STORR, as it's only clear in the text.

RESPONSE:
We have revised the Figure legend of Figure 5a as "The heatmap of mean normalized expression level of genes at the donor and the recipient loci related with STORR in six tissues of three Papaver species".

Reviewer #2 (Remarks to the Author):
In this study, Yang et al. generated chromosome-level genomes of three different poppy species that display extremes in noscapine and morphinan productions to study the evolution of these biosynthetic pathways. Cutting edge technologies were used to sequence and assemble high-quality sequences. These allow for the investigations of different whole-genome duplication events from around 7 million years ago.

The authors have used different pieces of evidence to support their claims on the ancestral genomes and their evolution, rearrangements and recruitments of new genes to BIA gene cluster locus. Intriguingly, 15 genes involved in the two distinct pathways were assembled into such a compact gene cluster on one single chromosome 4 in opium poppy. Of particular interest is the investigation of the formation of the fusion STORR enzyme through fusion translocation events, possibly thanks to transposable elements. This article could be of great interest to audiences who are interested in natural product pathway evolution.
RESPONSE: Thanks for your positive comments. We are also proud of our work on the genomes of three Papaver species. We applied the latest sequencing technologies to obtain high-quality chromosomal-level Papaver genomes, and found two rounds of whole genome duplications during about 7 million years. We are excited by this work because it is the closest, we have come in the plant biology community to characterizing the evolutionary history of the largest known metabolic gene cluster in plants producing two important medicinal compounds: morphine and noscapine. We believe our work is fundamentally important to understanding trait innovation, pathway evolution and the role of whole genome duplication during speciation and innovation of new traits.

Some minor comments:
-This manuscript lacks an abstract and introduction.

RESPONSE:
We have added the "Abstract" and "Introduction" in the manuscript to indicate these two parts.  Rev Microbiol., 1976), and is a model of pathway evolution, with ancestral enzymes, with broad specificities and catalyzing classes of reaction, forming a large network of possible pathways, and duplication and specialization of these enzymes accounting for extant pathways. We have revised our manuscript to describe the patchwork model at line 386-397 on Page 24.

RESPONSE:
The BIA gene cluster includes the noscapine gene branch and morphinan gene branch. Based on three species comparison, a "fusion, translocation" after WGD-1 but before the divergence of P. somniferum and P. setigerum (between 7.2Ma and 4.9Ma) led to the formation of STORR. And the formation of noscapine branch in P. somniferum was after the divergence of P. somniferum and P. setigerum (4.9 Ma to modern). That means formation of STORR happened prior to the formation of noscapine branch in BIA gene cluster.
There are five genes in the morphinan branch of BIA gene cluster, three of which (SALR, SALSYN, and THS) occurred before STORR formation in BIA gene cluster, while the order of STORR and SALAT remains unknown.

-It is interesting to see that new genes that arrived at the recipient locus might be pre-equipped with the stem-specific promoter or regulatory elements. With the high sequence quality that the authors have, it would be interesting to know what these promoters and regulatory elements would be. Can the authors provide more information regarding this?
RESPONSE: For the promoter sequences, we have added a Supplementary Data file of the promoter sequence of genes (2kb upstream region) in BIA gene cluster in P. somniferum. For the regulatory elements, we have tried to predict the transcription factor (TF) binding motifs in promoter regions of the BIA related genes by FIMO software based on JASPAR database (added as a Supplementary Table 17). However, to validate these motifs and to associate them with potential transcription factors in P. somniferum, a non-model plant, is quite challenging and will need more sequencing data such as Chip-seq for different TF, DNase-seq, DNA methylation sequencing data etc. The ChIP-seq is particularly tricky because the antibodies for specific TFs are unavailable for non-model organisms such as P. somniferum. Therefore, at this moment we could not make any solid conclusions about the regulatory elements for these BIA genes. We plan to combine various epigenomic sequencing, regulatory network inference and functional genomics to pinpoint the underlying regulatory mechanisms of the BIA related genes in our future work.

-
Can the authors elaborate more on the potential epigenetic regulatory mechanism for simultaneously regulating genes from the two branches? It would significantly increase the novelty of the work and would be of great interest for the audience who read this article to have deeper information regarding these mechanisms.

RESPONSE:
Thanks for your suggestion. We agree that more analysis of potential epigenetic regulatory mechanism of BIA related genes is very interesting. However, as our response to the last comment, it will require more epigenomic sequencing data, regulatory network inference and functional genomic analysis to pinpoint the epigenetic regulatory mechanisms for the BIA genes. For non-model organisms such as Papaver species, this is a challenging task and will require careful investigation in the future. In this work, we focus on comparative genome analysis of three Papaver species, and the epigenetic and transcriptional regulation of BIA genes will be our future work.

-
What are the levels of identity of the genes from these different species?

RESPONSE:
In Figure 4a and Supplementary Figures 22-30 in the revision manuscript, we labeled the identity as numbers of each gene compare to BIA related genes in P. somniferum. We made a table and a heatmap to show the identity distribution ( Figure R4, Table R1), and the dominate identity level is larger than 91%.  The genes involved in a pathway are not always clustered in a gene cluster, especially for Eukaryotes, e.g. plants. Although, many gene clusters have been reported to encode secondary metabolic pathways in plants, including opium poppy, none of them has been linked to the patchwork model. To our knowledge, our work on BIA gene cluster of Papaver genus revealed the first gene cluster explained by patchwork model in plants.
We have added the references and the discussions of patchwork related work in our manuscript at line 391-397 on Page 24.

Reviewers comment for Nature communication article-
In the present study under review, titled, "Three chromosome-scale Papaver genomes reveal punctuated patchwork pathway evolution leading to morphinan and noscapine biosynthesis", authors reported whole-genome assembly for two Papaver species, namely P. setigerum and P. rhoeas, and worked to improve previously published genome for P. somniferum. These three species represent extreme in the sense of morphine and noscapine synthesis. The authors then used these three genome assemblies to reconstruct the ancestral genome for Papaver species and derived a hypothesis for the evolution of BIAs in plants. The authors proposed the role of whole-genome duplication, followed by chromosomal rearrangement, fusion, and translocation as key events towards the innovation of secondary metabolites. The study is well written, and all figures are of good quality. I also need to mention here that overall, the manuscript is super-condensed and sometimes hard to follow, especially with the transition from ancestral genome construction towards the proposed hypothesis of the role of fusion and translocation towards the evolution of specialized metabolites. I highly recommend to consider expanding the explanation, simplification of text, and improving the flow if the argument with more detail in order to attract general and specialized readers towards this study. One of the weaknesses of this study is the generality, which is entirely missing. The offered method for ancestral genome construction, the whole hypothesis of evolution of specialized metabolites, and all discussions are Papaver specific. I feel that the authors have not attempted to draw comparisons with other studies or to draw a general point-of-view based on their results and other high-quality genomes that have also offered different approaches towards the evolution of specialized metabolites. This is one of the significant weaknesses of this study that I feel must be taken care of. While going through this study, I find several questions unanswered. I am listing these in the same flow as these are described in the manuscript.

RESPONSE:
Thanks for your comments. We intended to write in a way that keeps the main manuscript concise and provides abundant supporting evidence in the supplementary materials. In this way, the readers may grasp the main conclusions of the main text, while those who are interested in more technical details are directed to read the supplementary materials. We have revised our manuscript by adding some details.
As for the comment "the transition from ancestral genome construction towards the proposed hypothesis of the role of fusion and translocation towards the evolution of specialized metabolites", we added a few sentences for the transition from ancestral genome reconstruction to the STORR formation, and the transition from STORR formation to the BIA gene cluster evolution.
As for the comment "One of the weaknesses of this study is the generality, which is entirely missing", we included examples of the patchwork model, the biosynthesis pathway, and the gene clusters in bacterial and plants in the discussion section.
We have revised the manuscript as you suggested, and answered the questions point-to-point. . Given that different cultivars have various amount of indels and structural variations, which would cause artifacts and errors to genome rescaffolding in one cultivar using sequencing data from another cultivar. Therefore, for the best result of genome assembly, in this manuscript we generated Hi-C sequencing data of cultivar HN1 to improve the draft HN1 genome we published in 2018.
Secondly, we have doubts over the genome assembly statistics such as contig N50 reported in Li et al. paper. According to the methods in Li et al. 2020, their strategy for genome assembly is as the following "The sequences of the Guo assembly were broken into segments at gaps of 100 Ns, representing the inter-scaffold gaps of unknown size, and at large gaps of ≥1000Ns. The resulting split genome was comprised of 35,732 resulting segments, which were considered as contigs in the subsequent processing. The lengths of these contigs varied from 132 bp to 38.3 Mb, with the N50 length and the N50 number of the contigs of 7.6 Mb and 104, which provides a draft assembly with sufficient contiguity for making high confidence Hi-C based proximity-guided rescaffolding". What it says is that Li. et al did not assemble genome contigs from scratch, but simply took our published contigs to perform a Hi-C scaffolding. Therefore, the report that they have managed to produce a genome assembly with contig N50 improved from 1.74 Mb (Guo et al., Science, 2018) to a reported 7.6 Mb sounds incredible. In fact, we downloaded Li et al. genome assembly from NCBI genome repository and the global statistics demonstrates their contig N50 is 1.839 Mb (https://www.ncbi.nlm.nih.gov/assembly/GCA_010119995.1). We also recalculated the contig N50 by bbmap stats.sh, which showed that the contig N50 is 1.839 Mb, slightly larger than the 1.74Mb we reported in 2018, but far smaller than 7.6Mb Li et al. has reported. Actually, this slightly increased contig N50 from 1.74Mb could be down to the fact that some small contigs were removed (the number of contigs is reduced from 66,578 (Guo, et al. Science, 2018) to 9,646 (Li et al. Nature communications, 2020)), and which was also reflected from genome size reduction from 2.71Gb to 2.637Gb ( Figure R8).
Given these facts and observations, we opted to generate an improved P. somniferum HN1 genome assembly by ourselves and used it for genome comparisons of three Papaver species.  Figure 2). As the entire manuscript premise is about evolution of the specialized metabolites biosynthesis using these three plants, I feel that looking at overall metabolites and metabolite intermediates to morphine, noscapine across these species will be important. I am particularly curious as how about the intermediates from the earlier steps of biosynthesis pathways of key metabolites and their levels across these three species. Probably such quantification would allow to associate evolution of genes with presence/absence of the metabolites. If these metabolites/intermediates have already been reported and described for these plant species in the past, would be nicer to include as part of discussion and authors views on association of presence/absence of genes with the metabolites. If metabolites are not identified yet, authors should consider identifying and quantifying these intermediates.

RESPONSE:
We agree that it will be ideal to compare the three species for producing all the various intermediate metabolites. However, this is quite a daunting task and well beyond the objective and scope of this current manuscript. The most challenging part of it is to find trustworthy analytical standards for most of these intermediates. We are not a chemistry research group, and cannot syntheses these intermediate on our own. Due to the regulated substances and the influence of the COVID-19 pandemic, we can only obtain all available standards, i.e. Thebaine, Codeine, Morphine, and Noscapine, from domestic market. Therefore, we used HPLC-MS to quantify those four morphinans in the three species.

Regarding the question of "If these metabolites/intermediates have already been reported and described for these plant species in the past, would be nicer to include as part of discussion and authors views on association of presence/absence of genes with the metabolites. If metabolites are
not identified yet, authors should consider identifying and quantifying these intermediates.", yes, some metabolic profiling has been reported in literature for several Papaver species including P. rhoeas and P. setigerum. However, we focused on BIA pathways as these are the best understood pathways, for which our genomic analysis can reveal the genetic differences behind these metabolites. For other secondary metabolites, it is challenging to perform similar analysis or discussion as the genetic basis of their biosynthesis is largely unknown.

RESPONSE:
You are right that normally it is more difficult to achieve good contiguity for large and complex genomes than for small and simple genomes. The contiguity of genome assembly is usually affected by multiple factors, such as heterozygosity rates, polyploidy, raw reads quality, repeat content in genome etc. In our study, raw reads quality and genome heterozygosity are two main reasons on the rates of assignment to the scaffolds. For three Papaver species, we have high quality raw reads (Supplementary Table 1), e.g. the ONT raw reads N50 are about 30Kb in both P. setigerum and P. rhoeas, the Q30 of Hi-C data for P. setigerum, P. rhoeas, and P. somniferum are larger than 91%. The main difference of three Papaver species is the heterozygosity rate. P.
setigerum, despite its large genome size, has a low heterozygosity rate, as shown by a lack of clear heterozygosity peak in the k-mer frequency distribution of P. setigerum sequencing reads (Supplementary Figure 3). By contrast, despite the relatively smaller genome size, P. rhoeas has a high heterozygosity rate of 3.18% as shown by a clear heterozygosity peak in the k-mer frequency distribution of P. rhoeas sequencing reads (Supplementary Figure 5). Therefore, although the read length, quality of sequencing data and assembly methods for both genomes are similar, the genome assembly contiguity differed a lot. We have revised the supplementary method (Section 6.3 Assembly evaluation) by adding the discussion of the reason for different assignment rates in three Papaver species. Table 1, authors reported number of protein coding genes across these three species. The numbers are quite comparable between P. somniferum and P. rhoeas, while the former have undergone whole genome duplication while the latter has not. While we know that whole genome duplication is further followed by massive genome readjustment, including deletion and rejection of gene families and so on, I am curious to know as what genes specifically were gained and lost across these species. Are they somehow being favored for BIA biosynthesis? What kind of processes were lost in P. somniferum? Authors have briefly described it in section 10 of supplementary information, but I feel that a direct comparison between these three species in terms of understanding gene gain and loss would be interesting as they do have a contrasting levels of BIAs despite being closely related. I also feel that including or expanding this aspect in the discussion will make the discussion section better to understand evolution of BIAs as emerging through this study.

RESPONSE:
Thanks for your comments. We performed the syntenic analysis of three Papaver species, and found 28,660 genes in P. somniferum were syntenic with 19,512 genes in P. rhoeas with syntenic depth from one to five ( Figure R9a), indicating that 28,660 genes are kept in P.
somniferum following WGD-1 and diploidization. For any two-species comparison, it is difficult to differentiate gene "gain" and "loss" because gain for one species means loss for the other species, and vice versa. Alternatively, we found 21,958 and 26,654 genes are specific to P. rhoeas and P. somniferum respectively, by comparing P. rhoeas and P. somniferum genes. We performed the functional enrichment for species-specific genes to understand their functional roles. Based on the functional enrichment analysis, and the P. somniferum specific genes were significantly enriched in energy, photosynthesis, and metabolism related pathways, while P. rhoeas specific genes were significantly enriched in oxidative phosphorylation, ubiquitin system, ABC transporters related pathways ( Figure R9c).
Similarly, we compared P. somniferum with P. setigerum, and found 41,073 genes in P.
somniferum were syntenic with 71,398 genes in P. setigerum with synteny depth from one to 11 ( Figure R9b), indicating that 71,398 genes in P. setigerum were related with WGD-2, while 14,241 genes in P. somniferum were specific and 35,119 genes in P. setigerum were specific based on the comparison between P. somniferum and P. setigerum. The functional enrichment showed the P.
somniferum specific genes were significantly enriched in photosynthesis, ribosome, metabolism related pathways, while P. setigerum specific genes were significantly enriched in Spliceosome, metabolism, Endocytosis related pathways ( Figure R9d).
We have added a section "8.4 Protein coding gene number comparison based on synteny analysis" the Supplementary material and added Figure R9 as Supplementary Figure 14 in the revision manuscript to explain the difference of protein coding genes in three Papaver species. Figure R9. The synteny depth of P. rhoeas versus P. somniferum (a) and P. somniferum versus P. setigerum (b). c. The pathway enrichment of P. rhoeas specific genes and P. somniferum specific genes based on the comparison between P. rhoeas and P. somniferum. d. The pathway enrichment of P. setigerum specific genes and P. somniferum specific genes based on the comparison between P. setigerum and P. somniferum. We selected the top20 significantly enriched pathways in this figure.

Authors mentioned in the method and supplementary method section that they used purge_dup
to exclude any redundancies from the genome assembly. I think that its important to provide parameters used. As purge_dup may also exclude some of the genes that may be real and not artifacts, the used cutoff threshold needs to be mentioned. It may be a good idea if authors could provide all scripts and parameters for different stages of analysis/assembly/ comparative genomics as a github repository as they have provided for reconstruction of ancestral genome.

RESPONSE:
Thanks for your comments. In our assembly pipeline, we use the default parameters. For purge_dup used in P. rhoeas genome, the cutoffs parameter is automatically calculated by calculate module as "5 34 56 67 112 201". We totally agree your suggestions to provide all scripts and parameters at GitHub. We have uploaded the scripts and parameters used at https://github.com/xjtu-omics/Papaver-Genomics/tree/main/analysis_scripts. 6. This is just a minor comment. For phylogenetic tree, authors have used eight angiosperms. I feel that including more plant species for the phylogenetic tree by including plants from different lineages and from different time of evolution would be nicer to provide a relative view of the evolution of Papaver species with others. This is simply my preference, and I leave this up to the authors to agree or disagree.

RESPONSE:
Thanks for your kind suggestion. We have carefully selected plant species to compare with Papaver species in our phylogenetic analysis. In our opinion, it is critical to include plant species from lineages that are close to Papaver spp. in order to get an evolutionary picture of high resolution. However, available genome assemblies of high quality within Papaveraceae are scarce so far, which has limited our options to perform extensive intrafamily phylogeny analysis. Therefore, we included the Macleaya cordata (diverged from Papaver spp. around 60Ma), the first plant in Papaveraceae that was sequenced and assembled. On the other hand, we chose eight angiosperms from representative families that are commonly used in plant phylogenomic analysis, including monocots (rice), core eudicots (grape, Arabidopsis etc.) and an early-diverging dicot (Aquilegia). Sampling many species distant from Papaveraceae will probably not reveal any new insights into Papaveraceae evolution because such species have diverged a long time ago from Papaver spp. In addition, our selected angiosperms diverged from Papaver species at times ranging from 152Ma to 60Ma, covering a wide scale of evolutionary time. Accordingly, we believe our phylogenomic analysis (Figure 1e) provides a fairly good view of the evolutionary relations of Papaver species with others.  Figure 1d, we show the synonymous substitution rate (Ks) density distributions of syntenic paralogs and orthologs. It should be noticed that we split WGD-1 and WGD-2 in P. setigerum since the WGD-1 (at 7.2Ma, Ks peak value = 0.115) is very close to WGD-2 (at 4.0Ma, Ks peak value = 0.065). The two close WGD events are reflected by the mixture of two peaks ( Figure R10). In our work, we considered that only the reciprocal best matches among the syntenic gene pairs were considered as the pairs from WGD-2 while other syntenic gene pairs were grouped as the pairs from WGD-1, following the procedure introduced by Liu et al, Cell, 2020. For Arabidopsis, the WGD-1 (at 67Ma) is far away from WGD-2 (at 50Ma), which was reflected by two distinct peaks. Figure R10. The Ks distribution of P. setigerum without splitting WGD-1 and WGD-2

Much of this article, including hypothesis and interpretations for comparative genome analysis replies on the constructed ancestral genome of Papaver species.
RESPONSE: Thanks for your comments. We decoded the evolutionary history of BIA gene cluster based on the systematic gene syntenic analysis of three Papaver genomes, rather than the reconstructed ancestral genomes. Reconstruction of ancestral genome is the way to decode the karyotyping of Papaver species based on the syntenic blocks, while investigating the evolutionary history of BIA gene cluster require the syntenic gene pairs among three Papaver species. Therefore, the investigation of the evolutionary history of BIA gene cluster is based on the gene syntenic analysis among three modern Papaver species. Sorry for the unclear transition between "Ancestral genomes and accelerated nonrandom post-WGD rearrangements" section to "Recruitment of new genes to BIA gene cluster locus " section. We revised the manuscript at line 173-176 on Page 11 to make the transition clear.

At the beginning of section, "Ancestral genomes…", author says, "We developed a novel bottom-up framework to reconstruct pre-and post-WGD-1 ancestors of Papaver genomes". I have followed work mainly from Jerome Salse, and I wonder how this method is different from what his group have been using and have proposed before? A comparison in terms of what difference between this approach and what previously have been used and providing further discussion on advantages and novelty would be very important.
We have studied Jerome Salse's research (Murat et al, Nat Genet., 2017) in details, and found the methods is not able to handle the ancestral genome construction of three Papaver species for the following reasons:  In Jerome Salse's research, they used three representative species, including grape, cacao and peach, to build ancestral eudicot karyotype (AEK). These three species shared a whole-genome-triplication event (the γ event) at around 100 Ma and never happened another WGD event after that. Since the long time of γ event, the block ratio of three species become 1:1:1 which is suitable for modeling as GMP and solving by MGRA.
 However, the evolutionary scenario of Papaver species with a shared WGD by P. somniferum and P. setigerum at around 7.2Ma and a P. setigerum lineage specific WGD at around 4.0Ma. Such a close time scale of these two WGD makes the block ratio as 4:2:1.
To solve this problem, we developed our method for three Papaver species with two closed WGD events based on the work of Sankoff (Zheng et al, Evol Bioinform. 2007) and Pedro Feijao (Feijão et al, TCBB, 2011). Compared with previous work, we attempted to use block matching strategy to match block copies in related species first by minimizing single cut or join (SCoJ) distance (Feijão et al, TCBB, 2011) and then relabeled block copies to transform problems in Papaver species to traditional GMP and GGHP (Guided Genome Halving Problem). Moreover, we solved the ancestral genomes by integer programming frameworks.
We have revised our method section on reconstruction of ancestral genome in manuscript (Methods Section Ancestral genome reconstruction) and Supplementary Materials by adding the comparison as you suggested.

Authors provided github repository for the codes used in this study, and they say in the text as well as in the description of github repository that this code is suitable for reconstruction of ancestral genome of Papaver species. I wonder if the codes could be used for other plant species. Does this code be applicable for including and comparing different plant species and to derive ancestral genomes for a certain lineage? Authors need to clarify this in the text, and
if this is specific to Papaver species, why? This is a topic that would be of wider audience interest, and if addressed would be able to make this paper even stronger and would add general appeal.

RESPONSE:
Thanks for your comments. In this research, we designed and implemented the ancestral genome reconstruction according to the evolutionary scenarios with two recent WGDs, one shared and one lineage specific, of three Papaver species in this study. Our pipeline can also be applied to similar evolutionary scenarios. Besides the GMP and Papaver evolutionary scenarios, we are exploring the solutions for other complex evolutionary scenarios and eventually we will provide a generalized framework for all possible evolutionary scenarios. We have revised the method description in our manuscript and updated the readme at GitHub. , which means the ancestral genome is right before the WGD event, and we reconstructed the pre-WGD-1 ancestor at 7.2Ma. We compared pre-WGD-1 ancestor with AEK, and create a dotplot to show the comparison result ( Figure R11). We added the Figure R11 as the Supplementary Figure 16g, h in the revision manuscript.

RESPONSE:
We agree that assembly errors will affect the calculation of chromosomal fissions and fusions. In our work, we assembled three Papaver species based on the cutting-edge sequencing data and the wildly used assembly technologies. We evaluated our assembly from both completeness and base accuracy (Supplementary Section 6.3), and the results indicated the high-quality chromosome level genome assemblies of three Papaver species. We are confident of our data, methods, analysis, and the conclusions about the "non random distribution of chromosomal fissions and fusions" and the "superliner correlation between number of WGDs and number of chromosomal rearrangements". However, we admit that even with the cutting-edge sequencing data and widely used assembly methods, assembly errors are inevitable, especially for genomic regions with large and repetitive segments. Thus, we tune down the conclusion in our result (Section Ancestral genomes and accelerated nonrandom post-WGD rearrangements). If the editor and the Reviewer #3 insist that the genome assembly errors may compromise the conclusions about the number of chromosomal fissions and fusions, we are willing to remove the content related to chromosomal fissions and fusions.  Figure R6, Supplementary Table 18 in the revision manuscript), suggesting that the Hi-C loop between morphine branch genes and noscapine branch genes lacks robustness. Therefore, we removed the Figure 5b into Supplementary Figure 32 in the revision manuscript and tune down our conclusion accordingly.

Authors have very elegantly described the origin of STORR in the section
We also did the Hi-C interaction comparison between morphinan gene copies on chr15 and that on chr8 of each P. setigerum replicate (replicate 1 coverage is about 58x, and replicate 2 coverage is about 85x). We found that the Hi-C interactions on chr15 copy was significantly larger than that on chr8 copy on both replicates ( Figure R7), indicating the conclusion about different Hi-C contacts between two copies of morphine branch genes in P. setigerum is robust. Therefore, we added Figure R7 as Supplementary Figure 33 in the revision manuscript.

Since expression of genes involved in the BIAs biosynthesis are known to be tissue species,
and this is what authors also reported by saying that the epigenetic factor played a role in getting these genes highly expressed in a tissue specific manner. Do authors have any reference to support their point other than the comparative HiC data for these three species? What I mean is that it would be clearer if authors could show Hi-C data for tissues know to have no expression of BIA biosynthesis associated genes and metabolites, and tissues that show highest accumulation of BIAs and expression of genes. Those comparisons will be better to identify interactions at chromatin level (even if this is done for P. somniferum), which then could be explored across these three species to derive the conclusion that authors have made here.

RESPONSE:
Thanks for your comments. Since Papaver species are not model organism, epigenetic datasets, such as DNA methylation, histone modification, and tissue-specific Hi-C data, are still unavailable. In our work, we generated the Hi-C data from leaf successfully. We attempted several times to construct the Hi-C libraries of P. somniferum roots and stems. However, libraries construction failed due to the high lignin content in these two tissues. Therefore, we can only analyze leaf Hi-C data in this study, which is rather common in plant Hi-C studies. We are still attempting to construct the root and stem Hi-C libraries. In addition, we are planning to investigate epigenetic factors, such as DNA methylation, histone modifications, to decode the underlying regulatory mechanisms of BIA genes in the future. performed intra-species comparison by including resequencing data for nine P. somniferum cultivars and indicated the copy number variants of BIA genes are positively correlated with morphinan accumulation in these cultivars. By contrast, we did comparison among three Papaver species and revealed that novel insights into how BIA gene cluster was formed and evolved among different species as well as accumulated gene expression besides copy numbers contributed to the morphinan production in Papaver species. We have revised our discussion section by adding the comparison between our work and previously published work on P. somniferum on Page 22-23.
In addition, we revised the Results section (on Page 14) and Discussion section (on Page 23) to illustrate the evolutionary scenarios of available metabolic pathways, and other genes in the BIA gene cluster on Page 23. the BIA gene cluster among Papaver species, we de novo assembled two chromosome-level genomes of P. rhoeas (common poppy) and P. setigerum (Troy poppy) while also improved the previous draft P. somniferum HN1 genome.".
Line 52: Chapter "Genome assembly and annotation" -describing the pattern of morphine and noscapine accumulation in the model species P. somniferum, P. rhoeas and P. setigerum better fits to the introduction (see above). Also, check the sentence and cite references! And the cross reference to Fig. 1a is a bit misleading at this point (line 55), as the Fig. 1. does not illustrate alkaloid accumulation. Which, by the way, could be helpful. For example, the "alkaloid-type" could be included in Fig. 1e.

RESPONSE:
Although it is known that the production levels of BIAs in three Papaver species are different, the precise amount was determined in this study (Supplementary Figure 2). Thus, we decided to describe the patterns of morphinan and noscapine accumulation in three Papaver species in Results rather than Introduction.
Based on the divergence time and Ks distributions, we estimated the WGD time ( Supplementary Section 8.3). The method is like "Given the mean Ks value (0.12) of P. somniferum-P. rhoeas and their divergence date T (7.7 Ma), we calculated the synonymous substitutions per site per year (r) for Papaveraceae as 8.08e-9 (T = Ks / 2r) (Supplementary Table  6) which was applied to time the WGDs of P. somniferum and P. setigerum. We dated the WGD in P. somniferum (Ks = 0.116 ± 0.028) around 7.2 ± 1.7Ma, the first WGD in P. setigerum (Ks = 0.115 ± 0.018) around 7.1 ± 1.1 Ma, and the second WGD in P. setigerum (Ks = 0.065 ± 0.017) around 4.0 ± 1.0 Ma (Fig. 1d, Supplementary Table 6)." Page 10 -associated to Fig.2: In a pervious study ("The opium poppy genome and morphinan production"), the BIA gene cluster was described to be located on Chromosome 11. As far as I understood, the sequence data from this former analysis was re-analyzed in the actual study. Please elaborate why chromosome of BIA gene cluster changed. Have there been large re-assignments? RESPONSE: Thanks for your comments. After analysis of the Hi-C data of P. somniferum, we found our previous published genome assembly (Guo, et al, Science, 2018) has several mis-scaffolding ( Figure R12a). We used Hi-C data and 3d-DNA software to rescaffold our published contigs to obtain current assembly. We compared the rescaffolded genome assembly to the published genome assembly, and reassigned chromosome IDs based on the genome dotplot ( Figure R12b). After reassignment of chromosome IDs, we found the BIA gene cluster were located at chr4 rather than chr11. Based on the genome wide Hi-C contact maps, we can see that the rescaffolded genome assembly has indeed been improved by correcting the mis-assembled contigs in the previous published genome.  Fig. 4a, it seems to me, that Fig.4b is correct. Elaborate on the "dispersed" duplications. Are these small-scale duplications? Not result of WGD? How did they "arrive" in the gene cluster?

RESPONSE:
We have revised that reference of Fig. 4b. Dispersed duplication means the duplicated genes and their original copies are not adjacent. Yes, the dispersed duplications are small-scale duplications compared to WGD and not result of WGD. We tried to investigate the mechanisms of "how did they 'arrive' in the gene cluster", and labeled the annotated transposon elements (TEs) in related with these genes in Supplementary Figures 22-30 Fig.1c: Please include color code/legend that was used to color the syntenic blocks. I assume its Ks?

RESPONSE:
The color in the dotplot of Figure 1c is generated by MCScanX automatically. Different colors indicate different synteny blocks.     Fig. 4: Please include color code/legend that was used to color the syntenic block. And could you elaborate it bit on the findings, that we can get from the dotplots in the figure legend? E.g. in 4b -P. somniferum plotted against the post WGD-1 ancestor -chromosome 1 shows a syntenic block ration of 1:2, or in other words: a syntenic blocks in chromosome 1 P. somniferum "matches" chromosome 1 and chromosome 6 in the post-WGD-1 ancestor. What is the explanation? Small-scale duplication? RESPONSE: These dotplots were generated by MCScanX, and the color was generated by MCScanX automatically. Different colors indicate different synteny blocks. We have revised the figure legend to indicate this. Comparisons between three Papaver genomes and two reconstructed ancestors (pre-WGD-1 ancestor and post-WGD-1 ancestor) reveal the differences between the ancestor genomes and the modern Papaver genomes. For example, the comparison between post-WGD-1 ancestor and P. somniferum show a 2:2 block ratio, indicating both P.

Extended Data
somniferum and post-WGD-1 ancestor have WGD event. A syntenic block in chromosome 1 of P. somniferum "matches" both chromosome 1 and chromosome 6 of the post-WGD-1 ancestor indicating chromosome 1 and chromosome 6 in post-WGD-1 ancestor were a duplicated pair resulted from WGD-1 rather than small-scale duplication.
The authors have made a revision that addresses some of the issues raised in the last review. The new figures in the supplementary materials now provide useful information that goes part of the way to assessing whether tandem vs. dispersed duplication is responsible for cluster evolution, but the study really has not clearly demonstrated this VERY important point. It is not suitable for publication until this aspect is clarified and actually tested rigorously (comment #2). In several places, I still feel that the manuscript is putting forward one possible explanation for how things evolved, but not really doing a good job to test among different alternative explanations, as described in the comments below. I think this paper has a lot of REALLY interesting data --I really like the comparison with the other two poppy species --but it is still being hampered by these aspects to the point where it is not yet suitable for publication.
LINE NUMBERS refer to the line numbers in the word document included as track changes.
Major comments: 1. The writing needs a lot of work. For example, the first few sentences of the abstract are unclear and halting in their style and not precise in their meaning: -line 22: selection operates on things, not selection of, and this implies that evolution only operates through selection on this variation, which is not true/general (selection also acts on SNPs, indels, rearrangements, etc). -line 23: what is "single genomic variation"? This is not standard terminology -line 25: what does "forwardly or backwardly" mean? What is the patchwork model? This is not selfevident and everything in the abstract should be crystal clear. These problems extend throughout the paper and go far beyond what I can reasonably edit as a reviewer --the science is interest but it is not presented at the standard of a journal like Nature Communications, and this is not just about whether the language is polished but also whether it is clear in its presentation of ideas. This should not be a barrier to publication when the science is sound, but it is important that it is rendered in the proper style. As an example of how this hampers communication --the introduction is too short and fails to review the main ideas being discussed. As an example, the patchwork model is only introduced on line 605 and then it simply uses the name of the model as the definition, which is not clear and this is not something well-known to a general audience. One line 608 it is said that the BIA cluster is the first example in plants of the patchwork model, but what about the extensive work done by Nutzman and Osborn (http://dx.doi.org/10.1016/j.copbio.2013.10.009) on clustering plant secondary metabolic pathways? Why does previous work not conform? I don't know what the patchwork model is, and I work in this general area of the field. This is just one example but in general I find the scholarship aspect of this paper (clarity of terminology, clarity of concepts, linking with existing literature) is still lacking and this is an important part of academic publication.
2. While the authors have responded to my comments about STORR, they haven't really changed the presentation in a substantial way. I don't think it will be evident to any reader why this shows evidence for STORR evolution as the authors explain. Figure 3 is a cartoon showing a hypothesis, it doesn't present any evidence. Figure S20 shows evidence, but it is so convoluted it is hard to see anything besides "it's complicated". I can accept the author's story about how they think STORR evolved as one possible explanation, but Figure S20 doesn't really make that particularly clear, and other less-parsimonious explanations are possible (parsimony isn't always == truth when it comes to phylogenetics).
More importantly, the authors have largely ignored suggestions I made to improve Figure 4 and to more clearly identify whether genes arose from tandem duplication or dispersed duplication. I think it is REALLY important to nail this down because tandem duplication is very common and immediately leads the the evolution of a cluster, whereas the chance of a dispersed duplicate fortuitously landing right next to other genes to form a cluster is small ---so seeing evidence of this is very interesting but requires a high burden of proof. So the burden of proof here lies on showing that a gene is a dispersed duplicate (tandem duplicate is the null expectation). While I like this figure 4 as an illustration, I don't think this qualifies as comprehensive evidence about whether these genes arose from old tandem duplications vs. "dispersed duplication" --the blastp results added to Figure 4 don't really allow us to assess whether a gene arose from a tandem vs. dispersed duplication --they just show one line of "best blast hit", but this best blast hit could just be due to the WGD if that happened after an older tandem duplication. This can be assessed using Figure S31 but it's still hard to disentangle which copies are homologous because of the WGD vs. because of a dispersed duplication. I REALLY think the authors need to do a better job of discriminating between these hypotheses. As an example, consider PSMT3. The authors present a cartoon showing one possible way it evolved in Figure S28. However, examination of Figure S31 reveals other explanations could also be consistent with the data. The BIG problem to consider here is that PSMT3 could have evolved as an old tandem duplication of PSMT1 (or vice versa), long before the WGD. Figure 31 clearly shows there is homology between PSMT3 and PSMT1, but less than within the copies of PSMT3. Within the cluster labelled PSMT3, there are 3 Psom genes: -Pso04G00300 -PSo01G11570 -PSo05G43680 Within the cluster labelled PSMT1, there is only one Psom gene: -Pso04G00330 Thus, a hypothesis consistent with the data is that either PSMT1 or PSMT3 arose from a tandemduplication (of the other) prior to the WGD, then after the WGD, the copy of PSMT1 was deleted from the other WGD segment, while another dispersed duplication of PSMT3 gave rise to the other copy (either Pso01G11670 or Pso05G43680, whichever is not part of the WGD) but not in any cluster. This seems more parsimonious than assuming that PSMT3 arose from a dispersed duplication and just happened to land right next to PSMT1, in the proper cluster. This simply cannot be evaluated using the gene tree and dot plot shown in Figure S28. We need to see analysis of WGD regions vs. putative dispersed duplicates vs. putative tandem duplicates (assessing all nearby genes, not just those that are immediate neighbours --tandem duplication can duplicate several genes at once, that can subsequently be deleted).
Similar arguments apply for the other genes where the authors are claiming that a dispersed duplicate occurred and landed within the cluster (PSAT1, PSCXE1, THS, and SALAT). In each case, the authors need to conclusively disprove the hypothesis that an old tandem duplication gave rise to the putative "dispersed duplicate" by evaluating any cases where other genes nearby on the same chromosome have detectable/substantial %ID, as this may indicate an old tandem duplication. I think in many cases this simply won't be possible to confidently reconstruct what happened, but that's ok --better to be clear about the uncertainty than be wrong. This is vitally important as it helps us understand how these clusters formed and the current analysis is simply not up to a high standard.
3. Given the ambiguities in how these genes have evolved (comment #2), I have doubts about the validity of "donor" and "recipient" locus definitions. This affects all usages of these terms (e.g. lines 363-393). For example, if PSMT1/PSMT3 arose from a tandem duplication, then the labelling of other genes as "donor" would be wrong. 4. It is unclear why it is surprising that morphinan production is lower in setigerum than somniferum, given that somniferum has been actively under artificial selection for this trait. Also, the previous section just (nicely) showed that the promoter regions of one of the two clusters in setigerum have degraded, so it's even less surprising that setigerum has lower production. Unless I'm mistaken, the Li (2020) paper didn't show extensive correlation between copy number and production, they showed that deletions of particular genes corresponded to loss of alkaloid production, but not that increases in copy number were explanatory (?). Finding that setigerum has "twice as many copies of STORR..." (line 512) doesn't seem relevant after the previous section showing lower expression, but I see the point that is being made from line 520-532. 5. It is claimed that the role of transposable elements in gene duplication was investigated, but nothing conclusive is shown here, all that is shown is that there are TEs near these genes, this does little to nothing to clarify how these genes evolved. Claims must stand on evidence, and again, it feels like the authors are telling a story that goes along with our current understanding of how genomes evolve, rather than testing whether this is in fact what happened.
Minor comments: -line 260-267: the concepts of donor and recipient loci are not clear as described.
-line 354: add reference to other places this has been reported -line 568: this is not an example of punctuated equilibrium. Please reconsider the wording here.
-line 350: repeating a previous comment from my last review (!): please change all usage of "concerted" as this term already has another meaning (see https://en.wikipedia.org/wiki/Concerted_evolution). There are many other words that are more suitable, such as correlated or coordinated.
-line 528: seeing a difference in copy number of genes is not evidence of artificial selection. Again, this is another case of telling a story using data that is consistent with what we already know (somniferum experienced artificial selection), rather than using data to test among alternative hypotheses.
Reviewer #2 (Remarks to the Author): The authors have explained and addressed the issues and questions raised by the reviewer.
Reviewer #3 (Remarks to the Author): In the revision of the manuscript, titled, "Three chromosome-scale Papaver 1 genomes reveal punctuated patchwork pathway evolution leading to morphinan and noscapine biosynthesis", authors have considerably improved the content, especially by including content that appeals to a broader audience. The authors managed to address almost all my comments. The authors have now provided all scripts, including parameters that were used for genome assembly and analysis. After going through the responses to my queries, there are few points that remain unresolved, and I am listing them here.
 We tuned down the duplication types as "putative dispersed duplication" and "putative tandem duplication" in manuscript and revised the Figure 4. For your comments on STORR, we are delighted about your acceptance of our interpretation. In addition, reviewer #3 also appreciates our evolutionary model of STORR by saying "Authors have a very elegantly description of the origin of STORR" (see the comments #12 of reviewer #3 in previous review comments). Moreover, we think that our parsimonious explanation of STORR formation based on current three Papaver genomes is one of the many possible explanations. We agree with you that "parsimony isn't always == truth when it comes to phylogenetics", so we tuned down our conclusions about STORR evolutionary model and revised the results as well as discussion sections.
For comments on Figure 4, in our last response to your previous comments, we thought that we have addressed your previous comments as follows. Would you please let us know what additional modifications are needed?
 We added the Supplementary Figure 31 to show the identity of all genes in noscapine and morphinan clusters as well as their close paralogs.
 We added the detail explanation of the evidence to infer the ancestral origin and duplication type of each novel gene in BIA gene cluster.
 For "This would be SO much clearer if the paralogs to the pre-fusion modules within somniferum were shown on Figure 4. Where are they?" from your last comment about Figure 4, we have shown the paralogs of pre-fusion modules in three Papaver species in Figure 3a and Supplementary Figure 20, and these two figures showed the evidence of our STORR story. In our manuscript, Figure 3 is the story of STORR and Figure 4 is the story of whole BIA cluster evolution. The flow is from one gene to whole cluster, and we think this flow is smooth.
For comments on other genes (PSMT3, THS, SALAT, PSAT1, PSCXE1), thanks for the alternative explanation you pointed out. We revised the manuscript to include the alternative explanation and added Supplementary Figure 32 (Figure R1). 3. Given the ambiguities in how these genes have evolved (comment #2), I have doubts about the validity of "donor" and "recipient" locus definitions. This affects all usages of these terms (e.g. lines 363-393). For example, if PSMT1/PSMT3 arose from a tandem duplication, then the labelling of other genes as "donor" would be wrong.