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

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

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

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

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

1 材料與方法

1.1 試驗(yàn)材料

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

轉(zhuǎn)錄組數(shù)據(jù)：名優(yōu)品種晉谷21（JG21）種植于山西農(nóng)業(yè)大學(xué)試驗(yàn)田，在其生長周期中，對18個(gè)組織時(shí)期取樣，每個(gè)樣設(shè)置3個(gè)重復(fù)，進(jìn)行轉(zhuǎn)錄組學(xué)分析。具體組織時(shí)期有：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 試驗(yàn)方法

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

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

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

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

1.2.3 Camoco 分析及流程

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

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

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

1.2.5 模塊的GO 富集分析

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

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

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

2 結(jié)果與分析

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 討論

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