利用全基因組重測序技術研究182份大麥和青稞的基因組結構變異

徐齊君，王玉林，楊春葆，扎 桑，于明寨，原紅軍

(1.省部共建青稞和牦牛種質資源與遺傳改良國家重點實驗室,拉薩 850002；2.西藏自治區農牧科學院 農業研究所，拉薩 850032)

青稞(Hordeumvulgarevar.nudum)作為主要糧食作物和牲畜飼料主要種植于青藏高原，青稞生產分別占西藏糧食和農區飼草總產的80%和60%以上，對保障西藏糧食安全、促進畜牧發展具有重要意義。其栽培歷史可追溯到3 500年前。在過去的幾年里，研究人員已經陸續公布青稞[1-3]和大麥(HordeumvulgareL.)[4]的細胞核和細胞器基因組，為基于參考基因組的結構變異分析提供了重要基礎。結構變異(Structural Variation，SV)是一種基因組序列的重排現象，影響基因組片段大小超過50個堿基，不同于單核苷酸多態性(Single Nucleotide Polymorphisms，SNPs)和短插入/缺失(Insertions and Deletions,INDELs)。SV包括大片段的缺失(Deletion,DEL)、大片段的插入(Insertion,INS)、倒位(Inversion,INV)、易位和復制(Duplication,DUP)、拷貝數變異(Copy Number Variant,CNV)，CNV為不平衡SV[5-6]。SV最初的檢測方法是基于顯微鏡下的大規模染色體變化的核型觀察方法，例如非整倍體[7-8]、染色體重排[9]和CNVs[10-11]。二代測序技術(Next Generation Sequencing,NGS)為SV的檢測提供了一個前所未有的機會，可以從全基因組水平上闡明表型差異和復雜遺傳變異的分子機制。與傳統的實驗技術相比，NGS技術已被廣泛應用于動物和植物的遺傳學研究，特別是SVs的研究。對人類的研究表明，SVs與精神分裂癥[12]、癌癥[13]和復雜的遺傳紊亂疾病[14]有關。據估計，在人類基因組的遺傳變異中，CNVs占比17.7%[15]。對植物基因組的廣泛研究表明，SVs在植物適應進化和表型多樣性中起著重要作用，如擬南芥[16]、黃瓜[17]、水稻[18-19]、高粱[20]和大豆[21]。基于NGS技術，正朝著遺傳變異發現的新階段前進，該階段的重點是以單堿基分辨率來識別SVs在基因組中的分布式樣。在過去的幾年里，許多在種群水平上的重測序研究極大地豐富了人們對SVs機制和影響的認識。有研究表明，CNVs影響了大麥14個不同品種9.5%的編碼基因[22]。西藏自治區農牧科學院農業研究所青稞遺傳育種課題組之前的研究已對177份大麥材料進行了全基因組重測序，并進行了SNPs和INDELs分析，結果表明青稞起源于東方馴化大麥，最有可能在4 500年前至3 500年前，西藏青稞通過巴基斯坦北部-印度-尼泊爾傳入西藏東南部；由于高原極端生態條件導致的奠基者效應，使傳入后青稞遺傳多樣性快速降低，最終形成適應青藏高原生境的獨特的青稞類群[2]。但之前的研究只對177份大麥和青稞的群體樣本進行了SNP和INDELs分析，未對這些重測序數據進行基因組結構變異的分析。基于項目組之前已有的部分項目數據，結合從NCBI下載的5份數據，本研究分析182份大麥和青稞樣品的重測序數據，平均測序深度約為12 X。基于高覆蓋度基因組重測序數據，已經為大麥和青稞群體建立了高分辨率的基因組SV圖譜，包括DEL、INS、INV和ITX(染色體內易位)。這些數據將為進一步的基因組學和遺傳學研究，特別是探索相關基因的生物學功能起到重要作用。此外，作為一個有價值的補充，本試驗產出的數據將促進并豐富大麥屬物種的遺傳變異研究。

1 材料與方法

1.1 試驗材料

