基于溫帶和熱帶玉米群體全基因組選擇和雜種優勢候選位點的鑒定

周玲　熊威　胡俏強　戴惠學　趙涵

摘要：基于187份種質資源材料的全基因組重測序數據開發了120583個高質量SNP變異位點，通過這些位點可以將具有不同遺傳背景的187份玉米種質劃分為兩大類群，分別為包含100份材料的溫帶亞群和包含87份材料的熱帶亞群。通過對溫帶和熱帶玉米群體選擇信號的遺傳分化分析，檢測到3664個受到選擇的位點。選取187份材料中已報道能在溫/熱亞群形成雜種優勢的135份代表性自交系，其中75份來自溫帶玉米自交系，60份為熱帶玉米系，基于兩大雜種優勢群進行雜種優勢性狀的全基因組關聯分析（GWAS），[JP2]結果鑒定出2407個雜種優勢候選位點不均勻分布在玉米10條染色體上。整合選擇信號檢測和GWAS分析結果，共識別出1153個受到選擇的雜種優勢相關位點，其中，619個位點與26個已報道的雜種優勢相關QTLs一致。功能注釋發現與候選位點緊密連鎖的324個候選基因大部分都具有功能，其中包含61個重要的轉錄因子。根據GO富集分析發現這些候選基因主要參與了很多對雜種優勢形成有貢獻的關鍵生化代謝途徑，包括氮化合物代謝、葉酸代謝、糖酵解、發育過程的負調控及轉錄調控等重要生物學途徑。

關鍵詞：溫熱帶玉米;全基因組選擇;雜種優勢;關聯分析;候選位點

中圖分類號：S513.032文獻標志碼：A

文章編號：1002-1302（2021）04-0019-07

作者簡介：周玲（1987—），女，安徽合肥人，博士，副研究員，從事玉米遺傳育種學研究。E-mail：zlingxiaoyao@163.com。

通信作者：趙涵，博士，研究員，從事玉米遺傳育種學研究。E-mail：zhaohan@jaas.ac.cn。

玉米在我國是重要的糧食、飼料、經濟作物，盡管目前培育和審定的玉米品種很多，但是雜交種親本趨同化嚴重，所用親本自交系主要集中在少數幾個骨干自交系或其衍生系，狹窄的玉米育種遺傳資源已成為制約當今玉米育種的瓶頸因素，且生產上又存在病蟲害大面積暴發的風險，擴增、改良和創新玉米種質資源已經迫在眉睫[1]。

雜種優勢利用是玉米育種和生產的根本途徑，因此在種質創新中應優先考慮創新種質的雜種優勢利用問題。根據生長環境的不同，玉米可劃分為溫帶和熱帶兩大類群。在溫帶最主要的雜優利用模式是瑞德×蘭卡斯特，在熱帶、亞熱帶地區則是Tuxpeno×ETO。溫帶種質主要表現出莖稈強韌、出籽率高、產量潛力大、農藝性狀優良等特點。熱帶、亞熱帶玉米種質則表現出抗逆性強、根系發達、葉片濃綠、持綠期長等優點[2]。由于熱帶和亞熱帶玉米種質資源在病蟲、干旱等生物和非生物脅迫耐受性方面具有廣泛的遺傳變異，[JP2]對于玉米群體改良具有重要應用價值，特別是對于溫帶玉米。因此利用熱帶、亞熱帶×溫帶雜優模式，將有利于拓寬我國溫帶玉米種質基礎，在我國玉米育種中發揮重要作用。

由于熱帶玉米在經歷人工馴化和改良過程中逐漸適應溫帶環境，因此，有必要利用選擇信號檢測法揭示在改良過程中發生選擇的位點，挖掘候選馴化基因，進一步從基因層面上探究玉米的群體改良。例如Liu等利用溫帶、熱帶和亞熱帶玉米中代表性的260個玉米自交系，結合固定系數（FST）的計算方法，發現熱帶玉米相比于溫帶玉米具有更高的遺傳多樣性和更多的等位基因位點[3]。楊宇昕等通過對溫熱帶玉米群體的高測序深度的SNP進行分析，在熱帶玉米群體和溫帶玉米群體中鑒定到204752個SNP標記，通過FST選擇信號識別到557候選馴化基因，并鑒定到多個候選基因與玉米的開花調控密切相關[4]。[JP2]因此，利用選擇信號法探究溫熱帶玉米群體的基因組變化，可以鑒定在熱帶玉米適應性改良過程中受到選擇的基因區段，并且深入挖掘更多與雜種優勢關聯的基因，對于解析雜種優勢群的遺傳機制和指導雜種優勢群利用具有重要意義。

