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

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

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

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

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

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

1 材料與方法

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

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

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

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

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

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

為獲得全面的基因功能信息，對得到的基因序列進行7 個數(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 號，通過Novofinder 軟件在測序數(shù)據(jù)中檢索獲得SSR 位點所在的轉(zhuǎn)錄組第2 代測序基因序列。9 個基因序列平均長度為1 867 bp，其中PCZ023所在序列最短（910 bp），PCZ129所在序列最長（3 225 bp）（表2）。

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

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

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

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

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

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

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