基于轉(zhuǎn)錄組的池蝶蚌系統(tǒng)發(fā)育及其遺傳分析

周 偉 張萬昌 周 葉 王小敏 胡蓓娟 王軍花 洪一江

(南昌大學(xué)生命科學(xué)學(xué)院, 江西省水產(chǎn)動物資源與利用重點實驗室, 南昌 330031)

池蝶蚌(Hyriopsis schlegelii)和三角帆蚌(Hyriopsis cumingii)是目前在我國養(yǎng)殖的優(yōu)質(zhì)淡水珍珠育珠蚌, 它們具有相似的外形和生理功能, 但人們在分類學(xué)上對它們有不同的認識。有的認為它們是同屬不同種[1], 有的認為是同種不同亞種[2]。隨著分子生物學(xué)和生物化學(xué)的發(fā)展, 有研究者對池蝶蚌和三角帆蚌的同工酶和ITS1基因等進行了分析, 這些研究表明池蝶蚌和三角帆蚌的親緣關(guān)系很近, 達到了種內(nèi)分化的水平[2,3]。此外, 已有的研究表明: 相比于三角帆蚌, 池蝶蚌在生長速度、育珠性能和抗病性等方面更具優(yōu)勢[4,5]。然而, 鮮有基于基因組或轉(zhuǎn)錄組對二者進行比較分析的報道。

雖然同工酶和ITS1基因等分析表明池蝶蚌和三角帆蚌具有更近的親緣關(guān)系, 但是它們的種間關(guān)系還存在爭議, 因為這些傳統(tǒng)的分子進化研究方法幾乎都是選取一個或多個大分子的序列為標準以研究其在各物種同源序列之間的差異, 如16S rRNA和同工酶等。然而, 一個物種的基因編碼了數(shù)以萬計的序列, 以其中一個序列的差異代表整個生物體的差異是不全面的。因此, 有必要從生物體基因組或轉(zhuǎn)錄組水平來探討池蝶蚌和三角帆蚌的親緣關(guān)系。此外, 池蝶蚌生長、育珠和抗病等優(yōu)勢性能相關(guān)的分子機制也尚不清楚。然而, 池蝶蚌和三角帆蚌的基因組和轉(zhuǎn)錄組信息是未知的。

比較轉(zhuǎn)錄組學(xué)是對不同物種或近緣物種的轉(zhuǎn)錄組進行比較分析。在缺乏基因組信息時, 比較轉(zhuǎn)錄組學(xué)是分析近緣物種系統(tǒng)發(fā)育和同源基因家族的有力手段, 基于比較轉(zhuǎn)錄組學(xué)的分析主要應(yīng)用于近緣物種的進化關(guān)系和一些動植物適應(yīng)性進化相關(guān)的研究[6,7]。在本研究中, 我們對包括池蝶蚌和三角帆蚌在內(nèi)的4種淡水珠蚌轉(zhuǎn)錄組進行了二代測序(Next Generation Sequencing, NGS)和de novo組裝。在此基礎(chǔ)上, 通過比較轉(zhuǎn)錄組學(xué)的方法對池蝶蚌及其近緣物種進行了系統(tǒng)發(fā)育分析, 這旨在更深入地理解池蝶蚌和三角帆蚌的親緣關(guān)系。隨后, 我們對池蝶蚌中的特有基因進行了富集分析, 這有助于我們進一步研究池蝶蚌在生長發(fā)育、育珠和抗病性等方面的分子機制。

1 材料與方法

1.1 樣品收集、總RNA提取和轉(zhuǎn)錄組測序

池蝶蚌取自江西省撫州市池蝶蚌良種場, 三角帆蚌取自浙江省金華三角帆蚌良種場, 褶紋冠蚌(Cristaria plicata)和背角無齒蚌(Anodonta woodiana)分別采集于江西鄱陽湖和山東省微山湖。對這4種蚌每種隨機選取健康的三齡蚌10只, 置于充分曝氣的水族箱內(nèi)暫養(yǎng)一周, 水溫控制在18—25℃。

