基于SLAF-seq 的華山松球果數(shù)量性狀全基因組關(guān)聯(lián)分析*

李啟少，葉 鵬，趙文植，馬路遙，王正德，曹正英，王 飛，辛培堯

(1.西南林業(yè)大學(xué)，國家林業(yè)和草原局西南風(fēng)景園林工程技術(shù)研究中心，云南 昆明 650224；2.西南林業(yè)大學(xué)，西南地區(qū)生物多樣性保育國家林業(yè)和草原局重點實驗室，云南 昆明 650224；3.內(nèi)蒙古渾善達(dá)克規(guī)模化林場，內(nèi)蒙古 呼和浩特 010020)

華山松(Pinus armandiiFranch.)是松科(Pinaceae)松屬(Pinus)大喬木，也是中國特有的五針?biāo)蓸浞N，其耐寒能力強(qiáng)，喜好溫和涼爽、濕潤的氣候，適應(yīng)性較廣[1]。華山松是優(yōu)良的建筑和家具材料，同時也可以用于制漿造紙，種子可供食用和榨油，利用價值較高。近年來，由于種質(zhì)資源退化、良種選育滯后，導(dǎo)致華山松人工林陸續(xù)出現(xiàn)生長量下降、長勢衰弱、病蟲害嚴(yán)重、利用水平低和產(chǎn)品附加值低等問題[2-4]。雖然國內(nèi)相關(guān)科研單位已開展了對華山松的遺傳改良工作[5-9]，但由于其育種周期長和主要性狀的優(yōu)劣難以在早期進(jìn)行選擇，因此，找到準(zhǔn)確、可靠且有效的早期選擇標(biāo)記是解決良種問題的關(guān)鍵技術(shù)之一。

隨著分子生物學(xué)的不斷進(jìn)步，林木育種的方法與理論也在不斷發(fā)展，尤其是分子標(biāo)記輔助選擇技術(shù)的不斷更新，涌現(xiàn)出多種分子標(biāo)記方法，如RAPD、AFLP、SRAP、SSR 和SNP 等標(biāo)記，為植物遺傳改良提供了更多的選擇，也提高了選擇的可靠性與準(zhǔn)確度[10]。其中，單核苷酸多態(tài)性(single-nucleotide polymorphism，SNP)標(biāo)記是由單堿基突變(缺失、插入、轉(zhuǎn)換和顛倒)所引起的基因組水平上的DNA 序列多態(tài)性，具有穩(wěn)定性高、位點豐富、分布較廣和檢測迅速等特點，目前已廣泛應(yīng)用于植物育種[11]。全基因組關(guān)聯(lián)分析(genome-wide association study，GWAS)由RISCH等[12]首先提出。該技術(shù)能夠在全基因水平上發(fā)掘與復(fù)雜性狀相關(guān)聯(lián)的遺傳變異，以連鎖不平衡為基礎(chǔ)，通過識別和定位目標(biāo)群體中的高密度分子標(biāo)記獲得與復(fù)雜性狀表現(xiàn)型變異有關(guān)聯(lián)的SNP 標(biāo)記[13]。GWAS 最早應(yīng)用于人類視網(wǎng)膜黃斑變性與人類年齡的相關(guān)性分析上[14]。近年來，由于引入混合線性模型和Bonferroni 校正統(tǒng)計方法，能夠?qū)⒆匀蝗后w特征性狀和非家系群體作為研究對象開展GWAS 研究。因此，GWAS 也被大量應(yīng)用于植物的關(guān)聯(lián)分析與遺傳育種工作。目前，不僅在水稻[15]、小麥[16]和玉米[17]的產(chǎn)量、株高、穗長、成花時間及抽穗期等性狀關(guān)聯(lián)出多個相關(guān)的基因位點，還在植物抗病基因的篩選[18]和代謝物合成[19]等方面得到了應(yīng)用。GWAS 相關(guān)研究結(jié)果為闡明生物復(fù)雜性狀的遺傳結(jié)構(gòu)提供了理論基礎(chǔ)。生產(chǎn)中可將檢測到的關(guān)聯(lián)位點運用于分子標(biāo)記輔助選擇，對植物的遺傳改良工作具有重要的理論和實踐意義。

