棉花苗期耐鹽相關性狀QTL元分析

關鍵詞：棉花；耐鹽；元分析；一致性QTL；候選基因

doi：10.13304/j.nykjdb.2022.0283

中圖分類號：S562 文獻標志碼：A 文章編號：10080864（2023）12002609

棉花（Gossypium）是重要的經濟和纖維作物，同時也是鹽堿地的先鋒作物[1]，但隨著土壤鹽漬化愈發嚴重，棉花生產也受到一定危害。2021年，聯合國糧食及農業組織（Food and AgricultureOrganization of the United Nations ， FAO）發布，全球表層土壤（0—30 cm）鹽漬化面積逾4.24 億hm²，占全球土壤總面積的3%；底層土壤（30—100 cm）鹽漬化面積超8.33億hm²，占全球土壤總面積的6%，各大洲均有20%～50%不等的灌溉土壤鹽漬化[2]。新疆是我國棉花主產區，也是土壤鹽漬化嚴重地區，其鹽漬土壤面積為約占新疆總面積的6.61%，其中耕地土壤鹽漬化面積達12.6萬hm²，占耕地總面積的31%[3]。鹽害使棉花的營養和生殖生長均受到抑制，導致棉株減產甚至死亡[4]。

在棉花生長周期中，棉株對鹽分的敏感度和耐受性呈動態變化。苗期對鹽分的變化極為敏感且耐受性最低，在高鹽環境下，棉苗生長勢減弱甚至死亡。彭振等[5]研究表明，棉花幼苗在三葉期耐鹽性最差，之后耐鹽性隨植株生長發育不斷提高。Ashraf等[6]從組織、細胞、分子水平研究顯示耐鹽棉種在苗期具有更強大的Na+區室化能力（Na+ compartmentalization），可以減輕Na+過多累積對細胞質的毒害作用。

耐鹽性是典型的由多基因控制的復雜數量性狀。目前，數量性狀位點（quantitative trait loci，QTL）定位常用于對復雜數量性狀的研究，但是對棉花耐鹽性相關性狀的研究相對較少。劉晨晨[7]利用農大棉13號×農大601構建了包含588個株系的重組自交系（recombinant inbred lines， RIL）群體，采用復合區間作圖法（complex intervalmapping， CIM）共定位出14個控制相對發芽率的QTLs，分布在10條染色體上。Igbal等[8]以魯棉研28×新陸早24構建了RIL群體，采用CIM 方法共定位到與16個耐鹽性狀相關的262個QTLs。朱協飛等[9]以148份導入系為研究對象，共鑒定出23個與植株鮮物質量相關的QTLs，分布在9條染色體上，且單個QTL可解釋1.11%～3.18%的表型變異。但是，由于群體遺傳背景、標記類型等不同導致定位到控制棉花耐鹽性的QTLs數目、遺傳效應及其在染色體上的位置等存在較大差異，難以獲得真實、有效的QTLs，無法直接應用于育種實踐。因此，通過整合已報道的QTLs來提高棉花耐鹽性定位精確度具有重要意義。

元分析是一種整合已有研究的統計分析方法，通過增加樣本量使所獲得的結果更加穩健[10]。目前，該技術已經在水稻、小麥、玉米等作物中成功運用。通過對不同研究中的QTLs進行比較、整合，達到縮小置信區間、提高QTL定位精確度的目的。Arcade 等[11] 開發了用于元分析的軟件BioMercator，通過整合遺傳圖譜，并將相關QTLs都映射到一張遺傳圖譜上，再經過Overview等分析，預測一致性QTLs（meta quantitative trait loci，MQTL）。馮世超等[12]以Cornell 2001作為參考圖譜，整合了15 個作圖群體的分子圖譜，映射了263個與水稻耐鹽性相關的QTLs，共挖掘到14個MQTLs。王鵬等[10]收集了113個小麥粒長和86個粒寬QTLs，共挖掘出29個MQTLs，分布于8條染色體上。王曉麗等[13]收集了441個玉米產量及其構成因子QTLs，涉及10個產量性狀，并整合22張遺傳圖譜，共挖掘出5個MQTLs，篩選出7個產量候選基因。

