“三臺”核桃不同成熟期果實轉錄組分析

賀娜　耿樹香　寧德魯　肖良俊

摘要 以云南優(yōu)良品種“三臺”核桃4個不同成熟期種仁為材料，從分子水平上探索云南核桃（Juglans sigillata）種子不同成熟期種仁合成相關基因的表達模式。共有53 868條unigenes得到了注釋，其中33 407條unigenes注釋在Nr數(shù)據(jù)庫中;31 725條unigenes注釋在KEGG數(shù)據(jù)庫中;Swissprot數(shù)據(jù)庫中注釋到19 842條unigenes;COG數(shù)據(jù)庫中注釋到16 655條unigenes。四大數(shù)據(jù)庫共注釋到33 753條unigenes，占總unigenes的62.66%。CK-vs-HT1、CK-vs-HT2、CK-vs-HT3、HT1-vs-HT2、HT1-vs-HT3、HT2-vs-HT3分別有 2 417、2 703、4 166、1 256、5 347、3 283條上調表達的差異基因和9 515、12 455、9 498、1 751、3 132、1 365條下調表達的差異基因。不同時期基因差異表達富集度最高的為脂肪酸生物合成代謝通路。核桃果實在成熟中后期，油脂及蛋白大量積累，相關基因表達量增加。

關鍵詞 “三臺”核桃 ;成熟度;轉錄組;基因

中圖分類號 S 664.1 文獻標識碼 A 文章編號 0517-6611（2020）23-0165-05

doi：10.3969/j.issn.0517-6611.2020.23.041

Transcriptome Analysis of Santai Walnut Fruits at Different Ripening Stages

HE Na， GENG Shu-xiang， NING De-lu et al

（Yunnan Academy of Forestry and Grassland Sciences，Kunming， Yunnan 650201）

Abstract The expression patterns of genes related to seed synthesis at different mature stages of Juglans sigillata were explored at the molecular level by using 4 different mature kernels of Santai walnut as materials. A total of 53 868 unigenes were annotated， of which 33 407 were in the NR database， 31 725 in the KEGG database， 19 842 in the Swissprot database， and 16 655 in the COG database. A total of 33 753 unigenes were annotated in the four databases， accounting for 62.66% of the total. There were 2 417， 2 703， 4 166， 1 256， 5 347， 3 283 up regulated genes and 9 515， 1 2455， 9 498， 1 751， 3 132， 1 365 down regulated genes in CK-vs-HT1， CK-vs-HT2， CK-vs-HT3， HT1-vs-HT2， HT1-vs-HT3 and HT2-vs-HT3， respectively. Fatty acid biosynthesis pathway was the highest gene expression enrichment in different stages. In the middle and late stage of walnut fruit ripening， a large amount of oil and protein were accumulated， and the expression of related genes was increased.

Key words Santai walnut;Maturity;Transcriptome;Gene

云南“三臺”核桃是云南眾多深紋核桃（Juglans sigillata）中的一個主要栽培品種，產于云南賓川縣、大姚縣三臺鄉(xiāng)等地，無性系品種，脂肪含量44.42%～72.74%，蛋白質含量17.26%[1]。普通核桃（Juglans regia）為二倍體植物（2n=32），全基因組約606 Mb[2]。核桃是藥食兼用植物，隨高通量測序技術的飛速發(fā)展，基于高通量測序的轉錄組分析成為核桃代謝組研究的有力工具。轉錄組學是核桃功能基因組學一個重要分支，在發(fā)現(xiàn)核桃次生代謝產物生物合成關鍵基因、闡明次生代謝調控、篩選分子標記等方面具有重要的應用價值[3]。核桃代謝產物的研究主要集中于普通核桃脂肪代謝機制的分子生物學研究上，利用轉錄組測序技術構建普通核桃花后60、90和120 d脂肪積累期胚組織的轉錄組、表達譜，同時分析核桃胚脂肪積累期的基因表達模式及其代謝通路及核桃胚脂肪積累期脂肪代謝相關基因的表達模式，并利用轉錄組數(shù)據(jù)，克隆核桃異質型乙酰輔酶A羧化酶People acetyl coenzyme A carboxylase（ACCase）4個亞基的基因ORF[4-5]。

