馬尾松生長性狀顯著關(guān)聯(lián)標(biāo)記PCZ090所在基因的挖掘與生物信息學(xué)分析

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

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

馬尾松（Pinus massoniana）廣泛分布于我國南方17個(gè)?。▍^(qū)），是我國南方特有的多用途重要工業(yè)用材樹種[1-2]，具有分布面積廣、生長迅速、蓄積量大、經(jīng)濟(jì)價(jià)值高、用途廣和適應(yīng)性強(qiáng)等特點(diǎn)[3]。最新森林資源清查結(jié)果顯示，我國馬尾松分布區(qū)域達(dá)200萬km2，其中林分面積為804．30萬hm2，占全國喬木林面積的4．47%，林分總面積居全國喬木樹種首位，蓄積量居第六位[4]。馬尾松生長迅速，可在較短的經(jīng)營期內(nèi)提供大量木材，是許多木材工業(yè)、林產(chǎn)工業(yè)和造紙工業(yè)的支柱，同時(shí)也是荒山綠化和生態(tài)林建設(shè)的先鋒樹種[5]。通過選育和推廣馬尾松良種，可提高單位面積木材產(chǎn)量和質(zhì)量，緩解木材供需矛盾、保障生態(tài)安全及增加林業(yè)產(chǎn)業(yè)效益[6]。因此加快高材性狀的選育已成為馬尾松遺傳改良的重要內(nèi)容。

從20 世紀(jì)80年代開始，廣西壯族自治區(qū)林業(yè)科學(xué)研究院松樹課題組就開始了馬尾松生長性狀遺傳改良研究工作。至今已選育高材馬尾松將近2 000 份，建立了5 個(gè)馬尾松2 代種子園，3 個(gè)馬尾松高世代育種園。馬尾松的育種周期相對(duì)較長，通過傳統(tǒng)的遺傳育種方法，短時(shí)間內(nèi)無法做到準(zhǔn)確快速選育，難以實(shí)現(xiàn)良種的快速更新，限制了馬尾松的遺傳改良進(jìn)程，制約了森林生產(chǎn)力的提高速率。因此，為更快選育大量馬尾松高材良種，滿足林業(yè)生產(chǎn)的需求，開發(fā)研究與生長性狀具有較強(qiáng)相關(guān)性的遺傳標(biāo)記，結(jié)合分子標(biāo)記輔助選擇，是加快馬尾松生長性狀改良進(jìn)程的一種可行方法。通過分子標(biāo)記輔助選擇開展早期選擇，可有效縮短育種周期，提高選育效率。建立馬尾松分子標(biāo)記輔助選擇技術(shù)體系對(duì)于高效開展馬尾松遺傳改良意義重大。

課題組前期進(jìn)行了馬尾松主要經(jīng)濟(jì)性狀候選基因組關(guān)聯(lián)分析研究，已開發(fā)獲得與生長性狀顯著關(guān)聯(lián)SSR 標(biāo)記14 個(gè)。本研究以其中一個(gè)位點(diǎn)（PCZ090）為目標(biāo)，基于高通量轉(zhuǎn)錄組測(cè)序結(jié)果對(duì)其所在基因序列進(jìn)行追溯，測(cè)序驗(yàn)證獲得其基因全長，并通過生物信息學(xué)軟件分析PCZ090 基因序列，包括同源序列比對(duì)、蛋白結(jié)構(gòu)（一級(jí)、二級(jí)和三級(jí)結(jié)構(gòu)）及結(jié)構(gòu)域和功能位點(diǎn)預(yù)測(cè)等。對(duì)生長性狀顯著關(guān)聯(lián)標(biāo)記所在基因進(jìn)行深度挖掘，是發(fā)現(xiàn)控制目標(biāo)性狀重要基因的有效途徑，對(duì)于深層次探究馬尾松生長性狀遺傳調(diào)控機(jī)理具有重要意義。本研究可為進(jìn)一步揭示馬尾松生長性狀相關(guān)基因的作用效應(yīng)、深入研究馬尾松主要生長性狀遺傳機(jī)理、分子標(biāo)記輔助育種技術(shù)在其他林木上的應(yīng)用提供參考。

1 材料與方法

1．1 馬尾松生長性狀顯著關(guān)聯(lián)SSR位點(diǎn)的獲得