從水族箱中分別取出4種淡水珠蚌, 每種蚌為一組, 每組各3只。分別取每組每只蚌的多種組織(包括閉殼肌、心臟、腎臟、鰓、斧足、性腺、肝臟、腸胃和外套膜)各100 mg, 并分別將取出的三份相同組織混合。隨后采用Trizol試劑法提取4種淡水珠蚌的總RNA。以每種蚌各組織提取到的總RNA濃度為依據(jù),取等量的各組織總RNA進行混合, 最后分別得到4種淡水珠蚌混合組織的總RNA樣品。本研究采用的建庫類型為真核鏈特異性文庫(250—300 bp), 之后在Illumina HiseqXten平臺上進行雙端測序。

1.2 De novo組裝和功能注釋

用Trimmomatic[8]移除reads中的Illumina測序接頭和引物序列, 去除質(zhì)量過低的序列(Phred數(shù)值小于20), 過濾過短(序列長度小于50 bp)和不合格的序列。采用Trinity軟件[9]分別對4種淡水珠蚌轉(zhuǎn)錄組測序數(shù)據(jù)進行de novo組裝, 所有參數(shù)設(shè)置為默認參數(shù)。首先用CD-HIT軟件去除冗余轉(zhuǎn)錄本序列[10],再用Corset軟件進行分層聚類[11], 將分層聚類后得到的轉(zhuǎn)錄本序列進行后續(xù)的分析。對轉(zhuǎn)錄本序列進行統(tǒng)計并對其拼接質(zhì)量進行評估。

首先用TransDecoder軟件(Find Coding Region Within Transcripts, https://transdecoder.github.io/)對轉(zhuǎn)錄本序列進行ORF框預(yù)測并翻譯成蛋白序列, 然后用EggNOG-mapper將4種淡水珠蚌的蛋白序列分別比對到EggNOG數(shù)據(jù)庫, 依據(jù)同源序列的相似性進行基因功能注釋[12—14]。

1.3 直系同源基因鑒定和基因家族聚類

用Orthofinder對4種淡水育珠蚌的蛋白序列進行同源基因家族分析[15], 參數(shù)設(shè)置使用Diamond[14]以默認的Evalue值1e-3進行all-vs-all搜索, 尋找潛在的同源基因。使用馬爾可夫聚類算法(Markov Cluster Algorithm, MCL)對基因進行聚類, 劃分直系同源簇。提取4種淡水珠蚌中各自特有的基因家族, 結(jié)合之前得到的注釋信息獲取4種淡水珠蚌特有基因?qū)?yīng)的功能?；蚣易寰垲惖慕Y(jié)果可以直接從orthofinder分析結(jié)果中獲取。

1.4 系統(tǒng)發(fā)育分析

