芋全長(zhǎng)轉(zhuǎn)錄組測(cè)序分析及淀粉生物合成相關(guān)基因挖掘

何芳練，劉莉莉，蔣慧萍，韋紹龍，邱祖楊，董偉清

(1.廣西農(nóng)業(yè)科學(xué)院生物技術(shù)研究所，南寧 530007；2.荔浦市農(nóng)業(yè)農(nóng)村局，廣西 荔浦 546600；3.廣西亞熱帶作物研究所，南寧 530001)

【研究意義】芋[Colocasiaesculenta(L.) Schott]俗稱(chēng)芋頭，為天南星科單子葉植物，其利用歷史可追溯到28000年前的所羅門(mén)島[1]。芋的主要食用部位為富含淀粉和蛋白及其他營(yíng)養(yǎng)物質(zhì)的地下球莖，是亞太、非洲和美洲地區(qū)5億多人的主食和蔬菜[2-3]。我國(guó)有2000多年的芋栽培歷史，是世界上第二大芋生產(chǎn)國(guó)和最大出口國(guó)[2]。目前，雖然芋的基因組已被成功測(cè)序[4-5]，但其在轉(zhuǎn)錄水平的調(diào)控是一個(gè)復(fù)雜過(guò)程，包括可變剪接(AS)和可變多聚腺苷酸化(APA)等。AS和APA在植物生長(zhǎng)發(fā)育和脅迫應(yīng)答方面發(fā)揮重要作用[6-7]，但二代測(cè)序讀長(zhǎng)較短，無(wú)法準(zhǔn)確預(yù)測(cè)完整的全長(zhǎng)序列，對(duì)轉(zhuǎn)錄本的結(jié)構(gòu)分析存在較大困難，而基于單分子長(zhǎng)讀數(shù)測(cè)序技術(shù)(SMRT)的三代全長(zhǎng)轉(zhuǎn)錄組測(cè)序具有讀長(zhǎng)超長(zhǎng)優(yōu)勢(shì)，可直接獲得完整的全長(zhǎng)轉(zhuǎn)錄本，能準(zhǔn)確識(shí)別轉(zhuǎn)錄本的同源異構(gòu)體，在分析轉(zhuǎn)錄本AS、APA、融合基因和等位基因等方面具有非常大的優(yōu)勢(shì)[8]。因此，開(kāi)展芋三代全長(zhǎng)轉(zhuǎn)錄組測(cè)序，獲取芋全長(zhǎng)轉(zhuǎn)錄本信息并進(jìn)行功能注釋?zhuān)治鲛D(zhuǎn)錄本的AS、APA、長(zhǎng)鏈非編碼RNA(lncRNA)、轉(zhuǎn)錄因子(TF)及簡(jiǎn)單重復(fù)序列(SSR)等結(jié)構(gòu)信息，同時(shí)挖掘淀粉生物合成相關(guān)的基因及轉(zhuǎn)錄本信息，對(duì)豐富芋基因序列及結(jié)構(gòu)信息，為后續(xù)闡明淀粉生物合成的分子機(jī)制及深入挖掘芋基因資源具有重要意義。【前人研究進(jìn)展】測(cè)序技術(shù)已廣泛應(yīng)用于開(kāi)展芋種質(zhì)資源、基因組、球莖發(fā)育機(jī)理、顏色形成及分子標(biāo)記等方面的研究。在基因組研究方面，Bellinger等[4]完成了夏威夷地方芋品種Moi的基因組測(cè)序，組裝的基因組大小為2450 Mb，隨后Yin等[5]組裝高質(zhì)量染色體級(jí)別的芋基因組(龍香芋)，基因組大小為2405 Mb，高質(zhì)量基因組的公布為深入研究芋遺傳進(jìn)化及重要農(nóng)藝性狀形成的分子機(jī)制打下了基礎(chǔ)。在球莖發(fā)育和顏色形成方面，Dong等[3]對(duì)芋球莖發(fā)育過(guò)程進(jìn)行全轉(zhuǎn)錄組測(cè)序，富集到與淀粉和蔗糖代謝途徑相關(guān)的mRNA、CircRNA和miRNA分別為139、99和46個(gè)；He等[9]對(duì)芋球莖肉質(zhì)纖維顏色形成過(guò)程進(jìn)行轉(zhuǎn)錄組測(cè)序，鑒定出41和12個(gè)分別與類(lèi)黃酮和花青素相關(guān)的差異轉(zhuǎn)錄本。在種質(zhì)資源研究和分子標(biāo)記開(kāi)發(fā)方面，You等[10]從轉(zhuǎn)錄組測(cè)序中鑒定5278個(gè)SSR位點(diǎn)，并將68份芋種質(zhì)資源分為三大類(lèi)群；Wang等[11]從轉(zhuǎn)錄組測(cè)序中鑒定11363個(gè)SSR位點(diǎn)，使用18份芋種質(zhì)資源對(duì)隨機(jī)選取的150對(duì)引物進(jìn)行驗(yàn)證，結(jié)果顯示100對(duì)引物存在多態(tài)性信息含量值為0.042～0.778；Dong等[12]使用限制性位點(diǎn)關(guān)聯(lián)DNA測(cè)序(RAD-seq)鑒定4438個(gè)SSR位點(diǎn)，并將30份芋種質(zhì)資源分為三個(gè)類(lèi)群。Wang等[2]對(duì)234份芋種質(zhì)資源開(kāi)展特定長(zhǎng)度擴(kuò)增片段測(cè)序(SLAF-seq)，共獲得132 869個(gè)單核苷酸多態(tài)性(SNP)，并基于SNP標(biāo)記篩選出一套包含41份種質(zhì)資源的核心種質(zhì)資源庫(kù)。現(xiàn)階段，三代全長(zhǎng)轉(zhuǎn)錄組測(cè)序技術(shù)已在許多植物上成功應(yīng)用，如Wang等[8]對(duì)毛竹、Chen等[13]對(duì)紅花、潘敏等[14]對(duì)菠蘿蜜、尚驍堯等[15]對(duì)蒺藜苜蓿開(kāi)展全長(zhǎng)轉(zhuǎn)錄組測(cè)序，鑒定了大量AS事件、APA位點(diǎn)等轉(zhuǎn)錄本結(jié)構(gòu)信息，為深入研究基因的功能和轉(zhuǎn)錄調(diào)控機(jī)制提供了依據(jù)。【本研究切入點(diǎn)】雖然芋的基因組已被成功測(cè)序，但目前鮮見(jiàn)關(guān)于利用三代全長(zhǎng)轉(zhuǎn)錄組測(cè)序?qū)τ筠D(zhuǎn)錄本結(jié)構(gòu)進(jìn)行分析的報(bào)道。【擬解決的關(guān)鍵問(wèn)題】利用SMRT技術(shù)對(duì)芋不同組織(葉片、葉柄、球莖、匍匐莖和根)的混合樣開(kāi)展三代全長(zhǎng)轉(zhuǎn)錄組測(cè)序，鑒定新基因和轉(zhuǎn)錄本，分析AS、APA、lncRNA、TF及SSR位點(diǎn)等轉(zhuǎn)錄本的結(jié)構(gòu)信息，并挖掘淀粉生物合成相關(guān)基因和轉(zhuǎn)錄本，為后續(xù)闡明淀粉生物合成的分子機(jī)制及深入利用芋基因資源提供參考依據(jù)。

