楊 赟 陳夢嬌 杜慶鑫 朱景樂 杜紅巖 楊紹彬*

(1.國家林業(yè)和草原局泡桐研究開發(fā)中心,經(jīng)濟林種質創(chuàng)新與利用國家林業(yè)和草原局重點實驗室，鄭州 450003； 2.南京林業(yè)大學林學院，南京 210000； 3.河南農(nóng)業(yè)大學林學院，鄭州 450002)

杜仲(Eucommiaulmoides)隸屬于杜仲科(Eucommiaceae)，單屬單種，染色體2n=2X=34，是我國特有的經(jīng)濟林樹種，在我國北緯22°～42°，東經(jīng)100°～121°的廣大范圍內(nèi)均可栽培[1]。杜仲因其適應性強，干型通直，冠型優(yōu)美，主根較深，側根發(fā)達，具有較好的水土保持能力，被廣泛用于園林綠化和防止水土流失[2～3]。當前，對于杜仲的研究主要集中在藥理功效、成分分析、栽培模式，干旱脅迫響應等方面[4～6]，對選育觀賞型杜仲良種報道較少?！t葉’杜仲(Eucommiaulmoides‘Hongye’)是中國林業(yè)科學研究院選育的觀賞型良種，葉片鮮紅至紫紅色，是理想的綠化樹種。‘紅葉’杜仲葉片呈紅色的主要物質是花色苷[7]，其葉片花色苷含量顯著高于普通杜仲，但是‘紅葉’杜仲葉片的呈色機理尚不清楚。

隨著測序技術的不斷發(fā)展，擁有參考基因組的物種越來越多，全基因組重測序技術在發(fā)掘與植物表型相關的單核苷酸多態(tài)性(SNP)位點方面發(fā)揮著重要作用。研究者利用該技術，快速得到這一物種某個體的基因組信息和重要功能基因的序列信息，初步定位和鑒定差異位點，為從分子水平研究植物基因調(diào)控提供了依據(jù)[8～10]。植物葉片的呈紅色的物質是花色苷，花色苷代謝途徑是目前研究的最透徹的通路之一?；ㄉ沼杀奖彼峤?jīng)過一系列的酶催化合成，由上游基因群和下游基因群控制酶的合成，上游基因群通常編碼參與黃酮、黃酮醇和糅紅合成相關的酶，下游基因群編碼合成花青苷和原花青素的酶[11]。本研究通過重測序技術檢測并分析控制‘紅葉’杜仲與普通綠葉杜仲葉色差異的基因位點，比較與葉色形成相關的酶基因位點差異，開發(fā)杜仲紅葉性狀SNP標記。

1 研究方法

1.1 研究材料

‘紅葉’杜仲(Eucommiaulmoides‘Hongye’)和‘小葉’杜仲(Eucommiaulmoides‘Xiaoye’)均是在河南洛陽地區(qū)通過實生選育方法獲得的無性系良種。本研究以3年生‘紅葉’杜仲和‘小葉’杜仲嫁接苗為對象(砧木均為普通杜仲)，每個無性系選擇3個單株，每個單株選擇頂端中上部的3片葉片，按無性系將葉片分組后置于液氮中保存。

1.2 DNA提取

按照Qiagen公司的DNeasy Plant Mini Kit試劑盒的操作步驟提取杜仲DNA，使用Nano Drop 2000進行DNA純度檢測，Picogreen方法檢測DNA濃度。1%的瓊脂糖凝膠電泳檢測DNA完整性，檢測前將樣品放在冰上融化后，充分混勻并離心。每個無性系3個單株，從中取等量的質量較好的DNA混合建庫，用于二代測序。

1.3 全基因組重測序

將檢驗合格的基因組DNA樣品(樣品量>2 μg，樣品濃度>50 ng·μL-1，無色素污染，無明顯RNA、蛋白質等雜質污染)，送至上海美吉生物醫(yī)藥科技有限公司進行文庫構建。利用超聲波將DNA序列片段形成隨機片段，經(jīng)過處理后用磁珠吸附富集基因組長度為400 bp左右的片段，經(jīng)過PCR擴增形成測序文庫。質檢合格的文庫用Illumina HiSeqTM平臺進行基因組重測序，測序策略為Illumina PE150，總測序讀長為300 bp。

1.3.1 原始測序數(shù)據(jù)質控與過濾

