基于簡化基因組數據開發嶺南青岡微衛星引物

寧 馨 姜小龍 鄧 敏 徐剛標

(1.中南林業科技大學林木遺傳育種實驗室，長沙 410004；2.中國科學院上海辰山植物科學研究中心,上海辰山植物園,上海 201602)

微衛星又稱為短串聯重復序列(short tandem repeats,STRs)或簡單重復序列(simple sequence repeat,SSR)，通常由1～6個核苷酸的串聯重復片段構成，一般重復幾次到幾十次[1]。跟核苷酸的遺傳多樣性及突變率相比，微衛星具有更高的多態性及突變速率。因此，自20世紀90年代開始，微衛星標記就一直被廣泛用于法醫、親子鑒定、保護遺傳學、譜系地理學及群體遺傳學中[2～3]。不過由于微衛星的側翼序列在種間也可能存在變異，因此，針對某一物種設計的引物通常只適用于這一物種及其近源種。隨著近幾年新一代測序技術的快速發展，基于高通量測序技術可在短時期內以較低成本獲得大量SNP(Single nucleotide polymorphism)[4]。例如，近幾年出現的簡化基因組測序技術(restriction-site associated DNA sequence,RAD-seq)現已在群體遺傳學、保護生物學及譜系地理學中開始廣泛使用[5]。目前微衛星在保護生物學及群體遺傳學中還具有一定的優勢。首先，微衛星標記能有效地擴增DNA降解嚴重的樣本。雖然RAD-seq能以較低的成本得到大量的直系同源SNPs，但是其對提取DNA質量的要求較高，很多長期硅膠干燥保存的分子樣和標本上收集的樣本上提取的DNA及古DNA往往不能滿足二代測序的建庫要求；其次，RAD-seq中所得到SNPs的錯誤率跟DNA降解程度呈現正相關[6]；再次，微衛星的成本低廉。盡管利用微衛星標記得到的變異位點有限，但對于較少個體的瀕危物種群體遺傳學研究，微衛星標記仍是一個有效、快捷的手段。并且微衛星數據易于分析。使用普通的個人電腦就可以對其進行數據分析和處理，而RAD-seq產生的數據量較大，通常需要服務器或者工作站進行數據分析，并需要有熟練的數據分析技能。

經典微衛星引物開發的方法有文庫法、磁珠富集法和省略篩庫法。文庫法是先用超聲波或酶切方法將DNA切成小片段后建庫，然后用含有SSR的探針與文庫雜交，篩選出和探針有較強信號的克隆，此方法開發引物所得到的SSR陽性克隆比率低，不過操作簡單、容易掌握。富集法通過建立和篩選微衛星富集文庫，能極大地提高文庫法中SSR克隆得率低的問題。例如，磁珠富集法用生物素標記重復序列特異探針并與基因組DNA酶切片段雜交，再用包被鏈親和素的磁珠對生物素進行磁力吸附，從而完成重復序列目標片段的富集。省略篩庫法使用5′端帶7個簡并堿基的SSR引物(5′-KKVRVRV(CT)6-3′)對基因組SSR區域進行PCR富集后進行測序[7]。此方法不需要篩庫去識別含SSR的克隆，可以進一步提高有多態位點SSR的得率。近年來，隨著測序技術的快速發展，高通量測序技術被廣泛用來開發微衛星引物。例如，使用物種轉錄組序列[8～10]、基因組序列[11]和簡化基因組序列[12～13]來進行開發。這些方法直接利用測序序列設計SSR引物，簡化了實驗流程，并可得到更多的微衛星位點。先前使用簡化基因組序列開發SSR引物都是對雙端測序數據(2×150 bp)拼接后獲得較長的序列，然后再開發引物[14]。SSR擴增片段長度通常多為100 bp，因此，是否能直接使用簡化基因組測序中150 bp的R1端序列直接設計引物呢？如果可行，則能進一步簡化引物設計的流程。

本次研究我們基于前期獲得的嶺南青岡(QuercuschampioniiBenth)RAD-seq數據的R1端序列，進行SSR引物的開發，并作驗證。嶺南青岡是中國南方及中南半島亞熱帶中山云霧林的代表種，主要分布于海拔100～1 700 m的森林[15]。本研究旨在證明基于簡化基因組測序獲得的R1端序列(150 bp)是否能夠開發出足夠的SSR引物用于群體遺傳分析。

1 材料與方法

1.1 樣品來源

