国产日韩欧美一区二区三区三州_亚洲少妇熟女av_久久久久亚洲av国产精品_波多野结衣网站一区二区_亚洲欧美色片在线91_国产亚洲精品精品国产优播av_日本一区二区三区波多野结衣 _久久国产av不卡

?

分子動力學探究高壓對β-乳球蛋白與多酚結合的影響

2022-02-23 13:33:36孔慶新羅麗梅黃業(yè)傳徐偉平張徑舟劉碧林任麗萍
食品與發(fā)酵工業(yè) 2022年3期
關鍵詞:殘基氫鍵位點

孔慶新,羅麗梅,黃業(yè)傳,徐偉平,張徑舟,劉碧林,任麗萍

1(重慶化工職業(yè)學院 制藥工程學院,重慶,401228) 2(荊楚理工學院 生物工程學院,荊門,448000) 3(西南政法大學醫(yī)院,重慶,401120)

牛乳是一種營養(yǎng)豐富、老少皆宜的食物,除直接飲用外,現也將一些植物原料如茶、果汁加入牛奶中,使產品形式和口味多樣、營養(yǎng)更全面。植物原料中大都含有多酚,其具有保健作用,但也易與牛乳中蛋白形成復合物,進而影響多酚的生物活性、牛乳蛋白的性質和產品品質。很多學者都對蛋白與多酚的結合進行了研究,如兩者間的結合類型、結合位點、結合機理等[1-4]。多酚與蛋白的交聯主要靠非共價鍵,如氫鍵、疏水相互作用、范德華力、靜電相互作用,在一定條件下也發(fā)生共價交聯[5-6]。乳制品的殺菌仍以熱殺菌為主,熱殺菌對產品品質有負面影響。超高壓是一種冷殺菌技術,其具有良好殺菌效果且?guī)缀跄軌蛲旰帽A羰澄镏行》肿訝I養(yǎng)和風味物質,因此在包括乳制品的食品中應用越來越廣[7-8]。超高壓應用于牛乳時會改變乳蛋白的結構,進而影響乳制品的品質。一些學者報道了超高壓處理對牛乳蛋白結構的影響,RUSSO等[9]表明壓力處理時β-乳球蛋白的結構較穩(wěn)定,僅300 MPa以上的壓力會使蛋白逐漸展開;CONSIDINE等[10]發(fā)現超高壓條件下β-乳球蛋白的變性或聚集主要是蛋白中二硫鍵重排引起的,小于100 MPa為天然構象階段,150~350 MPa為二硫鍵發(fā)生重排的單體及二聚體,高于500 MPa主要為多聚體。超高壓用于植物乳飲料處理后是否影響到多酚與蛋白的結合目前報道很少,CHEN等[11]報道了大豆蛋白與茶多酚在400 MPa高壓下主要靠疏水作用和氫鍵交聯,高壓下蛋白對茶多酚的生物活性有保護作用。

光譜技術、波譜技術、量熱技術、原子力顯微鏡等都被用于研究蛋白質結構,但測得的只是靜態(tài)結構;生物體內蛋白質一直處于動態(tài)變化,現分子動力學模擬已成為研究蛋白結構的一種有效手段,成為繼實驗和理論手段后從分子水平了解和認識蛋白質世界的第3種手段,是一種具有足夠小的時間尺度和空間尺度的模擬技術[12]。牛乳蛋白與許多小分子的結合都用分子動力學進行了模擬,如柑橘黃酮[13]、姜黃素[14]、芹菜素[15]、茶多酚[4]等。而超高壓下的分子模擬較少,一些學者對胰島素[16]和脂肪酶[12]在高壓下的結構變化進行了模擬。

β-乳球蛋白是乳清中含量最為豐富的一種蛋白質,由于其良好的功能和營養(yǎng)特性在食品工業(yè)中有廣泛的應用,其每個單體含有162個氨基酸,分子質量為180 kDa,是由8條反向平行的β-折疊構成的桶狀結構,其外側含有α-螺旋結構,對疏水及兩性配體均具有良好的親和力[17],而表沒食子兒茶素沒食子酸酯(epigallocatechin gallate,EGCG)是茶多酚中最主要的成分之一。因此本研究以β-乳球蛋白和EGCG為研究對象,用分子對接和分子模擬技術探討蛋白經0.1、200、500、800 MPa的壓力處理后,其與小分子結合能和結合機理的變化,以期為進一步提高茶乳飲料的質量和超高壓殺菌技術在乳飲料中的應用奠定一定的理論基礎。

