基于剛毛檉柳轉錄組測序的EST-SSR標記識別與開發

李佳彬，黃 蕾，張雅楠，賈媛媛，田蕓蕓，張 雷，黨振華

（1. 內蒙古大學生態與環境學院 / 蒙古高原生態與資源利用教育部重點實驗室 / 內蒙古草地生態學重點實驗室，內蒙古 呼和浩特010021；2. 內蒙古大學生命科學學院 / 牧草與特色作物生物技術教育部重點實驗室，內蒙古 呼和浩特 010021；3. 內蒙古自治區林業科學研究院，內蒙古 呼和浩特 010010；4. 內蒙古大青山森林生態系統定位觀測研究站，內蒙古 呼和浩特 010010）

剛毛檉柳(Tamarix hispida)是中亞的廣布木質鹽生植物，是檉柳科(Tamaricaceae)檉柳屬(Tamarix)最耐鹽堿的種類之一，常分布于我國荒漠或半荒漠地帶的鹽土或低濕沙區[1]。它不僅在水土保持、防風固沙和綠化造林方面具有重要的生態價值，豐富的營養與化學成分還賦予其很高的藥用和經濟價值[2-3]。目前對剛毛檉柳的研究主要集中在形態解剖[4]、群落生態[5]、生理生態[6-7]、植物化學[8]、分子生態[9-10]等方面。然而，作為荒漠地區的主要資源樹種，對其適應生境的內在分子基礎還知之甚少，亟需深入了解和研究，以進一步挖掘蘊藏在該物種中的特有遺傳資源。

SSR (simple sequence repeat)，即簡單序列重復，指重復基序為1～6 個核苷酸的串聯重復序列，廣泛分布于真核生物基因組中，每隔10～15 kb 就存在1 個SSR 位點[11]。目前，SSR 分子標記已被用于遺傳圖譜的構建、遺傳多樣性分析、品種鑒定及分子育種等研究領域[12-14]。依據SSR 在基因組中的分布，可將其分為基因組SSR 和表達序列標簽SSR(EST-SSR)[15]。基因組SSR 位于基因組非編碼區，通常情況下，等位基因數目及多態性均相對較高，而EST-SSR 由于位于基因編碼序列內部，具有相對保守的特點。大量研究表明，EST-SSR 的變異能改變基因的活性或調節基因表達，故常與基因的功能密切相關，可影響相關的生物和細胞過程，如蛋白質結構、傳感和信號傳遞及基因轉錄等[16-17]。例如，編碼區(CAG/CTG)n 重復的擴增或減縮可導致蛋白質錯誤折疊或基因表達異常[18]；(CCTG/CAGG)n 重復與復制起點的距離可能在決定這些重復序列的遺傳不穩定性方面發揮重要作用[19]；(GAA/TTC)n 重復可形成三聚體，相關結構被稱為“粘性DNA”，在模板上捕獲RNA 聚合酶，從而阻斷轉錄延伸[20-21]。因此，EST-SSR 可能是生物體適應環境變化的分子裝置[22]，鑒定和深入分析多態EST-SSRs 的功能可為探索植物適應性進化機制提供新見解[23-24]。

CandiSSR 是一款近年開發的多態性微衛星識別軟件，可同時比較相同或近緣物種多個樣本的基因組或轉錄組序列，進行大規模多態性SSR 的識別及分析[25]。本研究通過轉錄組測序構建了5 個不同樣地剛毛檉柳的基因序列集，利用CandiSSR 軟件系統識別了該物種多態EST-SSR 位點；隨后對所識別多態SSR 的數目、類型、出現頻率、基因中的位置及關聯基因的功能進行了全面分析；經驗證，從15 個隨機挑選的SSR 中共開發出13 個多態性SSR 標記。本研究可為進一步鑒定由EST-SSRs 變異驅動的該物種生境適應機制提供素材，為開展其遺傳多樣性評價和種植資源開發奠定基礎。

1 材料及方法

1.1 植物材料