1 材料與方法

1.1 試驗(yàn)材料

供試材料為本課題組選育的芋新品種荔浦芋1號(hào)，2020年3月種植于廣西農(nóng)業(yè)科學(xué)院武鳴里建科學(xué)研究基地，農(nóng)事操作參考本課題組總結(jié)的檳榔芋水田輕簡(jiǎn)化高效栽培技術(shù)[16]。對(duì)播種90 d的植株(3月齡)進(jìn)行取樣，取樣部位為葉片、葉柄、球莖、匍匐莖和根，液氮速凍后保存于-80 ℃冰箱備用。

1.2 試驗(yàn)方法

1.2.1 總RNA提取與文庫(kù)構(gòu)建 參照天根生化科技(北京)有限公司的植物多糖多酚總RNA提取試劑盒(DP441)說(shuō)明提取所有樣品的總RNA。使用Nanodrop 2000(Thermo Fisher)檢測(cè)總RNA濃度和純度(OD260/280)，使用Agilent 2100(Agilent Technologies)檢測(cè)RNA的完整度(RIN值和28S/18S)。使用SMARTer PCR cDNA Synthesis Kit(TaKaRa，634926)將總RNA反轉(zhuǎn)錄成cDNA，然后將cDNA進(jìn)行PCR擴(kuò)增富集，使用AMpure PB(Pacbio，100-265-900)對(duì)擴(kuò)增產(chǎn)物進(jìn)行純化回收。將不同組織(葉片、葉柄、球莖、匍匐莖和根)的純化產(chǎn)物等量混合，然后使用SMRTbell?Express Template Prep Kit 2.0(PacBio，100-938-900)對(duì)混合產(chǎn)物進(jìn)行損傷修復(fù)、末端修復(fù)及連接接頭，構(gòu)建測(cè)序文庫(kù)。測(cè)序文庫(kù)質(zhì)檢合格后置于Sequel II(PacBio)測(cè)序儀上進(jìn)行三代全長(zhǎng)轉(zhuǎn)錄組測(cè)序。