使用QTL 元分析對棉花數量性狀進行QTL整合的研究較少，楊鑫雷[14]通過整合由中棉所8號（陸地棉）和Pima90-53（海島棉）構建的BC1、BC1F2群體定位到的92個纖維品質性狀QTLs，構建了一張包含其中63個QTLs、599個標記、覆蓋度達3 571.90 cM的整合圖譜，共在15條染色體上挖掘出30個MQTLs；Said等[15]基于棉花標記數據庫，收集來自42個報道中1 223個與棉花產量、形態建成、耐旱性、抗病性及纖維品質相關的QTLs，挖掘出8個與纖維品質、4個纖維顏色、1個葉片形態、5 個黃萎病抗性和4 個線蟲抗性相關的MQTLs；Zhang 等[16]收集401 個棉花抗性相關的QTLs，挖掘出13 個黃萎病抗性、8 個根結線蟲（Meloidogyne incognita）抗性和3 個腎形線蟲（Rotylenchulus reniformis）抗性相關的MQTLs。然而，暫無對棉花苗期耐鹽性相關性狀QTL元分析的相關報道。本研究將已報道的棉花苗期耐鹽相關性狀的QTLs進行統計整理，選取棉花高密度遺傳圖譜作為參考圖譜，利用BioMercator V4.2.3軟件，結合原始圖譜構建整合圖譜，映射原始QTLs至整合圖譜以構建一致性圖譜，并預測MQTL及其區間內與棉花耐鹽性相關的候選基因，以期為深入研究棉花苗期耐鹽機制提供理論依據。

1 材料與方法

1.1 材料

以棉花QTL 數據庫（http：//www.cottonqtldb.org）和已發表文獻中控制棉花苗期耐鹽性的QTLs信息為研究對象。

1.2 方法

1.2.1 耐鹽相關性狀QTLs數據收集 收集已報道的棉花苗期耐鹽相關性狀QTL定位信息，統計QTLs名稱、所在染色體及位置（最大可能性位置及置信區間）、表型變異率等信息。將在特定環境、群體、目標性狀條件下的QTL分析視為1次實驗，利用多點數據均值進行QTL分析亦視為1次實驗。當QTLs置信區間未知時，可通過下列公式求得95%置信區間（confidence interval，CI）[17]。

CI=163/（N×R2） （1）

CI=530/（N×R2） （2）

式中，Ｎ為群體大小；R²為貢獻率，%。

RIL群體采用公式（1）；F₂和回交（backcross，BC）群體使用公式（2）。

1.2.2 耐鹽相關性狀QTLs的映射 以Diouf等[18]發布的遺傳圖譜作為參考圖譜，使用BioMercatorV4.2.3軟件，將Oluoch等[19]和Abdelraheem等[20]發布遺傳圖譜中的分子標記映射至參考圖譜。確定收集到QTLs所涉及的染色體，將目標QTLs依據標記點數據，利用齊序函數（homothetic function）映射到整合的高密度遺傳圖譜上[21]。

1.2.3 耐鹽相關性狀QTL 元分析 利用BioMercator V4.2.3 軟件，通過染色體步移方法（chromosome walking method）對在同一染色體相同或部分重疊座位且具有獨立來源的QTLs計算出1 個MQTL[22]。依據赤池信息量準則值（akaikeinformation criterion，AIC），選取每個MQTL具有的5 個模型中AIC 值最小的作為真實QTL（real-QTLs）模型，并采用高斯定理最大似然比對MQTL存在的最大可能性位置和置信區間進行估算[23]。

2 結果與分析

2.1 QTL 信息收集和一致性圖譜整合

從已報道的研究中收集到194個與棉花苗期耐鹽性相關的QTLs，分別來自于CRI-12×AD₃-00、TM-1×NM24016 和CCRI35×Nan Dan（NH）共3 個不同遺傳背景的作圖群體，涉及發芽率、主根長、莖長、株高、地上部鮮重、地上部干重、根鮮重、根干重、葉鮮重、葉干重、鮮重、飽和葉重共12個形態指標及丙二醛含量、電導率、相對含水量、光合速率、葉綠素含量、氣孔導度、蒸騰作用共7個生理指標（表1）。

