馬尾松樹高生長(zhǎng)性狀顯著關(guān)聯(lián)基因挖掘

龔桂芳，馮源恒，羅群鳳，楊章旗

（廣西壯族自治區(qū)林業(yè)科學(xué)研究院 廣西優(yōu)良用材林資源培育重點(diǎn)實(shí)驗(yàn)室 國家林業(yè)和草原局馬尾松工程技術(shù)研究中心 廣西馬尾松工程技術(shù)研究中心，廣西南寧 530002）

馬尾松（Pinus massoniana）是我國南方主要的用材樹種之一，其自然分布廣泛，遺傳變異豐富，環(huán)境適應(yīng)性強(qiáng)，曾作為荒山綠化的首選樹種在南方各省區(qū)被大量種植。第八次森林資源清查數(shù)據(jù)表明，我國馬尾松林總面積為1 001萬hm2，蓄積量5．91億m3，其中人工林面積307 萬hm2，蓄積量1．72 億m3[1]。馬尾松也是我國最早一批開展遺傳改良研究的樹種，早在1958年就開始進(jìn)行馬尾松種源試驗(yàn)[2]。1980年后，馬尾松遺傳育種被正式列入國家重點(diǎn)科研項(xiàng)目[3]。進(jìn)入21 世紀(jì)后，各主要產(chǎn)區(qū)的馬尾松遺傳改良研究先后進(jìn)入第2 或第3 個(gè)輪回育種階段[4-7]，但較長(zhǎng)的育種周期制約了馬尾松的遺傳改良進(jìn)程。通過分子標(biāo)記進(jìn)行輔助育種成為當(dāng)前馬尾松育種研究的重點(diǎn)之一。

隨著馬尾松分子遺傳學(xué)研究的深入，可將分子標(biāo)記輔助育種分為兩個(gè)階段[8]。第一階段，以利用SSR 分子標(biāo)記技術(shù)分析育種群體的遺傳結(jié)構(gòu)、估算親緣關(guān)系、對(duì)選自半同胞家系的優(yōu)樹進(jìn)行親本分析為主要內(nèi)容，并將之作為劃分育種群體、研究選擇方法的重要依據(jù)[4,8-10]。第二階段，以基因關(guān)聯(lián)分析為主要技術(shù)手段，開發(fā)SSR、SNPs 等標(biāo)記，與表型性狀進(jìn)行關(guān)聯(lián)分析[11-12]，從而得到高度關(guān)聯(lián)的分子標(biāo)記，開展早期選育和基因功能分析。

本研究基于馬尾松候選基因組關(guān)聯(lián)分析獲得的9 個(gè)與樹高生長(zhǎng)性狀顯著相關(guān)的SSR 分子標(biāo)記，對(duì)其所在的基因組序列區(qū)域進(jìn)行分析，以獲得與馬尾松樹高生長(zhǎng)性狀相關(guān)的重要基因。

1 材料與方法

1．1 顯著關(guān)聯(lián)SSR位點(diǎn)的獲得

研究采用的9個(gè)與馬尾松樹高生長(zhǎng)性狀顯著相關(guān)的SSR分子位點(diǎn)來自研究團(tuán)隊(duì)前期進(jìn)行的候選基因組關(guān)聯(lián)分析研究結(jié)果。前期研究對(duì)兩個(gè)生長(zhǎng)及產(chǎn)脂量性狀存在差異的馬尾松無性系的頂芽組織、針葉及樹干韌皮部進(jìn)行轉(zhuǎn)錄組差異分析，得到差異表達(dá)基因。其目的是從頂芽組織中獲得與抽梢及縱向生長(zhǎng)相關(guān)的候選基因，從針葉中獲得與光合效率相關(guān)的候選基因，從樹干韌皮部獲得與徑向生長(zhǎng)及樹脂分泌相關(guān)的候選基因。在獲得的差異表達(dá)序列中，設(shè)計(jì)開發(fā)259 對(duì)EST-SSR 引物[11]，從中選出65對(duì)SSR引物進(jìn)行候選基因組關(guān)聯(lián)分析[13]。關(guān)聯(lián)分析試驗(yàn)群體由320株1994年造林的馬尾松1代種子園自由授粉獲得的子代組成，來自106個(gè)家系，在同一個(gè)隨機(jī)交配系統(tǒng)下產(chǎn)生，最大母本貢獻(xiàn)率為2．5%，不存在單一親本貢獻(xiàn)率過大的問題。2016年12月進(jìn)行生長(zhǎng)量測(cè)定，獲得樹高表型數(shù)據(jù)，與65 對(duì)SSR 引物進(jìn)行候選基因組關(guān)聯(lián)分析。以P＜0．05 作為標(biāo)記與樹高性狀存在連鎖不平衡的標(biāo)準(zhǔn)，共計(jì)獲得9 個(gè)與馬尾松樹高生長(zhǎng)性狀顯著相關(guān)的SSR 分子位點(diǎn)，平均表型變異解釋率1．42%（表1）。

