燕山矿区苜蓿恢复过程中土壤养分与微生物的演变特征

时间:2023-06-21 18:00:02 来源:网友投稿

马建军,姚虹,刘辉,田美荣,3*

1.廊坊师范学院生命科学学院

2.廊坊泽通林业工程设计有限公司

3.环境基准与风险评估国家重点实验室, 中国环境科学研究院

矿区是在人为因素严重干扰下形成的一种极度退化的生态系统,矿产资源开采时会因原有的表土及植被剥离而导致水土流失加剧、生态系统功能降低。植被恢复是控制水土流失和恢复生态系统功能的唯一策略,是矿区土壤生态系统重新发育和演替的基础。自然恢复通常需要较长时间,生态威胁将可能长期存在。因此,在矿区建立一个高效、可持续的人工植被恢复系统已逐渐引起人们的重视和广泛关注[1-2]。

植被恢复对土壤生态系统具有深远影响,而土壤对植被恢复的响应存在滞后效应[3]。因此,有必要评估植被恢复对土壤养分及微生物的长期影响[4]。苜蓿(Medicago sativa)是京津冀地区分布广泛的人工植被类型,其抗碱能力强、生态适应性强,具有增加生态系统氮输入的潜力,同时可促进植被和土壤中碳、氮的积累[5-6]。苜蓿的土壤改良效果尤其明显,在加速植被恢复[7]、防止水土流失[8]、修复退化土壤等方面效果显著[9]。然而,种植年限直接影响着苜蓿的生长情况,在北方干旱、半干旱地区,苜蓿生长盛期持续时间一般为5~6年,6~7年产量迅速下降,10年后逐渐退化演变为自然植物群落[10],苜蓿的这种阶段性生长对土壤水分、理化性质、养分变化和土壤酶等方面造成影响。

目前,苜蓿对土壤的影响研究多集中在土壤理化性质[11-12]、修复重金属污染[13]、土壤微生物[14-15]等方面,而对于在极端恶劣矿区环境中,苜蓿地土壤养分与微生物演变规律及其土壤改良的时效性和持续性仍缺乏足够的研究[16],但这方面的研究对于持续恢复矿区脆弱的生态系统有着积极的指导作用。

燕山山脉是中国北部著名山脉之一,是首都北京的绿色屏障和水源地,是环首都生态建设的重点地区,也是我国重要的矿产基地。河北省三河市东部矿区是廊坊市唯一的山区,系燕山余脉南麓,盛产白云岩。自20世纪70年代至2015年作为京、津建筑原料的供应基地,长期石材的开采破坏了山体的外貌和植被,导致土地资源和环境资源丧失,使土地失去原有的利用价值。本矿区在开采的同时进行了植被恢复工程,其中以种植豆科牧草苜蓿为主,在矿区形成了不同种植年限的苜蓿地恢复序列。

笔者以该矿区采石场平台不同种植年限的苜蓿地为研究对象,通过分析苜蓿地土壤养分和土壤微生物的动态变化,探索苜蓿地恢复过程中土壤养分和微生物的演变规律,旨在阐明以引种苜蓿作为土壤改良的措施在燕山矿区生态恢复中是否是高效和可持续的,即与未恢复对照样地相比,苜蓿对土壤的改良效果是否是高效的;
与原生境对照样地相比,苜蓿对土壤的改良效果是否是可持续的,是否有明显的时效性。

1.1 研究区概况

研究区为三河市东部白云石露天矿区(117°05′13″E~117°11′10″E,40°00′23″N~40°02′34″N)(图1),行政区划属三河市黄土庄镇和段甲岭镇。三河市东部矿区地处燕山山脉东段南缘与华北平原接壤部位的低山丘陵区。三河市属典型暖温带大陆性气候,年均气温11.1 ℃,年均降水量617.4 mm,年蒸发量1 681.9 mm。矿区土壤类型以石灰性盐碱褐土为主,土质疏松,抗蚀抗冲性差,水土流失严重,有机质浓度低。矿区采石场原生表土已经全无,现多为缺乏养分与水分的新生石砾土和沙质土。

图1 研究区及样地设置Fig.1 Schematic diagram of study area and sample plots

1.2 样地设置与采样