原始測序數(shù)據(jù)(Raw Data)經(jīng)過Fastp(v0.11.5)進行質控，包括堿基含量分布和堿基錯誤率分布。堿基含量分布檢查有無AT、GC分離現(xiàn)象。測序堿基序列錯誤率結果由于Illumina HiseqTM測序的技術特點，測序片段前端幾個Cycles和末端的錯誤率會偏高。使用NGSQCtoolkit(v2.3.3)套件對raw reads去除接頭等對reads進行trimming，去除低質量read并剪切掉5′端及3′端3 bp堿基并過濾掉長度小于50 bp的reads。最終得到質量較好的Clean Data。測序數(shù)據(jù)過濾采用Fastp默認參數(shù)流程[12]。

1.3.2 變異位點檢測及注釋

對于‘紅葉’杜仲和‘小葉’杜仲測序獲得的Clean Data，利用BWA(v0.7.16)基于MEM算法[13]將其比對到杜仲參考基因組序列上(杜仲參考基因組葉片為綠色)，利用GATK(v3.8.0)對BAM文件進行校正[14]，并進行SNP檢測。利用GATK軟件的Haplotyper算法進行變異位點檢測。利用SnpEff(v4.3)進行SNP變異位點功能注釋[15]，包括變異位點位置和功能注釋；SnpEff會根據(jù)變異位點對蛋白編碼的影響，從高到低依次分為High、Moderate、Low和Modifier 4個等級，High代表嚴重影響功能的SNP個數(shù)，Moderate代表中度影響功能突變的SNP；Low代表低影響程度；Modifier表示與功能無關的SNP。因此，只關注High和Moderate的個數(shù)。

1.3.3 變異基因的功能注釋和代謝通路分析

利用Diamond軟件比對到NR、UniProt和Gene Ontology(GO)數(shù)據(jù)庫，進行基因功能預測。使用KOBAS軟件比對到KEGG數(shù)據(jù)庫，注釋基因所涉及的代謝通路。利用emappe軟件比對EggNog數(shù)據(jù)庫，注釋基因功能。

1.3.4 變異位點的篩選

首先根據(jù)變異位點對蛋白編碼的影響，構建導致蛋白發(fā)生High、Moderate變化的基因集。根據(jù)變異基因的功能注釋，篩選與類黃酮和花色苷合成途徑相關的基因集。為了保證SNP位點的準確性，過濾掉低質量的SNP位點，要求SNP深度≥10X。選取‘紅葉’杜仲與參考基因組完全不一致的純合位點，‘小葉’杜仲與參考基因組完全一致的純合位點，這些位點被用于SNP位點開發(fā)。

1.4 基于Sanger測序SNP驗證

參考杜仲基因組序列信息，采用Primer Premier 5設計引物，退火溫度在60°左右，引物長度18～30 bp。PCR擴增產(chǎn)物長度：外顯檢測引物離外顯子上下游150 bp左右；目的基因測序的PCR產(chǎn)物條帶一般不超過1 200 bp。PCR反應體系：50 μL：1～2 μL的50 ng·μL-1cDNA，上下游引物、dNTP(mix)均2 μL(10 mmol·L-1), 10×Taq Buffer(with MgCl2) 5 μL，Taq酶(5 U·μL-1)0.5 μL，加雙蒸水至50 μL。PCR反應條件：95℃預變性3 min; 然后94℃變性30 s，60℃退火30 s，72℃延伸90 s，共35個循環(huán);最后72℃延伸8 min，4℃結束。Sanger測序由上海生工生物工程公司完成。

2 結果與分析

2.1 原始測序數(shù)據(jù)統(tǒng)計及質控

本次測序質量較好。‘紅葉’杜仲獲得14.16 G數(shù)據(jù)量，47 032 423個Clean reads，其中堿基質量Q30比例達到91.72%，GC含量為35.54%?！∪~’杜仲獲得14.29 G數(shù)據(jù)量，47 476 746個Clean reads，其中堿基質量Q30比例達到91.93%，GC含量為35.43%。數(shù)據(jù)質控分析顯示，除了reads的前幾個堿基存在正常的不平衡現(xiàn)象外，堿基質量分布基本無AT、GC分離現(xiàn)象，測序堿基錯誤率分布統(tǒng)計低于0.001%。

2.2 變異位點檢測結果

2.2.1 比對參考基因組數(shù)據(jù)統(tǒng)計

杜仲參考基因組大小為1.23 Gb，組裝水平是Scafflod(Scafflod N50=1.88 Mb)。比對參考基因組，‘紅葉’杜仲和‘小葉’杜仲測序數(shù)據(jù)的Mapped Ratio均為99.90%。本研究兩樣本對參考基因組(排除N區(qū))的平均覆蓋深度為10×左右，1×覆蓋深度(至少有一個堿基的覆蓋度)為87.77%以上。本研究覆蓋到參考基因組的區(qū)域較多，可檢測到的變異位點較多，基因組上堿基的覆蓋深度分布較均勻，測序的隨機性好。

