楊宇琭，羅韶凡，郄倩茹，張雅坤，段明，李亞軍，孫蓉，韓淵懷，李旭凱*，馬芳芳*

（1.山西農(nóng)業(yè)大學 農(nóng)學院，山西 太谷 030801；2.山西農(nóng)業(yè)大學 生命科學學院，山西太谷 030801；3.山西農(nóng)業(yè)大學 實驗教學中心，山西 太谷 030801）

谷子（Setaria italic），禾本科狗尾草屬，古稱稷、粟、粱禾，為中國古代的五谷之首。其具有耐干旱、耐貧瘠、高光效等特點，是典型的環(huán)境友好作物，完全有可能重新成為主栽作物和主要糧食［1］。但是谷子產(chǎn)量沒有較大的突破，經(jīng)濟效益低，種植面積增長幅度小［2］。影響作物產(chǎn)量的因素之一是光合作用，而葉綠體、碳代謝途徑、以及多種酶類等都會影響光合速率［3-4］。其中葉綠素是光合作用的主要色素之一。在高等植物的葉綠體中，形成葉綠素a 和葉綠素b 的過程中涉及到多種酶、多個基因［5-6］。葉綠素含量減少或者過多都會影響作物的生長發(fā)育，有些甚至會導致植株死亡。研究表明葉鞘中合成的光合產(chǎn)物主要運輸?shù)阶魑锼氲炔课唬瑢Ξa(chǎn)量約有15%的貢獻［7］。研究人員通過研究不同施氮條件下水稻葉鞘中與衰老相關(guān)的非結(jié)構(gòu)性碳水化合物（non-structural carbohydratem，NSC），發(fā)現(xiàn)葉鞘中NSC 的積累與轉(zhuǎn)運對水稻增產(chǎn)穩(wěn)產(chǎn)具有重要意義［8］。因此，研究控制谷子葉鞘顏色和葉枕顏色的遺傳機制作為挖掘植株顏色的切入點，對進一步提高谷子光合效率，提高谷子產(chǎn)量具有現(xiàn)實意義。

顏色性狀是一種重要的標記性狀，早在1980年就有人利用葉片顏色篩選谷子的有性雜交種［9］。一些研究人員按照葉鞘顏色對谷子品種類型進行了初步區(qū)分［10-11］。葉鞘顏色作為形態(tài)標記在水稻和小麥育種中的應(yīng)用更廣［12-13］。在高粱中結(jié)合葉鞘顏色開發(fā)了SRR（Simple Sequence Repeats）標記用于劃分類群［14］。另外，谷子葉鞘顏色和葉枕顏色相對豐富，因此研究控制谷子葉鞘顏色和葉枕顏色的基因，對開發(fā)谷子形態(tài)標記和篩選谷子品種具有重要意義。前期，刁現(xiàn)民等對916 份谷子種質(zhì)資源進行了重測序，對谷子的多種性狀進行了統(tǒng)計，包括葉鞘顏色和葉枕顏色，為本研究奠定了基礎(chǔ)［15］。

近年來，高通量測序技術(shù)快速發(fā)展，轉(zhuǎn)錄組測序在植物研究上也逐漸被普及［16］，與此同時，生物信息分析工具和方法也在不斷的發(fā)展與更新［17］。全基因組關(guān)聯(lián)分析（Genome-wide association studies，GWAS）已在作物上廣泛應(yīng)用，但傳統(tǒng)的GWAS 由于連鎖不平衡的特點，經(jīng)常會出現(xiàn)數(shù)十個甚至上百個基因座與性狀相關(guān)聯(lián)，而這些基因座中往往包含數(shù)百個SNP 位點，其中許多SNP 位點距離真正導致表型變異的位點還比較遠，這些多余的信息會干擾確定候選基因［18］。另外很多物種的基因注釋信息并不完善，這也成為鑒定和確認功能基因的障礙。因此本研究使用了美國明尼蘇達大學的研究團隊開發(fā)的Camoco（Co-analysis of molecular components）的算法，其原理是將全基因組關(guān)聯(lián)分析（GWAS）和共表達網(wǎng)絡(luò)整合起來，通過識別GWAS 篩選出的顯著SNP 位點附近的基因并對判斷優(yōu)先級，確定共表達較強的基因，將結(jié)果與隨機網(wǎng)絡(luò)進行比較，把P值分配給候選基因，篩選出高置信度的候選基因。已有研究使用此方法對玉米17個性狀進行分析，用3個基因表達數(shù)據(jù)集和兩種共表達網(wǎng)絡(luò)驗證此方法，并對篩選出的高置信候選基因使用突變體進行了功能驗證［19-20］。

