全基因組重測序分析4個毛竹變型的稈形和稈色變異

牟少華，李雪平，李 娟，高 健

(國際竹藤中心，國家林業和草原局/北京市共建竹藤科學與技術重點實驗室，北京 100102)

全基因組重測序是對已完成基因組序列的同一物種或者近緣種的其他個體進行全基因組測序，并在此基礎上對個體或群體進行差異性分析。新個體基因組序列通過映射參考基因組獲得同源數據集，然后檢測基因組范圍內的單核苷酸多態性(Single nucleotide polymorphism，SNP)、插入缺失(Insertion/deletion，InDel)及結構變異(Structural variation，SV)等多態位點[1]。運用高通量的重測序技術進行遺傳變異分析已經在谷子[2]、大豆[3]、水稻[4]、花生[5]、鼠尾草[6]、葡萄[7]、木蘭[8]、煙草[9]、竹黃菌[10]和蘑菇[11]等植物上得到了廣泛應用。

毛竹(Phyllostachysedulis)是我國特有的、利用歷史悠久的重要經濟竹種，也是栽培面積最廣的竹種，其生長速度快，材性優良，繁殖能力強，一次造林可以永續利用。毛竹筍是一種全天然營養食品，含有豐富的蛋白質、脂肪、糖類、胡蘿卜素、維生素、鈣、磷、鐵以及18種氨基酸等成分。因此，毛竹是優良的筍材兩用竹種，有著巨大的開發和應用潛力[12]。毛竹廣泛分布于我國南方16個省區，跨越北、中、南3種亞熱帶氣候，在長期適應不同的生長環境和栽培經營措施下，加上自身發生遺傳漂移等因素，發生了一系列的形態變異。目前，在天然毛竹林中已經發現21個毛竹變型[13]，這些變型在形態上表現出豐富的多態性，尤其是稈形、稈色等形態性狀差異顯著。為研究這些形態差異究竟是因為基因變異而導致的可遺傳變異還是由于生境不同而產生的表型可塑性，非常有必要在基因水平上揭示毛竹的遺傳多樣性和變異程度。

近年來，隨著分子生物學技術的突飛猛進以及組學技術在竹類植物上的應用，特別是毛竹全基因組測序工作的完成[14]，竹類植物基因組研究進入后基因組時代。目前，毛竹已得到一個完整的可變剪接圖譜[15]，對P450s[16]、AP2/ERF[17]、NAC[18]、SAUR[19]、AQP[20]、IQD[21]、SBP-like[22]、WRKY[23]、HD-Zip[24]和TCP[25]、Hsp[26]、CO-Like[27]等基因家族已采用生信方法進行了鑒定和功能驗證，但尚未發現有關毛竹變型基因組序列變異的研究。本研究采用成熟且經濟的二代測序技術，以毛竹全基因組為參考，對4個有代表性的毛竹變型進行重測序，對其SNP、Indel、結構變異SV等進行深度挖掘，從而揭示其全基因組的突變類型，分析其廣適性狀的變異基因，為從全基因組尺度分析毛竹的“種性”和變型種質“品性”獲得海量基因組序列數據。

1 材料和方法

1.1 試驗材料

試驗材料為毛竹的4個變型(表1)：圣音竹、黃槽毛竹、黃皮花毛竹和厚壁毛竹，編號分別為R01、R02、R03和R04，均采自國際竹藤中心安徽太平試驗中心。選取剛長出且尚未展開的嫩葉，液氮速凍后超低溫保存。

表1 4個毛竹變型形態特征Tab.1 Morphological characteristics of 4 variations of Moso bamboo

1.2 試驗方法

1.2.1 DNA提取與基因組測序 基因組DNA采用改良CTAB法提取[28]。將檢測合格的樣品基因組DNA，用超聲波方法將DNA進行片段化，然后進行DNA片段的純化、末端修復、在3′端加A、連測序接頭，進行瓊脂糖凝膠電泳完成片段大小選擇，通過PCR擴增建成測序文庫，進行文庫的質檢，合格文庫采用Illunima Hiseq 2500平臺進行測序，得到原始測序序列(Raw reads)，對其進行數據過濾，去掉帶接頭的reads，過濾N含量高(超過10%)和質量值低(低于10的堿基數超過50%)的reads，得到Clean reads，用作后續的信息分析。