項目在前期的研究中發(fā)現(xiàn)云南“三臺”核桃4個不同成熟期中的代謝產物含量差異較大，特別是次生代謝產物角鯊烯及磷脂類。選擇4個不同成熟期“三臺”核桃種仁構建云南核桃不同成熟期轉錄組文庫，并利用生物信息學方法對轉錄的基因結構進行分析，對不同發(fā)育階段差異表達基因進行GO功能富集分析，初步闡明不同成熟期云南核桃不同發(fā)育階段差異表達基因功能及轉錄水平。

1 材料與方法

1.1 植物材料

以云南主栽核桃 “三臺”為對象，依據(jù)方文亮等[1]對云南核桃果實成熟動態(tài)的研究結果和Li等[6]對果實發(fā)育動態(tài)的研究結果，分別在所選3株“三臺”核桃樣樹樹冠中部東南西北4個方位隨機采摘30個果實，立取種仁置于液氮中，于-80 ℃冰箱保存，備用。4個成熟期分別為硬核期（6月CK樣）、油脂轉化期（7月HT1樣）、油脂積累期（8月HT2樣）、果實成熟期（9月HT3樣）。

1.2 主要設備

冷凍高速離心機、核酸定量儀、-80 ℃超低溫冰箱（德國Thermo公司）、全自動凝膠圖像分析系統(tǒng)、水平電泳儀（美國 Bio-rad公司）、PCR擴增儀（美國ABI公司）、NanoPhotometer spectrophotometer、QuantStudio 12K Flex Real-Time PCR System（美國Thermo Fisher公司）、全套Eppendorf微量移液器（德國Eppendorf公司）、Agilent 2100 Bioanalyzer（Agilent RNA 6000 Nano Kit）。

1.3 主要試劑

Trizol reagent（Promega，Madison，USA）、NanoPhotometer spectrophotometer（IMPLEN，CA，USA）、NEBNext Ultra RNA Library Prep Kit for Illumina（Illumina Inc，San Diego，USA）、TruSeq PE Cluster Kit v3-cBot-HS（Illumia，San Diego，USA）、 SYBR Premix Ex Taq（TaKaRa，RR420）、DNA Marker-DL2000（TaKaRa 公司）、PCR 2×mix（上海碧云天生物技術有限公司）、SYBR Green Realtime PCR Master Mix（日本 Toyobo）。

1.4 試驗及數(shù)據(jù)處理方法

1.4.1 “三臺”核桃果實RNA的提取。

使用Trizol試劑提取“三臺”核桃果實總RNA20 μg，再用1%瓊脂糖凝膠電泳檢測RNA是否污染及降解，用Agilent2100 Bioanalyzer檢測總 RNA的濃度、RIN值等。使用Nano Drop超微量紫外分光光度計進行總RNA純度檢測。

1.4.2 cDNA文庫構建和測序。

“三臺”核桃的轉錄組測序由廣州基迪奧生物科技有限公司完成。構建流程如下：通過帶有Oligo（dT）的磁珠富集具有polyA尾巴的真核mRNA后，用超聲波把mRNA打斷。以片段化的mRNA為模板，隨機寡核苷酸為引物，在M-MuLV逆轉錄酶體系中合成cDNA第一條鏈，隨后用RNase H降解RNA鏈，并在DNA polymerase I 體系下，以dNTPs為原料合成cDNA第二條鏈。純化后的雙鏈cDNA經過末端修復、加A尾并連接測序接頭，用AMPure XP beads篩選200 bp左右的cDNA，進行PCR擴增并再次使用AMPure XP beads純化PCR產物，最終獲得文庫[7-8]。

1.4.3 測序數(shù)據(jù)的組裝。