1.2.2 全長(zhǎng)轉(zhuǎn)錄本序列獲取 全長(zhǎng)轉(zhuǎn)錄組序列獲取的過(guò)程主要包括全長(zhǎng)序列識(shí)別、全長(zhǎng)序列聚類(lèi)獲得一致序列和一致序列校正3個(gè)階段[17]。使用SMRT Link v7.0.0對(duì)下機(jī)數(shù)據(jù)進(jìn)行過(guò)濾、去除接頭獲得Subread序列，以Full passes≥3且序列準(zhǔn)確性>0.9的標(biāo)準(zhǔn)從Subread序列獲得環(huán)狀一致序列(CCS)，根據(jù)CCS中是否存在3′引物、5′引物和poly(A)獲得全長(zhǎng)非嵌合序列(FLNC)；隨后將與FLNC相似的序列聚成一簇(Cluster)，每簇得到一條一致序列；最后，對(duì)得到的一致序列進(jìn)行校正，獲得高質(zhì)量序列。將得到的高質(zhì)量序列通過(guò)GMAP v2017-11-15與芋參考基因組(Niue 2，https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA3 28799)進(jìn)行序列比對(duì)[18]，使用cDNA_Cupcake v6.1對(duì)比對(duì)結(jié)果去冗余，最終獲得非冗余全長(zhǎng)轉(zhuǎn)錄本序列。

1.2.3 新轉(zhuǎn)錄本的功能注釋和編碼區(qū)預(yù)測(cè) 將得到的新轉(zhuǎn)錄本序列與NR、eggNOG、SwissProt，GO、COG、KOG、Pfam和京都基因和基因組百科全書(shū)(KEGG)數(shù)據(jù)庫(kù)進(jìn)行BLAST比對(duì)，獲得新轉(zhuǎn)錄本的注釋信息。使用TransDecoder v5.0.0[19]預(yù)測(cè)新轉(zhuǎn)錄本的編碼區(qū)序列。

1.2.4 全長(zhǎng)轉(zhuǎn)錄本序列結(jié)構(gòu)分析 使用Astalavista v3.2[20]進(jìn)行AS分析；使用TAPIS pipeline v1.2.1[21]識(shí)別APA位點(diǎn)，通過(guò)MEME v4.9.1[22]對(duì)轉(zhuǎn)錄本poly(A)位點(diǎn)上游50 bp進(jìn)行motif分析；使用MISA v1.0進(jìn)行SSR分析；使用iTAK v1.7a[23]進(jìn)行TF預(yù)測(cè)；使用CPC2 v0.1[24]、CNCI v2.0[25]、PfamScan v1.6和CPAT v1.2.2[26]4種軟件對(duì)lncRNA進(jìn)行預(yù)測(cè)，根據(jù)lncRNA在參考基因組注釋信息(gff)上的位置，對(duì)lncRNA進(jìn)行分類(lèi)，并基于位置關(guān)系和互補(bǔ)序列2種方式對(duì)lncRNA的靶基因進(jìn)行預(yù)測(cè)；使用去冗余前的轉(zhuǎn)錄本進(jìn)行基因組跨區(qū)域預(yù)測(cè)，鑒定融合轉(zhuǎn)錄本。

1.2.5 淀粉生物合成相關(guān)基因挖掘 根據(jù)轉(zhuǎn)錄本在KEGG中淀粉與蔗糖代謝途徑(ko00500)的注釋信息，挖掘與芋淀粉生物合成相關(guān)的基因和轉(zhuǎn)錄本，并對(duì)其進(jìn)行AS和APA分析。

2 結(jié)果與分析

2.1 芋全長(zhǎng)轉(zhuǎn)錄組測(cè)序數(shù)據(jù)分析

為了鑒定盡可能多的轉(zhuǎn)錄本，從芋3月齡植株的葉片、葉柄、球莖、匍匐莖和根等組織部位提取高質(zhì)量的總RNA構(gòu)建混合樣測(cè)序文庫(kù)，在PacBio Sequel II平臺(tái)上進(jìn)行全長(zhǎng)轉(zhuǎn)錄組測(cè)序。下機(jī)數(shù)據(jù)經(jīng)過(guò)濾后獲得27.64 Gb測(cè)序數(shù)據(jù)，根據(jù)Full passes≥3且序列準(zhǔn)確性>0.9的條件從原始數(shù)據(jù)中提取CCS，共獲得CCS 275 469條，CCS總長(zhǎng)度為535 289 684 bp，平均為1943 bp，平均Full passes數(shù)為44。通過(guò)篩選含有5′引物、3′引物和poly(A)的CCS，共獲得FLNC 209 160條，占CCS總數(shù)的75.93%。將FLNC進(jìn)行聚類(lèi)得到一致序列85 053條，對(duì)一致序列進(jìn)行校正得到高質(zhì)量一致序列84 028條，將高質(zhì)量一致序列去冗余并與芋參考基因組進(jìn)行比對(duì)，最終得到38 043條全長(zhǎng)轉(zhuǎn)錄本序列。