2.2.2 SNP的檢測

在‘小葉’杜仲中共檢測到5 711 141個SNP，其中發(fā)生轉換(Ti)的SNP有4 139 309個，發(fā)生顛換(Tv)的SNP有1 565 071個，Ti/Tv比值為2.64，在‘紅葉’杜仲中共檢測到5 407 823個SNP，其中發(fā)生轉換(Ti)的SNP有3 915 050個，發(fā)生顛換(Tv)的SNP有1 487 022個，Ti/Tv比值為2.63。

2.2.3 SNP差異位點對蛋白編碼影響的評估

在‘小葉’杜仲中存在嚴重影響蛋白質功能的SNP為1 640個，中度影響功能的SNP為47 192個。在‘紅葉’杜仲中存在嚴重影響功能的有1 516個SNP，含有中度影響的41 328個SNP。

2.2.4 SNP差異位點功能注釋

SnpEff每個變異類型的功能統(tǒng)計如表1所示。某一基因在‘紅葉’杜仲和‘小葉’杜仲間存在堿基差異?！t葉’杜仲和‘小葉’杜仲每個SNP類型在數(shù)量上差異不大，可能它們分布在不同基因上，導致了‘紅葉’杜仲和‘小葉’杜仲在表型上的差異。

2.3 變異基因的功能分類

2.3.1 變異基因的NR、Uniprot、GO、KEGG、EggNOG注釋

含有變異位點的26 722基因，在Nr數(shù)據(jù)庫注釋到13 408條基因，在Uniport中注釋到14 418條基因，在GO數(shù)據(jù)庫中注釋到21 453條基因，在EggNOG數(shù)據(jù)庫中注釋22 660條基因，在KEGG數(shù)據(jù)庫中注釋到5 915條基因。為了進一步研究含有突變位點的基因功能，向GO數(shù)據(jù)庫的各個term進行映射(圖1)，發(fā)現(xiàn)在細胞組分(cellular component)中注釋到核、細胞質、膜的有機組成部分，在生物進程(biological process)中注釋到等離子體膜，在分子功能(molecular function)中注釋到蛋白結合、線粒體和葉綠體、ATP結合及轉錄過程。根據(jù)KEGG注釋信息，如表2所示，含有差異位點的基因可注釋到286條代謝途徑，其中較多的是淀粉和蔗糖代謝，苯丙素生物合成途徑等，與植物葉片生長及顏色密切相關。利用EggNOG數(shù)據(jù)庫對比分析，如圖2所示，有12 977個位點注釋到25個類別，其中信號轉導機制(1 220,9.38%)、轉錄(868,6.68%)、蛋白質翻譯后修飾與轉運、分子伴侶(838,6.45%)碳水化合物的運輸和代謝(801,6.16%)等類別注釋到的變異位點多，其次是能量產(chǎn)生和轉換(527,4.05%)，次生代謝產(chǎn)物的生物合成、轉運和分解代謝(456,3.51%)。

圖1 GO功能富集分類前20位統(tǒng)計圖 橫坐標表示位于功能區(qū)屬于High和Moderate類型的變異位點的個數(shù)，縱坐標表示GO Term的ID。Fig.1 GO function enrichment classification chart The abscissa indicates the number of variant sites located in the functional area belonging to the High and Moderate types,and the ordinate indicates the ID of the GO Term.

圖2 COG功能分類圖 A. RNA加工和修飾；B.染色質結構和動力學；C.能源生產(chǎn)和轉換；D.細胞周期控制，細胞分裂，染色體分配；E.氨基酸轉運和代謝；F.核苷酸轉運和代謝；G.碳水化合物轉運和代謝；H.輔酶轉運和代謝；I.脂質轉運和代謝；J.翻譯，核糖體結構和生物發(fā)生；K.轉錄；L.復制，重組和修復；M.細胞壁/膜/包膜生物發(fā)生；N.細胞運動性；O.翻譯后修飾，蛋白質周轉，分子伴侶；P.無機離子轉運和代謝；Q.次生代謝產(chǎn)物的生物合成，分解代謝；S.功能未知；T.信號轉導機制；U.細胞內(nèi)運輸，分泌和囊泡運輸；V.防御機制；W.細胞外結構；Y.核結構；Z.細胞骨架Fig.2 COG function classification chart A. RNA processing and modification; B. Chromatin structure and dynamics; C. Energy production and conversion; D. Cell cycle control,cell division,chromosome partitioning; E. Amino acid transport and metabolism; F. Nucleotide transport and metabolism; G. Carbohydrate transport and metabolism; H. Coenzyme transport and metabolism; I. Lipid transport and metabolism; J. Translation,ribosomal structure and biogenesis; K. Transcription; L. Replication,recombination and repair; M. Cell wall/membrane/envelope biogenesis; N. Cell motility; O. Posttranslational modification,protein turnover,chaperones; P. Inorganic ion transport and metabolism; Q. Secondary metabolites biosynthesis,and catabolism; S. Function unknown; T. Signal transduction mechanisms; U. Intracellular trafficking,secretion,and vesicular transport; V. Defense mechanisms; W. Extracellular structures; Y. Nuclear structure; Z. Cytoskeleton

