Genome-wide analysis of Cushion willow provides insights into alpine plant divergence in a biodiversity hotspot

The Hengduan Mountains (HDM) biodiversity hotspot exhibits exceptional alpine plant diversity. Here, we investigate factors driving intraspecific divergence within a HDM alpine species Salix brachista (Cushion willow), a common component of subnival assemblages. We produce a high-quality genome assembly for this species and characterize its genetic diversity, population structure and pattern of evolution by resequencing individuals collected across its distribution. We detect population divergence that has been shaped by a landscape of isolated sky island-like habitats displaying strong environmental heterogeneity across elevational gradients, combined with population size fluctuations that have occurred since approximately the late Miocene. These factors are likely important drivers of intraspecific divergence within Cushion willow and possibly other alpine plants with a similar distribution. Since intraspecific divergence is often the first step toward speciation, the same factors can be important contributors to the high alpine species diversity in the HDM.

1. The GO enrichment analyses were not convincing. The discussion centers around phenotypes that can be interpreted in terms of adaptation to high elevation, but the cited categories are hand-selected from among many categories with significant enrichment, and the cited categories also have marginal statistical significance and do not involve a large number of genes. This is not a robust analysis and it inappropriately layers interpretations onto the data that are not well supported by statistics. At the very least, the authors should find a way to support their interpretations more directly, rather than referring the reader to six supplementary tables in which the cited GO categories are buried among a multitude of categories that are irrelevant to their interpretations. 2. The evidence for whole genome duplications (WGD) within the willow lineage, as indicated in Fig.  2C, is very weak and poorly supported by the presented analyses. If the two Salix species share a very recent whole genome duplication, the peak should be higher than the older Salicacae whole genome duplication. It should also be very obvious from the dotplot comparison to P. trichocarpa. Rather than representing a true WGD, this likely reveals large-scale segmental duplications, or the presence of divergent haplotypes that have been included in the main genome assembly erroneously. I think the latter is the most likely explanation, as explained above. 3. There does not appear to be any robust protection from Type I error in the Fst outlier analysis, which was apparently done without protection for multiple testing or any of the other artifacts that plague such analyses. 4. The outlier Fst analysis also suffers from the same limitation as the GO enrichment gene family analysis. The highlighted genes and categories are the ones that support the premise, not the ones with the most robust statistical significance. At the very least, some phenotypic data supporting these hypotheses would be necessary to demonstrate the robustness of the conclusions. 5. There is no justification provided for the filtering of the SNP data from GATK, and this is a critical issue since it has a dramatic impact on the LD estimates and coalescent simulations that were performed with these data. Hard filtering of SNPs, as applied here, is generally not recommended, particularly when downstream analyses depend upon low false positive and false negative rates in SNP identification. 6. The estimated genome size is 486 Mb based on flow cytometry, and this closely matches the total length of the assembly, which is good. However, this result conflicts drastically with the kmer-based size estimate of 401 Mb. As stated in the manuscript, this discrepancy could be due in part to Illumina bias, but this would mean that large portions of the assembly would be derived from PacBio only 7. More details should be given concerning the molecular dating analysis (use of priors, age of fossils, setting of constraints).
8. Concerning environmental heterogeneity between high and lower elevation sites, more information should be given about the rationale behind using the climatic parameters finally chosen, and not other abiotic parameters, like soil, or other measures of geodiversity.
Minor issues (and details) Abstract: 29-31: "The Hengduan Mountains (HDM) is a biodiversity hotspot exhibiting exceptional plant diversity, particularly for alpine plants. However, the origins of this high alpine plant richness remain largely unknown". This is not entirely true (N-S orientation of HDM enabling easier migration during climatic changes, less glaciation impact, monsoonal influence, steep gradients with shorter migration distance during climatic fluctuations in the Quaternary, etc.). See Muellner-Riehl 2019 and references therein. Apply necessary text changes in all other relevant parts of the ms as well. 35: genome "sequence"? re-phrase 38: "topographic features of the HDM are associated with…" Which features exactly? Be more precise. Testable? 40: use of "sky islands" here and elsewhere -Please provide information which definition you apply for the term "sky islands". Not all readers will be familiar with this term, and it is applied in differently stricty ways of its definition. 41: "heterogeneous environments across different elevations" Be more precise (heterogeneity measured by? Which elevations?). 44: Use of "subnival" versus "alpine" -please check for consistency of the use of these terms throughout the ms 47-50: "These findings indicate that the HDM served as an evolutionary cradle and driver of diversification of alpine species during its history, which contributed to the very high level of plant diversity characterizing this region today." Strictly speaking, this was not investigated here in the study, the latter which focused on a single species. See my comment under "major issues". 165: You may want to consult Peter Stevens´s Angiosperm Phylogeny Website (go to "Salicaceae") at http://www.mobot.org/MOBOT/Research/APweb/welcome.html to take a look at other studies which have been conducted on Salicaceae phylogenetics and molecular clock dating. 179-198: How do patterns compare to other species of Salix included in this study? 207: alpine -see my previous comment about the use of "subnival" and "alpine" 208: provide information and citation for this high temperature 220: WGD events -at this stage is is unclear to the reader where (which specimen/species) this was investigated 224-225: unclear to which level this refers to (ancestor of ?) 231: genetic structures -which ones? 234: citation 28 missing here 238: recent burst of further gene duplication -what is exactly meant by this? 240-241: "footprints of adaptive evolution in the Cushion willow genome." -anything known about the other Salix species included in this study? How do they compare in terms of ecology/biology? 243-244: "may have facilitated species-specific adaptations driven by climate change during QTP uplift and the onset of harsh subnival environmental conditions" -what kind of climate change are you referring to? (attention: uplift of the QTP and adjacent mts. was not the cause of the origination of the monsoon system, and the QTP had an average elevation of 3-4 km already since the Eocene; see reviews by Favre et al. 2015, Renner et al., 2016, Spicer, 2017Botsyun et al., 2019). Favre, A., Päckert, M., Pauls, S.U., Jähnig, S.C., Uhl, D., Michalak, I., et al. (2015). The role of the uplift of the Qinghai-Tibetan Plateau for the evolution of Tibetan biotas. Biol. Rev. Camb. Philos. Soc. 90, 236-53. doi:10.1111/brv.12107. Renner, S.S. (2016. Available data point to a 4-km-high Tibetan Plateau by 40 Ma, but 100 molecular-clock papers have linked supposed recent uplift to young node ages. J. Biogeogr. 43, 1479-1487. Spicer, R. A. (2017 [328][329]: "This overlaps with the confidence interval for the divergence time of Cushion willow, suggesting that the species origin was associated with tectonic changes in the HDM." This is an ad hoc assumption, do you have real evidence for this? How could this be tested, what are alternative explanations? See Muellner-Riehl 2019 ("clearly defined phylogeographic hypotheses may be tested rigorously, instead of applying explanations that fit the observed results ad hoc (i.e., after analysis of the data), as has been done in the past."). 336-338: See some of my previous comments about the connection between population divergence within species and species richness. As presented now, the direct link (if it does exist at all -not proven yet for this region, but some studies pointing towards this) is not convincingly presented. The study by Xing and Ree investigates species richness / is based on species-level and above, not popuation divergence, but they provide interesting thoughts concerning the HDM. See also Muellner-Riehl 2019. 353: some unique terms -meaning unclear 355: Where appropriate in the text, it would be interesting to learn more about the potentially different impact of competition (from other co-occurring plant species) that the more lowland versus subnival populations of S. brachista might be subjected to, also given the fact that the selective pressure observed seems weaker in high elevation sites. 374 and following: The situation described is the one of nowadays, formerly it was colder, populations moved down and were better connected. See Muellner-Riehl 2019. 384: elevational gradients 386: diverse heterogeneous habitats -meaning unclear, specify 388: and following: Again, please see the problematic with this statement as outlined previously. It is not merely populations being uplifted with the mountain as it rises, it is climatic fluctuations making plant populations move and being separated or connected at various times on and across mountains.