研究采用與馬尾松生長性狀顯著相關(guān)的SSR分子位點(diǎn)來自研究團(tuán)隊(duì)前期進(jìn)行的候選基因組關(guān)聯(lián)分析研究結(jié)果。該研究對(duì)兩個(gè)生長及產(chǎn)脂量性狀存在差異的馬尾松無性系的頂芽、針葉及樹干韌皮部3 個(gè)組織進(jìn)行轉(zhuǎn)錄組差異分析，從而得到差異表達(dá)基因。在獲得差異表達(dá)序列中設(shè)計(jì)開發(fā)259 對(duì)EST-SSR 引物[7]。從中選擇65 對(duì)SSR 引物進(jìn)行候選基因組關(guān)聯(lián)分析。關(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月測(cè)定生長量，獲得樹高、胸徑和材積3 組生長性狀表型數(shù)據(jù)與65 對(duì)SSR 引物，進(jìn)行候選基因組關(guān)聯(lián)分析。以P＜0．05 作為標(biāo)記與生長性狀存在連鎖不平衡的標(biāo)準(zhǔn)。共計(jì)獲得14 個(gè)與馬尾松生長性狀顯著關(guān)聯(lián)的SSR 分子位點(diǎn)，其中PCZ090 與PCZ187 兩個(gè)標(biāo)記與樹高、胸徑和材積3組生長性狀均顯著關(guān)聯(lián)。

1．2 馬尾松PCZ090標(biāo)記所在基因序列的追溯

胸徑和材積均通過3 種線性模型進(jìn)行分析，分別為EMMMAX、GAPIT 和GLM，計(jì)算得到P值，胸徑P值分別為0．012 368、0．002 769和0．000 139；材積P值分別為0．006 429、0．002 796 和6．27E-07；樹高只通過GLM 模型進(jìn)行分析，計(jì)算得到P值為0．043 690（表1）。3 種生長性狀的P值都小于0．05，說明該標(biāo)記與生長性狀顯著關(guān)聯(lián)，符合要求。

表1 SSR關(guān)聯(lián)分析結(jié)果Tab.1 SSR association analysis results

PCZ090 標(biāo)記是根據(jù)馬尾松轉(zhuǎn)錄組測(cè)序結(jié)果開發(fā)所得的EST-SSR 標(biāo)記。從測(cè)序結(jié)果中可追溯檢索到PCZ090標(biāo)記的原序列Cluster-7694．28796。

1．3 引物來源及PCR反應(yīng)條件

根據(jù)PCZ090標(biāo)記的原序列Cluster-7694．28796，設(shè)計(jì)全長擴(kuò)增引物PCZ090-F（5′-TTTCCTCCCGAAGTAATTTGC-3′）、PCZ090-R（5′-CACCATGATGGATAAAACCGAA-3′），擴(kuò)增PCZ090標(biāo)記的cDNA全長序列，最后將目的片段進(jìn)行克隆測(cè)序。

反應(yīng)采用50 μL 體系進(jìn)行：TaKaRa LATaq（5 U/μL）0．5 μL，10 × LA PCR Buffer Ⅱ（Mg2+Free）5．0 μL，MgCl2（25 mmol/L）5．0 μL，dNTP Mixture（2．5 mmol/L each）8．0 μL，F(xiàn)orward Primer（10 μmol/L）2．0 μL，Reverse Primer（10 μmol/L）2．0 μL，cDNA template 1．0 μL，Milli-Q Water 26．5 μL。PCR 反應(yīng)條件為94 ℃3 min；94 ℃30 s，55 ℃30 s，72 ℃2 min，35 cycles；72 ℃延伸10 min，結(jié)束PCR 反應(yīng)。將獲得的目的條帶送至上海生工生物工程公司進(jìn)行測(cè)序。

1．4 生物信息學(xué)分析

采用ORF Finding（www．ncbi．nlm．nih．gov/orffinder）在線軟件進(jìn)行ORF 確定及氨基酸序列預(yù)測(cè)；采用在線工具Expasy ProtParam（https://web．expasy．org/protparam/）分析蛋白的分子量、理論等電點(diǎn)和氨基酸組成等；采用InterProScan（https://www．ebi．a(chǎn)c．uk/interpro/）、NCBI CDD（http://www．ncbi．nlm．nih．gov/Structure/cdd/wrpsb．cgi）在線工具軟件進(jìn)行結(jié)構(gòu)域和功能位點(diǎn)預(yù)測(cè)；采用ProtScale（https://web．expasy．org/protscale/）在線軟件繪制蛋白質(zhì)的親疏水性系列譜；采用SignalP（https://services．healthtech．dtu．dk/service．php? SignalP-5．0）在線軟件預(yù)測(cè)信號(hào)肽；采用TMHMM Server v．2．0（https://services．healthtech．dtu．dk/service．php?TMHMM-2．0）程序進(jìn)行蛋白序列跨膜區(qū)分析；采用SOPMA（https://npsa-prabi．ibcp．fr/cgibin/npsa_automat．pl?page=npsa_sopma．html）在線 軟件預(yù)測(cè)分析蛋白質(zhì)二級(jí)結(jié)構(gòu)[8]；采用SWISS-MODEL（https://swissmodel．expasy．org/interactive）在 線軟 件對(duì)蛋白質(zhì)三級(jí)結(jié)構(gòu)進(jìn)行建模[9-11]；采用NCBI BLAST（https://blast．ncbi．nlm．nih．gov/Blast．cgi）中 的Nucleotide BLAST 進(jìn)行同源比對(duì)；采用MEGA 5．0 構(gòu)建系統(tǒng)進(jìn)化樹，選擇鄰接法，Bootstrap 檢驗(yàn)使用1 000 次重復(fù)。

2 結(jié)果與分析

2．1 馬尾松PCZ090標(biāo)記相關(guān)基因序列分析

基因全長1 636 bp，其中5′非編碼區(qū)有432 bp，3′非編碼區(qū)有409 bp，開放閱讀框（ORF）序列為795 bp（433 ～1 227 bp），編碼264 個(gè)氨基酸，起始密碼子為ATG（黃色標(biāo)注），終止密碼子為TGA（黃色標(biāo)注）（圖1）。

圖1 PmGSTTCHQD1基因ORF結(jié)果Fig.1 ORF result of PmGSTTCHQD1 gene

PCZ090 所在基因編碼的蛋白質(zhì)分子式為C1411H2183N377O387S8，編碼蛋白的分子量約為30．88 kDa，理論等電點(diǎn)為9．27，屬于堿性蛋白，原子總數(shù)為4 366。20 種氨基酸組成中，亮氨酸（Leu）數(shù)量最多（13．6%）；半胱氨酸（Cys）數(shù)量最少（0．4%）。酸性氨基酸（Asp+Glu）總數(shù)為29 個(gè)，堿性氨基酸（Arg+Lys）總數(shù)為36 個(gè)，說明該蛋白帶弱正電荷。不穩(wěn)定指數(shù)為46．23，該蛋白為不穩(wěn)定蛋白；脂肪指數(shù)為94．28，親水性總平均值為-0．333。

2．2 編碼蛋白產(chǎn)物的功能分析和親疏水性比較

根據(jù)InterProScan 在線軟件分析基因編碼產(chǎn)物的結(jié)構(gòu)域和功能位點(diǎn)，發(fā)現(xiàn)PCZ090 所在基因具有典型的谷胱甘肽S-轉(zhuǎn)移酶的保守結(jié)構(gòu)域：GST-N 末端結(jié)構(gòu)域（從Met1-Gly80）、GST-C 末端結(jié)構(gòu)域（從Thr148-Pro229）（圖2），參與谷胱甘肽代謝過程，具有谷胱甘肽轉(zhuǎn)移酶活性。通過NCBI CDD 進(jìn)一步確定PCZ090 所在基因是谷胱甘肽S-轉(zhuǎn)移酶家族蛋白成員，位于1 ～234 位氨基酸，包含N 端結(jié)構(gòu)域和C 端結(jié)構(gòu)域，可催化還原谷胱甘肽與多種內(nèi)源性和外源性烷化劑的結(jié)合，包括致癌物、治療藥物、環(huán)境毒素和氧化應(yīng)激產(chǎn)物。依據(jù)GST 的系統(tǒng)命名體系，將該基因命名為PmGSTTCHQD1。

圖2 PmGSTTCHQD1保守結(jié)構(gòu)域Fig.2 Conserved domain of PmGSTTCHQD1

PmGSTTCHQD1 蛋白的親水性氨基酸殘基明顯多于疏水性氨基酸殘基，屬于親水性蛋白質(zhì)，親水指數(shù)為-2．467 ～2．211；親水區(qū)域的最高值（Score = -2．467）出現(xiàn)在第153 個(gè)氨基酸殘基，疏水區(qū)域的最高值（Score=2．211）出現(xiàn)在第65 個(gè)氨基酸殘基（圖3）。

圖3 PmGSTTCHQD1基因編碼蛋白的親疏水性序列譜Fig.3 Hydropathy profile of protein coded by PmGSTTCHQD1 gene

通過SignalP 進(jìn)行信號(hào)肽預(yù)測(cè)。結(jié)果表明，馬尾松PmGSTTCHQD1基因不存在信號(hào)肽酶切位點(diǎn)，是非分泌型蛋白。TMHMM 2．0蛋白質(zhì)跨膜結(jié)構(gòu)分析結(jié)果顯示，該基因所編碼的蛋白沒有跨膜結(jié)構(gòu)（圖4）。

圖4 PmGSTTCHQD1蛋白跨膜結(jié)構(gòu)Fig.4 Transmembrane structure of PmGSTTCHQD1 protein

2．3 編碼蛋白質(zhì)的二級(jí)結(jié)構(gòu)預(yù)測(cè)

SOPMA 分析表明，PmGSTTCHQD1基因編碼的蛋白質(zhì)二級(jí)結(jié)構(gòu)含有57．20%（151 個(gè)）的α 螺旋（Alpha helix）、5．68%（15 個(gè)）的β 轉(zhuǎn) 角（Beta turn）、10．98%（29個(gè)）的延伸鏈（Extended strand）和26．14%（69個(gè)）的無規(guī)則卷曲（Random coil）（圖5）。其中，α螺旋和無規(guī)則卷曲是主要成分。

圖5 PmGSTTCHQD1蛋白二級(jí)結(jié)構(gòu)預(yù)測(cè)Fig.5 Prediction of secondary structure of PmGSTTCHQD1 protein

2．4 編碼蛋白質(zhì)的三級(jí)結(jié)構(gòu)預(yù)測(cè)

蛋白質(zhì)三級(jí)結(jié)構(gòu)預(yù)測(cè)結(jié)果顯示，PmGSTTCHQD1基因編碼的蛋白質(zhì)以SWISS-MODEL 上來自毛果楊（Populus trichocarpa）的谷胱甘肽轉(zhuǎn)移酶F2 的晶體結(jié)構(gòu)（編號(hào)5ey6．1．B）模型作為模板，與模板有較高相似性的是1 ～231位氨基酸，模板覆蓋率為81%，兩者的序列具有33%的一致性（圖6）。

圖6 PmGSTTCHQD1蛋白三級(jí)結(jié)構(gòu)模型Fig.6 Predicted tertiary structure model of PmGSTTCHQD1 protein

2．5 同源性與系統(tǒng)進(jìn)化分析

由NCBI 核酸數(shù)據(jù)庫Blast 比對(duì)結(jié)果可知，PmGSTTCHQD1基因與油松（Pinus tabuliformis）TCHQD類谷胱甘肽S-轉(zhuǎn)移酶（JX962828．1）和日本落葉松（Larix kaempferi）TCHQD 類谷胱甘肽S-轉(zhuǎn)移酶（KF438050．1）的完整mRNA 同源性分別為99．75%和93．09%，覆蓋率均為48%；與白云杉（Picea glauca）GQ02824_L21（BT105585．1）的mRNA 同源性為91．01%，覆蓋率為49%；與火炬松（Pinus taeda）兩個(gè)分 離 株5639 和5637（FJ066498．1 和FJ066514．1）的基因組序列同源性分別為99．31%和99．14%，覆蓋率分別為26%和21%。

MEGA5．0 系統(tǒng)進(jìn)化結(jié)果表明，PmGSTTCHQD1蛋白與油松和落葉松TCHQD 類谷胱甘肽S-轉(zhuǎn)移酶蛋白關(guān)系最近，推測(cè)PmGSTTCHQD1 蛋白與油松和落葉松TCHQD 類谷胱甘肽S-轉(zhuǎn)移酶蛋白功能類似（圖7）。

圖7 PmGSTTCHQD1蛋白系統(tǒng)進(jìn)化分析Fig.7 Phylogenetic analysis of PmGSTTCHQD1 protein

3 結(jié)論與討論

本文基于高通量轉(zhuǎn)錄組測(cè)序結(jié)果對(duì)PCZ090 標(biāo)記所在基因序列進(jìn)行追溯，獲得了馬尾松PCZ090基因的全長。并通過生物信息學(xué)軟件分析SSRPCZ090 標(biāo)記所在的mRNA 序列，包括同源序列比對(duì)、蛋白質(zhì)一級(jí)、二級(jí)和三級(jí)結(jié)構(gòu)及功能基因的預(yù)測(cè)。同源序列比對(duì)結(jié)果表明，該基因是馬尾松谷胱甘肽S-轉(zhuǎn)移酶家族蛋白編碼基因成員。

植物谷胱甘肽S-轉(zhuǎn)移酶（GSTs）是一類種類豐富且普遍存在的酶[12]。其參與多種細(xì)胞功能，如初生代謝、次生代謝、解毒作用、抗低溫、抗氧化損傷及異源物質(zhì)隔離等[13]。該類酶還能作為配體蛋白參與植物激素代謝途徑[14]。在高等植物中，GST 被分為八種亞類：phi、tau、theta、zeta、lambda、EF1Bγ（elongation factor 1gamma）、DHAR（Dehydroascorbate reductase，谷胱甘肽依賴的抗壞血酸還原酶）和TCHQD（Tetrachlorohydroquinone dehalogenase，四氯代氫醌脫鹵素酶）[14-15]。本研究挖掘得到馬尾松PmGSTTCHQD1基因?qū)儆赥CHQD 亞家族中的一員。該類谷胱甘肽S-轉(zhuǎn)移酶的研究尚不全面，在模式植物擬南芥（Arabidopsis thaliana）中僅鑒定到1 個(gè)TCHQD 類GST，因其與原核酶序列同源，推測(cè)其可以催化氯代異源物質(zhì)的代謝[16]。