本研究利用來源不同的187份玉米自交系的全基因組重測序信息開發高質量的SNP位點，通過這些遺傳變異位點研究自交系群體的遺傳結構，劃分雜種優勢群，通過選擇信號分析和關聯分析策略主要研究以下內容：（1）鑒定溫帶和熱帶玉米在適應環境和遺傳改良過程中，基因組上由于選擇遺留的選擇信號。（2）篩選出溫熱帶兩大雜種優勢群之間的遺傳變異分化關聯的候選位點。（3）確定受到選擇的雜種優勢關聯候選位點及其緊密連鎖的候選基因，并分析潛在選擇候選基因顯著富集途徑，為挖掘和利用控制雜種優勢性狀的位點和候選基因提供信息參考。

1材料與方法

1.1研究材料與基因組測序數據

供試材料為187份玉米自交系，種質基因組重測序數據下載于NCBI網站（www.ncbi.nlm.nih.gov）的SRA數據庫（PRJNA389800），對其進行過濾后利用BWA軟件比對到玉米B73（v3版本）基因組上并使用SAMtool對結果進行整理，最后利用GATK軟件進行SNP變異檢測，獲得變異位點后利用Plink軟件剔除缺失率大于10%和最小等位基因頻率低于0.05的SNP位點以及覆蓋率較低的樣本，采用Beagle軟件對缺失基因型進行填補后，最終獲得基因組的變異信息，SNP標記密度分布利用R軟件包CMplot繪制。

1.2雜種優勢類群劃分

利用TASSELV5.0（www.maizegenetics.net/tassel）的鄰接算法（neighbor-joiningmethod，簡稱NJ）[5]計算187份玉米自交系之間的遺傳距離構建遺傳樹，并進行雜種優勢類群劃分。

1.3FST選擇信號檢測

基于不同群體的遺傳多樣性數據，可以利用FST進行溫熱帶玉米間選擇信號的檢測，其計算原理是依據染色體等位基因頻率變化。在中性進化條件下，F統計量的大小主要取決于遺傳漂變和遷移等因素的影響，若群體中一個等位基因因為對于特定環境的適合度較高而經歷適應性選擇，那么其頻率的升高會增大種群分化水平，反映在F統計量上就是有較高的FST值（0≤FST≤1，FST<0.05表示種群間沒有遺傳分化，FST>0.25則表示種群間分化程度非常高）[6]。為了降低因單位點SNP掃描受到遺傳漂變等因素影響而造成的假陽性，筆者利用VCFtools軟件計算滑動窗口10kb內群體間的FST值來增加選擇信號的靈敏度[7]，并運用R包CMplot繪制全基因組水平上的FST結果圖。為了鑒定FST值的受選擇位點，選擇FST值高于0.25作為顯著閾值線，高于閾值線的SNP位點定義為“受選擇位點”。

1.4全基因組關聯分析（GWAS）

運用R軟件GAPIT軟件包中的壓縮混合線性模型（compressedmixedlinearmodel，簡稱CMLM）程序[8]，計算前5個主成分控制群體結構，并根據VanRaden描述的方法[9]計算親緣關系K矩陣。在考慮研究材料群體結構和親緣關系的情況下，對溫熱帶兩大類群進行全基因組關聯分析，溫帶的雜種優勢表型值記作1，熱帶記作0。采用Bonferroni方法矯正P值，鑒定溫熱帶群體之間的特異性位點。

1.5功能注釋和GO富集分析

對基于FST和GWAS共識別的特異位點進行全基因組掃描，利用SnpEffv4.3軟件對變異信息進行功能注釋[10]，獲得與其緊密連鎖的候選基因，經MaizeGDB（http：//www.maizegdb.org/blast.php）[JP2]數據庫比對尋找有同源序列的功能基因。為了更深刻解析選擇清除區域候選基因的生物學功能，利用在線平臺AgriGO的單一富集分析（SEA）功能（http：//systemsbiology.cau.edu.cn/agriGOv2/index.php）[11]對候選基因進行生物學過程（biologicalprocess）富集分析，篩選P值<005且錯誤發現率（FDR）<005的GO富集術語作為選擇候選基因。

2結果與分析

2.1基因型鑒定和進化樹構建

本研究利用187份種質資源材料的高通量重測序數據開發了120583個SNP位點，這些變異位點基本均勻分布于10條染色體上（圖1）。采用TASSEL5.0軟件的NJ法構建了系統發育樹，結果（圖2）顯示，187份玉米材料被劃分為2個不同的類群上：100份材料被聚類為溫帶玉米亞群，另外87份材料被聚類為熱帶玉米亞群。

2.2基于全基因組選擇信號的檢測和關聯分析

