四個竹稈變異毛竹變型的全基因組序列分析

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

( 國際竹藤中心 竹藤科學與技術重點實驗室， 北京 100102 )

毛竹 () 是我國獨有的傳統經濟竹種(江澤慧，2002)，分布范圍較廣，現有毛竹變型20多種(馬乃訓等，2014)。毛竹變型在形態上表現出豐富的多態性，尤其是稈色性狀方面差異表現顯著，例如花毛竹稈為黃色，有寬窄不等的綠色縱條紋；而綠皮花毛竹稈為綠色，但節間有淡黃色細縱條紋。稈色的變異大大豐富了園林觀賞種類，提高了園林觀賞價值。毛竹竹稈性狀的遺傳變異是遺傳育種工作關注的重點。

目前，高通量測序技術可以分析一個物種的基因組的全貌，已經在谷子 (Bai et al., 2013; 賈小平等，2019)、水稻 (Takagi et al., 2013)、大豆 (Qi et al., 2014; Zhou et al., 2015; 張彥威等，2016)、黃秋葵(張少平等，2017)、厚樸(尹彥棚等，2020)、菜豆 (Jeremy et al., 2010)和番茄 (Lin et al., 2014)等植物上得到了廣泛應用。

近年來，隨著分子生物學和組學技術的發展，毛竹全基因組序列獲得公布 (Peng et al., 2013)，一些與毛竹性狀相關的基因家族，如2(Wu et al., 2015)、( Bai et al., 2016)、(Sun et al., 2016)、-like(Pan et al., 2016)、-(Chen et al., 2017)、(Xie et al., 2019)、-(Liu et al., 2016)等已進行鑒定和功能驗證，但是在基因組層面上的研究較少，特別是毛竹變型稈色相關的研究薄弱，僅對黃槽毛竹和黃皮花毛竹兩個毛竹變型基因組序列變異進行初步探索 (牟少華等，2020)，這在一定程度上限制了毛竹遺傳育種的應用發展。從基因水平上揭示毛竹的變異程度，是分析毛竹變型形態差異產生原因的重要手段之一。因此，開展毛竹變型基因組研究，揭示毛竹變型全基因組突變類型，探究黃酮類和硝酸還原酶等代謝途徑相關基因，對解析毛竹豐富的遺傳多樣性以及性狀相關的遺傳變異具有重要意義。

本研究以黃皮毛竹等4個有代表性的竹稈變異毛竹變型為研究對象，毛竹全基因組作為參考基因組，采用高通量測序技術，構建全基因組數據庫，并利用生物信息學的方法對獲得的核酸序列組裝，檢測并注釋其單核苷酸多態性(single nuclotide polymorphisms，SNP)、結構變異(structural variation，SV)和小片段插入缺失(insertion and deletion，InDel)等，注釋變異基因功能，積累基因組序列數據，以便為從全基因組水平上深入地分析毛竹的遺傳變異，為遺傳育種提供遺傳基礎。

1 材料與方法

1.1 材料

供試材料選自國際竹藤中心安徽太平試驗中心種質資源圃。選取4個毛竹變型(表1)適量新鮮幼嫩的葉片，經液氮罐中速凍后，放入-80 ℃冰箱凍存備用。

表 1 四個毛竹變型樣品簡表Table 1 Brief introduction of four variant samples of Moso bamboo

1.2 方法

1.2.1 基因組測序 毛竹變型葉片DNA提取(Zidani et al., 2005)后，經過打斷、損傷修復及連接接頭、PCR富集、文庫質量檢測，建成測序文庫，在Illunima Hiseq 2500測序平臺上運行獲得原始數據，將數據過濾后得到高質量數據。

1.2.2 比對統計 使用BWA軟件 (Li & Durbin, 2009) 將測序數據比對定位到已測序的毛竹基因組的位置，統計測序深度和基因組覆蓋度等信息。

1.2.3 檢測SNP、InDel和SV 使用Picard軟件 (Gordon et al., 2012) 去重復和GATK軟件 (Mckenna et al., 2010) 預處理后，檢測SNP和InDel變異。使用BreakDancer軟件 (Chen et al., 2009) 檢測SV變異，具體方法參照牟少華等(2020)。

1.2.4 注釋SNP、InDel和SV 運用SnpEff軟件 (Cingolani et al., 2012)注釋SNP、InDel和SV，具體方法參照牟少華等(2020)。