Methods:
393: Please provide voucher information on the specimen investigated (sequenced), collector and collector number, date, location, deposited at which herbarium, etc. 442: Here and elsewhere in the ms, please take care of properly setting spaces 455: artefact (though artifact also possible) 518: with the GTRGAMMA 519: Here and elsewhere in the ms, please take care of using hyphens in compound words 526: generated a maximum; which evolutionary model was used? 549-550: This sentence is not clear, please clarify and re-phrase. Certainly, you would have first run the actual NJ analysis, and then the independent boostrap analysis (the latter which is based not on the original data matrix, but the pseudo-replicates). 560 and following: Something is wrong with this sentence; phrasing and meaning of sentence unclear. Please provide the age of the fossils, based on the original primary literature. 564: Unclear -was a maximum root age set? Was a range of between 48 and 52 tested? Please specify in the text the rationale behind the setting of temporal constraints, and provide more detailed information about the priors. 602 and following: Provide information about the rationale behind considering only climatic parameters (and the ones you have chosen). What about soil or other geodiversity parameters? What would be the potentially most important determinants of occurrence for S. brachista? 609: Shannon, Simpson References: Need to be carefully checked for typos (diacritics in author names), formatting, use of upper and lower case, etc. throughout. 1f) Scheme or schematic? (both terms possible) This picture is very simplified, and does not consider mountain morphology (which is different for HDM versus Himalayas or QTP). You may consider using a real cross-section (compare Muellner-Riehl 2019) and then simplify it (smoothing the lines). Be aware that in the Himalayas/QTP, area size for alpine plants will increase at higher elevation. Figure 2: c) Please check whether the legend in the right-hand upper corner is correct. Insert "events" after "(WGD)". This is a very interesting and impressive study in which the authors conduct an ambitious population genomic analysis for a plant species native to an understudied biodiversity hotspot. The authors generated the first reference genome for Salix brachista and it important to note that they have been successful in scaffolding their reference genome into what appear to be the 19 chromosome of this species. To do so they have used state-of.the-art methodology such combining long and short reads for the assembly and HiC for the scaffolding. They also propose a very interesting set of 77 resequenced accessions sampled in localities with very different environmental and climatic conditions. I would like to felicitate the authors for having undertaken this comprehensive analysis of an interesting non-model plant species. I am convinced that these new genomic resources will be of great interest to the community. However, I have been disappointed by a series of imprecisions, inaccuracies through my reading of the analyses of this new data. But I am convinced that if the authors can address the issues listed in my review to meet the standard of the field, they could finalize a great contribution to the plant population genomic literature.
l. 114: The authors state here that given the measured heterozygosity of 0.2%, the genome must be highly repetitive. It seems that this is deduced from the joint observation of the k-mer distribution and the observed heterozygosity but it is unclear from the text what justifies making this statement. Please describe the causal relations implied here. Figure 2c: The annotation of this figure suggest that each peak is associated to a whole-genome duplication. I wonder, however, what effect the duplication of a single or two chromosomes would have on such distributions. For example, could it be that the most recent peaks in the S. brachista distribution be caused by the duplication of the Hic groups 8/19 and 9/12 only? Is the term wholegenome duplication applied loosely to refer to a duplication of a substantial part of the genome? Because the locations of each pairs of paralogues used to calculate the Ks distribution are known, I assume it should be possible to find out which chromosomes are associated with each of the three peaks.

