羅通宇,曾俊源,呂奕琳,周峻崗,余 垚,呂 紅
(1.復旦大學 遺傳工程國家重點實驗室,上海 200438; 2.復旦大學 生命科學學院 上海工業(yè)菌株工程技術研究中心,上海 200438)
馬克斯克魯維酵母(Kluyveromycesmarxianus)在分類上位于真菌界-子囊菌門-酵母綱-酵母目-酵母科-克魯維屬。屬名Kluyveromyces是為紀念荷蘭微生物學家Albert Jan Kluyver(1888—1956);種名marxianus是為紀念其發(fā)現者Marx,他于1888年首次從葡萄中分離得到該菌種[1]。之后,人們發(fā)現K.marxianus在發(fā)酵乳制品如奶酪和開菲爾粒(Kefir grains)中大量存在,也因此為人們所熟知[2]。除了安全可食用外,K.marxianus還具備耐高溫、生長快速、代謝底物廣等多種優(yōu)良特性,有廣泛的工業(yè)應用潛力。但是目前,由于人們對K.marxianus的遺傳學和生理學的認識有限,阻礙了其進一步應用[3]。因此,研究人員需要從基因組入手,更細致地研究其生理表型與基因突變、基因表達的內在聯系,增加對K.marxianus的認識。
現有對K.marxianus基因組的研究主要集中在基因組測序組裝與注釋方面,目前仍然缺乏對K.marxianus種內分化、核心基因鑒定的研究。截至2021年3月,美國國家生物信息中心(NCBI)上已有17株K.marxianus的基因組信息,其中5株組裝出完整的染色體,在此基礎上,有3株K.marxianus完成了基因組注釋工作。2015年,Lertwattanasakul等人使用鳥槍法完成了K.marxianusDMKU3-1042(簡稱DMKU3-1042)的基因組測序、組裝與注釋工作[4]。同年4月,Inokuma等人完成了對K.marxianusNBRC1777(簡稱NBRC1777)的基因組測序、組裝與注釋工作[5]。2019年10月,本實驗室完成了K.marxianusFIM1(簡稱為FIM1)的基因組測序、組裝與注釋工作[6],并將結果提交到了NCBI,編號為GCA_001854445.2。但后續(xù)應用FIM1的基因組信息進行深入地遺傳學研究時發(fā)現,對FIM1部分位點進行Sanger法測序,重測序結果與提交到NCBI上的對應序列不符。這種問題的存在不利于后續(xù)對其進行遺傳學研究與遺傳工程改造,所以需要對FIM1的基因組序列及注釋進行校驗更新。
綜上,研究從FIM1出發(fā),使用重測序的DNA-seq數據校驗了FIM1的基因組序列,隨后補充了FIM1的注釋,并使用比較基因組學的方法進一步完善了FIM1的基因組信息。將DNA-seq比對到之前已完成組裝的FIM1基因組上,根據比對結果校驗了FIM1的基因組序列。隨后,研究評估了FIM1基因組質量,補充了FIM1的非編碼RNA與次級代謝產物基因簇注釋信息。進一步,研究對FIM1進行了比較基因組學研究。對FIM1與DMKU3-1042、乳酸克魯維酵母KluyveromyceslactisNRRL Y-1140、釀酒酵母SaccharomycescerevisiaeS288C共3株酵母進行基因組共線性分析,劃定了FIM1在酵母分化過程中的基因排列順序保守的區(qū)域。對FIM1與NCBI來源的11株K.marxianus進行種內基因組比較,分析了FIM1與其他完成組裝的K.marxianus的全基因組序列差異,構建了12株酵母的同源基因家族,并依此繪制了種內系統進化樹,構建了K.marxianus的泛基因組,鑒定了FIM1菌株特性基因及其涉及的代謝通路,為FIM1遺傳學的研究與遺傳工程的改造提供了數據支持。
本文中所有基因組與注釋數據均下載自NCBI。在NCBI中的“Genome”條目下輸入“Kluyveromycesmarxianus”查找K.marxianus所有菌株序列信息,下載12株基因組組裝到Contig水平的K.marxianus基因組與注釋信息。乳酸克魯維酵母(K.lactis)、畢赤酵母(Pichiapastoris)、釀酒酵母(S.cerevisiae)的基因組與蛋白信息均源自NCBI上對應物種參考基因組。詳細信息如表1所示。
表1 本研究使用菌株的信息Tab.1 Information on strains used in this study
1.2.1 FIM1基因組序列校驗
研究使用抽提好的DNA基于Illumina Truseq Nano試劑構建文庫,隨后使用Illumina Hiseq X測序平臺進行重新測序,測序深度為200×。測序結果經FastQC(https:∥www.bioinformatics.babraham.ac.uk)質量檢驗后,使用HISAT2[7](version 2.1)將DNA-seq數據比對到基因組上,將比對結果與FIM1基因組序列輸入到inGAP[8](version 1)中,根據比對結果,刪除無reads比對的區(qū)域,生成新的FIM1序列。
1.2.2 FIM1基因組組裝質量評估
使用BUSCO[9](version 5.0)中的酵母目保守基因庫(saccharomycetales_odb10)與真核生物保守基因庫(eukaryota_odb10)分別評估FIM1的基因組質量,運行基因組評估默認參數。用于評估的酵母目保守基因庫包含2137個基因,真核生物保守基因庫包含255個基因。
1.2.3 FIM1的補充注釋
將Rfam[10](verison 14)及其依賴軟件下載到本地后,使用Rfam的默認參數預測FIM1的非編碼RNA序列位置與功能。
將antiSMASH[11](version 5.1.1)及其依賴軟件下載到本地后,運行antiSMASH,輸入FIM1基因組序列,分類學的運行參數選擇“fungi”,其他參數為全功能運行模式默認參數。
1.2.4 基因組共線性分析
使用Python版的MCscan[12](jcvi工具包)進行共線性分析,比較FIM1與DMKU3-1042、K.lactisNRRL Y-1140、S.cerevisiaeS288C 3株酵母的基因組共線性情況。使用jcvi.formats.gff模塊,輸入基因組序列與注釋文件獲得基因組bed格式的注釋信息。之后輸入蛋白序列與bed文件進行共線性分析,將分析結果整理后,使用jcvi.graphics.karyotype工具包生成共線性分析結果。
1.2.5K.marxianus種內基因組序列比較
使用Mummer[13](version 4)將FIM1與其他完成組裝的K.marxianus基因組序列進行比較。運行Mummer的dnadiff參數,輸入FIM1與待比較的K.marxianus基因組序列,在report文件中獲取基因組差異比較結果。
1.2.6 系統進化樹繪制及泛基因組分析
使用Orthofinder[14]進行同源基因家族分析,輸入菌株的蛋白序列到Orthofinder中,然后使用DIAMOND[15]方法比對蛋白序列并構建同源基因家族,極大似然法構建系統發(fā)育樹,建樹的推理方法使用默認的fasttree軟件,其他運行參數均為Orthofinder的默認參數。由于12株K.marxianus物種的分化時間尚不清楚,所以研究中增加K.lactis、P.pastoris、S.cerevisiae3株已經確定分化時間的酵母對進化樹進行標定校正,3個物種參考分化時間通過Timetree[16]網站獲取。然后使用CAFE[17](version 4.2.1)分析12株酵母在各個時間節(jié)點的基因收縮與擴增情況。隨后使用r8s軟件[18]繪制K.marxianus種內系統進化樹,在iTOL[19]中將用于時間標定的3株酵母的信息刪除。最后使用PanGP[20]預測并構建K.marxianus的泛基因組。
除FIM1蛋白序列外,本研究用到的蛋白序列由以下3種方法獲得。(1) 從NCBI上下載NBRC1777、DMKU3-1042、K.lactisNRRL Y-1140、S.cerevisiaeS288C、ATCC 28485共5株酵母的蛋白序列。(2) 輸入基因組序列與注釋信息到Cufflinks[21]中,使用gffreads模塊提取Olga1、Olga2的蛋白序列。(3) 使用AUGUSTUS[22](version 3.3)軟件預測基因,運行getAnnoFasta.pl腳本提取CBS6556、NRRL Y-6860、B0399、DMB1、IIPE453、100656-19、KCTC 17555共7株K.marxianus的蛋白序列。為了優(yōu)化預測結果,運行AUGUSTUS前,使用DMKU3-1042的基因組信息對AUGUSTUS進行訓練。
1.2.7 FIM1特有基因的GO富集分析
通過Orthofinder中生成的Orthogroups_UnassignedGenes.csv結果文件,獲得FIM1的特有基因信息。隨后使用R語言中的clusterprofiler[23]程輯包配合org.Sc.sgd.db物種數據庫(http:∥www.bioconductor.org/packages/release/data/annotation/html/org.Sc.sgd.db.html),輸入FIM1中有功能注釋的特有基因,進行GO富集分析。
本研究使用HISAT2將重測序的DNA-seq數據比對到FIM1基因組上,inGAP可視化比對結果(圖1,見第306頁),根據比對結果刪除DNA-seq無法覆蓋的序列片段61處,其中基因外部片段44處,基因內部片段17處,無基因啟動子或終止子被刪除,刪除序列的詳細信息見表2(見第306頁)。FIM1與校驗前相比,共刪除序列4 910 bp,新序列長度為10 909 543 bp。為了評估校驗是否提高了FIM1基因序列的精度,使用BUSCO評估校驗前后的FIM1基因組質量。采用BUSCO構建的真核生物保守基因庫(eukaryota_odb10)與酵母目保守基因庫(saccharomycetales_odb10)對FIM1基因組質量進行了評估(表3,見第307頁),結果顯示校驗后基因組質量未發(fā)生改變。評估顯示,FIM1基因組序列可覆蓋真核生物保守基因庫中91.4%的基因,與酵母目保守基因庫中99.5%的基因。2019年,McCarthy使用BUSCO的酵母目保守基因庫(saccharomycetales_odb9)評估了NCBI上100株S.cerevisiae酵母的基因組組裝質量,其中得分最高的酵母的基因組序列可覆蓋酵母目保守基因庫中98.8%的基因[24]。FIM1基因組組裝質量比該值高0.7%,表明FIM1基因組完整度較好,可以用于后續(xù)實驗與分析。
圖1 FIM1基因組序列校驗Fig.1 Genome sequence verification of FIM1(a) 圖中紅色和藍色線條分別代表正義鏈和負義鏈雙端測序的DNA-seq匹配到FIM1的CP015054.1號染色體上的結果,下方坐標指標明了染色體位置,灰色表示對應位點的DNA-seq比對結果不佳。(b) 將灰色部分的DNA-seq比對結果進行可視化,展示了具體的reads與對應位點的比對關系。
表2 FIM1刪除序列的位置信息Tab.2 Position of FIM1 deletion sequence
表3 FIM1基因組質量評估Tab.3 Quality assession of the FIM1 genome
為了進一步評估校驗結果,研究以基因內部有序列的刪除的HHT1(Histone H3)基因為例,對校驗結果進行再評估。HHT1基因編碼的蛋白為H3組蛋白,是一種高度保守的蛋白,在FIM1與DMKU3-1042中均有兩個拷貝,校驗前FIM1中H3組蛋白的長度為158 aa(amino acids),而校驗后長度為136 aa。將校驗前后FIM1的H3組蛋白序列與DMKU3-1042的H3組蛋白序列進行比對,結果顯示校驗后FIM1的H3序列與DMKU3-1042的H3序列匹配率為100%,而校驗前FIM1中的H3組蛋白序列與DMKU3-1042相比有22 aa重復(圖2),說明DNA-seq數據校驗提高了FIM1基因組序列的質量。
圖2 FIM1與DMKU3-1042的H3蛋白序列比對Fig.2 Sequence alignment of H3 between FIM1 and DMKU3-1042KLMA_60129、KLMA_80119為DMKU3-1042的H3蛋白序列,FIM1_4163、FIM1_4900為FIM1的H3蛋白序列。
使用Rfam注釋FIM1基因組的非編碼RNA,結果顯示FIM1基因組中含有tRNA 154個,snRNA(Small nucleolar RNA)55個,rRNA 8個。將Rfam預測的結果與FIM1原有注釋合并,整合后預測FIM1含有tRNA 172個,其他非編碼RNA數目與Rfam預測結果一致。
使用antiSMASH注釋FIM1次級代謝產物基因簇信息,結果顯示FIM1有一個次級代謝產物基因簇,為非核糖體肽合成酶簇片段(NRPS-like),基因簇總長44 161 bp(表4)。
表4 非核糖體肽合成酶簇片段信息Tab.4 Information of non-ribosomal peptide synthetase cluster fragment
為了探究FIM1在物種分化過程中基因組的變化,本文分別選取了與FIM1同種的DMKU3-1042、同屬的K.lactisNRRL Y-1140、同科的S.cerevisiaeS288C與FIM1進行基因組共線性分析。結果顯示FIM1與DMKU3-1042在基因組共線性程度極高(圖3(a)),共97.8%的FIM1序列與DMKU3-1042有同源性。FIM1與DMKU3-1042相比,1號與6號染色體的序列反向,序列無大片段缺失、重復等現象。FIM1與K.lactisNRRL Y-1140(KL)進行基因組共線性分析,結果顯示兩者基因組共線性程度較高(圖3(b)),共92.00%的FIM1序列與KL有同源性,兩者分化引起了大量的片段重排,未發(fā)現大規(guī)?;驈椭片F象,結果提示了基因組重排可能是克魯維屬酵母的一種重要進化方式。FIM1與S.cerevisiaeS288C(SC)的共線性分析結果顯示兩者基因組共線性程度不高(圖3(b)),共18.91%的FIM1序列與SC有同源性,FIM1的1號和8號染色體上有兩個基因片段,在SC上發(fā)生了基因片段復制現象,這一結果是對釀酒酵母曾發(fā)生過全基因組復制現象的支持。共線性分析的結果劃定了FIM1基因組上基因排序保守的區(qū)域,結合其他3株已有的基因功能注釋,可以對FIM1的基因功能注釋進行再校驗。
圖3 FIM1基因組共線性分析Fig.3 Genome collinearity analysis of FIM1(a) FIM1與K.marxianus DMKU3-1042(DMKU)的共線性分析;(b) FIM1分別與K.lactis NRRL Y-1140(KL)、S.cerevisiae S288C(SC)的共線性分析。圖中數字代表酵母的染色體編號。
除了FIM1外,NCBI上有DMKU3-1042、NBRC1777、CBS6556、NRRL Y-6860共4株K.marxianus組裝出完整染色體信息。本文使用Mummer將這4株酵母與FIM進行基因組間序列比對分析(表5)。
表5 FIM1與其他4株K.marxianus的基因組差異Tab.5 Genomic variation of 4 strains of K.marxianus compared with FIM1
分析結果顯示,與FIM1相比,4株酵母基因組差異情況總體比較接近,基因組上產生的單堿基插入/缺失(Indels)的數目在1 786~1 996之間;易位(Translocations)的數目在12~18之間;倒位(Inversions)的數目在6~8之間;大片段插入(Insertions)的數目在194~271之間。其中只有單堿基變異(SNP)數目出現了明顯差異,與FIM1相比,NRRL Y-6860出現了41 418個SNP變異,而其余3株酵母的SNP數目在22 339~26 498之間。
使用Orthofinder分析并構建K.marxianus的同源基因家族,分析得到5 110個同源基因家族,其中4 009個為12株酵母共有的同源基因家族。生成的結果中SpeciesTree_rooted.txt文件是軟件從所有同源基因家族中推斷出的STAG(Species Tree from All Genes)物種樹[25],并以STRIDE(Species Tree Root Inference from Duplication Events)為根[26]。分析結果顯示,12株K.marxianus早在約180萬年前分化為兩個亞型,第一個亞型的包含4株酵母在分化時主要發(fā)生基因擴增現象,共擴增112個基因,收縮48個基因;第二個亞型的包含8株酵母在分化時主要發(fā)生基因收縮現象,共擴增30個基因,收縮222個基因(圖4)。與未分化時相比,FIM1共擴增141個基因,收縮536個基因。
圖4 K.marxianus系統進化樹Fig.4 Phylogenetic tree of K.marxianus系統進化樹共包含12株K.marxianus的發(fā)育信息,分化時間由3株(K.lactis、P.pastoris、S.cerevisiae)已知分化時間的酵母標定校正。獲得系統發(fā)育信息后,使用iTOL將用于時間校正的3株酵母的發(fā)育信息刪除。系統進化樹樹干上,紅色代表基因擴增,藍色代表基因收縮,進化樹各個節(jié)點的support value由fasttree的SH test算法計算。Tree scale中0.1指10萬年。
根據Orthofinder生成的結果(Orthogroups_UnassignedGenes.csv),獲得FIM1的特有基因信息,結果顯示與其他K.marxianus相比,FIM1有129個特有基因,其中有功能注釋的基因有41個。將有功能注釋的41個基因進行GO富集分析,未有通路被顯著富集,但結果顯示,與其他K.marxianus相比,FIM1中NCE101、YPT31、KSH1共3個基因參與了蛋白分泌過程(GO: 0009306),該結果提示FIM1在蛋白分泌領域可能具有良好的應用潛力。
根據同源基因家族分析結果初步構建了12株酵母的泛基因組(表6)。結果顯示12株酵母共有的同源基因家族有4 009個,在12個基因組中平均占比為79.88%,附屬基因占比20.12%。附屬基因中,12株酵母多拷貝核心基因(2 733個)的平均占比為4.54%,其他非必需基因(9 380個)的平均占比為15.58%。在12株酵母中,100656-19含有最多的附屬基因,共有2 180個附屬基因,其中多拷貝核心基因741個,其他非必需基因1 439個;IIPE453的附屬基因最少,共有790個附屬基因,其中多拷貝核心基因221個,其他非必需基因569個(表6)。
表6 K.marxianus泛基因組分析Tab.6 Pan-genome analysis of K.marxianus
(續(xù)表)
基于同源基因家族構建結果,使用PanGP計算了新同源基因家族(N)與基因組數目(G)之間的關系(圖5(a)),計算出擬合方程N=1 734.8G-2.77(R2=0.816 7),當K.marxianus基因組數目增加到12時,只有2個新的同源基因家族,同源基因家族總數趨于穩(wěn)定。隨后使用PanGP軟件計算泛基因數、核心基因數和基因組數目之間的關系(圖5(b))。該軟件將泛基因定義為所有同源基因家族的集合,據此擬合得到泛基因的個數(P)與基因組數目(G)關系的為P=128 270G0-123 445(R2=1,精確到小數點后4位),根據方程預測K.marxianus泛基因組中的基因數為4 825個。核心基因的個數(C)與基因組數目(G)關系的擬合方程為C=1 113.91e-0.12G+3 760.1(R2=0.980 9)。根據方程預測K.marxianus的基因組中包含3 760個核心基因。根據新同源基因家族(N)與基因組數目(G)的擬合方程,預測K.marxianus基因組總數增加到15個以上時,不會有新的同源基因家族出現,則泛基因數不會再增加。因此本研究推測K.marxianus的泛基因組是閉合的,泛基因組由4 825個基因組成,核心基因有3 760個。
圖5 K.marxianus泛基因組分析Fig.5 Pan-genome analysis of K.marxianus(a) K.marxianus基因組數目與新基因家族數目的關系;(b) K.marxianus泛基因組分析,藍色曲線為泛基因擬合曲線,綠色為核心基因擬合曲線。
隨著測序技術的不斷進步,PacBio HiFi測序模式的面世后,3代測序已經可以同時實現重復序列的有效覆蓋與測序單堿基的高準確度[27]。這極大地提高了新物種或個體的基因組測序與組裝準確度。與此同時,我們也需要對已完成組裝的基因組序列進行校驗更新,避免它們被淘汰。本文首先使用重測序的DNA-seq數據校驗了FIM1基因組序列,隨后對校驗前后FIM1的基因組質量進行了評估。BUSCO評估顯示校驗后FIM1的基因組質量未發(fā)生改變,查詢默認參數后發(fā)現,BUSCO默認blastp的E-value值為1.0×10-3,比對閾值寬松,所以刪減基因中的短序列片段可能不會影響評估結果[9]。隨后,研究以有序列的刪除的HHT1(Histone H3)基因為例,對校驗結果進行再評估(圖2)。結果顯示校驗提高了HHT1基因的序列精度,說明DNA-seq數據校驗有助于提高FIM1基因組信息的質量。在未來的研究中,獲得新的基因組測序數據后,可以將新的數據用于校驗已完成組裝的基因組序列,從而不斷提高序列的精度。不僅是基因組序列,基因注釋也不可能一蹴而就。隨著各種注釋方法與基因功能數據庫的不斷完善,我們同樣需要持續(xù)補充與優(yōu)化注釋,這樣才能更好的將基因組信息用于實驗研究中。
隨后,本文構建了K.marxianus泛基因組,擬合方程后預測核心基因數為3 760個左右,K.marxianus泛基因組是閉合的,菌株間基因數目較穩(wěn)定,較少存在基因水平轉移(Horizontal Gene Transfer, HGT)的物種間基因交流現象[28]。HGT在原核生物中是一種重要的增加遺傳多樣性的方法[29],但在現有報道中,原核生物很少能直接通過HGT事件將基因轉移到真核生物基因組中[30]。McCarthy分析了真核生物核心基因與附屬基因的起源后認為,對于真核生物,能體現遺傳多樣性的非必需基因,主要通過基因復制或者基因收縮/擴增現象產生,而不是HGT[22]。結合本文結果,可推測K.marxianus主要通過基因組內部的片段重排、基因復制或者基因收縮擴增的方式,增加自身的遺傳多樣性,提高進化速率。
最后,NCBI現有的注釋預測K.marxianus基因的數目在4 912~5 081之間[4-6],均大于本研究預測的泛基因數4 825。該結果提示我們,現有K.marxianus注釋中可能存在較多的多拷貝核心基因以及其他非必需基因。所以下一步,我們可以從其他的非必需基因著手,通過實驗與信息學分析結合的方式,健全假基因的注釋,鑒定它們的功能?;蚪M信息的分析與研究,這始終是我們推動K.marxianus遺傳學研究深入、遺傳工程改造加快的不可或缺的一環(huán)。