基于家蠶中腸RNA-seq數據的新基因發掘及初步分析

周凱月, 唐 健, 李玉霞, 郝長富, 徐安英, 2*

(1.江蘇科技大學生物技術學院，江蘇鎮江 212018；2.中國農業科學院蠶業研究所，江蘇鎮江 212018)

我國是繭絲綢生產和出口大國，家蠶作為重要的吐絲昆蟲，具有很高的經濟價值，在生物反應器和模式生物中也顯現出廣闊的應用前景。隨著家蠶研究的深入，越來越多的新基因被發掘出來，一些新基因往往發揮著極為重要的作用，甚至是某些疑難問題的突破點，人們對新基因功能的研究也越來越重視。

轉錄組研究是一個發掘功能基因的重要途徑，與基因組學相比，轉錄組學只研究被轉錄的基因，研究范圍縮小，針對性更強[1]，越來越多地被運用到基因功能的相關研究中。轉錄組是特定組織或細胞在某一發育階段或功能狀態下轉錄出來的所有RNA的總和[2]，主要包括mRNA和非編碼RNA(ncRNA)。轉錄組測序即通過第二代高通量測序技術對特定組織或細胞的轉錄產物(主要是全部mRNA)反轉錄后測序并對其進行生物信息學分析的技術[3]，是當前在全基因組水平上研究基因表達模式的主要技術[4]，目前該技術已被廣泛應用于生物信息學研究的多個領域[5]。筆者擬采用RNA-seq技術對所構建的家蠶中腸的轉錄組進行測定，并在基因組水平上進行轉錄組分析，對新基因進行初步發掘，以期為新基因的功能鑒定基礎。

1 材料與方法

1.1材料試驗所用家蠶品種為抗BmNPV家蠶品種QFN和常規品種QF，來源于中國農業科學院蠶業研究所蠶資源中心課題組。

1.2方法

1.2.1轉錄組測序。首先進行樣品檢測、構建RNA文庫，文庫質控合格后，采用HiSeq2500進行高通量測序， Illumina HiSeq2500高通量測序獲得Reads或堿基信息，篩選除去冗余后得到Clean Reads，通過solexa QA軟件對其進行質量檢測可得到高質量的Clean Reads[6]。

1.2.2轉錄組數據比對。對于Clean Reads需要用高效的序列比對軟件TopHat2將其與參考基因組進行序列比對，得到Mapped Reads。比對效率可以直接反映出轉錄組數據的利用率[6]。

1.2.3新基因分析。通過Cufflinks軟件對Mapped Reads進行組裝，將得到的序列與參考基因組注釋信息進行比對，尋找未知的新基因。再利用Blast軟件對新基因進行功能注釋，然后利用各個數據庫分別對新基因的NR注釋信息、COG功能注釋及其分類、KEGG注釋通路進行分析，獲得新基因的相關注釋信息。

1.2.4部分新基因的克隆分析。以家蠶抗性品種QFN的中腸組織cDNA為模板，隨機挑選3個新基因(Silkworm New Gene 1、Silkworm New Gene 5、Silkworm New Gene 20)設計引物，進行RT-PCR擴增。PCR反應體系為(25 μL)：1 μL cDNA模板，2.5 μL 10×PCR Buffer，2 μL 10 mmol/mL dNTP，1 μL 20 pmol/mL上下游引物，0.3 μL 5 U/μL ExTaqDNA聚合酶，加雙蒸水補充至25 μL。PCR反應后檢測擴增片段的大小，與目的片段大小是否相符。所需引物序列見表1。

采用SanPrep柱式DNA膠回收試劑盒進行目的片段的回收和純化。試驗前確認Wash Solution中是否加入乙醇，將膠回收純化的PCR產物連接到pMDTM18-T Vector上，用冷熱刺激的方法轉化感受態細胞，過夜培養后挑取獨立菌落進行培養，再進行菌液PCR擴增判斷是否為目的條帶，挑選陽性菌液送生工生物工程(上海)股份有限公司測序鑒定。