本研究基于前期通過SLAF-seq 測序獲得的高通量SNP[6]構(gòu)建華山松樣品的系統(tǒng)發(fā)育樹與主成分聚類圖，繼而進(jìn)行全基因組關(guān)聯(lián)分析，篩選與華山松球果數(shù)量性狀顯著關(guān)聯(lián)的SNP 位點，并對顯著關(guān)聯(lián)位點區(qū)域的基因進(jìn)行分析，挖掘與球果數(shù)量性狀相關(guān)的候選基因，分析候選基因在華山松6 個部位的qPCR 基因組織特異性表達(dá)。研究結(jié)果可為華山松球果數(shù)量相關(guān)性狀的早期選擇提供理論依據(jù)。

1 材料與方法

1.1 采樣地概況

采樣地位于云南省楚雄市紫溪山華山松無性系種子園(N24°58′58″～25°4′，E101°22′29″～101°26′7″，海拔2 200～2 400 m)，于1990 年建成，園區(qū)內(nèi)80%的林地坡度小于10°。種子園地處北亞熱帶溫涼濕潤氣候，年降水量約1 000 mm，相對濕度80%～85%，年平均氣溫12 ℃，平均最低氣溫-6 ℃，平均最高氣溫20 ℃。

1.2 供試材料

試驗材料選自種子園內(nèi)6 個種源地的1 453株華山松單株，以連續(xù)多年的年份內(nèi)球果數(shù)量為基礎(chǔ)數(shù)據(jù)，共篩選出210 株華山松作為試驗材料，其中高、低球果產(chǎn)量的華山松單株各30株，球果數(shù)量及生長量均處中上水平的華山松單株150 株。另選3 株不同華山松優(yōu)良單株的幼嫩且無病害的球果、樹皮、韌皮部、木質(zhì)部、根和針葉，用于提取華山松RNA 以進(jìn)行q-PCR 驗證。

1.3 試驗方法

1.3.1 系統(tǒng)發(fā)育樹與主成分分析

利用前期已開發(fā)的SNP 標(biāo)記，采用MGEA X軟件[20]，基于鄰接法的Kimura2-parameter 模型構(gòu)建華山松樣本的系統(tǒng)發(fā)育樹，bootstrap 重復(fù)1 000次。利用EIGENSOFT 軟件[21]將SNP 標(biāo)記用于主成分分析，獲得華山松的主成分聚類圖，每1 種顏色代表1 個種源地，每1 個點代表1 個樣品，樣品間存在廣泛的差異，樣品距離越近，則代表親緣關(guān)系越近，主成分分析的結(jié)果可以輔助進(jìn)化分析，充分體現(xiàn)群體中個體親緣關(guān)系的離散程度。

表1 qPCR 引物序列Tab.1 qPCR primer sequence

1.3.2 球果數(shù)量的GWAS 分析

(1) GWAS 分析：將球果數(shù)量作為鑒定指標(biāo)，利 用TASSEL[22]、FaST-LMM[23]和EMMAX[24]軟件的光譜變換線性混合模型(factored spectrally transformed linear mixed models，F(xiàn)aST-LMM)、一般線性模型(general linear models，GLM)、混合線性模型(mixed linear models，MLM)和壓縮混合線性模型(compressed linear mixed model，CMLM)，在充分考慮210 份華山松單株的總體結(jié)構(gòu)與親緣關(guān)系的情況下進(jìn)行關(guān)聯(lián)分析。利用admixture 軟件[25]獲得Q值(樣本總體結(jié)構(gòu))；樣品間親緣關(guān)系K值則是利用SPAGeDi 軟件[26]獲得；X表示基因型，y表示表現(xiàn)型。所有SNP 位點均可得到1 個關(guān)聯(lián)結(jié)果。其中，TASSEL 軟件的混合線性模型公式為：y=Xα＋Qβ＋Kμ＋e。式中：α 為固定效應(yīng)向量；β 為SNP 替代效應(yīng)；μ 為服從分布N (0，Gσ2)的隨機(jī)加性遺傳效應(yīng)的向量，G 為從SNP 標(biāo)記導(dǎo)出的親屬關(guān)系矩陣，σ2為加性方差；e 為殘差。GWAS 結(jié)果基于R v4.2.0 中的CMplot 程序包繪制，Q-Q 圖(quantile-quantile plot)用于評估關(guān)聯(lián)分析結(jié)果的準(zhǔn)確性。

