橡膠樹(shù)種子性狀相關(guān)基因的預(yù)測(cè)與功能研究

劉嘉智，劉長(zhǎng)寧

（1.中國(guó)科學(xué)院西雙版納熱帶植物園，云南勐臘 666303；2.中國(guó)科學(xué)院大學(xué)，北京 100049）

橡膠籽是天然橡膠種植生產(chǎn)中的重要副產(chǎn)品，橡膠籽的含油量為35%～50%，年產(chǎn)量可達(dá)2 060 kg/hm2［1］，與其 他 非 食用 植 物 油相 比，橡膠籽油燃燒的CO2和NO2排放率更低，是理想的生物柴油來(lái)源［2-3］。然而，目前對(duì)橡膠樹(shù)的利用主要集中于獲取天然橡膠，忽略了橡膠籽作為生物柴油的巨大潛力，與之相關(guān)的研究也并不充足。

通常，具有同源序列的蛋白編碼基因具有相似的功能，因而基于同源分析可以在近緣物種中預(yù)測(cè)相關(guān)功能基因。本研究基于共線性分析和同源分析，利用目前研究中較為深入的蓖麻全基因組關(guān)聯(lián)分析數(shù)據(jù)，將蓖麻野生型品種中種子農(nóng)藝性狀相關(guān)的基因映射到橡膠樹(shù)reyan7-33-97 品種基因組上，使用生物信息學(xué)手段預(yù)測(cè)橡膠樹(shù)中種子性狀相關(guān)基因并進(jìn)行基因功能分析，為橡膠籽的合理開(kāi)發(fā)利用提供理論參考。

1 材料和方法

1.1 數(shù)據(jù)來(lái)源

蓖麻野生型品種的基因組文件以及基因組注釋文件下載于http://oilplants.iflora.cn/Download/castor_download.html，橡膠樹(shù)reyan7-33-97 品種的基因組文件以及基因組注釋文件從https://www.ncbi.nlm.nih.gov/data-hub/genome/GCF_001654055.1/下載。從上述文件分別提取基因ID、序列信息、基因組位置信息和功能描述等內(nèi)容。從NCBI SRA 數(shù)據(jù)庫(kù)中收集具有文獻(xiàn)來(lái)源且為雙端測(cè)序的蓖麻和橡膠樹(shù)根、莖、種子3 種組織的轉(zhuǎn)錄組數(shù)據(jù)，樣本編號(hào)分別為SRR3136156、SRR6127585、SRR3136168、SRR17218208、SRR17218209、SRR17218216。

1.2 研究方法

1.2.1 橡膠樹(shù)種子性狀相關(guān)基因預(yù)測(cè)

根據(jù)文獻(xiàn)調(diào)研［4］，在蓖麻野生型品種中收集基于全基因組關(guān)聯(lián)分析得到的種子特征性狀相關(guān)的SNP 位點(diǎn)。根據(jù)基因組注釋信息，在上述SNP 位點(diǎn)前后50 kb 范圍內(nèi)提取基因，并認(rèn)為這些基因可能與種子特征有關(guān)。

橡膠樹(shù)與蓖麻同屬大戟科，具有較近的親緣關(guān)系，因此，首先判斷蓖麻種子性狀相關(guān)基因所在區(qū)塊在橡膠樹(shù)中是否存在共線區(qū)塊。由于橡膠樹(shù)reyan7-33-97 品種基因組只組裝到scaffold 級(jí)別，使用BLAST 進(jìn)行共線性分析后認(rèn)為存在3 個(gè)及以上基因，且基因之間相互間隔不超過(guò)50 kb的片段具有共線性。然后，使用Orthofinder 軟件［5］對(duì)這兩個(gè)物種進(jìn)行同源組鑒定，得到的結(jié)果再根據(jù)雙向BLAST 最優(yōu)策略進(jìn)一步確定同源基因，并根據(jù)基因結(jié)構(gòu)驗(yàn)證預(yù)測(cè)的同源基因是否可靠，最終將收集到的蓖麻野生型種子性狀相關(guān)基因比對(duì)到橡膠樹(shù)基因組上，來(lái)挖掘橡膠樹(shù)中與種子特征性狀相關(guān)的基因。