2020年8 月,在坡位、坡向、海拔等立地条件基本一致的采石场平台上选择恢复年限分别为3年(Me3)、6 年(Me6)、10 年(Me10)和 15 年(Me15)的4个苜蓿恢复样地,苜蓿地种植后仅在前2年进行浇灌和2次刈割,以后未进行人工干预,处于自然演替状态。在矿区内选择已经开采但未进行恢复的区域(CK1)和有代表性且未经开采的自然区域(CK2)设置2个对照样地(图1)。样地基本情况见表1。

表1 样地基本情况Table 1 Basic properties of the sample plots

4个苜蓿恢复样地和对照样地CK1的表层覆盖客土,来源于研究区附近建筑工程地基开挖过程中所产生的开槽土和工程弃土。恢复期间所有样地均未施肥。

在每个样地中,设置3个10 m×10 m取样小区,小区之间及其与样地边界之间的距离均大于10 m。在每个采样小区中,按照S型多点采样法,沿S形路线均匀设置12个取样点(取样点距离植株基部3~5 cm左右),去除土壤表面凋落物层后,用土钻采集0~20 cm土壤样品,并进行混合,用无菌塑封袋密封后,用冰盒冷藏并迅速带回实验室。采集的土壤样品分为2份:一份过2 mm筛后室内自然风干,用于测定土壤养分;
一份于-20 ℃冷冻保存,用于土壤微生物群落的测定。

1.3 土壤养分的测定

速效氮(AN)浓度采用碱解扩散法测定;
速效磷(AP)浓度采用盐酸和硫酸溶液浸提法测定;
速效钾(AK)浓度用乙酸铵浸提-火焰光度法测定[17];
土壤有机质(OM)浓度采用重铬酸钾容量法-外加热法测定。每个样品做7个重复。

1.4 土壤微生物群落分析

采用CTAB法对不同样品的总DNA进行提取,经1%琼脂糖凝胶电泳测定抽提的基因组DNA。扩增细菌16S rRNA基因的V3+V4_b高变区段,引物为 338F(5"-ACTCCTACGGGAGCAGCA-3")和 806R(5"GGAC TACHVGGGTWTCTAAT-3")。真菌多样性对18S rRNA的ITS1_f区段进行测序,引物为ITS1F(5"-CTTGGTC ATTTAGAGGAAGTA A-3")和2043R(5"-GCTGCGTTCTTCATCGATGC-3")。扩增条件为 98 ℃ 预变形 1 min,98 ℃ 变形 10 s,50 ℃ 退火 30 s,72 ℃ 延伸 60 s,30 个循环,72 ℃ 延伸形5 min[18]。土壤样品完成DNA提取后,送北京百迈客生物科技有限公司进行Illumina高通量测序。每个样品做7个重复。

1.5 数据处理与分析

所有数据采用Excel及R-3.6.3中的Vegan分析包进行统计分析。根据方差分析和Tukey检验分析样品间的显著性差异。分析土壤细菌、真菌群落的Beta多样性时,首先对操作分类单元(operational taxonomic units, OTUs)表进行Hellinger标准化,采用R-3.6.3的Vegan包计算Bray-Curtis差异性距离(Bray-Curtis dissimilarity,dBCD),并进行相似性分析(analysis of similarities, ANOSIM)及置换多元方差分析(permutation-based multivariate ANOVA,PerMANOVA),分析各分组之间是否存在显著差异,并进一步通过主坐标分析(principal co-ordinates analysis, PCoA)进行可视化组间差异分析[19-20]。其中,dBCD计算公式如下:

式中:SA,i和SB,i为第i个OTU分别在A群落和B群落中的计数。dBCD为0~1时,数值越大,则2个样本之间的相似性就越低;
如果2个样本相同,则dBCD=0。

2.1 土壤养分演变特征

种植苜蓿可以有效地改善土壤养分[21-22]。由图2可见,苜蓿对于土壤AN浓度、AP浓度、AK浓度和OM浓度的改良效果有一定差异。

图2 不同样地土壤AN、AP、AK及OM浓度Fig.2 Concentrations of soil AN, AP, AK and OM in different sample plots