(2) 候選基因篩選：關(guān)聯(lián)分析得到與球果數(shù)量性狀相關(guān)的SNP 位點，并在SLAF 標(biāo)簽上找出對應(yīng)的位置以及堿基序列，利用SNP 序列在NCBI數(shù)據(jù)庫進(jìn)行BLAST 比對，篩選與球果數(shù)量性狀相關(guān)的基因。

1.3.3 熒光定量PCR 反應(yīng)

(1) 華山松總RNA 提取

使用植物RNA 提取試劑盒(OMEGA R6827)提取華山松針葉、樹皮、韌皮部、木質(zhì)部、球果和根的總RNA。

(2) mRNA 反轉(zhuǎn)錄為cDNA

參照周延清等[27]的方法得到第1 鏈cDNA 并儲存于-20 ℃冰箱中，備用。

(3) 引物設(shè)計

從NCBI 中檢索相關(guān)的基因序列，用Primer 3設(shè)計熒光定量PCR 引物(表1)。引物條件為：引物長度在18～24 bp 之間；GC 含量在40%～60%之間；5′端或中間區(qū)域為G 或C；擴(kuò)增片段的長度在100～300 bp 之間，退火溫度在55～60 ℃之間。從設(shè)計的7 對引物中篩選1 對能讓候選基因在各個組織中穩(wěn)定表達(dá)的引物，對檢測到的目標(biāo)基因進(jìn)行q-PCR 分析。

(4) 熒光定量PCR 反應(yīng)

使 用EvaGreen 2×qPCR MasterMix 試劑盒(上海愛必夢生物科技有限公司)進(jìn)行熒光定量PCR 反應(yīng)，反應(yīng)體系為20 μL，包括EvaGreen2×qPCR MasterMix 10 μL，上、下游引物各1 μL，模板DNA 1 μL，無酶水7 μL。q-PCR 擴(kuò)增程序為：95 ℃預(yù)變性 10 min；95 ℃ PCR 反應(yīng)15 s，60 ℃ 1 min，共40 個循環(huán)；95 ℃溶解曲線15 s，60 ℃ 1 min，95 ℃ 15 s。選擇和松樹相關(guān)的AT5G42190作為內(nèi)參基因，其引物序列為：F：ATGCTGGACAGGCTTTGAAC，R：GAGTTGCTCCGAGATCTTTACA。

1.3.4 數(shù)據(jù)統(tǒng)計與分析

試驗結(jié)束后，統(tǒng)計Ct 值于Excel 表中，利用公式ΔΔCt=(Ct處理組-Ct內(nèi)參)-(Ct對照組-Ct內(nèi)參)]計算各組織的Ct 值[28]，之后利用Excel 軟件作圖；采用Excel 2010 和SPASS 20.0 軟件進(jìn)行相關(guān)數(shù)據(jù)的描述統(tǒng)計與相關(guān)性分析。

2 結(jié)果與分析

2.1 系統(tǒng)發(fā)育樹與主成分分析

系統(tǒng)發(fā)育樹(圖1)顯示：華山松各材料之間存在著較為豐富的遺傳變異，來自相近地理位置的材料都基本在相近的位置，但是來源于楚雄和巍山的材料分布相當(dāng)廣泛，在各群體中均有分布。將SNP 標(biāo)記用于主成分分析(圖2)發(fā)現(xiàn)：同一種源地的樣品基本聚在一起，說明它們的親緣關(guān)系較近，其中來自巍山的樣品分布范圍較廣泛，表明巍山各材料間遺傳差異較大，與系統(tǒng)發(fā)育樹的分析結(jié)果基本一致。

圖1 華山松系統(tǒng)發(fā)育樹Fig.1 Phylogenetic tree of Pinus armandii

圖2 華山松PCA 聚類圖Fig.2 PCA cluster graph of P.armandii

2.2 華山松球果數(shù)量性狀的GWAS 結(jié)果

2.2.1 華山松全基因組關(guān)聯(lián)分析

GWAS 結(jié)果表明：球果數(shù)量性狀能夠在FaSTLMM、GLM 和MLM 模型中被定位。由圖3 可知：FaST-LMM 和MLM 模型的觀測值與期望值僅在最右側(cè)發(fā)生了偏離現(xiàn)象，說明球果數(shù)量性狀的差異并非因群體分層所造成，所選的模型為本試驗適合的模型；而GLM 模型中SNP 位點出現(xiàn)在預(yù)測線上方，且不與對角線重合，說明觀測值大大超過了預(yù)測值，表明該模型不合理。

