基于基因組重測序的栗屬種質(zhì)資源遺傳特性分析

中圖分類號：S664.2 文獻(xiàn)標(biāo)志碼：A 文章編號：1002-2910（2025）02-0011-08

Analysis of the genetic characteristics of Castanea germplasms based on genome resequencing

WANG Jinping1， TIAN Shoule1， CAI Gongzhan2， SUN Xiaoli1， TAO Diyang3， SHEN Guangning1* （1.Shandong Institute of Pomology，Tai'an，Shandong 2710oo， China; 2.Shandong Shangnong Agricultural Technology Co.，Ltd.，Heze，Shandong，China;3.Culaishan Forest Farm ofTai'an City，Tai'an， Shandong ， China）

Abstract： In this study， genome resequencing technology was employed to conduct an analysis of the genetic characteristics of a total of 180 Castanea germplasms，which included Chinese chestnut （Castanea molissima） cultivars from 16 provinces in China， as wellas European chestnut （Castanea sativa） and Japanese chestnut （Castanea crenata）. As a result， 48 374 108 genetic variation sites were successfully identified. According to the analysis of these variation sites，the l80 cultivars could be classified into eight major groups. There was a phenomenon of gene introgresson among different groups. The population dynamic changes of Chinese chestnut， European chestnut， and Japanese chestnut during the evolutionary process were significantly different. The study also discovered the gene flow of the Chinese chestnut populations from Beijing to Shandong. However， the study has limitations such as uneven distribution of the sample size and limited analysis dimensions. In the future， it is necessary to expand the sample scope and conduct in-depth research by integrating multi-omics approaches.

Key words： chestnut; germplasms; genome re-sequencing; genetic characteristics; population structure; gene flow

板栗（CastaneamollissimaBlume）在中國農(nóng)林產(chǎn)業(yè)與生態(tài)建設(shè)中占據(jù)重要地位[1]。深入剖析其種質(zhì)資源的遺傳特性，是推動品種優(yōu)化、培育高產(chǎn)優(yōu)質(zhì)且適應(yīng)性強(qiáng)新品種的關(guān)鍵。以往關(guān)于板栗的遺傳分析多采用形態(tài)學(xué)標(biāo)記，通過觀察植株的形態(tài)特征，如樹形、葉形、果實(shí)形狀等，以及基于簡單序列重復(fù)（SSR）等分子標(biāo)記技術(shù)。這些方法雖能初步揭示板栗種群間的遺傳差異，但由于易受環(huán)境因素干擾，且所能檢測到的遺傳變異位點(diǎn)有限，無法全面、精確地展現(xiàn)板栗復(fù)雜的遺傳全貌，滿足現(xiàn)代遺傳育種研究的需求[2]。

基因組重測序能夠?qū)σ延袇⒖蓟蚪M的物種進(jìn)行深度測序，通過與參考基因組細(xì)致比對，能夠高效、精準(zhǔn)地識別出單核昔酸多態(tài)性（SNP）、短片段插入缺失（InDeI）以及大片段結(jié)構(gòu)變異（SV）等豐富的遺傳變異信息[3]，助力解析群體遺傳結(jié)構(gòu)、挖掘關(guān)鍵基因，極大地推動了植物遺傳育種進(jìn)程[4]。

盡管基因組重測序技術(shù)在植物研究中成果豐碩，但栗屬種質(zhì)資源遺傳特性分析方面的應(yīng)用研究報(bào)道較少。本研究運(yùn)用基因組重測序技術(shù)，對廣泛收集的栗屬種質(zhì)資源展開全面深入的遺傳分析，旨在明晰栗屬種質(zhì)資源的遺傳多樣性、明確群體結(jié)構(gòu)劃分以及梳理親緣關(guān)系脈絡(luò)，為栗屬種質(zhì)資源的保護(hù)利用和遺傳改良提供理論與數(shù)據(jù)支撐。

1材料與方法

1.1 試驗(yàn)材料