2.2 新轉(zhuǎn)錄本的功能注釋及編碼區(qū)預(yù)測(cè)

通過(guò)與參考基因組進(jìn)行比對(duì)，共鑒定出新基因1878個(gè)，新發(fā)現(xiàn)轉(zhuǎn)錄本31 058條。將新發(fā)現(xiàn)的轉(zhuǎn)錄本與NR、eggNOG、Swiss-Prot、Pfam、KOG、KEGG、COG和GO等數(shù)據(jù)庫(kù)進(jìn)行比對(duì)，獲得各數(shù)據(jù)庫(kù)注釋的轉(zhuǎn)錄本數(shù)量在10 109～28 512條，共有28 785條轉(zhuǎn)錄本獲得功能注釋?zhuān)急?3.68%(表1)。在獲得注釋的新轉(zhuǎn)錄本中，300 bp<長(zhǎng)度<1000 bp的轉(zhuǎn)錄本數(shù)為1500條，長(zhǎng)度>1000 bp的轉(zhuǎn)錄本數(shù)為27 278條。

表1 芋新轉(zhuǎn)錄本的功能注釋統(tǒng)計(jì)Table 1 Functional annotation statistics of taro new transcripts

將新轉(zhuǎn)錄本序列比對(duì)到NR數(shù)據(jù)庫(kù)，獲取相似性最高的同源序列，統(tǒng)計(jì)比對(duì)到不同物種的序列數(shù)量和比例，結(jié)果顯示，排在前三位的同源物種為在KOG數(shù)據(jù)庫(kù)中，共有18 531條轉(zhuǎn)錄本被注釋?zhuān)⒏鶕?jù)功能將其歸為25類(lèi)，其中被注釋較多的功能分類(lèi)有：一般功能預(yù)測(cè)(General function prediction only，3700條)及翻譯后修飾、蛋白質(zhì)翻轉(zhuǎn)和分子伴侶(Posttranslational modification，protein turnover，chaperones，2150條)；信號(hào)轉(zhuǎn)導(dǎo)機(jī)制(Signal transduction mechanisms，2149條)和碳水化合物的運(yùn)輸和代謝(Carbohydrate transport and metabolism，1225條)。

Acanthamoebacastellanii、Capsasporaowczarzaki和Nematostellavectensis，比對(duì)到的同源序列分別有1920條(占比9.53%)、1303條(占比6.47%)和868條(占比4.31%)，此外，還比對(duì)到其他物種的序列11 940條(占比59.24%)。

在GO數(shù)據(jù)庫(kù)中，共有10 109條轉(zhuǎn)錄本被注釋?zhuān)譃榧?xì)胞組分(Cellular component)、分子功能(Molecular function)和生物過(guò)程(Biological process)三大類(lèi)，每大類(lèi)又細(xì)分為18、13和22個(gè)二級(jí)分類(lèi)。在細(xì)胞組分中，排在前三位的分別是細(xì)胞部分(Cell part，3830條)、細(xì)胞(Cell，3731條)和細(xì)胞器(Organelle，2447條)；在分子功能分類(lèi)中，催化活性(Catalytic activity)最多(6156條)，其次為連接(Binding，3941條)和轉(zhuǎn)運(yùn)器活性(Transporter activity，406條)；在生物過(guò)程中，細(xì)胞過(guò)程(Cellular process)最多(5383條)，其次為代謝進(jìn)程(Metabolic process，4989條)和單一有機(jī)體過(guò)程(Single-organism process，3008條)。

對(duì)新轉(zhuǎn)錄本的編碼區(qū)序列進(jìn)行預(yù)測(cè)，共鑒定出30 267個(gè)開(kāi)放閱讀框(ORF)，其中完整的ORF有20 046個(gè)。預(yù)測(cè)的ORF編碼蛋白序列為0～1000個(gè)氨基酸，其中以編碼100～200個(gè)氨基酸的ORF最多，數(shù)量為5530個(gè)，占比18.27%。

2.3 全長(zhǎng)轉(zhuǎn)錄本的結(jié)構(gòu)分析