基于全基因組檢測到的SNP標記利用VCFtools軟件計算溫熱帶亞群體之間的FST值，結果表明，全基因組水平上FST值高于閾值線0.25的受選擇位[CM（21][KG*8]點共有[KG*8]3[KG*3]664[KG*8]個（圖3-a），占總變異位點數的3%。其中染色體4受到選擇的顯著性位點最多（913個），染色體6（188個）和染色體7（189個）最少。第4染色體46380000～46390000區間內含有最高的FST值，其值為0.66。

選取187份材料中具有廣泛代表性的135份玉米自交系（表1），其中75份是選自育種程序的溫帶玉米自交系，代表經歷顯著改良并適應溫帶長日照環境的溫帶玉米（如B73、Mo17等），而60份主要來源于國際玉米小麥改良中心（CIMMYT）的熱帶、亞熱帶玉米自交系（如CML103、NC350等），代表具有豐富的遺傳多樣性、改良程度較低、更適應玉米馴化原始環境的熱帶玉米。在控制群體結構和材料間親緣關系條件下，基于溫熱帶兩大雜種優勢群體利用GAPIT軟件CMLM方法對其雜種優勢性狀的GWAS，采用Bonferroni方法矯正P值，顯著水平閾值為α=0.05/120583=4.15×10-7，在-lgP>683的水平下，共檢測到2407個特異位點與雜種優勢群劃分顯著關聯（圖3），其中顯著的P值的范圍是6.84～20.67，解釋了4.31%～11.97%的表型變異。在10條染色體的分布依次為214、180、568、576、41、51、97、274、122、284個，其中在染色體3（568個）和染色體4（576個）上面存在顯著富集。

整合全基因組選擇信號的檢測和GWAS分析的結果，發現受到選擇的雜種優勢相關SNPs共有1153個，其中在染色體3（233個）、染色體4（433個）和染色體10（201個）上有顯著富集，在染色體1、2、5、6、7、8、9染色體上分別包含34、37、21、17、44、114、19個顯著關聯SNPs。與已報道的雜種優勢相關QTL進行比較分析[12-15]，發現619個位點存在于24個已報道的QTL區段（表2），其中染色體3中4個QTL區段覆蓋了140個候選位點，染色體4中3個QTL區段覆蓋了146個位點，[JP2]染色體10中3個QTL區段覆蓋了178個位點，表明本研究識別的候選位點與已報道的雜種優勢QTL結果高度一致。

2.3候選基因的富集分析和功能注釋

利用Perl腳本和SnpEff4.3軟件對受到選擇的雜種優勢關聯的顯著SNP標記與玉米B73參考基因組進行比對，并進行候選基因預測及功能注釋，最終筆者獲得了324個候選基因，其中包含61個轉錄調控功能的基因（表3）。利用AgriGO在線軟件SEA對候選基因進行富集分析發現候選基因主要富集在55個GO富集術語生物過程（圖4）。這些GO富集術語主要涉及到氮化合物代謝（GO：0006807）、葉酸代謝（GO：0046655）、糖酵解（GO：0006096）、發育過程的負調控（GO：0051093）和轉錄調控（GO：0006355）等重要生物學途徑（圖4）。

3結論與討論

對玉米種質遺傳多樣性深入認識是合理利用種質的前提，隨著分子生物學的發展，分子標記已成為作物遺傳育種研究的重要手段，同時為評價玉米的遺傳變異提供了方便、快捷的研究方法。利用分子標記的方法對玉米種質進行類群劃分，方法易行、結果可靠，已經大規模用于玉米育種研究中。本研究以B73的參考基因組和187份玉米自交系的基因組信息為背景開發了120583個高質量的SNP分子標記均勻分布于10條染色體上，基于這些標記將187份自交系劃分為2個雜種優勢類群，分別為溫帶和熱帶兩大類群。合理的雜種優勢群，是把遺傳關系較近的自交系劃分為同一雜種優勢群，在雜交種選育過程中避免遺傳關系較近的自交系之間試組配工作，可以有效減少自交系之間的組配工作量，大大提高育種效率。

解析不同玉米類群遺傳差異，鑒定與雜種優勢相關的基因有助于深入了解玉米分化過程，同時為玉米種質的擴增和改良提供重要遺傳資源。在本研究中，筆者結合選擇信號分析和GWAS方法共鑒定出1153個受到選擇的雜種優勢相關候選位點，前人已報道的26個雜種優勢相關QTL區段覆蓋了本研究獲得的619個候選位點，表明本研究的結果與已報道的雜種優勢QTL結果高度一致，這為溫熱帶種質資源的雜交和改良奠定了基礎。功能注釋發現與候選位點緊密連鎖的324個候選基因中大部分都是具有功能的基因，利用GO富集分析發現，這些候選基因主要參與了很多對雜種優勢形成有貢獻的關鍵生化代謝途徑，包括氮化合物代謝、葉酸代謝、糖酵解、發育過程的負調控及轉錄調控等重要生物學途徑。