表1 SNP突變位點功能注釋統(tǒng)計表

Table 1 SNP mutation site function annotation statistics table

類型Type‘紅葉’杜仲E.ulmoides‘Hongye’‘小葉’杜仲E.ulmoides‘Xiaoye’5′ 端非翻譯區(qū)5_prime_UTR_variant53下游Downstream707591753545基因間區(qū)Intergenic region48268425097746內(nèi)含子區(qū)Intron variant503340532407外顯子區(qū)域的錯義突變Missense variant4132843402非編碼轉錄本區(qū)Non coding transcript variant1810內(nèi)含子左側連接區(qū)Splice acceptor variant188197內(nèi)含子右側連接區(qū)Splice donor variant124149距離內(nèi)含子或外顯子2bp的剪切區(qū)Splice region variant55275952突變導致啟動子缺失Start lost165181突變導致終止子獲得Stop gained753807突變導致終止子突變Stop lost287307突變導致終止子的編碼改變Stop retained variant8288同義突變Synonymous variant4051742898上游Upstream gene variant759898809593

2.3.4 用于分子標記開發(fā)的SNP位點篩選

根據(jù)基因功能注釋的5大數(shù)據(jù)庫結果，26 722條含有差異位點的基因中共有228條基因注釋到與花色苷和類黃酮合成相關。228條基因內(nèi)部共有829個差異位點，根據(jù)SnpEff注釋，46個位點嚴重影響蛋白質結構，783個位點中度影響蛋白質功能。與參考基因組比，‘紅葉’杜仲中存在的這228條基因，含有27個位點嚴重影響蛋白質結構，其中有15個純合位點；有556個位點中度影響蛋白質結構，269個純合位點。經(jīng)過最終篩選，得到12個位點在‘紅葉’杜仲中完全差異于參考基因組，且‘小葉’杜仲中的這12個位點與參考基因組一致，測序深度均≥10×(表3)。這12個位點分布在6條基因的內(nèi)部(表4)。

表2 杜仲重測序的KEGG代謝途徑分類圖

Table 2 Classification of KEGG metabolic pathways forEucommiaulmoidesresequencing

編號Number途徑PathwayDifferent sites with KEGG annotation(12997)ko ID1淀粉和蔗糖代謝Starch and sucrose metabolism183(1.41%)ko005002植物—病原互作Plant-pathogen interaction164(1.26%)ko046263苯丙素的生物合成Phenylpropanoid biosynthesis128(0.98%)ko009404碳代謝Carbon metabolism125(0.96%)ko012005植物激素信號轉導Plant hormone signal transduction105(0.81%)ko040756氨基酸的生物合成Biosynthesis of amino acids103(0.79%)ko012307核糖體Ribosome98(0.75%)ko030108內(nèi)質網(wǎng)蛋白加工Protein processing in endoplasmic reticulum98(0.75%)ko041419戊糖和葡萄糖醛酸酯互變Pentose and glucuronate interconversions92(0.71%)ko0004010氨基糖和核苷酸糖代謝Amino sugar and nucleotide sugar metabolism92(0.71%)ko0052011RNA運輸RNA transport73(0.56%)ko0301312糖酵解和糖質新生Glycolysis/Gluconeogenesis72(0.55%)ko0001013嘌呤代謝Purine metabolism67(0.52%)ko0023014剪接體Spliceosome64(0.49%)ko0304015真核生物核糖體生物起源Ribosome biogenesis in eukaryotes60(0.46%)ko0300816半乳糖代謝Galactose metabolism59(0.45%)ko0005217內(nèi)吞作用Endocytosis59(0.45%)ko0414418氰氨基酸代謝Cyanoamino acid metabolism57(0.44%)ko0046019細胞周期Cell cycle55(0.42%)ko0411020抗壞血酸和新陳代謝Ascorbate and aldarate metabolism54(0.42%)ko00053