經測序儀產生的原始數(shù)據(jù)經base caling轉化為序列數(shù)據(jù)，稱為raw reads，結果以fastq文件存儲，某些raw reads 帶有adaptor序列，含有少量低質量序列包括N比例大于10%的reads、質量值Q≤5的堿基數(shù)占整個read的50%以上，經數(shù)據(jù)處理，得clean reads。使用短reads組裝軟件Trinity重頭組裝，先將具有一定長度overlap的reads連成更長片斷，通過reads overlap 關系得到的組裝片斷為Contig，再將reads比對回Contig，通過paired-end reads 確定來自同一轉錄本不同Contig及Contig之間的距離，Trinity將這些Contig 連一起，得到兩端不能再延長的序列稱為Unigene，使用Tgicl去冗余和拼接，再對這些序列進行同源轉錄本聚類，得最終Unigene[9]，該Unigene分為兩部分，一是clusters（CL開頭），另一是singletons（以Unigene開頭），再將Unigene序列與Blastx比對，取最好比對結果蛋白確定Unigene序列方向。

1.4.4 轉錄組基因功能注釋和表達分析。

為了獲得基因功能注釋，將所得的unigenes與NCBI的Nt數(shù)據(jù)庫（NCBI non-redundant protein sequences database）、Nr數(shù)據(jù)庫（NCBI nucleotide sequences database）、KOG 數(shù)據(jù)庫（euKaryotic Ortholog Groups）以及Swiss-Prot數(shù)據(jù)庫等進行比對，過濾 E-value 的閾值設置為較低 1×10-5。在獲得所有 unigenes 的GO注釋信息后，使用WEGO軟件進行GO功能分類。此外，使用Blastx，過濾 E-value 的閾值設置為1×10-10，將所有unigenes與KEGG數(shù)據(jù)（the Kyoto Encyclopedia of Genes and Genomes Pathway Database）進行比對，分析核桃相關的代謝途徑，挖掘與核桃次生代謝產物合成相關的基因。

1.4.5 基因差異表達分析。

使用 Bowtie 2[10]將 clean reads 比對到參考序列以統(tǒng)計基因比對率，之后再使用 RSEM[11]計算基因和轉錄本的表達水平。RSEM是用于轉錄組reads計算基因以及轉錄本表達量的軟件包。RSEM 用 Paired-end的關系、reads的長度、fragment的長度分布、質量值等，基于最大期望的算法建立最大似然的豐度估計模型，用以區(qū)分哪些轉錄本是同一個基因的不同亞型。表達定量的結果以FPKM 為單位。差異表達基因（differentially expressed gene，DEG）檢測基于泊松分布。為了提高 DEG的準確性，將差異倍數(shù)為2倍以上并且 Q-value≤0.001 的基因定義為顯著差異表達基因。

2 結果與分析

2.1 轉錄組質量評估

使用 BGISEQ-500平臺進行測序，測序結果如表1所示，共得到 7 306 Mb 的原始數(shù)據(jù)，48 707 124條raw reads。通過去除低質量 reads，獲得48 647 921條clean reads，且clean reads的Q20比例為97.71%，GC含量為42.91%。以上數(shù)據(jù)證明本次測序質量較好，可用于后續(xù)拼接與分析。利用轉錄組拼接軟件 Trinity對clean reads進行從頭拼接，拼接獲得的 unigene 數(shù)目為53 868個，總長度為 52 736 772 bp，平均長度為979 bp，N50長度為1 710 bp。組裝結果質量評估可從N50來評估。將所有unigene從長到短排序，并依次累加長度，當累加片段長度達到總片段長度（所有unigene的長度）的50%時，對應那個片段的長度和數(shù)量，即為unigene N50長度和數(shù)量。unigene N50越長，數(shù)量越少，說明組裝質量越好[12]。此項分析統(tǒng)計結果中包含了N50、GC含量等詳細統(tǒng)計結果（表1）。

2.2 基因注釋與分類