1．2 顯著關(guān)聯(lián)基因的挖掘方法

1．2．1 顯著關(guān)聯(lián)位點(diǎn)所在基因序列的挖掘

根據(jù)該批次馬尾松EST-SSR 引物[11]的設(shè)計(jì)檔案，查找出9 個(gè)SSR 位點(diǎn)所在的序列ID 號(hào)，獲得9 個(gè)SSR 位點(diǎn)所在的轉(zhuǎn)錄組第2 代測(cè)序基因序列。開展第3 代全長(zhǎng)轉(zhuǎn)錄組測(cè)序，通過序列檢索、比對(duì)，獲得上述基因的全長(zhǎng)序列。

1．2．2 顯著關(guān)聯(lián)基因的功能注釋

為獲得全面的基因功能信息，對(duì)得到的基因序列進(jìn)行7 個(gè)數(shù)據(jù)庫的基因功能注釋，包括：Nr（NCBI non-redundant protein sequences）、Nt（NCBI non-redundant nucleotide sequences）、Pfam（Protein family）、KOG/COG（Clusters of Orthologous Groups of proteins）、Swiss-Prot（A manually annotated and reviewed protein sequence database）、KO（KEGG Ortholog database）和GO（Gene Ontology）。

2 結(jié)果與分析

2．1 顯著關(guān)聯(lián)基因序列挖掘

根據(jù)序列ID 號(hào)，通過Novofinder 軟件在測(cè)序數(shù)據(jù)中檢索獲得SSR 位點(diǎn)所在的轉(zhuǎn)錄組第2 代測(cè)序基因序列。9 個(gè)基因序列平均長(zhǎng)度為1 867 bp，其中PCZ023所在序列最短（910 bp），PCZ129所在序列最長(zhǎng)（3 225 bp）（表2）。

將獲得的基因序列進(jìn)行檢索和比對(duì)，以獲得完整的基因序列。分析結(jié)果表明，9個(gè)基因均在第3代全長(zhǎng)轉(zhuǎn)錄組測(cè)序結(jié)果中比對(duì)到相應(yīng)序列（表3）。9個(gè)基因序列平均長(zhǎng)度為1 706 bp，比第2代測(cè)序結(jié)果略短。其中PCZ023 所在基因序列最短（761 bp），PCZ157 所在基因序列最長(zhǎng)（3 831 bp）。9 個(gè)基因的第2 代測(cè)序結(jié)果與第3 代全長(zhǎng)轉(zhuǎn)錄組測(cè)序結(jié)果比對(duì)一致性平均為99．69%，均具有高度的一致性。

表3 馬尾松樹高性狀顯著關(guān)聯(lián)基因序列比對(duì)情況Tab.3 Sequence alignment of genes significantly associated with tree height traits of P.massoniana

2．2 顯著關(guān)聯(lián)基因的功能注釋

為進(jìn)一步分析與樹高性狀顯著關(guān)聯(lián)基因的功能，將得到的9 個(gè)基因在7 個(gè)數(shù)據(jù)庫中進(jìn)行基因功能注釋分析。功能注釋結(jié)果為PCZ002 所在的基因在云杉（Picea asperata）基因組中存在相似序列，其功能屬于轉(zhuǎn)錄因子類編碼基因；PCZ023所在的基因在云杉基因組中存在相似序列，其功能屬于60S 核糖體蛋白大亞基編碼基因；PCZ090所在的基因在油松（Picea tabuliformis）基因組中存在相似序列，其功能屬于TCHQD 類谷胱甘肽S-轉(zhuǎn)移酶編碼基因；PCZ099所在的基因在云杉基因組中存在相似序列，其功能屬于ATP酶編碼基因；PCZ129所在的基因在白云杉（Picea glauca）基因組中存在相似序列，其功能屬于泛核蛋白編碼基因；PCZ142所在基因在白云杉基因組中存在相似序列，其功能屬于60S 核糖體蛋白大亞基編碼基因；PCZ157所在基因在白云杉基因組中存在相似序列，其功能屬于轉(zhuǎn)錄因子類編碼基因；PCZ187所在的基因在云杉基因組中存在相似序列，其功能屬于氧化還原電子傳遞鏈酶類編碼基因；PCZ187所在的基因在云杉基因組中存在相似序列，其功能未知。

3 討論與結(jié)論