本研究收集了180份來自不同地區(qū)的栗屬種質(zhì)資源（表1）。這些種質(zhì)資源在形態(tài)特征、生態(tài)適應(yīng)性等方面具有一定的差異，能夠較好地代表栗屬種質(zhì)資源的多樣性。

1.2 試驗(yàn)方法

1.2.1DNA提取與質(zhì)量檢測采用改良的CTAB法提取板栗葉片的基因組DNA，通過 0 . 7 5 % 瓊脂糖凝膠電泳檢測DNA片段大小和DNA降解程度；

NanoDropOne分光光度計(jì)檢測DNA純度，檢測OD260/280比值在 1 . 8 ～ 2 . 2 ，無蛋白質(zhì)和肉眼可見雜質(zhì)污染；Qubit3.0熒光儀（LifeTechnologies，Carlsbad，CA，USA）檢測DNA濃度，檢測濃度大于 ，總量大于 2 μ g 。

1.2.2基因組重測序文庫構(gòu)建及庫檢DNA樣品檢測合格后，使用Covaris超聲波破碎儀將DNA隨機(jī)打斷成 3 0 0 ～ 5 0 0 b p 的片段，打斷后的樣品磁珠進(jìn)行片段選擇，使得樣品條帶集中在 2 0 0 ～ 4 0 0 b p 左右。再經(jīng)末端修復(fù)、加A尾、加測序接頭、純化、PCR擴(kuò)增、PCR產(chǎn)物環(huán)化等步驟完成整個(gè)文庫制備工作。

文庫構(gòu)建完成后，先使用Qubit2.0熒光定量儀進(jìn)行初步定量，隨后使用Agilent2100生物分析儀對文庫的插入片段進(jìn)行檢測，插入片段大小符合預(yù)期后，使用Q-PCR方法對文庫的有效濃度進(jìn)行準(zhǔn)確定量，以保證文庫質(zhì)量。測序時(shí)，采用雙端150bp的測序策略，測序深度達(dá)到10X以上，以保證數(shù)據(jù)的準(zhǔn)確性和可靠性。

DNBSEQ上機(jī)測序。檢測合格文庫安排上機(jī)測序（DNBSEQ）：單鏈環(huán)狀DNA分子通過滾環(huán)復(fù)制，形成一個(gè)包含300多個(gè)拷貝的DNA納米球（DNB）。將得到的DNBs采用高密度DNA納米芯片技術(shù)，加到芯片上的網(wǎng)狀小孔內(nèi)，通過聯(lián)合探針錨定聚合技術(shù)（cPAS）進(jìn)行測序。

1.2.3測序數(shù)據(jù)篩選測序平臺得到的原始圖像數(shù)據(jù)文件，經(jīng)過堿基識別分析即可轉(zhuǎn)化為原始測序序列，稱為RawData或RawReads，對Rawreads過濾，得到Cleanreads。使用fastp[5]對原始數(shù)據(jù)進(jìn)行質(zhì)控與過濾，具體過濾標(biāo)準(zhǔn)如下：剔除序列中的接頭序列，即adapter；去除reads尾部的polyG和polyX；以滑窗的方式統(tǒng)計(jì)窗口內(nèi)的堿基計(jì)算平均質(zhì)量值，將低質(zhì)量的滑窗剪裁掉；剔除N個(gè)數(shù)大于5的reads；剔除質(zhì)量低于15的堿基占比高于 40 % 的reads；剔除過濾后長度低于 1 5 b p 的reads。