1.2.2 與參考基因組進行比對統計 將重測序得到的測序reads重新定位到毛竹參考基因組上，進行后續的變異分析[14]。運用Bwa軟件將測序獲得的短序列與毛竹的參考基因組(毛竹基因組數據庫 BambooGDB下載網址：http://www.bamboogdb.org/page/download.jsp)進行比對，通過Samtools flagstat命令對Clean reads比對定位到參考基因組上的位置，然后統計4個樣品的測序深度和基因組覆蓋度等方面的信息，進行序列變異的檢測。

1.2.3 檢測基因組變異 使用GATK軟件工具包對樣品的SNP和InDel進行檢測[29]。根據參考基因組上Clean reads的定位結果，使用Picard軟件去重復、GATK軟件進行局部重比對和堿基質量值的校正等預處理，然后進行SNP和InDel的變異檢測。對SNP進行嚴格的過濾，包括：SNP cluster過濾(5 bp內若有2個SNP過濾掉)，InDel附近5 bp內的SNP過濾掉；相鄰2個InDel距離小于10 bp過濾掉，從而獲得SNP位點集。使用BreakDancer軟件檢測SV，主要是基于序列的比對結果，獲得測序文庫的插入片段大小、方差，再查找序列和參考基因組間比對的異常結果，尋找可能的結構變異。

1.2.4 注釋SNP、InDel和SV 通過SnpEff軟件進行SNP、InDel和SV的注釋[30]。根據檢測獲得的變異位點比對到參考基因組的位置基因、CDS位置等，獲得變異位點在基因組的發生區域(基因區、基因間區或CDS區等)，以及產生的同義非同義突變等。

1.2.5 挖掘變異基因和功能注釋 通過尋找樣品間與參考基因組發生非同義突變SNP基因、CDS區發生的InDel基因與SV基因，挖掘可能存在的功能變異基因。通過Blast將變異基因與NR[31]、GO[32]、COG[33]、KEGG[34]等功能數據庫進行比對，獲得基因的注釋，以便分析基因的功能。

2 結果與分析

2.1 與參考基因組比對統計

4個竹種通過高通量測序(Illunima Hiseq 2500等測序平臺)得到原始圖像數據文件，經堿基識別分析轉化為原始測序序列(Raw reads)，然后去除帶接頭的reads，過濾得到Clean reads，分別為180 506 436，184 142 480，213 567 146，188 998 942 bp，定位到毛竹參考基因組的占所有Clean reads數的百分比分別為96.13%，96.32%，96.33%和97.44%，雙端均定位到毛竹參考基因組上并且距離符合測序片段的長度分布占所有Clean reads數的百分比分別為88.18%，88.42%，88.57%和89.65%。樣品平均覆蓋深度在10×以上，各深度對應的基因組覆蓋比例如表2所示。

表2 4個毛竹變型的測序數據統計Tab.2 Summary of sequencing data of 4 variants of Moso bamboo

2.2 SNP的檢測與注釋

2.2.1 SNP檢測 對4個毛竹變型的基因組DNA進行SNP檢測，獲得SNP位點集(表3)，4個毛竹變型樣品分別檢測到的SNP數量為1 479 567～1 616 020，SNP總數為2 035 983，Ti/Tv值約為2.90。樣品之間的SNP數統計結果見表4。

表3 4個毛竹變型SNP數據統計Tab.3 Summary of SNPs detected in 4 variants of Moso bamboo

表4 毛竹變型間的SNP統計Tab.4 Summary of SNPs detected between variants of Moso bamboo

2.2.2 SNP注釋 運用SnpEff軟件對4個毛竹變型全基因組SNP進行注釋，得到其變異位點在基因組發生的區域或類型(表5)。圣音竹發生在CDS區域內的SNP共計24 426個，其中，同義突變占比為47.28%，非同義突變占比為51.56%。黃槽毛竹發生在CDS區域內的SNP共計25 319個，其中，同義突變占47.26%，非同義突變占51.50%。黃皮花毛竹發生在CDS區域內的SNP共計25 628個，其中，同義突變占47.26%，非同義突變占51.52%。厚壁毛竹發生在CDS區域內的SNP共計24 738個，其中，同義突變占47.49%，非同義突變占51.31%。

表5 4個毛竹變型SNP注釋結果Tab.5 SNP annotations in 4 variants of Moso bamboo

2.3 InDel檢測與注釋

2.3.1 InDel檢測 從4個樣品InDel檢測結果(表6)可以看出，圣音竹全基因組范圍與編碼區分別檢測出174 657，3 173個InDel，CDS區有插入突變2 296個和缺失突變877個，其中純合突變為578個。黃槽毛竹全基因組范圍與CDS區分別檢測出185 096，3 262個InDel，CDS區有插入突變2 323個和缺失突變939個，其中629個為純合突變。黃皮花毛竹全基因組范圍與CDS區分別檢測出195 518，3 296個InDel，CDS區有插入突變2 371個和缺失突變925個，其中純合突變為639個。厚壁毛竹全基因組范圍與CDS區分別檢測出176 393，3 121個InDel，CDS區有插入突變2 246個和缺失突變875個，其中純合突變為611個。對樣品在全基因組范圍和CDS區的InDel長度進行統計(圖1)，基因組范圍存在較多的+1、-1、+2、-2類型突變，在CDS區域存在較多的+1、-1、+3、-3類型突變。

表6 4個毛竹變型全基因組和編碼區InDel統計Tab.6 Summary of InDels detected in 4 variants of Moso bamboo

橫坐標小于0為缺失，大于0為插入。

2.3.2 InDel注釋 根據4個樣品與參考基因組的InDel檢測結果，進行兩兩比較，獲得樣品間的InDel變異位點(表7)。黃皮花毛竹與厚壁毛竹間的InDel數最少，為45 967；圣音竹與黃槽毛竹間InDel數最多，為48 551。

表7 毛竹變型間的InDel統計Tab.7 Summary of InDels detected between variants of Moso bambooo

2.4 SV檢測與注釋

2.4.1 SV檢測 使用BreakDancer進行SV的檢測(表8)，圣音竹共檢測到148 682個SV，其中插入3 701個，占總變異2.49%；缺失31 439個，占總變異21.15%；染色體間易位107 945個，占總變異72.60%；染色體內部易位4 991個，反轉321個，其他285個。黃槽毛竹共檢測到176 960個SV，其中插入17 377個，占總變異9.82%；缺失36 496個，占總變異20.62%；染色體間易位116 941個，占總變異66.08%。黃皮花毛竹共檢測到181 687個SV，其中插入9 498個，占總變異5.23%；缺失37 880個，占總變異20.85%；染色體間易位127 646個，占總變異70.26%。厚壁毛竹共檢測到148 968個SV，其中插入2 481個，占總變異1.67%；缺失30 971個，占總變異20.79%；染色體間易位110 054個，占總變異73.88%。

表8 4個毛竹變型結構變異數量統計Tab.8 Summary of SVs detected in 4 variants of Moso bamboo

2.4.2 SV注釋 對結構變異深入分析表明(表9)，圣音竹外顯子區變異總數為1 894個，內含子區變異總數為789個，基因間區變異總數為32 778。黃槽毛竹外顯子區變異總數為2 732個，內含子區變異總數為1 209個，基因間區變異總數為50 241個。黃皮花毛竹外顯子區存在2 498個變異，內含子區存在1 064個變異，基因間區存在44 150個變異。厚壁毛竹外顯子區分別存在1 731，內含子區存在773個，基因間區存在31 237個。

表9 4個毛竹變型結構變異注釋結果Tab.9 SV annotations in 4 variants of Moso bamboo

2.5 DNA水平的變異基因挖掘與功能注釋

2.5.1 變異基因挖掘 分別提取毛竹4個變型與參考基因組間發生非同義突變的SNP以及CDS區發生InDel與SV的基因，各樣品的變異如表10所示。稈形變異的2個竹種圣音竹和厚壁毛竹基因組分別存在9 307，9 148個基因變異，其中SNP突變基因分別為4 973，4 981個，InDel基因分別為2 735，2 699個，SV突變的基因分別為1 599,1 468個，多個基因同時存在多種類型突變。稈色變異的2個竹種黃槽毛竹和黃皮花毛竹基因組分別存在10 194,9 977個基因變異，其中，SNP突變基因分別為5 094,5 095個，InDel基因分別為2 799,2 816個，SV突變基因2 301,2 066個，且多個基因同時存在多種類型突變。

表10 4個毛竹變型的變異基因統計Tab.10 Summary of gene variations in 4 variants

2.5.2 變異基因的功能注釋 運用Blast軟件，將挖掘的變異基因分別與GO、KEGG和COG等功能數據庫進行比對，獲得這些變異基因的注釋，進而分析基因功能。4個毛竹變型注釋到數據庫的變異基因數分別為8 351，9 195，8 966，8 197個。

通過GO數據庫，獲得變異基因GO分類統計結果(圖2)。從圖中可以看出，細胞組件、分子功能和生物過程3個基因功能分類體系中，GO各分類內容對應的基因數以及基因數目所占百分比。對GO數據進一步分析發現，在細胞組分方面，葉綠素合成相關基因有1 848個，細胞壁相關基因750個，細胞骨架相關基因194個；在生物過程方面，參與類胡蘿卜素合成過程的基因有60個，參與花青素合成過程中的調控以及紫外光下組織中花青素積累的相關基因有55個；在功能基因方面，參與信號轉導的基因283個，轉錄因子475個，激素相關基因75個。

圖2 變異基因的GO注釋分類圖Fig.2 Classification of gene variations compared with GO database by Blast

COG數據庫是對基因產物進行直系同源分類的數據庫。從變異基因COG分類統計結果(圖3)中可以看出，在不同的功能分類中，基因數目相差很大，其中參與功能預測、轉錄、復制重組修復和信號轉導機制的基因數目較多。基因所占比例多少，反映對應時期和環境下代謝或者生理偏向等內容。COG注釋獲得的參與功能預測的基因1 131個，參與糖代謝途徑的基因355個，細胞壁、細胞膜合成相關的基因有276個，信號轉導機制的基因數為228個，細胞周期調控、細胞分裂和染色的相關基因145個。

通過KEGG數據庫，變異基因注釋到114個KEGG代謝通路中，分析共獲得了1 310個轉錄因子及813個信號轉導相關單基因；在功能基因方面，共獲得618個與細胞壁構建相關的代謝基因。KEGG數據庫能夠系統地分析基因產物在細胞中的代謝途徑和功能，有助于把基因及表達信息作為一個整體的網絡進行研究。從變異基因的糖代謝通路(圖4)可以看出，整個通路有46種不同的酶通過復雜的生化反應形成，框內的數字分別代表酶的號碼。其中的4種酶：1.2.4.1(丙酮酸脫氫酶)、1.8.1.4(二氫硫辛酸脫氫酶)、5.3.1.1(磷酸丙糖異構酶)和5.4.2.2(葡萄糖磷酸變位酶)對應的框代表變異基因中與此通路相關。

圖3 變異基因的COG注釋分類圖Fig.3 Classification of gene variations compared with COG database

將變異基因與GO、COG、KEGG等3個功能數據庫進行比對、注釋，共獲得1 739個參與細胞壁和細胞膜合成的基因、501個細胞骨架構建相關基因、1 785個轉錄因子、1 324個信號轉導、104個赤霉素及76個激素代謝相關基因、1 938個葉綠素合成相關的基因、73個類胡蘿卜素合成代謝相關基因和68個花青素合成調控基因。

3 結論與討論

3.1 4個毛竹變型的稈形、稈色變異表征

經過長期的自然演化和人工栽培后，毛竹種內產生豐富的遺傳變異，從而形成不同的結構形態。圣音竹、黃皮花毛竹、黃槽毛竹和厚竹是毛竹種的4個變型，具有明顯的稈形或稈色特征。圣音竹竹稈從上面向基部逐漸膨大呈喇叭狀，同時節間縮短，且上粗下細或下部緊靠節處有不規則凹陷，節部呈波浪狀歪斜。厚竹稈略呈四方形或橢圓形，壁厚，上部近實心，分枝以下節隆起較高。黃皮花毛竹稈黃色或黃綠色基本各半，有寬窄不等的綠色縱條紋。黃槽毛竹稈主要為綠色，但節間分枝一側縱溝槽為黃色或淡黃色。這些在竹稈形態和顏色等方面的不同變異使其具有很高的觀賞價值。

3.2 4個毛竹變型的基因組變異分析

通過對4個毛竹變型進行重測序，初步統計分析了這些形態變異發生的分子數據。樣品檢測到的SNP總數為2 035 983，Ti/Tv值約為2.90，說明同種類型的堿基之間突變比不同類型間更容易發生。分別提取毛竹4個變型與參考基因組間發生非同義突變的SNP、CDS區發生InDel與SV的基因，可以發現其基因組基因的變異數目和變異類型呈現一定的規律，其中，2個稈形變異竹種圣音竹和厚壁毛竹基因組基因的變異數目和變異類型相近，而2個稈色變異竹種黃槽毛竹和黃皮花毛竹基因組的變異數目和變異類型相近，但稈形和稈色變異的竹種兩兩之間差異較大。

3.3 變異基因的功能注釋與變異性狀的關聯分析

DNA水平的變異基因挖掘和功能注釋，可以分析基因產物在細胞中的代謝途徑以及這些基因產物的功能。將變異基因通過多個數據庫注釋，共獲得1 739個基因參與細胞壁、細胞膜合成，501個細胞骨架構建相關基因，以及1 785個轉錄因子、1 324個信號轉導、104個赤霉素和76個激素代謝相關基因。參考矢竹的研究結果，其地下莖生長發育受到從各類激素、轉錄因子到其下游功能基因等組成的復雜分子網絡的協同調控，細胞壁、細胞骨架下游功能相關基因在矢竹地下莖節間快速生長階段表達最高，而赤霉素合成及信號傳導相關基因則在節間伸長起始階段最高[35]。因此，進一步對毛竹變型細胞組分、生物過程和分子功能相關變異基因的研究，有利于推測毛竹稈形變異的分子調控機制。

竹類植物中的色素主要有葉綠素、類胡蘿卜素和花青素3大類，色素的變化會影響植株的呈色變異。通過數據庫注釋，發現的葉綠素合成相關基因、類胡蘿卜素合成代謝相關基因和花青素合成調控基因，深入研究這些差異基因的調控途徑利于從DNA水平上解釋稈色的變異。研究表明[36]，毛竹不同變異類型在凈光合速率、葉綠素含量、β胡蘿卜素含量、硝酸還原酶活性、游離氨基酸含量這5個生理指標中存在著顯著性差異，其中，黃皮花毛竹最高，其次為厚壁毛竹，然后是圣音毛竹、黃槽毛竹。對毛竹的栽培變種進行ISSR和AFLP分子標記分析，分子遺傳的分析結果與阮曉賽[37]進行的形態學研究結果相似，變型之間的遺傳變異程度較小。進一步分析這4個毛竹變型基因產物的功能及其代謝途徑，也有利于從基因水平上對其生物學和生理特性加以解釋和分析。數據的挖掘與深入分析工作正在進行，后續分析將解析毛竹眾多種質的分子遺傳基礎，闡析不同變異類型的分子品性、基因家族、功能基因等遺傳根基。