基因在轉(zhuǎn)錄后可通過(guò)AS產(chǎn)生豐富的轉(zhuǎn)錄本，利用Astalavista v3.2將去冗余后的轉(zhuǎn)錄本進(jìn)行AS分析，共鑒定9360個(gè)AS事件。對(duì)AS類(lèi)型進(jìn)行分析，結(jié)果顯示AS事件被分為5種類(lèi)型，其中內(nèi)含子保留(IR)類(lèi)型數(shù)量最多，AS事件4791個(gè)，占比51.19%；其次為3′端可變剪接(A3SS)類(lèi)型，AS事件1889個(gè)，占比20.18%；第三種類(lèi)型為外顯子跳躍(ES)，AS事件1598個(gè)，占比17.07%；第四種類(lèi)型為5′端可變剪接(A5SS)，AS事件966個(gè)，占比10.32%；而互斥可變外顯子(MEE)類(lèi)型數(shù)量最少，AS事件116個(gè)，占比1.24%。使用去冗余前的一致序列進(jìn)行融合轉(zhuǎn)錄本預(yù)測(cè)，共鑒定1911個(gè)融合轉(zhuǎn)錄本。使用TAPIS pipeline v1.2.1對(duì)轉(zhuǎn)錄本進(jìn)行APA分析，共鑒定6436個(gè)基因存在poly(A)位點(diǎn)，其中3283個(gè)基因具有2個(gè)或多個(gè)poly(A)位點(diǎn)(圖1)。利用MEME v4.9.1對(duì)poly(A)位點(diǎn)上游50 bp的序列進(jìn)行分析，結(jié)果顯示在poly(A)剪切位點(diǎn)上游存在3個(gè)motif元件(CCCCC/TCCCCC/CCCTCC)。

圖1 poly(A)位點(diǎn)分析結(jié)果Fig.1 Results of poly(A) sites analysis

2.4 轉(zhuǎn)錄本的SSR檢測(cè)

使用MISA v1.0對(duì)長(zhǎng)度500 bp以上的轉(zhuǎn)錄本進(jìn)行SSR分析，共檢測(cè)37 839條轉(zhuǎn)錄本序列，其中包含SSR位點(diǎn)的序列16 114條，占檢測(cè)序列數(shù)的42.58%。在含有SSR位點(diǎn)的序列中，含2個(gè)及以上SSR位點(diǎn)的序列有5821條，含混合SSR位點(diǎn)的序列有3121條。在所有序列中共檢測(cè)到25 081個(gè)SSR位點(diǎn)，其中，單核苷酸重復(fù)(Mononucleotide)最多，SSR位點(diǎn)10275個(gè)，平均分布密度為106.06個(gè)/Mb；其次為二核苷酸重復(fù)(Dinucleotide)，SSR位點(diǎn)9431個(gè)，平均分布密度為78.95個(gè)/Mb；第三為三核苷酸重復(fù)(Trinucleotide)，SSR位點(diǎn)5155個(gè)，平均分布密度為53.38個(gè)/Mb，其他SSR位點(diǎn)類(lèi)型數(shù)量較少(圖2)。在重復(fù)單元中，二核苷酸重復(fù)排在前三位的重復(fù)單元為CT/AG、TC/GA和GA/TC，SSR位點(diǎn)數(shù)量分別為1517、1359和1041個(gè)，最少的為CG/CG，SSR位點(diǎn)數(shù)量為14個(gè)。三核苷酸重復(fù)排在前三位的重復(fù)單元為CAG/CTG、GGC/GCC和GCC/GGC，SSR位點(diǎn)數(shù)量分別為256、249和225個(gè)，最少的為T(mén)AC/GTA和AGT/ACT，SSR位點(diǎn)數(shù)量均僅各有1個(gè)。

c為混合SSR，2個(gè)SSR位點(diǎn)距離小于100 bp；c*為混合SSR，2個(gè)SSR位點(diǎn)無(wú)間隔：p1為單核苷酸重復(fù)；p2為二核苷酸重復(fù)；p3為三核苷酸重復(fù)；p4為四核苷酸重復(fù)；p5為五核苷酸重復(fù)；p6為六核苷酸重復(fù)c represented compound SSR，2 SSR loci less than 100 bp apart；c* represented compound SSR with no spacing between two SSR loci；p1 represented mono nucleotide repeats；p2 represented dinucleotide repeats；p3 represented trinucleotide repeats；p4 represented tetranucleotide repeats；p5 represented pentanucleotide repeats；p6 represented hexanucleotide repeats圖2 轉(zhuǎn)錄本的SSR類(lèi)型密度分布比較Fig.2 Comparison of the SSR type density distribution of transcripts

2.5 轉(zhuǎn)錄本的lncRNA預(yù)測(cè)