基因組關(guān)聯(lián)分析是基于基因的連鎖不平衡原理，將基因型與觀測(cè)表型進(jìn)行群體水平的統(tǒng)計(jì)學(xué)分析，根據(jù)統(tǒng)計(jì)量或顯著性P值篩選出最有可能影響該性狀的分子標(biāo)記位點(diǎn)，挖掘與性狀變異相關(guān)基因的一種研究方法。所得的標(biāo)記位點(diǎn)極有可能與性狀變異相關(guān)的基因是強(qiáng)度連鎖，甚至處于該基因序列上。本研究基于候選基因組關(guān)聯(lián)分析研究結(jié)果，采用的分子標(biāo)記來自候選基因。通過該方法獲得的與樹高性狀顯著關(guān)聯(lián)的標(biāo)記位點(diǎn)極有可能處在控制該性狀變異的基因上。通過挖掘標(biāo)記位點(diǎn)所在的基因，得到控制樹高生長(zhǎng)主效基因的幾率比采用全基因組關(guān)聯(lián)分析與簡(jiǎn)化基因組關(guān)聯(lián)分析方法更高。

本研究中，挖掘的9 個(gè)與樹高性狀顯著關(guān)聯(lián)的基因均在7個(gè)數(shù)據(jù)庫中檢索到高度同源的基因。其中，5 個(gè)在云杉基因組中發(fā)現(xiàn)同源基因，3 個(gè)在白云杉基因組中發(fā)現(xiàn)同源基因，1 個(gè)在油松基因組中發(fā)現(xiàn)同源基因，說明上述基因可能是松屬植物特有的基因家族類型。

對(duì)馬尾松半雙列雜交家系遺傳測(cè)定試驗(yàn)的分析結(jié)果表明，馬尾松樹高性狀加性效應(yīng)高于顯性效應(yīng)[14]，說明馬尾松樹高性狀是典型的數(shù)量性狀，由眾多的基因甚至基因家族共同控制。樹木的光合作用、呼吸作用、激素調(diào)控、水分及營(yíng)養(yǎng)元素吸收和抗逆性等都會(huì)對(duì)樹高生長(zhǎng)產(chǎn)生重要影響。本研究中，有8 個(gè)基因在數(shù)據(jù)庫中獲得功能注釋。其中PCZ002 與PCZ157 所在基因?qū)儆谵D(zhuǎn)錄因子類，PCZ023 與PCZ142 所在基因?qū)儆?0S 核糖體蛋白大亞基編碼基因類，PCZ099 與PCZ187 所在基因可能參與了能量代謝及跨膜運(yùn)輸，PCZ090所在基因?qū)儆赥CHQD 類谷胱甘肽S-轉(zhuǎn)移酶編碼基因，廣泛參與植物體內(nèi)解毒及抗逆境脅迫等功能[15-17]，說明參與樹高生長(zhǎng)過程的基因種類多樣。值得關(guān)注的是，這些基因多與維持植物細(xì)胞基本功能的基因表達(dá)與代謝功能相關(guān)，而非預(yù)想的與細(xì)胞分裂、植物激素合成相關(guān)。由此推測(cè)，植物細(xì)胞基本代謝功能旺盛是保證植物高生長(zhǎng)的原動(dòng)力。

本研究采用通過第3代測(cè)序技術(shù)開展轉(zhuǎn)錄組測(cè)序獲得基因全長(zhǎng)的方法。第3代測(cè)序技術(shù)實(shí)現(xiàn)DNA聚合酶內(nèi)在自身的延續(xù)性，一個(gè)反應(yīng)就可以測(cè)非常長(zhǎng)的序列。第2代測(cè)序僅可測(cè)上百個(gè)堿基，但第3代測(cè)序可測(cè)幾千個(gè)堿基，并可對(duì)RNA 進(jìn)行直接測(cè)序，大大降低體外逆轉(zhuǎn)錄產(chǎn)生的系統(tǒng)誤差，且精度非常高，達(dá)到99．999 9%，克服了第2 代測(cè)序中因拼接過多造成的錯(cuò)誤。與傳統(tǒng)的通過RACE 技術(shù)等克隆獲得基因全長(zhǎng)的方法相比較，具有一次性獲得海量基因的完整序列明顯優(yōu)勢(shì)。本研究中的9 個(gè)基因，其第2代測(cè)序結(jié)果與第3代測(cè)序結(jié)果均存在長(zhǎng)度差異，這可能是因?yàn)榈?代測(cè)序以其他模式物種為模板進(jìn)行拼接造成的拼接錯(cuò)誤，也可能是因?yàn)镽NA 序列自身在轉(zhuǎn)錄過程中發(fā)生剪切、拼接等。