1 材料與方法

1.1 試驗材料

GWAS 數(shù)據(jù)：刁現(xiàn)民研究員團隊916 份谷子種質(zhì)資源重測序數(shù)據(jù)［15］。

轉(zhuǎn)錄組數(shù)據(jù)：名優(yōu)品種晉谷21（JG21）種植于山西農(nóng)業(yè)大學試驗田，在其生長周期中，對18個組織時期取樣，每個樣設(shè)置3個重復，進行轉(zhuǎn)錄組學分析。具體組織時期有：JG21-發(fā)芽種子-3 天、JG21-幼苗-2 周、JG21-植株-兩葉一心、JG21-未成熟小穗-S2、JG21-未成熟小穗-S4、JG21-旗葉-灌漿期S3、JG21-旗葉葉鞘-灌漿期S3、JG21-倒4葉-灌漿期S3、JG21-倒4 葉葉鞘-灌漿期S3、JG21-根-灌漿期、JG21-莖桿-灌漿期S3、JG21-倒2、3 葉-30 天、JG21-倒2、3 葉-抽穗、JG21-未成熟種子-S1、JG21-未成熟種子-S2、JG21-未成熟種子-S3、JG21-未成熟種子-S4、JG21-未成熟種子-S5。（S1：籽粒顏色為鮮綠色，質(zhì)地為乳狀；S2：籽粒顏色為暗綠色，質(zhì)地為乳狀；S3：籽粒顏色為黃綠色，質(zhì)地中出現(xiàn)粉狀物；S4：籽粒顏色為暗黃色，質(zhì)地為固態(tài)）。

1.2 試驗方法

1.2.1 GWAS 數(shù)據(jù)分析

以豫谷1號基因組為參考，依據(jù)上述916 份谷子種質(zhì)資源GWAS 結(jié)果選取候選SNP 位點上下游500 bp 的序列，通過Blast 比對到谷子高質(zhì)量基因組xiaomi上［21］，定位到新基因組位置。

1.2.2 轉(zhuǎn)錄組數(shù)據(jù)分析

基于Illumina 技術(shù)得到的雙端測序數(shù)據(jù)，利用FastQ［22］軟件對數(shù)據(jù)進行質(zhì)控，并使用Trimmonatic［23］軟件依據(jù)堿基質(zhì)量對fastq 文件進行修剪，過濾低質(zhì)量的reads，得到clean data。利用Hisat2［24］軟件將上一步數(shù)據(jù)比對到xiaomi基因組［21］。用featureCounts［25］獲得數(shù)據(jù)的reads 計數(shù)，之后計算TPM（transcripts per million）值。

1.2.3 Camoco 分析及流程

利用Camoco 將基因表達量數(shù)據(jù)和GWAS 分析獲得的SNP 位點信息結(jié)合［19］，定位到基因。具體參數(shù)設(shè)定可參谷子的設(shè)定，指定50 kb（SNP 位點向上下50 kb）范圍和1個側(cè)翼基因，以及范圍外的1個側(cè)翼基因。如果指定的范圍不覆蓋任何基因，側(cè)翼基因允許使用最近的基因。candidatewindow-size 自行根據(jù)所研究物質(zhì)的連鎖不平衡（LD）范圍改變。