使用CPC、CNCI、CPAT和Pfam蛋白結(jié)構(gòu)域4種方法對(duì)lncRNA進(jìn)行鑒定，對(duì)4種分析結(jié)果取交集，共鑒定304個(gè)lncRNA(圖3)。根據(jù)lncRNA在參考基因組上的位置，將lncRNA分為4種類(lèi)型，其中，基因間區(qū)lncRNA(lincRNA)最多，有147個(gè)，占比48.40%，其次為正義lncRNA(sense-lncRNA)，有110個(gè)，占比36.20%，再次為反義lncRNA(antisense-lncRNA)，有39個(gè)，占比12.80%，而內(nèi)含子型lncRNA(intronic-lncRNA)數(shù)量最少，僅有8個(gè)，占比2.60%。基于物理位置(lncRNA與mRNA的位置關(guān)系)和互補(bǔ)序列(lncRNA與mRNA的堿基互補(bǔ)配對(duì))的方法對(duì)lncRNA的靶基因進(jìn)行預(yù)測(cè)，共預(yù)測(cè)靶基因2712個(gè)。

圖3 轉(zhuǎn)錄本的lncRNA鑒定Fig.3 lncRNA identification

2.6 轉(zhuǎn)錄因子分析

使用iTAK v1.7a對(duì)TF進(jìn)行預(yù)測(cè)，共預(yù)測(cè)到1608個(gè)TF，這些TF可分為28個(gè)家族，其中，MYB家族的TF最多，有395個(gè)，占比24.56%，其次為bHLH家族，有205個(gè)，占比12.75%(圖4)。

圖4 轉(zhuǎn)錄因子類(lèi)型分布情況比較Fig.4 Comparison of transcription factor type distribution

2.7 淀粉生物合成相關(guān)基因挖掘

通過(guò)KEGG代謝通路富集分析，共挖掘到淀粉生物合成相關(guān)基因14個(gè)，其中ADP-葡萄糖焦磷酸化酶(AGPase)6個(gè)，淀粉分支酶(SBE)3個(gè)，淀粉合成酶(SS)和淀粉磷酸化酶(SP)各2個(gè)，ADP-糖焦磷酸化酶(AspP)1個(gè)(表2)。從表2還可看出，不同的基因均檢測(cè)到1～12條轉(zhuǎn)錄本，其中，基因PB.9743檢測(cè)到的轉(zhuǎn)錄本數(shù)量最多，為12條，其次為基因PB.7121，檢測(cè)到轉(zhuǎn)錄本10條，基因PB.13928、PB.7537、PB.14152和PB.1105檢測(cè)到的轉(zhuǎn)錄本數(shù)量最少，均為1條。

表2 芋淀粉生物合成基因挖掘Table 2 Identification results of starch biosynthesis genes in taro

續(xù)表2 Continued table 2

對(duì)淀粉生物合成的基因進(jìn)行AS分析，結(jié)果(圖5和表3)顯示，6個(gè)AGPase基因中有3個(gè)發(fā)生AS事件，包括5種AS類(lèi)型，其中，以基因PB.7121發(fā)生的AS類(lèi)型最多，包含IR、A5SS、A3SS和MEE 4種類(lèi)型；2個(gè)SP基因均發(fā)生AS事件，其中基因PB.9743含有4種AS類(lèi)型，分別為A5SS、A3SS、IR和MEE，而基因PB.8912只有A3SS類(lèi)型的AS事件；SBE基因中，只有基因PB.9363發(fā)生A3SS和IR 2種類(lèi)型的AS事件；SS和AspP基因不發(fā)生AS事件。此外，對(duì)淀粉生物合成的基因進(jìn)行APA分析結(jié)果(表2)表明，14個(gè)基因中有10個(gè)基因存在poly(A)位點(diǎn)，其中7個(gè)基因具有2個(gè)及2個(gè)以上poly(A)位點(diǎn)，尤其以基因PB.729的poly(A)位點(diǎn)數(shù)量最多(6個(gè))，其次為基因PB.7121和PB.9743(5個(gè))，基因PB.11557、PB.7516和PB.1105的poly(A)位點(diǎn)數(shù)量最少，均為1個(gè)。

表3 芋淀粉生物合成相關(guān)基因AS事件統(tǒng)計(jì)Table 3 Statistical of AS events of genes related to taro starch biosynthesis

3 討 論

芋的基因組雖已公布，但并未深入分析其轉(zhuǎn)錄本的結(jié)構(gòu)特征。Wang等[8]利用三代全長(zhǎng)轉(zhuǎn)錄組測(cè)序技術(shù)糾正毛竹基因組錯(cuò)誤注釋基因2241個(gè)，鑒定新轉(zhuǎn)錄本35 447條，尚驍堯等[15]通過(guò)三代全長(zhǎng)轉(zhuǎn)錄組測(cè)序鑒定蒺藜苜蓿新基因7209個(gè)，新轉(zhuǎn)錄本52 636條。本研究結(jié)果與上述研究結(jié)果相似，通過(guò)三代全長(zhǎng)轉(zhuǎn)錄組測(cè)序獲得芋全長(zhǎng)轉(zhuǎn)錄本序列38 043條，鑒定新基因1878個(gè)，新轉(zhuǎn)錄本31 058條。說(shuō)明三代全長(zhǎng)轉(zhuǎn)錄組測(cè)序技術(shù)在完善植物基因組功能注釋中可發(fā)揮重要作用。

