中國山楊群體CDS 遺傳多樣性格局與譜系地理分析

趙丹丹，張建國，惠 楠，陳為濤，王兆山＊

(1. 中國林業科學研究院林業研究所，國家林業局森林培育重點實驗室，北京 100091；2. 莒南縣自然保護地管理服務中心，山東 莒南 276600)

楊屬（Populus）物種是北半球十分重要的林木樹種及生態樹種， 中國山楊（Populus davidianaDode）是白楊派（Sect.Populus）的代表性物種之一，在我國山地森林中屬于非常重要的建群樹種。該樹種自我國東北地區至西南高山山區連續分布[1]，西南地區以亞熱帶熱帶季風氣候為主，冬季溫暖少雨，夏季高溫多雨，而北方-中部地區以溫帶季風氣候為主，冬天寒冷干燥，夏天炎熱多雨，不同的地理環境條件對中國山楊不同地區群體內遺傳多樣性的影響是巨大的。目前從表型性狀到分子水平，中國山楊物種內部不同自然群體之間已開展了遺傳多樣性和遺傳結構等譜系地理分析[2-5]，從局部以及整體分布區域的角度揭示了中國山楊的遺傳多樣性以及遺傳分化水平，然而針對中國山楊生長因子的分析鮮有研究。植物開花在植物生長甚至物種進化中處于核心地位，從植物開花基因的角度對中國山楊群體進行譜系地理分析，例如中國山楊遺傳多樣性及遷移路線，有助于為林木分子育種提供新的候選開花基因。

目前，學者們通過對植物成花調控關鍵因子的研究已總結出6 條調控途徑，即春化途徑、自主途徑、年齡途徑、光周期途徑、溫度途徑和赤霉素(GA)途徑，近年研究還發現2 條成花調控新途徑，即脫落酸(ABA)途徑和油菜素甾醇(BR)途徑，這8 條途徑共同構成了復雜的成花調控分子網絡，其中，FLC基因作為成花調控的關鍵因子，涉及多條途徑[6]。FLC基因是開花的主要抑制因子，正常情況下，植物通過環境低溫抑制FLC的表達[7]，促進開花。以往研究發現，REF6基因可以通過抑制FLC的表達，同樣可以促進植物開花[8]，突變型的BRR2a基因也會導致FLC的轉錄水平大大降低[9]。因此在低溫條件不足或低溫時長不夠的情況下REF6和BRR2a基因可以代替春化作用抑制FLC基因，從而促進植物適時開花。

本研究利用不同分布區的90 個中國山楊樣本的REF6和BRR2a基因的全長CDS 序列，從CDS 水平進行遺傳多態性、中性檢驗、單倍型網絡、系統發育分析，評估中國山楊不同群體內的遺傳多樣性，探討中國山楊在進化過程中的遷移路線，以及群體遺傳變異與環境變化的關系，解析相關變異的適生和進化意義，為后續研究中國山楊控制開花基因工作提供科學理論依據。

1 材料與方法

1.1 溫度數據獲取

根據中國山楊24 個采樣點經緯度信息，在worldClim 網站下載1970—2000 年氣候數據（ http://www.worldclim.org/）， 通過 DIVA-GIS（v7.3）軟件提取1970—2000 年間3—4 月份北方樣本、中部樣本和西南地區樣本點所對應經緯度的平均最高溫和最低溫（表1）。

表1 中國山楊采樣點信息及其1970—2000 年3—4 月平均氣溫最高及最低值Table 1 The information of sampling points of P. davidiana and the maximum and minimum average temperature values in the sample collection region from March to April, 1970—2000

1.2 基因CDS 序列數據獲取與分析