由图2(a)可知,随着恢复年限的增加,土壤AN浓度呈先增加后降低的变化规律。种植苜蓿的初期,样地Me3与未恢复对照样地CK1间土壤AN浓度无显著差异;
随着恢复年限的增加,样地Me6中土壤AN浓度达到峰值,其浓度显著高于除原生境对照样地CK2外的所有样地。这是由于苜蓿根部形成了大量根瘤菌,其具有固氮作用;
另外一部分根系和地上茎叶腐烂进入土壤,增加了土壤养分浓度。随着恢复年限的进一步增加,在恢复的第10年和第15年,土壤AN浓度持续降低。研究发现,本矿区苜蓿的生长盛期为6年,之后苜蓿进入衰退期,且随着种植年限的继续增加,根系活力下降,生物量减少,对氮素利用强度减弱,而导致的氮素归还减少,苜蓿地土壤肥力并不会继续保持增长[23-24]。

由图2(b)可知,随着恢复年限的增加,土壤AP浓度呈先增加后降低再增加的变化规律。尽管在样地Me6中土壤AP浓度达到峰值,但是所有样地间土壤AP浓度均无显著性差异,表明种植苜蓿对土壤AP浓度的影响并不明显。研究区各样地土壤AP浓度整体偏低(均在10 mg/kg以下),苜蓿对土壤AP的消耗甚微,甚至可以积累AP。随着种植苜蓿年限的延长,土壤AP浓度呈下降趋势[25]。

由图2(c)可知,随着恢复年限的增加,土壤AK浓度呈现与土壤AN浓度相同的变化规律。苜蓿对于土壤AK浓度的影响较为明显,样地Me6中土壤AK浓度达到峰值,其浓度显著高于所有样地。苜蓿在生长过程中对钾元素的需求量也很高,在本研究中样地Me3和Me6土壤AK浓度显著增加,这是由于在研究区盐碱性土壤中钾元素浓度较低,且主要以矿物态钾和非交换态钾形式存在,交换态钾和水溶态钾浓度较低[26]。不同形态间的钾存在着动态的平衡,苜蓿生长过程中需要大量的AK,而各苜蓿地中AK浓度依然高于对照样地CK1和CK2,说明长势旺盛的苜蓿(Me3和Me6)可以大大促进土壤AK的转化[27]。随着苜蓿生长的衰退(Me10和Me15),苜蓿对钾的转化能力逐渐减弱,同时大量本土植物在样地中出现,植物在生长过程中的吸收及淋溶作用会造成土壤AK的减少,因此随着种植年限的延长,土壤AK浓度降低[14]。由图2(d)可知,随着恢复年限的增加,土壤OM浓度总体呈增加趋势,且均高于对照样地CK1,但均显著低于对照样地CK2,这与样地中枯落物的淋溶腐解、腐殖质含量增加有关[28-29]。

有研究显示,苜蓿在不同地区有着7年[10]、10年[30]、11年[31]等不同的生长盛期,本研究区苜蓿的生长盛期为6年,这与采石场土壤现状(开槽土和工程弃土)、利用方式(保土固沙、2次刈割)和管理措施(未施肥、长期自然演替状态)等有关。因此,应加强苜蓿地的科学管理和合理利用以延长苜蓿的生长时间[32]。

2.2 土壤细菌群落结构演变特征

不同样地优势细菌的相对丰度如图3所示。本研究区细菌类群中相对丰度水平前5的优势门依次为变形菌门(Proteobacteria)、酸杆菌门(Acidobacteria)、放线菌门(Actinobacteria)、厚壁菌门(Firmicutes)、拟杆菌门(Bacteroidetes)〔图3(a1)~图3(a5)〕;
优势科依次为鞘氨醇单胞菌科(Sphingomonadaceae)、梭菌科(Pyrinomonadaceae)、伯克氏菌科(Burkholderiaceae)、根瘤菌科(Rhizobiaceae)、芽孢杆菌科(Bacillaceae)〔图3(b1)~图3(b5)〕。