2019 年9 月，在內蒙古阿拉善盟額濟納旗采集了5 個樣地剛毛檉柳(Tamarix hispida)的幼嫩葉片，詳細采樣地點的地理信息如表1 所列。每個樣地中每隔5～10 m 采集20～24 個個體。樣本采集后裝入凍存管，并迅速置于液氮速凍，帶回實驗室備用。

1.2 試驗方法

1.2.1 總RNA 提取及質量檢測

每個樣地選擇1 株植物用于總RNA 的提取。使用RNA plant plus 植物總RNA 提取試劑(DP437，天根生化科技有限公司，北京)進行。利用1%瓊脂糖電泳和超微量紫外分光光度計(NanoDrop 2000c，賽默飛世爾，美國)檢測總RNA 的質量。檢測合格后的RNA 樣本置于干冰中寄送至安諾基因(安諾優達基因科技有限公司，北京)。

表 1 5 個剛毛檉柳種群的地理位置信息Table 1 Location information for five Tamarix hispida populations

1.2.2 cDNA 文庫構建及轉錄組測序

利用安捷倫2100 RNA Nano 6000 Assay Kit (Agilent Technologies, 美國) 對RNA 樣本進行完整性檢測，RNA 完整值(RNA integrity number，RIN)達到7.0 以上。然后由測序公司進行測序文庫構建，簡要流程為：用帶有Oligo(dT) 的磁珠富集mRNA，向得到的mRNA 中加入片段緩沖液使其成為短片段，再以片斷后的mRNA 為模板，用六堿基隨機引物合成cDNA第一鏈。構建好的文庫利用Illumina HiSeqTM4000測序平臺進行測序。

1.2.3 轉錄組de novo 組裝及注釋

對原始序列進行過濾得到質量較高的Clean Reads，隨后采用Trinity (v2.4.0)分別對5 個樣本的Clean Reads 進行de novo 組裝，得到每個樣本轉錄組的Unigenes，然后進一步將各轉錄組Unigenes 進行拼接和去冗余，得到非冗余All-Unigenes。通過Blast、HmmScan、SignalP、TmHMMP 等 工 具 對All-Unigenes 進行功能注釋，采用Trinotate (v3.0.2)整合功能注釋信息。利用TransDecoder (v3.0.1)對unigene的編碼區進行鑒定，確定蛋白質編碼序列(CDS)。

1.2.4 多態SSR 位點識別與定位