1 實驗方法

1.1 β-乳球蛋白與茶多酚的分子對接

β-乳球蛋白數據從RCSB網站下載(ID號:3npo),用PyMOL除去結晶水。EGCG的結構用ChemDraw Professional 17.1和 Chem3D 17.1準備,并用自帶的MM2力場進行結構優(yōu)化,通過AutoDock Tools (1.5.6版本)分別對蛋白和小分子進行處理,其中小分子加氫、計算Gasteiger電荷,蛋白加氫、計算Gasteiger電荷、合并非極性氫,生成β-乳球蛋白和EGCG的pdbqt文件。采用盲對接,小分子設置為全柔性,盒子大小設定為40 ?×40 ?×40 ?,中心坐標為默認值,設置間隔為1 ?。對接過程中采用AutoDock Vina搜索前20個最優(yōu)的構象,對接參數energy_range設置為5,exhaustiveness設置為100。

1.2 不同壓力條件下EGCG與乳球蛋白結合的分子動力學模擬

以1.1中的最優(yōu)構象作為分子動力學模擬的初始結構,剝離出小分子的pdb結構,并在ATB網站生成小分子的top文件。分子動力學模擬采用Gromacs(2019.6)軟件,首先將小分子合并到蛋白中,修改復合物的top文件;模擬時選用GROMOS54a7力場,將蛋白放入立方體水盒子中,水模型采用SPC模型,使蛋白離盒子邊緣最短距離為1 nm并添加Na+使體系達到電中性,并使用最速下降法對體系進行能量最小化[13]。然后在NVT和NPT系綜下分別進行400 ps的平衡,平衡后溫度為300 K,壓力達到預先設定的壓力值。最后進行150 ns的分子動力學模擬,使用LINCS算法約束所有鍵,使用PME計算靜電作用,范德華相互作用使用截斷半徑為1.4 nm進行計算。選用Parrinello-Rahman壓浴方式,Isotropic為控壓方式,使用V-rescale的控溫方式。模擬步長為2 fs,每10 ps貯存1次數據,每個處理平行模擬2次。