1.2.4變異鑒定利用生物信息學(xué)軟件將測序數(shù)據(jù)與板栗參考基因組[（表2）進(jìn)行比對。首先，使用BWA[7軟件將測序reads比對到參考基因組上（參數(shù)：mem-R，其余參數(shù)采用軟件默認(rèn)參數(shù)），然后利用SAMtools[8軟件對比對結(jié)果進(jìn)行排序（參數(shù)：sort）、去重（參數(shù)：markdup-r）等處理，將比對結(jié)果從sam（SequenceAlignment/MAP）文件轉(zhuǎn)為排序后的bam文件（binaryAlignment/Map），統(tǒng)計(jì)比對情況。為了更具體更細(xì)致的評估各樣本比對參考基因組上的情況，采用滑窗的方式將基因組劃分成大小相等的子區(qū)間，統(tǒng)計(jì)各個(gè)子區(qū)間內(nèi)的平均測序深度（比對上該區(qū)間總的堿基數(shù)除以區(qū)間大?。┮约案采w度（區(qū)間內(nèi)有reads覆蓋的區(qū)域占區(qū)間大小的比例），從而更全面的了解各個(gè)樣本的比對情況。

由于二代測序數(shù)據(jù)reads長度較短且基因組上存在著重復(fù)序列從而導(dǎo)致reads的錯(cuò)誤比對，同時(shí)測序數(shù)據(jù)的分布不均勻。這些原因可能會導(dǎo)致變異鑒定結(jié)果的假陽性存在。因此為了降低這一比例，獲得高質(zhì)量的SNP和INDEL，使用GATK[9]（版本：4.2.5.0；參數(shù)：VariantFiltration），根據(jù)GATK官方推薦的hard-filtering標(biāo)準(zhǔn)對鑒定到的SNP和Indel分別進(jìn)行過濾，具體過濾標(biāo)準(zhǔn)如下：

基于GATKhard-filtering過濾后的剩下的變異結(jié)果文件，使用vcftools（版本：0.1.16；參數(shù)：-maf，-max-missing），剔除掉maf（次要等位基因頻率）小于0.03，基因組型缺失比例大于 20 % 的位點(diǎn)。最后使用剩余的變異位點(diǎn)進(jìn)行群體結(jié)構(gòu)分析。

1.2.5群體結(jié)構(gòu)分析使用PHYLIP[10]（版本：3.696;參數(shù)：neighbor）中的鄰接法（neighbor－joiningmethods，簡稱NJ），構(gòu)建進(jìn)化樹。后續(xù)基于樹文件（newick格式）使用ggtree進(jìn)行可視化。

使用fastStructure[11，2]（版本：1.0；參數(shù)：-K- seed - c v= 1 0 ）進(jìn)行群體遺傳結(jié)構(gòu)分析，K值取2～10，基于不同的seed值針對每個(gè)K值獨(dú)立重復(fù)10次（最佳K值根據(jù) fastStructure 提供的chooseK.py程序以及具體的群體背景進(jìn)行確定）。

使用PSMC[13] （pairwise sequentially Markoviancoalescent）（版本：latest；參數(shù)： - N2 5 - t1 5 - r5 -b-p ，堿基突變速率： 每代年數(shù)：15），利用個(gè)體重測序數(shù)據(jù)推測該個(gè)體所屬的種群在歷史上各個(gè)時(shí)期的有效群體大小。

利用TreeMix[14]軟件通過從多個(gè)種群中獲得等位基因頻率，返回該種群的最大似然（ML）樹，并推斷基因流。

2 結(jié)果與分析

2.1栗屬種質(zhì)資源統(tǒng)計(jì)

根據(jù)來源對重測序栗屬資源進(jìn)行了分析，180份栗屬種質(zhì)資源分別來源于中國16個(gè)?。ㄊ校┮约皻W洲和日本國家，其中山東49份、河北13份、安徽11份、湖南15份、湖北16份、江蘇24份、浙江14份、河南9份、陜西8份、北京6份、云南3份、廣東3份、廣西3份、福建2份、江西1份、貴州1份、歐洲栗1份、日本栗1份（表1）。

2.2 測序數(shù)據(jù)質(zhì)控

通過對180份栗屬資源的測序，每個(gè)樣本得到的Rawreads數(shù)量 6 3 0 1 0 8 6 6 ～ 1 8 8 5 7 2 2 0 4 ，通過進(jìn)行質(zhì)控與過濾獲得Cleanreads的數(shù)量范圍為： ratio在 9 7 . 2 3 % 以上，Q30ratio在 9 0 . 2 0 % 以上，GC含量 34 . 3 9 %～3 6 . 8 0 % 。以2020年發(fā)表的板栗全基因組組裝結(jié)果為參考基因組，基因組大小為6 8 8 . 9 9 M b 。

通過對比對結(jié)果分析，180個(gè)樣本的比對率在9 7 . 6 2 % ～ 9 9 . 4 6 % ，對參考基因組的覆蓋度在

，平均測序深度在13.11X～3 2 . 4 3 X 。

2.3等位基因類型及統(tǒng)計(jì)

基于堿基類型信息，對48374108個(gè)變異位點(diǎn)進(jìn)行了分類，其中二等位變異、三等位變異、四等位變異、五等位變異、六等位變異、七等位變異數(shù)量分別為43564584、4015291、405120、130080、91839、167194。對每個(gè)品種具體的變異情況（INDEL變異、堿基轉(zhuǎn)換、顛換、顛換轉(zhuǎn)換比、純合基因型、雜合基因型以及雜合比等）進(jìn)行了統(tǒng)計(jì)。整體變異中各類型顛換以及轉(zhuǎn)換變異的比例進(jìn)行了可視化（圖2），針對變異在基因組上的分布進(jìn)行可視化（圖3）。

2.4群體結(jié)構(gòu)分析

群體進(jìn)化樹構(gòu)建。根據(jù)所測180份資源的重測序結(jié)果，基于鄰接法構(gòu)建了系統(tǒng)進(jìn)化樹，180份資源可以分成八大類群，大的聚類結(jié)果與品種地理來源沒有明顯相關(guān)。

如來自于湖南的z-305（沅陵2號）、z-365（沅陵1號），江蘇的z-280（句容短刺）、z-385（青毛軟扎）、z-373（炮車7號）、Z-101（青扎）、z-103（焦扎），湖北的 z - 0 9 3 （鄂栗1號），山東的z-251（輻3）、z-239（青毛軟刺）、z-284（魯35）、z-341（沭河7號），河南的z-353（磚橋處暑紅），安徽的z-042（大紅袍）、z-094（舒城大紅袍），河北的 z - 2 9 0 （西溝2號）等16個(gè)品種聚為一類。有部分小規(guī)模聚類值得關(guān)注，如來自廣東的z-271（韶栗18號）、Z-363（早香1號）聚為一類；來自江蘇的z-019（優(yōu)選處暑紅）與來自安徽的z-055（處暑紅）聚到一起，來自湖北的z-046（紅毛早）和z-374（淺刺紅毛早）聚到一起；河北的 z - 0 0 7 （東密塢無花）與z-201（貴州野毛栗）聚到一起；z-198（廣西14-5）、z-011（高店10號）、z-324（合肥大紅袍）z-320（白毛栗）與Shenyinou（沈引歐）聚為一類。

群體Structure分析?；跈z測過濾后的品種SNP數(shù)據(jù)，使用fastStructure進(jìn)行群體遺傳結(jié)構(gòu)分析，當(dāng) K=8 時(shí)Cv-Error較小，與系統(tǒng)發(fā)育樹分群結(jié)果相一致，各來源地表現(xiàn)出基因結(jié)構(gòu)相互滲透的現(xiàn)象（圖4）。因此認(rèn)為收集的180份樣本分為8個(gè)類群比較合適。

2.5 有效群體大小分析

本研究從山東、河北、河南、江蘇、安徽、日本栗、歐洲栗等群體中分別選擇1個(gè)樣本的重測序數(shù)據(jù)推測個(gè)體所屬的種群在歷史上各個(gè)時(shí)期的有效群體大小。如圖5所示，除廣東外，中國各個(gè)產(chǎn)區(qū)的有效群體變化趨勢相近，而日本栗有效群體變化趨勢則與中國栗差異較大。由圖可知，大約在1 Ma.B.P.（Million-anniversary Before Present）中國栗、日本栗、歐洲栗三個(gè)群體大小大致重合在一起。中國栗在 0 . 4 ～ 0 . 5 M a . B . P 有效群體達(dá)到最大值，日本栗在0.2～0.3Ma.B.P.有效群體達(dá)到最大值，且有效群體數(shù)量明顯高于中國栗。此后經(jīng)歷了長期的衰退。大約在 2 0 ～ 2 5 K a .B.P.（kilion -anniversaryBeforePresent）又開始逐漸擴(kuò)張。

2.6 基因流分析

根據(jù)各群體中每個(gè)SNP類型變異位點(diǎn)等位基因數(shù)量統(tǒng)計(jì)，以每?。ㄊ校┑钠贩N為一個(gè)種群，使用Treemix進(jìn)行基因流預(yù)測。由于廣西、福建、廣東等省份的品種數(shù)量較少，不進(jìn)行分析。對湖北、湖南、山東、陜西、河北、北京、河南、安徽、浙江、江蘇等10個(gè)?。ㄊ校┑钠贩N進(jìn)行基因流分析。如圖6所示，河南、安徽、浙江、江蘇、湖北、湖南漂移參數(shù)較低，群體間遺傳差異?。粰z測出從北京板栗種群向山東板栗種群的基因流，反映出山東早期從北京引種交流活動比較頻繁。

3小結(jié)與討論

3.1群體結(jié)構(gòu)分析在育種中的利用

本研究首次運(yùn)用基因組重測序技術(shù)對來自中國、歐洲和日本共180份栗屬種質(zhì)資源展開了全面的遺傳特性分析，成功獲取了高質(zhì)量的測序數(shù)據(jù)，并精準(zhǔn)鑒定出豐富的遺傳變異位點(diǎn)。

通過系統(tǒng)進(jìn)化樹和Structure等群體結(jié)構(gòu)分析表明，180份資源可劃分為八大類群，且各地區(qū)品種間存在明顯的基因滲透現(xiàn)象，這充分證實(shí)了中國板栗資源遺傳背景的多樣性以及不同地區(qū)品種間頻繁的基因交流。有效群體大小分析揭示了中國栗、日本栗和歐洲栗在進(jìn)化歷程中的群體動態(tài)變化，為深入理解板栗的進(jìn)化歷史提供了關(guān)鍵線索?；蛄鞣治鲞M(jìn)一步解釋了不同省份板栗群體間的遺傳關(guān)系，為板栗品種的遺傳改良提供了信息。

研究鑒定出的大量遺傳變異位點(diǎn)，為后續(xù)開展板栗功能基因挖掘、分子標(biāo)記開發(fā)、進(jìn)化生物學(xué)研究以及深入解析板栗重要農(nóng)藝性狀的遺傳機(jī)制提供了豐富的素材。

板栗在相互引種和選育過程中出現(xiàn)了同名異物或者同物異名的現(xiàn)象，如本研究通過分析發(fā)現(xiàn)優(yōu)選處暑紅與處暑紅，紅毛早與淺刺紅毛早等很可能為同物異名。明確的群體結(jié)構(gòu)和遺傳關(guān)系，能夠指導(dǎo)種質(zhì)資源的合理收集、保存和管理，避免遺傳資源的重復(fù)收集和流失。育種工作者可依據(jù)遺傳距離和基因流信息，精準(zhǔn)選擇具有較大遺傳差異的親本進(jìn)行雜交，從而提高育種效率，培育出更具優(yōu)良性狀的新品種，有力推動板栗產(chǎn)業(yè)的高質(zhì)量發(fā)展。

3.2 研究展望

盡管本研究取得了一系列重要成果，但仍存在一定的局限性。在試驗(yàn)材料方面，雖然收集了來自多個(gè)地區(qū)的栗屬種質(zhì)資源，但歐洲栗、日本栗以及中國的廣西、云南等部分地區(qū)的樣本數(shù)量相對較少，可能對群體遺傳結(jié)構(gòu)分析的準(zhǔn)確性產(chǎn)生一定影響。未來研究可進(jìn)一步擴(kuò)大樣本收集范圍，增加樣本數(shù)量，以更全面地揭示栗屬種質(zhì)資源的遺傳多樣性。

在數(shù)據(jù)分析方面，本研究主要聚焦于群體結(jié)構(gòu)、有效群體大小和基因流等常規(guī)分析，后續(xù)研究將結(jié)合表型數(shù)據(jù)，運(yùn)用全基因組關(guān)聯(lián)分析（GWAS）等方法，深入挖掘與板栗產(chǎn)量、品質(zhì)、抗病性等重要性狀緊密相關(guān)的遺傳標(biāo)記和功能基因，為分子標(biāo)記輔助育種提供更精準(zhǔn)的靶點(diǎn)。

此外，本研究僅涉及基因組層面的分析，而植物的生長發(fā)育和性狀表現(xiàn)是一個(gè)復(fù)雜的調(diào)控網(wǎng)絡(luò)，受到轉(zhuǎn)錄組、蛋白質(zhì)組和代謝組等多個(gè)層面的協(xié)同調(diào)控。未來將結(jié)合開展多組學(xué)聯(lián)合分析，多個(gè)維度解析板栗的遺傳調(diào)控機(jī)制，為栗屬種質(zhì)資源的創(chuàng)新利用和遺傳改良拓寬路徑。

參考文獻(xiàn)：

[1］闞黎娜，李倩，謝爽爽，等．我國板栗種質(zhì)資源分布及營養(yǎng)成分比較［J]．食品工業(yè)科技，2016，37（20）：396-400.

［2］李沛，何治霖，談月霞，等．基于重測序數(shù)據(jù)與表型性狀的寬皮柑橘遺傳多樣性分析與優(yōu)異種質(zhì)篩選［J].中國農(nóng)業(yè)科學(xué)，2024，57（23）：4761-4795.

［3］姚遠(yuǎn)，鄧?yán)?，胡娟，等．‘脆紅李’及其早熟芽變?nèi)蚪M重測序分析［J]．園藝學(xué)報(bào)，2024，51（10）：2255-2266.

