馬尾松生長性狀顯著關聯標記PCZ090所在基因的挖掘與生物信息學分析

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

（廣西壯族自治區林業科學研究院 廣西優良用材林資源培育重點實驗室 廣西馬尾松工程技術研究中心，廣西南寧 530002）

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

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

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

1 材料與方法

1．1 馬尾松生長性狀顯著關聯SSR位點的獲得

研究采用與馬尾松生長性狀顯著相關的SSR分子位點來自研究團隊前期進行的候選基因組關聯分析研究結果。該研究對兩個生長及產脂量性狀存在差異的馬尾松無性系的頂芽、針葉及樹干韌皮部3 個組織進行轉錄組差異分析，從而得到差異表達基因。在獲得差異表達序列中設計開發259 對EST-SSR 引物[7]。從中選擇65 對SSR 引物進行候選基因組關聯分析。關聯分析試驗群體由320株1994年造林的馬尾松1 代種子園自由授粉子代組成，來自106個家系，在同一個隨機交配系統下產生，最大母本貢獻率為2．5%，不存在單一親本貢獻率過大的問題。2016年12月測定生長量，獲得樹高、胸徑和材積3 組生長性狀表型數據與65 對SSR 引物，進行候選基因組關聯分析。以P＜0．05 作為標記與生長性狀存在連鎖不平衡的標準。共計獲得14 個與馬尾松生長性狀顯著關聯的SSR 分子位點，其中PCZ090 與PCZ187 兩個標記與樹高、胸徑和材積3組生長性狀均顯著關聯。

1．2 馬尾松PCZ090標記所在基因序列的追溯

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

表1 SSR關聯分析結果Tab.1 SSR association analysis results

PCZ090 標記是根據馬尾松轉錄組測序結果開發所得的EST-SSR 標記。從測序結果中可追溯檢索到PCZ090標記的原序列Cluster-7694．28796。

1．3 引物來源及PCR反應條件

根據PCZ090標記的原序列Cluster-7694．28796，設計全長擴增引物PCZ090-F（5′-TTTCCTCCCGAAGTAATTTGC-3′）、PCZ090-R（5′-CACCATGATGGATAAAACCGAA-3′），擴增PCZ090標記的cDNA全長序列，最后將目的片段進行克隆測序。

反應采用50 μL 體系進行：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，Forward 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 反應條件為94 ℃3 min；94 ℃30 s，55 ℃30 s，72 ℃2 min，35 cycles；72 ℃延伸10 min，結束PCR 反應。將獲得的目的條帶送至上海生工生物工程公司進行測序。

1．4 生物信息學分析

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

2 結果與分析

2．1 馬尾松PCZ090標記相關基因序列分析

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

圖1 PmGSTTCHQD1基因ORF結果Fig.1 ORF result of PmGSTTCHQD1 gene

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

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

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

圖2 PmGSTTCHQD1保守結構域Fig.2 Conserved domain of PmGSTTCHQD1

PmGSTTCHQD1 蛋白的親水性氨基酸殘基明顯多于疏水性氨基酸殘基，屬于親水性蛋白質，親水指數為-2．467 ～2．211；親水區域的最高值（Score = -2．467）出現在第153 個氨基酸殘基，疏水區域的最高值（Score=2．211）出現在第65 個氨基酸殘基（圖3）。

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

通過SignalP 進行信號肽預測。結果表明，馬尾松PmGSTTCHQD1基因不存在信號肽酶切位點，是非分泌型蛋白。TMHMM 2．0蛋白質跨膜結構分析結果顯示，該基因所編碼的蛋白沒有跨膜結構（圖4）。

圖4 PmGSTTCHQD1蛋白跨膜結構Fig.4 Transmembrane structure of PmGSTTCHQD1 protein

2．3 編碼蛋白質的二級結構預測

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

圖5 PmGSTTCHQD1蛋白二級結構預測Fig.5 Prediction of secondary structure of PmGSTTCHQD1 protein

2．4 編碼蛋白質的三級結構預測

蛋白質三級結構預測結果顯示，PmGSTTCHQD1基因編碼的蛋白質以SWISS-MODEL 上來自毛果楊（Populus trichocarpa）的谷胱甘肽轉移酶F2 的晶體結構（編號5ey6．1．B）模型作為模板，與模板有較高相似性的是1 ～231位氨基酸，模板覆蓋率為81%，兩者的序列具有33%的一致性（圖6）。

圖6 PmGSTTCHQD1蛋白三級結構模型Fig.6 Predicted tertiary structure model of PmGSTTCHQD1 protein

2．5 同源性與系統進化分析

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

MEGA5．0 系統進化結果表明，PmGSTTCHQD1蛋白與油松和落葉松TCHQD 類谷胱甘肽S-轉移酶蛋白關系最近，推測PmGSTTCHQD1 蛋白與油松和落葉松TCHQD 類谷胱甘肽S-轉移酶蛋白功能類似（圖7）。

圖7 PmGSTTCHQD1蛋白系統進化分析Fig.7 Phylogenetic analysis of PmGSTTCHQD1 protein

3 結論與討論

本文基于高通量轉錄組測序結果對PCZ090 標記所在基因序列進行追溯，獲得了馬尾松PCZ090基因的全長。并通過生物信息學軟件分析SSRPCZ090 標記所在的mRNA 序列，包括同源序列比對、蛋白質一級、二級和三級結構及功能基因的預測。同源序列比對結果表明，該基因是馬尾松谷胱甘肽S-轉移酶家族蛋白編碼基因成員。

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

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

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