圖3 華山松球果數(shù)量Q-Q 圖Fig.3 Quantile-quantile (Q-Q) plot of cones of P.armandii

共檢測到12 個與球果數(shù)量性狀顯著相關(guān)的SNP 位點，獲得12 個SNP 所在SLAF 標(biāo)簽上的位置和P值(表2)。12 個SNP 位點分別為Marker170924、Marker193307、Marker357784、Marker6337143、Marker17681543、Marker201128、Marker215196、Marker250033、Marker251795、Marker-6508087、Marker292650、Marker18895864。

表2 與華山松球果數(shù)量顯著相關(guān)的SNP 位點Tab.2 SNP loci significantly related to the number of cones in Pinus armandii

2.2.2 候選基因的篩選

將12 個SNP 所在的序列片段在NCBI 中進(jìn)行基因檢索(表3)，在Marker193307 附近檢測到1 個緊密關(guān)聯(lián)的基因Korrigan(KOR1，編碼跨膜內(nèi)切β-1,4-D 葡聚糖酶)，其序列相似度較高(為74.74%)，NCBI 基因登錄號為AC241332.1，該基因?qū)χ参锢w維素合成具有重要作用。

表3 與球果數(shù)量相關(guān)的單核苷酸多態(tài)性(SNP)序列Tab.3 Single-nucleotide polymorphism (SNP) sequence related to the number of cones

2.3 基因組織特異性分析

由圖4 可知：KOR1在韌皮部的表達(dá)量最高，在幼嫩球果中的表達(dá)量次之，在根中的表達(dá)量最低。已知KOR1基因是與纖維素合成相關(guān)的基因，該基因在韌皮部和幼嫩球果中的表達(dá)量較高，說明其參與了韌皮部與幼嫩球果的生長代謝。因此，推測球果數(shù)量性狀與該基因相關(guān)。

圖4 KOR1 在華山松各個部位中的表達(dá)量變化Fig.4 Expression changes of KOR1 gene in different tissus of P.armandii

3 討論

3.1 華山松不同分子標(biāo)記方法比較

近年來，分子標(biāo)記輔助選擇技術(shù)不斷更新，越來越多的分子標(biāo)記方法應(yīng)用于華山松育種。趙楊等[29]對2 個不同年份的華山松子代材料進(jìn)行ISSR 遺傳多樣性分析，結(jié)果顯示：該種子園的遺傳多樣性處于較高水平，遺傳變異多處于種源內(nèi)，與種子園內(nèi)親代群體相比較，初期子代的遺傳多樣性指數(shù)略有降低，而后期子代的遺傳多樣性指數(shù)則為最低，可見，對于提高華山松球果產(chǎn)量和品質(zhì)刻不容緩。辛靜等[30]對紫溪山華山松無性系種子園進(jìn)行SSR 遺傳多樣性分析，結(jié)果表明：種源間遺傳變異不高，其主要遺傳變異存在于種源內(nèi)。徐劍等[31]利用SCoT 分子標(biāo)記對紫溪山華山松種子園進(jìn)行遺傳多樣性分析，表明種子園內(nèi)不同無性系擁有較高水平的遺傳多樣性，且遺傳分化主要存在于種源內(nèi)。劉成等[32]同樣利用紫溪山華山松無性系種子園的材料，通過SRAP分子標(biāo)記技術(shù)對其進(jìn)行遺傳多樣性分析，結(jié)果顯示：該園的華山松群體遺傳多樣性處于相對較高的水平，遺傳變異在種源內(nèi)。因此，利用該種子園進(jìn)行華山松雜交育種時，應(yīng)多進(jìn)行種源內(nèi)雜交，兼顧遺傳距離較遠(yuǎn)的種源間雜交。同時，需盡量開發(fā)較多的早期選擇標(biāo)記，以實現(xiàn)早期選擇。