1.2.4 加權(quán)共表達網(wǎng)絡(luò)分析

使用R 軟件（R version 4.0.2）中的WGCNA（R version 1.6.6）包［20］，構(gòu)建加權(quán)基因共表達網(wǎng)絡(luò)并劃分相關(guān)模塊。采用一步法構(gòu)建共表達矩陣，使用WGCNA 包中的pickSoftThreashold 計算權(quán)重值，根據(jù)表達量的相關(guān)性構(gòu)建基因聚類樹。設(shè)置基因模塊中最少基因數(shù)為30，按照混合動態(tài)剪切的標準劃分模塊，并計算每個模塊的特征向量值，對模塊進行聚類分析，將距離較近的模塊合并為新的模塊。

1.2.5 模塊的GO 富集分析

提取模塊中的基因，利用TBtools 進行GO 富集分析［26］，確定候選基因的主要生物學功能。使用Cytoscape（version 3.6.1）對模塊中的基因進行可視化［27］。

1.2.6 谷子葉鞘顏色和葉枕顏色候選基因的表達模式分析

在山西農(nóng)業(yè)大學谷子多組學數(shù)據(jù)庫（http://sky. sxau. edu. cn/MDSi. htm）中輸入候選基因的信息得到候選基因在谷子各組織時期的表達模式［21］。

2 結(jié)果與分析

2.1 基于GWAS 定位谷子葉鞘顏色和葉枕顏色相關(guān)基因

通過篩選GWAS 結(jié)果，發(fā)現(xiàn)葉鞘顏色和葉枕顏色與多個種植地區(qū)GWAS 結(jié)果關(guān)聯(lián)位點顯著且一致。2種性狀在不同環(huán)境下的GWAS 曼哈頓圖［15］結(jié)果顯示：在5種環(huán)境下谷子葉鞘顏色和葉枕顏色的GWAS 結(jié)果中，位于7號染色體和4號染色體上存在顯著大于閾值的SNP 位點（圖1）。篩選鑒定出與谷子葉鞘顏色和葉枕顏色顯著的SNP 位點，初步獲得顏色性狀相關(guān)的候選基因（增強出版附件1）。

圖1 谷子葉鞘顏色和葉枕顏色在不同環(huán)境下的GWAS 曼哈頓結(jié)果圖Fig.1 GWAS Manhattan Result Map of leaf sheath color and pulvinus color in Different Environment

使用Camoco 算法中SNP 到基因位點定位（SNP to gene mapping）的功能，對上述基因進行篩選，得到一個高置信度的候選基因列表（表1）。在葉鞘和葉枕中分別篩選到了8個和5個與顏色相關(guān)的候選基因。利用CandiHap［28］，繪制候選基因的結(jié)構(gòu)圖，并繪制出基因所在的染色體位置，非編 碼 區(qū)UTR（Untranslated Region）、編 碼 序 列CDS（Coding sequence）的位置和大小，包括相應(yīng)發(fā)生在外顯子區(qū)域的突變位置（增強出版附件2）。

表1 Camoco 篩選結(jié)果Table 1 Camoco screening results

2.2 基因共表達網(wǎng)絡(luò)模塊的構(gòu)建及篩選結(jié)果

對表達矩陣中的基因進行過濾，獲得31 980個高表達基因。通過WGCNA 算法，選擇權(quán)重值β=10 構(gòu)建無尺度網(wǎng)絡(luò)，按照混合動態(tài)剪切的標準劃分模塊，最終獲得37個模塊（圖2）。將Camoco分析中得到的候選基因與模塊中的基因進行比較，確定候選基因在模塊中的分布情況（圖3）。

圖2 基因聚類樹和樣品分割Fig.2 Gene cluster dendrograms and module detecting

2.3 相關(guān)模塊及候選基因分析結(jié)果