2 結果與分析

2.1測序數據為研究家蠶感染BmNPV后蠶體內基因表達調控情況，以家蠶抗性品種QFN和常規品種QF中腸組織為材料進行轉錄組測序分析(表2)。參考基因組組裝能否滿足信息分析的需求，可以通過轉錄組數據與參考基因組序列的比對結果評估(表3)。

經篩選，2個文庫共獲得12.8 Gb Clean Data，QF、QFN的≥Q30的堿基百分比分別為86.03%、87.01%。

表1 部分新基因RT-PCR引物序列

表2 樣品測序數據統計

表3 Clean Data與參考基因組比對結果統計

文庫中比對到參考基因組上的Reads在Clean Reads的效率達69.36%和70.24%，其中比對到參考基因組唯一位置的Reads在Clean Reads中所占比例分別為58.66%、60.75%。

2.2新基因分析

2.2.1新基因發掘及基因結構的分析。對測序得到的序列進行拼接和組裝，與原有的一些基因組注釋信息進行比對，尋找未被注釋的新基因。該研究過濾掉編碼的肽鏈過短(少于50個氨基酸殘基)或只包含單個外顯子的序列，得到了788個新基因。如Silkworm New Gene1007位于nscaf2794基因序列172～4 311的正鏈上，包含5個外顯子；Silkworm New Gene1008位于nscaf2795基因序列2 402 828～2 404 107的正鏈上，包含3個外顯子，部分新基因的文件見表4。

表4 部分新基因的文件

注： SNG為Silkworm New Gene簡寫

Note: SNG is short for silkworm new gene

2.2.2新基因功能注釋。使用Blast軟件將發掘的新基因分別同NR[7]，Swiss-Prot[8]、GO[9]、COG[10]、KEGG[11]數據庫進行序列比對[12]，結果發現788個新基因中746個得到了注釋，各數據庫得到注釋的基因數分別為746、418、443、187、224個。新基因功能注釋結果統計詳見表5。

表5 新基因功能注釋結果統計