[4]VARSHNEYRK，GRANERA，SORRELLSME.Genomicsequencing and discovery of markers for crop improvement[J]．Trends in biotechnology，2009，27（7）：442-450.

[5]CHEN S， ZHOU Y，CHEN Y，et al. Fastp： an ultra- fast all- in-oneFASTQ preprocessor[J].Bioinformatics，2018，34（17）：1884 -1890.

[6]WANG J， TIAN S，SUN X，et al. Construction of Pseudomoleculesfor the Chinese Chestnut （Castanea mollissima） Genome[J].G3-Genes GenomesGenetics，2020，10（10）：3565-3574.

[7]HENG L， RICHARD D. Fast and accurate short read alignment withBurrows-Wheeler transform[J].Bioinformatics，2010，14（25）：1754-1760.

[8]LI H，HANDSAKER B，WYSOKER A，et al.The SequenceAlignment/Map format and SAMtools[J].Bioinformatics，2009，25（16）：2078-2079.

[9]DANECEK P，AUTONA，ABECASIS G，et al.The variant allformatandVCFtools[J].Bioinformatics，201l，27（15）：2156-2158.

[10］RAJA，STEPHENS M，PRITCHARD JK. Fast STRUCTURE：Variational Inference of Population Structurein Large SNPDataSets[J].Genetics，2014，197（2）：573-589.

[11]LAM HM， XU X， LIU X，et al. Addendum： Resequencing of 31wild and cultivated soybean genomes identifies paterns of geneticdiversity and selection[J].Nature Genetics，2011，43（4）：387.

[12］LAYER R M，CHIANGC，QUINLAN AR，et al. LUMPY：aprobabilisticframeworkforstructuralvariantdiscovery[J]．GENOMEBIOLOGY，2014，15（6）：1-19.

[13]SHENGLIN，LIU，MICHAEL，et al. PSMC （pairwise sequentiallyMarkovian coalescent） analysis of RAD （restriction site associatedDNA） sequencing data[J]. Molecular Ecology Resources，2017，（17）：631-641.

[14]FITAK R R. Opt M： estimating the optimal number of migrationedges on population trees using Treemix[J].Biology Methods andProtocols，2021，6（1）：1-10.