表3 12個與紅葉性狀相關的SNP位點

表4 篩選出的6條基因的功能注釋

表5 SNP的驗證結果

圖3 基于SNP標記的紅葉性狀的初步分子鑒定Fig.3 Preliminary molecular identification of red leaf traits by SNP markers

2.4 Sanger測序結果

基于重測序數(shù)據(jù)和SNP位置設計了6對引物，一代測序的驗證結果與重測序篩選結果，匹配率100%。共篩選出12個分子標記用于‘紅葉’杜仲紅葉性狀的早期鑒定(表5)。圖3給出了表5中序號1和2的SNP驗證結果。

3 討論

杜仲是單科單屬單種，‘紅葉’杜仲和‘小葉’杜仲均是在河南洛陽地區(qū)通過實生選育得到的無性系良種，遺傳背景相對較小。本研究在參考杜仲基因組大小(1.18 G)的基礎上，前期按10×的測序深度，后期在分析中設定嚴格的閾值，同時再進行部分sanger測序驗證，以期盡可能準確地得到SNP位點。根據(jù)重測序結果，在‘紅葉’杜仲中共檢測到5 407 823個SNP，583個位點能夠嚴重或者中度影響蛋白質結構，且與花色苷合成相關。目前也有一些學者采用同樣的方法得到類似的結論，如Ni等[16]對日本杏紅綠皮果實2個品種的全基因組重測序分析，對CDS區(qū)域6個序列數(shù)據(jù)集的比對分析發(fā)現(xiàn)181 331個SNPs，3 318個InDels和51 427個SVs，篩選出與紅色果皮顏色形成相關的13個候選基因。Xu等[17]對兩種栽培葡萄株系的幼苗葉片進行基因組重測序，篩選出635個sbc特異性基因和665個sbbm特異性基因，分析了非同義突變與控制顏色形成的基因的關系，發(fā)現(xiàn)NCED6表達與NCED1有關，并可能影響花青素積累。此次研究可能存在一些未覆蓋的位點，得到的數(shù)據(jù)準確性也可能有一定偏差，今后我們會加大群體和測序規(guī)模，繼續(xù)深入開展相關研究。

‘紅葉’杜仲和‘小葉’杜仲除了葉片顏色外葉形也存在一定差異，我們在分析中鎖定的是與花色苷合成相關的代謝通路，以盡量避免其他因素的干擾。目前多種研究方法共同揭示花色苷代謝途徑中SNP位點差異被廣泛使用。Jiao等[18]在14種中國野生葡萄品種的樣本中克隆了調(diào)節(jié)花色素苷生物合成的mybA的兩個轉錄因子VvmybA1和VlmybA2，在測序的片段中檢測到總共121個SNP，這些SNPs分布在啟動子區(qū)、內(nèi)含子區(qū)和編碼區(qū)。Liu等[19]利用MISA、GATK3和注釋工具，根據(jù)轉錄組數(shù)據(jù)從毛茛花瓣151,136個單基因中，鑒定和定位SSR和SNP標記。在8個花青素合成關鍵酶基因家族的116個單基因中，鑒定出127個SNP。ANS酶基因中沒有SNP位點，而DFR酶基因具有最大量的SNP位點。Oh等[20]利用全基因組重測序檢測SNP和轉錄組聯(lián)合分析黑米花青素的合成調(diào)控，發(fā)掘到9個候選基因中的3個基因與花色苷的合成有很好的相關性，尤其是基因Os06t0192100(UDP糖基轉移酶)與之前報道的Os04g0557500基因相似性很高。這樣聯(lián)合分析的方法為我們進一步深入研究‘紅葉’杜仲花青素的系統(tǒng)進化和基因通路提供了思路。在本研究結合‘小葉’杜仲的重測序結果與一代測序，發(fā)現(xiàn)12個純合差異位點是‘紅葉’杜仲特有的，更驗證了12個與‘紅葉’杜仲紅葉性狀相關的SNP位點。下一步將選擇葉色有差別的杜仲，通過PCR擴增測序的方法對上述SNP位點開展進一步驗證。

導致‘紅葉’杜仲葉片特殊呈色的原因可能會有很多，如基因變異、基因表達量的差異及表觀遺傳修飾的不同等，本研究僅從基因層面通過SNP位點差異初步探討‘紅葉’杜仲呈現(xiàn)紅色的可能原因，至于確切的原因，還需進一步的探索。