低温条件下松材线虫在中国的风险分布区预测
2017-08-29赵捷，韩骁，石娟
赵 捷， 韩 骁， 石 娟
北京林业大学林学院，北京 100083
低温条件下松材线虫在中国的风险分布区预测
赵 捷， 韩 骁， 石 娟*
北京林业大学林学院，北京 100083
【目的】松材线虫是我国一种重要的检疫性外来入侵物种，给我国林业造成巨大的经济损失。近年来在松材线虫疫区的调查发现，一些年平均温度8 ℃左右的高海拔地区也开始有松材线虫病的发生，说明其分布线可能开始向北和高海拔地区移动。因此,在松材线虫低温适应性驯化环境下，对松材线虫在我国的最新风险分布区进行了预测。【方法】将松材线虫分为长期(30 d)暴露和短期(24 h)暴露2组进行低温(3、1、0、 -1、 -3、 -5、-6、 -7、 -9、 -12 ℃)暴露实验，得出致死温度LT50、LT75和LT99。利用CLIMEX 1.0适生性分析软件与ArcGIS 10.2地理信息系统软件以LT50、LT75和LT99为界限划分风险分布区。利用CLIMEX软件导入新的温度数据对松材线虫进行适生性分析，将结果导入ArcGIS中，进行IDW插值，以EI值划分松材线虫风险分布区。对有无低温胁迫及不同低温胁迫时间下的中度风险分布北线的变化趋势进行分析。【结果】松材线虫长期冷胁迫(30 d)条件下的中度风险分布区界线比短期冷胁迫(24 h)条件下和无胁迫条件下(以EI值划分)的中度风险分布区界线明显靠北。【结论】随着低温胁迫时间的延长，松材线虫分布区有向北扩散的趋势，即产生一定的低温适应性进化。
松材线虫； ArcGIS； CLIMEX； 风险分析； 低温适应性
松材线虫Bursaphelenchusxylophilus(Steiner & Buhrer) Nickle是我国一种重要的检疫性外来入侵物种，松材线虫病发病致死速度快,传播蔓延迅速,难以预先诊断，且松树一旦感染就无法治愈,防治难度大，因而被称为“松树的癌症”(朱丽华等，2011)。该病原产于北美(陈守常，2010)，在美国较为普遍，但对本地松树种类危害并不严重(潘沧桑，2011)，只危害几种由亚洲和欧洲引进的树种。在加拿大和韩国，只有几种松树受到危害(徐亮，2013)。松材线虫病在日本危害更严重，自发现一个多世纪以来，已经扩散到除北海道地区之外的45个县级地区，造成了巨大的经济损失(Kazuyoshi, 2013)。1982年,我国南京地区首次发现该病，现已扩散到我国15个省(自治区/直辖市)的210个县级行政区(国家林业局，2016)，导致松林大面积死亡，对我国森林生态系统造成了巨大危害。
日本的研究发现，松材线虫病在年平均气温低于10 ℃的地区不发生；在年平均气温10～12 ℃的地区能够生存并零星发生；在年平均气温12～14 ℃的地区可以流行；在年平均气温高于14 ℃的地区可暴发流行。故松材线虫对低温的适应能力是决定其在自然条件下存活和分布的重要因素之一。研究者们利用不同模型对松材线虫在中国的风险分布区做过预测(程功等，2015； 韩阳阳等，2015； 潘红伟，2009； 张海军，2007)。松材线虫病已入侵我国30余年，其分布区逐渐扩大，但因其主要寄主是马尾松PinusmassonianaLamb.和黑松PinusthunbergiiParl.，以及气候限制(如松材线虫病在年平均温度低于10 ℃的地区不发生)等方面的原因，主要发病区一直受到局限。2011年，在较高纬度的陕西高海拔(1700 m)区域也发生松材线虫病危害，意味着松材线虫分布区已突破以前预测的范围。因此，研究松材线虫在我国向北、向高纬度地区扩散蔓延的趋势，预测最新风险分布区十分必要，能够为新形势下更加客观有效地防控松材线虫提供指导，对北方油松林和高海拔地区松材线虫的防控工作也具有重要的生态学意义。
1 材料与方法
1.1 松材线虫低温暴露实验
低温存活率为衡量线虫耐寒性强弱(Jagdale & Grewal,2003; Leduc & Wharton,2008； Wharton & Block,1993)的普遍做法。松材线虫的生活史可分为2个阶段(Kikuchietal.,2011)，包括繁殖周期和扩散周期(张伟，2014)，且繁殖型的2龄幼虫蜕皮成的分散型3龄幼虫(dispersal third-stage juvenile)是松材线虫可以长期忍受恶劣环境的休眠阶段(黄瑞芬，2015)，所以，测定松材线虫分散型3龄幼虫的低温耐受能力更适合评价其耐寒性。
1.1.1 供试虫源 松材线虫样品采自湖北省宜昌市夷陵区疫区，用“贝尔曼漏斗法”从马尾松越冬病木中分离收集线虫，置于显微镜下观察，经形态鉴定后确定为分散型3龄松材线虫。
1.1.2 实验方法 (1)短期暴露：设置5个温度梯度0、-3、-6、-9、-12 ℃，每组从显微镜下分别挑取分散型3龄松材线虫幼虫100条，转移至200 μL的离心管中，离心管中注入40%甘油溶液，将离心管置于5个温度的恒温培养箱(上海智城ZWY-211C；一恒LRH-250CA)中静置24 h。将低温处理后的松材线虫在室温下放置24 h后，经物理刺激(用移液枪吹打数次)后，在显微镜下观察，能进行较强地卷曲和伸展运动，头部有轻微卷曲或尾部微弱摆动的线虫均判为存活线虫，不动的线虫判为死亡线虫。分别统计各处理线虫的存活数与死亡数，计算存活率，统计3次，取平均值。
(2)长期暴露：设置6个温度梯度3、1、-1、-3、-5、-7 ℃，每组从显微镜下分别挑取分散型3龄松材线虫幼虫100条，转移至200 μL的离心管中，将离心管置于5个温度的恒温培养箱(同上)中静置30 d。低温处理后物理刺激及判别统计方法同上。
1.1.3 数据处理 使用SPSS将松材线虫死亡例数作为响应频率，观察例数作为观测值汇总，温度作为协变量，进行数据分析。分别计算得出松材线虫短期(24 h)暴露和长期(30 d)暴露的半致死温度LT50、75%致死温度LT75、99%致死温度LT99。
1.2 分析软件及数据处理
采用ArcGIS v10.2以及CLIMEX for Windows 1.1软件的单物种地区比较[compare locations (1 species)]功能进行数据处理分析。
根据ArcGIS所需的气候数据格式，通过中国气象数据网(http:∥data.cma.cn/)提供的气象站点数据，收集了2010—2015年全国97个站点的1月平均温度和年最低日平均温。下载中国矢量地图(http:∥down.51cto.com/data/967123)，气象站点分布如图1所示。利用反距离权重插值法(inverse distance weighted, IDW)进行插值运算，分别以1月平均气温和日平均最低温为Z值，修改1月平均气温IDW插值结果的属性，根据松材线虫长期(30 d)低温暴露实验的结果，以LT50、LT75、LT99为划分值，分别生成3条等温线，并以此作为松材线虫风险分布区的划分界限，将松材线虫的风险分布区划分为高风险分布区(
CLIMEX拥有全球2031个地区的1961—1990年的气候资料数据，包括85个中国气象站点(宋红敏等，2004)。结合中国气象科学数据共享服务网中国地面气候标准值数据集(1981—2010年)中包含的885个可用气象站点数据，加上CLIMEX中包含的台湾省的2个站点(台北和花莲)和香港的1个站点(香港)，一共888个站点数据组成了中国气候数据，用于风险分布区预测(图2)。参考湿地松粉蚧适生性分析(陈燕婷和尤民生，2014)及橘小实蝇适生性研究(王涛等，2015)，结合低温暴露实验结果以及何善勇等(2012)的研究，根据松材线虫在全球分布的模拟结果与实际分布达到最大吻合后，确定研究松材线虫的最佳参数，如表1所示。
图1 气象站点分布图Fig.1 Locations of the meteorological stations
CLIMEX参数描述CLIMEXparameters设定值ValueCLIMEX参数描述CLIMEXparameters设定值Value发育起点温度Lowerthresholdtemperature(DV0)10冷胁迫开始积累的阈值Coldstresstemperaturethreshold(TTCS)-12．67适宜温度下限Loweroptimumtemperature(DV1)15冷胁迫积累速率Coldstressaccumulationrate(THCS)0．15适宜温度上限Upperoptimumtemperature(DV2)28热胁迫开始积累的阈值Heatstresstemperaturethreshold(TTHS)35发育最高温度Upperthresholdtemperature(DV3)33热胁迫积累速率Heatstressaccumulationrate(THHS)0．002有效积温Degree⁃daysnecessarytocompleteonegeneration(PDD)150干胁迫开始积累的阈值Drystresssoilmoisturethreshold(SMDS)0．02发育需要的最低土壤湿度Lowerthresholdofsoilmoisture(SM0)0．2干胁迫积累速率Drystressaccumulationrate(HDS)0．05适合发育的土壤湿度下限Lowerlimitofoptimumsoilmoisture(SM1)1湿胁迫开始积累的阈值Wetstresssoilmoisturethreshold(SMWS)5适合发育的土壤湿度上限Upperlimitofoptimumsoilmoisture(SM2)2湿胁迫积累速率Wetstressaccumulationrate(HWS)0．0015发育需要的最高土壤湿度Upperthresholdofsoilmoisture(SM3)3
松材线虫适生区分析结果如图2，并将得出的各站点的坐标和生态气候参数(ecoclimatic index, EI)值导入ArcMap 10.2中，利用IDW修改对于EI值的IDW插值结果属性，依照Chejaraetal.(2010)的界定，以0≤EI<0.5的地区为无风险区，0.5≤EI<10的地区为低风险区，10≤EI<20的地区为中度风险区，EI≥20的地区则为高风险区。
图2 松材线虫适生性分析结果Fig.2 Potential distribution of B. xylophilus considering adaptation to cold conditions
2 结果与分析
2.1 松材线虫低温暴露实验结果
经过SPSS分析，分别计算得出松材线虫短期(24 h)暴露和长期(30 d)暴露的半致死温度LT50、75%致死温度LT75、99%致死温度LT99(表2)。
表2 分散型松材线虫在不同胁迫时间下的致死温度Table 2 Lethal temperatures of B. xylophilus when exposed to low temperatures
2.2 以等温线划分的松材线虫风险分布区
2.2.1 松材线虫长期(30 d)冷胁迫条件下的风险分布区 图3显示，在1月份平均温度气候条件下，松材线虫病的风险分布区分布广泛，除了黑龙江、吉林、辽宁和新疆无松材线虫风险分布区外，其余省市均有风险区域。其中，高风险分布区与当前松材线虫在我国的实际分布区基本符合，包括海南、广西、广东、福建、台湾、云南、贵州、湖南、江西、浙江、上海、重庆、湖北、安徽、江苏、河南、陕西南部、河北邯郸、甘肃陇南、山西南部边缘地区、四川除西北部以外的所有地区、山东南部及泰山地区,以及西藏日喀则、拉萨和山南地区；中度风险分布区包括北京南部、天津大部分地区、西藏那曲和昌都地区、四川阿坝和甘孜藏族自治区中部、甘肃东南部分地区、陕西中部地区、山东除南部区域外所有地区、山西南部地区、宁夏南部地区及河北南部；低风险分布区包括西藏(双湖县、尼玛县及萨疆县)、青海省与西藏交界地区、四川西北部、甘肃兰州和白银地区、陕西中北部、宁夏除石嘴山外所有地区、内蒙古与甘肃和宁夏交界地区、山西中部、山东泰山地区、天津北部地区、河北保定和唐山地区。
2.2.2 松材线虫短期(24 h)冷胁迫条件下的风险分布区 图4显示，在年最低日平均温度气候条件下，松材线虫的风险分布区没有1月平均温度下的风险分布区广，除了黑龙江、吉林、内蒙古和新疆无风险区外，其余省市均有风险区域。其中，高风险分布区主要集中在我国南方地区，包括海南、广西、广东、福建、台湾、贵州、湖南、江西、浙江、上海、重庆、湖北、安徽、江苏、云南除西北地区外的所有地区、西藏日喀则地区、四川东部地区、陕西安康地区、河南黄河以南地区及山东南部沿海地区；中度风险分布区包括云南大理和保山地区、西藏日喀则周围地区、四川(成都、广元、雅安及凉山彝族自治州部分地区)、陕西西安和汉中地区、甘肃与陕西交界局部地区、山西南部边缘地区、山东除泰山和南部沿海外的南部地区；低风险分布区主要包括西藏(日喀则边缘地区)、青海西宁东部地区，四川中部部分地区，甘肃兰州和庆阳地区、陕西中部、山西东南地区、北京南部地区、天津南部地区、河北石家庄地区、山东泰山地区。
图3 松材线虫在长期(30 d)胁迫下的风险分布区Fig.3 Regions for potential invasion of B. xylophilus using the long-term (30 d) cold stress data
图4 松材线虫在短期(24 h)胁迫下的风险分布区Fig.4 Regions for potential invasion of B. xylophilus using the short-term (24 h) cold stress data
2.3 松材线虫无胁迫条件下的风险分布区
如图5所示，根据CLIMEX计算得出的松材线虫在全国各地的EI值，松材线虫风险分布区分布于除黑龙江、吉林、新疆以外的所有省份。其中，高风险分布区主要集中在我国南方地区，包括海南、广西、广东、福建、台湾、贵州、湖南、江西、浙江、上海、重庆、湖北、安徽、江苏、云南除普洱地区外的所有地区、西藏山南地区、四川东部和南部地区、陕西南部地区、河南南部地区及山东与江苏交界地区。中度风险分布区包括云南普洱地区，西藏左贡和灵芝地区、四川(成都、广元、雅安及凉山彝族自治州北部地区)、陕西西安和汉中地区、甘肃与陕西交界地区、山西南部边缘地区、山东除泰山和江苏边界外的南部地区、天津沿海地区、河南中部地区。低风险分布区包括西藏(拉萨、日喀则、那曲地区)、青海西宁东部地区，四川中部部分地区，甘肃兰州以南地区、陕西延安地区、宁夏南部地区、山西除北部外大部分地区、北京大部分地区、天津除沿海外地区、河南北部地区、河北除张家口以外地区、山东泰山及以北地区、辽宁大部分地区(包括大连)。
图5 以EI值划分的松材线虫风险分布区Fig.5 Regions for potential invasion of B. xylophilus data based on EI value
2.4 有无低温胁迫及不同低温胁迫时间下的中度风险分布北线变化趋势
将长期胁迫和短期胁迫下(以等温线划分)的松材线虫风险分布区(图3、图4)与无胁迫下(以EI值划分)的松材线虫风险分布区(图5)做比较。以中度风险分布区的界线为例做对比(图6)，松材线虫长期(30 d)冷胁迫条件下的中度风险分布区界线比短期(24 h)冷胁迫条件下和无胁迫条件下(以EI值划分)的中度风险分布区界线明显靠北。这表明从无胁迫到长期低温胁迫，以及从短期低温胁迫到长期低温胁迫，松材线虫的风险分布线有向北即向寒冷地区扩散的趋势。其中，短期(24 h)冷胁迫条件下和无胁迫条件下(以EI值划分)的松材线虫中度风险分布区界线差异很小。
图6 松材线虫中度风险分布区界线对比Fig.6 Comparison of the limits of distribution of B. xylophilus based on moderate risk
3 讨论
松材线虫长期(30 d)冷胁迫条件下的中度风险分布区界线比短期(24 h)冷胁迫条件下和无胁迫条件下(以EI值划分)的中度风险分布区界线明显靠北，表明随着低温胁迫时间的延长，松材线虫分布区有向北扩散的趋势，即产生一定的低温适应性进化。然而，短期(24 h)冷胁迫条件下和无胁迫条件下(以EI值划分)的松材线虫中度风险分布区的北界线差异很小，且部分区域短期胁迫条件下的分布线还向南移动。原因可能是24 h的低温胁迫时间不仅不足以使松材线虫产生一定程度的低温适应性驯化，而且有可能因为骤然的低温限制其向北扩散。
从图3和图4看，以高风险分布区为例，松材线虫在长期低温胁迫下的风险分布区要广于短期低温胁迫下的风险分布区，表明松材线虫分散型3龄幼虫受短时间内(24 h)的低温胁迫的限制要大。因此，对于松材线虫种群来说，在有些地区，长期且正常的自然冷胁迫条件并不是影响和限制松材线虫风险分布区的重要因素，但若出现冬季极端低温异常增加，将会对该种群的稳定繁衍和暴发及其风险分布区的扩散产生较大的制约作用。短期冷胁迫条件下，高风险分布区与最新松材线虫病疫区界线公告(国家林业局，2017)最为吻合。在无胁迫、短期低温和长期低温胁迫3种情况下，高风险分布区差异较小，中度风险分布区和低风险分布区的差异均较大，特别是在西藏、青海、甘肃、陕西、四川、河北、辽宁、山西这些风险分布边缘省份。其原因可能在于这些省份独特的地貌，这些地区大多位于高海拔、高山地区，因此，单独从温度因素判断并不能准确研究松材线虫在这些地区的风险分布。
除温度因素外，还要考虑媒介昆虫和寄主因素。就媒介昆虫而言，对松材线虫媒介昆虫——松褐天牛MonochamusalternatusHope在全国适生区的研究表明：松褐天牛在全国各地都有分布，且其适生区要广于松材线虫风险分布区(张海军，2007)。就寄主而言，松材线虫寄主种类多，目前已知自然条件下感病的寄主有47种，其中，松树38种，其他针叶树9种；人工接种条件下，另有17种松树、4种非松属针叶树能够感病。累计属于松属的寄主树种55种，非松属针叶树13种，共达68种。Shietal.(2013)在我国陕西油松林分离出松材线虫，表明在北方广布的油松树种已能在自然条件下感病。
综合本研究结果，随着松材线虫的低温适应性进化的产生，除了国家林业局公告的松材线虫病疫区外，松材线虫风险程度较高的分布区主要包括四川(成都、雅安地区)、西藏(日喀则地区)、青海西宁东部地区、甘肃陇南地区、陕西西安地区、宁夏南部地区、山西晋城地区、河南大部分地区、山东大部分地区、云南所有地区、北京和天津南部地区以及河北(邯郸、石家庄地区)。我国2017年松材线虫病疫区增加了大连市沙河口区(国家林业局，2017)，成为辽宁省第一个发现松材线虫病的地区。本研究中亦将辽宁省大部列为松材线虫无胁迫条件下的低风险分布区。为此，控制松材线虫病的扩散要加大植物检疫力度，且考虑到辽宁为风险分布区，以及广布寄主树种，需要在辽宁地区严防其扩散传播，以防威胁东三省的生态安全。
程功, 吕全, 冯益明, 理永霞, 王艳丽, 张星耀, 2015. 气候变化背景下松材线虫在中国分布的时空变化预测. 林业科学, 51(6): 119-126.
陈守常, 2010. 松材线虫病病原与致病机理研究进展. 四川林业科技, 31(1): 18-25.陈燕婷， 尤民生， 2014. 基于CLIMEX和ArcGIS的湿地松粉蚧在中国的适生性分析. 环境昆虫学报， 36(5)： 697-704.
国家林业局， 2016. 国家林业局公告(2016年第6号). (2016-04-14)[2017-02-10]. http:∥www.forestry.gov.cn/portal/main/govfile/13/govfile_2251.htm.
国家林业局， 2017. 公告： 国家林业局公告(2017年第4号)(2017年松材线虫病疫区公告). (2017-01-20)[2017-02-10]. http:∥www.forestry.gov.cn/main/4461/content-941973.html.
韩阳阳, 王焱, 项杨, 叶建仁, 2015. 基于Maxent生态位模型的松材线虫在中国的适生区预测分析. 南京林业大学学报(自然科学版), 39(1): 6-10.
何善勇, 温俊宝, 骆有庆, 宗世祥, 赵宇翔, 韩婧, 2012. 气候变暖情境下松材线虫在我国的适生区范围. 应用昆虫学报, 49(1): 236-243.
黄瑞芬, 2015. 中国不同地理种群松材线虫的耐寒性研究. 硕士学位论文. 北京: 北京林业大学.
潘沧桑, 2011. 松材线虫病研究进展. 厦门大学学报(自然科学版), 50(2): 476-483.
潘红伟, 2009. 松材线虫(Bursaphelenchusxylophilus)在我国的潜在分布区研究. 博士学位论文. 北京: 中国林业科学研究院.
宋红敏, 张清芬, 韩雪梅， 徐岩， 徐汝梅， 2004. CLIMEX: 预测物种分布区的软件. 昆虫知识, 41(4): 379-386.
王涛, 任艳玲, 何善勇, 张润志， 杨茂发, 2015 . 基于CLIMEX和ArcGIS的橘小实蝇在贵州适生性研究. 环境昆虫学报, 37(3): 534-542.
徐亮, 2013. 马尾松响应松材线虫侵染的CaM和CaMBP钙信号转导特征. 博士学位论文. 北京: 北京林业大学.
张海军, 2007. 中国松材线虫病潜在分布区研究. 硕士学位论文. 呼和浩特: 内蒙古农业大学.
张伟, 2014. 松材线虫繁殖力和致病性的差异形成及其分子基础. 硕士学位论文. 合肥: 安徽大学.
朱丽华, 季锦衣, 张文莉, 叶建仁, 2011. 松材线虫的低温冷冻保存研究. 南京林业大学学报, 35(6): 63-66.
CHEJARA V K, KRITICOS D J, KRISTIANSEN P, SINDEL B M, WHALLEY R D B, NADOLNY C, 2010. The current and future potential geographical distribution ofHyparrheniahirta.WeedResearch, 50(2): 174-184.
JAGDALE G B, GREWAL P S, 2003. Acclimation of entomopathogenic nematodes to novel temperatures: trehalose accumulation and the acquisition of the thermotolerance.InternationalJournalforParasitology, 33(2): 145-152.
KAZUYOSHI F T, 2013. Pine wood nematode,Bursaphelenchusxylophilus.AnnualReviewofPhytopathology, 51(2): 61-83.
KIKUCHI T, COTTON J A, DALZELL J J, HASEGAWA K, KANZAKI N, 2011. Genomic insights into the origin of parasitism in the emerging plant pathogenBursaphelenchusxylophilus.PlosPathogens, 7(9): e1002219.
LEDUC D, WHARTON D A, 2008. Three new species of free-living nematodes from inter-tidal sedimentsin southern New Zealand.Nematology, 10(5): 743-755.
SHI J, CHEN F, LUO Y Q, WANG Z, XIE B Y, 2013. First isolation of pine wood nematode fromPinustabuliformisforests in China.ForestPathology, 43: 59-66.
WHARTON D A, BLOCK W, 1993. Freezing tolerance in some antarctic nematode.FunctionalEcology, 7(5): 578-584.
(责任编辑:郭莹)
Potential distribution ofBursaphelenchusxylophilusin China due to adaptation cold conditions
ZHAO Jie, HAN Xiao, SHI Juan*
ForestryCollege,BeijingForestryUniversity,Beijing100083，China
【Aim】Bursaphelenchusxylophilus(pine wood nematode) is an important invasive species, causing huge economic losses to China. According to a recent survey ofB.xylophilusin the country, pine wilt disease occurred in high altitude where mean annual temperature is around 8 ℃. This suggests that its range probably moved to north and high altitude area. In this paper, we discuss the potential distribution ofB.xylophilusin the situation of adaptation cold conditions. 【Method】 The experiment included two treatments where nematodes were exposed to low temperatures (3, 1, 0, -1, -3, -5, -6, -7, -9, -12 ℃) for 24 h (short) or 30 d (long) periods. The lethal temperature (LT50, LT75, LT99) was then measured. The potential distribution ofB.xylophiluswas built using CLIMEX and ArcGIS and based on LT50, LT75, LT99results. CLIMEX was used to analyze the adaptive capacity ofB.xylophilusand the results were integrated into ArcGIS to calculate the inverse weighted distance, based on EI (ecological index) value. The trend of the Northern Line of pine wood nematode distribution at moderate risk with or without low temperature stress was analyzed. 【Result】 Considering the long-term (30 d) cold stress conditions, the distribution is significantly farther north than under short-term (24 h) cold stress and non-stress (based on EI value) conditions. 【Conclusion】 The potential distribution of the pine wood nematode can move northward over time under low temperature stress, suggesting pine wood nematode capacity to adapt to cold conditions.
Bursaphelenchusxylophilus; ArcGIS; CLIMEX; risk analysis; cold adaptation
2017-02-17 接受日期(Accepted)： 2017-04-06
中央高校基本科研业务费专项(2016ZCQ07)； 林业公益性行业科研专项(201504304)
赵捷， 女， 硕士研究生。 研究方向: 植物检疫。 E-mail: zhao_jie_symphony@163.com
*通信作者(Author for correspondence), E-mail: shi_juan@263.net
10.3969/j.issn.2095-1787.2017.03.003