利 用 CandiSSR (https://github.com/xiaenhua/CandiSSR)識別5 個剛毛檉柳轉錄組中的多態性SSRs[25]。篩選標準：包含2、3、4、5、6 重復單元的SSR 至少出現次數分別是6、5、5、4 和4 次。對于具有全長CDSs 的Unigenes，可根據SSRs 與相應基因起始(ATG) 和終止(TAA、TAG、TGA) 密碼子的相對位置來分析SSRs 的位置。對于組裝為非全長CDSs 的Unigenes，通 過BLAST(http://blast.ncbi.nlm.nih.gov/Blast.cgi)識別Genbank 中的全長同源基因，然后根據其在查詢序列中的位置預測SSRs 的位置。對于識別出的SSR，采用Primer 3.0 進行引物批量設計，設計原則：①引物長度18～24 bp；②退火溫度在53～62 ℃；③PCR 產物長度100～200 bp；④GC含量在40%～60%。

1.2.5 DNA 提取及多態SSR 標記開發

利用植物基因組DNA 試劑盒(DP305，天根生化科技有限公司，北京)從30 個剛毛檉柳樣本中提取基因組DNA，并用1% 瓊脂糖凝膠電泳和NanoDrop 2000c (賽默飛世爾，美國)進行DNA 質量及濃度的檢測。隨機挑選15 個多態SSR 位點進行驗證，PCR 擴增引物由擎科新業生物技術有限公司(北京)合成。

使用ABI2720 PCR 儀進行PCR 擴增，選擇25 μL反應體系：50 ng·μL-1DNA 模板1 μL，12.5 μL Premix Taq (寶生物工程有限公司，大連)，上、下游引物各0.5 μL(10 μmol·L-1)，10.5 μL 的ddH2O。PCR 擴增程序：94 ℃預變性5 min；94 ℃變性30 s，退火30 s，72 ℃延伸30 s，30 個循環；72 ℃ 10 min 終止反應。擴增產物用2% 瓊脂糖檢測。對于擴增成功的引物，在5′端加入3 種熒光染料(6-carboxy-fluorescine，hexachloro-6-carboxy-fluorescine，6-carboxy-X-rhodamine)中的一種進行標記，然后用相同的PCR 程序對所有檢測個體再次擴增，結果在ABI 3 730 DNA 分析儀(Applied Biosystems，北京) 上分析，毛細管電泳的內標為GS500LIZ。從上述驗證成功的引物中，選擇5 對引物(TR25868、TR19384、TR18283、TR23634、TR23597)的PCR 產物，依次進行2% 瓊脂糖凝膠電泳、凝膠回收及DNA 片段純化(愛思進生物技術有限公司，杭州)，將回收產物與pMD19-T 載體(寶生物工程有限公司，大連)連接并轉化大腸桿菌(全式金生物技術有限公司，北京)，經菌液PCR 鑒定后，選取陽性克隆送至擎科新業生物技術有限公司(北京)進行測序。

1.2.6 數據處理

采用GeneMarker (v2.6)對毛細管電泳峰圖進行基因型判讀及基因型統計，每個位點用GenAlEx(v6.5)對等位基因數(Na)、觀測雜合度(Ho)和期望雜合度(He)進行計算，然后用PowerMarker (v3.0)測量多態性信息含量(PIC)。通過GenePop (v4.7)查看各位點是否偏離哈迪-溫伯格平衡和連鎖不平衡。

2 數據分析

2.1 測序和de novo 組裝結果

經測序，5 個樣本的Clean Reads 最小為44 180 102條，最大47 575 708 條，平均Q30為93.13%。對上述Clean Reads 進行組裝后，每個個體分別獲得超過31 000 條Unigenes，平均長度范圍在1 043.24～1 139.67 bp，平均N50為1 913～1 970 bp (表2)。進一步組裝和去冗余后，共得到72 661 條All-Unigenes，總長度為65 674 878 個核苷酸，平均長度為903.85 bp，平均GC 含量為40.62%，N50長度為1 578 bp。

2.2 多態SSRs 的識別

在5 個剛毛檉柳轉錄組中，共鑒定到1 187 個多態性SSRs 位于1 123 個基因序列中，共有154 個SSR 重復單元類型，二、三、四、五、六核苷酸重復出現頻率有較大差異(圖1)。最常見的SSRs 是三核苷酸重復類型，共646 個(54.42%)，最豐富的重復序列類型有AGC/GCT(105，16.25%)、AAT/ATT(89，13.78%)、AAG/CTT (85，13.16%)等；其次是二核苷酸重復，有447 個(37.66%)，主要基元以AG/CT (217，48.55%)和AT/AT(177，39.60%)為主；四、五、六核苷酸重復類型的數量較少，占總SSRs 的7.92%，3 種重復類型分別以ATAA/TTAT(8，22.86%)、AAATA/TATTT(8，33.33%)、CCCTCT/AGAGGG (5，14.29%)為主。

表 2 測序數據產出與組裝結果Table 2 Sequencing outputs and assembly results

圖 1 SSR 重復類型及頻率示意圖Figure 1 Diagram of SSRs' motif types and frequencies

2.3 SSR 位置預測及頻率分析

在多態SSR 對應的Unigene 中，有856 個與公共數據庫中的已知蛋白質有同源比對結果。887 個SSRs 分別位于829 個Unigenes 的編碼區(CDSs)和非翻譯區(UTRs) 內。其中500 個位于CDSs 中，176 個位于3′UTRs，211 個位于5′UTRs。三核苷酸重復序列(91.40%)多位于CDSs 區，AGC/GCT 相對豐富；二核苷酸多位于UTRs，在3′UTRs 中占61.93%，其中AT/AT 相對豐富；5′UTRs 中占62.56%，AG/CT相對豐富(圖2)。

圖 2 SSR 分布分析圖Figure 2 Simple sequence repeat distribution analysis diagram

2.4 多態SSR 關聯基因的GO 及KEGG 功能注釋

GO 功 能 注 釋 顯 示，含685 個SSRs 的635 條Unigenes 與GO 數據庫中的5 582 個功能基因有比對結果(圖3)。它們被劃分為3 大類，分別為2 428個生物過程(43.50%)、1 253 個細胞組分(22.44%)和1 901 個分子功能(34.06%)。進一步將三大功能細分為1 248 種亞類。其中，生物過程包括743 個功能亞類。“調節轉錄”(regulation of transcription)類別在3′UTRs 與5′UTRs 內含有SSR 的Unigene 中所占比例最大，分別占GO 注釋總Unigene 的5.82%和8.27%；“轉 錄”(transcription) 類 別 在CDSs 含 有SSR 的Unigene 中最多(占10.17%)；細胞組分包括161 個功能亞類，代表性功能為“細胞核”(nucleus)，分別在3′UTRs、5′UTRs、CDSs 中 含 有SSR 的Unigene 中占20.59%、19.01%和29.77%；分子功能包括344 個功能亞類。“轉錄因子活性”(transcription factor activity)類別在3′UTRs (10.63%)及CDSs (13.99%)含有SSR的Unigene 中 最 多，“ATP 結 合”(ATP binding)在5′UTRs (8.80%)含有SSR 的Unigene 中最多。

KEGG 注釋分析顯示，120 條含有多態SSR 的Unigenes 參與到110 個通路中。對于包含多態SSR的5′UTRs、CDSs 和3′UTRs，參與最多的是“代謝途徑”(metabolic pathways)，分別在5′UTRs、CDSs 和3′UTRs 含有SSR 的Unigene 中占13.08%、12.60%和11.76%。其次是“植物激素信號轉導”(plant hormone signal transduction)，分別在5′UTRs、CDSs 和3′UTRs含有SSR 的Unigene 中占8.41%、11.02%和7.06%。

2.5 多態性SSR 引物篩選與標記開發

利用1 187 個多態SSRs 的側翼序列，成功為1 182 個SSRs 設計出了PCR 擴增引物。在隨機選取的15 個多態SSRs 引物中，13 個成功擴增出多態位點(表3)。為進一步驗證多態SSR 的可信度，選擇上述成功擴增的5 個位點，分別進行了SSR 片段的回收、克隆及測序。結果如圖4 所示，測序結果(圖4C)與基于高通量測序組裝的基因序列(圖4A)完全吻合，且利用CandiSSR 識別到的SSR 重復單元變異與毛細管電泳檢測到的多態性一致(圖4B)。同時發現，個別SSR重復內存在堿基替換事件，例如：TH7_DN13421、TR_25868_G7：A→C(圖4A)；TR_25868_G16：G→C(圖4C)。

圖 3 含有多態SSR 位點基因序列的GO 功能注釋結果Figure 3 Results of the gene ontology function annotation of the polymorphic SSR-containing sequence

在這些SSR 中，共檢測到73 個等位基因，每個位點得到3～8 個等位基因(平均5.615 個)。期望雜合度(He) 的范圍為0.301～0.786，觀測雜合度(Ho)的范圍為0.133～0.800，均值分別為0.619 和0.476。多態性信息含量(PIC)值范圍為0.283～0.744，平均值為0.565 (表3)，所有位點均符合哈迪-溫伯格平衡(P ＞ 0.05)，且未檢測到連鎖不平衡現象。

3 討論

高通量測序技術極大推進了獲得非模式物種基因序列的速度。利用這些數據，再結合一些生物信息學軟件(MISA[26]和SSR Finder[27]等)，研究人員可快速識別出相應物種數以千計的SSR 位點。無論開展種群遺傳學研究還是對特定SSR 的功能進行分析，首先都要得到一組多態性的SSR 標記。然而，對基于高通量測序數據的海量SSR 鑒定多態性，仍然是該標記技術應用和研究的瓶頸[24]。本研究通過構建相同物種不同分布區個體的轉錄組數據集，利用CandiSSR 在5 個剛毛檉柳轉錄組中共識別出1 187個多態性SSRs，隨機挑選的15 對引物中有13 個為多態性SSR 標記，成功率近87%，大大提高了該物種多態性SSR 標記開發的效率。遺傳參數分析表明，這些SSR 的平均PIC 為0.565，平均He 為0.619，平均Ho 為0.476，呈現高度多態性，可用于該物種后續的種群遺傳學和適應進化研究。同時，本研究所識別的SSR 數目高于用相同軟件鑒定的南北兩地區銀縷梅(Parrotia subaequalis)的497 個SSRs[28]、四合木(Tetraena mongolica)的881 個SSRs[29]，低于茶樹(Camellia)鑒定的1 663 個SSRs[30]，一定程度反映了該物種SSR 變異相對豐富及不同物種SSR 含量的差異。當然，這一差異也可能是因檢測個體數、識別參數不同造成。在剛毛檉柳多態SSR 中，三核苷酸出現頻率最高(46.93%)，且91.40%的三核苷酸重復位于CDSs 中，這與荔枝(Litchi chinensis)和刺梨(Rosa roxburghii)的研究結果相符[31-32]。這可能是因為密碼子由3 個核苷酸構成，若在編碼區發生突變，將引起最輕的蛋白質序列突變[33]。在本研究中，CDSs 中的SSRs 表現出對三核苷酸的強烈偏倚性，最豐富的三核苷酸重復為AGC/GCT，其次是AGG/

CCT、ACC/GGT、ATC/GAT 和AAG/CTT。這與Sonah等[34]認為的AAG、AAC、ATC、AGC、AGG 和ACG重復是雙子葉植物外顯子中最常見SSR 類型的觀點一致，一定程度反映了本研究SSR 頻率分析的正確性。

表 3 開發多態EST-SSR 標記的信息Table 3 Detailed information of the developed polymorphic EST-SSR markers

剛毛檉柳生態幅較廣，不同分布區的水分、溫度、地形、土壤屬性等均存在較大差異，不同生境下的種群很可能進化出不同的適應機制。張娟等[35]用隨機擴增多態性DNA (random amplified polymorphic DNA, RAPD)方法描述了分布在新疆的9 個剛毛檉柳居群的遺傳分化及遺傳結構；張道遠等[5]分析了剛毛檉柳在不同干旱情況下脯氨酸及可溶性糖含量的不同。EST-SSR 可能是生物快速適應外界環境變化的分子裝置[36]。本研究共鑒定出1 187 個多態SSR 分布在1 123 個轉錄本，它們可能在剛毛檉柳局部適應過程中發揮重要作用。在含有SSR 的序列中，267 個Unigenes 與公共數據庫中的已知基因沒有比對結果，可能代表了該物種獨有的遺傳資源。編碼區SSR 的變異可引起基因功能的獲得或喪失，5′UTRs 中SSR 的變異可影響基因的轉錄和翻譯，3′UTR 的SSR 則影響mRNA 的剪接[37]。對剛毛檉柳含EST-SSR 基因的GO 和KEGG 功能注釋發現，它們主要歸類于“調節轉錄”、“轉錄因子活性”、“序列特異性DNA 結合”等GO 條目，“代謝途徑”和“植物激素信號轉導”等KEGG 代謝通路，一定程度反映了該植物在適應不同生境時，眾多生物過程和代謝通路的基因很可能發生了較高頻率的變異或表達的調節，且EST-SSRs 的變異可能在這一過程中發揮重要作用。