此外，這些候選基因中還包括了61個轉錄因子，在馴化過程中，植物對各種環境、組織和發育信號做出反應，需要各種功能基因的表達，從而進行精確的調控。而植物感受外界環境如干旱、病害、鹽堿等體內細胞發育等信號時，需要一系列信號傳遞激發轉錄因子，從而啟動功能基因的轉錄表達，最后通過基因產物的作用對外界信號在生理生化等方面的變化做出適當的調節反應[16]?？梢?，植物的轉錄因子通過其功能域與DNA及其他蛋白質間的相互作用，可以調節植物的生長發育及環境脅迫應答等過程[17-18]。因此，這61個轉錄因子將成為研究玉米馴化的候選基因。

參考文獻：

[1]戴景瑞，鄂立柱.我國玉米育種科技創新問題的幾點思考[J].玉米科學，2010，18（1）：1-5.

[2]TallurySP，GoodmanMM.Experimentalevaluationofthepotentialoftropicalgermplasmfortemperatemaizeimprovement[J].TheoreticalandAppliedGenetics，1999，98（1）：54-61.

[3]LiuK，GoodmanM，MuseS，etal.GeneticstructureanddiversityamongmaizeinbredlinesasinferredfromDNAmicrosatellites[J].Genetics，2003，165（4）：2117-2128.

[4]楊宇昕，鄒棖.基于溫帶和熱帶玉米群體全基因組FST和XP-EHH的選擇信號檢測[J].中國農業科學，2019，52（4）：579-590.

[5]SaitouN，NeiM.Theneighbor-joiningmethod：anewmethodforreconstructingphylogenetictrees[J].MolecularBiologyandEvolution，1987，4（4）：406-425.

[6]WrightS.Evolutionandgeneticsofpopulations[M].Chicago：UniversityofChicagoPress，1978.

[7]MaY，DingX，QanbariS，etal.Propertiesofdifferentselectionsignaturestatisticsandanewstrategyforcombiningthem[J].Heredity，2015，115（5）：426-436.

[8]ZhangZ，ErsozE，LaiCQ，etal.Mixedlinearmodelapproachadaptedforgenome-wideassociationstudies[J].NatureGenetics，2010，42（4）：355-360.

[9]VanRadenPM.Efficientmethodstocomputegenomicpredictions[J].JournalofDairyScience，2008，91（11）：4414-4423.

[10]CingolaniP，PlattsA，WangLL，etal.Aprogramforannotatingand[HT][HT][LM]

[KH*4D]

[HT8.][KG2]predictingtheeffectsofsinglenucleotidepolymorphisms，SnpEff：SNPsinthegenomeofDrosophilamelanogasterstrainw1118;iso-2;iso-3[J].Fly，2012，6（2）：80-92.

[11]TianT，LiuY，YanH，etal.AgriGOv2.0：aGOanalysistoolkitfortheagriculturalcommunity，2017update[J].NucleicAcidsResearch，2017，45（1）：122-129.

[12]StuberCW，LincolnSE，WolffDW，etal.Identificationofgeneticfactorscontributingtoheterosisinahybridfromtwoelitemaizeinbredlinesusingmolecularmarkers[J].Genetics，1992，132（3）：823-839.

[13]LuH，Romero-SeversonJ，BernardoR.Geneticbasisofheterosisexploredbysimplesequencerepeatmarkersinarandom-matedmaizepopulation[J].TheoreticalandAppliedGenetics，2003，107（3）：494-502.

[14]FrascaroliE，CanèMA，LandiP，etal.Classicalgeneticandquantitativetraitlocianalysesofheterosisinamaizehybridbetweentwoeliteinbredlines[J].Genetics，2007，176（1）：625-644.

[15][JP2]ThiemannA，FuJ，SeifertF，etal.Genome-widemeta-analysisofmaizeheterosisrevealsthepotentialroleofadditivegeneexpressionatpericentromericloci[J].BMCPlantBiology，2014，14：88.

[16]LedentV，VervoortM.Thebasichelix-loop-helixproteinfamily：comparativegenomicsandphylogeneticanalysis[J].GenomeResearch，2001，11（5）：754-770.

[17]VaragonaMJ，SchmidtRJ，RaikhelNV.Nuclearlocalizationsignal（s）requiredfornucleartargetingofthemaizeregulatoryproteinOpaque-2[J].ThePlantCell，1992，4（10）：1213-1227.

[18]劉欣，李云.轉錄因子與植物抗逆性研究進展[J].中國農學通報，2006，22（4）：61-65.