從存放在西藏自治區農牧科學院農業研究所的177份大麥/青稞樣品中提取基因組DNA，這些樣本來源于省部共建青稞和耗牛種質資源與遺傳改良國家重點實驗室青稞項目[2]，數據已上傳到NCBI數據庫(BioProject accession number:PRJNA417220)。另外5份重測序數據從美國國立生物技術信息中心(NCBI)下載。本研究共獲得182份大麥樣品的數據，包括29份來自西藏(27份)、以色列(1份)和俄羅斯(1份)的半野生或野生大麥樣品；從中國西藏、四川、青海、云南和浙江等省收集了70份地方品種的青稞；從日本(1份)、加拿大(3份)、澳大利亞(2份)、墨西哥(1份)、埃及(1份)、日本(1份)、敘利亞(2 份)和埃塞俄比亞(1份)共收集了12份地方品種的青稞；總共從中國西藏、青海、甘肅和四川收集了31份青稞育成品種，從加拿大收集了兩份育成品種；其余34份大麥樣品是在世界范圍內收集的。

1.2 DNA提取和測序

使用快速植物基因組DNA分離試劑盒(Rapid Plant Genomic DNA Isolation Kit，中國上海,生工生物技術公司)從100 mg幼葉組織中提取基因組DNA，檢測合格的DNA送往華大基因公司進行全基因組重測序。Illumina HiSeq 2000文庫的構建按照制造商的說明書進行(Illumina，圣地亞哥，加利福尼亞州，美國)。使用Illumina cluster station進行合成，試驗流程如下：模板雜交、等溫擴增、線性化、封閉、變性和測序引物雜交。使用標準的Illumina base calling流程將熒光圖像轉換成序列。DNA文庫插入片段為(454±39)bp，進行雙端測序。Reads長度從 90 bp到150 bp不等。

1.3 基因組結構變異分析