1.2.5 注釋功能基因 運用BLAST軟件，對篩選得到的功能可能變異基因的基因序列與GO (Ashburner et al., 2000)、COG (Tatusov et al., 2000)和KEGG (Minoru et al., 2004)三大功能數據庫，進行BLAST比對，得到基因注釋。

2 結果與分析

2.1 與毛竹基因組比對

4個竹種通過高通量測序得到測序數據。金絲毛竹樣品(R02)過濾后的Clean Reads最少，為82 276 884 bp；花毛竹樣品 (R04) 的Clean Reads最多，為112 054 728 bp。定位到毛竹參考基因組的占所有Clean Reads數的百分比在99.45%以上，雙端均定位到毛竹參考基因組上并且距離符合測序片段的長度分布的占所有Clean Reads數的百分比在88%左右，說明參考基因組選擇合適，且相關實驗過程不存在污染，測序Reads的比對率會高于70%。另外，比對率的4個毛竹變型與毛竹參考基因組親緣關系較近、基因組組裝質量高，而且Reads測序質量高。4個樣品平均覆蓋深度均在10×左右(表2)。

表 2 四個樣品數據產出統計表Table 2 Output statistics among four samples

2.2 SNP的檢測與注釋

2.2.1 SNP檢測 4個毛竹樣品檢測后獲得SNP位點統計表(表3)。其中，花毛竹樣品(R04)的SNP數量最多，為1 691 715；綠皮花毛竹樣品(R03)的SNP數量最少，為1 534 648。4個樣品中，轉換類型(transition, Ti) SNP數量與顛換類型(transversion, Tv) SNP數量的比值Ti/Tv在3.05～3.10之間，說明這些毛竹變型轉換比顛換更容易發生。雜合類型 (heterozygosity, Het) SNP數量為純合類型 (homozygosity, Homo) SNP數量的10倍左右，雜合比率為88.53%～92.01%。其中，花毛竹樣品(R04)雜合比率最高，為92.01%，說明其雜合程度最高。綠皮花毛竹樣品(R03)雜合比率最低，為88.53%。

表 3 四個樣品SNP位點統計表Table 3 SNP loci statistics in four samples

根據4個毛竹樣品與參考基因組的比對結果，匯總樣品間SNP的統計結果見表4，表中各數值為對應的橫縱兩樣品之間的SNP數。從表中可以看出，金絲毛竹(R02)與綠皮花毛竹樣品(R03)間的SNP數最多。

表 4 四個樣品間的SNP統計表Table 4 Summary of SNPs detected between four samples

2.2.2 SNP注釋 對4個樣品SNP進行注釋，獲得其變異位點發生的區域或類型(圖1)。4個毛竹變型發生在編碼區(coding sequence，CDS)區域內的SNP數量占比均為2%左右，其中同義突變占48%左右，非同義突變占51%左右。非同義突變率與同義突變率的比值大于1，預示著有正向選擇效應。

圖 1 黃皮毛竹(R01)的SNP注釋圖Fig. 1 SNP annotations pie of Phyllostachys edulis f. holochrysa (R01)

2.3 InDel檢測與注釋

2.3.1 InDel檢測 對4個毛竹變型InDel進行統計(表5)，可以發現4個樣品全基因組范圍檢測出的InDel總數范圍為271 648～292 253，其中插入類型的突變總數略低于缺失突變總數；編碼區檢測出的InDel總數為4 711～4 877，其中插入突變總數為缺失突變的67%左右。各樣品中，全基因組范圍內純合突變數約為雜合突變數的2倍，編碼區純合突變數略低于雜合突變數。

表 5 四個樣品InDel統計表Table 5 Summary of InDels detected in four samples

對4個樣品各區域的InDel長度進行統計發現，編碼區存在較多的+1、-1、+ 3、-3 類型突變，而基因組范圍存在較多的+1、-1、+ 2、-2 類型突變。其中，數值代表InDel的長度(10 bp以內)；大于0為插入；小于0為缺失。

將4個樣品的InDel進行兩兩比較，統計結果見表6。表中各數值為對應的橫縱兩樣品之間的InDel數。

表 6 毛竹變型間的InDel統計表Table 6 Summary of InDels detected between variants of Moso bamboo

2.3.2 InDel注釋 對比毛竹參考基因組的基因、CDS位置等信息，注釋各樣品InDel位點的發生位置以及是否為移碼突變等，具體注釋結果如圖2所示。4個毛竹變型發生在編碼區的InDel數量均在1.7%左右。移碼突變的InDel有可能會引起基因功能的改變。