GST 作為一種重要的解毒酶類，其催化反應(yīng)主要分3 步進(jìn)行：（1）催化特異性底物進(jìn)行脫鹵水反應(yīng)；（2）已被活化的疏水底物與谷胱甘肽過氧化物酶（GSH）結(jié)合，增加底物的親水性；（3）復(fù)合物被轉(zhuǎn)運(yùn)出細(xì)胞質(zhì)。這一過程在植物中是由ATP 驅(qū)動(dòng)泵驅(qū)動(dòng)復(fù)合物進(jìn)入液泡或者非共質(zhì)體中[17]。在本研究團(tuán)隊(duì)通過候選基因組關(guān)聯(lián)分析得到的與生長性狀顯著關(guān)聯(lián)的其他位點(diǎn)中，PCZ099所在基因即為ATP酶基因。 由此推測(cè)PCZ090 所在基因PmGSTTCHQD1可能參與有毒物質(zhì)的脫鹵水反應(yīng)，而PCZ099 所在的ATP 酶基因則參與復(fù)合物被轉(zhuǎn)運(yùn)出細(xì)胞質(zhì)的過程，兩者共同參與細(xì)胞解毒功能。細(xì)胞解毒能力很有可能影響馬尾松的生長。

林木的生長性狀是典型的數(shù)量性狀，由多基因控制。本研究得到的PmGSTTCHQD1基因是通過關(guān)聯(lián)分析研究發(fā)現(xiàn)的對(duì)生長性狀表型變異解釋率較高的基因之一，其編碼蛋白所屬的GST 酶類主要通過參與細(xì)胞解毒機(jī)制以提高抗性。其對(duì)生長性狀的作用還需深入研究。可對(duì)關(guān)聯(lián)分析群體中生長性狀存在顯著差異的個(gè)體進(jìn)行PmGSTTCHQD1基因表達(dá)量分析、編碼蛋白活性分析等。并開展PCZ090位點(diǎn)不同基因型間表型效應(yīng)分析，研究其在不同林齡、不同地點(diǎn)的表型差異穩(wěn)定，為利用該標(biāo)記進(jìn)行分子標(biāo)記輔助選擇提供參考。