由图3(a1)可知,变形菌门的相对丰度呈先增加后降低的变化规律,样地Me10中变形菌门的相对丰度显著高于所有样地;
其他样地间,该门的相对丰度无显著差异;
该门中的鞘氨醇单胞菌科〔图3(b1)〕和伯克氏菌科〔图3(b3)〕呈现与该门一致的变化规律,而根瘤菌科〔图3(b4)〕则呈现较复杂的变化规律:峰值出现在样地Me10,而Me6、Me15及CK2样地中根瘤菌科的相对丰度均较低且无显著差异。变形菌门细菌在不同环境中广泛分布,适应能力强[33],是碱性土壤中的主要优势群落,广泛存在于以盐碱土壤为主的矿区土壤中[34]。与CK1相比,富营养型细菌类群变形菌门的相对丰度在所有苜蓿样地中有所增加,这与黄土高原草地植被自然演替过程中的细菌群落组成的变化规律表现出相同的趋势[35]。其中样地Me3、Me6和Me10中变形菌门的相对丰度均显著高于CK1和CK2。变形菌门在土壤中所占比例越大,在一定程度上代表了土壤越肥沃,在一些黑土地及半湿润地区的土壤中变形菌门往往是优势菌群[36]。本研究也证实,在种植苜蓿初期,可以改善土壤变形菌门分布状况。

酸杆菌门〔图3(a2)〕的相对丰度呈先增加后降低再增加的变化规律,第6年其相对丰度显著高于其他苜蓿样地及对照样地CK1,而显著低于样地CK2;
其中,隶属于该门的梭菌科〔图3(b2)〕的相对丰度随着苜蓿种植年限的增加而增加。酸杆菌门在自然环境中亦十分常见,可以降解植物纤维素等大分子聚合物[37]。本研究发现酸杆菌门主要存在于原生境样地CK2中,可能与样地中有大量的荆条、沙枣等灌木植物有关。

放线菌门〔图3(a3)〕和厚壁菌门〔图3(a4)〕呈现相同的变化规律,即二者的相对丰度先降低后增加。样地Me3、Me6中,放线菌门的相对丰度高于样地CK1(但无显著差异),表明种植初期苜蓿虽能够改善放线菌状况,但效果不明显。样地CK1中厚壁菌门的相对丰度最高,而该门的相对丰度在样地Me6、Me10及CK2间无显著性差异;
隶属于厚壁菌门的芽孢杆菌科〔图3(b5)〕的相对丰度的变化规律与厚壁菌门的变化规律一致。在干旱和寡营养的土壤中放线菌门和厚壁菌门的相对丰度较高[38]。厚壁菌门细菌细胞壁厚,结构简单,可以产生芽孢,它可以抵抗脱水和极端环境,因而能较好地适应矿区恶劣环境[39]。样地CK1几乎无植被覆盖,土壤水分含量低,故其厚壁菌门的相对丰度最高,在样地Me3和Me6中,植被覆盖显著增加,土壤水分含量有所增加,故厚壁菌门的相对丰度呈下降趋势。样地Me10(苜蓿呈现衰退迹象)和Me15(苜蓿基本消失)的土壤含水量减少,此时相对丰度出现增加的趋势,原生境样地CK2中厚壁菌门的相对丰度较低,这些差异说明厚壁菌门的相对丰度与植被恢复的不同阶段土壤含水量有关[40]。

图3 不同样地优势细菌的相对丰度Fig.3 Relative abundance of dominant bacteria in different sample plots

拟杆菌门〔图3(a5)〕的相对丰度则随着恢复时间的增加而增加,样地CK1、Me15的拟杆菌门的相对丰度显著最高,而其他样地间,该门的相对丰度则无显著差异。拟杆菌门在水体、土壤及沉积物中也是广泛分布的一个类群,在研究区该门也有较高的平均相对丰度。样地Me3、Me6、Me10和Me15中拟杆菌门的相对丰度逐渐增加,但是均低于样地CK1,这与变形菌门、酸杆菌门、放线菌门等优势菌群在上述样地中竞争性增强有关[41-42]。

2.3 土壤真菌群落结构演变特征

不同样地优势真菌相对丰度如图4所示。本研究区真菌类群中相对丰度水平前5的优势门依次为子囊菌门(Ascomycota)、担子菌门(Basidiomycota)、毛 霉 门 (Mucoromycota) 、 隐 真 菌 门(Rozellomycota) 、 被 孢 霉 门 (Mortierellomycota)〔图4(a1)~图4(a5)〕;
优势科及其平均相对丰度依次为曲霉科(Aspergillaceae)、毛壳菌科(Chaetomiaceae)、被孢霉科(Mortierellaceae)、粉褶菌科(Entolomataceae)、枝孢霉科(Cladosporiaceae)〔图4(b1)~图4(b5)〕。