本研究所用的4個樣本來自于Jiang等文章中4個不同群體的嶺南青岡簡化基因組數據，分別為廣東省龍門縣南昆山、廣州省佛岡縣觀音山、廣西省寧明縣公母山和廣東省羅博縣象頭山[16]。

1.2 數據過濾與引物設計

數據的過濾與引物設計按照如下流程進行：

(1)使用pyRAD[17]對序列進行質控過濾，去除低質量的reads。使用SciRoKo v3.4[18]提取包含有微衛星重復序列的reads。

(2)使用pyRAD對提取出來的序列進行個體內和個體間的聚類。由于我們關注的是微衛星中重復堿基數目的變異，這一變異類型在聚類時會包含大量的插入/缺失堿基(Indel base)。因此，為了在聚類中盡可能發現所有的微衛星多態性，將分析的聚類閾值(Clustering threshold)、樣本內和樣本間允許的最大插入缺失(Indel)堿基數分別設置為0.7，40和40。

(3)過濾掉序列中個體間存在堿基轉換/顛換的位點。樣本間聚類以后，除了微衛星堿基重復數具有差異以外，側翼序列還存在堿基轉換/顛換的位點。若在存在變異位點的側翼序列區設計引物會影響后續PCR擴增的成功率。因此，使用python腳本將存在堿基轉換/顛換的位點過濾掉。同時，設計微衛星引物時需要引物兩側都要具有一定長度的側翼序列，故只把微衛星重復堿基位于序列中26～129 bp的序列提取出來。這樣，微衛星序列的左側側翼序列至少有25 bp，右側側翼序列至少有21 bp可用于進行后期的引物設計，得出分析的流程圖見圖1。

(4)使用primer premier 5.0設計引物(Premier Biosoft International,CA,USA)。在設計引物時，只挑選個體間微衛星重復數存在差異的樣本進行設計(見圖2)。

圖1 基于RAD-seq開發SSR引物流程圖 reads中紅色表示微衛星重復序列Fig.1 Workflow of SSR Primer development based on the reads obtained from RAD-seq Red segment represents the microsatellite sequence in reads

圖2 微衛星引物RAD195設計區域Fig.2 Primer region for microsatellite RAD195

1.3 實驗驗證

為了減少熒光引物合成的成本，在正向引物的5′端添加一個接頭序列M13(序列為:5′-TGTAAAACGACGGCCAGT-3′)，先用接帶M13接頭的引物進行擴增，然后添加帶熒光的M13序列進行二次擴增[19]，PCR擴增體系如下：第一輪PCR擴增體系為：10×PCR Buffer(Mg2+free) 2 μL，dNTP(2.5 mmol·L-1) 0.5 μL，Mg2+(25 mmol·L-1) 2 μL，正向引物(10 μmol·L-1) 0.2 μL，反向引物(10 μmol·L-1) 0.7 μL，Taq酶(5 U·μL-1) 0.2 μL，基因組DNA模板(10 ng·μL-1) 2 μL，牛血清蛋白BSA(20 mg·mL-1) 0.2 μL，ddH2O 12.2 μL，PCR總體積共20 μL。PCR反應程序為:94℃預變性3 min、94℃變性30 s、55℃退火30 s、72℃延伸30 s，35個循環，72℃延伸7 min，12℃保存。第2輪擴增體系為：10×PCR Buffer(Mg2+free) 0.2 μL，Taq酶(5 U·μL-1)0.08 μL，熒光標記M13引物(10 ng·μL-1)0.4 μL，ddH2O 1.32 μL，PCR總體積共2 μL。PCR反應程序為:94℃預變性5 min、94℃變性30 s、53℃退火30 s、72℃延伸30 s，16個循環，72℃最后延伸10 min，12℃保存。

共對3個群體的36個樣本進行擴增。擴增結果電泳檢測后送美吉生物科技有限公司進行毛細管電泳分型。使用GeneMarker v2.2.0(SoftGenetics,PA USA)對SSR分型結果進行讀取和校對。

2 結果與分析

2.1 RAD-Seq測序數據質量評估

所選4個樣本的原始序列(Raw reads)數在240～640萬條。除去低質量的數據后序列數在180～520萬條(見表1)。這些序列中包含微衛星重復片段的序列在46 000～56 000條。

2.2 多態性SSR引物的開發