利用本地版Blast將獲得的unigenes與四大功能數(shù)據(jù)庫 Nr、KEGG、Swissprot、GO/COG進行比對，注釋情況見圖1 。共有53 868條unigenes得到了注釋，其中33 407條unigenes注釋在Nr數(shù)據(jù)庫中;31 725條unigenes注釋在KEGG數(shù)據(jù)庫中;Swissprot數(shù)據(jù)庫中注釋到19 842條unigenes;COG數(shù)據(jù)庫中注釋到16 655 unigenes。四大數(shù)據(jù)庫共注釋到33 753條unigenes，占總unigenes的62.66%。此外，仍有20 115條unigenes未獲得注釋信息，可能是由于數(shù)據(jù)庫中缺乏相關注釋信息或者為新的轉錄本。

2.3 Nr注釋分析

Nr數(shù)據(jù)庫屬于非冗余蛋白序列數(shù)據(jù)庫NR（non-redundant protein sequence database）非冗余蛋白庫，NCBI收錄，所有GenBank+EMBL+DDBJ+PDB中的非冗余蛋白序列，對于所有已知的或可能的編碼序列，NR記錄中都給出了相應的氨基酸序列（通過已知或可能的讀碼框推斷而來）以及專門蛋白數(shù)據(jù)庫中的序列號[13]?！叭_”核桃序列與Nr數(shù)據(jù)庫比對的結果如圖2。其中有25 505條unigenes比對到云南核桃即鐵核桃上，基因配比率達76.34%，其余依次為土瓶草（2.46%）、可可（2.12%）、葡萄（1.14%）、木豆（1.01%）、黃麻（0.85%）、粉掌（0.76%）、棗（0.74%）、烏梅（0.61%）等。而其余13.37%的序列沒有明確匹配到相應物種，可能是因為分散到了多種其他物種，或由于產生了一些數(shù)據(jù)庫里沒有的新基因。

2.4 KEGG pathway分析

基于 KEGG pathway的分析有助于進一步了解核桃基因的生物學功能。“三臺”核桃的53 868條unigenes中有31 725條注釋到了KEGG數(shù)據(jù)庫中（表2）。KEGG（kyoto encyclopedia of genes and genomes）的代謝途徑主要分為新陳代謝（metabolism）、遺傳信息過程（geneticInformation processing）、環(huán)境信息過程（environmental information processing）、細胞過程（cellular processes）和生物系統(tǒng)（organismal systems）等。其中，新陳代謝共注釋到 6 662條unigenes（20.99%），遺傳信息過程注釋到 2 298條 unigenes（7.87%），環(huán)境信息過程注釋到546條unigenes，細胞過程注釋到567條 unigenes，生物系統(tǒng)注釋到422條unigenes（圖3）。

核桃的主要生物活性物質為脂類和蛋白，因此著重對新陳代謝（metabolism）分類下的二級分類進行分析。在KEGG注釋中與新陳代謝相關的是global and overview maps（2 633條unigenes，39.52%），碳水化合代謝（936條unigenes，14.05%），能量代謝（589條unigenes，8.84%），氨基酸代謝（586條unigenes，8.79%），脂類代謝（510條unigenes，7.65%），

輔因子和維生素代謝（310條unigenes，4.65%），其他次級代謝生物合成（256條unigenes，3.84%），核苷酸代謝（254條unigenes，3.81%），其他氨基酸代謝（253條unigenes，3.79%），萜類及聚酮化合物代謝（189條unigenes，2.83%），多糖的生物合成和代謝（146條unigenes，2.19%）。這些注釋將為進一步研究與核桃油脂及蛋白合成有關的代謝途徑、基因的結構和功能提供參考信息。

2.5 COG功能分類

