Cao et al. BMC Biology (2022) 20:139 https://doi.org/10.1186/s12915-022-01342-y RESEARCH ARTICLE Open Access Chromosome‑level genome assemblies of four wild peach species provide insights into genome evolution and genetic basis of stress resistance Ke Cao1*†, Zhen Peng2†, Xing Zhao2†, Yong Li1†, Kuozhan Liu1, Pere Arus3, Weichao Fang1, Changwen Chen1, Xinwei Wang1, Jinlong Wu1, Zhangjun Fei4 and Lirong Wang1,5* Abstract Background: Peach (Prunus persica) is an economically important stone fruit crop in Rosaceae and widely cultivated in temperate and subtropical regions, emerging as an excellent material to study the interaction between plant and environment. During its genus, there are four wild species of peach, all living in harsh environments. For example, one of the wild species, P. mira, originates from the Qinghai-Tibet Plateau (QTP) and exhibits strong cold/ultraviolet ray environmental adaptations. Although remarkable progresses in the gene discovery of fruit quality-related traits in peach using previous assembled genome were obtained, genomic basis of the response of these wild species to dif- ferent geographical environments remains unclear. Results: To uncover key genes regulating adaptability in different species and analyze the role of genetic variations in resistance formation, we performed de novo genome assembling of four wild relatives of peach (P. persica), P. mira, P. davidiana, P. kansuensis, and P. ferganensis and resequenced 175 peach varieties. The phylogenetic tree showed that the divergence time of P. mira and other wild relatives of peach was 11.5 million years ago, which was consistent with the drastic crustal movement of QTP. Abundant genetic variations were identified in four wild species when com- pared to P. persica, and the results showed that plant-pathogen interaction pathways were enriched in genes contain- ing small insertions and deletions and copy number variations in all four wild relatives of peach. Then, the data were used to identify new genes and variations regulating resistance. For example, presence/absence variations which result from a hybridization event that occurred between P. mira and P. dulcis enhanced the resistance of their puta- tive hybrid, P. davidiana. Using bulked segregant analysis, we located the nematode resistance locus of P. kansuensis in chromosome 2. Within the mapping region, a deletion in the promoter of one NBS-LRR gene was found to involve † Ke Cao, Zhen Peng, Xing Zhao, and Yong Li contributed equally to this work. *Correspondence: wyandck@126.com; wanglirong@caas.cn 1 The Key Laboratory of Biology and Genetic Improvement of Horticultural Crops (Fruit Tree Breeding Technology), Ministry of Agriculture, Zhengzhou Fruit Research Institute, Chinese Academy of Agricultural Sciences, Zhengzhou 450009, People’s Republic of China 5 National Horticulture Germplasm Resources Center, Zhengzhou Fruit Research Institute, Chinese Academy of Agricultural Sciences, Zhengzhou 450009, People’s Republic of China Full list of author information is available at the end of the article © The Author(s) 2022. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://​creat​iveco​mmons.​org/​licen​ses/​by/4.​0/. The Creative Commons Public Domain Dedication waiver (http://​creat​iveco​ mmons.​org/​publi​cdoma​in/​zero/1.​0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data. Cao et al. BMC Biology (2022) 20:139 Page 2 of 17 the resistance by regulating gene expression. Furthermore, combined with RNA-seq and selective sweeps analysis, we proposed that a deletion in the promoter of one CBF gene was essential for high-altitude adaptation of P. mira through increasing its resistance to low temperature. Conclusions: In general, the reference genomes assembled in the study facilitate our understanding of resistance mechanism of perennial fruit crops, and provide valuable resources for future breeding and improvement. Keywords: Comparative genomics, Stress resistance, High-altitude adaptation, Peach Background characters in these studies. To dissect the genetic basis of Peach (Prunus persica) is the third most important of the adaptation to high plateau of P. mira, Wang et al. [13] de deciduous fruit crop in the world, and is widely cultivated novo assembled three high-quality genomes of Tibetan in temperate and subtropical regions, ranking only after Prunus species, including P. mira. The results showed apple and pear. Benefiting from the small genome size, that the expansion of SINE retrotransposons helped peach has been extensively used for comparative and Tibetan Prunus species adapt to the harsh environment functional genomic researches of the Rosaceae family [1]. of the Himalayan plateau by promoting the accumulation Nine years ago, the first high-quality reference genome of beneficial metabolites. sequence of peach was released [2]. Subsequently, the In this study, to dissect the mechanism by which reg- evolution [3, 4], domestication regions [3, 5, 6], and genes ulate peach evolution and species adaptation at the associated with important traits [7, 8] for peach had been genome-wide level, we de novo assembled the genomes successfully dissected. of four wild relatives of P. persica. Comparative genomic Adaptability is the basis of plant survival and repro- analysis showed that P. davidiana was evolved from P. duction, especially under extreme environmental con- mira directly and crossed P. dulcis during its evolutionary ditions. Among the four wild relatives of P. persica, P. history. Then, a nematode resistance gene was identified mira originate from the Qinghai-Tibet Plateau (QTP) in P. kansuensis using the identified genetic variations and show strong adaptability to cold and ultraviolet rays based on new assemblies. Finally, comprehensive char- (UV) resistance. P. davidiana is widely distributed in the acterization of selective sweeps revealed the new mecha- slopes and valleys of mountains in northern China and is nisms underlying the high-altitude adaptability in P. mira. resistant to cold and drought environments. P. kansuensis Compared with previous studies, the genomes of four and P. ferganensis grow in one specific region, Gansu and wild species assembled in the study are of higher quality, Xinjiang province of China, respectively [9]. However, providing more detailed basic data for other researchers genomic basis of the response of these wild species to dif- to study peach evolution and dissect the genetic mecha- ferent geographical environments remains unclear. nism of resistance traits. Genome assembly is a powerful tool to discover elite genes and the method have used frequently in peach in Results recently years. For example, Guan et al. [10] produced Assembly of the genomes of four wild peach species a comprehensive structure variations map for peach to The high-quality genome of P. mira, P. davidiana, P. kan- identify a large 1.67-Mb heterozygous inversion that seg- suensis, and P. ferganensis were assembled using a com- regates perfectly with flat-fruit shape after assembling a bination of PacBio, Illumina, and Hi-C (High-throughput flat-fruit peach cultivar “Rui You Pan 1.” Cao et al. [11] chromosome conformation capture) platforms. Using analyzed the molecular basis of the reduction of volatile the k-mer method (Additional file 1: Table S1; Additional compounds during the evolution of peaches and identi- file 2: Fig. S1), we estimated the size of four genomes to fied the key gene regulating linalool content through be 242.94, 237.29, 238.06, and 237.24 with a heterozy- the assembly of a cultivated peach variety “Chinese gous ratio of 0.76%, 1.10%, 0.56%, and 0.53% in P. mira, Cling” and the association analysis of 256 peach varie- P. davidiana, P. kansuensis, and P. ferganensis, respec- ties. Moreover, some wild related species of peach were tively. Then, a total of 600.37 ×, 574.84 ×, 343.99 ×, also assembled. Tan et al. [12] reported de novo genome and 387.90 × coverage of sequences were generated for assemblies for five species, including P. persica, P. mira, above species and used for genome assembly, respec- and P. davidiana. Next, the author performed an asso- tively (Additional file 1: Table S2). Among the ultimately ciation analysis of 417 peach accessions to identify loci obtained 253.63, 259.27, 253.17, and 261.3 Mb assemblies associated with fruit shape, fruit development period, (Table 1; Fig. S2), 92.29%, 89.46%, 94.03%, and 90.53% of and floral morphology, respectively. In general, research- them were anchored and allocated to eight pseudochro- ers mainly focused on gene identification of fruit quality mosomes for P. mira, P. davidiana, P. kansuensis, and P. Cao et al. BMC Biology (2022) 20:139 Page 3 of 17 Table 1 Genomic sequencing, assembly, and annotation statistic for four wild peach species Species P. mira P. davidiana P. kansuensis P. ferganensis Estimated genome size (Mb) 242.94 237.29 238.06 237.24 Sequencing depth (×) 600.37 574.84 343.99 387.9 Total sequencing data (Gb) 145.86 136.41 91.15 92.02 Assembled genome size (Mb) 253.63 259.27 253.17 261.30 Contig N50 (bp) 25,948,270 22,637,438 24,419,080 26,470,137 Contig N50 (number) 4 4 5 4 Scaffold N50 (bp) 27,933,902 28,549,251 28,885,292 28,247,033 Scaffold N50 (number) 4 4 4 4 GC content (%) 38.29 38.92 38.34 38.84 Percent of repeat sequence (%) 48.38 52.27 49.33 49.49 Predicted gene number 28,519 27,236 26,986 28,587 Annotated gene number 27,613 26,545 26,374 27,803 ferganensis, respectively (Additional file 1: Table S3). The microRNA genes in the four wild peach species (Addi- contig N50 sizes of the four final assemblies were 25.95, tional file 1: Table S10). 22.64, 24.42, and 26.47 Mb, respectively (Table 1), which Genome variations between species were also identified were higher than that of P. persica (250 kb), which pro- by anchoring the four assembled wild genomes onto the duced by the Sanger technology [14]. Compared with the reference genome of P. persica [14]. A total of 547,998- previously published genomes of P. mira and P. davidiana 2,971,987 single nucleotide polymorphisms (SNPs; Addi- with contig N50 of ~ 8.3 Mb and 2.3 Mb, respectively [12], tional file 1: Table S11), 106,436-490,421 small insertions our final assembly showed an ~ 3.1 and 10.0-fold increase and deletions (indels; Additional file 1: Table S12), 3870- in the length of the contig N50. In addition, compared 10,994 large structural variants (SVs; ≥ 50 bp in length; with another released P. mira assembly with a contig N50 Additional file 1: Table S13), and 3344-6840 copy number of 12.14 Mb [13], our P. mira assembly showed an ~ 2.14- variations (CNVs; Additional file 1: Table S14) were iden- fold increase in the length of the contig N50. We also per- tified in the four species (Fig. 1a; Additional file 2: Fig. formed Benchmarking Universal Single-Copy Orthologs S3). A total of 20 SVs that differed in P. mira and P. per- (BUSCO) analysis [15] (Additional file 1: Table S4) and sica was randomly selected (Additional file 1: Table S15) RNA-Seq read mapping rate calculation (Additional to design primers (Additional file 1: Table S16) for ampli- file 1: Table S5) to evaluate the quality of the assembly. fication in cultivar “2010-138” and “Shenzhou Li He Shui The results indicated that more than 95% of the highly Mi” to verify the accuracy of these variations. The results conserved embryophyta genes were completely present showed that among the 20 variants, only the electropho- in the four genomes, and more than 85% RNA-Seq reads resis results of two primers might be inconsistent with could be mapped to the corresponding genomes. Subse- the sequencing results (Additional file 2: Fig. S4). Mean- quently, repeat sequence were annotated and the results while, IGV software (www.​broad​insti​tute.​org/​igv) was showed that long terminal repeat (LTR) retrotranspo- used to observe the above 20 variant sites manually. We sons made up the majority of the transposable elements analyzed the function of the genes which comprised dif- (TEs), comprising 35.61%, 42.44%, 41.75%, and 41.86% of ferent variations (Additional file 2: Fig. S5-8) and found the P. mira, P. davidiana, P. kansuensis, and P. ferganensis that plant-pathogen interaction pathways were enriched genome, respectively (Additional file 1: Table S6). Gene in genes containing small indels and CNVs in all four prediction and annotation were performed resulting in wild relatives of peach (p = 0.01 ~ 1.74e−22), showing that 28,519, 27,236, 26,986, and 28,587 protein-coding genes responding with environmental stress was very impor- in P. mira, P. davidiana, P. kansuensis, and P. ferganensis, tant for the speciation of four species. respectively (Additional file 2: Fig. S2; Additional file 1: Table S7 and 8). About 96.82-97.73% of the protein-cod- Peach genome evolution and species divergence ing genes of the four wild species could be functionally Combining the genome of P. mira, P. davidiana, P. kan- annotated (Additional file 1: Table S9). In addition, we suensis, P. ferganensis, and P. persica, all of the genes also identified 8171-15,770 ribosomal RNA, 1226-1564 in the genomes of these five peach species were classi- transfer RNA, 514-753 small nuclear RNA, and 563-693 fied into 22,747 families on the basis of the homology Cao et al. BMC Biology (2022) 20:139 Page 4 of 17 Fig. 1. Variation screening of four wild species and evolutionary analysis of peach genome. a SNPs, small indels, SVs, and CNVs in P. mira, P. davidiana, P. kansuensis, and P. ferganensis compared to P. persica. b Core and dispensable gene families of four wild peaches and P. persica. c Gene Ontology annotation of genes specific in each species. d Estimation of divergence times of 15 species and identification of gene family expansions and contractions. Numbers on the nodes represent the divergence times from present (million years ago, Mya). MCRA, most recent common ancestor. e Distribution of Ks (synonymous mutation rate) values of orthologous genes between six genomes of the Prunus species and strawberry of their encoded proteins (Additional file 2: Fig. S9). pathways were enriched in specific gene families of P. Presence/absence variations (PAVs) analysis of gene davidiana (Fig. 1c). families revealed 8,138 (35.78%) dispensable gene-fam- A total of 830 single-copy gene families in the 15 spe- ilies distributed across all five genomes, and 368, 303, cies were identified (Additional file 2: Fig. S10) and used 303, 388, and 763 families specific to each of the above to construct a phylogenetic tree of P. persica and its species, respectively (Fig. 1b). Single-organism meta- wild related species as well as other representative plant bolic process (p = 1.36e−6) and oxidation-reduction species (Fig. 1d). Based on the known divergence time (p = 2.53e−15) were found to be enriched in gene fami- between Arabidopsis thaliana and strawberry, the age lies specific to P. mira. Cellular macromolecule biosyn- estimate for the split of the ancestor of P. mume, P. arme- thetic process (p = 0.05) and gene expression (p = 0.01) niaca, P. salicina, and the common ancestor of peach and almond was around 21.1 million years ago (Mya), later Cao et al. BMC Biology (2022) 20:139 Page 5 of 17 than that in a previous report, presumably 44.0 Mya [16]. that PAV may explain the multi-resistance of P. davidi- Hereafter, the divergence time of P. dulcis and P. mira was ana to diverse stresses and the deletions in the promoter about 13.3 Mya, which was obviously earlier than that of region could enhances the resistance of P. kansuensis to Yu et al. [4] and Alioto et al. [17] who found the diver- nematodes and P. mira to low temperature by regulating gence time of the two species was 4.99 and 5.88 Mya, the expression of target genes, respectively. Specific cases respectively. Furthermore, we found that P. mira split are described below. with the common ancestor of P. davidiana, P. kansuen- sis, and P. ferganensis approximately 11.5 Mya. The event Interspecific hybridization is essential for P. davidiana occurred around the drastic crustal movement of QTP to obtain the multi‑resistance traits [18] where P. mira originated. P. davidiana is a very important wild species of P. persica Lineage-specific gene family expansions may be associ- because of its resistance to various abiotic stresses such ated with the emergence of specific functions and physi- as drought and cold and biotic stresses such as aphids, ology [19]. With the genome evolution of wild peach which may be of great significance in the use of resist- relatives, the number of expanded gene families was ance breeding systems in the future [22]. Here, we first decreased from P. mira to P. kansuensis while increased screened the R genes from four wild species and found in P. ferganensis or P. persica compared with that in a total of 295, 341, 310, and 322 genes were identified the most recent common ancestor (MCRA; Fig. 1d). in P. mira, P. davidiana, P. kansuensis, and P. ferganensis We found that the expanded gene families were highly genomes, respectively (Additional file 1: Table S17). The enriched with those related to aminobenzoate degrada- largest number of R genes in P. davidiana might explain tion (or arginine and proline metabolism, styrene degra- its strong and multiple resistances to different pathogens, dation) in P. mira (p = 5.53e−21; Additional file 2: Fig. S11) such as aphid, Agrobacterium tumefaciuns, etc. We found which originated in QTP, same as in Polylepis tarapacana that most of the disease resistances QTLs/genes were [20] and Populus kangdingensis [21] where lived in high located in genome regions containing candidate R genes altitude regions. (Additional file 2: Fig. S13). To validate the speciation events of peach species, four- Next, we want to analyze the source of these R genes. fold synonymous third-codon transversion (4DTv) rates Considering that the k-mer frequency of P. davidi- were calculated and the results showed that all 4DTv val- ana presents a characteristic of a high heterozygosity ues (Additional file 2: Fig. S12) among paralogs in four genome (Additional file 2: Fig. S1b), we speculated that wild peach species peaked at around 0.50 to 0.55, consist- there was a hybridization event in the evolutionary his- ent with the whole-genome triplication event (γ event) tory of P. davidiana. To further verify this inference, shared by all eudicots and indicating that no recent one hundred and twenty-six accessions of wild peach whole-genome duplication occurred. A peak 4DTv value species (Additional file 1: Table S18; accessions 101- at around 0 for the orthologs between four wild related 126) were resequenced, and the obtained sequences species of peach and P. persica highlighted a very recent (Additional file 1: Table S19) were aligned to the P. diversification of Prunus species [16]. To estimate the mira genome to obtain a total of 839,431 high-quality time of species divergence of the four wild species, we SNPs. Using these SNPs, we constructed a phyloge- calculated the Ks (rate of synonymous mutation) values netic tree and found that cherries, plums, and apricots of orthologous genes between these species. As shown in had a close relationship with P. persica and its wild Fig. 1e, the peaks at a Ks mode for orthologs between P. related species. And during the latter category, P. dul- persica-P. dulcis was higher than that of other peaches. In cis, P. tangutica, P. davidiana, and P. mira formed a addition, among the wild related species of P. persica, the small subgroup (Fig. 2a). Principal component analysis similar Ks value of P. mira, P. davidiana, and P. kansuen- (PCA) also demonstrated that P. dulcis or P. tangutica sis showed significantly higher than that of P. ferganensis, is located between P. mira and P. davidiana in the dia- indicating similar results with phylogenetic tree. gram of PC2-PC3 (Fig. 2b). We then phenotyped the stone streak of different species and found that abun- Variations screening provided identification of new genes dant dot streaks were present in the P. davidiana but regulating resistance absent in the more recent species, P. kansuensis. How- To assess the advantage of genome assemblies of peach ever, the dot streaks could trace back to the ancient using five species, we analyzed the molecular mechanism species of this family, such as P. dulcis (Fig. 2c). Accord- for regulating the resistance traits of different species one ing to the above evidence, we proposed that P. davidi- by one. According to previous reports, P. fergenansis was ana showed intermediate genomic characteristics recognized as only a geographic type of P. persica [3], so between P. mira and P. dulcis and might originate from it was not discussed in the following analysis. We found the cross between these two species. A previous study Cao et al. BMC Biology (2022) 20:139 Page 6 of 17 Fig. 2. The evolution of P. davidiana. a Phylogenetic tree of 126 peach accessions which included P. persica and its wild relative species. b Principal component analysis (PCA) of above accessions. c Stone steak morphology of P. dulcis, P. mira, P. davidiana, P. kansuensis, and P. ferganensis. d Percent of P. davidiana-specific contigs covered by reads from different Prunus species. e Regional collinearity between P. davidiana with P. dulcis or with P. mira genomes. The detailed location of these three regions were at Chr. 1: 4.79 Mb, Chr. 1: 5.08 Mb, and Chr. 2: 7.50 Mb in P. davidiana genome. Green and blue boxes indicate the positive and negative direction of genes. “P” and “A” letters showed the gene was presence and absence in P. davidiana genome [4] reported that P. tangutica has a closer relationship to P. dulcis, and next was P. tangutica (Fig. 2d), which with P. davidiana. However, since no reference genome again proved that the introgression in P. davidiana is available for P. tangutica, a direct comparison of came from P. dulcis or P. tangutica. In addition, all the genes from P. davidiana and P. tangutica is not feasible. genes in P. davidiana were aligned with those from Therefore, we first aligned the sequences of P. davidi- the putative ancestors who had an assembled genome, ana to its putative ancestor, P. mira, and assembled including P. mira and P. dulcis (Additional file 2: Fig. the unmapped sequences to obtain a partial reference S14). The genes were then assigned as originating from genome. Genome resequencing data of different species a specific species if the highest score was observed were then aligned to this partial genome. We found that from the alignments of orthologous genes between P. the highest mapping rate was an accession belonging davidiana and the other species. The results showed Cao et al. BMC Biology (2022) 20:139 Page 7 of 17 Fig. 3. Identification of nematode resistance genes in P. kansuensis. a Bulked segregation analysis to locate the nematode resistance gene using a backcross population. b The Indels in the promoter and mRNA regions of one NBS-LRR genes on Chr. 2. “Ins” and “Del” indicate an insertion and deletion, respectively. And the number after the variation type indicates the length (bp) of Indels. c Expression of Prupe.2G053600 in two accessions (“Hong Gen Gan Su Tao 1#” and “Bailey”) inoculated with nematode. e Functional validation of Prupe.2G053600 through analysis of transgenic tomato plants expressing Prupe.2G053600 under nematode treatment. e Promoter activity assay. Promoters with different lengths were fused to the GUS gene in plasmid pBI101. GUS was dyed, and its activity was measured using protein extracts of tobacco. pBI101 and pBI121 indicate negative control and positive control, respectively that although about 50.74% of the genes (13,719) were followed by P. dulcis (23.05%). Therefore, we believed specific to P. mira, 26.87% (7,266) were thought to be that the cross between P. mira and P. dulcis or P. tangu- derived from P. dulcis, and no more than 8% from other tica enhanced the adaptation of P. davidiana when it was species. spread to new environments. To further study the evolutionary events leading to the genome structure of P. davidiana, we investigated the An indel in the promoter of an NBS‑LRR gene regulates chromosome-to-chromosome relationships based on the nematode resistance in P. kansuensis 16 (P. davidiana versus P. mira, including 21,186 genes) P. kansuensis shows high resistance to root-knot nem- and 28 (P. davidiana versus P. dulcis, including 18,480 atodes [23]. To discover candidate causal sequence genes) identified syntenic blocks. The mosaic syntenic variations involved in the resistance, we constructed patterns again demonstrated that P. davidiana might a backcross ­(BC1) population from the cross between have arisen during the evolutionary process of P. mira “Hong Gen Gan Su Tao 1#” (P. kansuensis) and a culti- but with a cross of other species (Additional file 2: Fig. vated peach “Bailey” (P. persica). “Hong Gen Gan Su S15). Regional collinearity between P. davidiana with P. Tao 1#” harbored high resistance to root-knot nema- dulcis or P. mira indicated that some regions contained todes (Meloidogyne incognita), whereas other acces- PAVs in P. davidiana genome showed a high structural sions including P. mira, P. davidiana, and P. ferganensis similarity with them in P. dulcis (Fig. 2e). Although P. dul- all showed low resistance to M. incognita [23]. Using this cis is mainly distributed in Georgia, Azerbaijan, Turkey, ­BC1 population, a nematode resistance locus was Syria, and Xinjiang province (China), it can also be found mapped at the top region of chromosome 2 (5.0-7.0 Mb) in Sichuan province of China where some P. mira grow in based on bulk segregation analysis (BSA) (Fig. 3a). The this region at the same time, indicating the hybridization most significant locus was located at 6,671,644 bp of event is highly possible (Additional file 2: Fig. S16). chromosome 2. In other species, the R genes to nema- We finally analyzed the origin of R genes in P. davidi- todes generally encode proteins containing the NBS- ana (Additional file 2: Fig. S12) and found 57.25% of LRR domain [24]. Moreover, previous reports showed all R genes were categorized to originate from P. mira, that the expression level of R gene which conferring Cao et al. BMC Biology (2022) 20:139 Page 8 of 17 root-knot nematode resistance was significantly higher difference of cis-elements in the promoter region of the in the roots of resistant plants than that of susceptible resistant accession “Hong Gen Gan Su Tao 1#” and the one and its value drastically increased when the plants susceptible one “Bailey.” The result showed that although were inoculated with nematodes [25, 26]. Therefore, the former accession was less than the latter by four transcriptome data [27] of two cultivars (resistant and TATA-box Elements, three more CAAAT-box elements susceptible to nematode infection, respectively) were in the former than that of the latter may be the main rea- analyzed again. We found that three of 19 expressed son for its higher activity in regulating the expression of NBS-LRR genes (Prupe.2G046000, Prupe.2G053600, the target gene (Additional file 1: Table S21). and Prupe.2G055500) have higher expression levels in “Hong Gen Gan Su Tao 1#” than those of “Bailey” dur- Combined RNA‑seq and selection sweep analysis to dissect ing most inoculation periods (Fig. S17). However, only the high‑altitude adaptability in P. mira Prupe.2G053600 was induced by root-knot nematodes P. mira is an attractive model for studying high-altitude in 6 or 12 h of resistant accession than that of 0 h (Addi- adaptability of perennial plants because it originated in tional file 2: Fig. S17; Fig. 3c). In the regions, we further the QTP in China. The region has an average elevation of focused on the genome variations mainly including ∼4000 m above the sea level, and the oxygen concentra- small indels and SVs which occurred in the 24 R genes tion is ∼40% lower and UV radiation is ∼30% stronger and their promoter. We found that there were 92 vari- than those at the sea level [28]. Up to date, knowledge ations which only present in P. kansuensis but not in on the mechanism of high-altitude adaptability has been other species (Additional file 1: Table S20). Among them, reported in pig [29], yak [30], human [28, 31], snake [32], Prupe.2G053600, which comprised a large deletion (Chr. hulless barley [33], and the herbaceous plant Crucihima- 2: 6,212,387..6212,422 bp) in the promoter was further laya himalaica [34]. However, little is known in perennial analyzed (Fig. 3b). crops about the genetic basis of response to harsh condi- Then, the coding sequence of this gene was inserted tions, such as low temperature and high UV radiation in into the plant expression vectors and transformed into high-altitude environments. tomato (cv. Micro-Tom). The transgenic lines were vali- Population genomic analyses were performed to ana- dated by genomic polymerase chain reaction (PCR) and lyze the high-altitude adaptability of P. mira. A total of qRT-PCR. We found the most transgenic lines showed 49 accessions (Additional file 1: Table S18; Accessions remarkable nematode resistance with less root knots 127-175) belonging to the species with an altitude rang- compared to control plants (Fig. 3d; Additional file 2: Fig. ing from 2293 to 4390 m were resequenced to an average S18). depth of 24.96×, and the sequencing reads were aligned We validated this promoter deletion (35 bp) and found to the P. mira genome to identify a total of 1,187,611 that it co-segregated with resistant phenotype of the SNPs. Based on the altitude of these samples (Additional seedlings in the ­BC1 population. To identify the active file 2: Fig. S20), the phylogenetic and structure analy- region of the promoter in Prupe.2G053600, we amplified ses using the identified SNPs, 12 accessions thought to 161, 282, 693, 1497, and 2063 bp of the 5′ flanking region have not corresponded to their altitude categories or of the gene in “Hong Gen Gan Su Tao 1#’ and linked the in an admixture subgroup between high and low alti- amplified products with the β-glucuronidase (GUS) cod- tude subgroups were excluded in the following analysis. ing sequence to transiently transform into Nicotiana Finally, ten accessions were classified into a high-altitude tabacum. Leaves from the transgenic lines were ana- subgroup and 27 into the low-altitude subgroup (Addi- lyzed for GUS activity by histochemical GUS staining tional file 2: Fig. S21). We calculated and compared the and GUS quantitative enzyme activity determination. nucleotide diversity (π; Fig. 4a), Tajima’s D (Fig. 4b), and The lines carrying the various Prupe.2G053600 promot- fixation index (FST, Fig. 4c) values using SNPs across the ers displayed remarkable but lesser GUS activity in com- genome of high- and low-altitude groups, resulting in the parison with the CaMV35S transformed one (pBI121 identification of selective sweeps of a total of 1.9 Mb and vector). An increase in GUS expression was observed containing 231 genes (Additional file 1: Table S22). These with promoters longer than 693 bp, indicating the dele- genes were mostly involved in resistance to a series of tion (310 bp ahead of start codon) could drive the expres- stresses, such as cold, UV light, and DNA damage (Addi- sion of Prupe.2G053600 (Fig. 3e). We also constructed a tional file 1: Table S23; Additional file 2: Fig. S22). mutation by recovery of 35 bp deletion to the promoter of Furthermore, using the young seedlings of P. mira “Hong Gen Gan Su Tao 1#” peach and linked it with GUS treated under low temperature (0 °C) and UV-light gene. The result found that the GUS activity of this muta- (~ 100 μW/cm2) for 10 h, we found that most genes tion was significantly inhibited compared with control presented stronger induction by UV than by low tem- (Additional file 2: Fig. S19). In addition, we analyzed the perature based on the RNA-Seq data (Fig. 4d). Gene Cao et al. BMC Biology (2022) 20:139 Page 9 of 17 Fig. 4. Selective regions associated with high-altitude adaptation in P. mira. a–c Domestication signals in accessions originating in high-altitude region compared to those in low-altitude. The signals were defined by the top 5% of πratio (a), Tajima’s D (b), and FST values (c). d Distribution of expression of genes induced by low temperature and UV of P. mira. Grey dots indicate the background genes and red dots indicate selective genes associated with high-altitude adaptation. e Detailed π ratio and FST values in the genome region of the candidate gene, Pmi02g3025 (pointed by the dashed line), which was substantially induced by low temperature. f Detailed Tajima’s D in the genome region of the candidate gene Pmi02g3025. g A. thaliana plants expressing Pmi02g3025 gene (OE) and the control (WT) treated with low temperature. h Genotypes (K indicates G/T) of a variation (Chr. 2: 28,612,439 bp) located at the promoter of Pmi02g3025 in accessions from different altitude regions expression analysis showed that genes encoding C-repeat under selection by altitude (Fig. 4e, f ). In addition, we binding transcription factor (CBF) and late embryogen- heterologously expressed the Pmi02g3025 gene in Arabi- esis abundant protein (LEA) are strongly expressed after dopsis and treated the ­T2 transgenic plants under 0 °C being induced by low temperature, while the enrich- for 48 h. The transgenic Arabidopsis seedlings showed ment of thaumatin-like protein plays a key role in the increased resistance to low temperature compared to the UV response in P. mira. One CBF gene, Pmi02g3025, wild type (Fig. 4g). We found that the gene expression of showed more than 240-fold induction of expression by Pmi02g3025 at 2 h after cold induction were higher than cold and about 6-fold induction by UV-light. According 0 h in the transgenic plants (Additional file 2: Fig. S23). to the FST and Tajima’s D values, this gene was indeed And the change of proline contents in wild-type and Cao et al. BMC Biology (2022) 20:139 Page 10 of 17 transgenic plants after cold induction also proved trans- genomic analysis using different accessions of P. mira genic ones have a higher resistance to low temperature originating from low- and high-altitude regions, we (Additional file 2: Fig. S24). found that the CBF genes play an important role in high- Based on resequencing data, we found seven SNPs and altitude adaptation of P. mira. However, analyzing the an indel in the promoter regions having a strong asso- distribution of the 19 bp deletion in varieties with differ- ciation with the phenotype, such as a SNP (Fig. 4 h) and ent low temperature resistance indicated that the varia- a 19 bp deletion (Additional file 2: Fig. S25) located at tion was specific in P. mira which lived in high altitude 1222 bp and 212 bp upstream of the start codon, respec- region (Additional file 1: Table S25). tively. The minor allele frequency of the former was 0.72 In addition, we found that gene families involved in and 0 in low and high altitude groups, respectively. And vitamin B6 metabolism were significantly expanded in the resequencing results showed that all high-altitude P. mira and P. ferganensis. Further analysis indicated that varieties had a 19 bp deletion in promoter of Pmi02g3025. among the above expanded gene families, a total of 148 All these results indicated that an obvious divergence genes encoding pyridoxal 5′-phosphate synthase were existed between high and low altitude groups. Subse- identified in P. mira and the corresponding gene numbers quent experiments also proved the relationship between decreased to 6, 19, and 92 in P. davidiana, P. kansuensis, the deletion and promoter activity. and P. ferganensis, respectively. Enzyme encoded by this gene family catalyzes the hydrolysis of glutamine to glu- Discussion tamate and ammonia as part of the biosynthesis of pyri- In peach, a high-quality reference genome of P. per- doxal 5’-phosphate [36]. It was reported that pyridoxal sica was released and has since widely used as a valu- 5’-phosphate was an essential compound to improve salt- able resource for effectively mining candidate genes for stress tolerance in Arabidopsis [37], wheat [38] and rice important traits [14]. However, this genome sequence [39]. The component also responded to several stresses alone is not adequate to uncover wild-specific sequences such as chilling, UV radiation, and drought in tobacco which might have been lost during domestication or arti- [40]. However, there were few reports on the involve- ficial selection [35]. In this study, we de novo assembled ment of the vitamin B6 metabolic pathway in altitude the genomes of four high-quality wild relatives of P. per- adaptation in plants. Interestingly, the researcher com- sica. We analyzed the divergence time of peach and its pared gene families in the genome of Triplophysa bleek- relatives and found that the divergence time of Prunus eri, an endemic fish inhabiting high-altitude regions, to genus from others was 56.6 Mya, later than that in other those fish species living in non-QTP regions and found studies, 66.2 Mya [16]. However, the divergence time that the expanded gene families in former species were between the wild relative of peach and almond was 13.3 primarily enriched in categories of vitamin B6 metabo- Mya, which was significantly earlier than previous report, lism [41]. Furthermore, after sequencing the 16S rRNA 4.99 Mya [4]. We believe that the more number of species genes of oral [42] or skin [43] microbiomes which sam- analyzed in our study and the differential software used pled from Tibetans at different altitude in Tibet, the vita- in previous two studies and our (BEAST vs PAML soft- min B6 metabolic pathway was all found to be enriched ware) was the main reasons for the differences between in the higher altitude group. Therefore, the expansion of different reports. gene involved in vitamin B6 metabolism in P. mira might Furthermore, using this large-scale comprehensive explain its high-altitude adaptability. dataset, millions of genomic variations including SNPs In recent years, we and other groups have made impor- and SVs were identified. In addition, hundreds of spe- tant progress in peach evolution, divergence, and excel- cific gene families in each of the wild peach species were lent genes identification. However, the main ideas of also identified. The above gene sets represent a useful these studies were different from that of this study. For source for in-depth functional genomic studies includ- example, by performing large-scale resequencing of 10 ing the identification of a nematode resistance gene from wild and 74 cultivated peach varieties, we show that P. P. kansuensis. We identified the genotype of the 35 bp mira is the most primitive wild species of peach, and P. deletion located in the promoter of Prupe. 2G053600 ferganensis can be classified as a subgroup of P. persica in other accessions and found that it has no correlation [3]. Obviously, gene introgression analysis was not per- with resistance (Additional file 1: Table S24). In addition, formed in the above study. However, in the study, we the putative interspecific hybridization of P. davidiana constructed the phylogenetic trees of different species identified in this study expands our understanding of and speculated that the specific genes in P. davidiana the evolutionary path of the species and the molecular may be derived from the hybridization of P. mira and mechanisms underlying multiple resistance traits. More- almond (Fig. 2). To verify this speculation, P. kansuen- over, based on expanded gene families and comparative sis and P. davidiana genomes were firstly aligned to P. Cao et al. BMC Biology (2022) 20:139 Page 11 of 17 dulcis. The results showed that the mapping rate of the (2010-138, Zhou Xing Shan Tao 1#, Hong Gen Gan Su latter (65.15%) was slight higher than that of the former Tao 1#, and Ka Shi 1#) were sequenced and de novo (63.73%). At the same time, these two genomes were also assembled. Second, 175 peach accessions were selected aligned with P. mira and found that the rate was 68.52% for genome resequencing. The above accessions were in the latter, which was also slight higher than the former conserved in the National Germplasm Resource Reposi- (66.91%), confirming that P. kansuensis was more evolu- tory of Peach at Zhengzhou Fruit Research Institute, Chi- tional than P. davidiana during peach evolution. Differ- nese Academy of Agricultural Sciences, China. Among ent from this result, Yu et al. [4] found that an ancient them, most accessions belonging to P. mira were sampled introgression that existed between P. mira and the com- from Tibet with different altitudes. Third, one B ­ C1 popu- mon ancestor of P. kansuensis and P. persica might lead lation was constructed between “Hong Gen Gan Su Tao to the development of fruit edibility in peach. Genome 1#” (P. kansuensis) and “Bailey” (P. persica) to identify assembly and identification and functional verification QTLs linked to nematode resistance. Resistance to nema- of several excellent genes are not involved in the above tode in this ­BC1 population was evaluated previously by studies. Here, we first completed the assembly of four our group [3]. Genomic DNA was extracted from young wild species of P. persica and identified the nematode- leaves using the Plant Genomic DNA kit (Tiangen, Bei- resistant genes of P. kansuensis and the altitude adapta- jing, China). bility-related genes of P. mira. However, there are still Moreover, different samples were used for RNA significant differences between the discovery of impor- sequencing (RNA-Seq). First, young leaves, mature fruits, tant genes in the study and that reported before [8] which seeds, phloem, and roots (obtained through asexual mainly focused on fruit quality traits. reproduction) of P. persica (Shang Hai Shui Mi), P. fer- Moreover, when analyzing the genetic variations in ganensis (Kashi 1#), P. kansuensis (Hong Gen Gan Su key genes controlling various adaptive traits of different Tao 1#), P. davidiana (Hong Hua Shan Tao), and P. mira peach species, we found that the source of variations was (2010-138) were collected. Second, roots of “Hong Gen diverse, including gene expansion, hybridization, and Gan Su Tao 1#” and “Bailey” infected with Meloidogyne mutation. In previous studies, genome expansion and incognita for 3, 6, 9, 12 h were collected for RNA-Seq contraction of repetitive elements within the genome analysis which was reported before [27]. was reckoned to fuel evolutionary diversity in many plant species [44]. And adaptive introgression or sex might Genome sequencing of wild peach species represent more efficient management option to con- The genomes of four wild species of peach were serve evolutionary potential in a changing environment sequenced using different platforms including PacBio [45, 46]. In view of the different origin of peach species, Sequel and Illumina according to the manufacturers’ pro- we speculated that the different sources of the variations tocols. Library construction and sequencing were per- were related to the intensity of the environmental change formed at Novogene Bioinformatics Technology Co., Ltd of the species. But, more evidences are needed to prove (Tianjin, China). For short-read sequencing, two short- the precise mechanism underlying plant adaptation to insert libraries (230 bp and 500 bp) and 4 large-insert changing environment during evolution. libraries (2 kb, 5 kb, 10 kb, and 20 kb) were constructed for P. mira and P. davidiana, while two short-insert Conclusions libraries and 2 large-insert libraries (2 kb and 5 kb) were In summary, the study reported four high-quality constructed for P. kansuensis and P. ferganensis. These genome assemblies and found different wild species libraries were sequenced on an Illumina HiSeq X Ten contained a large number of specific genetic variations, platform (Illumina, San Diego, CA). including PAV and indels. To demonstrate their utility for For single molecule real-time (SMRT) sequencing, basic genetics research and plant breeding applications, a 15-kb library was constructed and sequenced on the the role of some genetic variations in the resistance for- PacBio Sequel II platform (Pacific Biosciences, Menlo mation in wild peach species were dissected. We believe Park, CA). On the basis of subreads, the ccs software these genomic resources will increase our understanding (https://​github.​com/​Pacif​i cBio​scien​ces/​ccs) was used of the evolution of fruit crops and be available to acceler- to generate HiFi reads. The parameter was set as below: ating crop improvement in future. min-passes = 3, min-rq = 0.99. For chromosome conformation capture (Hi-C) Methods sequencing, leaves fixed in 1% (vol/vol) formaldehyde Plant materials were used for library construction. Cell lysis, chromatin In this study, different samples were used for DNA digestion, proximity-ligation treatments, DNA recovery, sequencing. First, genomes of four wild accessions and subsequent DNA manipulations were performed as Cao et al. BMC Biology (2022) 20:139 Page 12 of 17 previously described [47]. DpnII was used as the restric- Gene prediction and functional annotation tion enzyme in chromatin digestion. The Hi-C library Gene prediction was performed using a combination of was sequenced on the Illumina HiSeq X Ten platform homology, ab initio and transcriptome based approaches. (Illumina, San Diego, CA) to generate 150 bp paired-end For homology-based prediction, protein sequences from reads. P. persica, Pyrus bretschneideri, P. mume, Malus domes- tica, and Fragaria vesca (Genome Database for Rosaceae; https://​www.​rosac​eae.​org) and Vitis vinifera (http://​www.​ RNA‑Seq data generation genos​cope.​cns.​fr/​exter​ne/​Genom​eBrow​ser/​Vitis/) and To assist protein-coding gene predictions, we performed Arabidopsis thaliana (https://​www.​arabi​dopsis.​org) were RNA-Seq using five different tissues for each species downloaded and aligned to the peach assemblies. Augus- of 15 years old, and for each sample, three independ- tus [54], GlimmerHMM [55], and SNAP [56] were used ent biological replicates were generated. Total RNA was for ab initio predictions. For transcriptome-based pre- extracted with the RNA Extraction Kit (Aidlab, Beijing, diction, RNA-Seq data derived from root, phloem, leaf, China), following the manufacturer’s protocol. RNA- flower, and fruit were mapped to the assemblies using Seq libraries were prepared with the Illumina stand- HISAT2 software [57] and assembled into the transcripts ard mRNA-seq library preparation kit and sequenced using Cufflinks (version 2.1.1) with a reference-guided on a HiSeq 2500 system (Illumina, San Diego, CA) with approach [58]. Moreover, RNA-Seq data were also de paired-end mode. novo assembled using Trinity v2.0 [59] and open reading frames in the assembled transcripts were predicted using Genome assembly of wild peach species PASA [60]. Finally, gene models generated from all three The genome sizes of the four wild peach species were approaches were integrated using EvidenceModeler [60] estimated by K-mer analysis. The occurrences of K-mer to generate the final consensus gene models. with a peak depth were counted using Illumina paired- The predicted genes were functionally annotated by end reads, and genome sizes were calculated accord- comparing their protein sequences against the NCBI ing to the formula: total number of K-mers/depth at the non-redundant (nr), Swiss-Prot (http://​www.​unipr​ot.​ K-mer peak, using JELLYFISH 2.1.3 software [48] with K org/), TrEMBL (http://​www.​unipr​ot.​org/), Kyoto Ency- set to 17. Firstly, all HiFi reads which obtained from the clopedia of Genes and Genomes (KEGG, http://​www.​ above procedure were de novo assembled using Hifiasm genome.​jp/​kegg/), InterPro, and Gene ontology (GO) (v0.13-r308) [49]. Then, Hi-C data were used to correct databases. superscaffolds and cluster the scaffolds into pseudochro- tRNAscan-SE [61] was used with default parameters mosomes through ALLHiC [50] pipeline. to identify tRNA sequences in the genome assemblies. To evaluate the quality of the genome assemblies, rRNAs in the genomes were identified by aligning the ref- we performed BUSCO v4.0.0 [15] analysis on the four erence rRNA sequence of relative species to the assem- assembled genomes with the 1,614 conserved plant sin- blies using BLAST with E-values < 1 ­ e−10 and nucleotide gle-copy orthologs. Moreover, RNA-Seq data with dif- sequence identities > 95%. Finally, the INFERNAL v1.1 ferent tissues were used to validate the results of gene (http://​infer​nal.​janel​ia.​org/) software was used to com- prediction (Additional file 1: Table S5). pare the genome assemblies with the Rfam database (http://​rfam.​xfam.​org/) to predict miRNA and snRNA sequences. Repetitive element identification A combined strategy based on homology alignment and de novo search was used to identify repeat elements in the four wild peach genomes. For de novo prediction of Genome alignment and collinearity analysis TEs, we used RepeatModeler (http://​www.​repea​tmask​ Orthologous genes within the P. mira and P. persica er.​org/​Repea​tMode​ler.​html), RepeatScout (http://​www.​ genomes were identified using BLASTP (E value <1e−5), repea​tmask​er.​org/), Piler [51], and LTR-Finder [52] and MCScanX [62] was used to identify syntenic blocks with default parameters. For alignment of homologous between the two genomes. The collinearity of the two sequences to identify repeats in the assembled genomes, genomes was then plotted according to the identified we used RepeatProteinMask and RepeatMasker (http://​ syntenic blocks. www.​repea​tmask​er.​org) with the repbase library [53]. Four wild peach genomes were aligned to the P. per- TEs overlapping with the same type of repeats were inte- sica genome using LASTZ [63] with the parameters of grated, while those with low scores were removed if they “M=254K=4500 L=3000 Y=15000 --seed=match 12 overlapped more than 80 percent of their lengths and --step=20 --identity=85.” In order to avoid the inter- belonged to different types. ference caused by repetitive sequences in alignments, Cao et al. BMC Biology (2022) 20:139 Page 13 of 17 RepeatMasker (http://​www.​repea​tmask​er.​org) and amplified products was transformed into A. tunefaciens RepBase library [53] were used to mask the repetitive (GV1301) cells and collected and resuspended in infiltra- sequences in genomes of P. persica and four wild species. tion buffer and then transformed into 6-week-old tobacco The raw alignments were combined into larger blocks leaves using sterilized syringes. The transiently trans- using the ChainNet algorithm implemented in LASTZ formed tobacco plants were grown in a growth chamber [63]. for 48 h and the infection sites were cut to measure GUS activity as previously described in Jefferson et al. [70]. Variation identification The pBI121 vector was used as a positive control. We identified SNPs and small indels (< 50 bp) between the four wild and reference peach genomes using SAMtools Transgenic analysis (http://​samto​ols.​sourc​eforge.​net/) and LAST (https://​git- The full-length open reading frame of the lab.​com/​mcfri​th/​last) with parameters “-m20 -E0.05” and Prupe.2G053600 gene was amplified through PCR using SNP and indel filtering criteria “minimum quality = 20, cDNA synthesized from RNA that was isolated from minimum depth = 5, maximum depth = 200.” SVs were root of the “Hong Gen Gan Su Tao 1#” (P. kansuensis). identified from genome alignments by LAST with The amplified product was cloned into the pEASY vec- parameters “-m20 -E0.05.” The software has been used tor driven by the cauliflower mosaic virus (CaMV) 35 S in many studies [64, 65] for long read or genome align- promoter. The resulting vector was transformed into ments. CNVs were identified using CNVnator-0.3.3 [66]. Solanum lycopersicum cv. Micro-Tom by A. tumefaciens And observation by IGV software and PCR amplification C58. The T ­ 0 plants were generated and inoculated with (Additional file 2: Fig. S4; Additional file 1: Table S15 and M. incognita to observe resistance and measure gene 16) was used to confirm the accuracy of SVs. expression. Similarly, one candidate gene, Pmi02g3025, was cloned BSA to locate nematode resistance gene from the leaf of “2210-198” (P. mira) and ligated to the Using a previously constructed B ­ C1 population [24], we vector and transformed into A. thaliana “Columbia.” extracted genomic DNA from fresh leaves of its parents When the transformed plants were grown to about 5 and two bulked populations which mixed by 20 resist- leaves, low temperature (0 °C) treatment was applied and ance seedlings and 20 susceptible ones. The samples samples were collected at 0, 24, and 48 h post treatment. were used for resequencing on an Illumina HiSeq 2500 Phenotype was observed and the drooping leaves were platform (Illumina, San Diego, CA) to generate 125 bp used to indicate that the accession was susceptible to low paired-end reads. Then BWA software [67] was used to temperature. The proline content of leaves in Arabidopsis align the clean reads of each sample against the reference was measured according to previous reference as previ- genome of Prunus persica v2.0.a1 [14]. ously described [71]. Variants calling were performed for all samples by using GATK3.3 software [68] to identify a total of Comparative analysis 2,442,036 SNPs on all 8 chromosomes. The homozygous Protein sequences from 15 plant species including P. SNPs between two parents were extracted from the vcf persica (phytozomev10), P. ferganensis, P. kansuensis, P. files. The reads depth information for homozygous SNPs davidiana, P. mira, Prunus dulcis (https://​www.​rosac​eae.​ above in the offspring pools was gained to calculate the org/​speci​e s/​prunus/​prunus_​dulsis/​l aura​nne/​genome_​ SNP index [69]. The loci whose SNP index was less than v1.0), P. armeniaca (https://​www.​rosac​eae.​org/​speci​es/​ 0.3 or greater than 0.7 in both the pools were filtered. The prunus_​armen​iaca/​genome_​v1.0), Prunus mume (http://​ difference of the SNP index of the two pools was calcu- prunu​smume​genome.​bjfu.​edu.​cn/​index.​jsp), P. sali- lated as the Δ(SNP index). Finally, we made a graph of cina (NCBI accession: PRJNA574159), P. avium (NCBI Δ(SNP index) of the whole-genome by the sliding win- accession: PRJNA395588), M. domestica (NCBI acces- dow method. sion: PRJNA534520), P. bretschneideri (NCBI accession: PRJNA259338), Fragaria vesca (phytozome v10), A. Promoter activity measurement thaliana (phytozome v10), and Vitis vinifera (phytozome A total of five primers upstream and one downstream v10) were used to construct orthologous gene families. of the start codon of Prupe.2G053600 were synthesized To remove redundancy caused by alternative splicing, and used to amplify a series of 5 indel regions in the we retained only the gene model at each gene locus that Prupe.2G053600 promoter using PCR amplification. The encoded the longest protein. To exclude putative frag- amplified PCR products were ligated into pGEM-T easy mented genes, genes encoding protein sequences shorter vector and cloned into pBI101 binary vector after being than 50 amino acids were filtered out. All-against-all digested by XbaI and BamHI. Furthermore, each of the 5 BLASTp was performed for these protein sequences with Cao et al. BMC Biology (2022) 20:139 Page 14 of 17 an E-value cut-off of 1­ e−5. OrthoMCL V1.4 [72] was then Identification of an interspecific hybridization of P. used to cluster genes into gene families with the param- davidiana eter “-inflation 1.5.” Using SNPs identified from 126 peach accessions, we Protein sequences from 830 single-copy gene families performed PCA to evaluate the evolution path using the were used for phylogenetic tree construction. MUSCLE software GCTA [76]. [51] was used for multiple sequence alignment for pro- To trace the origin of genes in P. davidiana, each tein sequences in each single-copy family with default gene was aligned to other genomes to calculate the parameters. The alignments from all single-copy families alignment score. Genes were classified as putatively were then concatenated into a super alignment matrix, originating from the specie which had the high- which was used for phylogenetic tree construction using est alignment scores. Collinearity among P. dulcis, P. the maximum likelihood (ML) method implemented davidiana, and P. mira was plotted based on the iden- in the RAxML software (http://​cme.h-​its.​org/​exeli​xis/​ tified syntenic blocks. web/​softw​are/​raxml/​index.​html). Divergence times between the 15 species were estimated using MCMC- Resistance genes Tree in PAML software (http://​abacus.​gene.​ucl.​ac.​uk/​ Hidden Markov model search (HMMER; http://​ softw​are/​paml.​html) with the options “correlated rates” hmmer.​janel​ia.​org) was used to identify R genes in the and “JC69” model. A Markov Chain Monte Carlo analy- four wild peach genomes according to the nucleotide- sis was run for 10,000 generations, using a burn-in of binding site (NBS) domain (PF00931), Toll and inter- 10,000 iterations and sample-frequency of 2. Five cali- leukin-1 receptor-like (TIR) model (PF01582), and bration points were applied according to the TimeTree several leucine-rich repeat (LRR) models (PF00560, database (http://​www.​timet​ree.​org): A. thaliana and V. PF07723, PF07725, PF12799, PF13306, PF13516, vinifera (106.9-134.7 Mya), A. thaliana and the common PF13504, PF13855, and PF08263) in the Pfam database ancestor of F. vesca, M. domestica, and other Prunus spe- (http://​pfam.​sanger.​ac.​uk). Coiled-coil (CC) motifs cies (97.4-109.3 Mya), F. vesca and common ancestor of were detected using the COILS prediction program M. domestica, and other Prunus species (64.5-75.5 Mya), 2.2 (https://​embnet.​vital-​it.​ch/​softw​are/​COILS_​form.​ M. domestica, P. bretschneideri, and other Prunus spe- html) with a p score cut-off of 0.9. cies (51.7-61.8 Mya), M. domestica and P. bretschneideri (11.3-26.4 Mya). To detect the whole genome duplication events, we first Identification of selective sweeps associated identified collinearity blocks using paralogous gene pairs with high‑altitude adaptation in P. mira with software MCScanX [62]. Using the sum of trans- Raw genome reads of the 49 accessions of P. mira from version of fourfold degenerate site divided by the sum of Tibet, China with different altitudes were processed fourfold degenerate site, we then calculated 4dTv values to remove adaptor, contaminated and low quality of each block. In addition, Ks values of homologous gene sequences, and the cleaned reads were mapped to the pairs were also calculated using PAML [73] based on the assembled P. mira genome using BWA v0.7.8 [67]. Based sequence alignments by MUSCLE [51], to validate spe- on the alignments, the potential PCR duplicates were ciation times. removed using the SAMtools command “rmdup.” SNP calling at the population level was performed using SAM- tools [77]. The identified SNPs supported by at least five Gene family expansions and contractions mapped reads, mapping quality ≥ 20, and Phred-scaled Expansion and contractions of orthologous gene fami- genotype quality ≥ 5, and with less than 0.2 missing data lies were determined using CAFE [74], which uses a were considered high-quality SNPs (1,394,483) and used birth and death process to model gene gain and loss over for subsequent analyses. a phylogeny. Significance of changes in gene family size A phylogenetic tree constructed using TreeBeST 1.9.2 in a phylogeny was tested by calculating the p-value on [78] with 1000 bootstraps and structure analysis was per- each branch using the Viterbi method with a randomly formed for the 49 accessions of P. mira using the program generated likelihood distribution. This method calcu- frappe [79] with the number of assumed genetic clusters lates exact p-values for transitions between the parent (K) ranging from two to three, and 10,000 iterations for and child family sizes for all branches of the phylogenetic each run. tree. Enrichment of GO terms and KEGG pathways in To identify genome-wide selective sweeps associated the expanded gene families of each of the four wild peach with high-altitude adaptation, we scanned the genome in species were identified using the R package clusterPro- 50-kb sliding windows with a step size of 10 kb and cal- filer [75]. culated the reduction inπ based on the P. mira accessions Cao et al. BMC Biology (2022) 20:139 Page 15 of 17 originating in high-altitude to low-altitude regions (πhigh/ davidiana, P. kansuensis, and P. ferganensis. Fig. S3. Genome variations πlow). In addition, Tajima’s D and FST between the two across the pseudo-chromosomes of four wild peach species compared to groups were also calculated. Windows with the top 5% of the reference (Prunus persica). The circles from the outer to the inner (A-P) represent copy number variation (CNV) density in P. ferganensis (A), P. kan- the π ratios, Tajima’s D ratio, and FST values were consid- suensis (B), P. davidiana (C), and P. mira (D), and structure variations (SVs) in ered as selective sweeps. P. ferganensis (E), P. kansuensis (F), P. davidiana (G), and P. mira (H), and indels in P. ferganensis (I), P. kansuensis (J), P. davidiana (K), and P. mira (L), as well as SNPs in P. ferganensis (M), P. kansuensis (N), P. davidiana (O), and P. mira Abbreviations (P) in each sliding window of 0.1 Mb. Fig. S4. The electrophoresis results QTP: Qinghai-Tibet Plateau; Mya: Million years ago; CNVs: Copy number varia- of 20 primers which amplified in two varieties (2010-138 and Shen Zhou tions; BSA: Bulked segregant analysis; UV: Ultraviolet rays; BUSCO: Benchmark- Li He Shui Mi) belonging to P. mira (Pm) and P. persica (Pp), respectively. ing Universal Single-Copy Orthologs; LTR: Long terminal repeat; TEs: Transpos- Fig. S5. KEGG pathways enriched in genes comprising large-effect SNPs of able elements; SNPs: Single nucleotide polymorphisms; Indels: Small insertions between P. ferganensis and P. persica. Fig. S6. KEGG pathways enriched in and deletions; SVs: Structural variants; CNVs: Copy number variations; PAVs: genes comprising indels in four wild peach species compared to P. persica. Presence/absence variations; MCRA​: Most recent common ancestor; 4DTv: Fig. S7. KEGG pathways enriched in genes comprising structure variations Fourfold synonymous third-codon transversion; Ks: Rate of synonymous muta- in P. mira and P. kansuensiss compared to P. persica. Fig. S8. KEGG pathways tion; PCA: Principal component analysis; BC1: Backcross; GUS: β-Glucuronidase; enriched in genes comprising copy number variations in P. ferganensis π: Nucleotide diversity; FST: Fixation index; CBF: C-Repeat binding transcription (a), P. kansuensis (b), P. davidiana (c), and P. mira (d) compared to P. persica. factor; LEA: Late embryogenesis abundant protein; SMRT: Single molecule Fig. S9. Venn diagram of gene families identified from the five species real-time; Hi-C: High-throughput chromosome conformation capture; GO: of peach. Fig. S10. Statistics of single-copy orthologs, multiple-copy Gene ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; PCR: Poly- orthologs, and unique orthologs in 15 species. Fig. S11. KEGG pathways merase chain reaction; RNA-Seq: RNA sequencing; CaMV: Cauliflower mosaic enriched in expanded and contracted gene families of the four wild peach virus; ML: Maximum likelihood; NBS: Nucleotide-binding site; LRR: Leucine-rich species. Fig. S12. Whole-genome duplication and speciation events in repeat; TIR: Toll and interleukin-1 receptor-like; CC: Coiled-coil; SRA: Sequence peach as revealed by the distribution of 4DTv distance among paralogous read archive. and orthologs genes in different species. Fig. S13. Distribution of resist- ance (R) genes across the 8 chromosomes in five peach species of peach and their overlaps with disease resistance QTLs. Fig. S14. Percent of P. Supplementary Information davidiana-specific contigs covered by reads from different Prunus species. The online version contains supplementary material available at https://​doi.​ Fig. S15. Collinearity among P. davidiana, P. mira, P. dulcis, P. armeniaca, org/​10.​1186/​s12915-​022-​01342-y. and P. avium genomes. Fig. S16. Geographical distribution of P. mira (red circle), P. davidiana (yellow), and P. dulcis (orange) which originated in China. Fig. S17 .Expression of 19 R genes in two accessions (‘Hong Gen Additional file 1: Table S1. Genome survey of four wild peach species Gan Su Tao 1#’ and ‘Bailey’) inoculated with nematode after 0, 6, 12, 36, 60, (kmer = 17). Table S2. Summary of genome sequencing of four wild and 84 hours. Fig. S18. The evaluation of resistant to root-knot nematodes peach species. Table S3. Pseudochromosome lengths (bp) of the P. in transgenic tomato lines with Prupe.2G053600 gene and its wild type. mira, P. davidiana, P. kansuensis, and P. ferganensis assembly. Table S4. Fig. S19. The vector construction (a) and promoter activity assay (b) of BUSCO analysis of the genome assemblies of four wild peach species. Prupe.2G253600 promoter which recovered 35 bp deletion in ‘Hong Gen Table S5. Mapping statistics of RNA-Seq reads to the corresponding Gan Su Tao 1#’ peach. Fig. S20. The altitude distribution of accessions used genome assemblies of four wild peach species. Table S6. Statistic of to identify altitude adaptability related genes for genome resequencing. repeat sequences in the assemblies of four wild peach species. Table S7. Fig. S21. Phylogenetic tree of 32 accessions of P. mira (a) originating from Prediction of protein-coding genes in the genomes of four wild peach regions with different altitudes (b) and the population structure (c) when species. Table S8. Statistics of predicted protein-coding genes in four K=2 and 3. Fig. S22. Two genome regions associated with high-altitude wild peach species compared to other species. Table S9. Statistics of adaptation. Fig. S23. The gene expression of wild-type (WT) and three gene functional annotation in the four wild peach species. Table S10. transgenic plants (OE1-OE3) in Arabidopsis which induced by 4 °C for Non-coding RNAs identified in genomes of four wild wild peach species. different period. Fig. S24. The proline content of wild-type (WT) and Table S11. SNPs identified between genomes of each of the four wild transgenic plant (OE) in Arabidopsis which induced by 4 °C for 24 and 48 species and P. persica. Table S12. Small indels (<50 bp) identified between hours. Fig. S25. A deletion identified in the promoter of Pmi02g3025 gene genomes of each of the four wild species and P. persica. Table S13. in varieties with different altitude. Structural variants (≥ 50 bp) between the four wild species and P. persica. Table S14. Statistics of copy number variations between the four wild species and P. persica. Table S15. A total of 20 SVs which aligned P. mira Acknowledgements with P. persica genome. Table S16. The primers which designed to amplify We thank Dr. Xuwu Sun from the Henan University, Henan, China, for his the SVs between P. mira and P. persica genome. Table S17. Statistics of helpful suggestion. We thank the anonymous reviewers and editors for their resistance genes in the four wild peach species. Table S18. List of 175 constructive suggestions on our manuscript. peach samples used in the study. Table S19. The summary statistics of genome resequencing. Table S20. Variations in the promoter and mRNA Authors’ contributions regions of R genes on Chr. 2 (5-7 Mb) those were specific to P. kansuensis. K.C. and L.W. conceived the project. Y.L. and X.Z. contributed to the original con- Table S21. The cis-elements in promoter region of Prupe.2G053600 gene cept of the project. W.F., C.C., X.W., and J.W. collected the samples and performed in nematode-resistant and susceptible accessions. Table S22. Genes phenotyping. K.L. conducted the gene expression analysis and transgenic experi- selected between the two subgroups of P. mira which originated from ments. K.C., Z.P., and Z.F. analyzed the data. P.A. and Z.F. provided supervision. K.C. high- and low-altitude regions. Table S23. Selected genes in the KEGG wrote the paper. All authors read and approved the final manuscript. pathways associated with plateau adaptability. Table S24. The genotype of 35 bp deletion in the promoter of Prupe. 2G053600 in varieties with Funding different root-knot nematode resistance. Table S25. The genotype of 19 This work was supported by the National Key Research and Development Program bp deletion in the promoter of Pmi02g3025 in varieties with different cold (2019YFD1000203), the Agricultural Science and Technology Innovation Program resistance. Table S26. The primers used in the study. (CAAS-ASTIP-2019-ZFRI-01), and National Horticulture Germplasm Resources Center. Additional file 2: Fig. S1. Estimation of genome sizes of P. mira (a), P. davidiana (b), P. kansuensis (c), and P. ferganensis (d) based on K-mer Availability of data and materials analysis. Fig. S2. Genome size, annotated gene number, and the length Raw resequencing data for 175 peach accessions generated in this study of tandam repeat sequences, interpersed repeat sequences in P. mira, P. have been deposited into the NCBI database. And the accession numbers Cao et al. BMC Biology (2022) 20:139 Page 16 of 17 of Sequence Read Archive (SRA) were listed in Additional file 1: Table S19. In 10. Guan JT, Xu YG, Yu Y, Fu J, Ren F, Guo JY, et al. Genome structure vari- addition, the RNA-seq data of different tissues of four wild species and P. mira ation analyses of peach reveal population dynamics and a 1.67 Mb which were treated by UV or low temperature also have been deposited into causal inversion for fruit shape. Genome Biol. 2021;22:13. https://​doi.​ the NCBI database with BioProject under accession PRJNA753201 [80] and org/​10.​1186/​s13059-​020-​02239-1. PRJNA753549 [81]. The assemblies of four genomes have been uploaded to 11. Cao K, Yang XW, Li Y, Zhu GR, Fang WC, Chen CW, et al. New high-qual- Genome Database for Rosaceae (https://​www.​rosac​eae.​org/​Analy​sis/​12080​ ity peach (Prunus persica L. Batsch) genome assembly to analyze the 707; https://​www.​rosac​eae.​org/​Analy​sis/​11857​924; https://​www.​rosac​eae.​org/​ molecular evolutionary mechanism of volatile compounds in peach Analy​sis/​12080​706; https://​www.​rosac​eae.​org/​Analy​sis/​12080​705). The prim- fruits. The. Plant J. 2021;108:281–95. https://​doi.​org/​10.​1111/​tpj.​15439. ers used in the study can be found in Additional file 1: Table S26. 12. Tan QT, Li S, Zhang YZ, Chen M, Wen BB, Jiang S, et al. Chromosome- level genome assemblies of five Prunus species and genome-wide association studies for key agronomic traits in peach. Hortic Res. Declarations 2021;8:213. 13. Wang X, Liu SJ, Zuo H, Zheng WK, Zhang SS, Huang Y, et al. Genomic Ethics approval and consent to participate basis of high-altitude adaptation in Tibetan Prunus fruit trees. Curr Biol. Not applicable. 2021;31:3848–60. 14. Verde I, Jenkins J, Dondini L, Micali S, Pagliarani G, Vendramin E, et al. Consent for publication The Peach v2.0 release: high-resolution linkage mapping and deep Not applicable. resequencing improve chromosome-scale assembly and contiguity. BMC Genomics. 2017;18:225. Competing interests 15. Simao FA, Waterhouse RM, Ioannidis P, Kriventseva EV, Zdobnov EM. The authors declare that they have no competing interests. BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics. 2015;31:3210–2. Author details 16. Baek S, Choi K, Kim GB, Yu HJ, Cho A, Jang H, et al. Draft genome 1 The Key Laboratory of Biology and Genetic Improvement of Horticultural sequence of wild Prunus yedoensis reveals massive inter-specific Crops (Fruit Tree Breeding Technology), Ministry of Agriculture, Zhengzhou hybridization between sympatric flowering cherries. Genome Biol. Fruit Research Institute, Chinese Academy of Agricultural Sciences, Zheng- 2018;19:127. zhou 450009, People’s Republic of China. 2 Novogene Bioinformatics Institute, 17. Alioto T, Alexiou KG, Bardil A, Barteri F, Castanera R, Cruz F, et al. Trans- Beijing, People’s Republic of China. 3 IRTA, Centre de Recerca en Agrigenòmica, posons played a major role in the diversification between the closely CSIC-IRTA-UAB-UB, Campus UAB – Edifici CRAG, Cerdanyola del Vallès (Bel- related almond and peach genomes: results from the almond genome laterra), Barcelona, Spain. 4 Boyce Thompson Institute for Plant Research, sequence. Plant J. 2020;101:455–72. Cornell University, Ithaca, NY 14853, USA. 5 National Horticulture Germplasm 18. Chung SL, Lo CH, Lee TY, Zhang YQ, Xie YW, Li XH, et al. Diachronous uplift Resources Center, Zhengzhou Fruit Research Institute, Chinese Academy of the Tibetan plateau starting 40?Myr ago. Nature. 1998;394:769–73. of Agricultural Sciences, Zhengzhou 450009, People’s Republic of China. 19. Lespinet O, Wolf YI, Koonin EV, Aravind L. The role of lineage-specific gene family expansion in the evolution of Eukaryotes. Genome Res. Received: 19 December 2021 Accepted: 30 May 2022 2002;12:1048–58. 20. Rada F, García-Núñez C, Boero C, Gallardo M, Hilal M, González J, et al. Low-temperature resistance in Polylepis tarapacana, a tree growing at the highest altitudes in the world. Plant Cell Environ. 2001;24:377–81. 21. Yang F, Xu X, Xiao X, Li C. Responses to drought stress in two poplar spe- References cies originating from different altitudes. Biol Plantarum. 2009;53:511. 1. Abbott A, Georgi L, Yvergniaux D, Wang Y, Blenda A, Reighard G, et al. 22. Wang L, Zhu G, Fang W, Cao K, Wang X, Chen C, et al. Genetic diversity Peach: The model genome for Rosaceae. Acta Hortic. 2002;575:145–55. of peach wild relatives. In: Peach genetic resource in China (ed. Wang L, https://​doi.​org/​10.​17660/​ActaH​ortic.​2002.​575.​14. et al.). Beijing: China Agriculture Press; 2012. p. 16–98. 2. International Peach Genome Initiative. The high-quality draft genome 23. Zhu GR, Wang LR, Zuo QY, Zhang XW. Studies on the resistance of peach of peach (Prunus persica) identifies unique patterns of genetic diversity, rootstocks to root-knot nematodes. J Fruit Sci. 2000;S1:36–9. domestication and genome evolution. Nat Genet. 2013;45:487–94. 24. Cao K, Wang LR, Zhao P, Zhu GR, Fang WC, Chen CW, et al. Identification https://​doi.​org/​10.​1038/​ng.​2586. of a candidate gene for resistance to root-knot nematode in a wild peach 3. Cao K, Zheng ZJ, Wang LR, Liu X, Zhu GR, Fang WC, et al. Comparative and screening of its polymorphisms. Plant Breed. 2014;33:530–5. population genomics reveals the domestication history of the peach, 25. Zhai XG, Zhao T, Liu YH, Long H, Deng GB, Pan ZF, et al. Characterization Prunus persica, and human influences on perennial fruit crops. Genome and expression profiling of a novel cereal cyst nematode resistance gene Biol. 2014;15:415. https://​doi.​org/​10.​1186/​s13059-​014-​0415-1. analog in wheat. Mol Biol. 2008;42:960–5. 4. Yu Y, Fu J, Xu YG, Zhang JW, Ren F, Zhao HW, et al. Genome re-sequencing 26. Zhu X, Xiao K, Cui HY, Hu JF. Overexpression of the Prunus sogdiana NBS- reveals the evolutionary history of peach fruit edibility. Nat Commun. LRR subgroup gene PsoRPM2 promotes resistance to the root-knot nema- 2018;9:5404. https://​doi.​org/​10.​1038/​s41467-​018-​07744-3. tode Meloidogyne incognita in tobacco. Front Microbiol. 2017;8:2113. 5. Akagi T, Hanada T, Yaegaki H, Gradziel TM, Tao R. Genome-wide view of 27. Cao K, Li HY, Wang Q, Zhao P, Zhu GR, Fang WC, et al. Comparative tran- genetic diversity reveals paths of selection and cultivar differentiation in scriptome analysis of genes involved in the response of resistant and sus- peach domestication. DNA Res. 2016;23:271–82. https://​doi.​org/​10.​1093/​ ceptible peach cultivars to nematode infection. Sci Hortic. 2017;215:20–7. dnares/​dsw014. 28. Yang J, Jin ZB, Chen J, Huang XF, Li XM, Liang YB, et al. Genetic signa- 6. Li Y, Cao K, Zhu GR, Fang WC, Chen CW, Wang XW, et al. Genomic analyses tures of high-altitude adaptation in Tibetans. Proc Natl Acad Sci U S A. of an extensive collection of wild and cultivated accessions provide new 2017;114:4189–94. insights into peach breeding history. Genome Biol. 2019;20:36. https://​ 29. Li MZ, Tian SL, Jin L, Zhou GY, Li Y, Zhang Y, et al. Genomic analyses iden- doi.​org/​10.​1186/​s13059-​019-​1648-9. tify distinct patterns of selection in domesticated pigs and Tibetan wild 7. Cao K, Zhou ZK, Wang Q, Guo J, Zhao P, Zhu GR, et al. Genome-wide boars. Nat Genet. 2013;45:1431–8. association study of 12 agronomic traits in peach. Nat Commun. 30. Qu YH, Zhao HW, Han NJ, Zhou GY, Song G, Gao B, et al. Ground tit 2016;7:13246. https://​doi.​org/​10.​1038/​ncomm​s13246. genome reveals avian adaptation to living at high altitudes in the Tibetan 8. Cao K, Li Y, Deng CH, Gardiner SE, Zhu GR, Fang WC, et al. Comparative plateau. Nat Commun. 2013;4:2071. population genomics identified genomic regions and candidate genes 31. Huerta-Sánchez E, Jin X, Asan BZ, Peter BM, Vinckenbosch N, Liang Y, et al. associated with fruit domestication traits in peach. Plant Biotechnol J. Altitude adaptation in Tibetans caused by introgression of Denisovan-like 2019;17:1954–70. https://​doi.​org/​10.​1111/​pbi.​13112. DNA. Nature. 2014;512:194–7. 9. Layne DR, Bassi D. The peach: botany, production and uses. Cambridge: CAB; 2008. Cao et al. BMC Biology (2022) 20:139 Page 17 of 17 32. Li JT, Gao YD, Xie L, Deng C, Shi P, Guan ML, et al. Comparative genomic 58. Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, Van Baren MJ, et al. investigation of high-elevation adaptation in ectothermic snakes. Proc Transcript assembly and quantification by RNA-Seq reveals unannotated Natl Acad Sci U S A. 2018;115:8406–11. transcripts and isoform switching during cell differentiation. Nat Biotech. 33. Zeng XQ, Long H, Wang Z, Zhao SC, Tang YW, Huang ZY, et al. The draft 2010;28:511–5. genome of Tibetan hulless barley reveals adaptive patterns to the high 59. Grabherr M, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Full- stressful Tibetan Plateau. Proc Natl Acad Sci U S A. 2015;112:1095–100. length transcriptome assembly from RNA-Seq data without a reference 34. Zhang TC, Qiao Q, Novikova PY, Wang Q, Yue JP, Guan YL, et al. Genome of genome. Nat Biotech. 2011;29:644–52. Crucihimalaya himalaica, a close relative of Arabidopsis, shows ecological 60. Haas BJ, Salzberg SL, Zhu W, Pertea M, Allen JE, Orvis J, et al. Automated adaptation to high altitude. Proc Natl Acad Sci U S A. 2019;116:7137–46. eukaryotic gene structure annotation using EVidenceModeler and the 35. Xie M, Chung CYL, Li MW, Wong FL, Wang X, Liu AL, et al. A reference- Program to Assemble Spliced Alignments. Genome Biol. 2008;9:1–22. grade wild soybean genome. Nat Commun. 2019;10:1216. 61. Lowe TM, Eddy SR. tRNAscan-SE: a program for improved detec- 36. Tambasco-Studart M, Tews I, Amrhein N, Fitzpatrick TB. Functional tion of transfer RNA genes in genomic sequence. Nucleic Acids Res. analysis of PDX2 from Arabidopsis, a glutaminase involved in vitamin B6 1997;25:955–64. biosynthesis. Plant Physiol. 2007;144:915–25. 62. Wang Y, Tang HB, Debarry JD, Tan XF, Li JP, Wang XY, et al. MCScanX: a 37. Colinas M, Eisenhut M, Tohge T, Pesquera M, Fernie AR, Weber APM, et al. toolkit for detection and evolutionary analysis of gene synteny and collin- Balancing of B6 vitamers is essential for plant development and metabo- earity. Nucleic Acids Res. 2012;40:e49. lism in Arabidopsis. Plant Cell. 2016;28:439–53. 63. Harris RS. Improved pairwise alignment of genomic DNA. PhD thesis, The 38. Liu R, Zhang QN, Lu J, Zhang CH, Zhang L, Wu Y. The effects of exogenous Pennsylvania State University. 2007. pyridoxal-5-phosphate on seedling growth and development of wheat 64. Mitsuhashi S, Ohori S, Katoh K, Frith MC, Matsumoto N. A pipeline for under salt stress. Cereal Res Commun. 2019;47:442–54. complete characterization of complex germline rearrangements from 39. Hussain S, Huang J, Zhu CQ, Zhu LF, Cao XC, Hussain S, et al. Pyridoxal long DNA reads. Genome Med. 2020;12:67. 5’-phosphate enhances the growth and morpho-physiological charac- 65. Gong L, Wong CH, Cheng WC, Tjong H, Menghi F, Ngan CY, et al. Picky teristics of rice cultivars by mitigating the ethylene accumulation under comprehensively detects high-resolution structural variants in nanopore salinity stress. Plant Physiol Biochem. 2020;154:782–95. long reads. Nat Methods. 2018;15:455–60. 40. Huang SH, Zhang JY, Wang LH, Huang LQ. Effect of abiotic stress on the 66. Abyzov A, Urban AE, Snyder M, Gerstein M. Cnvnator: an approach to abundance of different vitamin B6 vitamers in tobacco plants. Plant discover, genotype and characterize typical and atypical cnvs from family Physiol Biochem. 2013;66:63–7. and population genome sequencing. Genome Res. 2011;21:974–84. 41. Yuan DY, Chen XH, Gu HR, Zou M, Zou Y, Fang J, et al. Chromosomal 67. Li H, Durbin R. Fast and accurate short read alignment with Burrows– genome of Triplophysa bleekeri provides insights into its evolution and Wheeler transform. Bioinformatics. 2009;25:1754–60. environmental adaptation. GigaScience. 2020;9:1–14. 68. McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, 42. Liu F, Liang T, Zhang ZY, Liu LJ, Li J, Dong WX, et al. Effects of altitude on et al. The Genome Analysis Toolkit: a MapReduce framework for analyzing human oral microbes. AMB Expr. 2021;11:41. next-generation DNA sequencing data. Genome Res. 2010;20:1297–303. 43. Li H, Wang YJ, Yu QL, Feng TS, Zhou R, Shao LY, et al. Elevation is associ- 69. Takagi H, Abe A, Yoshida K, Kosugi S, Natsume S, Mitsuoka C, et al. ated with human skin microbiomes. Microorganisms. 2019;7:611. QTL-seq: rapid mapping of quantitative trait loci in rice by whole 44. Hirsch CN, Foerster JM, Johnson JM, Sekhon RS, Muttoni G, Vaillancourt B, genome resequencing of DNA from two bulked populations. Plant J. et al. Insights into the maize pan-genome and pan-transcriptome. Plant 2013;74:174–83. Cell. 2014;26:121–35. 70. Jefferson RA, Kavanagh TA, Bevan MW. Gus fusion: β-glucuronidase as 45. Hamilton JA, Miller JM. Adaptive introgression as a resource for manage- a sensitive and versatile gene fusion maiker in higher plants. EMBO J. ment and genetic conservation in a changing climate. Conserv Biol. 1987;6:3901–7. 2016;30:33–41. 71. Zhao LN, Liu FX, Xu WQ, Di C, Zhou SX, Xue YB, et al. Increased expression 46. McDonald MJ, Rice DP, Desai MM. Sex speeds adaptation by altering the of OsSPX1 enhances cold/subfreezing tolerance in tobacco and Arabidop- dynamics of molecular evolution. Nature. 2016;531:233–6. sis thaliana. Plant Biotechnol J. 2009;7:550–61. 47. Lieberman-Aiden E, Van Berkum NL, Williams L, Imakaev M, Ragoczy 72. Li L, Stoeckert CJ, Roos DS. OrthoMCL: Identification of ortholog groups T, Telling A, et al. Comprehensive mapping of long-range interac- for eukaryotic genomes. Genome Res. 2003;13:2178–89. tions reveals folding principles of the human genome. Science. 73. Yang Z. PAML 4: phylogenetic analysis by maximum likelihood. Mol 2009;326:289–93. BiolEvol. 2007;24:1586–91. 48. Marcais G, Kingsford C. A fast, lock-free approach for efficient parallel 74. De Bie T, Cristianini N, Demuth JP, Hahn MW. CAFE: a computational tool counting of occurrences of k-mers. Bioinformatics. 2011;27:764–70. for the study of gene family evolution. Bioinformatics. 2006;22:1269–71. 49. Cheng HY, Concepcion GT, Feng XW, Zhang HW, Li H. Haplotype-resolved 75. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing de novo assembly using phased assembly graphs with hifiasm. Nat Meth- biological themes among gene clusters. OMICS. 2012;16:284–7. ods. 2021;18:170–5. 76. Yang JA, Lee SH, Goddard ME, Visscher PM. GCTA: a tool for genome-wide 50. Zhang X, Zhang S, Zhao Q, Ming R, Tang H. Assembly of allele-aware, complex trait analysis. Am J Hum Genet. 2011;88:76–82. chromosomal-scale autopolyploid genomes based on Hi-C data. Nat 77. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The Plants. 2019;5:833–45. sequence alignment/map format and SAMtools. Bioinformatics. 51. Edgar RC. MUSCLE: multiple sequence alignment with high accuracy and 2009;25:2078–9. high throughput. Nucleic Acids Res. 2004;32:1792–7. 78. Vilella AJ, Severin J, Uretavidal A, Heng L, Durbin R, Birney E. Ensem- 52. Xu Z, Wang H. LTR_FINDER: An efficient tool for the prediction of full- blCompara GeneTrees: complete, duplication-aware phylogenetic trees length LTR retrotransposons. Nucleic Acids Res. 2007;35:W265–8. in vertebrates. Genome Res. 2009;19:327–35. 53. Jurka J, Kapitonov VV, Pavlicek A, Klonowski P, Kohany O, Walichiewicz J. 79. Tang H, Peng J, Wang P, Risch NJ. Estimation of individual admix- Repbase Update, a database of eukaryotic repetitive elements. Cytogenet ture: analytical and study design considerations. Genet Epidemiol. Genome Res. 2005;110:462–7. 2005;28:289–301. 54. Stanke M, Steinkamp R, Waack S, Morgenstern B. AUGUSTUS: a web server 80. Different tissue of four wild species of peach. NCBI accession: PRJNA753201. for gene finding in eukaryotes. Nucleic Acids Res. 2004;32:W309–12. (2021). https://​www.​ncbi.​nlm.​nih.​gov/​biopr​oject/​PRJNA​753201. 55. Majoros WH, Pertea M, Salzberg SL. TigrScan and GlimmerHMM: 81. The cold and UV treatment of P. mira. NCBI accession: PRJNA753201. two open source ab initio eukaryotic gene-finders. Bioinformatics. (2021). https://​www.​ncbi.​nlm.​nih.​gov/​biopr​oject/​PRJNA​753549. 2004;20:2878–9. 56. Korf I. Gene finding in novel genomes. BMC Bioinformatics. 2004;5:59. 57. Kim D, Paggi JM, Park C, Bennett C, Salzerg SL. Graph-based genome Publisher’s Note alignment and genotyping with HISAT2 and HISAT-genotype. Nat Bio- Springer Nature remains neutral with regard to jurisdictional claims in pub- tech. 2019;37:907–15. lished maps and institutional affiliations.