1.2.2 表達(dá)譜分析和差異表達(dá)分析

下載蓖麻和橡膠樹(shù)根、莖、種子3 種組織對(duì)應(yīng)的轉(zhuǎn)錄組原始測(cè)序文件，使用Trim Galore 進(jìn)行質(zhì)控檢測(cè)去除低質(zhì)量堿基和接頭。利用HISAT2 軟件［6］對(duì)基因組序列構(gòu)建索引后進(jìn)行序列比對(duì)，使用Stringtie 工具［7］進(jìn)行定量分析來(lái)獲取基因水平上的Raw counts、TPM 值以及FPKM 值。將TPM 值進(jìn)行歸一化處理，使用R 包pheatmap 繪制表達(dá)量熱圖，并對(duì)上述兩個(gè)物種得到的表達(dá)量結(jié)果進(jìn)行比較。使用GFOLD 1.1.4 軟件［8］對(duì)橡膠樹(shù)基因進(jìn)行差異表達(dá)分析，進(jìn)一步挖掘橡膠樹(shù)中與種子性狀相關(guān)基因在種子與根、莖之間的表達(dá)差異。

1.2.3 miRNA 靶位點(diǎn)預(yù)測(cè)

根據(jù)文獻(xiàn)調(diào)研［9］以及pmiren 數(shù)據(jù)庫(kù)［10］中獲取橡膠樹(shù)的成熟miRNA 序列，使用gffread 軟件提取轉(zhuǎn)錄本序列，使用psRNATarget 在線工具［11］默認(rèn)參數(shù)進(jìn)行miRNA 靶基因預(yù)測(cè)，從預(yù)測(cè)結(jié)果中篩選期望值大于3.5 的相互作用的miRNA 和靶基因，使用Cytoscape 軟 件［12］將miRNA 與 靶 基因 間的 網(wǎng)絡(luò)關(guān)系進(jìn)行可視化展示。

1.2.4 蛋白互作網(wǎng)絡(luò)預(yù)測(cè)

基于同源組劃分，以擬南芥蛋白互作網(wǎng)絡(luò)作為參考映射到橡膠樹(shù)中，獲得橡膠樹(shù)reyan7-33-97 的蛋白互作網(wǎng)絡(luò)。從String 數(shù)據(jù)庫(kù)［13］中下載擬南芥蛋白互作物理子網(wǎng)，使用Orthofinder 軟件預(yù)測(cè)擬南芥與橡膠樹(shù)reyan7-33-97 之間的蛋白同源群組。使用Cytoscape 軟件對(duì)橡膠樹(shù)種子性狀相關(guān)候選基因進(jìn)行網(wǎng)絡(luò)可視化。

2 結(jié)果與分析

2.1 橡膠樹(shù)種子性狀相關(guān)基因預(yù)測(cè)

在野生型蓖麻中共發(fā)現(xiàn)18 個(gè)與種子特征性狀相關(guān)的SNP 位點(diǎn)，這些SNP 位點(diǎn)分別與種子面積（SA）、種子長(zhǎng)度（SL）、種子厚度（ST）、種子寬度（SW）以及單粒重（SSW）等5 個(gè)特征形狀相關(guān)。根據(jù)蓖麻基因組注釋文件，在上述SNP 位點(diǎn)上下50 kb 范圍內(nèi)共有67 個(gè)基因（圖1），分別位于蓖麻的1、2、3、5、6 以及9 號(hào)染色體上（表1）。其中Rc01G001599、Rc01G001600、Rc01G001601、Rc01 G001602、Rc01G001603、Rc01G001604、Rc01G001605、Rc01G001606、Rc01G001607、Rc01G001608、Rc01G 001609 對(duì)上述5 個(gè)性狀都產(chǎn)生影響。

表1 蓖麻野生型種子特征相關(guān)基因

圖1 蓖麻野生型種子性狀相關(guān)基因統(tǒng)計(jì)

基于共線性分析，在蓖麻野生型與橡膠樹(shù)reyan7-33-97 中發(fā)現(xiàn)4 個(gè)共線區(qū)塊，分別位于蓖麻的2、5、6 以及9 號(hào)染色體上，對(duì)應(yīng)橡膠樹(shù)的NW_018746954.1、NW_018745856.1、NW_018746236.1 以及NW_018746078.1 這4 條scaffold，共 覆 蓋19 個(gè)基因（表2）。根據(jù)Orthofinder 軟件，在蓖麻野生型與橡膠樹(shù)reyan7-33-97 中共鑒定到16 926 個(gè)同源組，根據(jù)雙向BLAST 最優(yōu)確定二者之間序列相似程度，上述67 個(gè)蓖麻基因最終在橡膠樹(shù)中比對(duì)得到49 個(gè)基因。對(duì)上述橡膠樹(shù)基因以及同源蓖麻基因根據(jù)基因結(jié)構(gòu)作圖進(jìn)行驗(yàn)證（圖2），發(fā)現(xiàn)基因結(jié)構(gòu)與最終得到的基因序列預(yù)期結(jié)果相符。

表2 蓖麻野生型與橡膠樹(shù)reyan7-33-97共線區(qū)塊基因統(tǒng)計(jì)

圖2 蓖麻野生型與橡膠樹(shù)reyan7-33-97同源基因結(jié)構(gòu)圖

2.2 種子性狀相關(guān)表達(dá)譜和差異表達(dá)分析

為進(jìn)一步探究蓖麻和橡膠樹(shù)種子性狀相關(guān)基因的功能，進(jìn)行了組織表達(dá)譜的分析，圖3 是蓖麻和橡膠樹(shù)種子性狀相關(guān)基因的組織表達(dá)譜。從中發(fā)現(xiàn)，有些同源基因在3 個(gè)組織中的表達(dá)量類似，但有些差別很大。對(duì)蓖麻和橡膠對(duì)應(yīng)的49 個(gè)基因在種子中表達(dá)結(jié)果進(jìn)行統(tǒng)計(jì)，發(fā)現(xiàn)在兩個(gè)物種中有5 個(gè)基因均為正表達(dá)，25 個(gè)基因均為負(fù)表達(dá)，19 個(gè)基因表達(dá)相反，這可能是由于物種在分化過(guò)程中部分基因功能改變導(dǎo)致的。

圖3 蓖麻和橡膠樹(shù)種子性狀相關(guān)基因組織表達(dá)譜

根據(jù)差異表達(dá)分析，選出16 個(gè)在種子中表達(dá)量遠(yuǎn)高于根莖的基因、17 個(gè)在種子中表達(dá)量遠(yuǎn)低于根莖的基因以及5 個(gè)在三者之中都沒(méi)有差異表達(dá)的基因，其余基因在3 個(gè)組織中無(wú)交集（表3）。

表3 種子與根莖差異表達(dá)基因

2.3 miRNA 靶位點(diǎn)預(yù)測(cè)

