Genetic diversity in temperate and boreal forests tree species has been strongly affected by late Pleistocene climate oscillations [2,3,5], but also by anthropogenic forces. Particularly in Europe, where a long history of human intervention has re-distributed species and populations, it can be difficult to know if a given forest arose through natural regeneration and gene flow or through some combination of natural and human-mediated processes. This uncertainty can confound inferences of the causes and consequences of standing genetic variation, which may impact our interpretation of demographic events that shaped species before humans became dominant on the landscape. In their paper entitled "Genomic data provides new insights on the demographic history and the extent of recent material transfers in Norway spruce", Chen et al.  used a genome-wide dataset of 400k SNPs to infer the demographic history of Picea abies (Norway spruce), the most widespread and abundant spruce species in Europe, and to understand its evolutionary relationship with two other spruces (Picea obovata [Siberian spruce] and P. omorika [Serbian spruce]). Three major Norway spruce clusters were identified, corresponding to central Europe, Russia and the Baltics, and Scandinavia, which agrees with previous studies. The density of the SNP data in the present paper enabled inference of previously uncharacterized admixture between these groups, which corresponds to the timing of postglacial recolonization following the last glacial maximum (LGM). This suggests that multiple migration routes gave rise to the extant distribution of the species, and may explain why Chen et al.'s estimates of divergence times among these major Norway spruce groups were older (15mya) than those of previous studies (5-6mya) – those previous studies may have unknowingly included admixed material . Treemix analysis also revealed extensive admixture between Norway and Siberian spruce over the last ~100k years, while the geographically-restricted Serbian spruce was both isolated from introgression and had a dramatically smaller effective population size (Ne) than either of the other two species. This small Ne resulted from a bottleneck associated with the onset of the iron age ~3000 years ago, which suggests that anthropogenic depletion of forest resources has severely impacted this species. Finally, ancestry of Norway spruce samples collected in Sweden and Denmark suggest their recent transfer from more southern areas of the species range. This northward movement of genotypes likely occurred because the trees performed well relative to local provenances, which is a common observation when trees from the south are planted in more northern locations (although at the potential cost of frost damage due to inappropriate phenology). While not the reason for the transfer, the incorporation of southern seed sources into the Swedish breeding and reforestation program may lead to more resilient forests under climate change. Taken together, the data and analysis presented in this paper allowed inference of the intra- and interspecific demographic histories of a tree species group at a very high resolution, and suggest caveats regarding sampling and interpretation of data from areas with a long history of occupancy by humans.
 Chen, J., Milesi, P., Jansson, G., Berlin, M., Karlsson, B., Aleksić, J. M., Vendramin, G. G., Lascoux, M. (2018). Genomic data provides new insights on the demographic history and the extent of recent material transfers in Norway spruce. BioRxiv, 402016. ver. 3 peer-reviewed and recommended by PCI Evol Biol. doi: 10.1101/402016
 Holliday, J. A., Yuen, M., Ritland, K., & Aitken, S. N. (2010). Postglacial history of a widespread conifer produces inverse clines in selective neutrality tests. Molecular Ecology, 19(18), 3857–3864. doi: 10.1111/j.1365-294X.2010.04767.x
 Ingvarsson, P. K. (2008). Multilocus patterns of nucleotide polymorphism and the demographic history of Populus tremula. Genetics, 180, 329-340. doi: 10.1534/genetics.108.090431
 Lockwood, J. D., Aleksić, J. M., Zou, J., Wang, J., Liu, J., & Renner, S. S. (2013). A new phylogeny for the genus Picea from plastid, mitochondrial, and nuclear sequences. Molecular Phylogenetics and Evolution, 69(3), 717–727. doi: 10.1016/j.ympev.2013.07.004
 Pyhäjärvi, T., Garcia-Gil, M. R., Knürr, T., Mikkonen, M., Wachowiak, W., & Savolainen, O. (2007). Demographic history has influenced nucleotide diversity in European Pinus sylvestris populations. Genetics, 177(3), 1713–1724. doi: 10.1534/genetics.107.077099 "