對谷子葉鞘顏色和葉枕顏色候選基因所在的模塊分別進行GO 功能富集分析，并使用R 語言中g(shù)gplot2 包，將其進行可視化，發(fā)現(xiàn)2 者共同富集到了與質(zhì)體（GO：0009536），葉綠體（GO：0044434），四吡咯生物合成過程（GO：0033014），色素的生物合成過程（GO：0046148），光反應(yīng)（GO：0019684）等生物學過程（圖4，增強出版附件2 圖3，附件3 表1）。以上的富集結(jié)果有力的證明了Camoco 挖掘的基因較為全面。

圖3 候選基因在模塊中的分布情況Fig.3 The distribution of candidate genes in the module

整合Camoco 分析得到的基因和GO 富集分析的結(jié)果（增強出版附件3），本研究共獲得了10個與顏色相關(guān)的候選基因，推測谷子中Si4g01380、Si4g01390、Si4g01300、Si7g20820、Si7g20840、Si7g20880對葉鞘顏色起調(diào)控作用，Si4g01380、Si4g01390、Si4g03530、Si4g03510對谷子葉枕顏色起調(diào)控作用。

2.4 谷子葉鞘顏色和葉枕顏色候選基因的表達模式分析

對谷子葉鞘顏色和葉枕顏色這2種性狀候選基因進行表達模式分析，發(fā)現(xiàn)這些基因在整個生育期都有表達，但普遍在苗期中表達量較高，說明在這些基因的調(diào)控下谷子表型在苗期開始出現(xiàn)差異（增強出版附件2，圖4～圖11）。如谷子葉鞘顏色的候選基因Si401300谷子全生育期都有表達但在苗期中表達量最高，在旗葉葉鞘和倒4 葉葉鞘中的表達量分別為90.225 和114.745（圖5）。

圖4 與谷子葉鞘顏色相關(guān)的3個模塊中共表達基因的GO 功能富集Fig.4 GO function enrichment of co-expressed genes in 3 modules related to the color of foxtail millet leaf sheath.

圖5 谷子Si4g01300 基因表達模式Fig.5 Expression pattern of Si4g01300 gene in foxtail millet

2.5 互作網(wǎng)絡(luò)的構(gòu)建及樞紐基因功能注釋

選取葉鞘顏色候選基因所在的模塊并將其作為核心基因使用Cytoscape 網(wǎng)絡(luò)圖工具，選擇軟閾值作為連通性計算方法，構(gòu)建局部共表達網(wǎng)絡(luò)圖（圖6）。同理構(gòu)建谷子葉枕顏色的共表達網(wǎng)絡(luò)圖（增強出版附件2 圖12）。同時，用WGCNA 深入挖掘相關(guān)模塊內(nèi)的樞紐基因，并對未知功能的基因進行功能注釋。本試驗對共表達網(wǎng)絡(luò)圖中的樞紐基因進行了功能注釋（增強出版附件4 表1），為以后驗證未知功能的基因提供參考。如在葉鞘顏色 的 green 模 塊 中Si1g38040、Si1g18850、Si2g07320等27個樞紐基因都注釋到與葉綠體相關(guān)；Si3g17790、Si3g17800直接注釋到細胞色素P450；Si3g31460、Si2g36540等基因注釋與水解酶相關(guān)，水解酶與植物的衰老密切相關(guān)。

圖6 谷子葉鞘顏色基因在相應(yīng)的3種模塊中的局部共表達網(wǎng)絡(luò)圖Fig.6 Local co-expression network of sheath color genes in the corresponding 3 modules

3 討論