模擬結束后,先去除周期性邊界條件,然后使用gmx的rms命令分析各壓力下小分子與蛋白的均方根誤差(root mean squared error,RMSD);并用gmx的相應命令分析各壓力下蛋白殘基的波動[均方根漲落(root mean square fluctuation,RMSF]、蛋白二級結構、可及表面積的變化;并用能量分析命令分析蛋白的體積變化。RMSD分析全部150 ns的數據,而其余指標分析結構穩(wěn)定后的數據,即100~150 ns。分析時每100 ps選取1個數據點,作圖用Origin 8.0軟件。

另外,將初始純蛋白(3npo)也按上述方法在不同壓力下進行150 ns的分子模擬,然后分別提取130、135、140、145、150 ns的蛋白結構,并按1.1的方法分別和小分子EGCG進行對接,以觀察高壓處理后小分子在蛋白的對接位點是否發(fā)生了變化。

1.3 MMPBSA計算小分子與蛋白的結合自由能

提取模擬過程中100~150 ns的軌跡(xtc文件),借助于g_mmpbsa工具,利用MMPBSA(molecular mechanics Poisson-Boltzmann surface area)方法計算小分子與蛋白之間的結合自由能[15]。

1.4 小分子與蛋白結合的二維平面圖和表面結構圖

提取各壓力條件下復合物100~150 ns的平均pdb結構,用LigPlot+(V 2.2)軟件分析小分子與蛋白在結合位點的疏水和氫鍵作用;并在PyMOL中繪制蛋白的表面結構圖及與小分子結合情況。

2 結果與分析

2.1 茶多酚在蛋白的結合位點

通過分子對接,發(fā)現EGCG在β-乳球蛋白前20個能量較高的結合位點主要分布在2個區(qū)域,如圖1所示,其中主要是結合在疏水腔上部的1號區(qū)域,另外還有位于蛋白表面的2號區(qū)域,在1號區(qū)域的平均對接能大于2號區(qū)域。關于多酚在β-乳球蛋白結合位點的報道較多,其中有研究認為在中性pH時小分子主要結合在疏水腔,而在酸性條件下則結合在蛋白表面[4,19],AL-SHABIBN等[20]也發(fā)現蘆丁既可以結合在疏水腔也可以結合在蛋白表面,其中疏水腔的結合力更大,這與本研究結論一致。另有學者發(fā)現小分子吸附在乳球蛋白的表面,未結合在疏水腔[1,21],這可能與小分子結構、對接的蛋白初始結構、模擬條件等差異有關。

蛋白在不同壓力條件下進行150 ns模擬后再進行分子對接,對接位點發(fā)生了明顯變化,500 MPa及以下的壓力處理后小分子仍結合到1、2號區(qū)域,其中常壓下(0.1 MPa)的前20個位點主要仍是位于1號區(qū)域;而200和500 MPa處理后則在2號區(qū)域為最優(yōu)區(qū)域,其結合能高于1號區(qū)域;800 MPa處理后,小分子的結合位點增多,除了1、2號區(qū)域外還增加了其他結合位點,如3、4和5號區(qū)域,沒有明顯的規(guī)律性,這種現象可能與蛋白的結構變化有關,800 MPa下蛋白結構發(fā)生變化后使可結合小分子的位點增加,或1、2號區(qū)域發(fā)生了不利于小分子結合的構象變化。

圖1 EGCG在β-乳球蛋白的主要結合位點Fig.1 Main binding sites of EGCG to β-lactoglobulin

2.2 EGCG與β-乳球蛋白結合的分子動力學

2.2.1 分子動力學模擬過程中RMSD的變化

RMSD是衡量特定時間蛋白結構與原始構象的平均偏差,是評價研究體系是否穩(wěn)定的重要指標。由圖2-a可知,500 MPa下小分子的波動較大,100 ns后才達到平衡,而其他壓力下基本在30 ns左右就達到了平衡,從平衡后(均取100 ns后的數據)的平均值來對比,常壓下為0.192 nm,與200 MPa下相當(0.194 nm),500 MPa(0.266 nm)時最高,800 MPa時為0.213 nm。就蛋白的RMSD而言(圖2-b),均是40 ns左右就達到了平衡,且整體波動很小,說明β-乳球蛋白與EGCG的結合過程對蛋白結構影響很小,兩者的結合總體平穩(wěn)。同樣,取100 ns后的平均值來衡量,常壓下為0.249 nm,200 MPa下蛋白的RMSD低于常壓(0.238 nm),隨著壓力的繼續(xù)增加,蛋白的RMSD也增加,且高于常壓,其中500和800 MPa下分別為0.272和0.259 nm。說明200 MPa的壓力處理能使蛋白結構更穩(wěn)定,而500和800 MPa的壓力處理能使β-乳球蛋白的穩(wěn)定性略為降低。SAHIHI等[13]報道柑橘黃酮與β-乳球蛋白間的結合很平穩(wěn),得到的RMSD值也較低;GHOLAMI等[21]研究柚皮苷和β-乳球蛋白的結合時發(fā)現柚皮苷的波動比蛋白小;ABDULATIF等[20]發(fā)現蘆丁與β-乳球蛋白結合時,小分子與蛋白的波動相當,本研究中小分子的波動略低于蛋白,與上面的報道基本一致。KURPIEWSKA等[16]報道胰島素在200或500 MPa處理時殘基位移會減小;HATA等[22]也報道高壓下蛋白的結構波動會減小,而本研究只發(fā)現200 MPa時蛋白波動小于常壓。

a-小分子;b-蛋白圖2 不同壓力處理模擬過程中小分子和蛋白的RMSD對比Fig.2 RMSD of small molecule and protein during different pressure simulations

2.2.2 分子動力學模擬過程中RMSF的變化

RMSF主要反映蛋白中各氨基酸殘基在不同壓力處理時的波動情況,波動越大,說明該殘基在壓力處理中柔性越大,平均柔性越高蛋白結構也就越不穩(wěn)定。從圖3可知,不同壓力處理體系中RMSF值相當,沒有明顯差異,說明高壓總體上對殘基波動影響不大。就小分子和蛋白結合口袋附近而言,殘基波動較大的有86號位的Ala,87號位的Leu和88號位的Asn。RMSF大小與模擬常壓下柑橘黃酮與β-乳球蛋白結合時得到的數據基本一致[13],但本研究中個別殘基的波動更大,可能是受高壓的影響。

圖3 分子模擬過程中不同壓力處理蛋白RMSF對比Fig.3 RMSF of different pressure treatment during molecular dynamics simulation

2.2.3 分子動力學模擬過程中蛋白體積的變化

如圖4所示,隨壓力的增加,蛋白體積顯著減少,且體積在150 ns的模擬過程中保持恒定。根據Le Chatelier原理,高壓下蛋白朝體積減小的方向進行,這可能是壓力下,水分子進入疏水腔使疏水腔變形,導致蛋白的二級和三級結構發(fā)生變化進而體積減小[22]。KURPIEWSKA等[16]也報道了胰島素在100 MPa下體積減少約1.4%。

圖4 不同壓力處理模擬過程中蛋白體積變化Fig.4 Volume comparison of different pressure treatment during molecular dynamics simulation

2.2.4 分子動力學模擬過程中蛋白溶劑可及表面積的變化

由表1可知,500 MPa以上的高壓處理可以明顯減少蛋白的溶劑可及表面積,其主要是由疏水表面積的減少所致,而親水表面積變化不大。這可能是由于高壓能導致蛋白總體積減少,從而使溶劑可及表面也減少,更多的疏水表面在壓力作用下被掩埋。但各壓力下親水表面積均大于疏水表面積,因此壓力作用不會改變蛋白的水溶性。疏水表面積的減少可能是高壓條件下水分子滲透進入蛋白分子內部使更多的疏水區(qū)域暴露在極性水溶液中[23]。其他一些研究者也報道了高壓下蛋白疏水性的變化,但結果并不一致,有的研究者認為隨壓力的增加蛋白的疏水作用減弱[11,22],而另一些學者報道高壓下蛋白的疏水性增加[24-25]。這些研究均是采用實驗手段研究表面疏水性,而本研究采用分子模擬研究疏水表面積的變化。

表1 不同壓力處理對蛋白溶劑可及表面積的影響 單位:nm2

2.2.5 分子動力學模擬過程中蛋白二級結構的變化

如表2所示,200 MPa時蛋白的二級結構變化不大,主要是彎曲減少而轉角增加;而500和800 MPa條件下,蛋白的二級結構發(fā)生了較大的變化,α-螺旋明顯減少,β-折疊也有所降低,而無規(guī)則卷曲顯著增加,500和800 MPa兩者間差異并不明顯。因此,壓力較低時(200 MPa)對β-乳球蛋白的二級結構影響不大,而壓力較高時(500 MPa以上)蛋白二級結構則變化較大。邱春江[24]研究高壓處理肌原纖維蛋白時也發(fā)現在高壓下蛋白的α-螺旋下降,α-螺旋對壓力比較敏感,可能是高壓下蛋白的氫鍵受到一定程度破壞,而氫鍵對α-螺旋的穩(wěn)定至關重要[26],本研究也發(fā)現500 MPa及以上壓力處理對α螺旋的破壞更大。白雨鑫[25]的研究表明500 MPa的壓力對β-乳球蛋白的二級結構影響很小,CHEN等[11]研究發(fā)現400 MPa的壓力處理能使大豆蛋白的α-螺旋減少而β-折疊增加,陳剛[12]報道脂肪酶在400 MPa處理時二級結構變化不大,而500 MPa以上時α-螺旋部分轉化為無規(guī)則卷曲或β-折疊,這些差異可能是由蛋白不同或測定方法不同引起,上述報道大都采用實驗手段測得而本研究是分子模擬。

表2 不同壓力處理對蛋白二級結構的影響Table 2 Effect of different pressure on the secondary structure of lactoglobulin

2.3 MMPBSA分析小分子與蛋白間的結合自由能

由表3可知,小分子和蛋白間的結合主要靠范德華力,靜電作用力和非極性溶劑化自由能也有一定的貢獻,而極性溶劑化自由能為正值,對蛋白的穩(wěn)定有負面作用。隨著壓力的不斷增加,靜電相互作用力有所增加,但范德華力相比常壓模擬體系均降低,另外極性溶劑化所消耗的能量也增加,而非極性溶劑化能基本不變。在以上幾種作用力的綜合作用下,總結合自由能隨壓力增加逐漸降低,因此在超高壓作用下,小分子與蛋白的結合不如常壓下穩(wěn)定。除了高壓對小分子結合的直接影響以外,這也可能與分子模擬只選擇了對接能最高的1號結合位點有關,在2.1的分析中可知,高壓處理后的蛋白更易結合在2號或其他位點,因此可以推斷,β-乳球蛋白高壓處理后在1號結合位點可能發(fā)生了不利于小分子結合的蛋白結構變化,具體通過后面的機制分析來闡明。

表3 不同壓力對蛋白與EGCG結合的MMPBSA自由能的影響 單位:kJ/mol

AL-SHABIBN等[20]研究蘆丁和β-乳球蛋白的結合時,用MMGBSA計算出兩者在蛋白表面的結合能為-68.011 2 kcal/mol,大于本研究的結果,這可能與小分子結構不同有關,也可能與計算方法不同有關,他們使用的是MMGBSA。ZHAN等[27]研究β-乳球蛋白與辣椒素結合時用MMPBSA計算得到的結合能是-142.744 kJ/mol,也是范德華作用為主要結合力,與本研究結果一致。陳剛[12]發(fā)現200 MPa處理時脂肪酶與乙酸乙酯的結合能增加,而400和600 MPa時結合能下降,其中600 MPa時下降更多,與本研究的結果稍有不一致。

從各個殘基對能量的貢獻來看,4個壓力下貢獻值均超過1 kJ/mol的有8個殘基,分別是Leu31、Leu39、Val41、Leu58、Ile71、Val92、Met107、Glu108,大部分為疏水氨基酸,說明蛋白與小分子的結合疏水作用也是一種主要結合力。特別值得注意的是,壓力增加時,Leu87、Met107和Ala118等殘基的貢獻顯著減弱,具體的機制有待于進一步研究。

2.4 蛋白與小分子可能的結合機制分析

由圖5-a可以看出0.1 MPa時小分子與蛋白的7個氨基酸存在疏水作用,分別是Leu39、Val41、Lys69、Ile71、Ala86、Met107、Glu108,小分子所處氨基酸環(huán)境與CHENG等[28]報道的矢車菊素-3-O-葡萄糖苷和KANAKIS等[17]報道的EGCG結合在β-乳球蛋白后所處的氨基酸環(huán)境基本一致。小分子與蛋白間形成6個氫鍵,分別與Asn88形成2個氫鍵、與Asn90間1個氫鍵、與Asn109間2個氫鍵,另外與Ser116間形成1個氫鍵。200 MPa時小分子所處的疏水環(huán)境和形成的氫鍵均和0.1 MPa沒有差異,500 MPa時所處的疏水環(huán)境和0.1 MPa相同,但形成的氫鍵少了1個,其中與Asn109間只形成了1個氫鍵(200和500 MPa條件下的二維平面圖未顯示)。而800 MPa時所處的疏水環(huán)境和形成的氫鍵均發(fā)生了比較明顯的變化,如圖5-b所示,雖然也是和周圍的7個氨基酸存在疏水相互作用,但氨基酸與前面3個壓力不同,其中多了Leu87和Asn90,而少了Lys69和Ala86;形成的氫鍵只有4個,分別與Lys69、Asn88、Asn109和Ser116形成氫鍵。

非配體鍵; 配體鍵; 氫鍵及長度; 參與疏水作用的非配體殘基; 參與疏水作用的相應原子 a-0.1 MPa;b-800 MPa圖5 0.1 和800 MPa下小分子與蛋白相互作用的二維平面圖Fig.5 The 2D plot for interaction between EGCG and β-lactoglobulin under 0.1 and 800 MPa注:B6e163在二維平面圖中代表小分子

由表2可知,隨著處理壓力的增加,小分子與蛋白間的靜電作用力也增加,可能的原因是除了氫鍵以外,蛋白與小分子間的其他弱靜電相互作用也很重要[29]。結合2.3的結果可推斷,壓力處理整體上對小分子與蛋白間的疏水作用影響不大,但800 MPa時小分子所處的疏水環(huán)境有微小改變;另外,隨著壓力的增加(500 MPa及以上),小分子與蛋白間的氫鍵會有一定程度的破壞,這可能是高壓下小分子與蛋白間結合能降低的重要原因之一。

賈晶晶[30]研究EGCG與β-乳球蛋白結合時發(fā)現,結合力以疏水作用為主,同時也有范德華力和氫鍵。CHENG等[28]報道矢車菊素-3-O-葡萄糖苷與β-乳球蛋白交聯時發(fā)現主要結合力為氫鍵和疏水作用,且形成的其中2個氫鍵是與Asn109和Ser116,與本研究基本一致。本研究中隨壓力增加,小分子與蛋白間結合能逐漸降低應該是氫鍵、范德華力和靜電相互作用力變化引起的綜合效果,其中范德華力和氫鍵均減弱,靜電作用有所增加;疏水作用雖然對兩者的結合也很重要,但500 MPa及以下的壓力沒有影響到疏水結合力和小分子的疏水環(huán)境。

為進一步了解高壓對蛋白表面結構及與小分子結合的影響,用PyMOL軟件繪制了蛋白分子表面結構圖,圖6-a和圖6-b分別為0.1和800 MPa下100~150 ns模擬過程中蛋白的平均表面結構及多酚與蛋白的結合情況,可以看出,蛋白在800 MPa時表面結構發(fā)生了顯著變化,特別是上部臨近疏水腔的口袋處,結合位點處蛋白結構的變化肯定會引起小分子結合姿勢與結合能的差異,表3的數據也證實了不同壓力下小分子與蛋白的結合能差異較大,說明壓力越大,小分子與蛋白的結合位點處發(fā)生了越不利于兩者結合的變化;前面的分子對接也表明壓力處理后小分子的最佳結合位點并不在1號位;另外,0.1 MPa時小分子結構更加伸展(圖6-a),能與蛋白更多殘基發(fā)生作用,從而結合力也更高;這也與前面的分析是一致的,0.1 MPa時小分子與蛋白間形成的氫鍵更多。

a-0.1 MPa;b-800 MPa圖6 0.1和800 MPa下蛋白分子表面結構及與多酚結合情況Fig.6 Surface structure and polyphenol binding of protein at 0.1 and 800 MPa

3 結論

不同壓力下茶多酚在β-乳球蛋白上的結合位點不一樣,其中常壓下主要結合在疏水腔的1號位點,蛋白在200和500 MPa處理后主要結合在蛋白表面的2號位點,而800 MPa處理后除了1號和2號外,還增加了其他結合位點,這可能與高壓下蛋白體積的收縮和蛋白變性引起的蛋白結構變化有關。以1號結合位點為目標進行分子動力學模擬,發(fā)現在150 ns的模擬過程中,200 MPa的壓力處理能使小分子和蛋白的波動都較小,而500和800 MPa下小分子和蛋白的波動都增加,但蛋白整體上在結合過程中結構穩(wěn)定。500 MPa及以上的高壓下蛋白的α-螺旋和β-折疊會受到一定程度的破壞、蛋白疏水表面積減少,蛋白的表面結構也會在高壓下發(fā)生較大的改變,特別是小分子與蛋白的結合口袋處,這些差異造成小分子與蛋白的結合自由能隨壓力的增加而降低,其中主要是范德華力和氫鍵的減少所致,而靜電作用力雖然高壓下有所增強,但影響有限。本研究為更好地了解蛋白與多酚在高壓下的結合,進而促進高壓在植物蛋白飲料殺菌的應用奠定了一定的理論基礎。

猜你喜歡
殘基氫鍵位點
教材和高考中的氫鍵
基于各向異性網絡模型研究δ阿片受體的動力學與關鍵殘基*
鎳基單晶高溫合金多組元置換的第一性原理研究
上海金屬(2021年6期)2021-12-02 10:47:20
CLOCK基因rs4580704多態(tài)性位點與2型糖尿病和睡眠質量的相關性
“殘基片段和排列組合法”在書寫限制條件的同分異構體中的應用
二項式通項公式在遺傳學計算中的運用*
生物學通報(2019年3期)2019-02-17 18:03:58
蛋白質二級結構序列與殘基種類間關聯的分析
基于支持向量機的蛋白質相互作用界面熱點殘基預測
二水合丙氨酸復合體內的質子遷移和氫鍵遷移
銥(Ⅲ)卟啉β-羥乙與基醛的碳氫鍵活化
大理市| 十堰市| 永定县| 大连市| 蓬溪县| 永昌县| 乾安县| 琼中| 小金县| 上饶市| 林口县| 遂溪县| 鸡泽县| 平和县| 乌鲁木齐县| 陕西省| 布尔津县| 长乐市| 乾安县| 长阳| 新和县| 天长市| 海兴县| 根河市| 固安县| 忻州市| 玛纳斯县| 利津县| 大宁县| 温宿县| 洪泽县| 砚山县| 五河县| 衡阳县| 阿合奇县| 盐城市| 禹城市| 平罗县| 乡宁县| 察哈| 淅川县|