在 NCBI 網站（ https://www.ncbi.nlm.nih.gov/）獲取擬南芥REF6與BRR2a蛋白質序列（ AT1G20960.1 和 AT3G48430.1） ， 使 用blastp 軟件對毛果楊（Populus trichocarpaTorr. &Gray） 蛋白質序列（ Ptrichocarpa_444_v3.1.protein.fa）進行比對。分別選取比對結果最佳的毛果楊基因號（Potri.015G095500.1 和Potri.015G089000.2，附件1），在Ptrichocarpa_444_v3.1.gene.gff3 文件[10]中查找兩個基因的染色體位置及相對應的CDS 區間。基于本研究組另一項研究的中間產物，即中國山楊27X 重測序數據，通過使用bwa（v0.7.17）軟件[11]將其比對到毛果楊參考基因組（Ptrichocarpa_444_v3.0.fa）[10]獲得bam 文件，使用samtools（v1.8）軟件[12]和bcftools 軟件獲取每個個體基因組全長序列文件（.fa），基于Ptrichocarpa_444_v3.1.gene.gff3文件[10]提供的CDS 的位置信息，使用gffread 軟件從每個中國山楊個體序列文件（.fa）中提取兩個基因的CDS 序列，最后對所有個體相同的CDS 合并在同一個文件， 用于后續分析。1）使用DnaSP（v6.12.03）軟件[13]對每個CDS 序列的合并文件進行Phase 處理后，計算不同CDS 序列每個位點的分子多樣性指數，包括多態核苷酸分離位點數目（S）、序列間每堿基平均差異值（π）、多態核苷酸位點數目S 的Watterson’s θw 參數[14]、重組事件（Rm）、單倍型數量（Nh）、3 種中性檢驗Tajima’sD[15]、Fu and Li’sD、Fu and Li’sF[16]，從CDS 水平檢驗群體是否處于中性狀態，及群體在適應環境的過程中是否受到自然選擇作用；2）以毛果楊為外類群，在MEGA5（v5.05）軟件采用鄰域連接法，對每個基因的每個CDS 的所有個體分別構建DNA 單倍型的系統進化樹，結合不同群體的地理分布，分析不同群體之間的進化關系及順序；3）為了確保足夠的變異位點并且使單倍型網絡圖不過于復雜，在每個基因中選擇一個具有適當變異位點數量的CDS（REF6基因第1 個CDS，BRR2a基因第3 個CDS），在Network（v10.2）軟件采用中間連接法構建了單倍型網絡圖，并根據樣本分布的地理位置、樣本數量以及不同單倍型數量，在中國地圖上繪制每個采樣點包含的單倍型比例。

2 結果分析

2.1 CDS 水平分子多樣性指數

將擬南芥REF6與BRR2a基因AT 號對應的蛋白質序列與毛果楊蛋白質序列進行比對，REF6獲得了37 個基因號，相似度最高的為Potri.015G089000.1 和 Potri.015G089000.2， 均為62.629%（詳細數據未列出）。Ptrichocarpa_444_v3.1.gene.gff3 文件[10]顯示Potri.015G089000.2 與Potri.015G089000.1 具有相同的6 個CDS（表2 中REF6基因前6 個CDS），Potri.015G089 000.1 剩余1 個CDS 為11 142 487—11 144 891，Potri.015G089000.2 剩余2 個CDS 為11 142 487—11 144 538 和11 144 758—11 144 891，本研究選用Potri.015G089000.2（共8 個CDS）進行后續分析。BRR2a獲得了25 個基因號，相似度最高的為Potri.015G095500.1（81.810%）（詳細數據未列出），Ptrichocarpa_444_v3.1.gene.gff3 文件[10]顯示存在3 個CDS。不同CDS 序列的分子多樣性指數如表2 所示。REF6基因位于第15 號染色體，8 個CDS 平均序列長度為594.75 bp；BRR2a基因同樣位于第15 號染色體，3 個CDS平均2 185.00 bp。在北方與中部群體中，兩個基因CDS 水平的Tajima’sD，Fu and Li’sD和Fu and Li’sF檢驗平均值分別為-0.69～0.08、-0.05～1.08 和0.01～0.53；而西南群體兩個基因3 種中性檢驗的平均值均呈現顯著的負值（REF6基因存在2 個CDS 顯示p＜ 0.05，BRR2a基因存在2 個CDS 顯示p＜ 0.01）， 分別為-0.91 和-1.83、-0.26 和-1.08、-0.42 和-1.63，顯著低于北方與中部群體。核苷酸多態性顯示西南群體均值（REF6與BRR2a基因π 值分別為0.001 7、0.000 6，θw值分別為0.002 2、0.002 0）均低于北方和中部群體均值，重組率和單倍型也表現出同樣的趨勢。