分析耐鹽相關QTLs在染色體上的分布（圖1）發現，D02染色體上最多，有25個QTLs，多數與形態性狀相關；A04、A05、A07和D09染色體上收集到的QTLs最少，均為2個。形態性狀QTLs多分布于D 染色體組（62.2%）；生理性狀QTLs 在A（53.1%）、D（46.9%）染色體組分布比較均等；出苗率QTLs 均分布于A 染色體組；根長、莖長、根干重、光合速率QTLs均定位于D染色體組。

參考圖譜總長4 768.10 cM，涉及5 718個SNPs標記[18]。整合分別具有1 295 個覆蓋3 328.24 cM的SSR標記圖譜[19]和具有499個標記（包含SSR、GBS-SNP、RGA-AFLP）總遺傳距離為2 221.28 cM的遺傳圖譜[20]。結果顯示，共整合5 322個分子標記于26 個連鎖群上，相鄰標記間的平均距離為1.15 cM；相鄰標記間最大間隔為105.73 cM，位于D02 染色體。分子標記映射最多的染色體為A02，共整合623 個標記，標記間平均距離為0.65 cM；分子標記映射最少的染色體為D09，共43 個標記，標記間平均距離為3.41 cM。共有182個QTLs被映射至整合圖譜，分布于26個連鎖群（部分連鎖群如圖1所示）。

2.2 苗期耐鹽相關性狀QTL 元分析結果

QTL元分析結果（表2）顯示，得到11個與棉花苗期耐鹽相關的MQTLs，每個MQTL至少包含3個原始QTLs，分布于A03、A06、A11、A12、A13、D01、D03、D06、D07和D08共10條染色體上，其中D01染色體上包含2個MQTLs；最小置信區間為0.92 cM，位于D07染色體上，最大為12.69 cM，位于D06染色體上。

2.3 候選基因預測結果

以TM-1 基因組[24] 作為參考基因組，依據MQTL在基因組物理圖譜上的位置，預測棉花苗期耐鹽相關性狀MQTL區間內的候選基因。因為較小的AIC值和圖距有利于提高定位的精確度，所以選取元分析結果中AIC 值和區間均較小的MQTL3作為目標區間，利用Jbrowse檢索，共獲得85個基因。使用模式作物擬南芥基因序列進行同源比對，利用GO（gene ontology）、KEGG（kyotoencyclopedia of genes and genomes）等數據庫注釋，結合表達量數據等方法，篩選出14 個候選基因（表3）。

由表3和圖2可知，Ghi_A11G02546 基因編碼的HD-ZIP蛋白和Ghi_A11G02636 編碼的MYB蛋白均為轉錄因子，且Ghi_A11G02636 基因在鹽脅迫1 h時高表達。Ghi_A11G02456、Ghi_A11G02496、Ghi_A11G02641、Ghi_A11G02631 均編碼蛋白激酶，其中Ghi_A11G02631 基因在鹽脅迫1 h 時表達量達對照組的2 倍。Ghi_A11G02386 編碼蛋白磷酸酶。Ghi_A11G02401、Ghi_A11G02611 和Ghi_A11G02416 均為編碼膜蛋白的基因，其中Ghi_A11G02401 及Ghi_A11G02611 基因編碼離子轉運蛋白；Ghi_A11G02416 基因編碼水通道蛋白。Ghi_A11G02376 基因參與生物合成途徑，編碼甜菜堿醛脫氫酶。這些基因通過編碼核糖核酸酶、MAPK激酶家族蛋白、H⁺-ATPase酶等，參與棉花苗期對鹽脅迫的響應，如信號轉導、離子及水分運輸、甜菜堿生物合成等，進而影響棉花幼苗抵御鹽脅迫的能力。

利用Jbrowse檢索，以候選基因上游2 000 bp核苷酸序列為目標序列，使用Plant CARE進行啟動子順式作用元件預測，結果（圖3）顯示，候選基因啟動子區域均具有鹽脅迫或水分脅迫誘導的順式作用元件。由此表明，候選基因可能參與鹽脅迫下的調控過程。

3 討論