同時發現，其中746(94.67%)個新基因與數據庫中匹配到的序列具有顯著的相似性(E

圖1 注釋基因E值分布Fig.1 Distribution of E-value for annotated genes

2.2.3與注釋基因匹配的物種分布。利用BlastX將組裝出來的unigene序列與NR數據庫進行比對后，共找到746個unigene與其他近緣生物的已知基因具有不同程度的同源性，746個注釋的基因中，有611條(81.90%)基因與家蠶(Bombyxmori)序列同源，99條(13.27%)與黑脈金斑蝶(Danausplexippus)序列同源、4條(0.54%)與玉帶鳳蝶(Papiliopolytes)序列同源、4條(0.54%)與赤擬谷盜(Triboliumcastaneum)序列同源、4條(0.54%)與致倦庫蚊(Culexquinquefasciatus)序列同源、3條(0.40%)與柑橘鳳蝶(Papilioxuthus)序列同源、2條(0.27%)與印度跳蟻(Harpegnathossaltator)序列同源、2條(0.27%)與佛羅里達弓背蟻(Camponotusfloridanus)序列同源、2條(0.27%)與畢氏粗角猛蟻(Cerapachysbiroi)序列同源，僅有15條(2.01%)與其他物種序列相匹配(圖2)。

2.2.4新基因GO富集分析。利用Blast2Go軟件對篩選到的基因進行GO富集分析，結果顯示，基因主要注釋到細胞組分、分子功能和生物學過程3個分支中，分別有947、551和1 777個(圖3)。在細胞組分模塊中(圖3A)，注釋到細胞(cell)、細胞部分(cell part)的基因數目較多，分別占19.6%和20.0%；在分子功能模塊中(圖3B)，注釋到黏合(binding)和催化活動(catalytic activity)的基因數目較多，分別占41.0%和37.7%。在生物學過程模塊中(圖3C)，注釋到細胞過程(cellular process)和單組織過程(single organismprocess)及代謝過程(metabolic process)的基因數目較多，分別占13.8%、13.6%和13.9%。

圖2 與注釋基因匹配的物種分布Fig.2 Distribution of species match to the annotated genes

2.2.5新基因在家蠶基因連鎖群上的分布。將新基因的Locus與家蠶基因連鎖群進行比對發現，有545個新基因分布在不同的染色體上，且在28條染色體上都有分布(圖4)，在18號染色體上分布數量最多，有58個，其中在nscaf2902上分布有54個；在15號染色體上分布有40個，其中在nscaf2888上分布有36個；在26號染色體上分布最少，僅有7個。

2.2.6新基因通路富集分析。通過KEGG分析，對新基因進行通路富集分析發現，共有 198個新基因注釋到KEGG數據庫中，分布于85條已知的通路中。映射基因最多的5個通路分別為剪接體(splicesome)(ko03040、12條)、RNA轉運(RNA transport)(ko03013、9條)、真核細胞核糖體合成(ribosome biogenesis in eukaryotes)(ko03008、6條)、過氧物酶體(peroxisome)(ko04146、5條)、內吞(endocytosis)(ko04144、5條)。映射到的信號通路見圖5。

2.2.7新基因COG數據庫功能注釋。將新基因與COG數據進行比對，并進行功能注釋與分類，結果發現共有22個類別里的258個新基因得到了注釋(圖6)，其中，一般功能(general function prediction only)的基因占總體的20.93%，所占比例最大；復制、重組和修復(replication, recombination and repair )、碳水化合物的運輸和代謝(carbohydrate transport and metabolism)、氨基酸的運輸和代謝(amino acid transport and metabolism)3個類別共居第2位，占總體的7.75%；轉錄(transcription)占總體的7.36%，其余分類的基因數較少，其中核結構(nuclear structure)、細胞運動(cell motility)、真核細胞的細胞外結構(extracellular structures)3個類別里均無新基因出現。

2.2.8部分新基因的克隆分析。以秋豐N中腸組織的cDNA為模板，隨機挑選3個新基因設計引物擴增目的基因，獲得特異性片段，與預期片段大小相符，經測序結果與參考序列比對后發現編碼序列高度相似。

3 討論

該研究基于所選參考基因組序列，共發掘788個新基因，通過生物信息學軟件將發掘的新基因與NR、Swiss-Prot、GO、COG及KEGG數據庫進行序列比對，共獲得746個新基因的注釋信息；746個注釋的基因中，其中611條(81.90%)基因與家蠶序列同源；利用Blast2Go軟件對篩選到的基因進行GO富集分析，結果顯示，基因主要注釋到細胞組分、分子功能和生物學過程3個分支中，分別有947、551和1 777個，注釋到細胞部分(cell part)、黏合(binding)及代謝過程(metabolic process)的基因數目最多；將新基因的Locus與家蠶基因連鎖群進行比對發現，有545個新基因分布在不同的染色體上，且在28條染色體上都有分布；通過KEGG分析對新基因進行通路富集分析發現，共有 198個新基因注釋到KEGG數據庫中，分布于85條已知的通路中；將新基因與COG數據進行比對，并進行功能注釋與分類，結果發現共有22個類別里的258個新基因得到了注釋。新基因與各個數據的序列比對結果，進一步證實了新基因的存在，該研究對新基因的功能做了初步分析，關于新基因具體的功能還需要進一步研究。

注：A.細胞組分；B.分子功能；C.生物學過程Note：A.Cellular components；B.Molecular function；C.Biological processes圖3 注釋基因在GO中的分類情況Fig.3 Classification of annotated genes in GO

圖4 新基因在家蠶基因連鎖群的分布Fig.4 Distribution of new gene in silkworm gene linkage group

圖5 KEGG通路富集分布Fig.5 KEGG pathways distribution

圖6 COG數據庫中功能注釋的unigenes分類Fig.6 COG function classification of unigenes