-
上饶早梨Pyrus pyrifolia是江西省上饶县花厅、田墩、五府山等11个乡镇的国家地理标志农产品，被认为是药食兼优的夏令佳品[1-2]，栽培历史长，以‘黄皮消’‘Huangpixiao’和‘六月雪’‘Liuyuexue’2个品种品质最为优良[3]。对上饶早梨脱毒快繁[4-5]和品种鉴定[6-7]的研究表明：来源上饶县花厅镇的上饶早梨道地性农家种‘六月雪’和‘黄皮消’经简单重复序列标记（SSR）和扩增片段长度多态性（AFLP）分子标记的聚类分析聚为1类，但生物学性状、农艺性状和果实品质均有一定差异[8]。因此，找出两者基因差异对地方品种的育种选种工作十分重要。自2000年模式植物拟南芥Arabidopsis thaliana的基因组被完全解析后，已有越来越多的种质被全基因组测序[9]；研究认为，野生种、农家种、栽培种的全基因组重测序是筛选重要性状关键基因的一个重要内容[10]。梨基因组测序的完成[11-13]和高通量测速技术的快速发展，使梨种质资源的全基因组变异分析成为可能，而全基因组单核苷酸多肽位点（SNP）和插入缺失位点（INDEL）相关标记的开发，对作物分子辅助育种和基因组学研究具有重要的作用。周贺等[14]对褐色砂梨‘黄花’Pyrus pyrifolia ‘Huanghua’色泽形成期的果皮转录组数据进行SNP分子标记开发，共筛选到1 178个可能与果皮色泽形成相关的SNP标记。李节法[15]以6年生砂梨品种‘翠冠’P. pyrifolia ‘Cuiguan’的果实膨大早期、中期和成熟期样品进行比较转录组学分析，鉴定到4 121个选择性剪切位点，30 560个SNP位点和7 443个SSR标记位点。WU等[16]、MONTANARI等[17]和TERAKAMI等[18]也通过SNP标记构建了梨遗传图谱。但梨的INDEL标记研究较少。本研究以上饶县花厅镇上饶早梨道地性农家种‘六月雪’和‘黄皮消’为材料，通过全基因组重测序，深度挖掘其基因组SNP，INDEL和结构变异位点（SV）等位点，为上饶早梨优质品种相关标记的开发、优异基因的挖掘提供参考。
HTML
-
上饶县花厅镇上饶早梨道地性农家品种‘六月雪’和‘黄皮消’试管苗由上饶师范学院生命科学学院植物组织培养室提供。
-
参照十六烷基三甲基溴化铵（CTAB）提取法提取样品DNA。
-
以吸光度比值D（260）/D（280）和Qubit 2.0（Life technologies, USA），Bioanalyzer 2100（Agilent，Germany）软件分析完成总DNA样品的质量控制。称取1 μg基因组总DNA片段化处理至300~400 bp，进行末端修复、末端加‘A’处理、接头（Adapters）连接反应，连接至Illumina公司测序平台的测序接头后进行聚合酶链式反应（PCR）扩增。连接产物用质量分数2%的琼脂糖凝胶电泳，选择400~500 bp的片段，随后用连接介导PCR（LM-PCR）进行扩增获得DNA文库。按照Illumina公司HiSeq 2500测序系统（Hiseq 2500）的操作说明对形成的DNA文库进行双端125 bp的高通量测序。
-
鉴于Illlumina Hiseq 2500错误率对结果的影响，对原始数据进行质量控制（QC）。参考数据库为西洋梨基因组Pyrus communis Genome v1.0数据库。分别对每个样本使用bowtie 2软件进行测序短序列匹配（reads mapping），并用UnifiedGenotyper软件进行SNP和INDEL的提取[19-21]。采用ANNOVAR软件对获得的SNP和INDEL进行功能注释[22]。SNP常见变异分析（common variation, CV）及差异INDEL分析首先获取2个样品相同位置的SNP/INDEL，再根据非同义SNP（nonsynonymous SNP，nsSNP）/INDEL获取相关基因[23]。将差异SNP和INDEL分别与转录组数据进行关联分析，分别考察差异SNP和INDEL相关的表达数据[24]。候选基因的富集分析递交至AgriGO软件用于富集基因本体术语（gene ontology terms）[25]。拷贝数变异位点（CNV）分析采用VarScan软件进行[26]。
1.1. 材料
1.2. 方法
1.2.1. DNA提取
1.2.2. 测序
1.2.3. 数据分析
-
以HiSeq 2500测序系统提供的起初测序数据为原始数据，即各样本测序得到的短序列（reads）数及碱基总数，共得到275 092 822个短序列和34 661 695 572个碱基，其中‘六月雪’中含短序列140 696 312个，碱基17 727 735 312个，‘黄皮消’中短序列134 396 510个，碱基16 933 960 260个。为剔除Illlumina Hiseq 2500错误率对结果的影响，需对原始数据进行质量控制，包括去除低质量序列，去除接头，以进行后续工作。质量控制后得到232 434 654个短序列和29 286 766 152个碱基，总有效数据比例为84.5%；其中‘六月雪’含短序列119 056 332个，碱基150 010 978 322个，有效数据比例为84.6%；‘黄皮消’含短序列113 378 320个，碱基14 285 668 320个，有效数据比例为84.4%。
-
经过测序将短序列匹配至参考基因组。‘六月雪’组中总匹配的短序列数为68 859 074个，占所有短序列数的57.8%，含唯一匹配的短序列数为31 889 494个，占总匹配数的26.8%；覆盖全基因的深度为24.35，覆盖全基因组的百分比为93.0%；当覆盖深度≥3时，覆盖全基因组的百分比为89.4%。‘黄皮消’组总匹配的短序列数为66 580 757个，占所有短序列数的58.7%，含唯一匹配的短序列数为32 165 247个，占总匹配数的28.4%；覆盖全基因的深度为23.89，覆盖全基因组的百分比为93.0%，覆盖深度≥3时，覆盖全基因组的百分比为89.5%。
-
‘六月雪’中共得到SNP 6 171 357个，在编码区的无义突变有335 659个，有义突变有281 871个；因SNP突变获得终止子6 001个，丢失终止子1 226个；在基因5'非翻译区（5'UTR内的SNP总数、在3'UTR内的SNP总数及位于5'UTR和3'UTR间的SNP总数均为0；在不同可变剪切的基因组区域内找到SNP 3 383个，内含子区域内找到1 298 966个，启动子区域内找到1 301 726个，基因间区域内找到2 942 525个。
‘黄皮消’中共得到SNP 6 140 603个，在编码区的无义突变有332 280个，有义突变有278 064个；因SNP突变获得终止子6 034个，丢失终止子1 210个；在基因5'非翻译区（5'UTR）内的SNP总数、在3'UTR内的SNP总数及位于5'UTR和3'UTR间的SNP总数均为0，在不同可变剪切的基因组区域内找到SNP 3 274个，内含子区域找到1 285 052个，启动子区域内找到1 291 598个，基因间区域内找到2 943 091个。
-
对获得的335 659个（‘六月雪’）和332 280个（‘黄皮消’）nsSNP进行常见变异分析发现，2个品种共有2 282个nsSNP关联了2 067个基因，nsSNP对基因的平均关联频率超过了1。其中，烟草花叶病毒耐药蛋白N（PCP017781），假定的抗病RPP13样蛋白1（PCP007357），可能的抗病RPP8样蛋白2（PCP030706），可能的LRR类受体丝氨酸/苏氨酸蛋白激酶At3g47570（PCP021305），未注释蛋白1（PCP008176），烟草花叶病毒耐药蛋白N（PCP007457），烟草花叶病毒耐药蛋白N（PCP018815），未注释蛋白2（PCP021753），含重复锚蛋白的蛋白质At5g02620（PCP022078），烟草花叶病毒耐药蛋白N（PCP030478），抗病蛋白RGA2（PCP014224），可能的LRR类受体丝氨酸/苏氨酸蛋白激酶At1g53420（PCP031574），假定的抗病蛋白RGA3（PCP023580）等蛋白编码的基因则关联了10个以上nsSNP。为进一步研究候选基因的选择压力，对得到的2 067个基因的nsSNP与同义SNP的比值（nsSNP/synonymous SNP，r）进行考察，发现r的对数呈现正态分布（图 1），其值约为2，说明进化有一定的正向选择压力。
Figure 1. Curve map of nsSNP/synonymous SNP ratio distribution of candidate genes
-
由于生物学定义混乱的原因，不同的生物学数据库可能会使用不同的术语。为了解决这个问题，基因本体联合会（Gene Onotology Consortium）建立了“基因本体论”（gene ontolog，GO）数据库，目的是通过利用统一化的、结构化的语言建立一个适用于不同物种的、对基因和蛋白质功能进行定义和描述，并且能够随着研究的不断深入而更新的语言词汇标准。GO数据库包含基因参与的生物过程、所处细胞位置及具有的分子功能3个方面信息，其注释信息可对基因功能进行预测。GO显著性富集分析以基因本体术语（GO term）为单位，确定差异表达基因行使的主要生物学功能。分析全局GO功能与候选基因所处的功能发现，刺激、结合反应等功能的基因相对于背景基因（1 306条）富集（图 2）。通过差异基因富集得到的GO terms（图 3及表 1）可知，全部25 698条信息中，ADP结合、蛋白酪氨酸激酶活性、腺苷核糖核酸结合、嘌呤核苷结合、核苷结合、腺苷核苷酸结合、防御反应、蛋白丝氨酸/苏氨酸激酶活性、转移酶活性、转运含磷基团、嘌呤核糖裂解键、核糖核酸结合、嘌呤核苷酸结合、蛋白激酶活性、RNA导向的DNA聚合酶活性、磷酸转移酶活性、醇基作为受体、DNA聚合酶活性、RNA依赖性DNA复制、翻译后蛋白质修饰、脯氨酸氨基酸、磷酸化、蛋白质改性过程等GO terms具有显著性意义。
Figure 2. Overall GO distribution of candidate gene (nsSNP)
Figure 3. GO enriched histogram of candidate gene (nsSNP)
GO数据库中唯一的标记信息 GO term类型 GO功能的描述信息 输入的具有GO term注释的候选基因数 该GO term中基因总数 富集分析统计学显著水平 错误发现率 GO:0043531 分子功能F ADP结合 60 415 9.20 E-13 9.60 E-10 GO:0004713 分子功能F 蛋白酪氨酸激酶活性 145 1 649 1.10 E-09 2.50 E-07 GO:0032559 分子功能F 腺苷核糖核酸结合 297 3 947 9.30 E-10 2.50 E-07 GO:0001883 分子功能F 嘌呤核苷结合 310 4 170 1.50 E-09 2.50 E-07 GO:0001882 分子功能F 核苷结合 313 4 181 6.50 E-10 2.50 E-07 GO:0030554 分子功能F 腺苷核苷酸结合 310 4 170 1.50 E-09 2.50 E-07 GO:0006952 生物学功能P 防御响应 52 402 1.50 E-09 2.10 E-06 GO:0004674 分子功能F 蛋白丝氨酸/苏氨酸激酶活性 118 1 331 1.90 E-08 2.80 E-06 GO:0016772 分子功能F 转移含磷基团的转移酶活性 213 2 779 2.90 E-08 3.80 E-06 GO:0032555 分子功能F 嘌呤核糖核酸结合 309 4 320 4.50 E-08 4.70 E-06 GO:0032553 分子功能F 核糖核酸结合 309 4 320 4.50 E-08 4.70 E-06 GO:0017076 分子功能F 嘌呤核苷酸结合 323 4 552 5.10 E-08 4.90 E-06 GO:0004672 分子功能F 蛋白激酶活性 153 1 995 1.50 E-06 1.30 E-04 GO:0003964 分子功能F RNA导向的DNA聚合酶活性 19 106 1.70 E-06 1.40 E-04 GO:0016773 分子功能F 以醇基为受体的磷酸转移酶活性 168 2 246 2.20 E-06 1.50 E-04 GO:0034061 分子功能F DNA聚合酶活性 22 138 2.20 E-06 1.50 E-04 GO:0006278 生物学功能P 依赖RNA的DNA复制 19 106 1.70 E-06 6.40 E-04 GO:0043687 生物学功能P 蛋白质翻译后修饰 170 2 265 1.60 E-06 6.40 E-04 GO:0006468 生物学功能P 蛋白氨基酸磷酸化 153 1 979 9.80 E-07 6.40 E-04 GO:0006464 生物学功能P 蛋白质修饰过程 178 2 429 3.80 E-06 1.10 E-03 GO:0016301 分子功能F 激酶活性 164 2 277 2.00 E-05 1.30 E-03 GO:0005524 分子功能F ATP结合 240 3 538 2.40 E-05 1.50 E-03 GO:0016779 分子功能F 核苷酸转移酶活性 36 335 2.70 E-05 1.60 E-03 GO:0005515 分子功能F 蛋白质结合 395 6 182 3.30 E-05 1.80 E-03 GO:0043412 生物学功能P 高分子修饰 180 2 513 1.20 E-05 2.30 E-03 GO:0006950 生物学功能P 胁迫响应 76 882 1.10 E-05 2.30 E-03 GO:0016310 生物学功能P 磷酸化 158 2 160 1.20 E-05 2.30 E-03 GO:0030246 分子功能F 碳水化合物结合 36 353 8.00 E-05 4.20 E-03 GO:0000166 分子功能F 核苷酸结合 362 5 741 1.40 E-04 6.90 E-03 GO:0016740 分子功能F 转移酶活性 303 4 740 1.70 E-04 8.00 E-03 GO:0006796 生物学功能P 磷酸盐代谢过程 163 2 319 6.60 E-05 9.60 E-03 GO:0006793 生物学功能P 磷代谢过程 163 2 319 6.60 E-05 9.60 E-03 GO:0051704 生物学功能P 多生物过程 18 129 1.10 E-04 1.40 E-02 GO:0050896 生物学功能P 刺激响应 87 1 124 1.30 E-04 1.60 E-02 GO:0009875 生物学功能P 花粉-雌蕊相互作用 16 122 5.00 E-04 4.60 E-02 GO:0008037 生物学功能P 细胞识别 16 122 5.00 E-04 4.60 E-02 GO:0048544 生物学功能P 花粉识别 16 122 5.00 E-04 4.60 E-02 GO:0009856 生物学功能P 授粉 16 122 5.00 E-04 4.60 E-02 Table 1. GO terms enriched by candidate gene (nsSNP)
-
‘六月雪’样本共得到INDEL 800 388个，编码区内移码插入总数为6 092个，移码缺失总数8 884个；编码区内因INDEL突变获得终止子426个，丢失终止子107个；在基因5'UTR内的INDEL总数为31个，在基因3'UTR内的INDEL总数和位于5'UTR和3'UTR间的INDEL总数均为0，在不同可变剪切的基因组区域内找到INDEL 786个，内含子区域内找到201 635个，启动子区域内找到198 924个，基因间区域内找到383 503个。
‘黄皮消’样本共得到INDEL 799 603个，编码区内移码插入总数为6 021个，移码缺失总数8 708个；编码区内因INDEL突变获得终止子440个，丢失终止子105个；基因5'UTR内找到INDEL 26个，基因3'UTR内、不同基因的5'UTR和3'UTR间则未找到；在不同可变剪切的基因组区域内找到INDEL 758个，内含子区域内找到199 949个，启动子区域内找到198 089个，基因间区域内找到385 507个。
-
2个样本共获得INDEL15 509和15 274个。对这些INDEL的差异分析发现，共有5 115个INDEL关联到了3 682个基因，其中24个终止子丢失（stop-loss）INDEL关联了24个基因，165个终止子获得（stop-gain）INDEL关联了160个基因，1 901个移码插入（frame-shift insertion）INDEL关联了1 629个基因，3 025个移码缺失（frame-shift deletion）INDEL关联了2 453个基因。分析发现1 276个基因内的INDEL数量超过1个；假定的抗病RPP13样蛋白1（PCP007357），烟草花叶病毒耐药蛋白N（PPCP015254），未注释蛋白1（PCP015680），烟草花叶病毒耐药蛋白N（PCP030478），假定的酰胺酶C869.01（PCP023678），ATP依赖的RNA解旋酶DHX36（PCP003694），甘露糖寡糖α-1, 2-甘露糖苷酶MNS1（PCP005093），未注释蛋白2（PCP017985），富有丝氨酸/精氨酸的分裂因子SC35（PCP000011），来自转座子逆转录病毒相关的Pol聚蛋白TNT 1-94（PCP032808），未注释蛋白3（PCP031973）等蛋白编码的基因均具有7个以上INDEL。
-
分析全局GO功能发现，INDEL关联的候选基因中，刺激响应、结合、催化活性等功能相对于背景基因（1 306条）来说更加富集；INDEL关联的差异基因的富集分析（图 5及表 2）结果与nsSNP关联的候选基因的GO富集分析一致，全部25 698条信息中。
Figure 4. Overall GO distribution of INDEL candidative gene
Figure 5. GO enrichment distribution map of INDEL candidative gene
GO数据库中唯一的标记信息 GO term类型 GO功能的描述信息 输入的具有GO term注释的候选基因数 该GO term中基因总数 富集分析统计学显著水平 错误发现率 GO:0043531 分子功能F ADP结合 60 415 9.20 E-13 9.60 E-10 GO:0004713 分子功能F 蛋白酪氨酸激酶活性 145 1 649 1.10 E-09 2.50 E-07 GO:0032559 分子功能F 腺苷核糖核酸结合 297 3 947 9.30 E-10 2.50 E-07 GO:0001883 分子功能F 嘌呤核苷结合 310 4 170 1.50 E-09 2.50 E-07 GO:0001882 分子功能F 核苷结合 313 4 181 6.50 E-10 2.50 E-07 GO:0030554 分子功能F 腺苷核苷酸结合 310 4 170 1.50 E-09 2.50 E-07 GO:0006952 生物学功能P 防御响应 52 402 1.50 E-09 2.10 E-06 GO:0004674 分子功能F 蛋白丝氨酸/苏氨酸激酶活性 118 1 331 1.90 E-08 2.80 E-06 GO:0016772 分子功能F 转移含磷基团的转移酶活性 213 2 779 2.90 E-08 3.80 E-06 GO:0032555 分子功能F 嘌呤核糖核酸结合 309 4 320 4.50 E-08 4.70 E-06 GO:0032553 分子功能F 核糖核酸结合 309 4 320 4.50 E-08 4.70 E-06 GO:0017076 分子功能F 嘌呤核苷酸结合 323 4 552 5.10 E-08 4.90 E-06 GO:0004672 分子功能F 蛋白激酶活性 153 1 995 1.50 E-06 1.30 E-04 GO:0003964 分子功能F RNA导向的DNA聚合酶活性 19 106 1.70 E-06 1.40 E-04 GO:0016773 分子功能F 以醇基为受体的磷酸转移酶活性 168 2 246 2.20 E-06 1.50 E-04 GO:0034061 分子功能F DNA聚合酶活性 22 138 2.20 E-06 1.50 E-04 GO:0006278 生物学功能P 依赖RNA的DNA复制 19 106 1.70 E-06 6.40 E-04 GO:0043687 生物学功能P 蛋白质翻译后修饰 170 2 265 1.60 E-06 6.40 E-04 GO:0006468 生物学功能P 蛋白氨基酸磷酸化 153 1 979 9.80 E-07 6.40 E-04 GO:0006464 生物学功能P 蛋白氨基酸磷酸化 178 2 429 3.80 E-06 1.10 E-03 GO:0016301 分子功能F 激酶活性 164 2 277 2.00 E-05 1.30 E-03 GO:0005524 分子功能F ATP结合 240 3 538 2.40 E-05 1.50 E-03 GO:0016779 分子功能F 核苷酸转移酶活性 36 335 2.70 E-05 1.60 E-03 GO:0005515 分子功能F 蛋白质结合 395 6 182 3.30 E-05 1.80 E-03 GO:0043412 生物学功能P 高分子修饰 180 2 513 1.20 E-05 2.30 E-03 GO:0006950 生物学功能P 胁迫响应 76 882 1.10 E-05 2.30 E-03 GO:0016310 生物学功能P 磷酸化 158 2 160 1.20 E-05 2.30 E-03 GO:0030246 分子功能F 碳水化合物结合 36 353 8.00 E-05 4.20 E-03 GO:0000166 分子功能F 核苷酸结合 362 5 741 1.40 E-04 6.90 E-03 GO:0016740 分子功能F 转移酶活性 303 4 740 1.70 E-04 8.00 E-03 GO:0006796 生物学功能P 磷酸盐代谢过程 163 2 319 6.60 E-05 9.60 E-03 GO:0006793 生物学功能P 磷代谢过程 163 2 319 6.60 E-05 9.60 E-03 GO:0051704 生物学功能P 多生物过程 18 129 1.10 E-04 1.40 E-02 GO:0050896 生物学功能P 刺激响应 87 1 124 1.30 E-04 1.60 E-02 GO:0009875 生物学功能P 花粉-雌蕊相互作用 16 122 5.00 E-03 4.60 E-02 GO:0008037 生物学功能P 细胞识别 16 122 5.00 E-03 4.60 E-02 GO:0048544 生物学功能P 花粉识别 16 122 5.00 E-03 4.60 E-02 GO:0009856 生物学功能P 授粉 16 122 5.00 E-03 4.60 E-02 Table 2. GO terms enriched by candidate gene (INDEL)
-
用VarScan软件对拷贝数变异（CNV）分析，共获得CNV 37 039个，其中缺失CNV（deletion CNV）20 629个，中性CNV（neutral CNV）7 577个，扩增CNV（amplification CNV）8 833个。
2.1. 测序数据预处理
2.2. 短序列匹配（reads mapping）统计
2.3. SNP统计分析
2.4. 常见变异（common variant，CV）分析
2.5. nsSNP关联的候选基因的基因本体论富集分析
2.6. INDEL统计分析
2.7. 相关INDEL差异分析
2.8. INDEL关联的候选基因的GO富集分析
2.9. 拷贝数变异（copy number variations，CNV）分析
-
近年来，迅猛发展的基因组测序技术已被广泛应用于植物基因组重测序等研究[27]。SNP是基因组DNA序列上广泛存在的最基本的变异形式[28]，植物基因组上平均每数百bp就存在一个SNP[29]。将SNP和传统分子标记相结合用于分子辅助育种，通过全基因组重测序（WGR）得到测序数据，与参考基因组进行序列比对，可分析SNP遗传变异信息，开发出数量较为丰富的分子标记，实现遗传资源的高效利用[30]。
全基因组重测序为从基因组水平开发SNP标记提供了新的技术条件，将SNP识别、验证和基因型分析与传统分子标记相结合，能快速挖掘到候选基因和获得导致表型的SNP位点[9]。中国6个玉米Zea mays优良自交系的非重复区大约存在1 272 134个SNP和30 178个INDEL，其中68 966个SNP和571个INDEL位于功能基因内[31]。对梁Setaria italica ‘SLX’品种的全基因组重测序发现，‘SLX’基因组存在762 082个SNP，26 802个INDEL，10 109个SV[32]；对包括野生种、培育种及改良品种在内的360份番茄Lycopersicon esculentum进行全基因组重测序，发现番茄的驯化和改良是2个相对独立的过程，影响果实颜色的主效基因是SlMYB12[33]。对302株大豆Glycine max进行全基因组重测序，检测到了162个受选择的拷贝数变异（CNV），并发现植株进化、发育性状与受选择区域相关[34]。大豆品种‘齐黄34’‘Qihuang34’经全基因组重测序检测到1 519 494个SNP，357 549个INDEL，4 506个SV，17 748个基因变异，其中转录、复制、重组、修复、信号传导机制等6个功能类序列存在较多的变异基因[35]。梨的全基因组重测序少见报道[11-13]，大多是通过转录组分析来发掘SNP位点并进行功能注释。编码区内的SNP一般可分为同义SNP（synonymous SNP）和非同义SNP（non-synonymous SNP），其中同义SNP所致的编码序列的改变不会引起氨基酸序列变化，而非同义SNP会使氨基酸序列发生改变，最终影响蛋白质序列[36]，因此认为非同义SNP是导致生物性状改变的直接原因，开发并研究这类SNP标记往往具有更为重要的生物学意义。本研究中，‘六月雪’和‘黄皮消’样本中总SNP数量分别为6 171 357和6 140 603个，编码区内无义突变内总数（nsSNP）分别为335 659和332 280个；对nsSNP的CV分析发现共有2 282个nsSNP关联了2 067个基因。2个样本的总INDEL数量分别为800 388和799 603个，位于编码区内的分别有15 509和15 274个；共5 115个差异INDEL关联到了3 682个基因，令烟草花叶病毒耐药蛋白N和抗病蛋白等关键基因发生变异。差异nsSNP基因和差异INDEL富集得到的GO terms一致。针对这些突变位点进行SNP和INDEL相关标记的开发、优异基因的挖掘将为分子标记辅助育种提供重要的标记资源，对上饶早梨‘六月雪’和‘黄皮消’育种研究具有重要的指导意义。
DownLoad: