王炯亮 趙韓生 高志民

(國際竹藤中心 北京 100102)

棕櫚藤屬于棕櫚科，多年生常綠攀援類植物，世界上共有600多個種，分別屬于13個屬，其中27個商業(yè)品種在熱帶地區(qū)被廣泛種植。棕櫚藤每年可產(chǎn)生大約70億美元的收益，超過500萬人在經(jīng)濟上依賴于棕櫚藤。棕櫚藤是熱帶森林生態(tài)系統(tǒng)的組成部分，其纖鞭帶有爪鉤，利于在森林環(huán)境中攀爬和生長；然而，這也使得棕櫚藤種植、管理和采收等經(jīng)營環(huán)節(jié)變得困難，導(dǎo)致棕櫚藤產(chǎn)品成本較高。因此，纖鞭生長發(fā)育的遺傳調(diào)控機制備受關(guān)注。作為竹藤基因組圖譜(GABR)項目的一部分，雖然國際竹藤中心在2018年發(fā)布了單葉省藤和黃藤2個代表性棕櫚藤的全基因組序列，但棕櫚藤轉(zhuǎn)錄組水平上的基因功能分析仍是空白。

基因協(xié)同功能網(wǎng)絡(luò)在構(gòu)建整體通路模型、細化基因注釋，以及在體內(nèi)模擬重要調(diào)控機制等方面發(fā)揮著關(guān)鍵作用。作為其中主要代表的共表達網(wǎng)絡(luò)(CEN)，它是基于基因轉(zhuǎn)錄表達水平以表征相關(guān)功能基因的強大工具，可在多種條件下批量注釋、聚類和探索海量基因。CEN分析和相應(yīng)的數(shù)據(jù)庫在植物中得到廣泛應(yīng)用，ATTED-II是一個提供了擬南芥、芥菜、大豆等9個物種的CEN數(shù)據(jù)集的數(shù)據(jù)庫，CcNET提供了基于CEN的二倍體和多倍體棉花功能基因比較分析，AraNet v2構(gòu)建了提供28個植物物種間同源關(guān)系的CEN，RiceNet v2是一個升級后的水稻CEN分析數(shù)據(jù)庫，Bamboo-NET提供了毛竹的CEN和優(yōu)化后的基因功能注釋。但是，棕櫚藤中的CEN研究仍處于空白。

一般CEN的構(gòu)建遵循以下3個步驟：計算基因間的相似性值；選定閾值提取基因?qū)π纬删W(wǎng)絡(luò)；在CEN中尋找功能模塊。每一步都有不同的方法可供選擇，所以構(gòu)建一個CEN可以有許多不同的組合方法。在本研究中，通過使用單葉省藤和黃藤的基因組和轉(zhuǎn)錄組數(shù)據(jù)，選擇成熟方法構(gòu)建CEN，最大化了基因的覆蓋率和模塊的數(shù)量，優(yōu)化對基因的功能注釋，在CEN中鑒定了覆蓋纖鞭不同發(fā)育階段的功能模塊，并使用基因家族分類、同源基因注釋、順式作用元件分析和基因本體分析等工具對優(yōu)化的基因功能注釋進行了評估。此外，構(gòu)建了棕櫚藤的CEN數(shù)據(jù)庫(Rattan-NET)，提供了潛在的功能模塊和基因功能注釋，方便研究人員使用。

1 材料與方法

1.1 樣品收集與測序

不同發(fā)育階段單葉省藤和黃藤纖鞭組織采自中國林業(yè)科學(xué)研究院廣州熱帶林業(yè)研究所(原文附表S1)。按照說明書，使用TRIzol試劑分離出RNA，并使用NanoDrop 2000分光光度計測定濃度和純度。所使用RNase-free DNase I試劑在37 ℃下處理30 min，以去除所提取RNA中殘余的DNA，之后將RNA在反轉(zhuǎn)錄系統(tǒng)下反轉(zhuǎn)錄成cDNA，并建庫測序。使用BGISEQ-500平臺對富集得到的文庫進行長度為100 bp的雙末端測序。同時，還收集了此前測序獲得的4個黃藤纖鞭組織的RNA測序數(shù)據(jù)，共使用了處于不同發(fā)育階段的24個單葉省藤和36個黃藤纖鞭轉(zhuǎn)錄組數(shù)據(jù)集進行分析。

1.2 數(shù)據(jù)處理和基因表達分析

使用FastQC v0.11.6，在默認參數(shù)下進行轉(zhuǎn)錄組數(shù)據(jù)的質(zhì)量分析。使用Trimmomatic v0.36在LEADING:3、TRAILING:3、SLIDINGWINDOWS 4:15、MINLEN:50、TOPPHRED 64等參數(shù)下對接頭序列和低質(zhì)量序列進行修剪和過濾。經(jīng)質(zhì)量控制后的序列被輸入到HISAT2 v2.1.0中進行轉(zhuǎn)錄組序列的回帖，參數(shù)如下：-min-intronlen 20、-max-intronlen 4000、-rna-strandness RF。HISAT2輸出的數(shù)據(jù)都被保留下來用于后續(xù)分析(原文附表S2)。使用Cufflinks v2.2.1進行FPKM值的計算，使用Cuffdiff v2.2.1評估差異表達基因。此外，此研究使用3σ準則(threshold = average(5% value) + 3*SD)計算FPKM閾值，用于篩除FPKM過低的數(shù)據(jù)(原文附圖S1)。

1.3 CEN的構(gòu)建

本研究采用的策略是將皮爾森相關(guān)系數(shù)(PCC)和互相排序(MR)整合在一起用于CEN的構(gòu)建(原文附圖S2)：1)使用FPKM計算基因?qū)χg的PCC值，并過濾掉弱相關(guān)的基因?qū)?，只保留強相關(guān)的基因?qū)?原文附圖S3)；2)使用PCC計算基因?qū)χg的MR值，并提取高置信度的基因?qū)?，剔除不可靠的基因?qū)Γ?)進行Receiver operating characteristic分析，用于篩選最優(yōu)的MR參數(shù)值。最后，本研究選定了單向MR<3+雙向MR<30的基因?qū)ψ鳛镃EN。

1.4 功能模塊的鑒定和基因注釋優(yōu)化

使用基于集團滲透算法的CFinder v2.0.6測定每個節(jié)點的密度，進而在網(wǎng)絡(luò)中鑒定模塊。為使基因的覆蓋度和模塊數(shù)量最大化，分別在單葉省藤和黃藤中選定k=6和k=5作為CFinder的參數(shù)(原文附圖S4)。此外，使用基因集富集分析(GSEA)工具對模塊進行功能注釋，不可靠的模塊在Fisher檢驗和多重假設(shè)檢驗過程中被剔除。然后，使用模塊和CEN數(shù)據(jù)對2個棕櫚藤物種的基因功能注釋進行優(yōu)化，即：用模塊的功能注釋對屬于該模塊的基因進行注釋優(yōu)化；對每個基因的共表達基因集進行GSEA，F(xiàn)DR<0.05的詞條也做為是該基因優(yōu)化后注釋。

1.5 棕櫚藤和擬南芥之間的直系同源基因鑒定

使用BLAST相互最佳匹配(RBH)方法在棕櫚藤和擬南芥之中鑒定直系同源基因，排位前三的RBH被認為是最優(yōu)的直系同源對，所有BLAST結(jié)果的E-value值分布的峰值作為閾值，小于該閾值的基因?qū)Ρ徽J定為次優(yōu)的直系同源基因?qū)Α?/p>

1.6 順式作用元件的顯著性分析

參照前人研究，使用基于Z-score和P-value的篩選策略用于順式作用元件的顯著性測試，在掃描棕櫚藤基因3 kb啟動子區(qū)域時，選擇P-value<0.05的元件模體作為顯著富集的調(diào)控元件。

1.7 捕光復(fù)合物基因家族的鑒定與系統(tǒng)發(fā)育分析

從TAIR網(wǎng)站(https://www.arabidopsis.org)下載擬南芥的捕光復(fù)合物(LHC)基因序列，作為查詢序列，在單葉省藤、黃藤和水稻基因序列組成的數(shù)據(jù)庫中進行相似性搜索。搜索工具為BLAST+，參數(shù)為E-value<1e-5。相似性搜索的結(jié)果進一步使用pfamscan.pl腳本進行檢驗，提取該基因家族的編碼序列，用于后續(xù)分析。

在系統(tǒng)發(fā)育分析中，使用MUSCLE v3.8.31在默認參數(shù)下進行多序列比對。使用Gblocks server提取保守區(qū)序列。使用jModeltest v2.1.6尋找最優(yōu)的替換模型，參數(shù)如下：-f、-g 4、-i、-s 203、-S BEST。最后，使用PhyML v20120412進行系統(tǒng)發(fā)育樹的構(gòu)建，其中模型參數(shù)由jModeltest得到，自檢參數(shù)設(shè)為1 000。

1.8 查詢與可視化平臺構(gòu)建

Rattan-NET的構(gòu)建基于Linux系統(tǒng)、Apache軟件、MySQL軟件和PHP語言。涉及計算相關(guān)的分析工具是基于Python、Perl和R腳本語言開發(fā)的。Cytoscape.js是Linux版本的Cytoscape，它用于模塊與CEN的動態(tài)展示。GBrowse和SequenceServer用于提供序列掃描和BLAST+服務(wù)。

2 結(jié)果和分析

2.1 CEN的構(gòu)建和功能模塊的鑒定

使用處于不同發(fā)育階段的24個單葉省藤和36個黃藤的纖鞭樣本進行CEN的構(gòu)建。分別在單葉省藤和黃藤確定了符合單向MR<3+雙向MR≤30條件的630 081和670 502個基因?qū)ψ鳛镃EN的邊(原文表1)。單葉省藤和黃藤的CEN分別包含了31 847和36 769個基因，它們各自覆蓋了62.16%和68.93%的基因(原文表1)，2個棕櫚藤CEN平均邊的數(shù)量相近。CEN基因所相連邊的數(shù)量密度圖表明只有很少一部分基因有著高連接度，這符合生物網(wǎng)絡(luò)的一般特征(原文附圖S5A)。

使用基于集團滲透算法的CFinder軟件和GSEA工具(閾值為FDR<0.05)在2個CEN中鑒定潛在的功能模塊，分別在單葉省藤和黃藤中鑒定得到了3 504和3 027個潛在的功能模塊(原文表1)。它們的分布圖表明隨著每個模塊包含基因數(shù)量的增大，該規(guī)模的模塊數(shù)量下降(原文附圖S5B)。單葉省藤和黃藤功能模塊平均包含的基因數(shù)為8.1和6.1(原文表1)，這些功能模塊使用包含了PO、GFam、GO和KEGG等基因集的GSEA工具進行功能注釋。注釋結(jié)果表明，這些功能模塊覆蓋了光合作用、植物次級細胞壁生物合成、木質(zhì)素生物合成、類黃酮生物合成和苯丙素生物合成等在發(fā)育中重要的生物過程。進一步分析模塊中顯著富集的GO詞條(Z-score>4)(原文附表S3)，發(fā)現(xiàn)了2個棕櫚藤物種之間保守的詞條，都包括泛素依賴蛋白代謝過程、質(zhì)子跨膜轉(zhuǎn)運、囊泡介導(dǎo)轉(zhuǎn)運、細胞分裂素響應(yīng)等。此外，Rattan-NET網(wǎng)站的Download頁面提供了由GSEA產(chǎn)生的功能注釋文件。共有54.6%的單葉省藤和62.4%黃藤的基因在模塊和CEN分析中得到注釋優(yōu)化。

2.2 棕櫚藤中類黃酮和木質(zhì)素生物合成通路的CEN分析

使用RBH方法在2個棕櫚藤物種中鑒定參與苯丙素生物合成通路的PAL、C4H和4CL基因(原文圖1)，結(jié)果表明，單葉省藤中有8個PAL基因、12個C4H基因、9個4CL基因，黃藤中有13個PAL基因、12個C4H基因和15個4CL基因(原文附圖S4)。通過檢測功能注釋，在這些基因中調(diào)查與類黃酮生物合成和木質(zhì)素生物合成相關(guān)的共表達基因(原文圖1)。結(jié)果表明，在單葉省藤和黃藤的上述3類基因家族中分別有12個和14個類黃酮生物合成相關(guān)的基因(原文圖1B與附表S5)，與發(fā)現(xiàn)的類黃酮生物合成相關(guān)基因相比，在單葉省藤和黃藤中分別只有2個和5個與木質(zhì)素生物合成相關(guān)的PAL、C4H和4CL基因(原文圖1)。由此表明，在棕櫚藤的PAL、C4H和4CL這3個基因家族中，類黃酮相關(guān)的基因比木質(zhì)素生物合成相關(guān)的基因更多。

2.3 棕櫚藤中次級細胞壁生物合成的CEN分析

選擇參與木質(zhì)素生物合成和次級細胞壁(SCW)轉(zhuǎn)錄調(diào)控的代表性基因，構(gòu)建它們的CEN(原文圖2A、B)。這些基因包括4CL1(Calsi_gene34733和Daeje_Gene51484)、CCR1(Calsi_gene01533和Daeje_Gene15213)和MYB103(Calsi_gene15542和Daeje_Gene04728)。通過分析，在網(wǎng)絡(luò)中鑒定到參與SCW生物合成的調(diào)控因子，在單葉省藤中有MYB20、MYB54和MYB52，而在在黃藤中有MYB55和MYB4，它們可以被回貼到SCW轉(zhuǎn)錄調(diào)控網(wǎng)絡(luò)的簡化圖中(原文圖2A、B、C)。因為木質(zhì)素是SCW中最重要的成分之一，本研究也在網(wǎng)絡(luò)中鑒定了參與木質(zhì)素生物合成的相關(guān)基因，并將它們回貼到簡化圖中(原文圖2A、B、C)。此外，作為纖維素生物合成酶復(fù)合物的核心成分的纖維素合酶基因也在網(wǎng)絡(luò)中被鑒定出來，包括IRX3(Calsi_gene01341和Calsi_gene30420)、IRX1(Calsi_gene32935)、CESA4(Calsi_gene08952)、IRX6(Daeje_Gene45742)、IRX3(Daeje_Gene59357)、CESA4(Daeje_Gene59357)(原文圖2A、B)。隨后通過在CCR1和4CL1基因上游3 kb區(qū)域的順式作用分析鑒定得到顯著富集(P-value<0.05)的“MYB61”元件(原文附表S6)。

使用包含GO、KEGG和GFam等數(shù)據(jù)集的GSEA工具，對2個CEN進行分析以檢測潛在的生物過程(原文圖2D)。SCW生物合成(GO:0009834)、木質(zhì)素生物合成(GO:0009809)和類黃酮生物合成(map00941)等GO詞條同時在單葉省藤和黃藤中被富集，這表明這些生物學(xué)過程在2個棕櫚藤物種之間具有保守性。此外，在單葉省藤中富集到了纖維素生物合成過程(GO:0030244)、黃酮和黃烷醇生物生物合成(map00944)等詞條，而黃藤中富集到了木質(zhì)素分解代謝過程(GO:0046274)、苯丙素生物合成(map00940)和苯丙素生物合成基因家族等詞條。

另外，還鑒定出與SCW生物合成相關(guān)的潛在功能模塊(原文附圖S6)。在GSEA中，SCW合成相關(guān)的詞條得到顯著富集(FDR<0.05)，包括纖維素生物合成過程(GO:0030244)、植物細胞壁生物合成家族、纖維素合成酶類似物，木質(zhì)素生物合成過程(GO:0009809)、細胞版合成等。

2.4 棕櫚藤中光合作用相關(guān)基因的分析

在棕櫚藤中鑒定得到了LHC基因，并重構(gòu)了系統(tǒng)發(fā)育樹以研究棕櫚藤中LHC基因家族的進化(原文圖3A和附表S7)。系統(tǒng)發(fā)育樹表明，LHC基因家族可分成13個亞家族(原文圖3B)。在單葉省藤和黃藤中分別有21和23個LHC基因家族成員，這與擬南芥中LHC基因的數(shù)量(21)類似，比水稻中(15)的多。雖然棕櫚藤與擬南芥的LHCB1亞家族成員數(shù)量相近，但是其它亞家族成員數(shù)量的分布卻并不相似。比如，在棕櫚藤LHCB6、LHCA2、LHCB5和LHCB3基因亞家族中各有著雙倍于擬南芥的成員數(shù)量。

通過GSEA分析LHC基因家族的共表達基因(原文附圖S7)，發(fā)現(xiàn)光合作用相關(guān)的GO詞條，如光合作用(GO:0015979)、葉綠體類囊體膜(GO:0009535)和光刺激響應(yīng)(GO:0009416)等顯著富集(FDR<0.05)，這與前人研究LHC蛋白的主要功能是通過捕獲光能供給光化學(xué)反應(yīng)的結(jié)果相符，也表明了2個棕櫚藤CEN具有可靠性。同時，通過對棕櫚藤LHCA1基因的PCC top300 CEN和擬南芥LHCA1基因的PCC top300 CEN(下載自ATTED-II和AraNet)之間的比較分析發(fā)現(xiàn)，2個網(wǎng)絡(luò)存在高度的一致性，而且還發(fā)現(xiàn)一些LHC直系同源基因?qū)Γ@同樣表明網(wǎng)絡(luò)的可靠性。本研究試著通過分析黃藤中不同LHC基因之間的共表達關(guān)系去研究它們之間可能的合作關(guān)系(原文圖3C)，發(fā)現(xiàn)大多數(shù)LHC基因關(guān)系較近，但屬于LHCB1亞家族的4個基因(Daeje_Gene30151、Daeje_Gene53180、Daeje_Gene67755和Dae_Gene18901)游離在其它LHC基因之外。

2.5 Rattan-NET：一個在線的棕櫚藤CEN數(shù)據(jù)庫

為提供方便直觀的棕櫚藤CEN分析功能，構(gòu)建了Rattan-NET數(shù)據(jù)庫。該數(shù)據(jù)庫整合了CEN分析、順式作用元件分析、GSEA、GBrowse和Sequenceserver等工具，有助于研究人員在轉(zhuǎn)錄水平上對棕櫚藤基因功能注釋進行精煉。其中，Cytoscape網(wǎng)頁版工具提供了CEN和模塊的交互式展示功能。在Search頁面，用戶可以通過輸入單個或多個基因的ID進行CEN搜索。Rattan-NET還提供了GO、KEGG和GFam等功能注釋的基因ID和關(guān)鍵詞搜索功能，在基因詳情頁提供了每個基因在不同組織下的FPKM值，GBrowse提供了基因組瀏覽和序列提取功能?？傊?，Rattan-NET為用戶提供了全面的數(shù)據(jù)和便捷的工具。此外，它會隨著分析工具的擴充、棕櫚藤物種數(shù)據(jù)的增加將持續(xù)得到維護和優(yōu)化。

3 討論

國際竹藤中心的研究團隊于2018年發(fā)布了2個棕櫚藤物種(單葉省藤和黃藤)的基因組序列和它們在全基因組水平上的基因功能注釋，這為棕櫚藤分子研究提供了重要的基礎(chǔ)數(shù)據(jù)。本研究通過使用上文描述的策略構(gòu)建CEN，將轉(zhuǎn)錄組層面的信息引入到基因功能注釋。2個CEN符合無標度網(wǎng)絡(luò)的特征，表明了它們的可靠性(原文附圖S5A)。在單葉省藤和黃藤中分別有3 504個和3 027個，覆蓋了多個在纖鞭生長發(fā)育過程中重要的生物學(xué)過程的模塊得到了GSEA鑒定。泛素依賴蛋白代謝過程、質(zhì)子跨膜轉(zhuǎn)運、囊泡介導(dǎo)轉(zhuǎn)運、細胞分裂素響應(yīng)等詞條是2個棕櫚藤物種之間功能模塊保守的GO注釋。然而，只有62.16%的單葉省藤和68.93%的黃藤基因被CEN覆蓋，這低于在毛竹中>90%的CEN覆蓋率，這可能是棕櫚藤的樣本數(shù)量過少和單葉省藤與黃藤的基因組scaffold過短導(dǎo)致的。隨著對竹藤基因組的維護和GABR計劃的推進，將得到棕櫚藤中更高質(zhì)量的基因組、基因注釋和源自更多組織的轉(zhuǎn)錄組數(shù)據(jù)集，這將為未來CEN的構(gòu)建和分析提供全面且高質(zhì)量的數(shù)據(jù)。同時，其他棕櫚藤物種的CEN將會被構(gòu)建并進行棕櫚藤物種之間的CEN比較分析。

4 結(jié)論

本研究利用單葉省藤和黃藤基因組數(shù)據(jù)和不同發(fā)育階段纖鞭的轉(zhuǎn)錄組數(shù)據(jù)，構(gòu)建了CEN并優(yōu)化了基因功能注釋，并通過功能富集分析、基因家族分類、順式作用元件分析和GO分析等證實了基因功能注釋和預(yù)測的可靠性。通過數(shù)據(jù)挖掘系統(tǒng)，豐富了對單葉省藤和黃藤重要農(nóng)藝性狀的遺傳基礎(chǔ)的理解，并對其他重要性狀研究提供了新的思路。構(gòu)建的在線的CEN數(shù)據(jù)庫(Rattan-NET)，將有助于棕櫚藤CEN的應(yīng)用，對促進棕櫚藤分子生物學(xué)發(fā)展，揭示棕櫚藤重要性狀分子調(diào)控機制具有重要科學(xué)價值。

原文出處

Wang J, Ma X, Yang J, Hui Y, She J, Tian T, Li Z, Xu W, Gao Z, Su Z, Zhao H. Coexpression analysis reveals dynamic modules regulating the growth and development of cirri in the rattans (CalamussimplicifoliusandDaemonoropsjenkinsiana).Frontiers in Genetics, 2020,11:378. DOI: 10.3389/fgene.2020.00378. eCollection 2020.