圖 2 黃皮毛竹(R01) InDel注釋圖Fig. 2 InDel annotations pie of Phyllostachys edulis f. holochrysa (R01)

2.4 SV檢測與注釋

2.4.1 SV檢測 檢測4個樣品與參考基因組間的插入(insertion，INS)、缺失(delection，DEL)、反轉(inversion，INV)、染色體內部易位(intra-chromosome translocation，ITX)、染色體間易位(inter-chromosomal translocation，CTX)，得到的各類型SV數量統計見表7。其中，4個毛竹變型都表現為缺失類型的SV數量最多，其次為染色體內易位類型。

表 7 四個樣品SV數量統計Table 7 Quantity statistics of SVs detected in four samples

2.4.2 SV注釋 檢測4個樣品SV發生位置信息，并對缺失、插入和反轉3種類型的結構變異進行注釋。結果表明，4個毛竹變型在各區域分布的SV總體情況一致，注釋到的變異基因數目以基因間區的缺失類型最多，其次為基因間區的插入類型(表8)。

表 8 四個毛竹變型SV注釋結果統計表Table 8 SV annotations in four variants of Moso bamboo

2.5 變異基因功能注釋與分析

2.5.1變異基因挖掘 分別統計4個樣品的非同義突變的SNP以及CDS區發生InDel和SV的基因(表9)，尋找可能存在功能變異的基因。在4個毛竹樣品中，花毛竹(R04)基因組存在12 555個基因變異，其中，非同義突變SNP基因為5 563個，InDel基因為4 006個，SV突變的基因為2 986個，差異基因總數和SV突變基因數最多。在3類變異基因中，非同義突變SNP基因數最多，InDel基因數量次之，SV突變的基因最少。

表 9 四個毛竹變型的變異基因統計表Table 9 Summary of mutant genes in four variants of Moso bamboo

2.5.2 變異基因的功能注釋 黃皮毛竹、金絲毛竹、綠皮花毛竹和花毛竹注釋到數據庫中的變異基因數分別為7 575、7 538、7 476和7 728。變異基因GO分類統計結果圖(圖3)中，顯示出在3大基因功能分類體系(分子功能、細胞組件和生物過程)的56個分類內容中所對應的基因數和基因占比。其中，在細胞組件分類中，與葉綠素合成相關的基因有2 431個；在生物過程分類中，參與類胡蘿卜素合成過程的基因有75個；在分子功能分類中，與花青素合成調控以及紫外光下組織中花青素積累的相關基因有80個。4個毛竹變型在相應的基因功能中的變異基因數目有差異，例如：花毛竹有21個與類胡蘿卜素合成相關的基因；綠皮花毛竹有17個相關基因；而黃皮毛竹有18個相關基因，基因數量和種類的差異，可能引起相應的功能變化。深入研究葉綠素、類胡蘿卜素和花青素合成相關基因以及這些差異基因的調控途徑，有利于從DNA水平上解釋稈色的變異。