由于三代全長(zhǎng)轉(zhuǎn)錄組測(cè)序技術(shù)無(wú)需對(duì)RNA進(jìn)行打斷和拼接，因此在轉(zhuǎn)錄本的結(jié)構(gòu)分析方面具有極大優(yōu)勢(shì)。AS在植物生長(zhǎng)發(fā)育、抗逆響應(yīng)方面具有重要作用[6]。本研究通過(guò)三代全長(zhǎng)轉(zhuǎn)錄組測(cè)序共鑒定9360個(gè)AS事件，與Wang等[8]對(duì)毛竹、尚驍堯等[15]對(duì)蒺藜苜蓿、Li等[27]對(duì)黃芪、Wang等[28]對(duì)海島棉的鑒定結(jié)果相似，說(shuō)明三代全長(zhǎng)轉(zhuǎn)錄組測(cè)序可有效鑒定轉(zhuǎn)錄本的AS。APA通過(guò)產(chǎn)生不同長(zhǎng)度3′UTR或不同編碼序列的轉(zhuǎn)錄本來(lái)提高轉(zhuǎn)錄的復(fù)雜性，從而通過(guò)多種機(jī)制調(diào)控植物的基因表達(dá)[21]。本研究結(jié)果表明，6436個(gè)基因存在poly(A)位點(diǎn)，與Wang等[8]對(duì)毛竹、尚驍堯等[15]對(duì)蒺藜苜蓿、Abdel-Ghany等[21]對(duì)高粱的研究結(jié)果相似，說(shuō)明三代全長(zhǎng)轉(zhuǎn)錄組測(cè)序也可有效鑒定轉(zhuǎn)錄本的APA，進(jìn)一步說(shuō)明三代全長(zhǎng)轉(zhuǎn)錄組測(cè)序技術(shù)在轉(zhuǎn)錄本結(jié)構(gòu)鑒定中發(fā)揮重要作用。

字體加粗代表發(fā)生AS事件的基因；IR：內(nèi)含子保留；ES：外顯子跳躍；A5SS：5′端可變剪接；A3SS： 3′端可變剪接；MEE：互斥可變外顯子；AGPase：ADP-葡萄糖焦磷酸化酶；SS：淀粉合成酶；SBE：淀粉分支酶；AspP：ADP-糖焦磷酸化酶；SP：淀粉磷酸化酶Bolded font represented the gene in which the AS event occurred；IR：Intron retention；ES：Exon skipping；A5SS：Alternative 5′splice site；A3SS：Alternative 3′ splice site；MEE：Mutually exclusive exon；AGPase：ADP-glucose pyrophosphorylase；SS：Starch synthase；SBE：Starch branching enzyme；AspP：ADP-sugar pyrophosphatase；SP：Starch phosphorylase圖5 芋淀粉生物合成途徑的AS事件示意圖Fig.5 Schematic diagram of AS events in the taro starch biosynthesis pathway

SSR標(biāo)記由于操作技術(shù)簡(jiǎn)單、穩(wěn)定性好及具有共顯性等特點(diǎn)，已廣泛應(yīng)用于開(kāi)展植物遺傳多樣性分析、DNA指紋圖譜構(gòu)建和遺傳作圖[29]。You等[10]、Wang等[11]分別從芋二代轉(zhuǎn)錄組測(cè)序中鑒定了5278和11 363個(gè)SSR位點(diǎn)，Dong等[12]使用RAD-seq技術(shù)從芋中鑒定了4438個(gè)SSR位點(diǎn)，而本研究通過(guò)三代全長(zhǎng)轉(zhuǎn)錄組測(cè)序可檢測(cè)到25 081個(gè)SSR位點(diǎn)，可見(jiàn)，利用二代轉(zhuǎn)錄組測(cè)序和RAD-seq技術(shù)2種方法檢測(cè)到SSR的位點(diǎn)數(shù)明顯少本研究檢測(cè)的SSR位點(diǎn)數(shù)量，說(shuō)明三代測(cè)序技術(shù)在鑒定SSR位點(diǎn)方面優(yōu)于上述2種方法，可鑒定更多的SSR位點(diǎn)。

