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

        ?

        基于酶約束的代謝網(wǎng)絡(luò)模型研究進(jìn)展及其應(yīng)用

        2019-10-31 12:35:52趙欣楊雪毛志濤馬紅武
        生物工程學(xué)報(bào) 2019年10期
        關(guān)鍵詞:成本方法模型

        趙欣,楊雪,毛志濤,馬紅武

        基于酶約束的代謝網(wǎng)絡(luò)模型研究進(jìn)展及其應(yīng)用

        趙欣1,2,楊雪1,2,毛志濤1,2,馬紅武1

        1 中國(guó)科學(xué)院天津工業(yè)生物技術(shù)研究所 中國(guó)科學(xué)院系統(tǒng)微生物技術(shù)重點(diǎn)實(shí)驗(yàn)室,天津 300308 2 中國(guó)科學(xué)院大學(xué),北京 100049

        基因組尺度代謝網(wǎng)絡(luò)模型已經(jīng)成功地應(yīng)用于指導(dǎo)代謝工程改造,但由于傳統(tǒng)通量平衡分析法僅考慮化學(xué)計(jì)量學(xué)和反應(yīng)方向約束,模擬得到的是理論最優(yōu)結(jié)果,對(duì)一些現(xiàn)象如代謝溢流、底物層級(jí)利用等無(wú)法準(zhǔn)確描述。近年來(lái)人們通過(guò)在代謝網(wǎng)絡(luò)模型中引入新的蛋白量、熱力學(xué)等約束發(fā)展了新的約束優(yōu)化計(jì)算方法,可以更準(zhǔn)確真實(shí)地模擬細(xì)胞在不同條件下的代謝行為。文中主要對(duì)近年來(lái)提出的多種酶約束模型進(jìn)行評(píng)述,對(duì)酶約束引入的基本思路、酶約束的數(shù)學(xué)方程表示及優(yōu)化目標(biāo)設(shè)定、引入酶約束后對(duì)代謝通量計(jì)算結(jié)果的影響及酶約束模型在代謝工程菌種改造中的應(yīng)用等進(jìn)行了全面深入的介紹,并提出了已有各種方法存在的主要問(wèn)題,展望了相關(guān)方法的未來(lái)發(fā)展方向。通過(guò)引入新的約束,代謝網(wǎng)絡(luò)模型能夠更精確模擬和預(yù)測(cè)細(xì)胞在環(huán)境和基因擾動(dòng)下的代謝行為,為代謝工程菌種改造提供更準(zhǔn)確可靠的指導(dǎo)。

        酶約束,代謝網(wǎng)絡(luò)模型,代謝工程,酶動(dòng)力學(xué),熱力學(xué)

        基于化學(xué)反應(yīng)計(jì)量學(xué)關(guān)系和反應(yīng)可逆性約束的約束優(yōu)化方法,如通量平衡分析(FBA)[1]、通量可變性分析(FVA)[2-3]等,是目前基因組尺度代謝網(wǎng)絡(luò)模型(GEM) 分析最常用的方法[4-5]。通過(guò)這些分析方法可以計(jì)算滿(mǎn)足特定優(yōu)化目標(biāo)(如最大化細(xì)胞生長(zhǎng)速率、ATP生成或產(chǎn)物生成速率等)的代謝通量分布,并基于結(jié)果分析基因必需性、確定產(chǎn)物合成最優(yōu)途徑、確定代謝工程改造靶點(diǎn)[6]等。雖然這類(lèi)方法已經(jīng)比較廣泛地應(yīng)用在指導(dǎo)代謝工程改造中,但在很多情況下菌體生長(zhǎng)時(shí)都處于一種非最優(yōu)代謝狀態(tài),因此模擬結(jié)果和實(shí)際結(jié)果可能有偏差[7]。而且由FBA等方法預(yù)測(cè)細(xì)胞生長(zhǎng)表型時(shí)需首先設(shè)定底物消耗速率才能計(jì)算最優(yōu)生長(zhǎng)速率,但實(shí)際上在實(shí)驗(yàn)測(cè)定底物輸入速率時(shí)可同時(shí)測(cè)得生長(zhǎng)速率,因此對(duì)于菌體生長(zhǎng)的模擬,這類(lèi)模型更適合于擬合實(shí)驗(yàn)結(jié)果而不是進(jìn)行預(yù)測(cè)。除此之外,由于模型中僅考慮計(jì)量學(xué)和反應(yīng)方向性約束,而真實(shí)細(xì)胞中的代謝流分布還受到許多其他因素影響,因此有很多生物學(xué)現(xiàn)象不能通過(guò)僅考慮計(jì)量學(xué)約束而不考慮其他約束的傳統(tǒng)FBA進(jìn)行準(zhǔn)確模擬預(yù)測(cè),舉例如下。

        1) 代謝溢流現(xiàn)象

        代謝溢流是在有氧條件下,微生物處于高生長(zhǎng)速率或高底物吸收速率時(shí),即使氧供應(yīng)充足仍有發(fā)酵副產(chǎn)物如乙酸、乙醇排出的現(xiàn)象。這種現(xiàn)象通常出現(xiàn)在單細(xì)胞生物當(dāng)中[8],最常見(jiàn)的為釀酒酵母產(chǎn)乙醇和大腸桿菌生成乙酸,被統(tǒng)稱(chēng)為Crabtree effect[8-9]。此外,這種現(xiàn)象在某些癌細(xì)胞中也有相關(guān)報(bào)道,被稱(chēng)為Warburg effect[10]。盡管代謝溢流現(xiàn)象已經(jīng)被廣泛研究,但傳統(tǒng)的代謝網(wǎng)絡(luò)模型由于僅以生長(zhǎng)最大化為目標(biāo),并不能準(zhǔn)確描述這種副產(chǎn)物生成的現(xiàn)象。

        2) 微生物在多種碳源上的層級(jí)利用現(xiàn)象

        細(xì)菌在兩種及以上混合碳源的培養(yǎng)基上培養(yǎng)時(shí)會(huì)表現(xiàn)出“先后利用” (Diauxie) 和“共利用” (Co-utilization) 的不同代謝行為[11-12]。傳統(tǒng)代謝網(wǎng)絡(luò)約束優(yōu)化方法因需要人為設(shè)定底物消耗速率,無(wú)法解釋細(xì)胞為何傾向于利用兩種底物中的一種或同時(shí)利用。需要考慮新的機(jī)理整合新的約束進(jìn)行模擬。

        在復(fù)雜的生物體內(nèi),細(xì)胞代謝受到很多因素的調(diào)控與限制,酶濃度的限制就是其中的一種,即在細(xì)胞內(nèi)催化代謝反應(yīng)的酶的總量是有限的,在1 g細(xì)胞干重中總蛋白的量一般不超過(guò)60%,而催化反應(yīng)的酶又只占總蛋白的一定比例。因此當(dāng)存在酶活很低的酶時(shí),細(xì)胞代謝轉(zhuǎn)化的速率就會(huì)受到限速酶能達(dá)到的最高濃度限制。并且,最近的研究結(jié)果表明,代謝溢流以及底物層級(jí)利用現(xiàn)象的發(fā)生與微生物體內(nèi)的酶濃度限制有關(guān)[13-15]?;谶@一考慮,近年來(lái)科學(xué)家們開(kāi)發(fā)了多種方法將細(xì)胞內(nèi)酶的約束整合到代謝網(wǎng)絡(luò)模型中,由此對(duì)與酶濃度限制相關(guān)的代謝現(xiàn)象進(jìn)行更準(zhǔn)確的模擬和預(yù)測(cè)。不同的方法從不同角度,如酶在細(xì)胞內(nèi)的空間擁擠度、酶的催化活性、代謝反應(yīng)熱力學(xué)等,對(duì)細(xì)胞內(nèi)的酶量約束進(jìn)行數(shù)學(xué)表征,進(jìn)而通過(guò)適當(dāng)?shù)募s束條件或優(yōu)化目標(biāo)求出滿(mǎn)足酶總量約束或者最小化酶成本的代謝通量分布及相應(yīng)酶量分布。這種模型目前已經(jīng)被成功地應(yīng)用到指導(dǎo)菌種代謝工程改造當(dāng)中[16]。

        文中將結(jié)合筆者的部分研究結(jié)果,對(duì)各種酶約束模型的基本原理、構(gòu)建方法及其在生物學(xué)研究和代謝工程領(lǐng)域的應(yīng)用情況進(jìn)行詳細(xì)介紹。

        1 整合動(dòng)力學(xué)信息對(duì)酶濃度水平進(jìn)行限制

        2007年Beg等基于大分子擁擠的原理開(kāi)發(fā)了FBAwMC (Flux balance analysis with molecular crowding) 的酶約束方法[17]。考慮到細(xì)胞的總體積是有限的,在有限的細(xì)胞質(zhì)體積內(nèi)大分子(催化各種反應(yīng)的酶) 的濃度受到物理和空間上的一個(gè)限制[18-20],即在細(xì)胞質(zhì)中酶的摩爾體積(v) 與摩爾數(shù)(n)滿(mǎn)足

        引入一個(gè)擁擠系數(shù)i,其中

        對(duì)于每一個(gè)酶都可以根據(jù)底物濃度、m值、cat值、酶的分子量、酶的質(zhì)量體積和細(xì)胞質(zhì)密度這6個(gè)參數(shù)求解出擁擠系數(shù)a,它反映了第個(gè)酶對(duì)總體積擁擠程度的貢獻(xiàn)。但是由于底物濃度和m值等實(shí)驗(yàn)數(shù)據(jù)量有限,因此FBAwMC方法在進(jìn)行模擬計(jì)算時(shí)通過(guò)使實(shí)驗(yàn)測(cè)量值與模擬值之間的方差最小,選取了一個(gè)統(tǒng)一的擁擠系數(shù)值=(0.004±0.000 5) h?g/mmol。這種在物理空間水平上對(duì)酶濃度進(jìn)行限制的方法能夠通過(guò)求解得到催化每個(gè)反應(yīng)的酶所占的空間體積以及每一個(gè)酶的濃度信息。

        在Beg等開(kāi)發(fā)出第一個(gè)在酶濃度水平上對(duì)代謝網(wǎng)絡(luò)模型進(jìn)行約束的方法后,人們又基于相似思路提出了更多類(lèi)似模型擴(kuò)展方法。例如2011年Shlomi等為了研究人類(lèi)癌細(xì)胞中Warburg effect現(xiàn)象產(chǎn)生的原因而開(kāi)發(fā)的方法[23],以及其在2012年為了探究微生物的實(shí)際生長(zhǎng)速率針對(duì)開(kāi)發(fā)的MOMENT (Metabolic modeling with enzyme kinetics)算法[24],都是在FBAwMC方法的基礎(chǔ)上進(jìn)行的改良和擴(kuò)展。與FBAwMC不同的是這兩種方法沒(méi)有在空間體積上對(duì)酶的濃度進(jìn)行限制,而是直接為細(xì)胞內(nèi)酶的總量施加了濃度上限,使細(xì)胞代謝時(shí)的總酶需求量不超過(guò)該設(shè)定值(,單位g/g DW)。根據(jù)反應(yīng)的最大速率與催化反應(yīng)的酶濃度之間的關(guān)系得到約束方程,

        兩種方法在酶約束原理上基本一致,不同的是在計(jì)算酶的總濃度時(shí),MOMENT直接將細(xì)胞內(nèi)蛋白質(zhì)的質(zhì)量分?jǐn)?shù)ptot作為酶濃度上限 (中為0.56 g蛋白/g DW[25]),對(duì)于微生物來(lái)說(shuō)這顯然是一個(gè)偏高的值,因?yàn)榧?xì)胞內(nèi)的蛋白不僅包括催化反應(yīng)的酶還包括一些結(jié)構(gòu)蛋白、轉(zhuǎn)運(yùn)蛋白等。而在對(duì)人類(lèi)癌細(xì)胞內(nèi)的Warburg effect現(xiàn)象進(jìn)行研究時(shí),Shlomi等還基于蛋白質(zhì)組學(xué)數(shù)據(jù)計(jì)算出了模型中催化代謝反應(yīng)的酶占細(xì)胞內(nèi)總蛋白質(zhì)的質(zhì)量分?jǐn)?shù)為0.1 g酶/g蛋白[26],在此基礎(chǔ)上對(duì)細(xì)胞內(nèi)的總酶濃度作了進(jìn)一步限制。兩種方法的共同問(wèn)題是假定酶均處于底物飽和狀態(tài),但實(shí)際細(xì)胞內(nèi)很多酶的底物濃度都未達(dá)到飽和濃度,因此這些方法計(jì)算得到的酶濃度值都偏低。

        2017年,Sánchez等提出GECKO[27]方法成功地構(gòu)建了帶有酶約束的釀酒酵母基因組尺度代謝網(wǎng)絡(luò)模型。與上述提到的兩種方法不同的是,在計(jì)算總酶量時(shí)考慮到了酶的平均飽和度σ[28-29]的影響(表1)。2019年Massaiu等利用GECKO的方法構(gòu)建出了枯草桿菌的酶約束模型ec_iYO844[16],該模型中僅運(yùn)用了17個(gè)中心代謝反應(yīng)的cat值,嚴(yán)格地講,并不是基因組尺度的酶約束模型,但運(yùn)用該模型進(jìn)行與中心代謝途徑相關(guān)的表型預(yù)測(cè)時(shí),相對(duì)于普通模型iYO844[30]的預(yù)測(cè)結(jié)果有了顯著提高。

        在上述酶約束模型構(gòu)建時(shí)均需要酶的cat值來(lái)確定滿(mǎn)足酶約束的最優(yōu)酶濃度分布。目前雖然有很多數(shù)據(jù)庫(kù)如BRENDA[31]、SABIO-RK[32]、Metacyc[33]等都已經(jīng)收集了很多酶的相關(guān)動(dòng)力學(xué)參數(shù),但是仍有大量酶的cat數(shù)據(jù)缺失,對(duì)同一個(gè)酶不同文獻(xiàn)報(bào)道的結(jié)果也常存在較大差異[34]。因此酶cat值的獲取和可靠性是影響這類(lèi)模型模擬預(yù)測(cè)結(jié)果準(zhǔn)確性的關(guān)鍵。最近Palsson課題組針對(duì)cat值缺失的問(wèn)題開(kāi)發(fā)了一種機(jī)器學(xué)習(xí)方法,基于酶結(jié)構(gòu)、通量取樣信息等成功地對(duì)cat值進(jìn)行預(yù)測(cè)[35],為解決這一問(wèn)題提供了新的思路。

        除了對(duì)FBAwMC擴(kuò)展的方法之外,人們還開(kāi)發(fā)了一些其他將動(dòng)力學(xué)的相關(guān)信息引入GEM來(lái)增加酶約束的方法。例如2010年Yizhak等開(kāi)發(fā)的IOMA算法[36]整合了代謝組學(xué)和蛋白質(zhì)組學(xué)數(shù)據(jù),通過(guò)這些測(cè)得數(shù)據(jù)對(duì)代謝網(wǎng)絡(luò)模型中部分反應(yīng)的通量進(jìn)行限制,使實(shí)驗(yàn)測(cè)量值與預(yù)測(cè)值之間的誤差平方和最小來(lái)求得代謝通量分布。它首先從經(jīng)典的米氏方程出發(fā)推導(dǎo)出反應(yīng)通量與代謝物濃度和酶濃度之間的關(guān)系,得到:

        其中e表示參考狀態(tài)下的酶濃度,e表示實(shí)驗(yàn)測(cè)量的酶濃度。使核心反應(yīng)(可以得到代謝物組學(xué)和蛋白質(zhì)組學(xué)數(shù)據(jù)的反應(yīng)) 的模擬通量盡可能與通過(guò)公式(8) 計(jì)算出來(lái)的通量一致,引入一個(gè)可變參數(shù),約束方程變?yōu)椋?/p>

        V表示模擬得到的通量,優(yōu)化目標(biāo)是使的誤差平方和最小。

        IOMA類(lèi)似于MOMA[37]的思想,以通量變化最小為優(yōu)化目標(biāo),可以預(yù)測(cè)在不同因素?cái)_動(dòng)下通量的改變,并且預(yù)測(cè)結(jié)果顯著優(yōu)于MOMA。但這種方法依賴(lài)大量的實(shí)驗(yàn)數(shù)據(jù)包括代謝物濃度、酶濃度和m值等,對(duì)基因組規(guī)模模型應(yīng)用時(shí)會(huì)因?yàn)榇罅繑?shù)據(jù)缺失而無(wú)法模擬。

        隨著生長(zhǎng)速率的變化,另外3個(gè)部分的蛋白成本會(huì)相應(yīng)調(diào)整以滿(mǎn)足方程(12)的約束,進(jìn)而引起碳源攝入的再調(diào)整。

        與CAFBA原理相似,專(zhuān)門(mén)針對(duì)的溢流代謝研究開(kāi)發(fā)的PAT (Proteome allocation theory) 方法[42]將細(xì)胞內(nèi)總酶量分為發(fā)酵相關(guān)的酶、呼吸相關(guān)的酶以及用于細(xì)胞其他活動(dòng)的酶這3部分,根據(jù)3部分的總和為1得出約束方程。這個(gè)方程只適合專(zhuān)注于對(duì)微生物的溢流代謝現(xiàn)象進(jìn)行分析,在對(duì)其他代謝現(xiàn)象分析時(shí)存在一定的局限性。

        2 整合動(dòng)力學(xué)與熱力學(xué)信息對(duì)酶濃度水平進(jìn)行限制

        前面提到的方法都是從動(dòng)力學(xué)角度出發(fā)經(jīng)過(guò)公式推導(dǎo)得出酶限制的方程,接下來(lái)將主要介紹結(jié)合熱力學(xué)和動(dòng)力學(xué)經(jīng)典公式來(lái)對(duì)GEM添加酶約束的方法。

        最早將熱力學(xué)信息引入GEM的思想是2004年Holzhütter[43]提出的,將熱力學(xué)平衡常數(shù)K作為可逆反應(yīng)通量的權(quán)重因子來(lái)求出通量和最小的解。之后就出現(xiàn)了將熱力學(xué)與動(dòng)力學(xué)信息結(jié)合來(lái)對(duì)代謝反應(yīng)的酶成本進(jìn)行表征進(jìn)而約束代謝網(wǎng)絡(luò)模型的方法。2013年開(kāi)發(fā)的mTOW (Metabolic tug-of-war)[44]能夠求解出穩(wěn)態(tài)時(shí)酶效率最高、總代謝物負(fù)載最小的途徑。在估計(jì)酶成本時(shí),引入了一個(gè)懲罰項(xiàng),假定需要維持一定反應(yīng)通量的細(xì)胞必須通過(guò)更高的酶成本來(lái)補(bǔ)償較小的熱力學(xué)驅(qū)動(dòng)力,熱力學(xué)驅(qū)動(dòng)力越小酶成本就會(huì)越高。由可逆速率定律=(+––) (表示催化反應(yīng)的酶濃度,+和–分別表示正向和逆向的反應(yīng)速率)出發(fā),得到酶成本的函數(shù)關(guān)系(方程14)。當(dāng)反應(yīng)的熱力學(xué)驅(qū)動(dòng)力小于下限的閾值時(shí),酶成本無(wú)窮大;當(dāng)大于一定閾值時(shí),酶成本可以忽略不計(jì);在閾值范圍內(nèi)酶成本是反應(yīng)速率與代謝物濃度的函數(shù)。

        代謝物濃度通過(guò)以下約束求出,其在滿(mǎn)足經(jīng)驗(yàn)值[45](17)約束的同時(shí)還需使反應(yīng)的吉布斯自由能變化大于閾值。

        mTOW可以求解出反應(yīng)的吉布斯自由能變化、代謝物濃度以及酶濃度信息,是第一個(gè)不需要利用動(dòng)力學(xué)模型就可以得到代謝物濃度的方法。在求解上,得到的是代謝物濃度與酶濃度之間的一個(gè)折中,它沒(méi)有考慮動(dòng)力學(xué)參數(shù)對(duì)酶濃度的影響,因此不需要?jiǎng)恿W(xué)參數(shù)的輸入。雖然對(duì)酶濃度的估計(jì)是一個(gè)近似值,但模擬生長(zhǎng)在以葡萄糖為碳源的基本鹽培養(yǎng)基上的酶濃度時(shí),預(yù)測(cè)結(jié)果與實(shí)驗(yàn)結(jié)果之間仍有比較好的相關(guān)性。

        與mTOW同期,F(xiàn)lamholz等[46]提出的方法可以求出反應(yīng)在滿(mǎn)足熱力學(xué)驅(qū)動(dòng)力為正值下的最小酶成本。類(lèi)似于mTOW,在求解最小酶成本時(shí),F(xiàn)lamholz的方法限制了熱力學(xué)驅(qū)動(dòng)力必須為正值,以保證反應(yīng)可以正向進(jìn)行。在對(duì)酶成本進(jìn)行表征時(shí),首先從米氏方程出發(fā)結(jié)合Haldane關(guān)系推導(dǎo)出一個(gè)簡(jiǎn)化的適合于任意反應(yīng)的速率方程。在推導(dǎo)時(shí)假設(shè)了酶全部與底物結(jié)合,對(duì)產(chǎn)物飽和度的影響忽略不計(jì),因此計(jì)算得到的酶濃度值偏低。對(duì)這個(gè)方法進(jìn)行完善和優(yōu)化以后的ECM方法[47]將cat、代謝物濃度、吉布斯自由能變化、飽和度相關(guān)因素以及與調(diào)控有關(guān)的因素等作為反應(yīng)通量的權(quán)重因子來(lái)對(duì)酶成本進(jìn)行評(píng)估,這個(gè)方法可以自由地選擇將哪些因素考慮到酶成 本的計(jì)算當(dāng)中。同樣的,考察的因素越多需要的實(shí) 驗(yàn)參數(shù)就會(huì)越多,但結(jié)果就會(huì)更加接近于真實(shí)情況。

        corsoFBA[48]也是求解的FBA解空間下酶成本最小的途徑,通過(guò)熱力學(xué)驅(qū)動(dòng)力來(lái)對(duì)酶成本進(jìn)行罰分,驅(qū)動(dòng)力越大酶成本越小。這種方法可以將熱力學(xué)的限制整合入模型中,計(jì)算出途徑的酶成本大小。但它在表征酶成本時(shí),僅定性地考慮了分子量和熱力學(xué)驅(qū)動(dòng)力這兩個(gè)因素的影響,因此計(jì)算結(jié)果僅能用于途徑之間酶成本的評(píng)價(jià)排序,數(shù)值本身并不代表真實(shí)的酶用量。各種酶約束模型的約束方程及相關(guān)信息見(jiàn)表1和表2。

        表2 縮略語(yǔ)與名稱(chēng)的對(duì)應(yīng)信息

        3 酶約束模型在生物表型分析和代謝工程中的應(yīng)用

        目前構(gòu)建出來(lái)的酶約束模型能夠顯著提高對(duì)細(xì)胞表型的預(yù)測(cè)能力,通過(guò)模擬預(yù)測(cè)各種生物學(xué)現(xiàn)象揭示微生物代謝策略,在基礎(chǔ)生物學(xué)研究和指導(dǎo)代謝工程改造當(dāng)中得到廣泛的應(yīng)用。

        3.1 模擬微生物在單一碳源下的生長(zhǎng)情況

        Beg等和Adadi等分別利用FBAwMC、MOMENT對(duì)在10種和24種不同單一碳源下的最大生長(zhǎng)速率進(jìn)行了預(yù)測(cè),結(jié)果顯示與實(shí)驗(yàn)測(cè)量值顯著相關(guān),其中MOMENT預(yù)測(cè)的最大生長(zhǎng)速率與實(shí)驗(yàn)測(cè)量值的皮爾遜相關(guān)系數(shù)為0.47 (值為0.02)[17,24]。通過(guò)GECKO方法構(gòu)建的的酶約束模型ecYeast7顯著降低了模型中60%的反應(yīng)的通量可變性,也成功預(yù)測(cè)了其在3種不同培養(yǎng)基、12種不同碳源下的最大生長(zhǎng)速率,同時(shí)還與未考慮酶約束的普通基因組尺度代謝網(wǎng)絡(luò)模型Yeast7的預(yù)測(cè)值進(jìn)行比較,ecYeast7預(yù)測(cè)結(jié)果的平均相對(duì)偏差僅為8%,而Yeast7預(yù)測(cè)值的平均相對(duì)偏差高達(dá)100%。證明了酶約束模型能夠在不需要人為設(shè)定底物輸入速率的情況下很好地對(duì)微生物的最大生長(zhǎng)速率進(jìn)行預(yù)測(cè)[27]。

        Flamholz等利用酶成本最小化的方法對(duì)EMP和ED途徑間的選擇進(jìn)行分析,揭示了微生物生長(zhǎng)是在能量和酶成本之間進(jìn)行權(quán)衡之后選擇了不同的糖酵解策略。EMP途徑能量產(chǎn)生效率高但酶成本也高,當(dāng)酶成為限制因素時(shí)就會(huì)切換到能量產(chǎn)生效率低但酶成本也低的ED途徑[46]。

        3.2 模擬并研究微生物體內(nèi)的代謝溢流(Overflow) 現(xiàn)象

        代謝溢流現(xiàn)象普遍存在于微生物體內(nèi)(圖1),針對(duì)普通基因組尺度代謝網(wǎng)絡(luò)模型不能預(yù)測(cè)該現(xiàn)象的問(wèn)題,Sánchez利用構(gòu)建的ecYeast7模型成功模擬出在無(wú)氧和有氧條件時(shí)高生長(zhǎng)速率和高壓生長(zhǎng)狀態(tài)下體內(nèi)產(chǎn)生的Crabtree effect[27],展示了釀酒酵母體內(nèi)存在的高效率代謝途徑(即經(jīng)過(guò)TCA和呼吸鏈產(chǎn)能) 與低效率代謝途徑(即經(jīng)過(guò)糖酵解后溢出發(fā)酵產(chǎn)物) 間的切換。但利用普通的Yeast7模型則無(wú)法模擬出釀酒酵母體內(nèi)的這種代謝溢流現(xiàn)象。

        圖1 不同菌株體內(nèi)存在的代謝溢流現(xiàn)象

        Vazquez等利用已經(jīng)開(kāi)發(fā)出來(lái)的FBAwMC方法也解釋了細(xì)胞內(nèi)酶濃度的限制是代謝溢流產(chǎn)生的主要原因[49],與目前對(duì)于代謝溢流現(xiàn)象的研究結(jié)論相統(tǒng)一。van Hoek還利用分子擁擠的方法探究了不同微生物采取不同溢流策略的現(xiàn)象(酵母在開(kāi)始利用低效率的溢流代謝途徑后高效率代謝途徑的比例隨之下降,而大腸桿菌在出現(xiàn)溢流時(shí)高效率途徑仍維持較高的水平),發(fā)現(xiàn)維持氧化還原平衡在微生物途徑調(diào)整策略中發(fā)揮關(guān)鍵作用[50]。

        今年Mori等利用CAFBA的方法研究了在生長(zhǎng)和酶成本之間進(jìn)行權(quán)衡之后選擇的呼吸和發(fā)酵兩種代謝行為。通過(guò)將攝入的碳源通量分為生成能量的碳和生成生物質(zhì)前體的碳兩部分,證明了ATP合成的效率是微生物在細(xì)胞得率和酶成本之間進(jìn)行權(quán)衡的關(guān)鍵驅(qū)動(dòng)力[51]。

        3.3 探究微生物的底物層級(jí)利用情況

        Beg等利用FBAwMC評(píng)估了酶濃度限制對(duì)細(xì)胞代謝的作用,檢測(cè)了在混合碳源限制培養(yǎng)基中的底物利用情況。FBAwMC模擬結(jié)果顯示在混合碳源下微生物會(huì)經(jīng)過(guò)首先利用葡萄糖,再同時(shí)利用半乳糖、乳酸鹽和麥芽糖,最后利用甘油和乙酸這3個(gè)階段[17],與實(shí)驗(yàn)觀(guān)察到的結(jié)果一致。但僅考慮計(jì)量學(xué)約束而不考慮其他約束的傳統(tǒng)FBA模擬結(jié)果表示菌體會(huì)同時(shí)利用所有的碳源,這與實(shí)驗(yàn)結(jié)果不相符。最近湯超研究組也利用粗粒模型對(duì)微生物在多碳源下的底物利用策略進(jìn)行模擬[15],但與酶約束模型相比粗粒模型并不能給出特定蛋白的濃度。由此可見(jiàn),基于酶約束的模型可更好地用來(lái)研究微生物的底物層級(jí)利用情況(圖2)。

        圖2 E. coli在23種單一碳源下單位質(zhì)量生物質(zhì)合成需要的酶成本(以核糖作為碳源時(shí)的酶成本進(jìn)行歸一化)

        3.4 酶約束模型在代謝工程中的應(yīng)用

        除了基礎(chǔ)的生物學(xué)現(xiàn)象和機(jī)理研究以外,酶約束模型還能夠成功地進(jìn)行靶點(diǎn)預(yù)測(cè)指導(dǎo)代謝工程改造。Sánchez等利用ecYeast7通過(guò)通量控制分析(FCC)[52]的方法找到了釀酒酵母產(chǎn)富馬酸的關(guān)鍵酶作為改造靶點(diǎn)[27]。

        Massaiu等通過(guò)利用MOMA算法對(duì)iYO844和ec_iYO844兩個(gè)模型進(jìn)行單基因和雙基因敲除,來(lái)預(yù)測(cè)提高多聚谷氨酸產(chǎn)量的靶點(diǎn)基因。兩個(gè)模型對(duì)于單基因敲除的預(yù)測(cè)結(jié)果相一致,雙基因敲除結(jié)果存在差異。后對(duì)雙基因敲除靶點(diǎn)的預(yù)測(cè)結(jié)果進(jìn)行驗(yàn)證,結(jié)果表明,的酶約束模型ec_iYO844預(yù)測(cè)的雙基因敲除靶點(diǎn)的敲除能夠明顯提高目標(biāo)產(chǎn)品多聚谷氨酸的產(chǎn)量,而按照傳統(tǒng)的iYO844模型預(yù)測(cè)結(jié)果進(jìn)行的雙基因敲除卻未能達(dá)到預(yù)期結(jié)果。進(jìn)一步證明基于酶約束的模型在指導(dǎo)代謝工程改造中相較于普通模型具有更廣泛的應(yīng)用前景[16]。

        雖然目前酶約束模型在指導(dǎo)代謝工程改造中的應(yīng)用實(shí)例較少,但其構(gòu)建原理和方法具有很好的生物學(xué)意義及可行性,可以預(yù)見(jiàn)這類(lèi)模型具備更為精確地指導(dǎo)代謝工程改造的潛力。如根據(jù)催化每一步反應(yīng)所需的酶量來(lái)確定限速酶,然后理性選擇過(guò)表達(dá)或酶改造等實(shí)驗(yàn)策略來(lái)優(yōu)化途徑,提高目標(biāo)產(chǎn)品得率。另外,此類(lèi)模型模擬的理論上非最優(yōu)但更接近于實(shí)際的途徑還可用于指導(dǎo)途徑中副產(chǎn)物的去除等,對(duì)于工業(yè)生產(chǎn)也具有重要的意義。此外酶約束模型還可以預(yù)測(cè)在環(huán)境擾動(dòng)下(如高溫條件),細(xì)胞內(nèi)發(fā)生的代謝行為變化,進(jìn)而指導(dǎo)過(guò)程控制和菌株優(yōu)化[27]。

        4 總結(jié)與展望

        本文主要總結(jié)了近年來(lái)開(kāi)發(fā)的各種基于酶約束的代謝網(wǎng)絡(luò)模型構(gòu)建方法,這些方法結(jié)合動(dòng)力學(xué)中的米氏方程或熱力學(xué)第二定律通過(guò)引入cat、m、飽和度、eq、吉布斯自由能以及代謝組學(xué)和蛋白質(zhì)組學(xué)數(shù)據(jù)來(lái)對(duì)反應(yīng)消耗的酶成本進(jìn)行表征,進(jìn)而通過(guò)適當(dāng)?shù)募s束條件或優(yōu)化目標(biāo)求出滿(mǎn)足酶總量約束或者最小化酶成本的代謝通量分布及相應(yīng)酶量分布等信息。經(jīng)研究發(fā)現(xiàn),這種整合了酶約束的代謝網(wǎng)絡(luò)模型在表型預(yù)測(cè)上顯著優(yōu)于普通FBA模型,能夠解決目前普通基因組尺度代謝網(wǎng)絡(luò)模型不能解釋的一些生物學(xué)現(xiàn)象,可以提供更準(zhǔn)確可靠的代謝工程改造指導(dǎo)。雖然目前這些方法還存在一些不足和可以改進(jìn)的地方,但其為開(kāi)發(fā)其他以酶約束為基礎(chǔ)的新的約束方法奠定了堅(jiān)實(shí)的基礎(chǔ)。

        首先,在給酶成本添加總濃度限制的方法中,考慮的是整個(gè)細(xì)胞內(nèi)代謝酶的總量,這個(gè)值只適合用于對(duì)微生物的生長(zhǎng)進(jìn)行預(yù)測(cè),當(dāng)對(duì)某一產(chǎn)品的合成途徑預(yù)測(cè)時(shí)這往往是一個(gè)過(guò)高的酶濃度值(因可用于某一特定合成途徑的酶占總蛋白量更小)。對(duì)此可以通過(guò)將總酶量的限制拆分成對(duì)幾個(gè)不同功能相關(guān)的酶的量限制來(lái)進(jìn)一步細(xì)化,這實(shí)際上反映了細(xì)胞在基因調(diào)控水平的一種約束。

        其次,在同工酶、多功能酶的處理上,目前的優(yōu)化結(jié)果都是傾向于利用其中酶活最高的酶,而在真實(shí)細(xì)胞中具有不同比酶活的同工酶會(huì)在不同條件下表達(dá),對(duì)此酶約束模型并不能很好模擬??梢钥紤]以酶約束模型作為基礎(chǔ),引入某種形式的調(diào)控約束,通過(guò)兩者結(jié)合開(kāi)發(fā)出能模擬細(xì)胞內(nèi)同工酶和多功能酶表達(dá)策略的方法,以實(shí)現(xiàn)更精確的細(xì)胞行為模擬。

        最后,這些不同的方法都僅僅是涉及到細(xì)胞內(nèi)影響酶成本的部分因素。對(duì)相同質(zhì)量的酶可能會(huì)因?yàn)榘被峤M成不同而具有不同合成成本,如何以酶約束方法為基礎(chǔ),再將這些氨基酸合成成本的約束因素引入到模型中開(kāi)發(fā)新的方法,也是這些酶約束模型進(jìn)一步完善的方向。

        [1] Orth JD, Thiele I, Palsson B?. What is flux balance analysis? Nat Biotechnol, 2010, 28(3): 245–248.

        [2] Mahadevan R, Schilling CH. The effects of alternate optimal solutions in constraint-based genome-scale metabolic models. Metab Eng, 2003, 5(4): 264–276.

        [3] Gudmundsson S, Thiele I. Computationally efficient flux variability analysis. BMC Bioinformatics, 2010, 11: 489.

        [4] Yuan QQ, Li FR, Luo H, et al. Discovery of new strain modification strategies by metabolic network analysis. Chem Ind Eng Prog, 2017, 36(12): 4592–4600 (in Chinese). 袁倩倩, 李斐然, 羅浩, 等. 由代謝網(wǎng)絡(luò)分析發(fā)現(xiàn)菌種代謝工程改造新策略. 化工進(jìn)展, 2017, 36(12): 4592–4600.

        [5] Hao T, Ma HW, Zhao XM. Progress in automatic reconstruction and analysis tools of genome-scale metabolic network. Chin J Biotech, 2012, 28(6): 661?670 (in Chinese).郝彤, 馬紅武, 趙學(xué)明. 基因組尺度代謝網(wǎng)絡(luò)自動(dòng)重構(gòu)及分析工具研究進(jìn)展. 生物工程學(xué)報(bào), 2012, 28(6): 661–670.

        [6] Feng J, Gu YY, Quan YF, et al. Improved poly-γ-glutamic acid production inby modular pathway engineering. Metab Eng, 2015, 32: 106–115.

        [7] Schuetz R, Kuepfer L, Sauer U. Systematic evaluation of objective functions for predicting intracellular fluxes in. Mol Syst Biol, 2007, 3: 119.

        [8] Molenaar D, van Berlo R, de Ridder D, et al. Shifts in growth strategies reflect tradeoffs in cellular economics. Mol Syst Biol, 2009, 5: 323.

        [9] Varma A, Palsson BO. Stoichiometric flux balance models quantitatively predict growth and metabolic by-product secretion in wild-typeW3110. Appl Environ Microbiol, 1994, 60(10): 3724–3731.

        [10] Hsu PP, Sabatini DM. Cancer cell metabolism: Warburg and beyond. Cell, 2008, 134(5): 703–707.

        [11] Hermsen R, Okano H, You CH, et al. A growth-rate composition formula for the growth of.on co-utilized carbon substrates. Mol Syst Biol, 2015, 11(4): 801.

        [12] Narang A, Pilyugin SS. Bacterial gene regulation in diauxic and non-diauxic growth. J Theor Biol, 2007, 244(2): 326–348.

        [13] Basan M, Hui S, Okano H, et al. Overflow metabolism inresults from efficient proteome allocation. Nature, 2015, 528(7580): 99–104.

        [14] Berkhout J, Bosdriesz E, Nikerel E, et al. How biochemical constraints of cellular growth shape evolutionary adaptations in metabolism. Genetics, 2013, 194(2): 505–512.

        [15] Wang X, Xia K, Yang XJ, et al. Growth strategy of microbes on mixed carbon sources. Nat Commun, 2019, 10(1): 1279.

        [16] Massaiu I, Pasotti L, Sonnenschein N, et al. Integration of enzymatic data ingenome-scale metabolic model improves phenotype predictions and enablesdesign of poly-γ-glutamic acid production strains. Microb Cell Fact, 2019, 18(1): 3.

        [17] Beg QK, Vazquez A, Ernst J, et al. Intracellular crowding defines the mode and sequence of substrate uptake byand constrains its metabolic activity. Proc Natl Acad Sci USA, 2007, 104(31): 12663–12668.

        [18] Minton AP. Influence of macromolecular crowding upon the stability and state of association of proteins: predictions and observations. J Pharm Sci, 2005, 94(8): 1668–1675.

        [19] Hall D, Minton AP. Macromolecular crowding: qualitative and semiquantitative successes, quantitative challenges. Biochim Biophys Acta, 2003, 1649(2): 127–139.

        [20] Ellis RJ. Macromolecular crowding: obvious but underappreciated. Trends Biochem Sci, 2001, 26(10): 597–604.

        [21] Zimmerman SB, Trach SO. Estimation of macromolecule concentrations and excluded volume effects for the cytoplasm of. J Mol Biol, 1991, 222(3): 599–620.

        [22] Lee B. Calculation of volume fluctuation for globular protein models. Proc Natl Acad Sci USA, 1983, 80(2): 622–626.

        [23] Shlomi T, Benyamini T, Gottlieb E, et al. Genome-scale metabolic modeling elucidates the role of proliferative adaptation in causing the Warburg effect. PLoS Comput Biol, 2011, 7(3): e1002018.

        [24] Adadi R, Volkmer B, Milo R, et al. Prediction of microbial growth rate versus biomass yield by a metabolic network with kinetic parameters. PLoS Comput Biol, 2012, 8(7): e1002575.

        [25] Bremer H, Dennis PP. Modulation of chemical composition and other parameters of the cell by growth rate//Neidhardt FC, Curtiss III R, Ingraham JL, et al, Eds.and Salmonella: Cellular and Molecular Biology. 2nd ed. Washington, DC: American Society for Microbiology, 1996: 1553–1569.

        [26] Lee JK, Havaleshko DM, Cho H, et al. A strategy for predicting the chemosensitivity of human cancers and its application to drug discovery. Proc Natl Acad Sci USA, 2007, 104(32): 13086–13091.

        [27] Sánchez BJ, Zhang C, Nilsson A, et al. Improving the phenotype predictions of a yeast genome-scale metabolic model by incorporating enzymatic constraints. Mol Syst Biol, 2017, 13(8): 935.

        [28] Bakker BM, Bro C, K?tter P, et al. The mitochondrial alcohol dehydrogenase Adh3p is involved in a redox shuttle in. J Bacteriol, 2000, 182(17): 4730–4737.

        [29] van Hoek P, van Dijken JP, Pronk JT. Effect of specific growth rate on fermentative capacity of baker’s yeast. Appl Environ Microbiol, 1998, 64(11): 4226–4233.

        [30] Oh YK, Palsson BO, Park SM, et al. Genome-scale reconstruction of metabolic network inbased on high-throughput phenotyping and gene essentiality data. J Biol Chem, 2007, 282(39): 28791–28799.

        [31] Placzek S, Schomburg I, Chang A, et al. BRENDA in 2017: new perspectives and new tools in BRENDA. Nucleic Acids Res, 2017, 45(D1): D380–D388.

        [32] Wittig U, Rey M, Weidemann A, et al. SABIO-RK: an updated resource for manually curated biochemical reaction kinetics. Nucleic Acids Res, 2018, 46(D1): D656–D660.

        [33] Caspi R, Billington R, Ferrer L, et al. The MetaCyc database of metabolic pathways and enzymes and the BioCyc collection of pathway/genome databases. Nucleic Acids Res, 2016, 44(D1): D471–D480.

        [34] Bar-Even A, Noor E, Savir Y, et al. The moderately efficient enzyme: evolutionary and physicochemical trends shaping enzyme parameters. Biochemistry, 2011, 50(21): 4402–4410.

        [35] Heckmann D, Lloyd CJ, Mih N, et al. Machine learning applied to enzyme turnover numbers reveals protein structural correlates and improves metabolic models. Nat Commun, 2018, 9(1): 5252.

        [36] Yizhak K, Benyamini T, Liebermeister W, et al. Integrating quantitative proteomics and metabolomics with a genome-scale metabolic network model. Bioinformatics, 2010, 26(12): i255–i260.

        [37] Segrè D, Vitkup D, Church GM. Analysis of optimality in natural and perturbed metabolic networks. Proc Natl Acad Sci USA, 2002, 99(23): 15112–15117.

        [38] Mori M, Hwa T, Martin OC, et al. Constrained allocation flux balance analysis. PLoS Comput Biol, 2016, 12(6): e1004913.

        [39] Scott M, Gunderson CW, Mateescu EM, et al. Interdependence of cell growth and gene expression: origins and consequences. Science, 2010, 330(6007): 1099–1102.

        [40] Hui S, Silverman JM, Chen SS, et al. Quantitative proteomic analysis reveals a simple strategy of global resource allocation in bacteria. Mol Syst Biol, 2015, 11(2): 784.

        [41] You CH, Okano H, Hui S, et al. Coordination of bacterial proteome with metabolism by cyclic AMP signalling. Nature, 2013, 500(7462): 301–306.

        [42] Zeng H, Yang AD. Modelling overflow metabolism inwith flux balance analysis incorporating differential proteomic efficiencies of energy pathways. BMC Syst Biol, 2019, 13(1): 3.

        [43] Holzhütter HG. The principle of flux minimization and its application to estimate stationary fluxes in metabolic networks. Eur J Biochem, 2004, 271(14): 2905–2922.

        [44] Tepper N, Noor E, Amador-Noguez D, et al. Steady-state metabolite concentrations reflect a balance between maximizing enzyme efficiency and minimizing total metabolite load. PLoS ONE, 2013, 8(9): e75370.

        [45] Bennett BD, Kimball EH, Gao M, et al. Absolute metabolite concentrations and implied enzyme active site occupancy in. Nat Chem Biol, 2009, 5(8): 593–599.

        [46] Flamholz A, Noor E, Bar-Even A, et al. Glycolytic strategy as a tradeoff between energy yield and protein cost. Proc Natl Acad Sci USA, 2013, 110(24): 10039–10044.

        [47] Noor E, Flamholz A, Bar-Even A, et al. The protein cost of metabolic fluxes: prediction from enzymatic rate laws and cost minimization. PLoS Comput Biol, 2016, 12(11): e1005167.

        [48] Schultz A, Qutub AA. Predicting internal cell fluxes at sub-optimal growth. BMC Syst Biol, 2015, 9: 18.

        [49] Vazquez A, Oltvai ZN. Macromolecular crowding explains overflow metabolism in cells. Sci Rep, 2016, 6: 31007.

        [50] van Hoek MJ, Merks RM. Redox balance is key to explaining full. partial switching to low-yield metabolism. BMC Syst Biol, 2012, 6: 22.

        [51] Mori M, Marinari E, de Martino A. A yield-cost tradeoff governs’s decision between fermentation and respiration in carbon-limited growth. npj Syst Biol Appl, 2019, 5(1): 16.

        [52] Nilsson A, Nielsen J. Metabolic trade-offs in yeast are caused by F1F0-ATP synthase. Sci Rep, 2016, 6: 22264.

        Progress and application of metabolic network model based on enzyme constraints

        Xin Zhao1,2, Xue Yang1,2, Zhitao Mao1,2, and Hongwu Ma1

        1 Key Laboratory of Systems Microbial Biotechnology, Tianjin Institute of Industrial Biotechnology, Chinese Academy of Sciences, Tianjin 300308, China 2 University of Chinese Academy of Sciences, Beijing 100049, China

        Genome-scale metabolic network models have been successfully applied to guide metabolic engineering. However, the conventional flux balance analysis only considers stoichiometry and reaction direction constraints, and the simulation results cannot accurately describe certain phenomena such as overflow metabolism and diauxie growth on two substrates. Recently, researchers proposed new constraint-based methods to simulate the cellular behavior under different conditions more precisely by introducing new constraints such as limited enzyme content and thermodynamics feasibility. Here we review several enzyme-constrained models, giving a comprehensive introduction on the biological basis and mathematical representation for the enzyme constraint, the optimization function, the impact on the calculated flux distribution and their application in identification of metabolic engineering targets. The main problems in these existing methods and the perspectives on this emerging research field are also discussed. By introducing new constraints, metabolic network models can simulate and predict cellular behavior under various environmental and genetic perturbations more accurately, and thus can provide more reliable guidance to strain engineering.

        enzyme constraints, metabolic network model,metabolic engineering, enzyme kinetics, thermodynamics

        10.13345/j.cjb.190220

        馬紅武 中國(guó)科學(xué)院天津工業(yè)生物技術(shù)研究所“百人計(jì)劃”研究員,中科院系統(tǒng)微生物技術(shù)重點(diǎn)實(shí)驗(yàn)室副主任。2001年在天津大學(xué)獲得生物化工博士學(xué)位。后分別在德國(guó)和英國(guó)從事研究工作。2011年底回國(guó)任中科院天津工業(yè)生物技術(shù)研究所研究員。主要研究方向包括代謝網(wǎng)絡(luò)分析和途徑設(shè)計(jì)、代謝工程、計(jì)算生物學(xué)軟件開(kāi)發(fā)等,在相關(guān)領(lǐng)域取得了在國(guó)際上具有重要影響的研究成果,在、等刊物發(fā)表學(xué)術(shù)論文40余篇,SCI引用千余次,2014年獲湯森路透高被引科學(xué)家獎(jiǎng)。兼任中國(guó)生物工程學(xué)會(huì)生物資源專(zhuān)業(yè)委員會(huì)副主任委員。

        趙欣, 楊雪, 毛志濤, 等. 基于酶約束的代謝網(wǎng)絡(luò)模型研究進(jìn)展及其應(yīng)用. 生物工程學(xué)報(bào), 2019, 35(10): 1914–1924.

        Zhao X, Yang X, Mao ZT, et al. Progress and application of metabolic network model based on enzyme constraints. Chin J Biotech, 2019, 35(10): 1914–1924.

        May 29, 2019;

        July 15, 2019

        Supported by: International Partnership Program of Chinese Academy of Sciences (No. 153D31KYSB20170121), National Key Research and Development Program of China (No. 2018YFA0900301).

        Hongwu Ma. Tel: +86-22-24828735; E-mail: ma_hw@tib.cas.cn

        中國(guó)科學(xué)院國(guó)際合作局對(duì)外合作重點(diǎn)項(xiàng)目(No. 153D31KYSB20170121),科技部重點(diǎn)專(zhuān)項(xiàng)(No. 2018YFA0900301) 資助。

        2019-07-23

        http://kns.cnki.net/kcms/detail/11.1998.Q.20190722.1442.004.html

        (本文責(zé)編 郝麗芳)

        猜你喜歡
        成本方法模型
        一半模型
        2021年最新酒駕成本清單
        河南電力(2021年5期)2021-05-29 02:10:00
        重要模型『一線(xiàn)三等角』
        重尾非線(xiàn)性自回歸模型自加權(quán)M-估計(jì)的漸近分布
        溫子仁,你還是適合拍小成本
        電影(2018年12期)2018-12-23 02:18:48
        可能是方法不對(duì)
        3D打印中的模型分割與打包
        用對(duì)方法才能瘦
        Coco薇(2016年2期)2016-03-22 02:42:52
        四大方法 教你不再“坐以待病”!
        Coco薇(2015年1期)2015-08-13 02:47:34
        捕魚(yú)
        天天综合久久| 国产精品无码dvd在线观看| 性高朝大尺度少妇大屁股| 国产又爽又黄的激情精品视频| 亚洲成AV人片无码不卡| 丝袜av乱码字幕三级人妻| 一本加勒比hezyo无码专区| 婷婷丁香五月中文字幕| 美女视频很黄很a免费国产| 成年人视频在线观看麻豆| 成人免费无码大片a毛片抽搐色欲| 亚洲18色成人网站www| 日韩国产精品无码一区二区三区| 亚洲学生妹高清av| 亚洲三级在线播放| 乳乱中文字幕熟女熟妇| 日本在线 | 中文| 在线观看av永久免费| 91久久久久无码精品露脸 | 亚洲人成电影在线观看天堂色| 乱码午夜-极品国产内射| 噜噜噜色97| 国产精品亚洲一级av第二区| 亚洲成在人线av品善网好看| 免费网站国产| 国产亚洲美女精品久久久2020| 91精品国产九色综合久久香蕉| 亚洲av无码国产精品草莓在线| 日韩精品大片在线观看| 国产一区二区三区色区| 国产一区二区三区视频网| 天天做天天摸天天爽天天爱| 亚洲色欲久久久久综合网 | 国产白嫩美女在线观看| 亚洲女同系列高清在线观看 | 久久综合亚洲鲁鲁五月天| 2019最新中文字幕在线观看| 好日子在线观看视频大全免费动漫| 99久久超碰中文字幕伊人| 男人一插就想射的原因 | 国产精品久久久久久亚洲av|