楊正釗 王梓豪 胡兆榮 辛明明 姚穎垠 彭惠茹 尤明山 宿振起 郭偉龍
小麥主栽品種濟麥22與良星99的基因組序列多態(tài)性比較分析
楊正釗 王梓豪 胡兆榮 辛明明 姚穎垠 彭惠茹 尤明山 宿振起*郭偉龍*
中國農(nóng)業(yè)大學農(nóng)學院/ 農(nóng)業(yè)生物技術(shù)國家重點實驗室/ 雜種優(yōu)勢研究與利用教育部重點實驗室, 北京 100193
濟麥22和良星99是我國黃淮冬麥區(qū)和北部冬麥區(qū)大面積推廣的高產(chǎn)小麥品種, 也是目前小麥雜交育種的重要親本。雖然濟麥22和良星99的來源和系譜不同, 但在重要農(nóng)藝、產(chǎn)量等性狀上存在較高的相似性。為了從全基因組水平研究其遺傳組成的異同, 本研究采用Illumina HiSeq2500測序平臺對上述兩個品種進行了全基因組測序(平均測序深度為5.8×), 并系統(tǒng)地比較了兩個品種拷貝數(shù)變異(CNV)、單核苷酸多態(tài)性(SNP)和插入/缺失(InDel)的序列差異。與中國春參考基因組序列相比, 兩個品種除了具有總長466 Mb的共有CNV變異區(qū)間外, 濟麥22和良星99的特有CNV變異區(qū)間的總長分別為91 Mb和45 Mb, 這些特有CNV區(qū)間主要集中在2B和4B染色體上; 濟麥22和良星99間存在1,547,371個SNP差異位點和137,817個InDel差異位點。以差異SNP分布規(guī)律為依據(jù), 在全基因組水平鑒定出濟麥22和良星99間存在14.2%的差異多態(tài)性熱點區(qū)間, 這些區(qū)間集中分布在1D、2B和4B染色體上。通過對5個控制小麥株高和穗長基因的序列分析, 發(fā)現(xiàn)有2個位于多態(tài)性熱點區(qū)間的基因在品種間存在移碼突變。本研究為利用重測序數(shù)據(jù)在基因組水平上比較小麥品種間遺傳差異提供了重要參考, 同時揭示了濟麥22和良星99在全基因組的遺傳相似區(qū)間和差異區(qū)間, 為今后小麥育種改良中更好利用濟麥22與良星99提供了重要遺傳信息。
小麥; 濟麥22; 良星99; 全基因組重測序; SNP; CNV
小麥是重要的口糧作物。我國是世界上最大的小麥生產(chǎn)和消費國, 提高單產(chǎn)是保障我國小麥生產(chǎn)可持續(xù)發(fā)展和糧食安全的重要途徑。新中國成立以來, 我國小麥主產(chǎn)區(qū)經(jīng)歷了5~6次的品種更換, 小麥品種的產(chǎn)量、抗性、品質(zhì)等重要性狀得到了大幅提高。小麥品種的遺傳改良對我國小麥單產(chǎn)和總產(chǎn)的提高發(fā)揮了重要作用[1]。不同小麥生態(tài)區(qū)中對照品種既是各時代育種水平的代表, 通常又是下一代品種改良的骨干親本[2]。例如, 黃淮冬麥區(qū)的對照品種周麥18[3]、魯麥14[4]、石4185[5]等在小麥的遺傳改良中做出重要貢獻。
常規(guī)育種是我國小麥遺傳改良的主要手段, 但在親本選配、后代選擇等關(guān)鍵育種過程中依賴于有限的表型性狀和育種者的經(jīng)驗。育種過程中存在不確定因素和盲目性, 影響了小麥遺傳改良的效率。分子標記的開發(fā)和應(yīng)用為小麥性狀的基因定位、標記輔助選擇、親本遺傳多樣性分析等方面提供了重要支撐, 可有效提升抗性、品質(zhì)等性狀遺傳的改良效率。繼水稻、玉米、大麥等主要作物基因組測序的完成, 六倍體普通小麥中國春的高質(zhì)量基因組參考基因組序列(IWGSCv1)的測序完成, 標志著小麥遺傳研究已進入“后基因組學時代”[6]。與基于表型和有限分子標記輔助的選擇相比, 利用基因組測序并從全基因組水平解析品種間的遺傳差異為小麥的遺傳研究和品種的改良提供高維度的組學數(shù)據(jù)支持[7]。
濟麥22是山東省農(nóng)業(yè)科學院作物研究所選育[8]。良星99是山東省德州市良星種子研究所選育[9]。濟麥22和良星99分別于2006年通過黃淮冬麥區(qū)北片審定, 此后相繼通過北部冬麥區(qū)和黃淮南片等麥區(qū)審定, 累計推廣面積超過6000萬公頃。由于綜合性狀突出, 濟麥22 (2015至今)和良星99 (2010—2015)先后作為黃淮冬麥區(qū)北片區(qū)域試驗的對照品種, 也是我國當前小麥品種改良的重要親本資源。育種家在長期育種和生產(chǎn)實踐中發(fā)現(xiàn), 雖然濟麥22和良星99來源及系譜不同, 但兩品種間除株高外, 田間總體形態(tài)相近; 主要農(nóng)藝、產(chǎn)量和抗病反應(yīng)等性狀表現(xiàn)一致, 兩者的雜交后代亦無明顯分離, 普遍認為兩個品種遺傳組成相似性較高, 親緣關(guān)系較近。近期研究發(fā)現(xiàn)麥22和良星99的2個衍生材料的白粉病抗性相近, 通過分子標記定位發(fā)現(xiàn)兩者具有相同基因[10]。但對其遺傳相似性的認識仍以有限的表型性狀描述和經(jīng)驗為主, 尚未從遺傳水平上對兩者進行系統(tǒng)解析。
本研究通過對濟麥22和良星99的重測序數(shù)據(jù)分析, 在全基因組水平系統(tǒng)解析了兩者的遺傳組成差異, 以期對今后小麥育種的親本利用、重要性狀基因區(qū)間追溯、新品種系譜分析、基因組水平上研究重要骨干種質(zhì)的遺傳組成等方面提供重要參考。
供試材料為小麥品種濟麥22 (系譜: 935024× 935106)與良星99 [系譜: (濟91102×魯麥14)×PH85-16]。2018年10月1日將上述兩個品種播種于中國農(nóng)業(yè)大學上莊實驗站(北京, 40°14′N, 116°19′E), 按照完全隨機試驗設(shè)計, 設(shè)3次重復, 2行區(qū), 行長1.5 m, 行距25.0 cm, 株距3.0 cm。收獲前每小區(qū)取中部10株調(diào)查株高、穗長、小穗數(shù)、穗粒數(shù)、旗葉長和寬等農(nóng)藝性狀; 收獲后利用考種儀(良田高拍儀S500A3B)考察千粒重、粒長、粒寬等籽粒性狀。品種間性狀差異檢驗用測驗進行統(tǒng)計分析, 數(shù)據(jù)分析通過R語言計算完成。
用CTAB的方法分別從濟麥22和良星99幼根提取DNA。用Illumina HiSeq2500測序平臺采用雙端測序的方法進行全基因組測序, 建庫和測序由北京諾禾致源科技股份有限公司完成。對重測序原始數(shù)據(jù)用Trimmomatic軟件(v0.36)[11]進行過濾后再使用BWA軟件(v0.7.15)[12]的BWA-MEM工具將過濾后的數(shù)據(jù)在小麥參考基因組(IWGSCv1)[13]上進行比對, 選擇保留存在“唯一最優(yōu)匹配(Unique Best Hit)”的讀段對; 最后利用與samtools (v1.4)[14]工具進行過濾。
利用GATK軟件(v3.8)[15]中的HaplotypeCaller、GenotypeGVCF、SelectVariants和VariantFiltration等功能計算并過濾獲得材料基因組上的高質(zhì)量SNP和InDel位點。其中, SNP位點過濾參數(shù)為“QD < 2.0, FS > 60.0, MQRankSum < –12.5, ReadPosRankSum < –8.0, SOR > 3.0, MQ < 40.0, DP > 30||DP < 3”; InDel過濾參數(shù)為“QD < 2.0, FS > 200.0, ReadPosRankSum < –20.0, DP > 30||DP < 3”。對過濾后的SNP與InDel位點通過SNPEff軟件(v4.3t)[16]進行注釋。使用的中國春基因注釋版本為IWGSC RefSeq v1.1 (https:// urgi.versailles.inra.fr/download/iwgsc/IWGSC_RefSeq_Annotations/v1.1/)。
拷貝數(shù)變異(copy-number variation, CNV)是小麥基因組研究的主要多態(tài)性類型之一。本研究為鑒定供試材料相對于中國春參照基因組的CNV變異區(qū)間, 在分析中以1 Mb為單位將全基因組劃分為小窗, 利用bedtools (v2.26.0)[17]計算兩個品種的重測序比對讀段在每個窗口內(nèi)的“平均覆蓋深度”(Depbin);并結(jié)合該材料的全基因組“平均讀段覆蓋深度”(Depave)進行歸一化, 得到每個小窗的“平均相對覆蓋深度”: Depbin/Depave。由于在CNV分析的結(jié)果中, 基因組序列丟失相比于序列插入更容易被檢測, 本研究計算得到的CNV變異區(qū)間僅考慮相對于參照基因組低覆蓋度(“丟失”)情況。下文中的CNV均指“丟失(低覆蓋度)變異”?!跋鄬ζ骄采w深度”低于0.5的小窗被視為“CNV變異區(qū)間”(附圖1)。CNV變異區(qū)間分析均由每個材料分別和“中國春”參照基因組的序列比較計算獲得。其中, 在兩個品種中同時存在CNV變異的區(qū)間被稱為“共有CNV變異區(qū)間”。
為鑒定兩個品種的基因組序列差異區(qū)間, 本研究以1 Mb為單位對全基因組進行劃分窗口, 對每個窗口中的單核苷酸多態(tài)性(single nucleotide polymorphism, SNP)的頻率進行分析, 在每個窗口中統(tǒng)計兩品種間存在差異的純合SNP的密度分布(附圖2)。這部分研究中僅考慮拷貝數(shù)正常的區(qū)間, CNV變異區(qū)間未被統(tǒng)計。考慮該密度分布具有明顯的高斯混合分布的特點, 我們利用EM算法擬合出兩個正態(tài)分布的模型, 并據(jù)此選擇閾值=1.5 (差異頻率約為3.1×10–5)作為2個分布的分隔邊界。其中, 高密度SNP分布區(qū)間被認為兩個材料的多態(tài)性熱點區(qū)間;低密度SNP分布區(qū)間被認為兩個材料的相似區(qū)間。
對濟麥22與良星99基因組內(nèi)的特有與共有CNV區(qū)間分別進行驗證, 選取區(qū)間內(nèi)2000 bp左右長度的序列設(shè)計染色體區(qū)間特異的引物對CNV區(qū)間進行驗證。擴增引物由在線軟件Primer 3.0設(shè)計, 其中引物長度范圍為18~24 bp, GC含量范圍在40%~60%之間, 退火溫度范圍為54~60℃, 擴增產(chǎn)物大小1~2 kb (表1)。以濟麥22、良星99與中國春的DNA為模板進行PCR擴增, 通過有無目標擴增產(chǎn)物判斷CNV是否存在。對與小麥株高相關(guān)的候選基因、的序列差異SNP位點, 設(shè)計出基因組特異引物(表1)進行PCR擴增并測序。PCR反應(yīng)體系為20 μL, 包括10 μL的2X M5 HiPer plusHiFi PCR mix, 正、反向引物(10 μmol L–1)各1 μL, 150 ng μL–1模板DNA 2 μL, 用ddH2O補至20 μL。PCR擴增程序為95℃ 3 min; 95℃30 s, 56~57.4℃ 30~60 s (依引物退火溫度以及目標序列而定), 72℃ 2 min, 35個循環(huán); 72℃5 min。
表1 本研究使用的引物編號和序列
田間觀察發(fā)現(xiàn)2個品種株形相近, 兩者除株高、穗長存在顯著差異外(< 0.01) (圖1), 旗葉大小、小穗數(shù)、千粒重、粒長、粒寬等主要農(nóng)藝和產(chǎn)量相關(guān)性狀差異均不顯著(表2)。這表明濟麥22和良星99在主要農(nóng)藝和產(chǎn)量性狀上存在較高的相似性。
利用Illumina HiSeq2500分別對濟麥22和良星99進行了全基因組測序, 分別獲得了88 Gb和128 Gb的重測序數(shù)據(jù)。通過和中國春參考基因組序列(IWGSCv1)進行測序數(shù)據(jù)的比對并去除重復序列, 選擇“唯一最優(yōu)匹配”的讀段進行后續(xù)分析, 最終兩品種的平均覆蓋深度均為5.8×。經(jīng)過SNP calling分析, 濟麥22和良星99中分別鑒定出15,326,314個與15,060,004個和中國春參照基因組不同的高質(zhì)量純合SNP位點, 約占全基因組的0.109%和0.107% (表3)。其中, A、B基因組中的SNP的頻率(0.130%~ 0.155%)明顯高于D基因組(0.0224%~0.0225%)。兩品種間基因型不一致的純合SNP位點僅占總SNP位點個數(shù)的9.8%, 即超過90%的SNP是兩個品種間共有的。同時, 在濟麥22和良星99中分別鑒定出1,143,512與1,118,564個的InDel位點(表4), 分別占全基因組的0.0813%和0.0795%。與D基因組, A、B基因組中鑒定出的InDel頻率更高。濟麥22和良星99間差異的InDel多態(tài)性位點占鑒定出的InDel位點總數(shù)的11.8%。
圖1 濟麥22與良星99的表型比較
A: 濟麥22 (左)與良星99 (右)的株高及節(jié)間長度比較; B: 濟麥22 (上)與良星99 (下)的粒型比較。
A: plant height and internode length comparison between Jimai 22 (left) and Liangxing 99 (right); B: grain morphological comparison between Jimai 22 (top) and Liangxing 99 (bottom).
表2 濟麥22與良星99的主要農(nóng)藝性狀差異顯著性分析
**表示在0.01水平上差異顯著。
**represents significantly different at the 0.01 probability level.
表3 濟麥22與良星99的純合SNP位點統(tǒng)計
表4 濟麥22與良星99的純合InDel位點統(tǒng)計
濟麥22的基因組相對于中國春參照基因組中存在557 Mb的CNV變異區(qū)間, 占全基因組的4.0%; 而在良星99基因組中檢測出511 Mb的CNV變異區(qū)間, 占全基因組的3.6% (表5)。其中, 兩品種共有的CNV變異區(qū)間的長度為466 Mb, 主要集中于2B、2D、5A和6A染色體; 兩品種特有的CNV區(qū)間分別為91 Mb和45 Mb, 共計約占總CNV區(qū)間的22.6%, 主要集中在2B和4B染色體上(圖2)。
為了驗證CNV區(qū)間的分析結(jié)果, 在2B和2D染色體上的CNV差異區(qū)間中設(shè)計了特異PCR引物。其中2B_LX99_CNV1與2B_LX99_CNV2位于2B染色體上良星99特有CNV區(qū)間(Chr.2B: 667 Mb~ 726 Mb), 2B_COM_CNV位于2B染色體上良星99與濟麥22共有的CNV區(qū)間(Chr. 2B: 738 Mb~769 Mb), 2D_COM_CNV位于2D染色體上濟麥22與良星99共有的CNV區(qū)間(Chr. 2D: 570 Mb~638 Mb)。利用這4對引物在濟麥22、良星99和中國春材料內(nèi)進行PCR擴增, 結(jié)果表明材料間的差異CNV區(qū)間可以通過擴增條帶的有無進行驗證(圖3); 其中2B_ LX99_CNV1與2B_LX99_CNV2在中國春與濟麥22中均擴增出目標條帶, 但未在良星99中擴增出條帶, 證明該區(qū)間僅在良星99中丟失, 2B_ COM_CNV與2D_COM_CNV僅在中國春中擴增出目標條帶, 證明該區(qū)間在濟麥22與良星99中均丟失。PCR驗證結(jié)果與全基因組分析結(jié)果吻合, 表明該分析方法計算得到的CNV變異準確可靠。
表5 濟麥22與良星99的全基因組CNV區(qū)間長度統(tǒng)計
圖2 濟麥22與良星99的CNV區(qū)間(相對于中國春參照基因組)在各染色體上的分布
綠色: 濟麥22的特有CNV區(qū)間; 粉色: 良星99的特有CNV區(qū)間; 紫色: 濟麥22與良星99的共有CNV區(qū)間。
Green: specific CNV regions in Jimai 22; Pink: specific CNV regions in Liangxing 99; Purple: common CNV regions in the two varieties.
圖3 利用PCR對濟麥22(JM22)、良星99(LX99)和中國春(CS)中檢測的CNV區(qū)間進行驗證
2D_COM_CNV和2B_COM_CV分別為2D和2B染色體上共有區(qū)間序列設(shè)計的引物; 2B_LX99_CNV1和2B_LX99_CNV2為在2B染色體上良星99特有的CNV區(qū)間中序列設(shè)計的引物。每對引物對應(yīng)的3個泳道從左到右分別為: 良星99(LX99)、濟麥22(JM22)和中國春(CS)。
2D_COM_CNV and 2B_COM_CNV were designed utilizing sequences in common CNV regions in 2D and 2B chromosomes, respectively; 2B_LX99_CNV1 and 2B_ LX99_CNV2 were designed utilizing sequences in Liangxing 99-specific CNV regions in chromosome 2B. The corresponding three lanes from left to right for each primer are Liangxing 99 (LX99), Jimai 22 (JM 22) and Chinese Spring (CS).
通過對濟麥22與良星99中特有的CNV區(qū)間內(nèi)高可信度(HC)基因進行GO富集分析發(fā)現(xiàn), 在濟麥22的特有CNV區(qū)間內(nèi)有982個基因, 其中294個基因顯著富集在25個GO條目, 主要包括“ADP 結(jié)合”(GO: 0043531)、“肽酶抑制劑活性”(GO: 0030414),輔酶結(jié)合”(GO: 0050662)、“酶抑制劑活性”(GO: 0004857)、“半胱氨酸型內(nèi)肽酶抑制劑活性”(GO: 0004869)等分子功能, 以及“肽酶活性負調(diào)控”(GO: 0010466)、“內(nèi)肽酶活性負調(diào)控”(GO: 0010951)、“催化活性負調(diào)控”(GO: 0043086)、“碳代謝過程”(GO: 0008152)等生物學過程中, 細胞組分分類內(nèi)未有富集(圖4-A)。在良星99特有CNV區(qū)間內(nèi)包括513個基因, 其中富集到的395個基因分布在25個GO條目中, 主要有“ADP結(jié)合”(GO: 0043531)、“血紅素結(jié)合”(GO: 0020037)、“過氧化物酶活性”(GO: 0004601)、“萜烯合酶活性”(GO: 0010333)、“O-甲基轉(zhuǎn)移酶活性”(GO: 0008171)、“腺苷脫氨酶活性”(GO: 0004000)等分子功能, 以及“抗氧化反應(yīng)”(GO: 0006979)、“過氧化氫分解過程”(GO: 0042744)、“細胞氧化解毒”(GO: 0098869)、“蛋白磷酸化”(GO: 0006468)、“糖原生物合成”(GO: 0005978)等生物學過程和“胞外區(qū)”(GO:0005576)這一細胞組分中(圖4-B)。濟麥22特有CNV區(qū)間內(nèi)富集在GO條目, 如“生長素反應(yīng)”(GO: 0009733), “色素合成過程”(GO: 0046148)下的基因, 可能是影響兩材料間的表型差異的候選基因。
圖4 濟麥22(A)與良星99(B)特有CNV區(qū)間中基因的GO富集分析結(jié)果
綠色: 生物過程; 紫色: 細胞組分; 橙色: 分子功能。
以中國春IWGSCv1參考基因組為模板, 利用剔除品種間大片段CNV后的13.5 G大小的基因組序列進行分析。以差異純合SNP的密度為3.1×10–5為閾值將染色體區(qū)分為“多態(tài)性熱點區(qū)間”(高密度差異SNP區(qū)間)和“序列相似區(qū)間”(低密度差異SNP區(qū)間)。濟麥22與良星99之間共有1915 Mb (14.2%)區(qū)間為多態(tài)性熱點區(qū)間(圖5-A)。在多態(tài)性熱點區(qū)間中, 兩品種間有1,375,981個純合差異SNP位點, 約占兩品種全部純合差異SNP位點總數(shù)的96.5% (圖5-B), 這些多態(tài)性熱點區(qū)間主要集中于1D、2B、3B和4B染色體(圖5-C)。
通過對兩品種多態(tài)性熱點區(qū)間內(nèi)的14,306個高可信度基因進行GO富集分析, 其中1461個基因顯著富集在20個GO條目中, 主要包括“5'–3' RNA聚合酶活性” (GO: 0003899)、“二磷酸核酮糖羧化酶活性” (GO: 0016984)、“光合作用的循環(huán)電子轉(zhuǎn)運通路內(nèi)轉(zhuǎn)移電子的電子轉(zhuǎn)運器活性” (GO: 0045156)、“醛糖1-表異構(gòu)酶活性” (GO: 0004034)、以NAD或NADP為受體提供CH-OH基團的氧化還原酶活性” (GO: 0016616)等分子功能, “光合作用” (GO: 0015979)、“己糖代謝過程” (GO: 0019318)和“光系統(tǒng)II中的光合電子傳輸” (GO: 0009772)等生物過程, 以及“光系統(tǒng)II” (GO: 0009523)、“光系統(tǒng)II反應(yīng)中心” (GO: 0009539)、“葉綠體” (GO: 0009507)、“質(zhì)體” (GO: 0009536)和“類囊體” (GO: 0007579)等細胞組分內(nèi)(圖6)。綜合看來, 該部分基因的GO富集分析結(jié)果主要集中在光合作用通路相關(guān)條目。
圖5 濟麥22與良星99全基因組內(nèi)多態(tài)性熱點區(qū)間分布
A: 兩品種間差異SNP位點中分布在多態(tài)性熱點區(qū)間與序列相似區(qū)間的比例; B: 多態(tài)性熱點區(qū)間與序列相似區(qū)間的長度在全基因組上的比例; C: 多態(tài)性熱點區(qū)間與序列相似區(qū)間在全基因組分布。深藍色: 多態(tài)性熱點區(qū)間; 淺藍色: 序列相似區(qū)間; 白色: CNV區(qū)間。
A: pie chart for the proportions of differential SNPs in the polymorphism hotspot regions (PHRs) and the genetic similar regions (GSRs); B: pie chart for the proportions of the PHRs and GSRs in length; C: distributions of the PHRs and the GSRs across chromosomes.
圖6 濟麥22與良星99間多態(tài)性熱點區(qū)間的基因GO富集分析結(jié)果
綠色: 生物過程; 紫色: 細胞組分; 橙色: 分子功能。
為比較兩品種間多態(tài)性熱點區(qū)間和序列相似區(qū)間中的突變類型的異同, 我們首先統(tǒng)計差異SNP位點的堿基變異類型的頻率(圖7-A和表6)。無論在多態(tài)性熱點區(qū)間還是序列相似區(qū)間, T?C和G?A類型的突變頻率均為最多, 在多態(tài)性熱點區(qū)間內(nèi)約占71.4%, 略高于在序列相似區(qū)間內(nèi)所占總SNP的比例(66.5%)。同時, 我們還統(tǒng)計了兩品種間差異InDel的長度分布(圖7-B和表7)。多態(tài)性熱點區(qū)間中的單堿基的InDel占該區(qū)間內(nèi)所有InDel位點的64.7%, 其比例高于單堿基的InDel在序列相似區(qū)間中所占的比例(49.0%)。
圖7 濟麥22與良星99間多態(tài)性熱點區(qū)間與序列相似區(qū)間中各自的InDel長度分布(A)和SNP突變類型(B)
深藍: 多態(tài)性熱點區(qū)間; 淺藍: 序列相似區(qū)間。
Dark blue: polymorphism hotspot regions; Light blue: genetic similar regions.
我們分別對多態(tài)性熱點區(qū)間與序列相似區(qū)間內(nèi)兩品種間差異SNP與InDel位點對于基因功能的影響進行了注釋與比較(表8), 結(jié)果表明多態(tài)性熱點區(qū)間與序列相似區(qū)間中基因序列上的SNP和InDel變異的類型主要以錯義突變(missense variation)和同義突變(synonymous variation)為主。在多態(tài)性熱點區(qū)間中錯義突變(51.3%)與終止子相關(guān)突變(stop-codon gain/lost)(1.5%)所占比例要分別略高于兩者在序列相似區(qū)間的44.7%和0.8%比例。
表6 濟麥22與良星99全基因組差異純合SNP的突變類型的統(tǒng)計
表7 濟麥22與良星99的純合InDel變異的長度分布
表8 濟麥22與良星99的純合點突變對蛋白編碼功能影響的統(tǒng)計
濟麥22與良星99在株高和穗長上存在顯著差異, 利用兩者的重測序數(shù)據(jù)對小麥目前已克隆的矮稈基因[18]、[18]、株高相關(guān)基因[19]和在2D染色體的內(nèi)定位到2個株高穗長候選基因[20]等進行序列分析, 發(fā)現(xiàn)區(qū)間內(nèi)的2個基因在品種間存在序列差異(表9)。
表9 小麥株高相關(guān)基因列表及在序列相似區(qū)間和多態(tài)性熱點區(qū)間的分布
利用基因特異性引物對QPht/Sl.cau-2D.1 QTL區(qū)間內(nèi)的和基因分別在濟麥22與良星99中進行擴增、測序發(fā)現(xiàn), 2個基因在品種間存在引起氨基酸改變的差異位點(圖8)。該QTL區(qū)間是矮稈基因所在區(qū)間, 該位點對小麥的株高和穗長均有影響[21], 但其遺傳功能仍需進一步試驗驗證。
圖8 TraesCS2D02G051500(A)與TraesCS2D02G055700(B)在中國春、濟麥22和良星99中功能改變位點附近的核苷酸序列
本研究通過對濟麥22和良星99兩個小麥品種的重測序數(shù)據(jù)進行分析, 系統(tǒng)研究了兩個品種在基因組水平上存在的CNV區(qū)間和SNP熱點區(qū)間。與中國春參照基因組比較發(fā)現(xiàn), 兩個品種均有大量CNV丟失變異區(qū)間, 約占全基因組的4.3%。除集中分布在2B染色體上148 Mb的CNV變異外, 大部分CNV變異區(qū)間零散地分布在不同染色體上。這一結(jié)果表明, 與中國春參考基因組序列相比, 濟麥22和良星99中存在大量的CNV變異(特別是CNV丟失變異)。通過比較2個材料中CNV變異區(qū)間的異同, 發(fā)現(xiàn)大部分CNV區(qū)間是共有的, 但仍有22.6%的差異CNV區(qū)間; 其中60.3%的差異CNV區(qū)間集中于2B染色體。需要指出的是, 本研究中發(fā)現(xiàn)的“CNV缺失區(qū)間”是通過讀段比對后的覆蓋度確定的, 不能排除該區(qū)域可能來源于遺傳距離較遠的基因組序列的可能性。根據(jù)品種間差異SNP的分布密度對染色體區(qū)間進行了劃分, 鑒定出基因組14.2%的區(qū)間是兩品種間的多態(tài)性熱點區(qū)間, 集中了全基因組96.5%差異SNP位點, 這些區(qū)間主要位于1D、2B和4B染色體上。通過鑒定相似區(qū)間和多態(tài)性熱點區(qū)間變異位點的特征, 發(fā)現(xiàn)兩類區(qū)間在SNP變異類型、InDel長度分布、突變影響編碼功能的比例上有一定區(qū)別。上述結(jié)果表明濟麥22和良星99在全基因組水平具有較高的相似性, 但在特定的染色體上確實存大片段的差異區(qū)間。
濟麥22和良星99的主要農(nóng)藝和產(chǎn)量性狀存在較高相似性, 只有在株高和穗長上存在差異。結(jié)合品種間序列差異位點分析, 我們推測控制品種間株高和穗長的差異基因可能位于差異SNP熱點區(qū)域或差異CNV區(qū)間。結(jié)合目前已經(jīng)在小麥中克隆的株高或穗長的5個相關(guān)基因, 并分析兩品種間基因序列, 發(fā)現(xiàn)了2個位于差異SNP熱點區(qū)且存在影響編碼功能的變異的候選基因。本研究在已有基因組變異和表型信息的基礎(chǔ)上對已知性狀相關(guān)基因和品種間差異序列區(qū)間進行了初步分析, 找到了位于多態(tài)性熱點區(qū)間的2個具有序列差異的基因, 但其生物學功能和關(guān)聯(lián)仍需進一步驗證。利用表型差異和全基因組序列多態(tài)性分析相結(jié)合的方法, 本研究的方法可為今后借助基因組測序方法快速定位和克隆候選基因提供了新的參考。
Zhao等[22]研究發(fā)現(xiàn), 良星99的2BL上有一個控制抗白粉病的主效基因, 該基因被命名為[23]。利用中國春參照基因組的信息, 將基因定位在chr2B的581 Mb~596 Mb的區(qū)間[24]。本研究的結(jié)果顯示, 濟麥22和良星99的2BL的絕大部分區(qū)間為多態(tài)性熱點區(qū)間, 但575 Mb~603 Mb區(qū)段為本研究中鑒定的“序列相似區(qū)間”, 該區(qū)段包含所在的物理區(qū)間; 本研究間接支持了Qu等[24]在濟麥22和良星99的衍生材料中均攜帶相同的基因的結(jié)論。該結(jié)果也表明本研究鑒定的序列相似區(qū)間可為抗病等關(guān)鍵性狀基因的定位研究提供參考。
小麥品種的重測序研究為探索小麥基因組序列的多態(tài)性提供高維度、高精細度、類型豐富的基因組變異信息。本研究也為比較小麥品種間的基因組序列差異分析提供了方法學參考。隨著小麥材料重測序數(shù)據(jù)的不斷積累, 小麥品種間基因組序列變異的特征和規(guī)律也將得到進一步揭示。有效分析和利用豐富的基因組變異信息, 將有助于對關(guān)鍵性狀基因功能進行闡釋, 輔助小麥育種對優(yōu)異基因區(qū)間的利用, 提高小麥品種選育的效率。
濟麥22和良星99兩個重要小麥品種的田間表型相似。本研究通過對兩個品種進行全基因序列差異分析發(fā)現(xiàn), 濟麥22和良星99間的遺傳組成相似性較高(序列相似區(qū)間占85.8%), 但仍存在大區(qū)段CNV變異(136 M)和占全基因組14.2%的差異多態(tài)性熱點區(qū)間; 其中差異SNP熱點區(qū)間主要集中在1D、2B和4B染色體上。本工作為進一步研究和利用濟麥22和良星99提供了重要的基因組變異數(shù)據(jù); 也為小麥品種間的基因組變異熱點區(qū)間的鑒定提供了方法參考。
說明: 本研究中的原始測序數(shù)據(jù)已提交至中國科學院北京基因組研究所BIG數(shù)據(jù)中心[25]的組學原始數(shù)據(jù)歸檔庫[26]中(https://bigd.big.ac.cn/gsa), 編號為CRA002333。為便于數(shù)據(jù)的共享, 我們利用SnpHub數(shù)據(jù)庫模型[27]搭建了2個材料的基因組變異信息查詢數(shù)據(jù)庫(http://wheat.cau.edu.cn/Wheat_SnpHub_ Portal/lx99_jm22/)。
[1] 何中虎, 莊巧生, 程順和, 于振文, 趙振東, 劉旭.中國小麥產(chǎn)業(yè)發(fā)展與科技進步. 農(nóng)學學報, 2018, 8(1): 99–106. He Z H, Zhuang Q S, Cheng S H, Yu Z W, Zhao Z D, Liu X. Wheat production and technology improvement in China., 2018, 8(1): 99–106 (in Chinese).
[2] 張婷, 逯臘虎, 張偉, 楊斌, 袁凱, 史曉芳. 黃淮麥區(qū)300個小麥種質(zhì)的主要農(nóng)藝性狀比較.小麥研究, 2018, 39(1): 13–20. Zhang T, Lu L H, Zhang W, Yang B, Yuan K, Shi X F. Comparison of main agronomic traits of 300 wheat germplasms in Huanghuai Wheat Region., 2018, 39(1): 13–20 (in Chinese).
[3] 曹廷杰, 王西成, 趙虹. 國審小麥新品種周麥18號豐產(chǎn)性、穩(wěn)產(chǎn)性及適應(yīng)性分析. 中國農(nóng)業(yè)科技導報, 2007, 9(1): 39–41. Cao T J, Wang X C, Zhao H. Analysis on the yield potential, stability and adaptability of national authorized wheat variety Zhoumai 18., 2007, 9(1): 39–41 (in Chinese with English abstract)
[4] 蓋紅梅, 李玉剛, 王瑞英, 李振清, 王圣健, 高峻嶺, 張學勇. 魯麥14對山東新選育小麥品種的遺傳貢獻. 作物學報, 2012, 38: 954–961. Gai H M, Li Y G, Wang R Y, Li Z Q, Wang S J, Gao J L, Zhang X Y. Genetic contribution of Lumai 14 to novel wheat varieties developed in Shandong province., 2012, 38: 954–961 (in Chinese with English abstract)
[5] 傅曉藝, 郭進考, 劉艷濱, 史占良, 單子龍, 韓然, 何明琦. 利用石4185為骨干親本培育高產(chǎn)小麥新品種. 種子, 2017, 36(12): 95–99. Fu X Y, Guo J K, Liu Y B, Shi Z L, Shan Z L, Han R, He M Q. New wheat varieties with high yield of core parent breeding by using Shi 4185., 2017, 36(12): 95–99 (in Chinese).
[6] 凌宏清. 小麥及其近緣種基因組測序研究進展與發(fā)展趨勢. 麥類作物學報, 2016, 36: 397–403. Ling H Q. Progress and perspectives of the genome sequencing in wheat and its relatives., 2016, 36: 397–403 (in Chinese with English abstract)
[7] Uauy C. Wheat genomics comes of age., 2017, 36: 142–148.
[8] 殷貴鴻, 李根英, 何中虎, 劉建軍, 王輝, 夏先春. 小麥新品種濟麥22抗白粉病基因的分子標記定位. 作物學報, 2009, 35: 1425–1431. Yin H G, Li G Y, He Z H, Liu J J, Wang H, Xia X C. Molecular mapping of powdery mildew resistance gene in wheat cultivar Jimai 22., 2009, 35: 1425–1431 (in Chinese with English abstract)
[9] 趙廣軍, 張英肖. 高產(chǎn)穩(wěn)產(chǎn)冬小麥品種良星99栽培技術(shù). 小麥研究, 2006, 27(3): 12–13. Zhao G J, Zhang Y X. Cultivation technique of Liangxing 99, a winter wheat variety with high and stable yield., 2006, 27(3): 12–13 (in Chinese)
[10] Qu Y F, Wu P P, Hu J H, Chen Y X, Shi Z L, Qiu D, Li Y H, Zhang H J, Zhou Y, Yang L, Liu H W, Zhu T Q, Liu Z Y, Zhang Y M, Li H J. Molecular detection of the powdery mildew resistance genes in winter wheats DH51302 and Shimai 26., 2020, 19: 931–940.
[11] Bolger A M, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data., 2014, 30: 2114–2120.
[12] Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform., 2009, 25: 1754–1760.
[13] Appels R, Eversole K, Feuillet C, Keller B, Rogers J, Stein N, Pozniak C J, Stein N, Choulet F, Distelfeld A, Eversole K, Poland J, Rogers J, Ronen G, Sharpe A G, Pozniak C, Ronen G, Stein N, Barad O, Baruch K, Choulet F, Keeble-Gagnere G, Mascher M, Sharpe A G, Ben-Zvi G, Josselin A A, Stein N, Mascher M, Himmelbach A, Choulet F, Keeble-Gagnere G, Mascher M, Rogers J, Balfourier F, Gutierrez-Gonzalez J, Hayden M, Josselin A A, Koh C, Muehlbauer G, Pasam R K, Paux E, Pozniak C J, Rigault P, Sharpe A G, Tibbits J, Tiwari V, Choulet F, Keeble-Gagnere G, Mascher M, Josselin A A, Rogers J, Spannagl M, Choulet F, Lang D, Gundlach H, Haberer G, Keeble-Gagnere G, Mayer KFX, Ormanbekova D, Paux E, Prade V, Simkova H, Wicker T, Choulet F, Spannagl M, Swarbreck D, Rimbert H, Felder M, Guilhot N, Gundlach H, Haberer G, Kaithakottil G, Keilwagen J, Lang D, Leroy P, Lux T, Mayer KFX, Twardziok S, Venturini L, Appels R, Rimbert H, Choulet F, Juhasz A, Keeble-Gagnere G, Choulet F, Spannagl M, Lang D, Abrouk M, Haberer G, Keeble-Gagnere G, Mayer KFX, Wicker T, Choulet F, Wicker T, Gundlach H, Lang D, Spannagl M, Lang D, Spannagl M, Appels R, Fischer I, Uauy C, Borrill P, Ramirez-Gonzalez R H, Appels R, Arnaud D, Chalabi S, Chalhoub B, Choulet F, Cory A, Datla R, Davey M W, Hayden M, Jacobs J, Lang D, Robinson S J, Spannagl M, Steuernagel B, Tibbits J, Tiwari V, van Ex F, Wulff BBH, Pozniak C J, Robinson S J, Sharpe A G, Cory A, Benhamed M, Paux E, Bendahmane A, Concia L, Latrasse D, Rogers J, Jacobs J, Alaux M, Appels R, Bartos J, Bellec A, Berges H, Dolezel J, Feuillet C, Frenkel Z, Gill B, Korol A, Letellier T, Olsen O A, Simkova H, Singh K, Valarik M, van der Vossen E, Vautrin S, Weining S, Korol A, Frenkel Z, Fahima T, Glikson V, Raats D, Rogers J, Tiwari V, Gill B, Paux E, Poland J, Dolezel J, Cihalikova J, Simkova H, Toegelova H, Vrana J, Sourdille P, Darrier B, Appels R, Spannagl M, Lang D, Fischer I, Ormanbekova D, Prade V, Barabaschi D, Cattivelli L, Hernandez P, Galvez S, Budak H, Steuernagel B, Jones J D G, Witek K, Wulff B B H, Yu G, Small I, Melonek J, Zhou R, Juhasz A, Belova T, Appels R, Olsen O A, Kanyuka K, King R, Nilsen K, Walkowiak S, Pozniak C J, Cuthbert R, Datla R, Knox R, Wiebe K, Xiang D, Rohde A, Golds T, Dolezel J, Cizkova J, Tibbits J, Budak H, Akpinar B A, Biyiklioglu S, Muehlbauer G, Poland J, Gao L, Gutierrez-Gonzalez J, N'Daiye A, Dolezel J, Simkova H, Cihalikova J, Kubalakova M, Safar J, Vrana J, Berges H, Bellec A, Vautrin S, Alaux M, Alfama F, Adam-Blondon A F, Flores R, Guerche C, Letellier T, Loaec M, Quesneville H, Pozniak C J, Sharpe A G, Walkowiak S, Budak H, Condie J, Ens J, Koh C, Maclachlan R, Tan Y, Wicker T, Choulet F, Paux E, Alberti A, Aury J M, Balfourier F, Barbe V, Couloux A, Cruaud C, Labadie K, Mangenot S, Wincker P, Gill B, Kaur G, Luo M, Sehgal S, Singh K, Chhuneja P, Gupta O P, Jindal S, Kaur P, Malik P, Sharma P, Yadav B, Singh N K, Khurana J, Chaudhary C, Khurana P, Kumar V, Mahato A, Mathur S, Sevanthi A, Sharma N, Tomar R S, Rogers J, Jacobs J, Alaux M, Bellec A, Berges H, Dolezel J, Feuillet C, Frenkel Z, Gill B, Korol A, van der Vossen E, Vautrin S, Gill B, Kaur G, Luo M, Sehgal S, Bartos J, Holusova K, Plihal O, Clark M D, Heavens D, Kettleborough G, Wright J, Valarik M, Abrouk M, Balcarkova B, Holusova K, Hu Y, Luo M, Salina E, Ravin N, Skryabin K, Beletsky A, Kadnikov V, Mardanov A, Nesterov M, Rakitin A, Sergeeva E, Handa H, Kanamori H, Katagiri S, Kobayashi F, Nasuda S, Tanaka T, Wu J, Appels R, Hayden M, Keeble-Gagnere G, Rigault P, Tibbits J, Olsen O A, Belova T, Cattonaro F, Jiumeng M, Kugler K, Mayer K F X, Pfeifer M, Sandve S, Xun X, Zhan B, Simkova H, Abrouk M, Batley J, Bayer P E, Edwards D, Hayashi S, Toegelova H, Tulpova Z, Visendi P, Weining S, Cui L, Du X, Feng K, Nie X, Tong W, Wang L, Borrill P, Gundlach H, Galvez S, Kaithakottil G, Lang D, Lux T, Mascher M, Ormanbekova D, Prade V, Ramirez-Gonzalez R H, Spannagl M, Stein N, Uauy C, Venturini L, Stein N, Appels R, Eversole K, Rogers J, Borrill P, Cattivelli L, Choulet F, Hernandez P, Kanyuka K, Lang D, Mascher M, Nilsen K, Paux E, Pozniak C J, Ramirez-Gonzalez R H, Simkova H, Small I, Spannagl M, Swarbreck D, Uauy C. Shifting the limits in wheat research and breeding using a fully annotated reference genome., 2018, 361: eaar7191.
[14] Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R. The sequence alignment/map format and SAMtools., 2009, 25: 2078–2079.
[15] van der Auwera G A, Carneiro M O, Hartl C, Poplin R, Del Angel G, Levy-Moonshine A, Jordan T, Shakir K, Roazen D, Thibault J, Banks E, Garimella K V, Altshuler D, Gabriel S, Depristo M A. From FastQ data to high-confidence variant calls: The genome analysis toolkit best practices pipeline., 2013, 43: 10–11.
[16] Cingolani P, Platts A, Wang L L, Coon M, Nguyen T, Wang L, Land S J, Lu X, Ruden D M. A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff., 2014, 6: 80–92.
[17] Quinlan A R, Hall I M. BEDTools: a flexible suite of utilities for comparing genomic features., 2010, 26: 841–842.
[18] Peng J, Richards D E, Hartley N M, Murphy G P, Devos K M, Flintham J E, Beales J, Fish L J, Worland A J, Pelica F, Sudhakar D, Christou P, Snape J W, Gale M D, Harberd N P. ‘Green revolution’ genes encode mutant gibberellin response modulators., 1999, 400: 256–261.
[19] Beales J, Turner A, Griffiths S, Snape J W, Laurie D A. A pseudo-response regulator is misexpressed in the photoperiod insensitivemutant of wheat (L.)., 2007, 115: 721–733.
[20] 柴嶺嶺. 普通小麥2D染色體株高和穗長QTL的精細定位與圖位克隆. 中國農(nóng)業(yè)大學博士學位論文, 北京, 2019. pp 38–39. Chai L L. Fine Mapping and Map-based Cloning of QTL Controlling Plant Height and Spike Length on Chromosome 2D in Common Wheat (L.). PhD Dissertation of China Agricultural University, Beijing, China, 2019. pp 38–39 (in Chinese with English abstract).
[21] Chai L L, Chen Z Y, Bian R L, Zhai H J, Cheng X J, Peng H R, Yao Y Y, Hu Z R, Xin M M, Guo W L, Sun Q X, Zhao A J, Ni Z F. Dissection of two quantitative trait loci with pleiotropic effects on plant height and spike length linked in coupling phase on the short arm of chromosome 2D of common wheat (L.).,2018, 131: 2621–2637.
[22] Zhao Z H, Sun H G, Song W, Lu M, Huang J, Wu L F, Wang X M, Li H J. Genetic analysis and detection of the gene MlLX99 on chromosome 2BL conferring resistance to powdery mildew in the wheat cultivar Liangxing 99., 2013, 126: 3081–3089.
[23] 鄒景偉, 邱丹, 孫艷玲, 鄭超星, 李靜婷, 吳培培, 武小菲, 王曉鳴, 周陽, 李洪杰.——小麥品種良星99抗白粉病基因的有效性. 作物學報, 2017, 43: 332–342.Zou J W, Qiu D, Sun Y L, Zheng C X, Li J T, Wu P P, Wu X F, Wang X M, Zhou Y, Li H J.: effectiveness of the gene conferring resistance to powdery mildew in wheat cultivar Liangxing 99., 2017, 43: 332–342 (in Chinese with English abstract).
[24] Wu P P, Hu J H, Zou J W, Qiu D, Qu Y F, Li Y H, Li T, Zhang H J, Yang L, Liu H W, Zhou Y, Zhang Z J, Li, J T, Liu Z Y, Li H J. Fine mapping of the wheat powdery mildew resistance geneusing comparative genomics analysis and the Chinese Spring reference genomic sequence., 2019, 132: 1451–1461.
[25] Wang Y, Song F, Zhu J, Zhang S, Yang Y, Chen T, Tang B, Dong L, Ding N, Zhang Q, Bai Z, Dong X, Chen H, Sun M, Zhai S, Sun Y, Yu L, Lan L, Xiao J, Fang X, Lei H, Zhang Z, Zhao W. GSA: genome sequence archive., 2017, 15: 14–18.
[26] National genomics data center members and partners., 2020, 48: D24–D33.
[27] Wang W X, Wang Z H, Li X T, Ni Z F, Hu Z R, Xin M M, Peng H R, Yao Y Y, Sun Q X, Guo W L. SnpHub: an easy-to-set-up web server framework for exploring large-scale genomic variation data in the post-genomic era with applications in wheat., 2020, 9: giaa060.
附圖1 濟麥22(A)與良星99(B)全基因組reads覆蓋度的密度分布直方圖
X軸為經(jīng)歸一化的每1 Mb內(nèi)的平均reads覆蓋度; Y軸為覆蓋度的密度。
X-axis, the averaged read coverage per 1 Mb after normalization. Y-axis, the distribution density. A and B represent density plot of bin-wise normalized average read coverage in Jimai 22 and Liangxing 99, respectively.
附圖2 濟麥22與良星99間的單區(qū)間差異SNP的密度分布圖
X軸, 每Mb區(qū)間內(nèi)差異SNP位點數(shù)的對數(shù)(以10為底); Y軸, 位點數(shù)的密度。
X-axis, 10-based logarithm of the counts of differential SNP sites per Mb. Y-axis, the distribution density.
Comparative analysis of the genomic sequences between commercial wheat varieties Jimai 22 and Liangxing 99
YANG Zheng-Zhao, WANG Zi-Hao, HU Zhao-Rong, XIN Ming-Ming, YAO Ying-Yin, PENG Hui-Ru, YOU Ming-Shan, SU Zhen-Qi*, and GUO Wei-Long*
College of Agronomy and Biotechnology / State Key Laboratory for Agrobiotechnology / Key Laboratory of Crop Heterosis and Utilization, Ministry of Education, China Agricultural University, Beijing 100193, China
Jimai 22 and Liangxing 99 are high-yield wheat varieties widely planted in the North Huang-Huai Rivers Valley Winter Wheat Zone and Northern Winter Wheat Zone in China, and are currently used as important parents in wheat breeding programs. Although the origins and pedigrees of Jimai 22 and Liangxing 99 are different, they are highly similar in many important agronomic traits, yield-associated traits, and so on. To identify the genomic differences between the two varieties, we performed whole-genome sequencing using the Illumina HiSeq2500 platform, with an average sequencing depth of 5.8×. We aligned the raw sequencing data against the Chinese Spring reference genome and identified the difference of copy-number variation (CNV) regions, single-nucleotide polymorphisms (SNPs) and InDels in sequence between the two varieties. Lengths of 466 Mb CNV intervals were shared by the two varieties. The lengths of cultivar-specific CNV intervals in Jimai 22 and Liangxing 99 were 91 Mb and 45 Mb, respectively, and these intervals are mainly located on chromosomes 2B and 4B. Beyond the CNV intervals, 1,547,371 SNPs and 137,817 InDels were different between the two cultivars. Based on the distribution of SNP densities in the intervals, we identified the polymorphic hotspot regions on chromosomes 1D, 2B, and 4B, making up 14.2% of the whole genome. The sequences of five previous cloned dwarf genes and spike length related genes were investigated, and two genes located in the polymorphic hotspot regions were detected with the frame shift variations. This study provides an important guidance for evaluating the genetic differences between two wheat varieties in the genomic level, and also identified both genetic similarity regions and polymorphic hotspot regions between Jimai 22 and Liangxing 99, which provided a valuable genetic information for future genetic improvement utilizing Jimai 22 and Liangxing 99 as parents.
wheat; Jimai 22; Liangxing 99; whole-genome resequencing; SNP; CNV
本研究由國家重點研發(fā)計劃項目(2018YFD0100803), 國家自然科學基金項目(31701415)和中央高?;究蒲袠I(yè)務(wù)費專項資助。
This study was supported by the National Key Research and Development Program of China (2018YFD0100803), the National Natural Science Foundation of China (31701415), and the Chinese Universities Scientific Fund.
郭偉龍, E-mail: guoweilong@cau.edu.cn; 宿振起, E-mail: suzhenqi80@163.com
E-mail: yangzhengzhao@cau.edu.cn
2020-01-15;
2020-06-02;
2020-08-10.
URL: https://kns.cnki.net/kcms/detail/11.1809.S.20200810.1525.002.html
10.3724/SP.J.1006.2020.01009