lncRNA是一類(lèi)長(zhǎng)度大于200 bp但缺乏編碼能力的轉(zhuǎn)錄本。已有許多研究表明，lncRNA在植物生長(zhǎng)發(fā)育及響應(yīng)生物和非生物脅迫方面發(fā)揮重要作用[30]。本研究鑒定了304個(gè)lncRNA，其中l(wèi)incRNA的數(shù)量最多，這與對(duì)其他植物的研究報(bào)道結(jié)果相似[15，31]，說(shuō)明三代全長(zhǎng)轉(zhuǎn)錄組測(cè)序可有效鑒定lncRNA。TF是一類(lèi)具有特殊結(jié)構(gòu)的蛋白，通過(guò)與靶基因上游啟動(dòng)子特異性結(jié)合而調(diào)控靶基因的表達(dá)，在植物生長(zhǎng)發(fā)育中發(fā)揮重要作用[32]。本研究共鑒定了1608個(gè)TF，其中MYB家族的TF最多，這些鑒定的TF將為后續(xù)研究芋基因功能提供充實(shí)的數(shù)據(jù)。

本研究從芋淀粉生物合成的基因中共挖掘到5種酶的14個(gè)基因和60條轉(zhuǎn)錄本。在其他研究中，如Liu等[33]通過(guò)二代轉(zhuǎn)錄組測(cè)序挖掘到淀粉生物合成的AGPase、SS、SBE和SP基因共9個(gè)，Dong等[3]通過(guò)比較轉(zhuǎn)錄組在球莖發(fā)育過(guò)程中富集到淀粉生物合成的基因(AGPase、SS、SBE和SP)共10個(gè)，說(shuō)明三代全長(zhǎng)轉(zhuǎn)錄組測(cè)序可挖掘到更多的基因數(shù)量。此外，本研究檢測(cè)到基因PB.9743的轉(zhuǎn)錄本數(shù)量多達(dá)12條，說(shuō)明三代全長(zhǎng)轉(zhuǎn)錄組測(cè)序可鑒定同一基因的不同轉(zhuǎn)錄本，而二代轉(zhuǎn)錄組測(cè)序較難實(shí)現(xiàn)，這為研究轉(zhuǎn)錄本結(jié)構(gòu)信息提供了可靠依據(jù)。基于此，本研究對(duì)淀粉生物合成通路的基因進(jìn)行AS和APA分析。AS分析結(jié)果顯示，AGPase基因、SBE基因和SP基因均發(fā)生了AS事件，包括IR、ES、A5SS、A3SS和MEE 5種類(lèi)型，說(shuō)明芋淀粉生物合成過(guò)程中存在豐富的轉(zhuǎn)錄調(diào)控。對(duì)其他植物的研究也獲得類(lèi)似的結(jié)果，如Li等[34]對(duì)草莓進(jìn)行AS分析，發(fā)現(xiàn)AS可影響草莓發(fā)育過(guò)程相關(guān)基因和TF的功能；孫銘陽(yáng)等[35]對(duì)穿心蓮內(nèi)酯前體合成途徑基因進(jìn)行AS分析，結(jié)果表明穿心蓮內(nèi)酯前體合成的兩條途徑(質(zhì)體MEP途徑和細(xì)胞質(zhì)MVA途徑)均發(fā)生了AS事件，其中有1個(gè)基因產(chǎn)生了6個(gè)可變啟動(dòng)子式的IR亞型。本研究對(duì)芋淀粉生物合成的基因進(jìn)行APA分析，結(jié)果發(fā)現(xiàn)14個(gè)基因中有10個(gè)基因存在poly(A)位點(diǎn)，其中7個(gè)基因具有2個(gè)及以上poly(A)位點(diǎn)，說(shuō)明APA參與調(diào)控芋淀粉的生物合成。在其他植物的研究中，Simpson等[36]證實(shí)APA在擬南芥開(kāi)花過(guò)程中具有重要的調(diào)控作用；Abdel-Ghany等[21]對(duì)高粱進(jìn)行干旱處理，結(jié)果顯示同一個(gè)基因產(chǎn)生了不同的APA；Wang等[8]在毛竹中鑒定了11個(gè)纖維素合酶基因(CesA)、11個(gè)纖維素合成酶相似基因(CsI)和2個(gè)木質(zhì)素基因由APA調(diào)控，說(shuō)明APA可能參與調(diào)控細(xì)胞壁結(jié)構(gòu)和次生細(xì)胞壁的形成。

4 結(jié) 論

通過(guò)三代全長(zhǎng)轉(zhuǎn)錄組測(cè)序分析能獲取芋全長(zhǎng)轉(zhuǎn)錄本的序列和結(jié)構(gòu)信息，并挖掘到參與芋淀粉生物合成相關(guān)的基因14個(gè)，轉(zhuǎn)錄本60條，可為后續(xù)闡明芋淀粉生物合成分子機(jī)制及深入利用芋基因資源提供科學(xué)依據(jù)。