橫坐標為GO的分類內容，縱坐標的左側為基因數占比，右側為基因數量。 1. 代謝過程; 2. 細胞過程; 3. 對刺激的反應; 4. 生物調節; 5. 定位; 6. 定位確立; 7. 細胞組成或生物形成; 8. 分化過程; 9. 多細胞生物過程; 10. 繁殖; 11. 生殖過程; 12. 信號; 13. 多組織過程; 14. 生長; 15. 免疫系統過程; 16. 死亡; 17. 細胞增殖; 18. 生物粘附; 19. 節律過程; 20. 病毒繁殖; 21. 色素沉著; 22. 運動; 23. 細胞死亡; 24. 碳利用; 25. 細胞部分; 26. 細胞; 27. 細胞器; 28. 膜; 29. 細胞器部分; 30. 膜部分; 31. 高分子復合物; 32. 細胞外區域; 33. 膜內腔; 34. 細胞連接; 35. 細胞外基質; 36. 類核; 37. 病毒粒子; 38. 細胞外基質部分; 39. 細胞外區部分; 40. 病毒粒子部分; 41. 綁定; 42. 催化活性; 43. 運輸活動; 44. 核酸結合轉錄因子活性; 45. 結構分子活性; 46. 電子載體活性; 47. 酶調節活性; 48. 活動分子傳感器; 49. 抗氧化活性; 50. 受體活性; 51. 蛋白結合轉錄因子活性; 52. 營養庫活性; 53. 翻譯調節活性; 54. 金屬伴侶活性; 55. 蛋白質標記; 56. 通道調節活性。Abscissa is GO classification, the left side of the ordinate is the percentage of genes, and the right side is the number of genes. 1. Metabolic process; 2. Celluar process; 3. Response to stimulus; 4. Biological regulation; 5. Localization; 6. Establishment of localization; 7. Cellular component organization or biogenesis; 8. Developmental process; 9. Multicellular organismal process; 10. Reproduction; 11. Reproductive process; 12. Signaling; 13. Multi-organism process; 14. Growth; 15. Immune system process; 16. Death; 17. Cell proliferation; 18. Biological adhesion; 19. Rhythmic process; 20. Viral reproduction; 21. Pigmentation; 22. Locomotion; 23. Cell killing; 24. Carbon utilization; 25. Cell part; 26. Cell; 27. Organelle; 28. Membrane; 29. Organelle part; 30. Membrane part; 31. Macromolecular complex; 32. Extracellular region; 33. Membrane-enclosed lumen; 34. Cell junction; 35. Extracellular matrix; 36. Nucleoid; 37. Virion; 38. Extracellular matrix part; 39. Extracellular region part; 40. Virion part; 41. Binding; 42. Catalytic activity; 43. Transporter activity; 44. Nucleic acid binding transcription factor activity; 45. Structural molecule activity; 46. Electron carrier activity; 47. Enzyme regulator activity; 48. Molecular transducer activity; 49. Antioxidant activity; 50. Receptor activity; 51. Protein binding transcription factor activity; 52. Nutrient reservoir activity; 53. Translation regulator activity; 54. Metallochaperone activity; 55. protein tag; 56. Channel regulator activity.圖 3 黃皮毛竹(R01)變異基因的GO注釋分類圖Fig. 3 Classification of Phyllostachys edulis f. holochrysa (R01) mutant genes compared with GO database by BLAST

變異基因COG注釋分類圖(圖4)直觀顯示出COG功能分類條目上分別對應的頻率，其中涉及到功能注釋、轉錄、復制重組修復和信號轉導機制的對應數值高。獲得的功能注釋基因為1 630個，參與復制、重組和修復的基因數為369個，信號轉導機制的基因數為291個，轉錄的相關基因為222個。

圖 4 黃皮毛竹(R01)變異基因的COG注釋分類圖Fig. 4 Classification of Phyllostachys edulis f. holochrysa (R01) mutant genes compared with COG database

KEGG數據庫系統地分析4個毛竹變型的基因產物在生物學過程中的功能。以黃皮毛竹(R01)的纈氨酸、亮氨酸和異亮氨酸生物合成通路為例(圖5)，注釋到57個基因參與該通路，其中23個變異基因。整個通路涉及不同的酶連接一系列生化反應形成，其中，框內的數字代表enzyme的號碼，紅色的框代表通路相關變異基因。

3 討論

全基因組重測序可以在已知植物的基因組序列基礎上，對其不同品種的基因組序列進行測序，從而找出個體與該物種間的差異性(Ley et al., 2008)。隨著毛竹全基因組序列的公開發表(Peng et al., 2013)，研究毛竹不同變種或變型基因組序列差異成為可能。全基因組重測序可以檢測個體的全部基因組序列，掃描出一些與該個體生長性狀密切相關的變異位點(宋志芳等,2017)。毛竹的地下莖中單軸散生，其不同變異類型也都是散生竹。在毛竹自身的遺傳基因漂移，以及長期的栽培措施和自然環境變遷等因素的影響下，毛竹種內產生了很多遺傳變異，產生各種獨特的結構形態，表現出豐富的園林觀賞性狀。其中，黃皮毛竹、花毛竹、綠皮花毛竹在竹稈顏色方面表現出不同程度的變異，使其具有更高的園林觀賞價值。

圖 5 黃皮毛竹(R01)變異基因的KEGG通路代謝圖Fig. 5 Pathway of Phyllostachys edulis f. holochrysa (R01) mutant genes compared with KEGG database by BLAST