基因的miRNA 靶位點(diǎn)預(yù)測(cè)也是研究基因調(diào)控的重要途徑。本研究共收集到橡膠樹(shù)的217 條成熟miRNA 序列，對(duì)橡膠樹(shù)種子性狀相關(guān)基因進(jìn)行miRNA 靶位點(diǎn)預(yù)測(cè)，結(jié)果發(fā)現(xiàn)LOC110671483、LOC110673189、 LOC110645416、 LOC110651200、LOC110672934、 LOC110641727、 LOC110653311、LOC110659520 受到不同miRNA 調(diào)控（圖4）。由于miRNA 在物種進(jìn)化中相當(dāng)保守，大部分miRNA 在特定的組織和發(fā)育階段表達(dá)，我們預(yù)測(cè)得到的miRNA 許多功能都與種子性狀相關(guān)。例如之前研究發(fā)現(xiàn)miR156 參與調(diào)控?cái)M南芥種子成熟過(guò)程［14］、miR408 在水稻種子發(fā)育中具有重要作用［15］等，因此推測(cè)上述8 個(gè)基因可能受到miRNA 調(diào)控影響種子發(fā)育，從而影響了種子特征性狀。

圖4 橡膠樹(shù)種子性狀相關(guān)基因miRNA靶位點(diǎn)預(yù)測(cè)

2.4 蛋白互作網(wǎng)絡(luò)預(yù)測(cè)

基于同源組的劃分，橡膠樹(shù)reyan7-33-97 品種中共有22 815 個(gè)基因映射到19 762 個(gè)擬南芥基因。根據(jù)預(yù)測(cè)得到的橡膠樹(shù)蛋白互作網(wǎng)絡(luò)，我們將橡膠樹(shù)reyan7-33-97 品種中與種子性狀相關(guān)的49 個(gè)基因放入蛋白互作網(wǎng)絡(luò)中，發(fā)現(xiàn)除LOC110654127 和LOC110668244 外，其 它 基 因 之 間不存在邊。這可能是影響種子性狀的橡膠樹(shù)候選基因與其他基因也存在相互作用導(dǎo)致的。為此，提取這49 個(gè)基因以及最近的相鄰基因重新進(jìn)行網(wǎng)絡(luò)構(gòu)建，結(jié)果發(fā)現(xiàn)其中16 個(gè)橡膠樹(shù)種子性狀候選基因通過(guò)一級(jí)相互連接構(gòu)成網(wǎng)絡(luò)（圖5），因此認(rèn)為這16 個(gè)基因存在間接蛋白相互作用從而影 響 種 子 性 狀，LOC110654127 和LOC110668244 同時(shí)存在直接蛋白相互作用影響種子性狀，而其它33 個(gè)基因可能單獨(dú)影響種子性狀或者存在不同于擬南芥的蛋白互作成員。

圖5 橡膠樹(shù)種子性狀候選基因蛋白互作網(wǎng)絡(luò)

3 展望

本研究使用全基因組關(guān)聯(lián)分析、基因組共線性分析和基因序列同源分析等多種技術(shù)手段，預(yù)測(cè)了49 個(gè)與橡膠樹(shù)種子性狀相關(guān)的基因并進(jìn)行了多種組學(xué)分析，研究結(jié)果為進(jìn)一步開(kāi)發(fā)利用橡膠籽提供了一定的工作基礎(chǔ)。然而，橡膠籽種子性狀基因與種子產(chǎn)油之間的關(guān)系是一個(gè)復(fù)雜的過(guò)程，植物激素的生物合成和信號(hào)傳遞對(duì)種子大小和形態(tài)的控制至關(guān)重要，例如赤霉素可以促進(jìn)種子的擴(kuò)大和發(fā)育［16］，而乙烯則可以抑制種子的生長(zhǎng)［17］。因此，后續(xù)研究可深入挖掘橡膠籽種子性狀基因和植物激素生物合成和信號(hào)傳遞通路之間的關(guān)系，通過(guò)對(duì)橡膠樹(shù)中這些激素合成和信號(hào)傳遞通路進(jìn)行深入分析，揭示橡膠籽種子性狀基因在種子性狀形成過(guò)程中的作用機(jī)制，以便更好地理解種子大小和形態(tài)的調(diào)控機(jī)制，為橡膠樹(shù)的育種和生產(chǎn)提供指導(dǎo)。