下载:
-
叶面积指数(LAI)为投影到水平基准面上的每单位地表面积上总绿叶面积的一半，是描述植被生长状况和开展陆面过程研究的重要参数[1]。传统获取LAI数据的地面测量方法只适合小范围区域的短期监测，而卫星遥感具有大面积和重复观测的特点，为获取大尺度长期LAI提供了有效的数据来源[2]。目前已经有多个全球范围、长时间序列LAI产品，如AVHRR LAI[3]、MODIS LAI[4]、GEOV2 LAI[5]、PROBA-V LAI[6]、VIIRS LAI[7]、GLASS LAI[8]、MUSES LAI[9]等。然而这些产品的空间分辨率大多较低，在反演过程中受到云雪覆盖、气溶胶、大气辐射以及传感器自身参数等因素影响，导致数据存在算法误差，数据集存在异常值、数据缺失和精度偏低等问题，限制了产品的应用研究[10]。
多源数据同化融合技术可以显著提高数据精度。基于集合卡尔曼滤波(EnKF)的数据同化算法具有整合不同来源和不同误差信息数据的能力，通过观测数据对模型模拟结果进行调整和校正，提高模拟数据精度[11]。LIU等[12]整合3种不同来源数据的时间、光谱和角度信息，利用EnKF技术将反演的LAI进行气候学的递归更新，可以填补缺失数据，获得更准确的LAI产品。李喜佳等[13]耦合辐射传输模型和动态模型，运用双集合卡尔曼滤波分别对林地、草地和农田的LAI进行同化，改善了模型预测能力。周红敏等[14]利用EnKF方法耦合神经网络模型和动态模型，将反演得到的Landsat高分辨率LAI数据与MODIS LAI融合生成30 m空间分辨率的LAI数据集，估算结果与地面实测数据一致性较高。然而，LAI产品在反演过程中输入数据和模型算法的不确定性导致其存在质量问题，如MODIS产品由主算法反演得到的高质量LAI在全球尺度上约占67%，其余为备用算法反演得到的低质量LAI[15]。主算法失效主要是由于残云或较差的大气校正，失效后采用备用算法，导致备用算法检索数据准确性较低[16]。
为了尽可能减少LAI产品由于观测不稳定造成的不确定性，本研究将多年LAI产品数据作为先验知识，利用产品的质量控制文件构建LAI背景库修正其低质量数据，利用混合像元分解方法对2017年的MODIS、VIIRS和PROBA-V的LAI数据进行降尺度处理，实现LAI时序数据的重建工作。利用3种LAI数据集和Sentinel-2反射率数据，结合EnKF技术耦合LAI动态模型和辐射传输模型生成同化LAI，最后对3种同化后的LAI数据进行赋权融合，达到减少不同来源数据差异、提高LAI产品估算精度和一致性的目的。
-
研究区位于黑龙江省尚志市帽儿山实验林场(45°14′~45°29′ N，127°29′~127°44′ E，图1)，平均海拔为428.0 m，属于温带季风气候，年均气温为2.8 ℃，土壤为典型暗棕壤。植被属于长白山植物区系，主要针叶树种有红松Pinus koraiensis、落叶松Larix gmelinii、樟子松Pinus sylvestris var. mongolica等，主要阔叶树种有白桦Betula platyphylla、山杨Populusdavidiana、色木槭Acer mono、水曲柳Fraxinus mandshurica、胡桃楸Juglans mandshurica等[17]。
图 1 研究区和样地分布示意图
Figure 1. Geographical location and plot distribution of the study area
-
在研究区内根据植被类型选取2块1 km×1 km大样地(图1)，每个大样地内包含5块100 m×100 m小样地。使用LAI-2200和TRAC冠层分析仪测量LAI[18−19]。使用GPS记录位置坐标。
-
考虑到研究区数据的可获取性，本研究选取了3种应用广泛的LAI产品(表1)。利用HEG软件，对MOD15A2H和VNP15A2H数据进行格式转换和投影变换。利用ArcGIS 10.2对CGLS 300m V1数据进行投影变换。
表 1 LAI产品数据及基本信息
Table 1. LAI product data and basic information
产品名称 传感器 空间分辨率/m 时间分辨率/d 投影方式 数据来源 MOD15A2H MODIS 500 8 正弦曲线 https://ladsweb.modaps.eosdis.nasa.gov/ VNP15A2H VIIRS 500 8 正弦曲线 https://ladsweb.modaps.eosdis.nasa.gov/ CGLS 300m V1 PROBA-V 300 10 简易圆柱 https://land.copernicus.eu/global/ -
选取2017年云层覆盖率低于10%的11幅Sentinel-2影像，利用Sen2cor工具进行大气校正得到L2A级数据，利用SNAP工具箱将数据重采样至10 m分辨率，并转换为与LAI产品一致的UTM WGS84投影。
-
采用2017年8月14日的Sentinel-2 MSIL2A影像，利用ENVI 5.3进行最大似然法分类[20]。定义了农田、阔叶林、针叶林、草地和其他5类训练样本，通过实地调查数据和30 m地表覆盖分类产品GLC_FCS30相结合的方式目视判读各样本点的类别，采用混淆矩阵的方法进行精度检验，计算得到的总体精度为79.34%，Kappa系数为0.7517，整体分类质量高。将分类后的影像掩膜提取生成10 m分辨率的阔叶林和针叶林土地覆盖图(图1)。
-
利用2014—2016年的MOD15A2H和VNP15A2H的质量控制数据集进行像元质量标识，采用时间加权插值算法建立两者的LAI背景库，基于背景库对2017年的MODIS LAI和VIIRS LAI进行数据重建，计算公式参考文献[15]。最后使用三次样条帽盖算法(LACC)[21]对初修正的2个数据集和PROBA-V LAI数据集进行平滑处理，去除异常值。
根据解混理论，假设每个粗分辨率像元的反射率是构成混合像元的每个土地覆盖类型反射率的线性组合，同理混合像元的LAI由不同类别的LAI及其在像元中所占的面积比例(即丰度)构成，基于Sentinel-2的10 m土地覆盖类型图来计算丰度，利用约束最小二乘法计算混合像元内不同类别的平均LAI值，并根据土地覆盖类型图分配到对应类别的像元上，从而得到10 m分辨率的高精度LAI数据[22]。
-
利用降尺度后的3种LAI数据集分别构建经验性的LAI动态模型，将模型预测的LAI作为参数，连同其他参数一起输入PROSAIL模型模拟森林冠层反射率。同时，将Sentinel-2的红光波段(Band 4)和近红外波段(Band 8)影像作为观测数据，通过EnKF算法比较观测数据集合和模拟反射率集合，对预测LAI进行校正。结合精度评价结果对同化后的3种LAI数据进行赋权融合。
-
LAI动态模型能模拟叶面积指数[23]，描述LAI随时间的动态变化过程，计算公式参考文献[24]。PROSAIL模型是由叶片光学特性模型(PROSPECT)和冠层二向反射率模型(SAIL)组合成的耦合模型，参数设置参考文献[25−27]。PROSPECT模型输入参数为：叶片结构参数(N=1.4)、叶绿素含量(Cab=58 μg·cm−2)、叶片含水量(Cw=0.025 g·cm−2)和干物质含量(Cm=0.009 g·cm−2)；SAIL模型输入参数为：叶面积指数(LAI)﹑土壤亮度参数(Psoil=0.25)、热点参数(Hspot=0.037)、太阳天顶角(tts=30)、观测天顶角(tto=0)和相对方位角(psi=0)，叶片反射率和透射率采用PROSPECT模型输出数据，tts、tto和psi可通过遥感影像头文件获取。
-
EnKF是基于蒙特卡罗的集合预报方法估计预报误差统计量的数据同化技术[11]，计算公式如下：
$$ {{\boldsymbol{A}}}_{t}^{\mathrm{f}}={\boldsymbol{M}}\left({{\boldsymbol{A}}}_{t-1}^{\mathrm{a}}\right)+\varepsilon \text{；} $$ (1) $$ {{\boldsymbol{A}}}_{t}^{\mathrm{a}}={{\boldsymbol{A}}}_{t}^{\mathrm{f}}+{{\boldsymbol{P}}}_{t}{{\boldsymbol{H}}}^{\mathrm{T}}({\boldsymbol{H}}{{\boldsymbol{P}}}_{t}{{\boldsymbol{H}}}^{\mathrm{T}}+{{\boldsymbol{R}}}_{t}{)}^{-1}({\boldsymbol{D}}-{\boldsymbol{H}}{{\boldsymbol{A}}}_{t}^{\mathrm{f}}) \text{；} $$ (2) $$ {{\boldsymbol{P}}}_{t}=\frac{1}{N-1}({{\boldsymbol{A}}}_{t}^{\mathrm{f}}-{\overline{{{\boldsymbol{A}}}_{t}^{\mathrm{f}}}})({{\boldsymbol{A}}}_{t}^{\mathrm{f}}-{\overline{{\boldsymbol{A}}_{t}^{\mathrm{f}}}}{)}^{\mathrm{T}} \text{；} $$ (3) $${\overline{{{\boldsymbol{A}}}_{t}^{\mathrm{a}}}}=\frac{1}{N}\sum _{i=1}^{N}{{\boldsymbol{A}}}_{{{it}}}^{\mathrm{a}} 。 $$ (4) 式(1)~(4)中：${\boldsymbol{A}}=({x}_{1}, {x}_{2}, \cdots , {x}_{N})\in {{{{{\boldsymbol{R}}}}}}^{n\times N}$，是包含N个成员的状态变量集合，x表示n维状态向量(LAI)；M代表LAI动态模型，$ \varepsilon $为模型误差；${{\boldsymbol{A}}}_{t}^{\mathrm{f}}$表示t时刻模型预测的状态变量集合；a表示分析状态变量，T表示转置矩阵。${\boldsymbol{D}}\in {{{{{\boldsymbol{R}}}}}}^{m\times N}$代表扰动观测集合矩阵，是通过给定的m维观测向量${\boldsymbol{d}}\in {{{{{\boldsymbol{R}}}}}}^{m}$加上N个扰动构成，本研究中m=2。${\boldsymbol{P}}\in {{{{{\boldsymbol{R}}}}}}^{n\times n}$代表预测误差协方差矩阵，${\boldsymbol{R}}\in {{{{{\boldsymbol{R}}}}}}^{m\times m}$代表观测误差协方差矩阵，其算法与P相似，如式(3)所示。式(4)中，$\overline{{{\boldsymbol{A}}}_{t}^{\mathrm{a}}}$为状态最优估计值。H代表PROSAIL模型，因为该模型是非线性观测算子，不适用于标准EnKF，所以构建了扩展集$\widehat{{\boldsymbol{A}}}=[{{\boldsymbol{A}}}_{t}^{\mathrm{f}},\; h({{\boldsymbol{A}}}_{t}^{\mathrm{f}}\left)\right]$，将H转化成线性算子$\widehat{{\boldsymbol{H}}}$，此时P为扩展集协方差矩阵，故式(2)改写如下：
$$ {{\boldsymbol{A}}}_{t}^{\mathrm{a}}={{\boldsymbol{A}}}_{t}^{\mathrm{f}}+{\boldsymbol{P}}{\widehat{{\boldsymbol{H}}}}^{\mathrm{T}}(\widehat{{\boldsymbol{H}}}{\boldsymbol{P}}{\widehat{{\boldsymbol{H}}}}^{\mathrm{T}}+{{\boldsymbol{R}}}_{t}{)}^{-1}({\boldsymbol{D}}-\widehat{{\boldsymbol{H}}}{{\boldsymbol{A}}}_{t}^{\mathrm{f}})。 $$ (5) 本研究观测集合和第1次的LAI预测集合是通过添加N个扰动(符合均值为0，方差为1的标准正态分布的随机数)生成，其中集合数N取值为100。
-
对3种同化LAI数据的精度检验结果归一化处理后分别赋予权重，将3种同化LAI数据分别与各自的权重相乘得到融合LAI数据。除了相关系数外的指标数值的绝对值越大，精度越差，需要在归一化处理前对它们取倒数，计算公式如下：
$$ {Y}_{i,j}=\left|1/{Y}_{i,j}\right|,j=\mathrm{1, \;2},\; \cdots, \;n-1 \text{；} $$ (6) $$ {Y}_{i}=\frac{1}{n}\sum _{j=1}^{n}\left({Y}_{i,j}/\sum _{i=1}^{m}{Y}_{i,j}\right) 。 $$ (7) 式(6)~(7)中：n表示评价指标的个数，m表示LAI同化数据的个数，i=1$, \cdots , $ m，$ {Y}_{i,j} $表示第i个LAI同化数据对应的第j个指标的值，$ {Y}_{i} $表示第i个LAI同化数据在所有指标进行归一化后的值。
-
利用实测数据对3种同化LAI数据和融合LAI数据进行精度评价，评价指标选用相关系数(R)、均方根误差(ERMSE)、平均绝对误差(EMAE)和相对偏差(ERE)，计算公式参考文献[28]。
-
由图2可见：原始MODIS和VIIRS LAI时序曲线存在较为剧烈的跳跃现象，影响了其对植被物候特征的表征，而PROBA-V LAI时序曲线较为平滑，非生长季出现较小的跳跃现象。各产品都清晰表现了植被LAI夏季高、冬季低的生长特点。在植被LAI关键点缺失的情况下，MODIS和VIIRS产品基于背景库修正的方法，能够更好地重建LAI时序曲线，修正后的LAI与PROBA-V LAI差异较小，且时序曲线在形态上与PROBA-V也更加接近。修正后的LAI使用LACC算法进行平滑处理，平滑后的非生长季LAI均值都有所提高。在生长季的阔叶林VIIRS和PROBA-V LAI变化趋势相似；在生长季的针叶林MODIS和PROBA-V LAI变化趋势相似。在生长季，PROBA-V阔叶林的LAI高于针叶林。总体上，PROBA-V LAI均值高于MODIS和VIIRS。
图 2 3种LAI时序曲线在阔叶林和针叶林区域的变化情况以及重建效果
Figure 2. Changes and reconstruction effects of 3 LAI time series curves in broad-leaved forest and coniferous forest areas
-
由表2可见：3种LAI产品与实测数据的相关系数均较低，均方根误差、平均绝对误差和相对偏差值均较高，说明3种LAI产品与实测数据之间的一致性不高，没有明显的线性相关关系。经过EnKF同化后，3种LAI产品的相关系数均提升到0.49以上，MODIS LAI、VIIRS LAI和PROBA-V LAI的相关系数在阔叶林分别提高了0.57、0.52和0.57，在针叶林分别提高了0.52、0.32和0.40。3种同化后LAI数据的均方根误差均在0.40以下。其中，PROBA-V LAI在阔叶林和针叶林的均方根误差最小。VIIRS LAI在阔叶林的均方根误差小于MODIS LAI，而MODIS LAI在针叶林的均方根误差小于VIIRS LAI。相比同化前，VIIRS LAI在针叶林的均方根误差减小程度最低，仅减小了0.67，而VIIRS LAI在阔叶林的均方根误差减小程度最高，减小了1.69。对比3种同化后LAI数据的平均绝对误差，MODIS和VIIRS LAI在阔叶林的平均绝对误差最小，均为0.60，其中，VIIRS LAI在阔叶林的平均绝对误差相比同化前减小了1.11，减小程度最高。对比相对偏差，同化后除了MODIS LAI在阔叶林的相对偏差为负值，其他相对偏差都为正值。相比同化前，除了MODIS LAI在针叶林的偏差增大(同化前相对偏差=−0.72%，同化后相对偏差=8.34%)，其他数据的偏差均有不同程度的降低。
表 2 3种LAI产品在阔叶林和针叶林数据同化前后的精度指标
Table 2. Accuracy indexes of 3 LAI products before and after data assimilation in broad-leaved and coniferous plots
LAI产品 样地 同化前 同化后 R ERMSE EMAE ERE/% R ERMSE EMAE ERE/% MODIS 阔叶林 0.02 1.60 1.37 −23.75 0.59 0.37 0.60 −2.56 针叶林 0.07 1.46 1.30 −0.72 0.59 0.24 0.64 8.34 VIIRS 阔叶林 0.04 2.00 1.71 −28.62 0.56 0.31 0.60 2.29 针叶林 0.17 0.95 0.77 11.55 0.49 0.28 0.72 10.53 PROBA-V 阔叶林 0.05 1.20 1.04 15.24 0.62 0.14 0.86 12.49 针叶林 0.16 1.54 1.39 28.39 0.56 0.19 1.03 20.53 由图3可见：同化后LAI数据在阔叶林和针叶林从大到小依次为PROBA-V、VIIRS、MODIS，PROBA-V LAI具有明显的高估现象，与表2的相对偏差情况一致，与PROBA-V LAI原始数据就具有高估情况有关。同化后MODIS和VIIRS LAI数值相近，与两者原始数据就具有较高一致性有关。
图 3 同化LAI与实测值的散点验证图
Figure 3. Scatter verification diagram of assimilated LAI and measured value
-
由图4可见：相比同化后LAI，融合后LAI与实测数据的相关性较高。融合LAI在阔叶林和针叶林的相关系数分别为0.83和0.76，均方根误差均在0.20以下，平均绝对误差和相对偏差趋势相同。在阔叶林，相比同化后的MODIS、VIIRS和PROBA-V LAI，融合后的LAI相关系数分别提高了0.24、0.27和0.21；对比同化后的MODIS和VIIRS LAI，融合后的LAI均方根误差分别降低了0.22和0.16；融合后的LAI平均绝对误差相比同化后的降低了0.02~0.28；融合后的LAI 相对偏差相比同化后的PROBA-V LAI，偏差降低程度最大，降低了8.59%。在针叶林，相比同化后的MODIS、VIIRS和PROBA-V LAI，融合后的LAI相关系数分别提高了0.17、0.27和0.20；均方根误差分别降低了0.11、0.15和0.06；对比同化后的PROBA-V LAI，融合后的LAI平均绝对误差降低了0.29；融合后的LAI相对偏差相比同化后的PROBA-V LAI，偏差降低程度最大，降低了9.88%。
图 4 融合LAI与实测值的散点验证图
Figure 4. Scatter verification diagram of fused LAI and measured values
-
一般而言，LAI时序曲线是光滑且连续的，因此，当严重锯齿状的变化曲线中出现与正常趋势差距较大的拐点时，认为不符合植被正常生长规律[29]。MODIS和VIIRS LAI时序曲线存在剧烈的跳跃现象，PROBA-V LAI即使利用时间平滑和间隙填充技术，使得时序曲线较为平滑，但仍存在较多缺失值[30]，因此在产品应用前进行重建至关重要。本研究对3种LAI产品进行重建，初步改善了数据集质量较差的问题，使得三者LAI时序曲线变化趋势趋于一致。由于LAI背景库值反映的是研究区LAI生长变化的平均趋势，受空间和时间因素的影响非常大，因此与实际值有一定的偏差。本研究利用混合像元分解法进行降尺度处理，由于粗分辨率像元存在异质性，会存在一定的不确定性，经过降尺度处理，样地内同质像元具有相同的LAI值。本研究表明：经过同化方法校正了LAI值，使得LAI产品与实测值具有较高的相关性，但存在低值高估和高值低估的现象。经过EnKF同化，3种LAI产品精度和一致性均有所提高，经过同化融合后的LAI数据，降低了单一LAI数据误差，可有效弥补地面实测数据缺乏的不足，但能否满足植被监测和地表模型输入参数的要求还需进一步验证。本研究通过将LAI数据降尺度后作为同化算法的输入数据，使融合LAI的精度得以提高，生成10 m空间分辨率LAI数据，然而由于缺少研究区长期有效的观测数据，生成数据的时间分辨率受限制，在后续研究中可以考虑加入其他反射率数据，以期得到高精度长时序的LAI产品。
-
本研究结果表明：3种原始LAI产品与实测数据的一致性均较差，精度较差；EnKF同化后的LAI更接近实测数据；相比针叶林，阔叶林LAI同化效果较好；总体上MODIS LAI同化效果最好，其次是PROBA-V LAI，VIIRS LAI在阔叶林和针叶林的同化效果存在较大差异。融合后的LAI精度优于同化后单一的LAI产品，其相关系数为0.83和0.76，比同化后LAI提升了0.17~0.27。同化融合后的LAI数据提高了单一LAI产品的精度，减小了不同源数据的差异，达到了充分利用多源遥感数据优势和提高产品精度的目的，可以为开展高精度LAI数据生成方法的研究提供参考。
Reconstruction and fusion correction method of multi-source LAI products in Maoershan forest based on ensemble Kalman filter
-
摘要:
目的 现有叶面积指数(LAI)产品大多存在分辨率低、数据异常和精度低等问题，难以满足某些应用需求。因此，本研究提出一种多源LAI数据的融合方法，以减少不同来源数据的差异并提高产品精度。 方法 以帽儿山实验林场的阔叶林和针叶林区域为研究区，基于2017年的MODIS、VIIRS和PROBA-V 的LAI产品，利用多年LAI数据作为先验知识建立LAI背景库修正低质量数据，对3种LAI数据集进行混合像元分解的降尺度处理，基于Sentinel-2反射率产品耦合集合卡尔曼滤波(EnKF)算法、LAI动态模型和辐射传输模型进行数据同化，最后对同化后的3种LAI数据进行赋权融合，使用实测数据进行精度评价。 结果 在阔叶林，同化后的MODIS、VIIRS和PROBA-V LAI与实测数据的相关系数分别为0.59、0.56和0.62，比原始数据提升了0.57、0.52和0.57；均方根误差分别为0.37、0.31和0.14，比原始数据减小了1.23、1.69和1.06。在针叶林，同化后的MODIS、VIIRS和PROBA-V LAI与实测数据的相关系数分别为0.59、0.49和0.56，比原始数据提升了0.52、0.30和0.40；均方根误差分别为0.24、0.28和0.19，比原始数据减小了1.22、0.67和1.35。通过融合方法，阔叶林LAI和针叶林LAI的相关系数分别为0.83和0.76，比同化后数据的相关性更高；均方根误差分别为0.15和0.13，比同化后数据的误差更小。 结论 通过数据同化提升了3种LAI产品精度，融合后LAI较同化后单一LAI具有更高的精度和可靠性。图4表2参30 -
关键词:
- 叶面积指数(LAI) /
- MODIS /
- VIIRS /
- PROBA-V /
- 重建 /
- 集合卡尔曼滤波(EnKF) /
- 数据融合
Abstract:Objective Most of the existing leaf area index (LAI) products have some problems, such as low resolution, abnormal data and low accuracy, which are difficult to meet the requirements of some applications. Therefore, this study proposes a method of fusing multi-source LAI data to reduce the differences of data from different sources and improve product accuracy. Method The broad-leaved forest and coniferous forest in Maoershan experimental forest farm were taken as the research area. Based on MODIS LAI, VIIRS LAI and PROBA-V LAI products in 2017, the LAI background database was established to correct low-quality data by using years of LAI data as prior knowledge, and 3 LAI data sets were downscaled by mixed pixel decomposition. Based on Sentinel-2 reflectivity product coupling ensemble Kalman filter (EnKF) algorithm, LAI dynamic model and radiative transfer model, data assimilation was carried out. Finally, 3 LAI data after assimilation were weighted and fused, and the accuracy was evaluated by using measured data. Result In broad-leaved forest, the correlation coefficients between the assimilated MODIS, VIIRS and PROBA-V LAI and the measured data were 0.59, 0.56 and 0.62, respectively, which were 0.57, 0.52 and 0.57 higher than the original data. The root mean square error (ERMSE) were 0.37, 0.31 and 0.14 respectively, which were 1.23, 1.69 and 1.06 lower than the original data. In coniferous forest, the correlation coefficients between the assimilated MODIS, VIIRS and PROBA-V LAI and the measured data were 0.59, 0.49 and 0.56, respectively, which were 0.52, 0.30 and 0.40 higher than the original data. ERMSE were 0.24, 0.28 and 0.19 respectively, which were 1.22, 0.67 and 1.35 lower than the original data. Through the fusion method, the correlation coefficients of LAI in broad-leaved forest and coniferous forest were 0.83 and 0.76 respectively, which were higher than the data after assimilation. ERMSE were 0.15 and 0.13, respectively, which were smaller than the error of the assimilated data. Conclusion Through data assimilation, the accuracy of 3 LAI products is improved, and the fused LAI data has higher accuracy and reliability than the single LAI data after assimilation. [Ch, 4 fig. 2 tab. 30 ref.] -
Key words:
- leaf area index (LAI) /
- MODIS /
- VIIRS /
- PROBA-V /
- reconstruction /
- ensemble Kalman filter (EnKF) /
- data fusion
-
密码子是识别和传递生物体遗传信息、联系蛋白质与DNA之间的重要桥梁，在生物体遗传和变异中起着至关重要的作用[1]。编码同一氨基酸的不同密码子被称为同义密码子。由于基因突变和自然选择的影响，某些同义密码子在蛋白质翻译过程中往往被高频使用，被称为密码子的使用偏好性[2−3]。物种的生物学功能与密码子偏好性密切相关，密码子偏好性不仅可以影响生物编码基因的蛋白质合成速率和翻译速率[4]，还会影响蛋白质结构、折叠程度和mRNA的合成[5]。研究表明：同一物种或亲缘关系相近的物种，具有相似的密码子偏好使用模式[6]，通过分析物种的密码子偏好性可以衡量物种之间的基因表达量，进而探究物种之间亲属关系[7]。通过密码子偏好性的研究，能够更好地阐明物种进化过程中基因的表达规律[8]，为利用基因工程技术改良物种目标基因提供参考依据[9]。
梁山慈竹Dendrocalamus farinosus属竹亚科Bambusoideae牡竹属Dendrocalamus，又名大叶竹和瓦灰竹，是中国西南地区重要的经济竹种[10]，生长速度快，适应性强，竹笋效益高，属于优良的笋竹两用竹种，与硬头黄竹Bambusa rigida都属于竹编和制浆造纸的优质原料[11]。针对梁山慈竹叶绿体基因组密码子使用偏好性的研究鲜见报道。为了更好地挖掘和利用梁山慈竹的潜在经济价值，本研究以梁山慈竹叶绿体基因组序列为研究对象，分析其密码子偏好性使用模式，探究并总结其相关表达基因的密码子偏好性，以期分析影响梁山慈竹叶绿体基因组密码子偏好性的主要因素，并筛选出最优密码子，为后续梁山慈竹叶绿体基因工程改造等研究提供理论基础。
1. 材料与方法
1.1 叶绿体基因组序列的获取
根据GenBank登录号MZ681865.156在美国国家生物技术信息中心(NCBI)数据库中搜索并下载梁山慈竹叶绿体基因组序列，共有85条编码序列(CDS)。序列重复或小于300 bp会对密码子偏好性指标的测定产生影响[12]。对基因序列进行筛选，剔除序列长度小于300 bp且重复的序列，获取起始密码子为ATG，终止密码子为TAG、TGA和TAA的序列，最终获得51条CDS序列作为后续分析的样本序列。
1.2 密码子组成分析
运用CodonW1.4.2 (
http://sourceforge.net/projects/codonw )和EMBOSS (http://imed.med.ucm.es/EMBOSS/ )计算有效密码子数(ENC)、适应指数(CAI)、密码子偏性指数(CBI)、最优密码子频率(FOP)以及密码子第3位核苷酸A、T、C、G的含量(分别记为A3、T3、C3、G3)。利用ENC判断密码子偏好性程度，ENC＞35说明密码子偏好性比较弱；反之，说明偏好性强[13]。通过CUSP软件分析并获得密码子鸟嘌呤(G)和胞嘧啶(C)所占的比率(GC比率)及GC平均比率(GCall)，使用SPSS 25.0软件对梁山慈竹密码子各位置的GC比率与ENC进行相关分析。1.3 相对同义密码子使用度分析
运用CodonW 1.4.2对同义密码子相对使用度(RSCU)进行分析，即该密码子的实际使用频率与其理论使用频率的比值[14]。当RSCU大于1时，同义密码子中偏好使用该密码子，被称为高频密码子；当RSCU等于1时，密码子无偏好性；当RSCU小于1时，密码子使用偏好性较弱[15]。
1.4 中性绘图分析
中性绘图分析是对影响密码子使用偏好性的关键因素进行分析，X轴为GC3，Y轴为GC1和GC2的平均值，绘制二维散点图对GC3和GC12 (各基因 GC1和GC2的平均值)的相关性进行分析(GC1、GC2、GC3分别代表第1、2、3位密码子的GC比例)。若回归系数接近1，代表GC3和GC12显著相关，碱基组成没有差异，说明突变是决定密码子偏好性的主要因素；若回归系数接近0，则代表自然选择是主要因素。
1.5 ENC-plot绘图分析
ENC-plot绘图分析表现密码子的使用偏好性受到突变和自然选择的影响程度。使用Python 3.7进行ENC-plot绘图分析，构建散点图，横纵坐标分别为GC3、ENC，并绘制ENC的标准曲线。基因位点靠近或在标准曲线上，表明突变是决定密码子偏好性的主要因素，若基因位点和标准曲线距离很大，则说明偏好性主要由自然选择决定。
1.6 PR2-plot偏倚分析
PR2-plot分析表明基因中密码子的第3位碱基的构成情况。计算密码子碱基中第3位上4种碱基A、T、C、G比例，G3/(G3+C3)为X轴，A3/(A3+T3)为Y轴，绘制PR2-plot散点图，中心点为碱基比例A=T、C=G时的值，代表处于此区域的密码子并无使用偏好性[16]。
1.7 最优密码子分析
将51条基因升序排列后的ENC前后两端10％的基因建立高、低表达基因库。通过CodonW软件计算2个表达库中密码子的RSCU和ΔRSCU，同时满足高频密码子(RSCU＞1)和高表达密码子(ΔRSCU≥0.08)的为最优密码子[17]。
1.8 梁山慈竹和其他几种生物密码子偏好性比较分析
在Codon Usage Database (
http://www.kazusa.or.jp/codon/ )下载异源表达宿主和植物代表类群，包括巨龙竹D. farinosus、粉麻竹D. sinicus、小叶龙竹D. pulverulentus、硬头黄竹、大肠埃希菌Escherichia coli、烟草Nicotiana tabacum、拟南芥Arabidopsis thaliana和酿酒酵母Saccharomyces cerevisiae等物种基因组密码子的使用频率，与梁山慈竹基因组密码子使用频率比值进行比较分析，当梁山慈竹密码子使用频率比其他生物的比值≥2.0或≤0.5时，说明该物种与梁山慈竹的同义密码子的使用偏好性差异较大，当比值不在上述范围内时，表明这2个物种对该密码子的偏好性较接近。1.9 对应分析
将叶绿体基因如表1所示进行功能分类，使用CodinW软件，选择对应分析计算样本中各个基因的RSCU，将分析结果分布在59维向量空间中，分析指标间的对应性。
表 1 梁山慈竹叶绿体基因结构分析Table 1 Structural analysis of the choroplast genome of D. farinosus基因分类 基因分组 基因名称 光合系统基因 光系统Ⅰ基因 psaA、psaB、psbA、psbC、psbD、psbB 光系统Ⅱ基因 petA、petB、petD 细胞色素b/f复合体基因 atpA、atpB、atpE、atpF、atpI 三磷酸腺苷合成酶基因 ndhA、ndhB、ndhC、ndhD、ndhE、ndhF、ndhG、ndhH、ndhI、ndhJ、ndhK 遗传系统基因 烟酰胺腺票吟二核甘酸氧化还原酶基因 rbcL 二磷酸核酮糖羧化酶大亚基基因 rpoA、rpoB、rpoC1、rpoC2 RNA聚合酶亚基基因 rps2、rps3、rps4、rps7、rps8、rps11、rps12、rps14、rps18 核糖体蛋白小亚基基因 rpl2、rpl14、rpl16、rpl20、rpl22 其他基因 成熟酶K基因 matK 膜蛋白基因 cemA 细胞色素合成基因 ccsA 酪蛋白分解蛋白酶基因 clpP 未知功能基因 假定叶绿体阅读框 ycf2、ycf3、infA 2. 结果与分析
2.1 密码子的碱基组成分析
分析梁山慈竹叶绿体基因组CDS序列的碱基组成：梁山慈竹的4种碱基所对应的同义密码子的第3位碱基比例 (T3s、A3s、C3s、G3s)分别为45.28%、42.07%、18.13%、17.96%，T3s和A3s远高于G3s和C3s，表明梁山慈竹叶绿体基因组密码子的第3位碱基以A/U结尾为主。梁山慈竹的ENC为50.40，CAI为16.6%，第3位同义密码子的GC比率 (GC3S)为28.1%，表明其叶绿体基因组密码子偏好性较弱。
梁山慈竹叶绿体基因组密码子的GC平均比率为39.48%，且GC1 (47.69%)＞GC2 (39.70%)＞GC3 (31.05%)。ENC为 39.04~61.00，均值为49.51，GC比率在基因密码子上并没有均匀分布(表2)。ENC和密码子3个位置GC比率的相关分析(表3)结果发现：ENC与GC3比率显著相关，与GC1、GC2不显著相关，说明密码子使用偏好性形成过程中GC3的影响作用大于GC1、GC2。
表 2 梁山慈竹叶绿体基因组各基因密码子相关参数统计Table 2 Statistics of codon related parameters of various genes in the chloroplast genome of D. farinosus基因 GC比率/% ENC CAI FOP 基因 GC比率/% ENC CAI FOP GC GC1 GC2 GC3 GC GC1 GC2 GC3 rps12 41.87 52.00 47.20 26.40 44.85 0.140 0.341 rps18 33.53 34.50 39.77 26.32 39.04 0.147 0.333 psbA 42.56 49.72 42.94 35.03 41.33 0.313 0.532 rpl20 36.11 38.33 40.83 29.17 50.97 0.112 0.298 matK 34.44 40.82 32.42 30.08 49.49 0.166 0.329 clpP 43.01 52.53 38.25 38.25 52.37 0.175 0.337 psbD 44.44 53.39 43.50 36.44 48.99 0.242 0.456 psbB 44.01 54.42 45.97 31.63 50.73 0.190 0.380 psbC 44.66 53.59 44.73 35.63 48.91 0.183 0.386 petB 41.06 48.93 41.20 33.05 47.31 0.191 0.333 rpoB 39.19 49.81 38.01 29.74 49.69 0.153 0.353 petD 40.37 50.93 39.13 31.06 49.46 0.161 0.305 rpoC1 39.87 49.93 38.07 31.63 52.77 0.156 0.347 rpoA 37.06 46.18 35.59 29.41 49.94 0.151 0.311 rpoC2 38.95 49.01 36.64 31.18 52.29 0.154 0.333 rps11 43.52 50.69 56.25 23.61 44.33 0.174 0.396 rps2 38.40 40.51 40.93 33.76 52.55 0.168 0.338 infA 40.35 43.86 35.96 41.23 61.00 0.181 0.409 atpI 38.84 47.58 36.29 32.66 50.55 0.163 0.353 rps8 36.50 41.61 41.61 26.28 46.62 0.122 0.374 atpF 38.27 47.62 35.45 31.75 53.17 0.147 0.353 rpl14 38.71 54.84 37.10 24.19 51.90 0.181 0.392 atpA 42.06 56.01 39.96 30.12 49.96 0.182 0.385 rpl16 44.76 52.14 53.57 28.57 39.41 0.115 0.354 rps14 39.42 39.42 46.15 32.69 41.73 0.135 0.384 rps3 33.47 43.75 31.67 25.00 48.03 0.193 0.402 psaB 41.81 48.71 43.13 33.61 49.34 0.172 0.350 rpl22 37.56 41.33 36.67 34.67 47.48 0.188 0.415 psaA 43.68 51.80 43.28 35.95 52.07 0.198 0.373 rpl2 44.56 51.77 48.58 33.33 53.33 0.143 0.361 ycf3 39.69 47.40 38.15 33.53 55.45 0.156 0.343 ndhB 38.16 42.07 39.33 33.07 46.71 0.156 0.348 rps4 37.13 47.52 37.13 26.73 49.59 0.169 0.386 rps7 39.49 49.68 45.22 23.57 48.31 0.164 0.373 ndhJ 39.38 49.38 36.88 31.88 51.48 0.176 0.356 ndhF 34.19 37.84 38.92 25.81 46.19 0.144 0.321 ndhK 38.60 41.70 43.72 30.36 51.91 0.159 0.329 ccsA 33.64 33.74 41.10 26.07 45.60 0.152 0.307 ndhC 39.67 50.41 36.36 32.33 48.75 0.177 0.345 ndhD 36.19 40.72 36.93 30.94 48.98 0.133 0.314 atpE 42.51 52.17 39.13 36.23 59.51 0.167 0.405 ndhE 33.33 41.18 32.35 26.47 59.06 0.144 0.316 atpB 42.62 53.91 41.68 32.26 47.43 0.192 0.381 ndhG 34.46 44.07 32.77 26.55 45.77 0.125 0.250 rbcL 44.14 57.11 43.93 31.38 50.19 0.271 0.454 ndhI 34.99 37.57 38.67 28.73 52.09 0.171 0.345 ycf4 41.22 48.39 39.78 35.48 47.14 0.162 0.385 ndhA 33.98 42.42 36.36 23.14 44.35 0.140 0.321 cemA 33.62 41.99 27.71 31.17 55.91 0.176 0.342 ndhH 37.82 50.76 34.77 27.92 49.95 0.155 0.322 petA 40.29 53.58 35.2 32.09 51.12 0.155 0.331 表 3 梁山慈竹叶绿体基因组中各基因参数的相关性分析Table 3 Correlation analysis of various gene parameters in the chloroplast genome of D. farinosus参数 GC1 GC2 GC3 ENC CAI CBI FOP GC3s GC GC1 1 GC2 0.300* 1 GC3 0.265 −0.009 1 ENC 0.142 −0.425** 0.389** 1 CAI 0.409** 0.076 0.370** 0.012 1 CBI 0.438** 0.272 0.322* −0.092 0.774** 1 FOP 0.402** 0.312* 0.341* −0.064 0.797** 0.965** 1 GC3s 0.271 −0.029 0.946** 0.445** 0.330* 0.330* 0.370** 1 GC 0.814** 0.673** 0.525** 0.010 0.407** 0.512** 0.518** 0.499** 1 说明： *表示显著相关 (P＜0.05)；**表示极显著相关 (P＜0.01)。 2.2 相对同义密码子使用度分析
梁山慈竹叶绿体基因组中共包含18110个密码子(表4)，总计编码20个氨基酸，密码子数为12~705个，其中密码子UGA共有12个，密码子含量最多的是编码谷氨酸的GAA，共有705个。梁山慈竹叶绿体基因组蛋白编码序列RSCU分析表明：氨基酸含量较高的有亮氨酸(Leu)和精氨酸(Arg)，均为6个密码子编码，编码精氨酸的是UUA、UUG、CUU、CUC、CUA和CUG；编码亮氨酸的有AGA、AGG、CGU、CGC、CGA和CGG；除此之外，蛋氨酸(Met)和色氨酸(Trp)均只有1个密码子编码，分别是AUG和UGG，其余氨基酸密码子编码个数分别为2~4个。
表 4 梁山慈竹叶绿体基因组蛋白编码序列RSCU分析Table 4 RSCU of protein coding region in the chloroplast of D. farinosus氨基酸 密码子 数量 RSCU 氨基酸 密码子 数量 RSCU 氨基酸 密码子 数量 RSCU 氨基酸 密码子 数量 RSCU Phe UUU* 644 1.29 Tyr UAU* 532 1.59 Ser UCU* 343 1.58 Cys UGU* 151 1.53 Phe UUC 351 0.71 Tyr UAC 137 0.41 Ser UCC* 260 1.19 Cys UGC 47 0.47 Leu UUA* 634 1.94 TER UAA* 28 1.56 Ser UCA* 222 1.02 Arg AGA* 322 1.75 Leu UUG* 362 1.11 TER UAG 14 0.78 Ser UCG 119 0.55 Arg AGG 119 0.64 Leu CUU* 420 1.29 TER UGA 12 0.67 Ser AGU* 273 1.25 Arg CGU* 261 1.41 Leu CUC 138 0.42 Trp UGG* 328 1.00 Ser AGC 89 0.41 Arg CGC 95 0.51 Leu CUA 295 0.90 Gln CAA* 477 1.53 Thr ACU* 403 1.68 Arg CGA* 234 1.27 Leu CUG 107 0.33 Gln CAG 148 0.47 Thr ACC 181 0.75 Arg CGG 76 0.41 Ile AUU* 740 1.48 Glu GAA* 705 1.46 Thr ACA* 259 1.08 Gly GGU* 421 1.24 Ile AUC 295 0.59 Glu GAG 263 0.54 Thr ACG 116 0.48 Gly GGC 145 0.43 Ile AUA 461 0.92 Lys AAA* 647 1.44 Ala GCU* 493 1.73 Gly GGA* 538 1.58 Met AUG* 416 1.00 Lys AAG 253 0.56 Ala GCC 172 0.60 Gly GGG 259 0.76 Val GUU* 382 1.47 Asp GAU* 522 1.54 Ala GCA* 343 1.20 Pro CCU* 286 1.48 Val GUC 126 0.49 Asp GAC 155 0.46 Ala GCG 135 0.47 Pro CCC* 196 1.01 Val GUA* 390 1.50 His CAU* 311 1.47 Asn AAU* 528 1.48 Pro CCA* 209 1.08 Val GUG 139 0.54 His CAC 112 0.53 Asn AAC 187 0.52 Pro CCG 84 0.43 说明：*表示RSCU大于1的高频密码子。 梁山慈竹叶绿体基因组RSCU大于1的密码子数目为34个(分别为UUU、UUA、UUG、CUU、AUU、AUG、GUU、GUA、UCU、UCC、UCA、AGU、ACU、ACA、GCU、GCA、AAU、UAU、UAA、UGG、CAA、GAA、AAA、GAU、CAU、UGU、AGA、CGU、CGA、GGU、GGA、CCU、CCC和CCA)，即筛选出了34个高频密码子，其中以A、U、C、G结尾的密码子分别有13、16、2和1个，这说明密码子偏好以A和U结尾，RSCU较高的3个密码子分别为UUU (1.94)、CUA (1.73)和UCU (1.75)。
2.3 中性绘图分析
中性绘图分析量化自然选择和突变压力之间的关系，阐明3个密码子位置之间的联系。结果表明：横坐标GC3的数值为23.14%~41.23%，纵坐标GC12的数值为39.04%~61.00% (图1)。梁山慈竹的Pearson相关系数为0.17，呈正相关关系，数据拟合后的回归系数为0.1868，决定系数(R2)较小，为0.0282，GC12和GC3的相关性不显著，说明其叶绿体基因组密码子偏好性受自然选择影响较大。
2.4 ENC-plot绘图分析
图2显示：ENC分布并不紧密，少量分布在标准曲线附近，还有个别分布在标准曲线上侧，位点的ENC均大于35，与预期ENC值有差距。说明梁山慈竹密码子偏好性较弱且自然选择和突变都对其偏好性有影响。由于落在标准曲线下方的基因点数量比较多，所以梁山慈竹基因组密码子使用偏好性主要受自然选择的影响。
2.5 PR2-plot绘图分析
图3显示：基因位点在平面图4个区域内分布并不均匀，在A3/(A3+T3)＜0.5和G3/(G3+C3)＞0.5区域范围内分布最多。表明第3位碱基使用频率为：T＞A、G＞C，梁山慈竹叶绿体基因组密码子的第3位碱基在选择上具有偏好性，同时说明其密码子使用偏好性主要受自然选择的影响。
2.6 最优密码子的确定
对梁山慈竹的ENC进行升序排列，前10%为高表达基因，即rps18、rpl16、psbA、rps14、rps11，后10%为低表达基因，即 ycf3、cemA、ndhE、atpE、infA。梁山慈竹的RSCU和ΔRSCU表明(表5)：梁山慈竹叶绿体基因组有32个高频密码子，筛选出GCA、GCU等25个高表达密码子，最终确定18个密码子作为梁山慈竹叶绿体基因组的最优密码子，分别为UAA、GCA、GCU、UUC、GGU、AAA、CUU、UUA、CCA、CCU、CAA、AGA、CGU、AGU、UCC、ACU、GUA、GUU。其中16个以A/U结尾，2个以C结尾。
表 5 梁山慈竹叶绿体基因组各氨基酸的RSCU分析及最优密码子分析Table 5 RSCU analysis and optimal codon analysis of amino acids in chloroplast genome of D. farinosus氨基酸 密码子 基因组
RSCU高表达
RSCU低表达
RSCUΔRSCU 氨基酸 密码子 基因组
RSCU高表达
RSCU低表达
RSCUΔRSCU Ter UAA*** 1.560 0 1.800 0 1.200 0 0.600 0 Met AUG 1.000 0 1.000 0 1.000 0 0 UAG 0.780 0 0.600 0 1.200 0 −0.600 0 Asn AAC* 0.520 0 0.893 6 0.625 0 0.268 6 UGA 0.670 0 0.600 0 0.600 0 0 AAU 1.480 0 1.106 4 1.375 0 −0.268 6 Ala GCA** 1.200 0 1.200 0 0.734 7 0.465 3 Pro CCA** 1.080 0 0.800 0 0.500 0 0.300 0 GCC 0.600 0 0.457 1 0.653 1 −0.196 0 CCC 1.010 0 0.800 0 1.166 7 −0.366 7 GCG 0.470 0 0.228 6 0.734 7 −0.506 1 CCG 0.430 0 0.444 4 1.000 0 −0.555 6 GCU* 1.730 0 2.114 3 1.877 6 0.236 7 CCU*** 1.480 0 1.955 6 1.333 3 0.622 3 Cys UGC** 0.470 0 0.400 0 0 0.400 0 Gln CAA* 1.530 0 1.500 0 1.368 4 0.131 6 UGU 1.530 0 1.600 0 2.000 0 −0.400 0 CAG 0.470 0 0.500 0 0.631 6 −0.131 6 Asp GAC* 0.460 0 0.500 0 0.411 8 0.088 2 Arg AGA* 1.750 0 1.723 4 1.534 9 0.188 5 GAU 1.540 0 1.500 0 1.588 2 -0.088 2 AGG 0.640 0 0.319 1 0.837 2 −0.518 1 Glu GAA 1.460 0 1.189 2 1.578 9 −0.389 7 CGA 1.270 0 1.276 6 1.395 3 −0.118 7 GAG** 0.540 0 0.810 8 0.421 1 0.389 7 CGC 0.510 0 0.319 1 0.837 2 −0.518 1 Phe UUC** 1.290 0 1.041 7 0.650 0 0.391 7 CGG 0.410 0 0.319 1 0.279 1 0.040 0 UUU 0.710 0 0.958 3 1.350 0 −0.391 7 CGU*** 1.410 0 2.042 6 1.116 3 0.926 3 Gly GGA 1.580 0 1.253 7 1.818 2 −0.564 5 Ser AGC 0.410 0 0.384 6 0.470 6 −0.086 0 GGC 0.430 0 0.417 9 0.484 8 −0.066 9 AGU** 1.250 0 1.846 2 1.411 8 0.434 4 GGG 0.760 0 0.119 4 0.363 6 −0.244 2 UCA 1.020 0 0.615 4 1.058 8 −0.443 4 GGU*** 1.240 0 2.209 0 1.333 3 0.875 7 UCC*** 1.190 0 1.769 2 0.941 2 0.828 0 His CAC** 0.530 0 0.941 2 0.571 4 0.369 8 UCG 0.550 0 0.153 8 0.705 9 −0.552 1 CAU 1.470 0 1.058 8 1.428 6 −0.369 8 UCU 1.580 0 1.230 8 1.411 8 −0.181 0 Ile AUA 0.920 0 0.850 7 0.949 4 −0.098 7 Thr ACA 1.080 0 1.181 8 1.176 5 0.005 3 AUC* 0.590 0 0.626 9 0.531 6 0.095 3 ACC 0.500 0 0.818 2 1.058 8 −0.240 6 AUU 1.480 0 1.522 4 1.519 0 0.003 4 ACG 0.480 0 0.363 6 0.588 2 −0.224 6 Lys AAA** 1.440 0 1.471 7 1.155 6 0.316 1 ACU** 1.680 0 1.636 4 1.176 5 0.459 9 AAG 0.560 0 0.528 3 0.844 4 −0.316 1 Val GUA*** 1.500 0 1.767 4 1.257 1 0.510 3 Leu CUA 0.900 0 0.833 3 1.295 5 −0.462 2 GUC 0.490 0 0 1.028 6 −1.028 6 CUC 0.420 0 0 0.545 5 −0.545 5 GUG 0.540 0 0.372 1 0.342 9 0.029 2 CUG 0.330 0 0.250 0 0.477 3 −0.227 3 GUU** 1.470 0 1.860 5 1.371 4 0.489 1 CUU* 1.290 0 1.333 3 1.227 3 0.106 0 Trp UGG 1.000 0 1.000 0 1.000 0 0 UUA*** 1.940 0 2.166 7 1.022 7 1.144 0 Tyr UAC** 0.410 0 0.521 7 0.166 7 0.355 0 UUG 1.110 0 1.416 7 1.431 8 −0.015 1 UAU 1.590 0 1.478 3 1.833 3 −0.355 0 说明: 高频密码子(RSCU＞1.00)带下划线；*. ΔRSCU≥0.08；**. ΔRSCU≥0.3；***. ΔRSCU≥0.5; 加粗的密码子表示最优密码子。 2.7 梁山慈竹和其他几种生物密码子偏好性比较分析
将梁山慈竹基因组密码子使用频率与巨龙竹、粉麻竹、小叶龙竹、硬头黄竹、大肠埃希菌、烟草、拟南芥和酿酒酵母等物种的基因组密码子使用频率进行比较(图4)。结果显示：梁山慈竹与巨龙竹、粉麻竹、小叶龙竹和硬头黄竹的密码子使用频率为0.5~2.0，说明它们的密码子使用偏好性相似，推测具有亲缘关系的禾本科Gramineae牡竹属植物叶绿体基因组密码子偏好性相似；在大肠埃希菌、烟草、拟南芥和酿酒酵母的密码子使用比值中筛选≥2.0或≤0.5的密码子，分别有28和15、15、14个，表明梁山慈竹与这些物种在同义密码子的偏好性上有一定差异。
图 4 梁山慈竹与其他物种密码子偏好性比较Figure 4 Comparison of codon preference between D. farinosus and other species2.8 对应分析
将梁山慈竹的51个叶绿体基因的基因功能分为光合系统基因、遗传系统基因、其他基因和未知功能基因四大类，在计算RSCU的基础上将各个基因分布到59维的向量空间。对应分析结果(图5)显示：前4个向量轴分别存在18.3%、16.8%、15.6%和15.4%的差异，前4向量轴累计差异为66.1%，4个轴对密码子均有不同程度的影响；第1轴的值大于其他轴，说明第1轴对梁山慈竹叶绿体基因组密码子偏好性的影响较大。对第1轴与CAI、CBI、FOP、ENC和GC3s等指数进行进一步的相关分析发现：梁山慈竹基因在第1轴上的坐标值与CAI (r=−0.001 7，P＜0.01)、CBI (r=0.099 0，P＜0.01)、FOP (r=0.083 0，P＜0.01)、ENC (r=0.112 0，P＜0.01)、GC3s (r=−0.145 0，P＜0.01)间具有极显著的相关关系，其中CAI和GC3s第1轴具有负相关关系，表明基因组密码子的偏好性不止受单一因素的影响，自然选择、基因突变均有可能影响梁山慈竹基因组密码子使用偏好性[18]。
3. 讨论
本研究对梁山慈竹叶绿体基因组密码子进行使用偏好性分析，筛选出51条CDS序列，分析表明：GC1＞GC2＞GC3，密码子在3个位置上的分布并不均匀，密码子偏好使用以A或U结尾的碱基，且梁山慈竹叶绿体基因组的ENC均值为49.51，表明其叶绿体基因组密码子使用偏好性较弱。这与乳油木Vitellaria paradoxa[19]和二乔玉兰Magnolia soulangeana[20]等植物叶绿体基因组密码子偏好性相似。
对梁山慈竹叶绿体基因组密码子进行中性绘图、ENC-plot分析、PR2-plot分析和对应分析。在中性绘图分析中，回归系数为0.412 8，说明密码子偏好性更多受到自然选择的影响；在ENC-plot分析中，多数基因离标准曲线距离较远，实际ENC和预期ENC有差距，表明该部分基因的密码子偏好性主要受自然选择的影响；在PR2-plot绘图分析中，大部分基因位于平面图的右下方，即T＞A、G＞C，表明其密码子的使用更多受自然选择的影响。综上所述，影响梁山慈竹叶绿体基因组密码子偏好性的主要原因是自然选择。该研究结果与巨桉Eucalyptus grandi[21]、灰毛浆果楝Cipadessa cinerascens、酸枣Ziziphus jujuba var. spinosa[22]和云南油杉Keteleeria evelyniana[23]等叶绿体基因组密码子偏好性研究结果基本一致；但在对4种蔷薇科 Rosaceae果树[24]和银白杨Populus alba[25]的研究中发现：突变是影响密码子偏好性的主要因素。这说明密码子的使用偏好性受自然选择或基因突变因素影响。基于RSCU的对应分析表明：梁山慈竹的密码子使用变异原因除了突变和自然选择之外，还有其他的因素，这其中光合系统基因和遗传系统基因分布相对集中，各类基因密码子使用偏好性较为接近。该结论与木薯Manihot esculenta[26]和高山松Pinus densata[27]的研究结果一致。密码子使用频率比较结果显示：梁山慈竹与禾本科牡竹属的植物密码子偏好性相似，在基因选择外源系统表达时，可以选择密码子偏好性差异相对较小的酿酒酵母，在选择大肠埃希菌、烟草和拟南芥作为外源表达宿主时，需要根据密码子使用偏好性进行碱基优化，从而使基因在宿主体内更好地表达。
最优密码子分析表明：梁山慈竹叶绿体基因组有GCU、GAU以及GGU等18个最优密码子，最优密码子大部分以A或U结尾。该结果与抽筒竹Gelidocalamus tessellatus[28]和毛竹Phyllostachys edulis[29]叶绿体基因组最优密码子分析结果一致，这可能与亲缘关系相近，但不同物种之间叶绿体基因组进化过程中的相对保守性有关系[21]。通过筛选获取梁山慈竹偏好使用密码子，可进一步对目标基因进行密码子优化，提高梁山慈竹的竹笋产量和造纸纤维含量，以及利用新一代精准基因编辑工具CRISPR/Cas9优化梁山慈竹密码子，从而改造梁山慈竹基因组编辑的Cas9基因，提高该基因在梁山慈竹中的表达水平[30]。
4. 结论
本研究通过分析梁山慈竹叶绿体基因组的CDS序列，对梁山慈竹的叶绿体基因组进行生物信息学分析，筛选出梁山慈竹叶绿体基因组有GCU、GAU以及GGU等18个最优密码子。研究结果表明：影响梁山慈竹密码子偏好性的主要因素是自然选择。研究结果为后续在分子层面上利用基因工程开发梁山慈竹优良资源提供参考。
-
图 1 研究区和样地分布示意图
Figure 1 Geographical location and plot distribution of the study area
图 2 3种LAI时序曲线在阔叶林和针叶林区域的变化情况以及重建效果
Figure 2 Changes and reconstruction effects of 3 LAI time series curves in broad-leaved forest and coniferous forest areas
图 3 同化LAI与实测值的散点验证图
Figure 3 Scatter verification diagram of assimilated LAI and measured value
图 4 融合LAI与实测值的散点验证图
Figure 4 Scatter verification diagram of fused LAI and measured values
表 1 LAI产品数据及基本信息
Table 1. LAI product data and basic information
产品名称 传感器 空间分辨率/m 时间分辨率/d 投影方式 数据来源 MOD15A2H MODIS 500 8 正弦曲线 https://ladsweb.modaps.eosdis.nasa.gov/ VNP15A2H VIIRS 500 8 正弦曲线 https://ladsweb.modaps.eosdis.nasa.gov/ CGLS 300m V1 PROBA-V 300 10 简易圆柱 https://land.copernicus.eu/global/ 
下载: 导出CSV
表 2 3种LAI产品在阔叶林和针叶林数据同化前后的精度指标
Table 2. Accuracy indexes of 3 LAI products before and after data assimilation in broad-leaved and coniferous plots
LAI产品 样地 同化前 同化后 R ERMSE EMAE ERE/% R ERMSE EMAE ERE/% MODIS 阔叶林 0.02 1.60 1.37 −23.75 0.59 0.37 0.60 −2.56 针叶林 0.07 1.46 1.30 −0.72 0.59 0.24 0.64 8.34 VIIRS 阔叶林 0.04 2.00 1.71 −28.62 0.56 0.31 0.60 2.29 针叶林 0.17 0.95 0.77 11.55 0.49 0.28 0.72 10.53 PROBA-V 阔叶林 0.05 1.20 1.04 15.24 0.62 0.14 0.86 12.49 针叶林 0.16 1.54 1.39 28.39 0.56 0.19 1.03 20.53
下载: 导出CSV
-
[1] CHEN J M, BLACK T A. Defining leaf area index for non-flat leaves [J]. Plant,Cell &Environment, 1992, 15(4): 421 − 429. [2] MYNENI R B, HOFFMAN S, KNYAZIKHIN Y, et al. Global products of vegetation leaf area and fraction absorbed PAR from year one of MODIS data [J]. Remote Sensing of Environment, 2002, 83(1/2): 214 − 231. [3] WU Wei, SUN Yin, XIAO Kun, et al. Development of a global annual land surface phenology dataset for 1982−2018 from the AVHRR data by implementing multiple phenology retrieving methods [J/OL]. International Journal of Applied Earth Observation and Geoinformation, 2021, 103: 102487[2023-12-10]. doi: 10.1016/j.jag.2021.102487. [4] ZOU Dongxiao, YAN Kai, PU Jiapin, et al. Revisit the performance of MODIS and VIIRS leaf area index products from the perspective of time-series stability [J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2022, 15: 8958 − 8973. [5] VERGER A, SÁNCHEZ-ZAPERO J, WEISS M, et al. GEOV2: improved smoothed and gap filled time series of LAI, FAPAR and FCover 1 km copernicus global land products [J/OL]. International Journal of Applied Earth Observation and Geoinformation, 2023, 123: 103479[2023-12-10]. doi: 10.1016/j.jag.2023.103479. [6] FUSTER B, SÁNCHEZ-ZAPERO J, CAMACHO F, et al. Quality assessment of PROBA-V LAI, fAPAR and fCOVER collection 300 m products of copernicus global land service [J/OL]. Remote Sensing, 2020, 12(6): 1017[2023-12-10]. doi: 10.3390/rs12061017. [7] YAN Kai, PU Jiapin, PARK T, et al. Performance stability of the MODIS and VIIRS LAI algorithms inferred from analysis of long time series of products [J/OL]. Remote Sensing of Environment, 2021, 260: 112438[2023-12-10]. doi: 10.1016/j.rse.2021.112438. [8] MA Han, LIANG Shunlin. Development of the GLASS 250-m leaf area index product (version 6) from MODIS data using the bidirectional LSTM deep learning model [J/OL]. Remote Sensing of Environment, 2022, 273: 112985[2023-12-10]. doi: 10.1016/j.rse.2022.112985. [9] XIAO Zhiqiang, SONG Jinling, YANG Hua, et al. A 250 m resolution global leaf area index product derived from MODIS surface reflectance data [J]. International Journal of Remote Sensing, 2022, 43(4): 1409 − 1429. [10] XING Luqi, LI Xuejian, DU Huaqiang, et al. Assimilating multiresolution leaf area index of Moso bamboo forest from MODIS time series data based on a hierarchical Bayesian network algorithm [J/OL]. Remote Sensing, 2019, 11(1): 56[2023-12-10]. doi: 10.3390/rs11010056. [11] 李雪建, 毛方杰, 杜华强, 等. 双集合卡尔曼滤波LAI同化结合BEPS模型的竹林生态系统碳通量模拟[J]. 应用生态学报, 2016, 27(12): 3797 − 3806. LI Xuejian, MAO Fangjie, DU Huaqiang, et al. Simulating of carbon fluxes in bamboo forest ecosystem using BEPS model based on the LAI assimilated with dual ensemble Kalman filter [J]. Chinese Journal of Applied Ecology, 2016, 27(12): 3797 − 3806. [12] LIU Qiang, LIANG Shunlin, XIAO Zhiqiang, et al. Retrieval of leaf area index using temporal, spectral, and angular information from multiple satellite data [J]. Remote Sensing of Environment, 2014, 145: 25 − 37. [13] 李喜佳, 肖志强, 王锦地, 等. 双集合卡尔曼滤波估算时间序列LAI[J]. 遥感学报, 2014, 18(1): 27 − 44. LI Xijia, XIAO Zhiqiang, WANG Jindi, et al. Dual ensemble Kalman filter assimilation method for estimating time series LAI [J]. Journal of Remote Sensing, 2014, 18(1): 27 − 44. [14] 周红敏, 张国东, 王昶景, 等. 塞罕坝地区高空间分辨率叶面积指数时序估算与变化检测[J]. 遥感学报, 2021, 25(4): 1000 − 1012. ZHOU Hongmin, ZHANG Guodong, WANG Changjing, et al. Time series high-resolution leaf area index estimation and change monitoring in the Saihanba area [J]. Journal of Remote Sensing, 2021, 25(4): 1000 − 1012. [15] 张慧芳, 高炜, 施润和. 基于背景库的高质量LAI时间序列数据重建[J]. 遥感学报, 2012, 16(5): 986 − 999. ZHANG Huifang, GAO Wei, SHI Runhe. Reconstruction of high-quality LAI time-series product based on long-term historical database [J]. Journal of Remote Sensing, 2012, 16(5): 986 − 999. [16] XU Bao, PARK T, YAN Kai, et al. Analysis of global LAI/FPAR products from VIIRS and MODIS sensors for spatio-temporal consistency and uncertainty from 2012−2016 [J/OL]. Forests, 2018, 9(2): 73[2023-12-10]. doi: 10.3390/f9020073. [17] 魏萌, 范文义, 张海军, 等. 4-Scale几何光学模型冠层反射率模拟的空间尺度适用性[J]. 应用生态学报, 2023, 34(3): 605 − 613. WEI Meng, FAN Wenyi, ZHANG Haijun, et al. Spatial scale applicability of canopy reflectance simulation of 4-Scale geometric optical model [J]. Chinese Journal of Applied Ecology, 2023, 34(3): 605 − 613. [18] ROGERS C, CHEN Jingming, CROFT H, et al. Daily leaf area index from photosynthetically active radiation for long term records of canopy structure and leaf phenology [J/OL]. Agricultural and Forest Meteorology, 2021, 304: 108407[2023-12-10]. doi: 10.1016/j.agrformet.2021.108407. [19] WEN Yibo, ZHUANG Linlan, WANG Hezhi, et al. An Automated Hemispherical Scanner for monitoring the leaf area index of forest canopies [J/OL]. Forests, 2022, 13(9): 1355[2023-12-10]. doi: 10.3390/f13091355. [20] 贾玉洁, 刘云根, 杨思林, 等. 面向Sentinel-2A影像的大理市土地利用分类方法适用性研究[J]. 浙江农林大学学报, 2022, 39(6): 1350 − 1358. JIA Yujie, LIU Yungen, YANG Silin, et al. Applicability of land use classification method in Dali City based on Sentinel-2A image [J]. Journal of Zhejiang A&F University, 2022, 39(6): 1350 − 1358. [21] CHEN Jingming, DENG Feng, CHEN Mingzhen. Locally adjusted cubic-spline capping for reconstructing seasonal trajectories of a satellite-derived surface parameter [J]. IEEE Transactions Geoscience and Remote Sensing, 2006, 44(8): 2230 − 2238. [22] ZHAI Huan, HUANG Fang, QI Hang. Generating high resolution LAI based on a modified FSDAF model [J/OL]. Remote Sensing, 2020, 12(1): 150[2023-12-10]. doi: 10.3390/rs12010150. [23] 杨绍钦, 王翔, 许澄, 等. 基于MODIS时间序列数据的竹林地上生物量估算[J]. 浙江农林大学学报, 2022, 39(4): 734 − 741. YANG Shaoqin, WANG Xiang, XU Cheng, et al. Estimation of aboveground biomass of bamboo forest based on MODIS time series data [J]. Journal of Zhejiang A&F University, 2022, 39(4): 734 − 741. [24] JI Jiayi, LI Xuejian, DU Huaqiang, et al. Multiscale leaf area index assimilation for Moso bamboo forest based on Sentinel-2 and MODIS data [J/OL]. International Journal of Applied Earth Observation and Geoinformation, 2021, 104: 102519[2023-12-10]. doi: 10.1016/j.jag.2021.102519. [25] 于颖, 刘敏, 范文义, 等. 基于PROSPECT和4-scale模型的光化学植被指数尺度转换[J]. 北京林业大学学报, 2020, 42(10): 27 − 35. YU Ying, LIU Min, FAN Wenyi, et al. Scale conversion of photochemical reflectance index based on PROSPECT and 4-scale models [J]. Journal of Beijing Forestry University, 2020, 42(10): 27 − 35. [26] 王晓冰, 范文义. 叶面积指数与植被指数饱和点分析及验证[J]. 东北林业大学学报, 2023, 51(9): 83 − 94, 111. WANG Xiaobing, FAN Wenyi. Analysis and verification of saturation point of leaf area index and vegetation index [J]. Journal of Northeast Forestry University, 2023, 51(9): 83 − 94, 111. [27] 周杨凡, 郑兴明, 孙源, 等. GF-1与Sentinel-2在典型农作物LAI估算中的比较[J]. 遥感技术与应用, 2023, 38(3): 599 − 613. ZHOU Yangfan, ZHENG Xingming, SUN Yuan, et al. Comparison of GF-1 and Sentinel-2 for estimation of leaf area index in typical crops [J]. Remote Sensing Technology and Application, 2023, 38(3): 599 − 613. [28] 孟庆博, 刘艳丽, 鞠琴, 等. 基于多源数据同化融合的尼洋河降水时空分布特征[J]. 南水北调与水利科技, 2020, 18(3): 110 − 118. MENG Qingbo, LIU Yanli, JU Qin, et al. Analysis of temporal and spatial distribution characteristics of precipitation based on multi-source data assimilation and fusion in Niyang River Basin [J]. South-to-North Water Transfers and Water Science &Technology, 2020, 18(3): 110 − 118. [29] 周旻悦, 沈润平, 陈俊, 等. 基于像元质量分析和异常值检测的LAI时序数据S-G滤波重建研究[J]. 遥感技术与应用, 2019, 34(2): 323 − 330. ZHOU Minyue, SHEN Runping, CHEN Jun, et al. Filtering and reconstruction of LAI time series data by S-G filter based on pixel quality analysis and outlier detection [J]. Remote Sensing Technology and Application, 2019, 34(2): 323 − 330. [30] BROWN L A, MEIER C, MORRIS H, et al. Evaluation of global leaf area index and fraction of absorbed photosynthetically active radiation products over North America using Copernicus Ground Based Observations for Validation data [J/OL]. Remote Sensing of Environment, 2020, 247: 111935[2023-12-10]. doi: 10.1016/j.rse.2020.111935. 期刊类型引用(1)
1. 何刀山，阳小强，秦雅林，何海燕，谢维，罗治国，李鹏. ‘湘辣14号’叶绿体基因组密码子偏好性分析. 中国果菜. 2025(01): 47-54+79 . 
百度学术其他类型引用(0)
-
-
链接本文:
https://zlxb.zafu.edu.cn/article/doi/10.11833/j.issn.2095-0756.20230601
点击查看大图
计量
- 文章访问数: 307
- HTML全文浏览量: 76
- PDF下载量: 21
- 被引次数: 1