Breakdancer[23]和SVDetect[24]是結構變異檢測常用軟件,專門針對雙端測序數據進行開發，Breakdancer 和SVDetect可用于檢測DEL、INS、INV、ITX(Intra-chromosomal translocation)和染色體間易位(Inter-chromosomal translocation)，此外，SVDetect還可以用于檢測DUP和CNV。Breakdancer運行速度較快，而SVDetect由于步驟較多，運行速度較慢。為保證結果的準確性，本研究使用Breakdancer和SVDetect對182份大麥和青稞群體中的主要SVs類型進行鑒定，篩選共有的SVs類型進行下一步分析。大麥的參考基因組序列從IBSC(http://webblast.ipk-gatersleben.de/barley_ibsc/downloads/150831_barley_pseudomolecules.fasta.gz；https://webblast.ipk-gatersleben.de/barley_ibsc/downloads/160517_Hv_IBSC_PGSB_r1_CDS_HighConf_REPR_annotation.fasta.gz)下載。此外，為消除未定位到染色體上的序列對SVs檢測的影響，大麥參考基因組中chrUn的數據被排除在外。首先使用NGS QC Toolkit過濾原始數據并進行質控[25]。使用Bowtie 2軟件采用默認參數，將質控過后的序列比對到大麥參考基因組[26]。SAMtools用于將SAM格式轉換為BAM格式[27]。然后使用Picard包(https://broadinstitute.git hub.io/Picard/)去除重復序列，并使用Breakdancer[21]和SVDetect[22]軟件執行SVs分析。Breakdancer采用默認參數。SVDetect采用以下參數:(1)ration filter參數為 0.35，coverage filter參數為6、0.05和2，depth filter參數為1.5和1.5，SVN filter參數為2，其它參數為默認值。本研究將上述2個軟件的結果合并為最終的SVs結果。合并SVs條件如下:(1)多于4條支持序列或不少于3條不一致序列對；(2)SVs片段長度在100 bp和1 Mb之間[28]；(3)重疊區長度大于等于25 bp并且2個SV位于同一染色體并屬于同一變異類型，符合上述篩選條件的SVs被合并成一個SV[29]。本研究中只保留了Breakdancer和SVDetect軟件檢測到的主要SVs類型，包括DEL、INS、INV和ITX 4種SVs類型，并提取位于SVs變異區的基因用于后續功能分析。

1.4 基因功能分析

通過Gene Ontology Consortium(GO)所建立的數據庫和參考基因組注釋數據庫(http://webblast.ipk-gatersleben.de/barley_ibsc/)對位于SVs變異區的基因進行功能注釋和富集分析。

2 結果與分析

2.1 重測序數據和SVs在染色體上分布分析 結果

全基因重測序經數據質控后，總共產生 10 515 Gb數據。每個樣本的重測序數據平均為 57.78 Gb，序列平均覆蓋度為11.95X。SVs分析結果顯示共獲得了74 262個SVs，包括48 078個DEL(65%)、13 461個INSs(18%)、7 012個INV(9%)和5 711個ITXs(8%)(圖1)。每個樣本的平均SVs數量為408。在染色體分布式樣中，7H的SVs的數量最多，占24.39%，其次是6H(14.07%)、2H(14.01%)、3H(13.07%)、4H(12.34%)、5H(11.50%)和1H(10.59%)。此外，與其他染色體相比，7H含有最高密度(28個/Mb)的SVs(圖1)。

圖1 SVs 在大麥染色體上的分布

長度分布顯示，SVs片段長度主要分布在100 bp到1 kb之間，占所有SVs的64.26%。其中只有8.76%大于100 kb。大多數DELs(35 395/47 959)和INSs(11 466/13 180)長度從100 bp到1 kb不等。有趣的是，大多數INVs(3 139/7 012)和ITXs(2 252/5 708)的片段長度大于100 kb。與大麥參考基因組相比，在青稞的進化和選育過程中，很多基因組區域發生丟失，占所有SVs的64.93%(表1)。

表1 182份樣本的SVs 長度分布

進一步為獲得與SVs區域相關的可靠的基因集，從182份樣本中選擇了145份具有明確品種類型的樣品，將其分成3組，包括野生大麥(30個)、青稞改良品種(33個)和青稞地方品種(82個)。對每份樣品，采用以下條件下進行SVs區域相關基因的篩選:(1)基因區完全位于SV區；(2)每組80%以上樣本具有相同的基因被保留。基于大麥參考基因的注釋結果，分析結果顯示18.72%的基因(7 440/39 734)位于SVs結構變異區。與大麥參考基因組相比，在野生大麥、地方品種和改良品種中共有380、357和391個基因在進化或選育過程中發生丟失。

分析還發現大量基因插入到野生大麥群體中，占INSs的99.73%(1 093)。1個基因插入青稞地方品種中，2個基因插入到青稞育成品種中。INV和ITX相關基因的平均值分別為3 461和 2 257(圖2)。進一步分析發現，基因的差異顯示青稞育成品種、青稞地方品種和野生大麥中分別有62、23和37個獨特的與DELs相關基因,280個與DELs相關基因在3個群體中共有(圖2)。INSs分析結果表明，在青稞育成良品種、青稞地方品種和野生大麥群體中有2、1和1 093個獨特的基因，在青稞地方品種中獲得1個與INSs相關的基因。INVs結果表明，青稞育成品種、青稞地方品種和野生大麥分別有207、94和177個獨有基因。與INVs相關基因數量最多。染色體分布結果顯示，2H與DELs、INVs和ITX相關的基因數量最多。INSs主要分布在3H(158)和6H(936)染色體上，主要分布在野生大麥群體中，表明野生大麥含有較多的遺傳資源，是遺傳育種的重要基因資源庫(圖3)。

圖2 與SVs相關基因的韋恩圖

DEL.刪除;INS.插入;INV.倒位;ITX.染色體內易位

2.2 位于SVs結構變異區基因的功能富集分析結果

與SVs相關基因的功能分析結果顯示，篩選出的7 440個基因富集到945個GO term中，主要包括蛋白質結合(GO:0005515)、膜(GO:0016020)和氧化還原過程(GO:0055114)。GO功能富集最多的為INV(695)，其次是ITX(583)、INS(331)和DEL(174)(圖4)。此外，13個GO term與光系統相關，如光合作用(GO:0015979)、光周期開花調節(GO:2000028)和光合作用光反應(GO:0019684)。共有14個GO與植物應激反應過程有關。有趣的是，和大麥參考基因組相比，與真菌(GO:0050832)和細菌(GO:0042742)防御反應相關基因只在野生大麥中插入，與之相反，防御反應相關基因(GO:0006952)只在青稞育成品種發生丟失。在青稞育成品種中共有62個基因發生丟失，這些基因富集在48個GO term中，其中最豐富的是與蛋白質結合相關(44個基因)、鋅離子結合相關(16個基因)和催化活性相關(16個基因)。大多數與DELs相關的基因屬于大的多基因家族，包括糖基轉移酶家族、WRKY、GRAS、鋅指蛋白和五肽重復(PPR-like)家族。

圖4 與SVs相關基因的GO terms 注釋結果

3 討 論

不斷涌現的生物信息學工具和NGS測序技術的進步極大地促進了核苷酸變異的研究，盡管如此，對結構變異的研究在大多數物種中仍難以捉摸。本研究通過分析182份大麥和青稞群體全基因組重測序數據，用已公布的大麥基因組作為參考基因組，平均覆蓋度為12 X。共獲得74 262個結構變異，包含48 078個缺失(65%)、13 461個插入(18%)、7 012個倒位(9%)和5 711個染色體內易位(8%)。與SVs相關基因分析結果發現18.72%(7 440/39 734)的基因位于結構變異區域內。基因的功能分析發現，許多與結構變異相關的基因與光合作用和抗逆性相關，13個GO term與光系統相關，包括光合作用(GO:0015979)、光周期的開花調節(GO:2000028)和光合作用中的光反應(GO:0019684)。此外，還有14個GO term與植物應激反應過程有關。相關研究表明，與小尺度的核苷酸(<50 bp)變異相比，大片段的SVs在核苷酸變異中所占的比例更大，在人類基因組中，平均有8.9 Mb的堿基序列受到SVs的影響，而受到SNPs 影響的堿基序列只有3.6 Mb[30]。基于先前的研究結果，大麥基因組序列受到1.55 百萬個SNPs變異的影響[2]，其堿基序列長度遠小于本研究中SVs 影響的堿基序列。結合與SVs相關基因的功能富集分析結果，說明SVs在青稞適應外界環境中重要性，亦表明物種基因組的多樣性是其適應環境多樣性所必需的遺傳基礎。

測序技術和生物信息學工具雖然在近幾年取得了長足的進步，但核苷酸變異的研究目前主要集中在人類相關疾病的研究[30-31]。這種限制的主要原因是SVs是大規模的脫氧核糖核酸重組，帶來了計算和生物信息學的挑戰[32]。每種SV檢測方法都有不同的優勢和劣勢，這取決于SVs的類型或SVs位點相關序列的性質[33]。本研究中，使用2種不同工具來分析SVs：Breakdancer[21]和SVDetect[22]。這些工具使用1種或2～3種基于參考基因組的比對方法，比如測序深度、雙端比對或單端比對來檢測SVs。這些方法本身很可能不足以揭示所有SVs[34]。本研究通過設置嚴格的篩選參數，合并以上2種軟件的結果，以期望克服每種方法的弱點。

圖2分析發現，INS相關基因與其余三種類型的分布模式不同，與INS相關的基因在野生大麥中大量聚集，這一結果與圖3的結果一致，特別是在6號和3號染色體上。在野生大麥中發生插入結構變異，從另一個角度可以理解為在大麥參考基因組中發生大量丟失，我們推測在大麥參考基因組中，6號和3號染色體收到的選擇壓力遠遠強于其余5條染色體，這一現象也表明人工馴化或自然選擇會引起基因組發生大片段的丟失，導致片段上基因資源的丟失，亦說明保護物種野生種質資源的重要性[35]。對SVs相關基因進行進一步分析發現，野生大麥、青稞育成品種和青稞地方品種之間存在如下差異，與DELs相關的基因中，青稞育成品種和青稞地方品種共有20個相同的基因，其中62個基因為青稞育成品種特有，23個基因為青稞地方品種特有，在INV和ITX結構變異類型中，青稞育成品種具有類似的規律，其特有基因數量都大于青稞地方品種和野生大麥，可能與青稞育成品種受到強烈的人工馴化作用相關，對玉米和水稻等馴化歷史的研究顯示，人工馴化會導致遺傳多樣性顯著下降，即瓶頸效應[36-37]，導致優異基因資源的丟失。育種工作者可以利用這些基因開展相關研究，特別是開展野生大麥中與SVs相關優異基因的研究，促進優良青稞品種的選育工作。

綜上所述，本試驗首次通過大規模全基因組重測序數據，對大麥基因組中的SVs進行了研究。由于分析工具和計算資源的局限，只選擇了DEL、INS、INV和ITX 4種SVs類型，未能對大麥基因組中的CNVs進行分析，但本研究從各種類型結構變異中獲得了新的生物學認識，并建立一個結構變異的綜合數據集，為后續研究這些SVs相關基因對大麥和青稞生物多樣性的影響提供支撐，可能有助于進一步理解大麥適應外界環境和非生物脅迫的分子機制。