The revisions made by the authors have added sufficient clarity to the manuscript to address the concerns about how SNPs were identified and filtered. The work described represents a significant technical achievement in application of genomic technologies to conifer population genetics.
I believe the authors have satisfactorily addressed all reviewers' questions and have made appropriate changes in the manuscript. Therefore, I don't have any more questions or comments.
Dear Dr. Lascoux,
Thank-you for your submission to PCI Evol Biol and apologies for the delay in getting back to you. I now have reviews from two experts in the field and invite you to revise your manuscript according to their recommendations. One reviewer mainly asked for clarification regarding how the data were partitioned for the various analyses, and also about the population groupings and their display (in addition to some suggestions for textual changes/clarification). The other reviewer was focused mainly on how the data were processed, the description of associated parameters, and the impact of using an incomplete draft genome on accurately disentangling paralogs.
Once these issues are addressed, I think your paper will make a nice addition to our understanding of the demographic history of trees, and of spruce in particular.
Review of manuscript “Genomic data provides new insights on the demographic history and the extent of recent material transfers in Norway Spruce” by Chen et al. 2018.
In this paper, Chen and co-authors describe a complex demographic history in Norway Spruce, Siberian spruce and Serbian spruce (P. omorika). Population structure in Norway spruce has been shaped by complex migration events, admixture events, and bottlenecks. They also estimated divergence times between the three species, which seem to have occurred earlier than reported in other studies. In general, I found the manuscript interesting and well-written, however I have some concerns in relation to the experimental design and other minor comments.
My main concern is with the experimental design of this work, particularly about the type of markers used. Exome capture requires a priori knowledge of target sequences, therefore the sequences obtained are not “random”. Also, there is generally a heavy emphasis on coding regions, if the SNPs are meant to be used for GWAS or linkage mapping. Of course, non-coding regions can also be targeted but that is not frequent. Alternatively, non-coding regions can be linked to coding regions in targeted sequencing, however these again would not be random. I understand that the 40,018 probes used in this study were designed to cover only coding regions. Then, my question is where did the non-coding SNPs come from? I am mainly worried about ascertainment bias and the absence of markers that could be used as null models in all admixture, SFS and population genomic analyses. Also, it is clear that the SNPs chosen were not “random” therefore may not accurately represent genome-wide patterns. A paragraph about potential limitations and sources of error in this study should be incorporated.
How much of the exome capture design and SNP calling (lines 39-69) was done for this study and how much was previously done (and already reported) in Vidalis et al. 2018; Bernhardsson et al. 2018; and Baison et al. 2018? Please make this distinction in the paper.
From what I understood, only 399k SNPs were used in the admixture and population genomic analyses. If this is the case then “1M SNPs” should be replaced by 399k or other number of SNPs used in the analyses. In the Results section (line 8, page 6), it is mentioned that only 30,000 SNPs were used, so this is confusing. It would be good to make clear the number of SNPs used in each of the analyses.
The admixture plot in Figure 2 is confusing. I wonder if it would look better with K=3 (abies, omorika, obovata). Including both plots (K=5, and K=3) would be more informative, I think. I was surprised to see very little signs of admixture between abies and obovata in the Admixture plot, where other analyses in this paper prove otherwise (lines 52-54, page 8). Any comments on that?
Page 9, line 17- Smaller mutation rates have been reported in more recent literature.
Page 11, line 51-52.- There is convincing evidence that came from several recent papers that trees have lower and not higher mutation rates than annual plants (Lanfear et al. 2013; Bromham et al. 2015; De La Torre at al. 2017). Minor comments
Page 5, Line 48- eliminate “of” after population sizes. Table 1- Kirov population is not in footnote. Figure 1- sample points are not easy to differentiate from each other, and from other group points. I suggest changing shape and/or color. Figure 1- EUFOGEN should read EUFORGEN Lines 15-20, page 7-Although levels of admixture in Kirov and Indigo are mentioned in these lines, these populations are not present in the Admixture plot in Figure 2. Table 3- Please explain what the different parameters names mean?
Chen et al present an extensive set of population genetic analyses of European spruce trees based on genotype data from over 1 million SNP loci, but the manuscript and supplementary material provide little detail on how the sequencing data were filtered and how quality control analyses were conducted on the genotype data, and no discussion of whether the overall conclusions would be different if different filtering and QC thresholds were imposed. A comparable analysis could be conducted with a much smaller genotype dataset, and the reader is left to speculate whether the outcomes would have been different if different criteria were used in the process of calling SNP genotypes and filtering to remove low-quality data.
Conifer genomes are large (typically >15 Gb) and full of repetitive sequences, including not only various classes of mobile elements but also processed pseudogenes that are very similar to existing functional genes in the genome. The use of exome capture sequencing as a method for genotyping therefore requires considerable care in filtering the sequencing data to avoid confounding reads derived from paralogous sequences with those derived from the intended target exons. The authors filtered the sequence data from their samples to minimize the likelihood of detecting paralogous sequences, but some additional information could be provided to allow readers to better judge the degree of rigor used. Some caution is also warranted due to the incomplete nature of the current genome assembly, which limits the ability of the authors to detect (and therefore to exclude) sequence reads derived from multi-copy paralogous sequences. An appropriate strategy to deal with this limitation is to filter the resulting SNP loci carefully to exclude those likely to represent confounded data from paralogous sequences rather than true genotype calls from single-copy loci.
The sequencing reads were aligned to the entire v1.0 draft genome assembly of Norway spruce rather than just the sequences of the target exons, which is good - this allows reads from diverged paralogous sequences that are represented in the draft assembly the opportunity to align to the copy of the sequence to which the read is most similar. The Chen et al manuscript does not point out, however, that the v1.0 draft genome assembly is estimated to include only 60% of the Norway spruce genome, although this fact is highlighted in the abstract of the cited reference by Bernhardsson et al (https://doi.org/10.1101/292151). The sequence reads are derived from the genomes of many individual trees (which may be different from each other as well as from the Z4006 reference individual used for genome sequencing), and are not limited to the subset of the genome represented in the v1.0 assembly. It is llikely, therefore, that paralogous sequences exist in the genome that are not represented in the assembly, and so the genotype data used as the basis for the rest of the manuscript may be a mixture, consisting of true genotypes derived from single-copy sequences and also confounded data derived from paralogous sequences. The relative proportion of this mixture cannot be determined from the data presented in this manuscript, nor even in the original complete dataset, given the fragmented and incomplete state of the v1.0 Norway spruce genome assembly.
One criterion commonly used to address this question is testing for genotype frequencies consistent with expectations based on the assumption of Hardy-Weinberg equilibrium - an excess of heterozygotes can be due to the presence of reads derived from paralogous sequences in the filtered sequence dataset used for SNP genotype callling. There is no mention of such a test in the manuscript, and no discussion of why such a test would (or would not) be suitable for filtering the genotype data prior to conducting the analyses described in the rest of the manuscript. There is also no discussion of the possible impacts on the conclusions drawn if confounded data are present in the genotype calls.
In the spirit of reproducible research, the authors could include the parameters used for read alignment by BWA and extraction of "uniquely-aligned pairs", as these steps are critical to the process of excluding reads from paralogous sequences. The supplementary material provided consists of material supporting the conclusions drawn by the authors regarding the genetic hypotheses, but does not include any details about the process of generating the genotype data on which those hypothesis tests are based. The authors could also provide additional information in supplementary material about steps taken to filter the putative SNP loci used for genotype calling, including the results of testing for consistency with HWE expectations for all SNP loci used in the analyses. Summary data could be included for all candidate loci, including the nature of the deviation from HWE expectations for those that exceed a reasonable threshold, recognizing that some correction for multiple testing will be required due to the large number of loci to be tested. The effects on the number of genotyped loci of imposing different filtering thresholds could also be summarized - such thresholds would include both the depth of read coverage per allele called (the authors used a minimum of two), the minimum fraction of individuals genotyped (the authors used 50%), but could also include deviation from HWE expectations at different FDR thresholds.