COG 是直系同源家族蛋白數(shù)據(jù)庫，按照COG的分類信息將“三臺”核桃的16 655條unigenes分配到 25個類別里，結果見圖4。其中，最多的是一般功能預測，包含4 870條unigenes，占總unigenes的29.24%;其余有代表性的包括信號轉導機制（2 260條 unigenes，13.56%），轉錄后修飾、蛋白轉換、伴侶蛋白（1 889條 unigenes，11.34%），翻譯（1 121條 unigenes，6.73%），翻譯、核糖體構成和起源（960條 unigenes，5.76%），胞內運輸、分泌及囊泡轉運（917條 unigenes，5.50%），未知功能（903條 unigenes，5.42%），碳水化合物的運輸和代謝（892條 unigenes，5.35%），RNA加工及修飾（784條 unigenes，4.70%），能量產生及轉化（704條 unigenes，4.22%），次生代謝產物合成、運輸及分解（685條 unigenes，4.11%）等，其中最少是細胞遷移（13條 unigenes，0.08%），其中未知功能（903條 unigenes，5.42%）表明“三臺”核桃中存在許多新基因，其生物學功能有待挖掘。

2.6 不同成熟期核桃差異基因分析

將差異倍數(shù)為2倍以上并且 Q-value≤0.001 的基因確定為顯著差異表達基因。對云南“三臺”核桃4個不同成熟期樣進行兩兩比較，進行差異基因分析（圖5）。其中7月與6月相比有2 417個基因上調表達，9 515個基因下調表達;8月與6月相比有2 703個基因上調表達，12 455個基因下調表達;9月與6月相比有4 166個基因上調表達，9 498個基因下調表達;8月與7月相比有1 256個基因上調表達，1 751個基因下調表達;9月與7月相比有5 347個基因上調表達，3 132個基因下調表達;9月與8月相比有3 283個基因上調表達，1 365個基因下調表達。核桃果實在成熟中后期，油脂及蛋白大量積累，相關基因表達量增加。在核桃發(fā)育期，果實油脂積累期與硬核期相比，下調表達基因最多，也是差異基因較多的一組，這與CK樣品有關，CK樣品為果實整個子房，其他3個樣品（HT1、HT2、HT3）均為果肉，另外，果實形成初期也是果實整個發(fā)育過程中代謝最活躍的時期，尤其是激素在子房中的含量，子房中幾乎包括果實發(fā)育期間所有必需的生物代謝途徑[14]。從HT1到HT3，后一個階段都比前一階段上調表達的基因多，這表明果實處于不斷發(fā)育變化中。HT1和HT3上調基因比例明顯增加，核桃發(fā)育最快，HT2和HT3相比差異基因減少，可能與后期果實成熟代謝減慢有關。

2.7 GO富集分析

Gene Ontology（簡稱GO）是一個國際標準化的基因功能分類體系，提供了一套動態(tài)更新的標準詞匯表（controlled vocabulary）來全面描述生物體中基因和基因產物的屬性。GO總共有3個ontology（本體），分別描述基因的分子功能（molecular function）、細胞組分（cellular component）、參與的生物過程（biological process）[15]。將不同成熟期“三臺”核桃的差異基因按照其功能進行GO分類，將每個分類數(shù)量最多的2種term和數(shù)量列入圖6中。在生物學過程（biological process）中，metabolic process（代謝過程）、cellular process（細胞過程）和single-organism process（單組織過程）所占的比例最多;細胞成分（cellular component）中，cell（細胞）、cell part（細胞部分）、organelle（細胞類脂質）所占的比例最多;分子功能（molecular function）中，數(shù)量最多的為catalytic activity（催化活性）和binding（結合）。

3 結論

（1）以4個不同成熟期（硬核期、油脂轉化期、油脂積累期、果實成熟期）云南主栽的“三臺”核桃為材料，轉錄組分析不同成熟期基因表達關系。共有53 868條unigenes得到了注釋，其中33 407條unigenes注釋在Nr數(shù)據(jù)庫中;31 725條unigenes注釋在KEGG數(shù)據(jù)庫中;Swissprot數(shù)據(jù)庫中注釋到19 842條unigenes;COG數(shù)據(jù)庫中注釋到16 655條unigenes。四大數(shù)據(jù)庫共注釋到33 753條unigenes，占總unigenes的62.66%。