遺傳圖譜的完善和加密提高了QTL 定位的可靠度和精準度，棉花基因組測序工作的推進促進了候選基因的篩選。本研究選取通過GBS高通量測序分型技術構建的SNP標記遺傳圖譜[18]作為參考，其總長為4 768.10 cM，涉及5 718 個標記。此圖譜密度較大，但標記種類較為單一，僅具有GBS-SNP類型標記。因此，本研究整合Oluoch等[19]和Abdelraheem等[20]發布的遺傳圖譜，整合后的遺傳圖譜標記類型增多，包含SSR、GBS-SNP、RGA-AFLP標記，且標記密度增大，但仍有部分目標性狀QTLs不能映射成功，這可能是由于原始圖譜標記和整合后圖譜標記的一致性較差。對于不能成功映射的QTLs，本研究將其剔除。最終本研究整合的遺傳圖譜具有高密度、高標記多樣性的特點，為高效QTLs映射、縮小側翼標記區間和精細定位奠定了基礎。

QTL定位的精確度和參考基因組的完善是篩選和克隆目標性狀重要候選基因的良好前提，利用QTL元分析可以有效地縮短目標性狀的置信區間，提升基因篩選和鑒定的效率。隨著越來越多目標性狀被精細定位，使用生物信息學、比較基因組學等方法，結合完備的基因注釋數據庫、模式作物基因測序和注釋結果、基因表達量數據庫等工具，通過同源比對，可以預測基因功能，更加可靠地篩選候選基因。本研究在調控棉花苗期耐鹽相關性狀的MQTL 區間內初步挖掘到了14個候選基因。其中，Ghi_A11G02456、Ghi_A11G02496、Ghi_A11G02641 基因在擬南芥中的同源基因均參與MAPK信號級聯通路，該通路在調控細胞分裂和抵御鹽害等非生物脅迫中發揮重要作用。研究表明，Ghi_A11G02456 基因編碼的蛋白質位于細胞器膜上，作為滲透傳感器響應鹽、干旱脅迫及脫落酸（abscisic acid， ABA）刺激，并將壓力信號傳遞至下游MAPK 級聯反應通路[25-27]；Ghi_A11G02496和Ghi_A11G02641 基因產物均屬于MAPK激酶家族，在MAPK 信號轉導通路中發揮作用[2829]；Ghi_A11G02291 是核糖核酸酶基因，在鹽脅迫期間充當 mRNA 脫帽的正調節劑，可作為調節因子參與鹽逆境下mRNA分解代謝過程[30]；Ghi_A11G02376基因可受鹽脅迫誘導編碼一種甜菜堿醛脫氫酶參與甘氨酸甜菜堿的合成，甘氨酸甜菜堿可作為滲透調節物質參與維持細胞滲透平衡[31]；Ghi_A11G02286是編碼鈣離子結合蛋白基因，可能參與Ca2+介導的信號傳導途徑，參與滲透調節，起到維持離子平衡的作用；Ghi_A11G02416 是在根部特異性表達、編碼水通道蛋白的基因，對維持根部細胞滲透勢起到重要作用；Ghi_A11G02401 編碼的單價陽離子轉運蛋白組成位于液泡膜上的H⁺-ATPase酶，在鹽逆境下促進離子和代謝物的轉運[32]；Ghi_A11G02631 基因編碼SNF1 蛋白相關蛋白激酶，受離子或非離子滲透脅迫激活，參與信號轉導過程[33]。

本研究收集了來自3個作圖群體、19個苗期耐鹽相關性狀、194個QTLs數據，構建了棉花苗期耐鹽相關性狀的一致性圖譜。通過QTL元分析共得到11個苗期耐鹽性相關MQTL位點，置信區間最短縮小至0.92 cM，降低了QTL定位分析誤差。挖掘到的MQTLs中涉及到同時包含生理指標和形態指標的原始QTLs，說明在棉花苗期耐鹽脅迫響應過程中，其生理調節和生長發育過程中發揮作用的部分基因分布于染色體上同一物理區域。在A11染色體上MQTLs區域進行候選基因預測，共挖掘到14個與棉花苗期耐鹽相關的基因，涉及鹽逆境下MAPK級聯通路等信號轉導途徑、水通道蛋白等滲透調節過程及質子泵等離子平衡活動。在鹽脅迫1、3、6 h均有不同基因高表達，說明此MQTL區間內基因可能在鹽脅迫不同時間段內起相應作用，為調控棉花苗期耐鹽性的重要位點。本研究為棉花苗期耐鹽相關性狀QTL精細定位和相關基因克隆鑒定提供理論依據，對利用分子標記輔助輔助棉花耐鹽性育種具有重要意義。

（責任編輯：張冬玲）