在真菌群落中,子囊菌门和担子菌门是自然界中分布最广、丰度最高的2个真菌类群,这与其代谢特性及在多种生境中较强的生存能力有关[43]。本研究区子囊菌门〔图4(a1)〕和担子菌门〔图4(a2)〕的相对丰度均较高,这与Yang等[44-45]在黄土高原半干旱区退耕草地中研究的结果相似。子囊菌门和担子菌门大多数为腐生菌,本研究中各样地土壤pH为7.5~8.5,碱性土壤最适合腐生真菌的生长,这可能是子囊菌门和担子菌门为优势菌门的一个原因[46]。

由图4(a1)可知,随着恢复年限的增加,子囊菌门的相对丰度呈先降低后增加趋势。子囊菌门是土壤腐生真菌,易受到植物种类和植物残茬的强烈影响,其功能是分解木质化植被碎屑[47]。本研究中各苜蓿样地中植被覆盖度及植被残茬均显著增加,使子囊菌门真菌能够更好地利用可降解的植被残茬,促进菌群的快速增长与繁殖[48],故各苜蓿地子囊菌门的相对丰度均显著高于CK1。本研究中,隶属于子囊菌门的优势科为曲霉科、毛壳菌科和枝孢霉科。曲霉科〔图4(b1)〕的相对丰度随着恢复年限增加而增加,而毛壳菌科〔图4(b2)〕与枝孢霉科〔图4(b5)〕呈现完全相反的变化规律,即随着恢复年限增加,毛壳菌科的相对丰度是先降后增,而枝孢霉科的相对丰度是先增后降。

图4 不同样地优势真菌相对丰度Fig.4 Relative abundance of dominant fungi in different sample plots

担子菌门〔图4(a2)〕呈现与子囊菌门相反的变化规律,这与子囊菌门真菌优势度增加有关[49-50]。自然界中超过98%的陆生真菌属于子囊菌门和担子菌门,而且前者的物种多样性明显多于后者,子囊菌门的物种数量是担子菌门的2倍多[51]。本研究中,隶属于担子菌门的粉褶菌科〔图4(b4)〕的相对丰度随着恢复年限增加而增加。

毛霉门〔图4(a3)〕和隐真菌门〔图4(a4)〕的变化规律较为一致,即二者的相对丰度随恢复年限的增加而呈现降低的趋势。毛霉门真菌能够吸收糖和简单的多糖以及N、P、K等营养物质,其孢子萌发快,菌丝生长速度快。因此,在植被恢复的初期最先形成优势菌,然而毛霉门真菌对自身代谢副产物的积累较为敏感,尤其是环境中CO2的积累,使其停止生长,产生休眠结构,进入休眠状态。毛霉门真菌的相对丰度在苜蓿种植初期(3~6年)之后逐渐较低,这一演替阶段随着毛霉门真菌自身可利用的营养物质的消耗和土壤中CO2的积累而消失,样地Me15中毛霉门真菌几乎消失。隐真菌门是在水体及极端环境中广泛分布的真菌,在土壤中相对较少[52-53]。本研究显示,各样地的隐真菌门无论从相对丰度(均较低),还是从变化规律(逐年递减)呈现出非常一致的变化规律。

随着苜蓿种植年限的增加,被孢霉门〔图4(a5)〕的相对丰度在样地Me6中达到峰值(与Me10间无差异);
在样地Me3、Me15和CK2中,其相对丰度均极低;
在样地CK1、Me6、Me10中,其相对丰度均显著高于样地CK2。被孢霉门是土壤养分丰富的标志类群,与土壤养分有密切联系[54]。本研究显示,尽管各样地中被孢霉门的相对丰度较低,但是呈现了一定的规律性,即随苜蓿种植年限的增加呈先升后降趋势,最高值出现在第6年,这与土壤AN和AK变化规律相一致。本研究中,随着恢复年限的增加,隶属于被孢霉门的被孢霉科〔图4(b3)〕的相对丰度呈现与该门较为一致的变化规律。

2.4 细菌/真菌群落Beta多样性演变特征

为了进一步阐明不同恢复年限苜蓿地间微生物群落组成上的差异性,本研究基于Bray-Cutis距离矩阵采用PCoA分析了6个样地土壤细菌及真菌群落物种组成的相似度,结果见图5、表2。

表2 基于Bray-Curtis距离算法的样地间土壤细菌/真菌的距离(dBCD)矩阵Table 2 Distance matrix of soil bacteria / fungi between sample plots based on Bray-Curtis distance(dBCD) algorithm