（2）不同成熟期“三臺”核桃差異基因兩兩比較中，CK-vs-HT1、CK-vs-HT2、CK-vs-HT3、HT1-vs-HT2、HT1-vs-HT3、HT2-vs-HT3分別有 2 417、2 703、4 166、1 256、5 347、3 283條上調表達的差異基因和9 515、12 455、9 498、1 751、3 132、1 365條下調表達的差異基因。

（3）在“三臺”核桃生物合成途徑中，以次生代謝產物生物合成途徑注釋到的兩端不能再延長的序列最多，在脂肪酸代謝途徑中富集到的差異表達基因最多，其中acc系列酶及FAB系列基因可能是“三臺”核桃油脂和脂肪酸生物合成關鍵基因，其后續(xù)通過表達載體構建和蛋白異源表達等對基因功能進行驗證。

參考文獻

[1] 方文亮，寧德魯.云南核桃[M].北京：科學出版社，2019.

[2] YOU F M，DEAL K R，WANG J R，et al.Genome-wide SNP discovery in walnut with an AGSNP pipeline updated for SNP discovery in allogamous organisms [J].BMC Genomics，2012，13：1-16.

[3] RAI A，KAMOCHI H，SUZUKI H，et al.De novo transcriptome assembly and characterization of nine tissues of Lonicera japonica to identify potential candidate genes involved in chlorogenic acid，luteolosides，and secoiridoid biosynthesis pathways[J].Journal of natural medicines，2017，71（1）：1-15.

[4] 張楠.核桃胚脂肪積累期轉錄組分析[D].泰安：山東農業(yè)大學，2014.

[5] 楊麗，陳虹，潘存德，等.核桃種子油脂轉化期轉錄組分析[J].果樹學報，2017，34（9）：1084-1094.

[6] LI Y T，MA S M，WANG Y F，et al.The dynamics of fat，protein and sugar metabolism during walnut（Juglans regia L.）fruit development [J].African journal of biotechnology，2012，11（5）：1267-1276.

[7] 潘教文，李臻，王慶國，等.NaCl處理谷子萌發(fā)期種子的轉錄組學分析[J].中國農業(yè)科學，2019，52（22）：3964-3976.

[8] 張坤.紅地球葡萄延后栽培生育后期樹體與果實的水分關系研究[D].蘭州：甘肅農業(yè)大學，2018.

[9] 王燕.蘿卜鉛（Pb）脅迫響應的分子機制研究[D].南京：南京農業(yè)大學，2014.

[10] LANGMEAD B，SALZBERG S L.Fast gapped-read alignment with Bowtie 2[J].Nature methods，2012，9（4）：357-359.

[11] LI B，DEWEY C N.RSEM：Accurate transcript quantification from RNA-Seq data with or without a reference genome[J].BMC Bioinformatics，2011，12：1-16.

[12] 魯強.三椏苦種子休眠的解除方法和初步機理研究[D].廣州： 廣州中醫(yī)藥大學，2019.

[13] 劉文浩.熏硫及脫硫處理對龍眼貯藏生理的影響及轉錄組研究[D].廣州：華南農業(yè)大學，2016.

[14] 李文彬.‘紅陽獼猴桃果實發(fā)育轉錄組及花青素累積機理研究[D].武漢：中國科學院研究生院（武漢植物園），2015.

[15] 白琳云.靈武長棗嫁接與根蘗繁殖植株生長、果實特性的比較[D].銀川：寧夏大學，2019.

基金項目 云南省重大專項科技計劃項目（2018ZG003）;云南省科技創(chuàng)新人才培養(yǎng)項目（2016HB004）。

作者簡介 賀娜（1982—），女，湖北京山人，副研究員，碩士，從事經濟林栽培及精深加工研究。

*通信作者，研究員，博士，從事植物資源評價及木本油料精深加工研究。

收稿日期 2020-05-07