對4個毛竹變異類型全基因組重測序，初步統計分析了其基因組數據，與毛竹參考基因組進行比對，檢測其SNP、InDel和SV。SNP類型的變異分為轉換和顛換兩種，4個毛竹樣品的轉換/顛換(Ti/Tv)的比值均約等于3，說明轉換類型比顛換類型更容易發生。SNP雜合比例約為90%，說明樣品有很高的雜合度，即同源染色體上SNP位點含不同類型的堿基比例高。InDel位點數同樣能反映不同樣品與毛竹基因組之間的差異，并且編碼區的InDel會引起移碼突變，影響基因功能。SV中缺失、插入、反轉、易位4種類型的數量，反映出基因組水平上大片段的缺失、插入、倒置、易位等序列差異。通過生物信息學分析，比較不同稈色的變異類型在全基因組水平上的結構差異，并進行差異注釋，從而為毛竹選育提供遺傳基礎，也為重要基因的功能研究提供有利依據。

顏色變異是植物中較常見的表型變異，其中關于水稻、擬南芥、菊花等多種植物均有葉色變異的報道。據不完全統計，水稻葉綠體含量基因超過140個(趙紹路等，2018)。竹子中的色素分為3大類：葉綠素、花青素和類胡蘿卜素。通過功能數據庫比對，對4個毛竹變型的變異基因，進行基因功能注釋和分析。GO數據庫注釋聚類反映了毛竹變型在不同功能組分類中基因數目和基因產物的屬性，其中與稈色變異有關的葉綠素、類胡蘿卜素和花青素等色素合成相關基因作為重點關注對象進行分析。COG數據庫注釋了基因產物的直系同源分類，不同分類對應的基因數目差別很大，反映了不同條件下的生理或者代謝偏好等。 KEGG數據庫將基因和多種酶形成通路，有氨基酸生物合成、類胡蘿卜素生物合成、類黃酮生物合成、 萜類化合物的生物合成、植物激素信號轉導、參與卟啉和葉綠素代謝等顯著富集。其中，葉綠體、類胡蘿卜素和花青素等色素合成相關通路，是與稈顏色相關的主要代謝通路。

結合不同稈色毛竹變型的生物學和生理學特性，研究全基因組序列色素合成相關基因，有助于從基因水平上解析其稈色變異原因。有研究表明，毛竹不同變異類型在葉綠素含量、β胡蘿卜素含量等生理指標中存在顯著性差異(陳建華等，2011)，株型較大的花毛竹生理指標值比較小型的龜甲竹、綠槽毛竹大(晏育存，2011)。毛竹變型ISSR和AFLP分子標記分析表明，變型間的遺傳變異程度較小(阮曉賽，2008)。在參考黃槽毛竹和黃皮花毛竹兩個稈色變異毛竹變型的研究結果(牟少華等，2020)基礎上，通過對黃皮毛竹等4個毛竹變異類型重測序， 進行DNA水平的變異基因功能注釋，可以分析基因產物在細胞中的代謝途徑及功能，尤其是對黃酮類、類胡蘿卜素、硝酸還原酶等合成通路的深入分析，為揭示相關代謝通路有關基因提供重要理論依據，對于探究毛竹變型稈色變異有重要意義。另外，顏色變異通常是一個不穩定的性狀。例如，花毛竹在不同的生境條件下有可能變回全部綠色或者變成綠皮花毛竹，這表明竹類植物的顏色變異在遺傳上不是一個穩定的性狀。因此，從分子機制上探索顏色變異有其復雜性，其代謝調控還需要進一步研究。

4 結論

采用第二代高通量重測序技術，對4個毛竹變型材料進行全基因組重測序研究，對其單核苷多態性、小片段插入缺失和結構變異進行分析和注釋，篩選可能發生功能變異的基因。將變異基因與GO、COG、KEGG等功能數據庫進行比對，每樣品都有7 000多個變異基因得到功能注釋。GO注釋分類包括細胞組件、分子功能和生物過程3個基因功能分類體系的56個功能組，在細胞組件分類中，葉綠素合成相關基因有2 431個；在生物過程分類中，參與類胡蘿卜素合成過程的基因有75個；在分子功能分類中，參與花青素合成調控以及紫外光下組織中花青素積累的相關基因有80個。COG分類表明參與復制、重組和修復的基因數為369個，信號轉導機制的基因數為291個，轉錄的相關基因為222個。通過KEGG數據庫系統地分析變異基因參與的黃酮類、類胡蘿卜素等物質代謝合成途徑。后續數據的深入分析將解析不同變異類型的基因家族和基因功能，初步闡析不同竹稈變異毛竹變型的分子遺傳基礎。