對樣本內進行聚類、過濾掉旁系同源序列后，每個樣本包含的位點數在5 546～15 844。樣本間聚類后，得到所有樣本都具有的位點(一致性位點)1 158個，并得出詳細的樣本序列信息(見表2)。使用python腳本過濾掉1 158個一致性位點中包含轉換/顛換的序列以后，共得到308個位點。進一步過濾掉微衛星重復序列位于reads前25或129 bp以后位置的位點后，得到位點數186個。在這186對引物中挑選出微衛星是2堿基或3堿基重復并最大重復數大于6次,而且在4個樣本中至少存在2種不同重復數的引物。共挑選出了40對，使用Primer Premier 5.0設計出了25對引物用于實驗驗證(見表2)。

表1 用于SSR引物開發的樣本信息

表2 開發用于微衛星多態性篩選的引物序列信息

使用設計好的25對引物對36個個體進行擴增。所有引物的PCR擴增條件和體系見1.3。其中，17對引物被成功擴增出來。這17對引物的等位基因數共為106個，每個引物的等位基因數在2～12，平均為6.2。引物的期望雜合度和觀測雜合度分別為0.19～0.88和0.11～0.76。得出引物的詳細信息(見表3)。

表3 研究成功擴增引物的多態性信息

Table 3 Polymorphism information of the amplified primers

引物名Primername個體數Individualnumber等位基因數Allelenumber期望雜合度Expectedheterozygosity觀測雜合度ObservedheterozygosityRAD131730.640.59RAD393160.610.68RAD543380.830.76RAD593440.410.21RAD693440.190.21RAD1012340.400.30RAD10932120.880.75RAD1122950.660.48RAD1232470.770.63RAD1803140.630.55RAD1903390.790.70RAD19134110.850.74RAD1953480.770.74RAD1983320.210.18RAD2023190.810.74RAD2113440.520.44RAD2133560.700.11

3 討論

本研究基于簡化基因組測序方法(RAD-seq)獲得的R1端序列成功開發了10多對具有遺傳多樣性的嶺南青岡微衛星引物。根據合成的25對SSR引物，一次性擴增出來17對具有多態性的引物。這17對引物的等位基因數共106個，SSR位點的等位基因數在2～12，平均數為6.2。引物的期望雜合度和觀測雜合度分別為0.19～0.88和0.11～0.76。雖然獲得的微衛星引物比使用轉錄組開發的少，但是研究獲得的等位基因數目能滿足進行物種的群體遺傳學、保護生物學及譜系地理學研究的需要。隨著測序技術的發展及費用的減低，可以使用PE250對簡化基因組數據進行測序。跟本研究采用的PE150測序模式相比，更長讀長的PE250測序模式能有更多位點符合微衛星引物的設計。不過對于超大基因組(大于3G)的物種，例如百合科植物的基因組通常達到幾十G[20]，為了滿足測序位點覆蓋度的要求可能需要測十幾甚至幾十G的數據，這極大地增加了測序的成本。因此，我們不建議使用簡化基因組測序的方法開發具有超大基因組物種的SSR引物。本次研究中設計的引物都是物種的直系同源序列，但還存在8對引物不能一次性擴增出來。由于正向引物都加了M13接頭序列，會對引物的結構和結合特異性有一定影響，可能導致擴增的失敗，雖然通過直接合成熒光引物可以消除這一因素的影響，但會增加引物合成的成本。

跟磁珠吸附法開發微衛星引物相比，本次SSR開發方法簡化了實驗的流程、減少了工作量，并且所需費用也更低。跟利用轉錄組序列開發微衛星引物相比，使用簡化基因組數據的R1端序列開發引物所需要的測序成本、數據的分析能力及計算機資源要求更低[21]。轉錄組數據的拼接需要計算機具有較大的內存，通常需要服務器才能滿足數據分析的需求。本研究的開發方法只需提取含有微衛星位點的序列進行比對即可，這些位點只占所有序列的1%左右。因此，普通的個人電腦即可快速地(1 d以內)完成整個分析。同時，轉錄組序列位于基因的編碼區，比較保守。為了獲得足夠的遺傳變異信息需要進行大量的實驗篩選。而簡化基因組序列是均勻分布于物種基因組中，具有更高的遺傳變異。再次，本研究使用了4個樣本進行樣本間的聚類，能有效提高篩選出具有多態性的微衛星引物的概率，極大地減少了微衛星引物后期篩選實驗的工作量。對于開發引物用于群體遺傳學、保護生物學或譜系地理學分析時，此方法可作為有效的候選方法之一。