由于分子標(biāo)記方法的不同，開發(fā)基因標(biāo)記的試驗結(jié)果會有一定差異。本研究所采用的分子標(biāo)記方法是SNP 標(biāo)記技術(shù)，與其他分子標(biāo)記方法相比較，SNP 標(biāo)記數(shù)量多、遺傳穩(wěn)定，且不需要凝膠電泳檢測，為試驗節(jié)省了大量的時間和成本。而利用如簡單重復(fù)間序列分子標(biāo)記(ISSR)、序列相關(guān)擴(kuò)增多態(tài)性分子標(biāo)記(SRAP)和簡單重復(fù)序列分子標(biāo)記(SSR)等方法開發(fā)華山松基因組的DNA 分子標(biāo)記難度相對較大、需要的時間較長，且標(biāo)記數(shù)量有限。趙楊等[29]利用ISSR 分子標(biāo)記共獲得13 條具有多態(tài)性的引物；劉成等[32]利用SRAP 分子標(biāo)記共篩選出15 對多態(tài)性引物組合；祝娟[33]利用SSR 標(biāo)記對華山松基因組DNA 進(jìn)行PCR 擴(kuò)增，篩選獲得12 對具有多態(tài)性的引物。而本研究利用SLAF-seq 測序技術(shù)共得到3 469 074個SNP 標(biāo)記。可見，SLAF-seq 測序技術(shù)開發(fā)分子標(biāo)記的效率較其他分子標(biāo)記技術(shù)效率較高。

3.2 華山松球果數(shù)量性狀的相關(guān)基因

本研究共檢測到與華山松球果數(shù)量顯著相關(guān)的SNP 位點12 個，并在Marker13307 位點中檢索到KOR1基因。KOR1基因是一種纖維素酶，它參與調(diào)控纖維素合成，對植物生長發(fā)育、組織脫落以及維持植物正常生理代謝有著極其重要的作用[34]。李萍[35]對新疆杏中早熟、中熟和晚熟品種的果實發(fā)育期進(jìn)行了生理生化機(jī)制研究，結(jié)果表明：伴隨著果實的發(fā)育，果實中纖維素以及半纖維素的含量先呈上升趨勢，14 d 后急劇下降，而纖維素酶在果實生長發(fā)育過程中則一直呈上升趨勢。另有學(xué)者以壺瓶棗作為試驗材料，測定其果實的果肩果皮、果肩果肉、果底果皮和果底果肉在盛花后40、55、70 和85 d 的纖維素含量變化，結(jié)果表明：壺瓶棗果實整體的纖維素含量呈先上升后下降的趨勢[36]，這與李萍[35]對新疆杏果實的研究基本一致?？梢?，纖維素在果實的生長發(fā)育中起到重要的作用，而參與調(diào)控纖維素合成的KOR1基因也在果實的生長發(fā)育中有著不可或缺的地位。KOR1是纖維素合成所必需的基因，在植物和其他生物(如細(xì)菌)質(zhì)膜—細(xì)胞壁的合成中也發(fā)揮著重要作用[37-38]。目前，已經(jīng)確定KOR1與纖維素合成酶復(fù)合物具有物理相互作用，并參與了纖維素伸長過程中甾醇連接引物的斷裂或纖維素纖維結(jié)晶度的降低，KOR1等位基因突變會導(dǎo)致纖維素合成受阻以及生長遲緩，在KOR1基因點突變的irx2 突變體中，束間纖維和木質(zhì)部的纖維素含量降低，導(dǎo)致導(dǎo)管塌陷，表明KOR1參與次生壁中纖維素的生物合成[39-41]。此外，KOR1的其他突變會導(dǎo)致細(xì)胞板形成缺陷[42]和初級細(xì)胞壁纖維素含量減少[43-44]，表明KOR1在初級細(xì)胞壁的纖維素生物合成中也起著重要作用。本研究的qPCR 結(jié)果顯示KOR1基因在華山松幼嫩球果中的表達(dá)量較高，因此認(rèn)為KOR1基因參與了幼嫩球果的生長發(fā)育。

4 結(jié)論

巍山種源的華山松較其他種源擁有較為廣泛的遺傳變異。結(jié)合華山松表型數(shù)據(jù)和SNP 標(biāo)記進(jìn)行全基因組關(guān)聯(lián)分析，得到12 個與球果數(shù)量性狀顯著關(guān)聯(lián)的SNP 標(biāo)記，并通過檢索得到1 個與植物纖維素合成相關(guān)的基因KOR1，該基因在幼嫩球果中的表達(dá)量較高，可能與球果數(shù)量性狀相關(guān)。