亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        不同基因型小麥凍害無人機(jī)遙感高通量表型

        2023-05-15 07:53:56劉易雪吳建輝韓德俊蘇寶峰
        關(guān)鍵詞:分析

        劉易雪,蔚 睿,吳建輝,韓德俊,蘇寶峰

        ·農(nóng)業(yè)信息與電氣技術(shù)·

        不同基因型小麥凍害無人機(jī)遙感高通量表型

        劉易雪1,2,3,蔚 睿4,5,吳建輝4,5,韓德俊4,5,蘇寶峰1,2,3※

        (1. 西北農(nóng)林科技大學(xué)機(jī)械與電子工程學(xué)院,楊凌 712100;2. 農(nóng)業(yè)農(nóng)村部農(nóng)業(yè)物聯(lián)網(wǎng)重點(diǎn)實(shí)驗(yàn)室,楊凌 712100; 3. 陜西省農(nóng)業(yè)信息感知與智能服務(wù)重點(diǎn)實(shí)驗(yàn)室;4. 西北農(nóng)林科技大學(xué)農(nóng)學(xué)院,楊凌 712100; 5. 西北農(nóng)林科技大學(xué)旱區(qū)作物逆境生物學(xué)國家重點(diǎn)實(shí)驗(yàn)室,楊凌 712100)

        為了實(shí)現(xiàn)田間條件下小麥抗凍性狀相關(guān)的數(shù)量性狀基因座(quantitative trait locus, QTL)分析,該研究針對4個(gè)試驗(yàn)地491份小麥核心種質(zhì)資源的抗凍性狀,基于無人機(jī)多光譜遙感提出了一種高通量表型方法。首先通過光譜植被指數(shù)對小麥抗凍性狀進(jìn)行評估,基于機(jī)器學(xué)習(xí)分類算法使用16個(gè)光譜植被指數(shù)特征構(gòu)建了小麥凍害評價(jià)模型,并完成了光譜特征相關(guān)性分析及對評價(jià)模型的貢獻(xiàn)率分析。對比隨機(jī)森林(random forests,RF)、分布式梯度增強(qiáng)(extreme gradient boosting,XGBoost)、梯度提升決策樹(gradient boosting decision tree,GBDT)及支持向量機(jī)(support vector machine,SVM)算法建立的小麥凍害等級(jí)評價(jià)模型,結(jié)果表明,使用XGBoost建立的評價(jià)模型準(zhǔn)確率最高,達(dá)67.94%;16個(gè)光譜特征相關(guān)性及其對評價(jià)模型的貢獻(xiàn)率分析表明,簡化冠層葉綠素含量指數(shù)(simplified canopy chlorophyii content index, SCCCI)對小麥抗凍表型鑒定的貢獻(xiàn)率最大。其次,使用SCCCI作為小麥抗凍表型,結(jié)合通過全基因組關(guān)聯(lián)分析檢測小麥抗凍相關(guān)QTL,檢測到3個(gè)已被證明與抗凍性狀相關(guān)的QTL,證明了基于無人機(jī)獲取的光譜特征可以作為小麥抗凍表型定性定量分析指標(biāo),可提供小麥抗凍性狀遺傳解析必需的表型信息。小麥凍害的無人機(jī)遙感高通量表型方法的提出促進(jìn)了小麥抗凍基因功能解鎖。

        無人機(jī);遙感;小麥凍害;多光譜;關(guān)聯(lián)分析;機(jī)器學(xué)習(xí)

        0 引 言

        小麥高通量表型方法的研究對加速育種進(jìn)程有重要意義,但小麥育種技術(shù)面臨著周期長、效率低、遺傳背景狹窄等境況[1-3]。育種水平的提高是作物產(chǎn)能提高的突破口,小麥的表型分析是挖掘性狀相關(guān)基因功能的必要環(huán)節(jié)[4]。作物的田間表型不僅是對育種前期的種質(zhì)篩選的指導(dǎo)數(shù)據(jù),同時(shí)也是在后期推廣種植中的評估數(shù)據(jù)[5-6]??焖?、準(zhǔn)確、實(shí)時(shí)地獲取大田環(huán)境下小麥育種群體的表型,并將表型參數(shù)與基因型、環(huán)境因子進(jìn)行關(guān)聯(lián)性分析,分析小麥田間表型與基因型的關(guān)系是加速育種進(jìn)程的重要突破口。

        低溫脅迫是小麥在越冬期遇到的最大生存阻力,小麥能否成功越冬主要受抗凍性基因調(diào)控[7-9]。小麥抗凍性相關(guān)基因位點(diǎn)的定位工作依賴于小麥田間受低溫脅迫后凍害表型分析,過去主要由育種專家長期的田間調(diào)查完成[10-12],數(shù)百種基因型數(shù)據(jù)的田間調(diào)查工作需要幾天的時(shí)間才能完成,并且表型數(shù)據(jù)會(huì)受到變化的環(huán)境條件(風(fēng),太陽角,溫度,濕度)以及晝夜節(jié)律的影響[13],存在效率低且數(shù)據(jù)質(zhì)量難以保證的問題[5,14-17]。田間表型平臺(tái)可以代替人工調(diào)查,提高效率[11]。盡管目前已有表型平臺(tái)投入到實(shí)際應(yīng)用中,仍無法保證數(shù)據(jù)質(zhì)量。例如,ANDRADE-SANCHEZ等[12]的研究發(fā)現(xiàn),在5 h的過程中經(jīng)過干旱脅迫的棉花的歸一化差異植被指數(shù)(normalized difference vegetation index, NDVI)下降了21%,而未經(jīng)處理的棉花也下降7%。

        盡管許多研究表明基于遙感的作物田間高通量表型方法可以為解決效率不足、質(zhì)量難以保證的問題提供方案[18-22],但大部分研究仍停留在表型分析階段,如何結(jié)合遺傳學(xué)全基因組關(guān)系分析方法處理表型分析結(jié)果仍是分析植物表型與基因型的關(guān)系追蹤遺傳變異的瓶頸。

        從2020年以來,已有研究開始嘗試結(jié)合基于遙感的植物高通量表型方法和全基因組關(guān)聯(lián)分析(genome-wide association study, GWAS)方法挖掘單核苷酸多態(tài)性標(biāo)記(single nucleotide polymorphism, SNP)[23-25]。例如,SANTINI等[26]基于無人機(jī)遙感開發(fā)了一種半自動(dòng)的地中海白松冠層分割方法,并利用分割結(jié)果進(jìn)行GWAS找到了相關(guān)基因組區(qū)域。COUPEL-LEDRU等[27]基于近地遙感獲取蘋果的冠層結(jié)構(gòu)相關(guān)參數(shù),基于GWAS篩選相關(guān)基因組區(qū)域找到了參與每個(gè)性狀相關(guān)途徑的候選基因。針對基于無人機(jī)遙感的小麥SNP挖掘,HASSAN等[24]基于無人機(jī)遙感獲取田間小麥株高,驗(yàn)證了以往通過人工測量方式獲取到的表型鑒定到的相關(guān)基因組區(qū)域;CHEN等[9]基于無人機(jī)遙感觀察冬小麥的越冬能力并根據(jù)視覺算法結(jié)果基于GWAS找到了相關(guān)基因組區(qū)域。雖然這些研究驗(yàn)證了無人機(jī)遙感可以替代人工田間調(diào)查,證明了基于無人機(jī)遙感的田間高通量表型方法服務(wù)全基因組關(guān)聯(lián)分析的可能性,但由于表型特征信息不足,導(dǎo)致在基因功能挖掘結(jié)果方面沒有突破進(jìn)展。而JIANG等[28]基于無人機(jī)遙感影像構(gòu)建水稻抗旱指數(shù)等3種動(dòng)態(tài)性狀,通過全基因組關(guān)聯(lián)分析檢測到111個(gè)與干旱脅迫顯著相關(guān)位點(diǎn),其中30.6%在之前研究未檢測到。這表明了結(jié)合無人機(jī)遙感以及機(jī)器學(xué)習(xí)算法在表型特征上的創(chuàng)新為突破目前植物基因功能解鎖貢獻(xiàn)了主要?jiǎng)恿?。然而,這部分工作在小麥抗凍性的基因功能挖掘方面還留有空白,如何分析并解釋表型特征與基因型的關(guān)系,以挖掘更多性狀相關(guān)功能位點(diǎn)仍然是目前亟待解決的問題。

        針對上述問題,本文提出一種基于無人機(jī)遙感的小麥抗凍高通量表型方法,以自然條件下受到低溫脅迫的小麥為研究對象,結(jié)合遙感技術(shù)、計(jì)算機(jī)視覺算法以及機(jī)器學(xué)習(xí)提出田間高通量小麥凍害等級(jí)評價(jià)方法,使用GWAS處理表型分析結(jié)果檢測與小麥抗凍性相關(guān)QTL(quantitative trait locus),擬為分析小麥抗凍性的關(guān)鍵遺傳因素提供參考。

        1 材料和方法

        1.1 試驗(yàn)材料和試驗(yàn)區(qū)概況

        小麥群體試驗(yàn)材料包括農(nóng)家種、主栽品種、新品系、核心種質(zhì)以及國外品種等,共491種不同基因型的小麥品種。試驗(yàn)材料種植于4個(gè)不同的試驗(yàn)地,分別為河南省洛陽市農(nóng)科院試驗(yàn)地(112°28.428′E,34°38.352′N),陜西省咸陽市楊陵區(qū)曹新莊試驗(yàn)農(nóng)場(108°5.641′E,34°18.267′N),河南省南陽市農(nóng)科院試驗(yàn)地(112°24.960′E,32°54.582′N)以及江蘇省宿遷市瑞華農(nóng)業(yè)科技有限公司試驗(yàn)地(118°22.650′E,34°2.040′N)。

        育種試驗(yàn)設(shè)計(jì)采取隨機(jī)增廣設(shè)計(jì)[29],每個(gè)試驗(yàn)地有560個(gè)小區(qū),包括16個(gè)區(qū)組,每個(gè)區(qū)組內(nèi)隨機(jī)分布5個(gè)不同的對照品種[30](周麥18、百農(nóng)207、濟(jì)麥22號(hào)、偃展4110、西農(nóng)511),其余小區(qū)為試驗(yàn)品種。4個(gè)試驗(yàn)地皆處于黃淮冬麥區(qū),氣象數(shù)據(jù)由田間微型氣象站獲取,氣候環(huán)境如表1所示。

        表1 試驗(yàn)地氣候環(huán)境對比

        1.2 數(shù)據(jù)采集

        1.2.1 無人機(jī)遙感影像

        無人機(jī)遙感影像采集系統(tǒng)由Matrice 100(深圳DJI)和Micasence Rededge Camera(美國西雅圖MicaSense Inc.)組成。Matrice 100搭載MicaSence Rededge相機(jī)用于快速獲取小麥冠層的多光譜圖像。MicaSence Rededge多光譜相機(jī)分辨率為1 280×960像素,包括藍(lán)光(Blue)、綠光(Green)、紅光(Red)、紅邊(Red edge)和近紅外(NIR)5個(gè)波段,它們的中心波長分別為475、560、668、717和840 nm,其帶寬分別為20、20、10、10和40 nm。

        在小麥越冬期(播種后90 d左右),在晴天10:00至12:00進(jìn)行小麥冠層的無人機(jī)多光譜圖像采集,飛行高度為15 m,飛行路徑覆蓋整個(gè)試驗(yàn)區(qū),同時(shí)將飛行路徑縱向重疊率設(shè)置為80%,橫向重疊率設(shè)置為75%。每次飛行任務(wù)執(zhí)行前后,使用Micasence校正白板獲取參照圖像用于反射率校準(zhǔn)(校正參數(shù)Blue為0.57,Green為0.58,Red為0.58,Red edge為0.58,NIR為0.53)。

        1.2.2 人工田間調(diào)查

        依據(jù)氣象情況,在2019年1月15日,中國北方地區(qū)大面積遭受一次寒流,于2019年1月19日到1月22日分別前往4個(gè)試驗(yàn)地進(jìn)行田間調(diào)查并記錄小麥?zhǔn)艿蜏孛{迫后凍害發(fā)生的具體情況。人工調(diào)查時(shí)間平均為2 d,每個(gè)試驗(yàn)地的調(diào)查間隔最長不超過2 d。在無人機(jī)采集小麥冠層多光譜圖像數(shù)據(jù)的同一時(shí)間,參照中華人民共和國農(nóng)業(yè)行業(yè)標(biāo)準(zhǔn)《NY/T 1301-2007》對4個(gè)試驗(yàn)地鑒定小麥?zhǔn)艿蜏氐拿{迫情況,并將小麥?zhǔn)軆龊η闆r分為5個(gè)等級(jí),分級(jí)標(biāo)準(zhǔn)如表2所示。

        表2 小麥凍害等級(jí)分級(jí)標(biāo)準(zhǔn)

        1.2.3 SNP獲取

        采用Affymetrix公司開發(fā)的小麥660 K SNP芯片對491份品種材料進(jìn)行基因型檢測,分型工作主要由博奧晶典生物技術(shù)有限公司(北京,http://www.capitalbiotech.com)完成?;谪惾~斯系統(tǒng)發(fā)育法對種群結(jié)構(gòu)進(jìn)行聚類分析,計(jì)算微等位基因頻率(micro allele frequency, MAF)、遺傳多樣性和多態(tài)性信息含量(polymorphism information content, PIC)。為保證基因分型數(shù)據(jù)的準(zhǔn)確性,利用R語言包對標(biāo)記進(jìn)行質(zhì)量控制,過濾缺失率大于20%的標(biāo)記、顯著偏分離的標(biāo)記以及微等位基因頻率低于0.05%的標(biāo)記[31]。

        1.3 研究方法

        本研究技術(shù)路線主要包括3個(gè)階段:1)無人機(jī)多光譜圖像數(shù)據(jù)預(yù)處理,為了提取每個(gè)小麥品種對應(yīng)的種植小區(qū)的冠層光譜特征,首先對無人機(jī)多光譜圖像進(jìn)行預(yù)處理生成光譜植被地圖,然后對單個(gè)小區(qū)冠層進(jìn)行聚類分析實(shí)現(xiàn)冠層分割;2)群體小麥凍害表型鑒定與分析,為了實(shí)現(xiàn)基于冠層光譜特征的小麥抗凍性評價(jià)分析,首先基于機(jī)器學(xué)習(xí)分類算法完成了小麥凍害評價(jià)模型,并對比不同算法的分類性能,然后分析基于無人機(jī)遙感的田間小麥高通量表型獲取的特征重要性;3)表型與基因型關(guān)系分析,為了實(shí)現(xiàn)小麥抗凍性相關(guān)QTL,包括了使用獲取的表型信息結(jié)合試驗(yàn)材料的基因型信息,通過全基因組關(guān)聯(lián)分析定位與小麥抗凍性相關(guān)位點(diǎn)。具體技術(shù)路線圖如圖1所示。

        圖1 技術(shù)路線圖

        1.3.1 單個(gè)品種冠層植被指數(shù)特征提取

        使用Pix4D Mapper軟件處理無人機(jī)多光譜圖像,經(jīng)幾何校正、輻射校正、圖像拼接及指數(shù)計(jì)算后,生成楊凌、洛陽、南陽和宿遷4個(gè)試驗(yàn)地對應(yīng)的光譜植被指數(shù)地圖。

        低溫脅迫對小麥生長的影響主要體現(xiàn)在含水量、光合作用以及滲透物質(zhì)含量等方面。小麥?zhǔn)艿降蜏孛{迫后,主要表現(xiàn)為葉尖發(fā)黃、光合作用停滯、葉片枯死、分蘗凍死[32],推測冠層植被覆蓋度或反應(yīng)葉綠素水平的光譜特征可能與小麥?zhǔn)艿蜏孛{迫后的表型相關(guān)。因此,小麥抗凍性鑒定模型的輸入特征包括5個(gè)波段(Blue、Green、Red、Red edge及NIR)和由這5個(gè)波段計(jì)算得出的與冠層覆蓋度、葉綠素相關(guān)的16個(gè)植被指數(shù)[33-46](表2)。

        表2 光譜特征及其計(jì)算式

        注:R、R、R、RR分別代表了Blue、Green、Red、Red edge、NIR正射影像的灰度值。

        Note:R,R,R,R,Rrepresent the gray value of blue, green, red, red edge and NIR orthophoto respectively.

        1.3.2 冠層光譜特征與低溫脅迫的相關(guān)性分析

        本研究獲取到表型數(shù)據(jù)包括人工田間調(diào)查的小麥?zhǔn)艿降蜏孛{迫后的凍害等級(jí)結(jié)果以及基于無人機(jī)遙感獲取的16個(gè)光譜特征結(jié)果。首先,基于機(jī)器學(xué)習(xí)分類算法證明基于無人機(jī)遙感的小麥凍害表型鑒定方法可以代替人工田間調(diào)查法,并分析光譜特征對小麥凍害等級(jí)鑒定的貢獻(xiàn)。

        整合4個(gè)試驗(yàn)地491種不同基因型的小麥品種的16個(gè)植被指數(shù)特征構(gòu)建數(shù)據(jù)集,一共包括1 964個(gè)樣本,每個(gè)樣本有16個(gè)輸入特征,同時(shí)以人工田間調(diào)查結(jié)果為參照作為實(shí)際標(biāo)簽。隨機(jī)抽取數(shù)據(jù)集的20%作為測試集,剩余的數(shù)據(jù)作為訓(xùn)練集,測試集不參與分類器訓(xùn)練。

        本文使用隨機(jī)森林[47](random forests,RF)、分布式梯度增強(qiáng)[48](extreme gradient boosting,XGBoost)、梯度提升決策樹[49](gradient boosting decision tree,GBDT)及支持向量機(jī)[50](support vector machine,SVM)分類算法進(jìn)行訓(xùn)練。XGBoost是在Boosting框架下的前向迭代模型,訓(xùn)練弱學(xué)習(xí)器的同時(shí)引入列采樣和學(xué)習(xí)率,使用節(jié)點(diǎn)遞歸分裂的貪心準(zhǔn)則來實(shí)現(xiàn)樹的生成,同時(shí)使用稀疏感知策略來應(yīng)對缺失值。XGBoost通過融合隨機(jī)森林算法中列采樣,對每次節(jié)點(diǎn)分裂前進(jìn)行隨機(jī)采樣提升抗噪能力,同時(shí)在樹模型的參數(shù)化過程中加入了正則項(xiàng)控制模型的復(fù)雜度防止過擬合[47-49]。設(shè)置RF分類器的弱學(xué)習(xí)器的最大迭代次數(shù)為10,子結(jié)點(diǎn)上的最小樣本數(shù)量為2;設(shè)置XGBoost分類器學(xué)習(xí)率為0.1以二分類邏輯回歸為目標(biāo)函數(shù)訓(xùn)練分類器;設(shè)置GBDT分類器學(xué)習(xí)率為0.1,每棵子樹的深度為2,子結(jié)點(diǎn)上的最小樣本數(shù)量為2;設(shè)置SVM分類器的核函數(shù)為徑向基函數(shù)。使用Python3.8.8在配備 Intel i7-9700 3.00 GHz CPU、16 GB 內(nèi)存和 Nvidia GeForce GTX 3090 顯卡、運(yùn)行 Win10 操作系統(tǒng)的工作站上實(shí)現(xiàn)小麥凍害評價(jià)模型的建立。

        以混淆矩陣評估4種分類器的預(yù)測性能。首先計(jì)算每一個(gè)類別對應(yīng)的召回率、精確度,同時(shí)考慮到數(shù)據(jù)集類間不平衡,使用準(zhǔn)確率、Weight-F1分?jǐn)?shù)(F1)、Micro-F1(F1)分?jǐn)?shù)評價(jià)小麥?zhǔn)艿蜏孛{迫嚴(yán)重度分類模型的精度[51-52],計(jì)算式如下:

        式中T表示預(yù)測值和標(biāo)簽真值均為真;F表示標(biāo)簽真值為真,預(yù)測值為假;F表示標(biāo)簽真值為假,預(yù)測值為真;T表示標(biāo)簽真值為假,預(yù)測值為假為準(zhǔn)確率,%,為精確度,為召回率,為類別數(shù)目權(quán)重,為類別數(shù)目。

        1.3.3 全基因組關(guān)聯(lián)分析

        利用GEMMA軟件采取合適的模型進(jìn)行全基因組關(guān)聯(lián)分析,值的閾值基于有效標(biāo)記的數(shù)量(=1/N,N為有效SNP標(biāo)記數(shù)量)。顯著性關(guān)聯(lián)標(biāo)記通過軟件R 3.0.3繪制的曼哈頓圖(Manhattan plot)來挖掘,通過GCTA軟件進(jìn)行顯著性關(guān)聯(lián)標(biāo)記對表型貢獻(xiàn)率的評估。使用分位數(shù)圖對關(guān)聯(lián)分析結(jié)果進(jìn)行可視化。

        2 結(jié)果與分析

        2.1 冠層光譜特征與低溫脅迫的相關(guān)性

        2.1.1 分類算法對比結(jié)果與分析

        以準(zhǔn)確率、Weight-F1分?jǐn)?shù)作為評價(jià)標(biāo)準(zhǔn),評估4種分類器在測試集上的預(yù)測結(jié)果,如表3所示。相比于經(jīng)典的機(jī)器學(xué)習(xí)分類算法SVM,XGBoost的預(yù)測性能更好。測試結(jié)果表明,XGBoost的準(zhǔn)確率比SVM高5.47個(gè)百分點(diǎn)?;贕BDT的機(jī)器學(xué)習(xí)方法是目前分析預(yù)測問題的高效方法,RF通過集成多棵決策樹實(shí)現(xiàn)分類,而GBDT也是基于提升樹實(shí)現(xiàn)分類,同樣基于提升樹的XGBoost相比于RF的準(zhǔn)確率高出3.16個(gè)百分點(diǎn),比GBDT高出1.87個(gè)百分點(diǎn)。

        由分類器SVM、RF、XGBoost及GBDT在測試集上的預(yù)測結(jié)果與實(shí)際標(biāo)簽生成的混淆矩陣[53](圖2)來看,數(shù)據(jù)集中一級(jí)和二級(jí)的數(shù)量遠(yuǎn)遠(yuǎn)超過其他3個(gè)類別的數(shù)量之和,不平衡的類別數(shù)量為分類器的訓(xùn)練帶來了負(fù)面影響,對比表3中的Weight-F1分?jǐn)?shù)和Micro-F1分?jǐn)?shù)結(jié)果,使用Micro-F1分?jǐn)?shù)評估分類器的預(yù)測性能時(shí),可能會(huì)受到數(shù)量較多的類別的結(jié)果的影響,而Weight-F1是根據(jù)類別數(shù)目所占權(quán)重來計(jì)算的,因此雖然Micro-F1分?jǐn)?shù)結(jié)果要比Weight-F1分?jǐn)?shù)的結(jié)果好,但Weight-F1的結(jié)果更接近實(shí)際情況。

        由表3可知,XGBoost分類器的分類性能優(yōu)于SVM、RF及GBDT,準(zhǔn)確率達(dá)67.94%,各個(gè)凍害等級(jí)的F1分?jǐn)?shù)均超過了0.5(表4),這個(gè)試驗(yàn)結(jié)果表明使用無人機(jī)多光譜圖像基于XGBoost建立的小麥凍害評價(jià)模型能夠替代人工田間調(diào)查作為田間高通量表型獲取方法。在本研究中,對于任一個(gè)試驗(yàn)地,專家田間調(diào)查560個(gè)種植小區(qū),評估小麥冠層受低溫脅迫情況,并鑒定凍害等級(jí),至少需要8 h以上,然而,無人機(jī)執(zhí)行一次飛行任務(wù)僅需要15 min,結(jié)合本文所提全自動(dòng)的表型分析方法分析,研究基于無人機(jī)遙感的小麥田間高通量表型方法將為育種技術(shù)的發(fā)展提供一大助力。

        表3 不同分類器預(yù)測性能

        圖2 不同分類器在測試集上的預(yù)測結(jié)果與實(shí)際標(biāo)簽的混淆矩陣

        表4 基于XGBoost不同凍害等級(jí)分類結(jié)果評估

        2.1.2 光譜特征重要性分析

        用XGBoost函數(shù)庫中對每個(gè)子節(jié)點(diǎn)分裂時(shí)16個(gè)特征帶來的信息增益進(jìn)行統(tǒng)計(jì),將其作為排序索引后,得到了小麥冠層的16個(gè)光譜特征對小麥凍害鑒定的貢獻(xiàn)率排序(圖3)。

        圖3 小麥冠層16個(gè)光譜特征對凍害鑒定的貢獻(xiàn)率排序

        不含抗凍基因的小麥品種在受到低溫脅迫后,隨著凍害程度的加重,葉片含水量、葉綠素a、光系統(tǒng)Ⅱ最大光化學(xué)量子產(chǎn)量均顯著下降[54],小麥的光合作用受到嚴(yán)重影響,光合作用相關(guān)參數(shù)發(fā)生變化,同時(shí),硬化組織的酶、蛋白也會(huì)發(fā)生變化[55]。而SCCCI常用于反演植物葉綠素[41-42,45],相比于圖3中其他與植被覆蓋度相關(guān)的植被指數(shù)更符合如小麥?zhǔn)艿蜏孛{迫后光合作用停滯葉綠素含量降低的生理生化反應(yīng)。

        2.2 全基因組關(guān)聯(lián)分析

        通過全基因組關(guān)聯(lián)分析發(fā)現(xiàn)在小麥的21條染色體上,共有3個(gè)位點(diǎn)出現(xiàn)了超出閾值(﹣lg=4)且連續(xù)分布的顯著SNP,分別位于位點(diǎn)2B、3A與5A(圖4a)。從圖4b中看出,當(dāng)- lg<4時(shí),值的分布和均勻分布的結(jié)果集中在一條直線上,表示確定與抗凍表型性狀不關(guān)聯(lián)的位點(diǎn),這些位點(diǎn)的值觀測值與期望值一致;當(dāng)-lg>4時(shí),值的分布和均勻分布的結(jié)果出現(xiàn)快速分離的情況,特別是值越低的時(shí)候分離程度就越高,說明這些位點(diǎn)的效應(yīng)超過了隨機(jī)效應(yīng),位點(diǎn)2B、3A與5A是潛在與抗凍表型性狀相關(guān)的候選位點(diǎn)。對2B短臂上顯著SNP注釋發(fā)現(xiàn)他們位于編碼一個(gè)Cor14b蛋白的基因上,該蛋白已被證明在大麥和二倍體小麥中參與對逆境脅迫的響應(yīng)[56];同樣在5A長臂的位點(diǎn)與前人定位研究的Fr1與Fr2重合,F(xiàn)r1與Fr2是兩個(gè)重要的抗凍(frost tolerance, FroT)位點(diǎn),即抗凍1(FR-A1)和抗凍2(FR-A2),擁有Fr1與Fr2的小麥具有耐凍性和冬季抗凍性,其中Fr1與VRN-A1基因緊密相連,他們共同作用參與到小麥耐凍調(diào)控的通路上[57-58]。3AL上的位點(diǎn)很可能與近期挖掘的一個(gè)小麥抗凍QTL一致[59]。

        圖4 SCCCI的曼哈頓圖與QQ Plot圖

        上述結(jié)果表明,結(jié)合基于無人機(jī)遙感的田間小麥高通量表型方法獲取的表型數(shù)據(jù)與基因型進(jìn)行全基因組關(guān)聯(lián)分析,QTL結(jié)果和前人的研究結(jié)果實(shí)現(xiàn)了一致,證明光譜植被指數(shù)可以很好地定性且定量小麥的受凍情況,可作為小麥抗凍種質(zhì)資源篩選與抗凍遺傳解析的表型評定指標(biāo)。

        3 結(jié) 論

        1)本文對4個(gè)試驗(yàn)地的491份小麥種質(zhì)材料進(jìn)行了人工田間調(diào)查,同時(shí)獲取4個(gè)試驗(yàn)地的無人機(jī)多光譜遙感影像并生成16個(gè)光譜植被指數(shù)特征,對比分析了XGBoost、GBDT、RF及SVM 4種方法建立的小麥凍害評價(jià)模型,結(jié)果表明,使用XGBoost建立的評價(jià)模型準(zhǔn)確率最高,可達(dá)67.94%。同時(shí),并以信息增益作為評估標(biāo)準(zhǔn)評估了16個(gè)不同光譜特征對評價(jià)模型分類性能的影響,結(jié)果表明簡化冠層葉綠素含量指數(shù)(simplified canopy chlorophyii content index, SCCCI)是對小麥凍害鑒定貢獻(xiàn)最大的特征。

        2)基于無人機(jī)遙感的小麥抗凍表型方法可用于抗凍性狀相關(guān)的遺傳解析。使用GWAS對小麥挖掘抗凍基因SNP,檢測結(jié)果表明,位于2B、3A與5A的3個(gè)位點(diǎn)出現(xiàn)了超出閾值且連續(xù)分布的顯著SNP,證明基于無人機(jī)獲取的光譜植被指數(shù)可以很好地定性且定量小麥?zhǔn)軆龊η闆r,可作為小麥抗凍種質(zhì)資源篩選與抗凍遺傳解析的表型評定指標(biāo)。

        [1] 鄒華文,王道杰,邊秀秀,等. 近15年國家自然科學(xué)基金稻、麥類作物遺傳育種領(lǐng)域項(xiàng)目申請情況分析[J]. 作物學(xué)報(bào),2015,41(5):820-828.

        ZOU Huawen, WANG Daojie, BIAN Xiuxiu, et al. Analysis of NSFC program applications in rice and wheat crops genetics and breeding fields in recent 15 years[J]. Acta Agronomica Sinica, 2015, 41(5): 820-828. (in Chinese with English abstract)

        [2] 翟俊鵬,李海霞,畢惠惠,等. 普通小麥主要農(nóng)藝性狀的全基因組關(guān)聯(lián)分析[J]. 作物學(xué)報(bào),2019,45(10):1488-1502.

        ZHAI Junpeng, LI Haixia, BI Huihui, et al. Genome-wide association study for main agronomic traits in common wheat[J]. Acta Agronomica Sinica, 2019, 45(10): 1488-1502. (in Chinese with English abstract)

        [3] 胡文靜,李東升,裔新,等. 小麥穗部性狀和株高的QTL定位及育種標(biāo)記開發(fā)和驗(yàn)證[J]. 作物學(xué)報(bào),2022,48(6):1346-1356.

        HU Wenjing, LI Dongsheng, YI Xin, et al. Molecular mapping and validation of quantitative trait loci for spike-related traits and plant height in wheat[J]. Acta Agronomica Sinica, 2022, 48(6): 1346-1356. (in Chinese with English abstract)

        [4] FURBANK R T, TESTER M. Phenomics-technologies to relieve the phenotyping bottleneck[J]. Trends in Plant Science, 2011, 16(12): 635-644.

        [5] 劉建剛,趙春江,楊貴軍,等. 無人機(jī)遙感解析田間作物表型信息研究進(jìn)展[J]. 農(nóng)業(yè)工程學(xué)報(bào),2016,32(24):98-106.

        LIU Jiangang, ZHAO Chunjiang, YANG Guijun, et al. Review of field-based phenotyping by unmanned aerial vehicle remote sensing platform[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2016, 32(24): 98-106. (in Chinese with English abstract)

        [6] ARAUS J L, KEFAUVER S C, ZAMAN-ALLAH M, et al. Translating high-throughput phenotyping into genetic gain[J]. Trends In Plant Science, 2018, 23(5): 451-466.

        [7] 孟繁圓,馮利平,張豐瑤,等. 北部冬麥區(qū)冬小麥越冬凍害時(shí)空變化特征[J]. 作物學(xué)報(bào),2019,45(10):1576-1585.

        MENG Fanyuan, FENG Liping, ZHANG Fengyao, et al. Temporal and spatial variations of winter wheat freezing injury in northern winter wheat region[J]. Acta Agronomica Sinica, 2019, 45(10): 1576-1585. (in Chinese with English abstract)

        [8] 王慧芳,顧曉鶴,董瑩瑩,等. 冬小麥凍害災(zāi)情及長勢恢復(fù)的變化向量分析[J]. 農(nóng)業(yè)工程學(xué)報(bào),2011,27(11):145-150.

        WANG Huifang, GU Xiaohe, DONG Yingying, et al. Monitoring freeze injury and growth recovery of winter wheat based on change vector analysis[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2011, 27(11): 145-150. (in Chinese with English abstract)

        [9] CHEN Y, SIDHU H S, KAVIANI M, et al. Application of image-based phenotyping tools to identify QTL for in-field winter survival of winter wheat (L.)[J]. Theoretical and Applied Genetics, 2019, 132(9): 2591-2604.

        [10] ARAUS J L, CAIRNS J E. Field high-throughput phenotyping: The new crop breeding frontier[J]. Trends in Plant Science, 2014, 19(1): 52-61.

        [11] BAO Y, TANG L, BREITZMAN M W, et al. Field‐based robotic phenotyping of sorghum plant architecture using stereo vision[J]. Journal of Field Robotics, 2019, 36(2): 397-415.

        [12] ANDRADE-SANCHEZ P, GORE M A, HEUN J T, et al. Development and evaluation of a field-based high-throughput phenotyping platform[J]. Functional Plant Biology, 2013, 41(1): 68-79.

        [13] FELDERHOFF T J, MURRAY S C, KLEIN P E, et al. QTLs for energy‐related traits in a sweet× grain sorghum [(L.) Moench] mapping population[J]. Crop Science, 2012, 52(5): 2040-2049.

        [14] 李繼宇,胡瀟丹,蘭玉彬,等. 基于文獻(xiàn)計(jì)量學(xué)的2001-2020全球農(nóng)用無人機(jī)研究進(jìn)展[J]. 農(nóng)業(yè)工程學(xué)報(bào),2021,37(9):328-339.

        LI Jiyu, HU Xiaodan, LAN Yubin, et al. Research advance on worldwide agricultural UAVs in 2001-2020 based on bibliometrics[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2021, 37(9): 328-339. (in Chinese with English abstract)

        [15] 劉忠,萬煒,黃晉宇,等. 基于無人機(jī)遙感的農(nóng)作物長勢關(guān)鍵參數(shù)反演研究進(jìn)展[J]. 農(nóng)業(yè)工程學(xué)報(bào),2018,34(24):60-71.

        LIU Zhong, WAN Wei, HUANG Jinyu, et al. Progress on key parameters inversion of crop growth based on unmanned aerial vehicle remote sensing[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2018, 34(24): 60-71. (in Chinese with English abstract)

        [16] 劉軻,黃平,任國業(yè),等.基于冠層反射率模型的作物參數(shù)多階段反演方法研究進(jìn)展[J]. 農(nóng)業(yè)工程學(xué)報(bào),2017,33(1):190-198.

        LIU Ke, HUANG Ping, REN Guoye, et al. Review on multi-stage inversion techniques of canopy reflectance models for retrieving crop variables[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2017, 33(1): 190-198. (in Chinese with English abstract)

        [17] XIE C, YANG C. A review on plant high-throughput phenotyping traits using UAV-based sensors[J]. Computers and Electronics in Agriculture, 2020, 178: 105731.

        [18] HU P, CHAPMAN S C, WANG X, et al. Estimation of plant height using a high throughput phenotyping platform based on unmanned aerial vehicle and self-calibration: Example for sorghum breeding[J]. European Journal of Agronomy, 2018, 95: 24-32.

        [19] AMPATZIDIS Y, PARTEL V. UAV-based high throughput phenotyping in citrus utilizing multispectral imaging and artificial intelligence[J]. Remote Sensing, 2019, 11(4): 410.

        [20] HERZIG P, BORRMANN P, KNAUER U, et al. Evaluation of RGB and multispectral unmanned aerial vehicle (UAV) imagery for high-throughput phenotyping and yield prediction in barley breeding[J]. Remote Sensing, 2021, 13(14): 2670.

        [21] CHIVASA W, MUTANGA O, BURGUE?O J. UAV-based high-throughput phenotyping to increase prediction and selection accuracy in maize varieties under artificial MSV inoculation[J]. Computers and Electronics in Agriculture, 2021, 184: 106128.

        [22] FULLANA-PERICàS M, CONESA M à, GAGO J, et al. High-throughput phenotyping of a large tomato collection under water deficit: Combining UAVs’ remote sensing with conventional leaf-level physiologic and agronomic measurements[J]. Agricultural Water Management, 2022, 260: 107283.

        [23] LYRA D H, VIRLET N, SADEGHI-TEHRAN P, et al. Functional QTL mapping and genomic prediction of canopy height in wheat measured using a robotic field phenotyping platform[J]. Journal of Experimental Botany, 2020, 71(6): 1885-1898.

        [24] HASSAN M A, YANG M, FU L, et al. Accuracy assessment of plant height using an unmanned aerial vehicle for quantitative genomic analysis in bread wheat[J]. Plant Methods, 2019, 15(1): 1-12.

        [25] BREITZMAN M W, BAO Y, TANG L, et al. Linkage disequilibrium mapping of high-throughput image-derived descriptors of plant architecture traits under field conditions[J]. Field Crops Research, 2019, 244: 107619.

        [26] SANTINI F, KEFAUVER S C, ARAUS J L, et al. Bridging the genotype–phenotype gap for a Mediterranean pine by semiautomatic crown identification and multispectral imagery[J]. New Phytologist, 2021, 229(1): 245-258.

        [27] COUPEL-LEDRU A, PALLAS B, DELALANDE M, et al. Tree architecture, light interception and water‐use related traits are controlled by different genomic regions in an apple tree core collection[J]. New Phytologist, 2022, 234(1): 209-226.

        [28] JIANG Z, TU H, BAI B, et al. Combining UAV-RGB high-throughput field phenotyping and genome‐wide association study to reveal genetic variation of rice germplasms in dynamic response to drought stress[J]. New Phytologist, 2021, 232(1): 440-455.

        [29] FEDERER W T, REYNOLDS M, CROSSA J. Combining results from augmented designs over sites[J]. Agronomy Journal, 2001, 93(2): 389-395.

        [30] 中華人民共和國農(nóng)業(yè)部. 農(nóng)作物品種區(qū)域試驗(yàn)技術(shù)規(guī)程小麥:NY/T 1301-2007 [S]. 北京:農(nóng)業(yè)農(nóng)村部,2007.

        [31] MENG L, LI H, ZHANG L, et al. QTL IciMapping: Integrated software for genetic linkage map construction and quantitative trait locus mapping in biparental populations[J]. The Crop Journal, 2015, 3(3): 269-283.

        [32] 王洋洋. 春季低溫對小麥生理和產(chǎn)量的影響及凍害評價(jià)[D]. 鄭州:河南農(nóng)業(yè)大學(xué),2021.

        WANG Yangyang. Effect of Spring Low Temperature on Wheat Physiology and Yield and Evaluation of Freezing Injury[D]. Zhengzhou: Henan Agricultural University, 2021. (in Chinese with English abstract)

        [33] TUCKER C J. Red and photographic infrared linear combinations for monitoring vegetation[J]. Remote Sensing of Environment, 1979, 8(2): 127-150.

        [34] HUETE A, DIDAN K, MIURA T, et al. Overview of the radiometric and biophysical performance of the MODIS vegetation indices[J]. Remote Sensing of Environment, 2002, 83(1/2): 195-213.

        [35] GITELSON A A, GRITZ Y, MERZLYAK M N. Relationships between leaf chlorophyll content and spectral reflectance and algorithms for non-destructive chlorophyll assessment in higher plant leaves[J]. Journal of Plant Physiology, 2003, 160(3): 271-282.

        [36] CRIPPEN R E. Calculating the vegetation index faster[J]. Remote Sensing of Environment, 1990, 34(1): 71-73.

        [37] MACCIONI A, AGATI G, MAZZINGHI P. New vegetation indices for remote measurement of chlorophylls based on leaf directional reflectance spectra[J]. Journal of Photochemistry and Photobiology B: Biology, 2001, 61(1/2): 52-61.

        [38] GITELSON A A, KAUFMAN Y J, MERZLYAK M N. Use of a green channel in remote sensing of global vegetation from EOS-MODIS[J]. Remote sensing of Environment, 1996, 58(3): 289-298.

        [39] HUETE A R. A soil-adjusted vegetation index (SAVI)[J]. Remote Sensing of Environment, 1988, 25(3): 295-309.

        [40] RONDEAUX G, STEVEN M, BARET F. Optimization of soil-adjusted vegetation indices[J]. Remote Sensing of Environment, 1996, 55(2): 95-107.

        [41] GITELSON A A, MERZLYAK M N, Chivkunova O B. Optical properties and nondestructive estimation of anthocyanin content in plant leaves[J]. Photochemistry and Photobiology, 2001, 74(1): 38-45.

        [42] WU C, NIU Z, TANG Q, et al. Estimating chlorophyll content from hyperspectral vegetation indices: Modeling and validation[J]. Agricultural and Forest Meteorology, 2008, 148(8/9): 1230-1241.

        [43] LOUHAICHI M, BORMAN M M, JOHNSON D E. Spatially located platform and aerial photography for documentation of grazing impacts on wheat[J]. Geocarto International, 2001, 16(1): 65-70.

        [44] FILELLA I, SERRANO L, SERRA J, et al. Evaluating wheat nitrogen status with canopy reflectance indices and discriminant analysis[J]. Crop Science, 1995, 35(5): 1400-1405.

        [45] RAPER T B, VARCO J J. Canopy-scale wavelength and vegetative index sensitivities to cotton growth parameters and nitrogen status[J]. Precision Agriculture, 2015, 16(1): 62-76.

        [46] BROGE N H, LEBLANC E. Comparing prediction power and stability of broadband and hyperspectral vegetation indices for estimation of green leaf area index and canopy chlorophyll density[J]. Remote Sensing of Environment, 2001, 76(2): 156-172.

        [47] GONZáLEZ A G, HERRADOR M á, ASUERO A G. Intra-laboratory assessment of method accuracy (trueness and precision) by using validation standards[J]. Talanta, 2010, 82(5): 1995-1998

        [48] CHEN T, GUESTRIN C. Xgboost: A scalable tree boosting system[C]. Proceedings of the 22nd Acm Sigkdd International Conference on Knowledge Discovery and Data Mining. New York: Association for Computing Machinery, 2016: 785-794.

        [49] PERRAS M, SARHAN F. Synthesis of freezing tolerance proteins in leaves, crown, and roots during cold acclimation of wheat[J]. Plant Physiology, 1989, 89(2): 577-585.

        [50] NOBLE W S. What is a support vector machine?[J]. Nature Biotechnology, 2006, 24(12): 1565-1567.

        [51] HABOUDANE D, MILLER J R, TREMBLAY N, et al. Integrated narrow-band vegetation indices for prediction of crop chlorophyll content for application to precision agriculture[J]. Remote Sensing of Environment, 2002, 81(2/3): 416-426.

        [52] HENDERSON C R. Best linear unbiased estimation and prediction under a selection model[J]. Biometrics, 1975, 31(2): 423-447.

        [53] QUINONERO-CANDELA J, RASMUSSEN C E. A unifying view of sparse approximate Gaussian process regression[J]. The Journal of Machine Learning Research, 2005, 6: 1939-1959.

        [54] SZALAI G, PAP M, JANDA T. Light-induced frost tolerance differs in winter and spring wheat plants[J]. Journal of Plant Physiology, 2009, 166(16): 1826-1831.

        [55] FRIEDMAN J H. Greedy function approximation: A gradient boosting machine[J]. Annals of Statistics, 2001, 29(5): 1189-1232.

        [56] CROSATTI C, POLVERINO de LAURETO P, BASSI R, et al. The interaction between cold and light controls the expression of the cold-regulated barley gene cor14b and the accumulation of the corresponding protein[J]. Plant Physiology, 1999, 119(2): 671-680.

        [57] MONROY A F, DRYANOVA A, MALETTE B, et al. Regulatory gene candidates and gene expression analysis of cold acclimation in winter and spring wheat[J]. Plant Molecular Biology, 2007, 64(4): 409-423.

        [58] GALIBA G, VáGúJFALVI A, LI C, et al. Regulatory genes involved in the determination of frost tolerance in temperate cereals[J]. Plant Science, 2009, 176(1): 12-19.

        [59] SOLEIMANI B, LEHNERT H, BABBEN S, et al. Genome wide association study of frost tolerance in wheat[J]. Scientific Reports, 2022, 12(1): 1-11.

        High-throughput phenotyping for different genotype wheat frost using UAV-based remote sensing

        LIU Yixue1,2,3, YU Rui4,5, WU Jianhui4,5, HAN Dejun4,5, SU Baofeng1,2,3

        (1.712100,; 2.712100,; 3.712100; 4.712100; 5.,712100,)

        Wheat (triticum aestivum l.) breeding technology can face a great challenge on the long cycle, low efficiency, and narrow genetic background. An important breakthrough can be combining the high-throughput phenotyping of in-field wheat and genome-wide association, thereby revealing the genetic variation in dynamic response to environmental stress. Fortunately, the unmanned aerial vehicle (UAV) remote sensing and machine learning can be expected to bridge the genotype–phenotype gap of the wheat in the breeding process. Among them, frost tolerance is an important phenotype target, particularly with the winter survival of wheat in various environments. It is a high demand for the rapid and cost-effective assessment of frost tolerance from the UAV multi-spectral imagery using machine learning. In this study, a genome-wide association study (GWAS) was assessed for the quantitative genomic analysis of wheat frost tolerance. A bi-parental wheat population consisting of 491 doubled haploid lines was also used in four study sites. 491 wheat core materials with a relatively consistent growth stage were selected to obtain their high-density genotype data with the 660 K single nucleotide polymorphism (SNP). The UAV-based multi-spectral imagery of the wheat canopy was collected at the overwintering stage at four experimental sites. At the same time, the wheat in-field phenotypes of frost tolerance were investigated by the wheat breeding experts at the same time. The image pre-processing was performed on the features generation of 16 spectral vegetation indices, including image mosaic, geometric correction, radiometric correction and index calculation. Image segmentation was utilized to obtain the features of the wheat canopy using unsupervised clustering. The features correlation analysis and importance analysis were implemented to compare with the in-field investigation, in order to identify quantitative trait loci (QTL) underlying frost tolerance. A comparison was then made on the evaluation models of wheat freezing injury established by random forests (RF), extreme gradient boosting (XGBoost), gradient boosting decision tree (GBDT), and support vector machine (SVM). The results showed that significantly high accuracy was achieved up to 67.94% of the classifier in the XGBoost, compared with the in-field investigation. The correlation and importance of features were also analyzed during this time. The importance of 22 spectral features to the prediction performance of the classifier was evaluated using the information gain brought by the feature, when the sub node of the classifier split. The results showed that there was the most important for the prediction performance of the classifier in the simplified Canopy Chlorophyll content index (SCCCI) among the 16 spectral features of the wheat canopy. Three QTLs were also closely related to the frost resistance detected by the genome-wide association analysis. The three loci of 2B, 3A, and 5A on chromosome 21 of wheat presented a significant SNP, even exceeding the threshold (-lg=4). The SNPs were continuously distributed. Therefore, the spectral features using UAV remote sensing can be expected to serve as the wheat frost resistance QTL. The UAV-enabled phenotyping can be an effective, high-throughput, and cost-effective approach to understanding the genetic basis of wheat frost tolerance in genetic studies and practical breeding. This finding can also provide a fast way for the high-throughput phenotyping of wheat frost tolerance for wheat winter survival in the field.

        UAV; remote sensing; wheat frost; multispectral; GWAS; machine learning

        10.11975/j.issn.1002-6819.202206279

        S127

        A

        1002-6819(2023)-05-0128-09

        劉易雪,蔚睿,吳建輝,等. 不同基因型小麥凍害無人機(jī)遙感高通量表型[J]. 農(nóng)業(yè)工程學(xué)報(bào),2023,39(5):128-136.doi:10.11975/j.issn.1002-6819.202206279 http://www.tcsae.org

        LIU Yixue, YU Rui, WU Jianhui, et al. High-throughput phenotyping for different genotype wheat frost using UAV-based remote sensing[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2023, 39(5): 128-136. (in Chinese with English abstract) doi:10.11975/j.issn.1002-6819.202206279 http://www.tcsae.org

        2022-06-27

        2023-02-22

        國家重點(diǎn)研發(fā)計(jì)劃(2021YFD1200600);楊凌種業(yè)創(chuàng)新中心重點(diǎn)研發(fā)項(xiàng)目(YLzy-xm-01)

        劉易雪,博士生,研究方向?yàn)樘镩g植物表型方法。Email:sunnyliu@nwafu.edu.cn

        蘇寶峰,副教授,博士生導(dǎo)師,研究方向?yàn)樘镩g植物表型快速獲取及應(yīng)用。Email:bfs@nwsuaf.edu.cn

        猜你喜歡
        分析
        禽大腸桿菌病的分析、診斷和防治
        隱蔽失效適航要求符合性驗(yàn)證分析
        電力系統(tǒng)不平衡分析
        電子制作(2018年18期)2018-11-14 01:48:24
        電力系統(tǒng)及其自動(dòng)化發(fā)展趨勢分析
        經(jīng)濟(jì)危機(jī)下的均衡與非均衡分析
        對計(jì)劃生育必要性以及其貫徹實(shí)施的分析
        GB/T 7714-2015 與GB/T 7714-2005對比分析
        出版與印刷(2016年3期)2016-02-02 01:20:11
        網(wǎng)購中不良現(xiàn)象分析與應(yīng)對
        中西醫(yī)結(jié)合治療抑郁癥100例分析
        偽造有價(jià)證券罪立法比較分析
        青草国产精品久久久久久| 国产一区二区三区在线大屁股| 国产亚洲精品久久午夜玫瑰园| 亚洲色成人网站www永久四虎| 欧美日韩综合网在线观看| 精品人妻免费看一区二区三区 | 色综合久久加勒比高清88| 高清亚洲精品一区二区三区| 中文字幕本久久精品一区| 久久国产精品99精品国产| 亚洲 欧美 综合 另类 中字| 激情综合网缴情五月天| 日本一区二区三区光视频| 网禁拗女稀缺资源在线观看| 国产免费破外女真实出血视频| 亚洲av永久无码精品成人| 日韩av一区二区三区精品久久| 国产亚洲一区二区在线观看| 亚洲综合久久成人a片| 无码中文字幕专区一二三| 蜜臀av在线一区二区尤物| 欧美成人看片一区二区三区尤物| 99精品视频在线观看| 国产好片日本一区二区三区四区| 久久蜜桃资源一区二区| 又大又紧又粉嫩18p少妇| 亚洲精品二区中文字幕| 国产一区二区三区涩涩涩| 亚洲免费观看视频| 丰满人妻av无码一区二区三区| www.久久av.com| 激情视频在线观看好大| 亚洲乱码国产乱码精品精| 国产成人精品日本亚洲18| 国产一级一片内射在线| 中文字幕一区二区三区视频 | 五月婷婷激情六月开心| 91久久精品色伊人6882| 少妇装睡让我滑了进去| 成人午夜免费福利| 青青草成人免费在线视频|