GWAS 分析與傳統(tǒng)的連鎖分析相比較具有耗時少、廣度大、精確度高等優(yōu)勢，但是GWAS 研究中很容易出遺傳性缺失現(xiàn)象，以往為解決此問題通常通過增加樣本量，提高檢測效率，或以整個代謝通路或其功能類型作為單位進行關(guān)聯(lián)以及富集分析，通常不關(guān)注個體的功能基因。另外從SNP位點到基因的定位通常是人工檢查篩選，因此增加了很多人為因素，很可能錯過真正對性狀起調(diào)控作用的基因。因此有人利用一些功能基因組數(shù)據(jù)的互補性來解釋GWAS，在人類和擬南芥的研究中建立了這樣的共表達網(wǎng)絡(luò)［29-30］，此方法被稱為GWAB（genome-wide association boosting），但大部分作物功能基因的基因注釋數(shù)據(jù)都不完整，此方法普遍應(yīng)用較為困難。因此本文使用Camoco算法利用基因表達量篩選候選基因，既大大減少了候選基因的數(shù)量，又可以將表達模式相似的基因篩選出來，減少GWAS 分析定位到假陽性基因的概率。本研究通過GWAS 分析，在葉鞘顏色中篩選到10個相關(guān)的SNP 位點，葉枕顏色中篩選到11個相關(guān)的SNP 位點，此區(qū)段共預測到相關(guān)基因278個。Camoco 算法計算后共得到12個基因，進一步結(jié)合WGCNA 分析共得到10個候選基因，且這些基因都直接或間接與植物體內(nèi)葉綠體的功能有關(guān)［31-36］。如谷子葉鞘顏色的候選基因Si4g01300編碼了一種DUF3754 域蛋白質(zhì)，其在谷子中的具體功能還不是很明確，但有研究人員發(fā)現(xiàn)一些DUF 家族基因與植物發(fā)育和細胞的生長、葉綠體運動等功能有關(guān)［37］，已有研究表明葉鞘顏色會受多對隱性基因或者單顯性基因控制，并證明這些基因調(diào)控的代謝途徑會受光照等因子的影響［38］。另外谷子葉枕顏色的候選基因Si4g03530的推定功能是泛素連接酶家族，E3 泛素連接酶，在玉米的穗位葉的葉綠體中檢測到該基因的表達［39］，并且該連接酶家族可能與植物細胞的自然凋亡相關(guān)。雙形狀共定位到的基因Si4g01380則與半胱氨酸合成酶在葉肉細胞葉綠體中對無機硫的催化有關(guān)，在基因GO 注釋中主要表達于質(zhì)體組織等細胞組分內(nèi)，其編碼GRAS 家族轉(zhuǎn)錄因子，該家族中有多個轉(zhuǎn)錄因子參與光敏色素（phytochrome）信號傳導［40～42］，光敏色素為色素蛋白結(jié)構(gòu)，是葉綠素等四吡咯物質(zhì)的基本類型。除此之外，通過共表達網(wǎng)絡(luò)挖掘到的樞紐基因也很重要，其中一些基因如Si3g31460、Si2g36540沒有直接富集到與色素合成相關(guān)的GO term，但富集到與蛋白質(zhì)水解，水解酶活性，轉(zhuǎn)移酶活性和激素相關(guān)的GO term，這些與植物衰老有密切關(guān)系，植物衰老對谷子的顏色性狀有顯著影響。

4 結(jié)論

本研究基于谷子葉鞘顏色和葉枕顏色的全基因組關(guān)聯(lián)分析，初步確定相關(guān)的SNP 位點，利用Camoco 軟件初步確定的SNP 位點進行篩選，挖掘到高置信度的候選基因，對其進行功能注釋，并構(gòu)建共表達網(wǎng)絡(luò)，對相關(guān)模塊進行GO 富集分析，得到了具有生物學意義的功能富集結(jié)果。總之，本研究利用Camoco 算法結(jié)合共表達網(wǎng)絡(luò)確定了10個候選基因，其中葉鞘顏色6個，葉枕顏色4個，Si4g01380和Si4g01390是兩個性狀共同定位到的候選基因。這些發(fā)現(xiàn)為深入探索調(diào)控谷子相關(guān)顏色的分子機制提供一定的理論依據(jù)，對培育高產(chǎn)谷子品種、發(fā)掘谷子形態(tài)標記以及開發(fā)各種功能性谷子具有重要意義。