表2 候選基因在CDS 水平的分子多樣性指數Table 2 Molecular diversity index of the candidate genes at the CDS level

2.2 CDS 水平系統發育樹

以毛果楊為外類群，針對90 個個體，對REF6的8 個CDS 和BRR2a的3 個CDS 構建了系統發育樹，從而評估每個CDS 的聚類程度和分支順序。如圖1 所示，不同的CDS 序列顯示出不同的樹狀拓撲結構，可以分為3 類。變異位點較多的系統發育樹顯示一部分北方和中部地區的個體均勻分布并與外類群相連，西南地區的個體主要聚集在末端，如圖1 中（A1）、（A6）和（A7），說明西南群體的CDS 存在進一步的分化；第二類為西南地區的個體與北方-中部地區的個體分為兩大枝，與外類群距離均等，如（B1）和（B3），或來自3 個地區的個體均隨機分布，如（A4）、（A5）和（B2），說明3 個群體的CDS 隨機分布，沒有明顯分化。第三類為變異位點少于5 且無明確規律的CDS。因此與北方和中部群體相比，西南群體存在較大分化，但3 個群體仍然共享基因庫。

圖1 中國山楊個體間REF6 與BRR2a 基因系統發育樹Fig. 1 Phylogenetic trees of the REF6 and BRR2a between individuals of P. davidiana

2.3 CDS 水平單倍型網絡圖

本研究選擇REF6基因第1 個CDS（Chr15：11 135 555—11 136 422）和BRR2a基因第3 個CDS（Chr15：11 629 168—11 632 128）分別構建了單倍型網絡圖，并在中國地圖上繪制每個采樣點包含的單倍型比例。對于REF6基因，在其CDS 序列（11 135 555—11 136 422）中發現了15 個單倍型（圖2A），3 種常見單倍型（頻率＞10%）為：H8（43.33%）、H11（15.56%）、H2（13.89%）。其中，北方群體存在單倍型H1～H10 共10 種，H8（21.67%）、H2（18.33%）、H3（16.67%）和H5（13.33%）為主要的單倍型。中部群體存在單倍型H2、H3、H5、H8、H9、H11、H12，其中H11（48.28%）、H2（24.14%）和H8（8.62%）占比較高。在西南群體中，H8（96.77%）為最常見的單倍型，H13 和H14 為西南群體特有的單倍型，僅各占1.61%。3 個群體共有的H8 單倍型在北方與中部群體占比較小，在西南群體占比非常大。在BRR2a基因的CDS 序列（11 629 168—11 632 128）中發現了28 個單倍型（圖2B），H22（27.78%）和H1（22.22%）是3 個群體中最常見的兩種單倍型。其中，北方群體存在單倍型H1～H12 共12 種，H1（31.67%）、H6（20.00%）和H5（10.00%）為主要的單倍型。中部群體存在單倍型H1、H2、H3、H8、H11～H21 共15 種， 其中H1（ 36.21%） 和H3（24.14%）較為常見。在西南群體中，存在H3、 H12、 H22～H27 共 8 種單倍型， H22（80.65%）為最常見的單倍型，H22～H27 均為西南群體特有的單倍型。兩種基因的單倍型在北方與中部多樣性均較為豐富，西南群體中單倍型種類相對單一，因此在西南地區的百分比相對較大。同時在西南群體中占比最大的H8 單倍型（REF6基因）和H22 單倍型（BRR2a基因）均分別與其外類群單倍型連接甚遠。

圖2 REF6 基因CDS（11 135 555-11 136 422）與BRR2a 基因CDS（11 629 168-11 632 128）單倍型網絡圖Fig. 2 Haplotype networks of the REF6 CDS (11 135 555-11 136 422) and BRR2a CDS (11 629 168-11 632 128)