對4種蚌的單拷貝同源基因進行多序列聯(lián)配(Multiple sequence alignment, MSA), 在MEGA軟件中采用pairwise uncorrected p-distance模型對多序列聯(lián)配得到的長片段序列的遺傳距離進行計算[16]。從Uniprot數(shù)據(jù)庫(https://www.uniprot.org/proteomes/)下載光滑雙臍螺(Biomphalaria glabrata)、太平洋牡蠣(Crassostrea gigas)、霸王蓮花青螺(Lottia gigantea)、蝦夷扇貝(Mizuhopecten yessoensis)和加利福尼亞雙斑章魚(Octopus bimaculoides)這5個物種的蛋白質(zhì)組數(shù)據(jù), 結(jié)合本研究中的4種淡水珠蚌蛋白序列數(shù)據(jù)共9個貝類物種以供進一步的分析。隨后利用Orthofinder對這9個物種蛋白質(zhì)組序列篩選出的1317個共有單拷貝同源基因重建了貝類的系統(tǒng)發(fā)育。默認使用mafft[17,18]進行多序列聯(lián)配, 用fasttree推斷系統(tǒng)發(fā)生樹[19]。

通過生命進化時間的公共信息庫TimeTree網(wǎng)站(http://www.timetree.org/)查找到霸王蓮花青螺與光滑雙臍螺分歧時間大約在471(423—519)百萬年前(Millions of years ago, Mya), 太平洋牡蠣與蝦夷扇貝分歧時間大約在455(421—497)百萬年前。以上述2個已知的分歧時間作為固定值和校正值, 根據(jù)上一步得到的系統(tǒng)發(fā)生樹用最大似然法(r8s軟件)預(yù)測它們的分歧時間[20]。

1.5 池蝶蚌特有基因分析

首先用AnnotationForge[21]依據(jù)池蝶蚌功能注釋信息構(gòu)建OrgDb數(shù)據(jù)包, OrgDb包用于存儲不同數(shù)據(jù)庫基因ID之間的對應(yīng)關(guān)系及基因與GO、KEGG等注釋的相對關(guān)系。根據(jù)同源基因家族的分析結(jié)果提取池蝶蚌中特有基因, 隨后基于本研究構(gòu)建的池蝶蚌OrgDb數(shù)據(jù)包用clusterprofiler[22]對池蝶蚌中特有基因進行KEGG和GO富集分析。

2 結(jié)果

2.1 轉(zhuǎn)錄組組裝和注釋

通過Illumina測序獲得了池蝶蚌、三角帆蚌、背角無齒蚌和褶紋冠蚌的轉(zhuǎn)錄組數(shù)據(jù), 質(zhì)控后分別得到12.97、13.58、13.41和11.36 Gb轉(zhuǎn)錄組測序數(shù)據(jù), 各含有clean reads約4323萬、4525萬、4471萬和3785萬條。對這4種蚌的轉(zhuǎn)錄組數(shù)據(jù)進行了de novo組裝, 其中池蝶蚌轉(zhuǎn)錄組組裝后得到138948個轉(zhuǎn)錄本和104803個Unigenes, reads回貼率(reads回貼到轉(zhuǎn)錄本的比率)達85.62%, Unigenes的平均長度為1373.20 bp, N50長度為1928 bp(表 1)。以后生動物單拷貝同源基因集作為比對數(shù)據(jù)庫, 使用Busco軟件對轉(zhuǎn)錄組完整性進行分析, 結(jié)果表明: 在4種淡水珠蚌的轉(zhuǎn)錄本中, 超過98%都是完整的, 只有不到2%的轉(zhuǎn)錄本是部分比對或比對不上, 這說明我們組裝結(jié)果完整度較好, 轉(zhuǎn)錄本序列可用于后續(xù)分析(表 2)。

將得到的池蝶蚌、三角帆蚌、背角無齒蚌和褶紋冠蚌基因翻譯成蛋白序列后比對到EggNOG數(shù)據(jù)庫, 基于相似性(E-value＜10–3), 這4種蚌分別有34143、28180、28609和24778個基因得到了注釋,在蛋白質(zhì)直系同源簇(Clusters of orthologous groups of proteins, COG)分類統(tǒng)計表明, 這4種育珠蚌的COG類別分布大致相同, 其中, 富集基因最多的2個分類是“功能未知(S. Function unknown)”和“信號轉(zhuǎn)導(dǎo)機制(T. Signal transduction mechanisms)”(圖 1)。

表 1 轉(zhuǎn)錄組組裝數(shù)據(jù)統(tǒng)計Tab. 1 The statistic of transcriptome assembled data

表 2 組裝轉(zhuǎn)錄組的BUSCO分析結(jié)果Tab. 2 Results of BUSCO analysis on the assembled transcriptomes

圖 1 池蝶蚌(a)、三角帆蚌(b)、背角無齒蚌(c)和褶紋冠蚌(d)的蛋白質(zhì)直系同源簇(COG)分類Fig. 1 Clusters of Orthologous Group of proteins (COG) classifications for H. schlegelii (a), H. cumingii (b), A. woodiana (c) and C.plicata (d)

2.2 直系同源基因的鑒定及基因家族聚類

通過Orthofinder在4種淡水珠蚌中共鑒定到29468個潛在的同源基因家族, 其中4個物種共享的16001個同源基因家族構(gòu)成了淡水珠蚌核心基因集,這些核心基因參與廣泛的基礎(chǔ)活動, 幾乎在每一個GO功能中都有分布。只存在于池蝶蚌和三角帆蚌中的同源基因家族有5175個(圖 2)。另外, 池蝶蚌、三角帆蚌、褶紋冠蚌和背角無齒蚌4種育珠蚌間各自特有的基因家族數(shù)量分別為12、11、9和15個。注釋結(jié)果表明, 池蝶蚌特有基因主要與生長發(fā)育(TRIM45、TRIM33和JAG1等)和免疫(ABCB9)等過程相關(guān), 三角帆蚌的特有基因與抗腫瘤(CUZD1)和甾醇代謝(HDLBP)等功能有關(guān), 褶紋冠蚌特有基因主要與離子平衡(SLC12A2)和免疫(CEACAM1)功能相關(guān), 背角無齒蚌特有基因主要與疾病(TMPRSS12和CLCA2)、細胞內(nèi)通訊和信號傳遞(FLNB)等功能相關(guān)(表 3)。

圖 2 池蝶蚌、三角帆蚌、褶紋冠蚌和背角無齒蚌中的基因家族聚類韋恩圖Fig. 2 Gene families clustering venn diagram in H. schlegelii, H.cumingii, C. plicata and A. woodiana

表 3 四種淡水珠蚌中物種特有基因列表Tab. 3 The list of species-specific genes in four freshwater mussels

2.3 系統(tǒng)發(fā)育分析

我們對池蝶蚌、三角帆蚌、背角無齒蚌和褶紋冠蚌的單拷貝同源基因進行了多序列聯(lián)配, 采用pairwise uncorrected p-distance模型對多序列聯(lián)配得到的長片段序列計算遺傳距離(E-value＜1e–3), 結(jié)果表明池蝶蚌和三角帆蚌的遺傳距離小于1%, 而池蝶蚌與其他物種的遺傳距離幾乎是它與三角帆蚌的10倍(表 4)。結(jié)合Uniprot數(shù)據(jù)庫中已發(fā)布的光滑雙臍螺、太平洋牡蠣和蝦夷扇貝等5個物種的蛋白質(zhì)組數(shù)據(jù)和本研究中獲得的4種淡水珠蚌的Unigenes對應(yīng)的蛋白質(zhì)組數(shù)據(jù), 我們篩選了包含1317個共有單拷貝直系同源基因家族的蛋白集, 隨后利用這些單拷貝同源蛋白集重建了貝類的系統(tǒng)發(fā)育。結(jié)果表明, 池蝶蚌與三角帆蚌的親緣關(guān)系最近, 2個物種聚為一支(帆蚌屬)后再與冠蚌屬、無齒蚌屬聚在一起, 瓣鰓綱的真瓣鰓目與異柱目聚為一支(瓣鰓綱), 然后它們再與腹足綱、頭足綱聚在一起(圖 3)。隨后, 我們對這9個物種進行了分歧時間評估, 結(jié)果表明, 池蝶蚌與三角帆蚌大約在424萬年前發(fā)生分歧, 帆蚌屬與冠蚌屬、無齒蚌屬發(fā)生分歧的時間大約在5144萬年前, 真瓣鰓目與異柱目大約在5.72億年前發(fā)生分歧, 瓣鰓綱、腹足綱和頭足綱大約在6.18億年前發(fā)生分歧(圖 4)。此外, 我們還對它們的同源基因家族擴張和收縮情況進行了分析, 發(fā)現(xiàn)池蝶蚌在淡水珠蚌中發(fā)生了最多的同源基因家族擴張。

表 4 基于4種淡水珠蚌單拷貝同源基因序列的遺傳分化距離Tab. 4 Genetic differentiation distance based on single copy ortholog sequences of four freshwater mussel

圖 3 九種貝類動物的系統(tǒng)發(fā)育樹Fig. 3 Phylogenetic tree of nine species in mollusca

圖 4 九種貝類動物的分歧時間及其基因家族擴張/收縮情況Fig. 4 Time of divergence of nine species in mollusca and their expansion/contraction of genefamilies

2.4 池蝶蚌相對于三角帆蚌的特有基因

通過比較池蝶蚌和三角帆蚌的轉(zhuǎn)錄組, 得到了池蝶蚌轉(zhuǎn)錄組中相對于三角帆蚌特異性表達的基因3721條, 其中345個特有基因得到了注釋。這些特有基因可能與池蝶蚌的生長發(fā)育、育珠及抗病能力相關(guān)。我們對這些特有基因進行了KEGG和GO富集分析。KEGG 通路富集分析結(jié)果表明, 池蝶蚌特有基因主要富集在免疫細胞膜分子(CD molecules)和蛋白消化吸收(Protein digestion and absorption)等通路上(圖 5)。此外, 池蝶蚌特有基因GO注釋結(jié)果表明: 有很多特有基因與池蝶蚌的優(yōu)良性狀相關(guān), GO類別包括生長發(fā)育過程(91個), 系統(tǒng)發(fā)育(69個), 動物器官發(fā)育(51個), 免疫系統(tǒng)過程(25個), 抗原的處理和呈遞(4個, 圖 6)。

圖 5 池蝶蚌特有基因的KEGG富集分析Fig. 5 KEGG enrichment analysis of species-specific genes in H. schlegelii

圖 6 池蝶蚌優(yōu)勢性能相關(guān)的特有基因GO富集Fig. 6 GO enrichment of specific genes related to ecological adaptation in H. schlegelii

3 討論

3.1 池蝶蚌與三角帆蚌的親緣關(guān)系

池蝶蚌和三角帆蚌的親緣關(guān)系一直存在爭議,這是由于我國的淡水雙殼類蚌科分類系統(tǒng)主要依據(jù)傳統(tǒng)的形態(tài)學(xué)分類, 而蚌科動物存在趨同現(xiàn)象且貝殼具有高度可塑性, 因此傳統(tǒng)的形態(tài)學(xué)分類并不能完全反映各類群之間的親緣關(guān)系。近年來, 分子生物學(xué)技術(shù)常用于解決一些分類學(xué)上的疑難問題。在淡水珠蚌中的分子進化研究方法, 大多采用某個大分子的序列為標準, 通過比較相應(yīng)物種的同源序列分析物種的進化關(guān)系。以往池蝶蚌和三角帆蚌的分子進化研究主要是基于同工酶[2]、ITS1[3]和線粒體DNA[23]的比較分析, 但這樣單個基因或少數(shù)基因的比較不能代表整體差異, 而且可能由于遺傳變異存在一定的偶然性, 不足以使人信服。比較基因組或比較轉(zhuǎn)錄組可以最大限度地上對不同物種的遺傳差異進行分析, 在缺乏基因組信息的非模式生物中, 往往可以基于轉(zhuǎn)錄組數(shù)據(jù)解決物種的進化關(guān)系[24,25]。我們基于轉(zhuǎn)錄組水平的系統(tǒng)發(fā)育分析表明池蝶蚌與三角帆蚌的遺傳距離不到1%, 大大低于其他淡水珠蚌間的遺傳距離, 這與傳統(tǒng)的形態(tài)學(xué)分類有所不同, 但與Sakai等[2]的研究結(jié)果一致,很大程度上支持了池蝶蚌與三角帆蚌可能是不同亞種或不同種群的假設(shè), 二者的形態(tài)學(xué)差異可能是由于生存環(huán)境差異造成的。進一步的分析表明, 小方蚌亞科與無齒蚌亞科大約在5144萬年前發(fā)生分歧, 而池蝶蚌和三角帆蚌發(fā)生分化的時間大約在424萬年前, 這表明可能約在424萬年前由于某些環(huán)境變化使池蝶蚌和三角帆蚌發(fā)生地理隔離導(dǎo)致它們產(chǎn)生一些遺傳上的差異, 但這些差異還處于種內(nèi)分化的水平。

3.2 池蝶蚌中的特有基因

特有基因的分析能幫助我們理解表型性狀產(chǎn)生的分子機制, 同時可以揭示基因表達模式與表型形成之間的關(guān)系。有研究表明, 池蝶蚌生長速度快,其殼寬、外套膜的厚度、晶桿的長度和直徑都大于三角帆蚌[4]。個體的生長及其形態(tài)與生長發(fā)育過程息息相關(guān), 如外套膜的生長與發(fā)育。本研究對池蝶蚌和三角帆蚌轉(zhuǎn)錄組中的基因進行比較分析, 結(jié)果表明, 根據(jù)GO分類, 多達51個池蝶蚌特有基因被標注為動物器官發(fā)育。其中, 基因α-Ⅰ型膠原質(zhì)(COL1A1)主要負責(zé)編碼Ⅰ型膠原的主要成分, 即大部分結(jié)締組織中的纖維膠原, 包括軟骨[26]。另有基因JAG1和Notch信號的調(diào)控有關(guān), 該基因在胚胎發(fā)育中起著重要作用, 可能涉及造血過程中的細胞命運決定和抑制成肌細胞的分化[27]。結(jié)締組織中的纖維膠原可能與池蝶蚌外套膜的厚度相關(guān), 而胚胎發(fā)育對有機體的生長發(fā)育至關(guān)重要。這表明池蝶蚌中可能存在一些特有基因或者調(diào)控這些特有基因表達的相關(guān)基因, 使得池蝶蚌在生長發(fā)育方面優(yōu)于三角帆蚌。

此外, 有研究表明池蝶蚌對蚌瘟病和嗜水氣單胞菌的抵抗力明顯強于三角帆蚌[5,28]。有機體的抗病能力主要依賴免疫系統(tǒng)發(fā)揮防御作用, 這個過程包括抗原的識別和免疫系統(tǒng)的調(diào)控。GO富集結(jié)果表明有25個池蝶蚌特有基因涉及免疫系統(tǒng)過程, 這個結(jié)果與富集在免疫細胞膜分子通路中的13個基因相一致。作為免疫系統(tǒng)過程的重要環(huán)節(jié), 抗原的處理和呈遞發(fā)揮著重要的作用。例如, 基因ABCB9是ATP依賴的低親和力肽轉(zhuǎn)運體, 它將廣譜的多肽從胞漿轉(zhuǎn)移到溶酶體腔[29]。類似地, 基因SEC24A編碼衣殼蛋白復(fù)合物Ⅱ的組成成分, 促進內(nèi)質(zhì)網(wǎng)轉(zhuǎn)運囊泡的形成。衣殼主要有兩種功能, 一是內(nèi)質(zhì)網(wǎng)膜在囊泡中的物理變形, 二是選擇轉(zhuǎn)運到高爾基體的貨物分子, 這與異質(zhì)蛋白在內(nèi)質(zhì)網(wǎng)-高爾基體界面被識別、分揀和捕獲有關(guān)。抗原的有效識別和處理決定了免疫系統(tǒng)能否有效地發(fā)揮作用[30]?；谝陨辖Y(jié)果, 我們推測池蝶蚌優(yōu)良的抗病能力可能與池蝶蚌中特有的免疫基因有關(guān)。

4 結(jié)論

本研究通過轉(zhuǎn)錄組學(xué)分析對池蝶蚌和三角帆蚌進行了比較, 為厘清池蝶蚌和三角帆蚌的親緣關(guān)系提供了重要的證據(jù)和幫助, 并為研究池蝶蚌優(yōu)良性狀相關(guān)的分子機制提供了新的視角和有力的支持。在本研究中, 我們的計算表明池蝶蚌和三角帆蚌的遺傳距離不到1%, 大約只有其他蚌種間遺傳距離的十分之一, 這與種內(nèi)分化水平相當。另外,本研究對包括池蝶蚌、三角帆蚌、背角無齒蚌和褶紋冠蚌在內(nèi)的9種貝類重建了貝類系統(tǒng)發(fā)育分析,結(jié)果表明, 池蝶蚌和三角帆蚌發(fā)生分歧的時間大約在424萬年前。池蝶蚌產(chǎn)自日本, 而三角帆蚌是中國本土物種, 可以推斷, 它們可能在那個時間點由于某些事件發(fā)生了地理上的隔離。池蝶蚌在一些性能上優(yōu)于三角帆蚌, 這很可能是遺傳差異導(dǎo)致的。池蝶蚌轉(zhuǎn)錄組中特有基因富集分析結(jié)果說明池蝶蚌具有潛在的優(yōu)勢性能, 如生長發(fā)育和免疫。這很可能是池蝶蚌在生長速度、外套膜厚度和抗病能力等方面具有優(yōu)勢的關(guān)鍵所在。由于轉(zhuǎn)錄組學(xué)分析難以獲得全部和完整的遺傳信息, 因此需要進一步的研究。盡管如此, 本研究依然為池蝶蚌和三角帆蚌的親緣關(guān)系提供了重要參考價值, 為進一步研究池蝶蚌生長發(fā)育和免疫等相關(guān)的分子機制提供了有益線索。