由图5可见,土壤细菌和真菌的7个重复样本聚在一起,这说明样本的重复性较好,组内变异相对较小。在细菌PCoA分析中PCoA第1轴和PCoA第2轴分别解释了变异信息的32.59%、16.36%,二者累计贡献率达48.95%。结果表明,在第1轴方向上,样地Me3与样地CK1的细菌群落组成较为相似,样地Me6和Me10相似,样地Me15与样地CK2相似性较低。在真菌PCoA分析中PCoA第1轴和PCoA第2轴分别解释变量方差的24.97%、14.69%,二者累计贡献率达39.66%。结果表明在第1轴方向上,随着苜蓿种植年限的增加,土壤真菌群落组成发生明显的阶段性变化,样地Me3和Me6与样地CK1相似,而样地Me15与CK2相似。

图5 土壤微生物PCoA分析图Fig.5 PCoA analysis of soil microorganisms

土壤细菌和真菌随着恢复时间的增加呈现逐渐本土化的演变趋势,但二者的演变进度存在差异。由Bray-Curtis距离矩阵(表2)可知,与未恢复样地CK1相比,在种植苜蓿的第10年,样地Me10与样地CK1的细菌群落的Bray-Curtis距离就达到了0.524 4,而此时两样地的真菌群落的Bray-Curtis距离为0.245 2;
在种植苜蓿的第15年,样地Me15与样地CK1的真菌群落的Bray-Curtis距离为0.572 2。各样地与样地CK2的真菌群落的Bray-Curtis距离也明显大于与细菌群落的Bray-Curtis距离。显然,苜蓿地土壤细菌群落演变较快,而真菌群落的演变相对缓慢。

Dangi等[55]研究了美国怀俄明州东北部露天煤矿复垦地经过14年的土地复垦,土壤真菌群落结构基本达到正常水平。本研究与Dangi等[55]研究结果存在差异的原因在于苜蓿对土壤细菌群落的碳源供应更充足,恢复能力很强[56],所以本研究中细菌的恢复进程更快,而真菌在竞争中处于劣势,恢复较为缓慢[57]。另有研究显示,随着生态系统的成熟,土壤微生物生物量的主导地位将由细菌向真菌转变[58]。

(1)与未恢复对照样地相比,苜蓿可以高效地提高土壤AN、OM、AK等养分的浓度以及变形菌门、酸杆菌门、子囊菌门、担子菌门、鞘氨醇单胞菌科、伯克氏菌科、根瘤菌科、芽孢杆菌科、毛壳菌科、粉褶菌科等优势微生物类群的相对丰度。

(2)与原生境对照样地相比,苜蓿的土壤改良效果有明显的时效性。苜蓿对土壤养分的恢复以第6年为最佳,之后随着恢复年限的延长,土壤AN浓度及AK浓度呈下降趋势;
在第6年或第10年,土壤优势微生物类群的相对丰度显著升高;
苜蓿对细菌的恢复进度较快,而对真菌的恢复较慢。苜蓿地恢复中应充分考虑气象、土壤等环境因素,根据不同区域特点对苜蓿地进行科学的管理与利用,以延长其土壤改良的时间。

猜你喜欢子囊菌门苜蓿苜蓿的种植及田间管理技术现代畜牧科技(2021年9期)2021-10-13破壁方式对冠突散囊菌有性繁殖体提取物抑菌活性的影响食品与生物技术学报(2021年9期)2021-09-28苜蓿的种植技术现代畜牧科技(2021年4期)2021-07-21冠突散囊菌繁殖体提取物抑菌活性研究食品与发酵工业(2020年14期)2020-07-29野生树鼩与人工饲养树鼩消化道不同部位微生物组成的比较研究中国比较医学杂志(2020年4期)2020-05-26粗糙脉孢菌7种子囊型归类教学探究遗传(2019年11期)2019-11-28饥饿与重摄食对河蟹肠道菌群结构的影响水生生物学报(2019年4期)2019-07-20昆虫体内微生物多样性的影响因素研究进展生物安全学报(2019年3期)2019-02-15妊娠期糖尿病肠道菌群变化及临床价值分析川北医学院学报(2019年6期)2019-02-10苜蓿:天马的食粮中国三峡(2017年4期)2017-06-06

推荐访问:燕山 苜蓿 矿区