3 討論

3.1 中國山楊遷移路線

在過去的三百萬年間，氣候振蕩已經嚴重地影響了北半球許多動植物物種的分布和遺傳結構[17-18]。目前已經有多項系統地理學的研究發現，中國東北部存在寒溫帶落葉喬木冰川避難所[19-20]。Hewitt 則指出另一種可能性，北方物種在冰期期間南遷，待冰期結束再回遷至高緯度地區；同時有一些物種的南方居群因無法適應南方高溫等因素而滅絕[21]。中國山楊在我國東北地區至西南高山山區連續分布，東北地區的群體可能是從南方遷移而來，也可能是東北地區的避難所幫助中國山楊在北方生存下來，冰期結束后逐漸向南遷移。基于BRR2a與REF6基因，本系統發育樹研究結果（圖1，A1，A6，A7）顯示外類群首先與北方-中部地區的中國山楊連接，西南群體往往傾向于聚集在末端，從CDS 水平支持了中國山楊從北方向南方遷移的結論[4]。同時，西南地區群體的核苷酸多態性指數π 與θw（表2）均低于北方和中部群體，單倍型網絡圖（圖2）也顯示外類群先與北方的單倍型相連，隨后與南方的單倍型相連，且北方和中部群體單倍型種類豐富，而西南群體相對來說種類單一，以上結果也均支持中國山楊從北向南遷移的觀點[4]。

3.2 群體遺傳變異與環境變化的關系

中國山楊耐寒冷、耐干旱瘠薄土壤，廣泛分布于東亞地區，分布區緯度跨度巨大，從寒溫帶延伸到亞熱帶，氣候差異很大。東亞西南部的中國山楊群體主要分布在中國四川省南部、貴州省西部和云南省，屬于南亞熱帶地區，而中國北方-中部地區的山楊群體主要分布在北亞熱帶、暖溫帶、溫帶和寒溫帶地區，因此西南地區的氣候特征，如海拔、溫度、濕度、光照、積溫等，與北方地區的氣候特征存在明顯差異，不同的地理環境條件對中國山楊個體的生長及種群間的遺傳結構存在一定的影響。另一方面，中國山楊在我國的花期為3—4 月[22-23]，西南地區與北方地區花期氣溫差異巨大（表1）。如果不同地區的中國山楊群體在同一溫度下開花，會造成北方花期太遲或者南方花期太早，導致其他氣候因子不匹配。因此一些學者認為，由于云貴高原氣候的獨特性，分布在西南地區的中國山楊分化嚴重[4,24]。本研究顯示在溫度較低的環境（北方和中部地區）單倍型豐富，各個單倍型比率相對較低，溫度較高的環境（西南地區）單倍型比率增高，例如西南群體中REF6基因H8 單倍型比例高達96.77%，BRR2a基因H22 單倍型比例高達80.65%，因此比例高的單倍型可能在適應西南地區的過程中受到了正選擇作用。同時REF6和BRR2a基因CDS 序列的各項遺傳多樣性指數顯示西南群體均明顯低于北方與中部群體，Tajima’sD，Fu and Li’sD和Fu and Li’sF中性檢驗指標在西南群體中均顯著偏離中性模型，呈顯著負值（表2），因此推測這兩個基因在西南群體分化過程中受到了強烈的自然選擇，與之前學者的研究結果一致[2,4-5]。這可能是因為在環境的選擇壓力下，中國山楊群體為了可以長期適應西南地區的環境特點，發生了局部適應。

4 結論

本研究REF6與BRR2a基因系統發育樹與單倍型網絡圖顯示外類群首先與中國山楊北方-中部群體連接，西南群體往往聚集在末端，西南群體遺傳多樣性水平低于北方與中部群體，支持了中國山楊從北方向南遷移的結論。REF6與BRR2a基因在西南地區均表現出正選擇的特征，推測西南地區的山楊群體在適應云貴高原獨特的氣候環境過程中發生了局部適應。未來需要更多的工作從基因克隆和表達分析等多個角度對控制植物開花的相關基因展開研究。