Phylogenetics
PAML analyses: the authors write in the M&Ms that they have estimated the Ka/Ks ratio for each branch in the phylogeny (note that there is more than one ratio per branch in both models). How does the number of positively selected gene in S. brachita (1582, l. 175) compares to the distribution of equivalent numbers on all other branches. This is important to help the reader evaluate whether there is a particular signal of adaptation on the external branch leading to S. brachita or whether this very large proportion (1582/11288) could be an artefact of the statistical method employed here. l.180: significant GO enrichment: what set is enriched, PSG? The authors are not being precise here.
l.182-183: in Sup Table 8 (line 252), the GO term "regulation of response to stimulus to DNA damage stimulus" is associated with a q-value of 0.96, that is non-significant. At line 180 the author claim this to be significant. Is there a mistake somewhere? l.183-185: same here, all these terms are associated with q-values > 0.05 according to Sup. Table 9). It looks like the authors have based their assessment of significance on p-values not corrected for multiple testing.
l.185-again, same problem for the PSG class. None of these GO terms seems to be significant at the 5% level after correcting for multiple testing. Furthermore, the classes reported by the authors in this paragraph are not even the one associated with lowest q-values (see Sup Table 8,9,10). Rather, it seems that the authors have cherry-picked these classes to be consistent with their a priori hypothesis of an adaptation to UV-radiation-induced damage. My conclusion for this paragraph is that either the authors manage to make sense of the list of significant GO terms (i.e. q-value < 0.05) or the whole paragraph should be removed. My personal recommendation is to remove this paragraph.
l-189-192: this result has no relevance given my comments above.
No single sentence in this paragraph is backed by the results presented in the manuscript.
l. 240-241: This assertion should be reformulated to clearly indicate that it is based on results from PAML only. Given the very large number of genes identified as significant and the poor treatment of the GO analysis, I am not sure how much confidence can be given to this statement. l-241->end of paragraph: no results are presented that show a significant enrichment of PSG genes in species-specific or species-specific-expanded gene families. The conclusions reached here a pure speculation.
Genetic structure and demographic analysis -) Please add a table with population specific values of theta, Taj D , Fst calculated by excluding genic regions from the data.
-) authors should report the cross-validation error plot for the structure results.
-) Model-based phylogenetic inference as presented in sup Figure 9 is not an acceptable representation of intraspecific variation data. This is because recombination can occur within each locus/gene and violate the assumption of a bifurcating tree assumed by the method. Based on this, I suggest to remove this analysis from the manuscript. Also it does not provide additional information to the structure and pca results. (I feel less strongly about the NJ tree because it can be seen as a clustering method and not an explicit model-based inference method) -) tree mix does not infer a phylogenetic tree because the units here are not species but populations. Rather it infers the patterns of split and mixture in the history of a set of populations.
-) line 286. I agree that the Fst can be influenced by the distance/migration between populations in this study. But it could also be affected by the history of recent population splits (see Figure 3d). See my comment below about model-based demographic inference .
-) stairway plot analysis: this is a very interesting analysis, which seems to have been conducted very carefully (SFS estimated with ANGSD, re-estimation of the mutation rate). How is this method affected by the significant population differentiation reported in this manuscript? The authors do not specify which populations/accessions have been used to create the SFS used for the stairway analysis. I would like to see a figure with one stairway plot for each identified cluster in the structure analysis ( Figure 3a Line 307-311: this sentence sounds like the HD mountains were adapted to subnival conditions. Is this what the authors meant to write? Did they mean that the Cushion willow was well adapted to the HD mountains during glacial periods?
Origin of cushion willow associated with HDM tectonic events.
It is not acceptable to use phylogenetic reconstruction to evaluate the demographic history of these populations. As mentioned before this is because there is not a single bifurcating tree underlying this genomic data but many (due to recombination), see https://www.nature.com/articles/s41559-018-0717-x for an example. Instead the authors should estimate the age of the splits between the different genetic clusters (or main clusters) using either the software dadi (https://bitbucket.org/gutenkunstlab/dadi) or fastsimcoal2 (http://cmpg.unibe.ch/software/fastsimcoal2/). Note that the authors should try to construct the SFS using neutral sites (i.e excluding genic regions for example) Selection on Cushion willow at high and low altitudes.
In its current form the genome scan for selection to higher altitudes does not meet the standards of the field. The authors should compare the distribution of the statistic(s) used for scanning the genome to a null distribution obtained from a neutral demographic model. Such a model can estimated using one of the two tools mentioned above (dadi or fastsimcoal2). The authors should use these models to define a statistical threshold indicating how extreme the genome scan statistics should be expected to be given the demographic history of the cushion willow. Also note that tests like HapFLK, EHH (his, XP-EHH, implemented in selscan) and the CLR test implemented in the software sweepfinder are likely to be more powerful to detect selective sweeps than the statistics used by the authors so far.
Please, also indicate what is known about the breeding system of the species. Is it an outcrosser or a selfer?
Minor 149: Where are the references to these datasets? Figure 2a: The tracks "iii" and "iV" of the circus plot can not be read from a printed version. Even by zooming in track "iv" does not seem to display any information.
We would like to thank the three reviewers for their positive and constructive comments and suggestions which have been very helpful in revising and improving our manuscript. We have carefully revised our manuscript and provide below point-by-point responses (in blue) to all comments received. A major change in our revised manuscript is the inclusion of an updated Cushion willow genome to address the problem of the large discrepancy between k-mer-based estimated genome size (~400 Mb) and assembled genome size (~480 Mb). We present a PCR-free based k-mer estimate of the genome size of Cushion willow which confirms it to be ~400 Mb. We further conducted Oxford Nanopore Technology (ONT) long read sequencing, and re-assembled the Cushion willow genome based on ONT long reads and polished by PacBio and Illumina reads. As a result, we found multiple small fragments of redundancy in the previous assembly, which might be caused by alternative haplotypes. The contiguity of the updated Cushion willow genome assembly is significantly enhanced, with contig N50 of ~9.5 Mb and scafford N50 of ~17.9 Mb, and in line with k-mer-based genome size estimation. We replaced all statistics and results of the updated genome assembly in the revised genome.
We then repeated all analyses based on the updated reference Cushion willow genome and yielding results similar to those reported previously, addressing concerns and comments of reviewers in the process (e.g., stricter variant and SNP calling criteria). All revised portions of the manuscript are marked in red.

Point-by-point response to reviewer comments: Reviewer #1 (Remarks to the Author):
This manuscript describes a whole genome sequence for a dwarf willow, Salix brachista, together with genome resequencing of 77 individuals from populations in the Hengduan Mountains in China. The genome assembly is impressive and highly contiguous, and some of the analyses are interesting and reveal some insights about the evolutionary history of this species.
However, I thought there were several quite serious flaws with the interpretation of the genomic data, as detailed below. In particular, I disagree with the premise that this manuscript demonstrates "alpine environmental responses, e.g., to hypoxia and DNA damage." (lines 47-48). Also, given that the manuscript focuses on a single species with potentially unusual life history, I disagree with the assertion that the manuscript provides new, unambiguous insights about species diversification in the Hengduan Mountains. I recommend that the authors restrict their range of inference to this single species, reduce their reliance on hand-selected categories of genes to weave adaptationist stories, and tighten up their interpretation of the genomic data. The basics of an excellent manuscript are here, but it needs further work.
Response: Thanks very much for your comments. Regarding the issue of alpine environmental response, Cushion willow is a typical alpine plant in the HDM and adjacent areas which encounters a subnival environment; it is reasonable to expect, therefore, that adaptive footprints related to alpine adaptation might be present in its genome. We therefore searched for such genomic signals of adaptation based on an analysis of gene families. In the former version of the manuscript, we included some GO or KEGG terms with significant p values (< 0.05) but not significant q values (< 0.05). In the revised version, we only discuss possible alpine adaptive footprints based on terms with significant p values (< 0.05) corrected for multiple comparisons using the method of Benjamini and Hochbergsignificant. Since there were no significant enriched GO/KEGG terms in positively selected genes identified by the branch-site model using CodeML in PAML package, we removed this analysis as suggested by reviewer #3.
Regarding the origin of HDM biodiversity richness, it is true that we only conducted a case study of one species -Salix brachista -a dominant component of many subnival assemblages in the HDM. Our results indicate that diverse isolated heterogeneous habitats on sky islands, together with climatic fluctuations that have possibly caused recurrent cycles of admixture and isolation between some populations (MIM cycles), were likely important drivers of intra-specific population diversification in Cushion willow, and by extension, other alpine species with similar distributions and histories in the HDM. Because intraspecific divergence is often a first step in speciation, we propose that such factors were also likely important contributors to the generation of high alpine species diversity in the HDM.

Major Issues
Comment 1: The GO enrichment analyses were not convincing. The discussion centers around phenotypes that can be interpreted in terms of adaptation to high elevation, but the cited categories are hand-selected from among many categories with significant enrichment, and the cited categories also have marginal statistical significance and do not involve a large number of genes. This is not a robust analysis and it inappropriately layers interpretations onto the data that are not well supported by statistics. At the very least, the authors should find a way to support their interpretations more directly, rather than referring the reader to six supplementary tables in which the cited GO categories are buried among a multitude of categories that are irrelevant to their interpretations.

Response:
In the revised manuscript, we only report and discuss alpine adaptation based on enriched GO/KEGG terms with significant p values (<0.05) corrected for multiple comparisons.
We consider it is logical to do this because Cushion willow is an alpine plant which grows in subnival habitats and is likely, therefore, to have evolved adaptations to such conditions with signatures of such adaptation in its genome. Morphological variation between Cushion willow populations is minor, except for redness of branchlets, which differs between elevations. We measured the anthocyanin concentration in lowland and highland Cushion willow, as well as greenhouse maintained individuals, and observed it to be higher in high elevation individuals., where it may function to protect plants against UV radiation.
Following the suggestion of Reviewer #2, we also compared GO/KEGG enrichment results of Cushion willow with two other Salix species (S. purpurea, S. suchowensis) analyzed in this study, and did not find significant enriched GO/KEGG terms related to alpine adaptation (i.e., DNA repair, response to hypoxia) in either of these species (supplementary tables 10-13).

Comment 2:
The evidence for whole genome duplications (WGD) within the willow lineage, as indicated in Fig. 2C, is very weak and poorly supported by the presented analyses. If the two Salix species share a very recent whole genome duplication, the peak should be higher than the older Salicacae whole genome duplication. It should also be very obvious from the dotplot comparison to P. trichocarpa. Rather than representing a true WGD, this likely reveals large-scale segmental duplications, or the presence of divergent haplotypes that have been included in the main genome assembly erroneously. I think the latter is the most likely explanation, as explained above.
Response: Using the updated genome assembly, we re-identified orthologous genes and repeated the WGD analysis. We only observed an obvious peak representing the common WGD event in Salicaceae, indicating that the additional peaks observed previously in Cushion willow were caused by alternative haploid or redundancy in the previous genome assembly (please see detail explanation in response below). Consequently, we modified the text accordingly.
Comment 3: There does not appear to be any robust protection from Type I error in the Fst outlier analysis, which was apparently done without protection for multiple testing or any of the other artifacts that plague such analyses.

Response:
In the revised version, we did not use Fst and θπ outliers to identifiy genomic regions under selection. As suggested by Reviewer #3, we used instead a composite likelihood ratio statistic (CLR) method implemented in Sweepfinder2.

Comment 4:
The outlier Fst analysis also suffers from the same limitation as the GO enrichment gene family analysis. The highlighted genes and categories are the ones that support the premise, not the ones with the most robust statistical significance. At the very least, some phenotypic data supporting these hypotheses would be necessary to demonstrate the robustness of the conclusions.

Response:
We did GO and KEGG enrichment and reported enriched terms with significant p values (<0.05) corrected for multiple comparisons in the revised manuscript. We also added some phenotypic data, i.e., anthocyanin concentration data for Cushion willow individuals from two different elevation sites (2950 and 3950 m), as well as individuals maintained in a greenhouse.

Comment 5:
There is no justification provided for the filtering of the SNP data from GATK, and this is a critical issue since it has a dramatic impact on the LD estimates and coalescent simulations that were performed with these data. Hard filtering of SNPs, as applied here, is generally not recommended, particularly when downstream analyses depend upon low false positive and false negative rates in SNP identification.

Response:
We agree that filtering of the SNP data has a dramatic impact on subsequent analyses. In the manuscript, we therefore did not use hard filtering anymore; instead, we used much stricter variant and SNP calling and filtering criteria in the manuscript. These criteria were used in a recent similar genome-wide analysis in Populus (Wang, J. et al. 2018. A major locus controls local adaptation and adaptive life history variation in a perennial plant. Genome Biology 19:72.). We found a similar number of SNPs (~6.5 M) by SNP calling as before, but using the stricter filtering criteria, we obtained ~1.6 M high quality SNPs instead of the previous ~2.9 M.

Comment 6:
The estimated genome size is 486 Mb based on flow cytometry, and this closely matches the total length of the assembly, which is good. However, this result conflicts drastically with the kmer-based size estimate of 401 Mb. As stated in the manuscript, this discrepancy could be due in part to Illumina bias, but this would mean that large portions of the assembly would be derived from PacBio only sequence. Given that the authors must have aligned the Illumina reads back to the assembly, they should be able to estimate how much of the genome was not covered by Illumina. This should be reported.

Comment 7:
If it is true that 80 Mb of the genome has no Illumina coverage, then the error correction would not have been applied as in the rest of the genome. This would result in an overestimate of polymorphisms for those parts of the genome. This could result in several artifacts, including failure to remove alternative haplotypes with the very stringent filtering that was applied (98% identity). I suspect that this genome assembly still contains alternative haplotypes, and that accounts for the discrepancy in kmer and flow cytometry genome size estimates.

Response to comments 6 and 7:
We also considered the genome size of our assembly might be a serious problem. We therefore conducted PCR-free Illumina sequencing that is not affected by PCR biases to obtain a more accurate estimate of genome size. This also yielded an estimate of ~ 400 Mb, as before, so we considered the flow cytometry estimation, which had a wide range (431-539 Mb), could be erroneous. We therefore repeated the flow cytometry analysis again using another machine and obtained a genome size estimate of ~421±8 Mb based on 17 replicates (supplementary table 1), which is closer to the k-mer based estimate.
Following this, we conducted Oxford Nanopore Technology (ONT) sequencing to obtain longer reads to try to solve the problem of alternative haplotypes and redundancy in our genome assembly. A new version of the Cushion willow genome was obtained based on ONT long reads and polished by PacBio long reads and Illumina short reads. The contiguity of the updated genome was significantly enhanced, with a size of 339.6 Mb, containing 78 contigs, with contig N50 = 9.52 Mb and scafford N50 = 17.92 Mb. We found there were large amounts of high frequency (> 1×10 5 ) k-mers in the PCR-free reads (1705 kinds of 17-k-mers sum to ~1.626 Gb), which might originate from highly repetitive regions such as centromeres and telomeres. These high frequency k-mers amount to ~47.8 Mb, meaning that the genome size is reduced to ~350 Mb if these high frequency k-mers are excluded. It is very difficult to successfully assemble the centromeres and telomers based on current sequence data, because the sequence read length is still not long enough, and since in subsequent analysis repeat regions are masked. Consequently, we did not attempt to assemble these highly repetitive sequences.
The reasons for the large difference in genome size between the updated genome assemby and the previous one might be: (1), when the there are bubbles in the assembly graph caused by heterozygosity, HGAP/CANU tend to assembly all allelic haplotigs, while SMATRTdenovo/CANU tend to combine allelic hplotigs into one; (2), long reads generate by ONT can successfully sequencing high teterozygous regions more effectively and therefore combine allelic haplotigs more effective.
When we mapped all types of sequence reads (ONT, PacBio and Illumina) back to the updated genome assembly, 98.49% of Illumina short reads were mapped successfully back to the genome assembly, and only 1.18% of the genome assembly was not covered by Illumina short read (supplementary Table 6), implying that the current assembly covered most unique genomic regions and high accuracy of the genome assembly. We report this result in the revised manuscript. Fig. S7 is quite interesting, and seems to show some large-scale duplications in cushion willow, and an apparent lack of the large discrepancy found between Chr1 and Chr16 in S. suchowensis. However, it is very hard to tell due to the high background (the vertical lines). Please produce a clearer version of this figure and discuss the genome comparisons more thoroughly. It would be quite interesting if cushion willow also has the Chr1 fusion that is present in Populus but absent in S. suchowensis and S. purpurea.

Comment 8: The dot plot in
Response: Sorry for the unclear figure. We present a clear version of dotplot (Fg. 2c, Supplementary Fig. 8) based on the updated genome assembly. We compared not only Cushion willow with Populus trichocarpa, but also Populus deltodies and Salix purpurea. We show that genome synteny within Salix and Populus species are both high, but some large rearrangements were evident in Salix chromosomes 4, 15 and 19. The fission and fushion events between choromosomes 1 and 16 were also observed (Fig 3d, Supplementary Fig. 8) when comparing Salix species with Populus species. We therefore added these results to the revised manuscript.

Minor Issues
Comment 9: Lines 112-114. The rationale for concluding that the genome is highly repetitive based on the heterozygosity rate of 0.2% is not clear. Are the authors concluding that this level of polymorphism is too low to account for the second peak in the kmer distribution? This has not been demonstrated, and in my experience this assumption would be untrue. In any case, the rationale for this conclusion needs to be more clearly stated.

Response:
We agree that the conclusion is not appropriate and therefore deleted it from the revised manuscript. We estimated heterozygosity rate based on the updated genome assembly, it is 0.71%, i.e. higher than the previously version.
Comment 10: Lines 137-138. The manuscript states that 97.85% of genes were annotated based on sequence similarity to annotated proteins in public databases. This is likely to be misleading, since many of the "annotations" in public databases are probably "gene of unknown function" or something similar. I think the authors mean that 97.85% of the predicted proteins matched a predicted protein in public databases. "Gene of Unknown Function" is not an annotation. The functional annotations should be quantified more precisely, or the statement should be clarified.
Response: Thanks for this suggestion, we therefore rephrased this as: "Most of the predicted proteins (97.4%) matched a predicted protein in public databases (Supplementary Table 9)". Comment 11: Lines 160-165. The analysis reported here is not a robust demonstration of monophyly of Salix and Populus, which is not really in doubt in any case. I recommend deleting this paragraph. However, I agree that the phylogeny is a useful framework for presenting the gene family results, so Fig. 2b is helpful.
Response: It's true that the monophyly of Salicaceae as well as Salix and Populus is solid, as proved by previous studies including ours, so we deleted this content from the revised manuscript.

Comments to the authors
The authors report a high-quality, chromosome-scale genome of Cushion willow (Salix brachista), which they use as a reference when resequencing 77 individuals from 14 populations across the species' distribution range (incl. sites in the Hengduan Mountains region and neighbouring regions), to characterize the species´ genetic diversity, population structure, demographic history, and footprints of adaptive evolution to the subnival environment. The authors interpret their results in the context of explaining high plant species richness within the HDM.
I would like to express my great appreciation to the authors for this extensive study, which aims at bridging the gap between genomics, population genetics/genomics, and species richness patterns, at the backbone of a geologically active mountainous region. I found the manuscript very interesting to read, but also came across several issues that will deserve some careful attention. A very recent paper (Muellner-Riehl, 2019), dealing with the geographic region which is also the focus of this manuscript, tackles many of these issues and provides suggestions for future research avenues, some of which ideally could be taken up by this study.
In the following, where appropriate, I will refer to the line numbers as given in the merged pdf file. I hope that my comments will be helpful for the authors.
Response: Thanks very much for your positive comments to our study and highlighting the recent paper by Muellner-Riehl .

Major issues
Comment 1: The authors interpret their population-level, infra-specific results in the context of explaining high plant species richness within the HDM. The connection between the two (which would be speciation) is currently not convincingly presented in the manuscript. There have been some recent attempts to view speciation and the origination of high levels of species richness in the same context. For example, the recently presented Mixing-Isolation-Mixing (MIM) model of speciation (He et al. 2018;compare comments by Abbott, 2018, andRieseberg, 2018), supports that speciation may involve repeated cycles of admixture, followed by geographic isolation until complete isolation is finally achieved. By permitting intermittent gene flow during speciation, MIM cycles can potentially generate species at higher rates than just by geographic isolation, which in turn means that speciation and the origination of high biodiversity levels in a region may be viewed in the same context framework (Heet al., 2018). Response: Thanks for these valuable suggestions and comments. First, we want to emphasize that we have updated the Cushion willow genome assembly (see above) and performed all analyses again, yielding results similar to those reported previously.
Our results indicate that S. brachista split from one of its close relatives around 13 Ma, and subsequently underwent intraspecific divergence more recently (around 0.18-0.42 Ma), i.e. after the intense uplift of the HDM. In the revised ms, we propose that our results are consistent with the hypothesis that diverse isolated heterogeneous habitats on sky islands, together with climatic fluctuations that have possibly caused recurrent cycles of admixture and isolation between some populations (MIM cycles), were important drivers of intra-specific population diversification in Cushion willow, and by extension, other alpine species with similar distributions and histories in the HDM. We further speculate from our results that geographic isolation and environmental variation associated with elevation resulting from the uplift of THR, together with the effects of Quaternary climate fluctuation, have contributed to the high species richness observed in the HDM by driving intraspecific population divergence, a common precursor of speciation.

Comment 2:
The authors argue for uplift-driven diversification (as many other authors have done before) in an ad hoc manner, but they don't test this. What could be alternative explanations for the patterns observed? How could different scenarios be compared and tested? What is the most likely explanation? (see Muellner-Riehl 2019 for some suggestions). Connected to the issue raised above, the relation of population-level, infra-specific differentiation and (species-level) diversification is not apparent.

Response:
In the revised manuscript, we estimated the divergence time of Cushion willow from one of its closest relatives. Crown age estimation of Cushion willow was deleted (Reviewer #3 pointed out it is not appropriate to estimate divergence time of a species using individuals of it in the context of phlygenetic tree construction). We also estimated the split times of Cushion willow populations as suggested by Reviewer #3 using fastisimcoal2. Our results indicate that S. brachista split from one of its closest relatives around 13 Ma, in accordance with the intense uplift of the QTP and Himalayas, whereas populations of Cushion willow diverged from each other much more recently (around 0.18-0.42 Ma), i.e. after the intense uplift of HDM. Therefore, the origin of S. brachista might be associated with tectonic events, while its population divergence was likely driven by habitat isolation (across sky islands) and possibly MIM cycles as mentioned previously.
Comment 6: Of a somewhat similar problematic is the discussion of population bottlenecks with regard to glaciations. Information needs to be provided about the regional extent of these glaciations (into the species distribution range), how these glaciations would have possibly impacted vegetation cover in the HDM, either directly by glacial cover, or a drop in temperature.
The age for the LGM given is too young, and it is known that the maximum glacial extent in the THR was not at the time of the LGM, but pre-dated the LGM (see Muellner-Riehl 2019 and references therein for details), thus rendering the time of the global LGM as not specifically relevant for the THR.

Response:
We have re-analyzed the demographic history of Cushion willow based on the genetic groups identified by genetic structure analysis as suggested by Reviewer #3. We found multiple bottlenecks to occur for most of these groups, but their timing and number differed. Such variation in demographic history between these groups is not unexpected given that glaciation of the THR region was highly variable, with strong microclimatic variation present within individual mountain ranges. Each bottleneck event was usually followed by a population expansion, reflecting major climatic oscillations during the Pleistocene. Since bottleneck events are population or mountain specific and do not accord with the five distinct glaciation periods in the THR, we deleted these glaciation periods from figures and present a δ 18 O curve (which is an indicator of global ice volume and temperature) instead (Fig. 3f, Supplementary Fig. 10).
Comment 7: More details should be given concerning the molecular dating analysis (use of priors, age of fossils, setting of constraints).

Response:
We have added these details in Methods. In the new version of the ms, we estimate the divergence time of Cushion willow from another dwarf willow in THR (Salix souliei, also placed in section Lindleyanae). Since the earliest reliable fossil of Salicaceae is in the early Middle Eocene, about 48 Ma, and another study estimated the divergence time of Populus and Salix lineage as ~52 Ma, we calibrated the root age as 48-52 Ma.

Minor issues (and details)
make the use of "subnival" and "alpine" appropriate depend on specific context. Comment 14: 47-50: "These findings indicate that the HDM served as an evolutionary cradle and driver of diversification of alpine species during its history, which contributed to the very high level of plant diversity characterizing this region today." Strictly speaking, this was not investigated here in the study, the latter which focused on a single species. See my comment under "major issues".

Response:
We re-wrote this part emphasizing intraspecific divergence, often a first step in speciation: "These factors, Response: We have capitalized the first letter for these words in the revised ms.

References:
Comment 61: Need to be carefully checked for typos (diacritics in author names), formatting, use of upper and lower case, etc. throughout.

Response:
We have checked the references section for these issues carefully and made corresponding corrections. Response: Thanks, we have made corresponding corrections as suggested. We also have changed to a map which shows mountain morphology. Figure 2: c) Please check whether the legend in the right-hand upper corner is correct. Insert "events" after "(WGD)".

Response:
Since we only found the common WGD event in Salicaceae, we do not show the Ks distribution in Fig. 2, and only present it in Supplement Fig. 7 without this legend. Response: Sorry for the typo, it is "time". Since bottleneck events are population or mountain specific and do not accord with the five distinct glaciation periods in THR, we deleted these glaciation periods from figures and present a δ 18 O curve (which is an indicator of global ice volume and temperature) instead (Fig. 3f, Supplementary Fig. 10). Related content in the manuscript was deleted.

Figures:
Comment 65: Improve readability of fonts used (e.g., Fig. 1e, river and mountain names) Response: We changed to a map with clear mountain terrain. The mountains are hard to distinguish, so we deleted their names from the map.

Final:
Comment 66: -Carefully read, check and correct for clarity and logic of English language and use of scientific terms throughout the manuscript.

Response:
We have checked throughout the manuscript for the issues listed. Response: Thank for this suggestion, we have read and compared this study with ours, it helped us in revising the manuscript.

Reviewer #3 (Remarks to the Author):
This is a very interesting and impressive study in which the authors conduct an ambitious population genomic analysis for a plant species native to an understudied biodiversity hotspot. The authors generated the first reference genome for Salix brachista and it important to note that they have been successful in scaffolding their reference genome into what appear to be the 19 chromosome of this species. To do so they have used state-of.the-art methodology such combining long and short reads for the assembly and HiC for the scaffolding. They also propose a very interesting set of 77 re-sequenced accessions sampled in localities with very different environmental and climatic conditions. I would like to felicitate the authors for having undertaken this comprehensive analysis of an interesting non-model plant species. I am convinced that these new genomic resources will be of great interest to the community. However, I have been disappointed by a series of imprecisions, inaccuracies through my reading of the analyses of this new data. But I am convinced that if the authors can address the issues listed in my review to meet the standard of the field, they could finalize a great contribution to the plant population genomic literature.
Response: Thanks very much for your comments on our study. We are sorry for the imprecisions and inaccuracies you encountered. We have addressed these issues when revising the ms.

Genome assembly.
Comment 1: 114: The authors state here that given the measured heterozygosity of 0.2%, the genome must be highly repetitive. It seems that this is deduced from the joint observation of the k-mer distribution and the observed heterozygosity but it is unclear from the text what justifies making this statement. Please describe the causal relations implied here.

Response:
We agree that this statement is not convincing, so we delete this. Moreover, we updated the Cushion willow genome assembly based on Nanopore Oxford Technology (ONT) long reads, contiguity of the genome assembly was significantly enhanced, and we also found that the previous genome assembly included many small fragments of redundancy (please also see our response to comment 2). Heterozygosity rate based on the updated genome assembly is 0.71%, higher than previous estimated. Figure 2c: The annotation of this figure suggest that each peak is associated to a whole-genome duplication. I wonder, however, what effect the duplication of a single or two chromosomes would have on such distributions. For example, could it be that the most recent peaks in the S. brachista distribution be caused by the duplication of the Hic groups 8/19 and 9/12 only? Is the term whole-genome duplication applied loosely to refer to a duplication of a substantial part of the genome? Because the locations of each pairs of paralogues used to calculate the Ks distribution are known, I assume it should be possible to find out which chromosomes are associated with each of the three peaks.

Comment 2:
Response: The S. brachista-specific peak and the peak shared by S. brachista and S. purpurea were very weak (especially in S. pupurea) (reviewer #1 also addressed this issue). In the WGD analysis based on the updated S. brachista genome assembly, we only observed the common WGD event in Salicaceae, indicating that the additional peaks in S. brachista resulted from redundant sequence (Supplementary Fig. 7).

Phylogenetics
Comment 3: PAML analyses: the authors write in the M&Ms that they have estimated the Ka/Ks ratio for each branch in the phylogeny (note that there is more than one ratio per branch in both models). How does the number of positively selected gene in S. brachita (1582, l. 175) compares to the distribution of equivalent numbers on all other branches. This is important to help the reader evaluate whether there is a particular signal of adaptation on the external branch leading to S. brachita or whether this very large proportion (1582/11288) could be an artefact of the statistical method employed here.

Response:
We have redone the PSG analysis, we found many PSGs in S. brachista and S. suchowensis, and much less PSGs in other species and branches (please see the below table for detail). But like previous results, we did find significant enriched GO/KEGG terms in S. brachista, so as you suggested in comment 7, we remove this analysis and related contents from the revised manuscript.

Salix brachista 1163
Ancestral Salix 73 Comment 4: l.180: significant GO enrichment: what set is enriched, PSG? The authors are not being precise here.

Response:
We re-wrote this paragraph using new gene family enrichment (including expanded, fast evolve and species-specific gene families, and PSG analysis was removed as mentioned in response to comment 3) results based on the updated S. brachista genome assembly.
See the revised manuscript for details (line 176-200).
Comment 5: l.182-183: in Sup Table 8 (line 252), the GO term "regulation of response to stimulus to DNA damage stimulus" is associated with a q-value of 0.96, that is non-significant. At line 180 the author claim this to be significant. Is there a mistake somewhere? (Fig. 3f, Supplementary Fig. 10). Moreover, since bottleneck events are population or mountain specific and do not accord with the five distinct glaciation periods, we deleted these glaciation periods from figures and text.

Comment 19: Line 295: there is no figure 3f
Response: Sorry for this mistake, it should have been Fig 3e. However, we added a panel to Figure 3, so there is Figure 3f in the revised manuscript. Response: In our new stairway plot analysis based on genetic groups, some have confidence intervals for Ne estimates shrinking to 0 for some period while others do not. We think this is a normal output from the analysis.
Comment 21: Line 307-311: this sentence sounds like the HD mountains were adapted to subnival conditions. Is this what the authors meant to write? Did they mean that the Cushion willow was well adapted to the HD mountains during glacial periods?
Response: We mean Cushion willow was well adapted to subnival conditions in HDM, so we rephrased it as: "... it was well adapted to subnival conditions...".

Origin of cushion willow associated with HDM tectonic events.
Comment 22: It is not acceptable to use phylogenetic reconstruction to evaluate the demographic history of these populations. As mentioned before this is because there is not a single bifurcating tree underlying this genomic data but many (due to recombination), see https://www.nature.com/articles/s41559-018-0717-x for an example. Instead the authors should estimate the age of the splits between the different genetic clusters (or main clusters) using either the software dadi (https://bitbucket.org/gutenkunstlab/dadi) or fastsimcoal2 (http://cmpg.unibe.ch/software/fastsimcoal2/). Note that the authors should try to construct the SFS using neutral sites (i.e excluding genic regions for example) Response: Thank you. In the revised manuscript, we have estimated split time and gene flow for six populations pairs (because this analysis is time consuming, we chose six representative pairs) by fastisimcoal2 based on SFS constructed using neutral sites.

Selection on Cushion willow at high and low altitudes.
Comment 23: In its current form the genome scan for selection to higher altitudes does not meet the standards of the field. The authors should compare the distribution of the statistic(s) used for scanning the genome to a null distribution obtained from a neutral demographic model. Such a model can estimated using one of the two tools mentioned above (dadi or fastsimcoal2). The authors should use these models to define a statistical threshold indicating how extreme the genome scan statistics should be expected to be given the demographic history of the cushion willow. Also note that tests like HapFLK, EHH (his, XP-EHH, implemented in selscan) and the CLR test implemented in the software sweepfinder are likely to be more powerful to detect selective sweeps than the statistics used by the authors so far.

Response:
In the revised manuscript, we have not used Fst and θπ outliers to scan for selective sweep regions. Instead, we scan for selective sweep region using CLR implemented in sweepfinder2 based on site ancestral state as suggested. We identify selected region using the top5% CLR values, this kind of criteria, that using genome-wide distribution as a null for identifying selective sweep was also used in some recent studies (see references below). We did not found significantly enriched GO/KEGG terms for the genes residue in the sweep region, we therefore discuss some genes with the very high CLR values. Comment 24: Please, also indicate what is known about the breeding system of the species. Is it an outcrosser or a selfer?

Response:
The Saliaceae family is dioecious, so S. brachista is an outcrosser. We have added this information in the Introduction (line 69).

Minor
Comment 25: 149: Where are the references to these datasets?

Response:
We have added the URL and reference for the genome datasets used in our study (Supplementary Table 23). Figure 2a: The tracks "iii" and "iV" of the circus plot can not be read from a printed version. Even by zooming in track "iv" does not seem to display any information.

Comment 26:
Response: Sorry, we have added these tracks (i-vi) in the circus plot (Fig. 2a) in the revised manuscript.

Reviewers' comments:
Reviewer #1 (Remarks to the Author): Alpine plant divergence in the Hengduan Mountain biodiversity hotspot: insights from a genome-wide evolutionary analysis of Cushion willow This manuscript is greatly improved over the original submission. The manuscript describes the genome of the cushion willow and makes inferences about adaptation to high elevation based on gene family expansion and population genomics analyses. They also infer that the sky island distribution of the populations has contributed to diversification of this and other species in the Hengduang Mountain region of China. The manuscript is well written and the results are quite interesting. The authors have also made many satisfactory responses to the original reviews. The genome assembly is muchimproved and the authors have satisfactorily addressed concerns about the previous analyses of genome size and heterozygosity. The selective sweep analysis is also greatly improved, and yielded some very interesting results. I also think that the authors have sufficiently reduced their range of inference and the statements they make about diversification in the THR and HDM are now much more consistent with the data that they present in the manuscript. Nevertheless, I have detected some persistent and additional flaws in the analyses and interpretations, and these greatly dimmed my enthusiasm for this work.
1. The authors do not specify the comparative database used for the GO and KEGG enrichment analyses. What is enrichment being compared to? These methods need to be expanded, because this is absolutely critical to proper interpretation of these tests. As presented, I cannot tell what is being compared to what. To be more specific, the authors initially report on gene orthology analyses involving 10 Malpighiales species, which should give a very robust picture of which genes are exclusive or enriched in cushion willow. They then report the GO enrichment results. What is the input gene list for this analysis, and what is the background? For the questions being posed, I think the optimal analysis would focus on genes that are unique or expanded (or rapidly evolved) in cushion willow compared to the other Malpighiales, and the background set should be the full complement of genes in cushion willow. Other gene lists or reference databases could also be appropriate, but the question being posed needs to be clearly stated, and the reference database should be clearly identified. Based on the legends of the supplementary tables I guess that the input list is correct, though this needs to be explicitly stated in the manuscript. I can't find any information about the background set(s). 2. Even if the GO enrichment analyses are conducted correctly, I still do not find them to be convincing evidence of adaptation to high altitude conditions. The authors are highlighting a small subset of the categories that are significantly enriched based on their a priori biological expectations that this enrichment will reflect high altitude adaptation. There are 54 and 26 over-represented GO categories reported in Supplemental Tables 11 and 12, respectively, and only 4 of these are picked out to show evidence of adaptation to high elevations, and reasonable alternative interpretations could be provided for each of these four categories. The analysis is slightly improved over the original version by eliminating categories that are not significant after correction for multiple testing, and the comparison to the other Salix species is certainly an improvement. However, this approach would be much stronger if there were a more direct comparison of the differences in gene family size between the closely-related species and a focus on expansion of genes for which there is a reasonable expectation that the selective force could be high elevation. The GO analysis lacks specificity in this regard, and therefore is insufficient for a journal of this prominence. 3. The authors have missed the opportunity to present some potentially interesting comparative analyses. For example, they report the repeat composition of the genome, but they do not put it into context. The genome consists of ~9% Gypsy and ~9% Copia elements. I believe this number is quite low compared to other sequenced Salix and Populus, but the authors do not make these comparisons. Are the repeats different or shared with other sequenced Salix? Since there are at least 3 public Salix genomes now (S. purpurea, S. suchowensis, S. viminalis) and at least 3 public Populus genomes (P. trichocarpa, P. euphratica, and P. tremula), there is really no reason not to do this, and it would add substantial value to the manuscript. 4. The selective sweep analysis is greatly improved over the original version, and it yielded some interesting results that lend more credence to the claim of high altitude adaptation than the GO/KEGG analyses. However, I think this analysis is flawed because it suffers from a severe Wahlund effect based on combining populations that are well differentiated into false groups based on their elevation. The neutral population structure analyses clearly show that the low elevation populations LJ and MN are strongly differentiated, yet the selective sweep analysis is performed on merged populations. This is likely to flatten the site frequency spectrum, which will affect the null distribution and have unpredictable impacts on individual loci. If this effect is differentially affecting the merged lowland versus highland populations, inferences about differential selection are likely to be flawed due to possible artifacts driven by different power in the two analyses. It would be preferable to run these analyses on undifferentiated populations with normalized sizes (i.e., the same number of individuals for the populations that are being compared). 5. The authors do not specify the criteria they used to define "rapid evolution" of gene families (lines 172-173) 6. Do the authors have permission to publish the Populus deltoides comparative analysis? As far as I know, this is still protected by the JGI data embargo. In fact, the authors do not specify where they acquired any of the comparative data that they used in this thesis. At the very least they should specify the source and version of the genomes used for the analyses. 7. The use of color for the population dots in Fig. 1e is potentially confusing because some of the same colors are also used in the FastStructure plots, where they have a different meaning. I suggest just using one color for the population dots, which is not reused for the FastStructure charts. The first two letters of the population names are sufficient to represent the groupings of the collections 8. The Supplementary Materials are not well organized within the provided files. It is not obvious from the file names where to find particular figures and tables, and figures and tables are mislabeled and mixed within Excel files. For example, the file 198041_1_related_ms_3888373_ptdwfy.xlsx contains a sheet labeled Figure 3d that is actually a table, and Supplementary Table 1 is actually an excel sheet with a bunch of figures. And Supplementary Figure 13 is actually a table. The file 198041_1_data_set_3888370_ptdw3c.xlsx is much better organized, but the supplemental figures are not labeled, so it is difficult to know if they are in the correct order. This may have been due to upload problems during review, but in any case I recommend that these be improved for the published manuscript. 9. Lines 260-261: The authors discuss long-distance dispersal by seed as an explanation for shared genome composition between distantly-related populations. I expect that long-distance pollen dispersal is more likely, but of course this depends on the pollination vector (wind vs insect, and which insects). 10. As pointed out by a previous reviewer, the assumption of a 3 year generation time is probably not justified. First, the citation provided is apparently not in the peer-reviewed literature, and it was provided by a willow breeder for the purpose of planning a breeding program. There is a big difference between flowering time in a nursery bed with abundant resources and flowering time in a harsh wild environment. Second, flowering time is not the same as generation time in the wild. In reality, seed establishment opportunities may only occur a handful of times over the lifetime of a successful plant, and for a clonal plant like cushion willow, this could mean an effective generation time of decades. A better estimate is needed, or at the very least a range should be tested. Nevertheless, since the authors are no longer trying to relate the timing of bottlenecks to geological events, this criticism is not fatal for the manuscript, since the relative comparisons among populations should mostly be valid. 11. Although it is generally well-written, some of the text, particularly the newer sections, have not undergone sufficient language editing. I have marked up the pdf of the manuscript to facilitate these revisions.
Reviewer #2 (Remarks to the Author): My review to R1 I am very pleased with the answers to my questions, comments and suggestions regarding the original version of the manuscript.
I am attaching an annotated pdf version of the main manuscript text, with some minor, mostly editorial corrections and comments.
With regard to the response to my "comment 45" in the first revision, I would like to suggest to the authors to also include one or two short sentences in their main ms text, providing some of the information they have given in their reply ("Thanks for this valuable suggestion. S. brachista only grows in open habitats: habitats in subnival scree slope are completely open (Fig. 1d), while in low elevation localities around 3000 m, S. brachista is only found on open river banks and never occurs in forest near the river (Fig. 1c). Therefore, it seems that competition of lowland and subnival populations of Cushion willow with co-occuring plant species are unlikely to be different.").

Alexandra Muellner-Riehl
Reviewer #3 (Remarks to the Author): I am pleased by how the authors addressed my comments. There is only one aspect that I would recommend to improve: In their selective sweep analysis, the authors consider as significant (i.e. putative sweeps) all regions with CLR values falling within the top 5% of the distribution.
The standard in the field, however, is too calculate the statistical threshold for the CLR statistic using neutral simulations obtained from the best demographic model (as obtained by fastsimcoal2 in this study for example). This is actually relatively easy to do because fastsimcoal2 creates a file called "project_maxL.par" that can be used as input file for fastsimcoal2 to SIMULATE DNA data under your best demographic mdoel for the high and low populations. A bit of scripting may be required to transform the simulated data to the correct input for Sweepfinder2. this procedure will allow the authors to identify "how high" the CLR value can be interesting. The authors have also made many satisfactory responses to the original reviews. The genome assembly is much-improved and the authors have satisfactorily addressed concerns about the previous analyses of genome size and heterozygosity. The selective sweep analysis is also greatly improved, and yielded some very interesting results. I also think that the authors have sufficiently reduced their range of inference and the statements they make about diversification in the THR and HDM are now much more consistent with the data that they present in the manuscript. Nevertheless, I have detected some persistent and additional flaws in the analyses and interpretations, and these greatly dimmed my enthusiasm for this work.
Response: Thanks for positive comments for our study. We accepted all the revisions in the attached annotated PDF file, which are marked as red in the revised manuscript. And we made revisions according to concerns addressed, please see responses below.
1. The authors do not specify the comparative database used for the GO and KEGG enrichment analyses. What is enrichment being compared to? These methods need to be expanded, because this is absolutely critical to proper interpretation of these tests. As presented, I cannot tell what is being compared to what. To be more specific, the authors initially report on gene orthology analyses involving 10 Malpighiales species, which should give a very robust picture of which genes are exclusive or enriched in cushion willow. They then report the GO enrichment results. What is the input gene list for this analysis, and what is the background? For the questions being posed, I think the optimal analysis would focus on genes that are unique or expanded (or rapidly evolved) in cushion willow compared to the other Malpighiales, and the background set should be the full complement of genes in cushion willow. Other gene lists or reference databases could also be appropriate, but the question being posed needs to be clearly stated, and the reference database should be clearly identified. Based on the legends of the supplementary tables I guess that the input list is correct, though this needs to be explicitly stated in the manuscript. I can't find any information about the background set(s).

Response:
We agree with Rev 1 that we did not describe the GO and KEGG enrichment analysis in sufficient detail in the previous version. We tested GO terms and KEGG pathway enrichments using expanded, rapidly-evolved or species-specific sets of genes for the species examined (Salix brachista, Salix purpurea and Salix suchowensis) as input gene lists, and using all the genes of the genome with GO or KEGG annotations of the examined species (Salix brachista, Salix purpurea and Salix suchowensis) as background. We explain our approach in the revised manuscript (Lines 232-235 and Lines 586-588).
2. Even if the GO enrichment analyses are conducted correctly, I still do not find them to be convincing evidence of adaptation to high altitude conditions. The authors are highlighting a small subset of the categories that are significantly enriched based on their a priori biological expectations that this enrichment will reflect high altitude adaptation. There are 54 and 26 overrepresented GO categories reported in Supplemental Tables 11 and 12, respectively, and only 4 of these are picked out to show evidence of adaptation to high elevations, and reasonable alternative interpretations could be provided for each of these four categories. The analysis is slightly improved over the original version by eliminating categories that are not significant after correction for multiple testing, and the comparison to the other Salix species is certainly an improvement. However, this approach would be much stronger if there were a more direct comparison of the differences in gene family size between the closely-related species and a focus on expansion of genes for which there is a reasonable expectation that the selective force could be high elevation. The GO analysis lacks specificity in this regard, and therefore is insufficient for a journal of this prominence.

Response:
Our expectation is that since Salix brachista is an alpine plant, we should find some hints for adaptation to alpine conditions, such as response to UV radiation (that can cause DNA damage), and low oxygen level. Some of the enriched terms may seem cherry-picked, and we agree there may be reasonable alternative interpretations for enrichment of these categories. So, we now only discuss the top two enriched KEGG terms in expanded gene families (base excision repair, flavonoid biosynthesis), and the GO terms related to flavonoid/anthocyanin biosynthesis.
We further compare the gene numbers of the gene families involved in these terms between the three Salix species examined. We found that Salix brachista has more genes for most of the gene families than the other two non-alpine Salix species (supplementary table 15). These modifications to the manuscript are presented in "Genomic footprints of alpine adaptationlines" start from Line 232 in the revised manuscript.
We think these results provide indications of genomic footprints of alpine adaptation in Salix brachista. Of course further in-depth studies are required to prove this, but our results could provide indications and candidate genes for such future research.
3. The authors have missed the opportunity to present some potentially interesting comparative analyses. For example, they report the repeat composition of the genome, but they do not put it into context. The genome consists of ~9% Gypsy and ~9% Copia elements. I believe this number is quite low compared to other sequenced Salix and Populus, but the authors do not make these comparisons. Are the repeats different or shared with other sequenced Salix? Since there are at least 3 public Salix genomes now (S. purpurea, S. suchowensis, S. viminalis) and at least 3 public Populus genomes (P. trichocarpa, P. euphratica, and P. tremula), there is really no reason not to do this, and it would add substantial value to the manuscript.
Response: Thanks for this suggestion. We did the LTR/Copia and LTR/Gypsy analysis for six Saliaceae species (we did not use S. viminalis because the genome assembly is too fragmented and still not publically available). Our results show that in comparison to genomes of closely related other Salicaceae species, Copia elements appear to have expanded in Cushion willow. In contrast, Gypsy elements have contracted in Salix species examined, although Cushion willow has the highest proportion among the Salix species. A phylogenetic analysis of Copia and Gypsy families across the six Salicaceae species revealed no evidence for any species-specific family members, demonstrating that Copia/Gypsy divergence predates the divergence of these species.
4. The selective sweep analysis is greatly improved over the original version, and it yielded some interesting results that lend more credence to the claim of high altitude adaptation than the GO/KEGG analyses. However, I think this analysis is flawed because it suffers from a severe Wahlund effect based on combining populations that are well differentiated into false groups based on their elevation. The neutral population structure analyses clearly show that the low elevation populations LJ and MN are strongly differentiated, yet the selective sweep analysis is performed on merged populations. This is likely to flatten the site frequency spectrum, which will affect the null distribution and have unpredictable impacts on individual loci. If this effect is differentially affecting the merged lowland versus highland populations, inferences about differential selection are likely to be flawed due to possible artifacts driven by different power in the two analyses. It would be preferable to run these analyses on undifferentiated populations with normalized sizes (i.e., the same number of individuals for the populations that are being compared).

Response:
We agree with these comments. So we use undifferentiated high (LJ_XF and LJ_LS) and low (LJ_BS and LJ_GH) elevation LJ populations, which are located on the same Yulong mountain in Lijiang county. Structure analysis and PCA analysis indicate these high and low populations are undifferentiated and therefore appropriate for selective sweep analysis. Nine individuals from each of the high and low elevation groups were included in the analysis. Since we use the default parameter of CAFE, we chose not to report it in detail in the Method section, and we just state that we used the default parameters.
6. Do the authors have permission to publish the Populus deltoides comparative analysis? As far as I know, this is still protected by the JGI data embargo. In fact, the authors do not specify where they acquired any of the comparative data that they used in this thesis. At the very least they should specify the source and version of the genomes used for the analyses.

Response:
We have contacted (by email) the principal investigator of the genome project of Populus deltoides, professor Gerald A. Tuskan. He told us that there would be no restrictions in the use of this data in a comparative gene family analysis focused on Salix and cite JGI and Phytozome as our source of data. If necessary, we can provide these emails.
7. The use of color for the population dots in Fig. 1e is potentially confusing because some of the same colors are also used in the FastStructure plots, where they have a different meaning. I suggest just using one color for the population dots, which is not reused for the FastStructure charts. The first two letters of the population names are sufficient to represent the groupings of the collections.   Figure 13 is actually a table. The file 198041_1_data_set_3888370_ptdw3c.xlsx is much better organized, but the supplemental figures are not labeled, so it is difficult to know if they are in the correct order. This may have been due to upload problems during review, but in any case I recommend that these be improved for the published manuscript.

Response:
We are sorry for the organization of the Supplementary Material. The name of these files seem to be automatically renamed by the submission system. We will carefully check this problem when we upload the revised manuscript files.
9. Lines 260-261: The authors discuss long-distance dispersal by seed as an explanation for shared genome composition between distantly-related populations. I expect that long-distance pollen dispersal is more likely, but of course this depends on the pollination vector (wind vs insect, and which insects). We think the use of 15 years as the maximum generation time is reasonable, because by using this generation time, time of demographic events of the DC population were estimated to start from around 16.4 Ma ( Supplementary Fig 13c), which is close to the upper bound of the estimated divergence time of Salix brachista with one of its close relative Salix souliei, i.e., 11.73 Ma (95% 95% confidence intervals: 6.83-18.34 Ma) ( Supplementary Fig 9).

Response
Our main result is that demographical history was affected by climate fluctuation, so although we used different generations times, our conclusions still stands.
11. Although it is generally well-written, some of the text, particularly the newer sections, have not undergone sufficient language editing. I have marked up the pdf of the manuscript to facilitate these revisions.
Response: Sorry for this, we have checked thoroughly for language improvement in the revised manuscript.
Reviewer #2 (Remarks to the Author): My review to R1 I am very pleased with the answers to my questions, comments and suggestions regarding the original version of the manuscript.
I am attaching an annotated pdf version of the main manuscript text, with some minor, mostly editorial corrections and comments.
With regard to the response to my "comment 45" in the first revision, I would like to suggest to the authors to also include one or two short sentences in their main ms text, providing some of the information they have given in their reply ("Thanks for this valuable suggestion. S. brachista only grows in open habitats: habitats in subnival scree slope are completely open (Fig. 1d), while in low elevation localities around 3000 m, S. brachista is only found on open river banks and never occurs in forest near the river (Fig. 1c). Therefore, it seems that competition of lowland and subnival populations of Cushion willow with co-occuring plant species are unlikely to be different.").

Alexandra Muellner-Riehl
Response: Thanks very much for positive comments for our study. We added sentences to describe the competition of high and low elevation Cushion willow plant in Lines 367-368. We also accepted all the revisions in the attached annotated PDF file, which are marked as red in the revised manuscript.
Reviewer #3 (Remarks to the Author): I am pleased by how the authors addressed my comments. There is only one aspect that I would recommend to improve: In their selective sweep analysis, the authors consider as significant (i.e. putative sweeps) all regions with CLR values falling within the top 5% of the distribution.
The standard in the field, however, is too calculate the statistical threshold for the CLR statistic using neutral simulations obtained from the best demographic model (as obtained by fastsimcoal2 in this study for example).
This is actually relatively easy to do because fastsimcoal2 creates a file called "project_maxL.par" that can be used as input file for fastsimcoal2 to SIMULATE DNA data under your best demographic mdoel for the high and low populations. A bit of scripting may be required to transform the simulated data to the correct input for Sweepfinder2. this procedure will allow the authors to identify "how high" the CLR value can be in their system in the absence of any positive selection/ adaptation minor: what do the author mean by "based on SAS" at line 614