2. 广东省林业科学研究院,广东 广州 510520
2. Guangdong Academy of Forestry, Guangzhou 510520, China
林木遗传改良通常需在田间试验基础上,估算参试样本(种源、家系、单株或无性系)目标性状的遗传参数,并运用一定的评价方法筛选出具有优异目标性状的种源、家系或无性系[1]。目前,林木遗传评估基本采用混合线性模型,即根据观测表型值和设计的关联矩阵,利用最佳线性无偏估计(Best linear unbiased estimation, BLUE)方法获得固定效应估计值,利用最佳线性无偏预测(Best linear unbiased prediction, BLUP)方法获得随机效应预测值。而育种值(Breeding value)即属随机效应预测值,且BLUP值依赖于混合线性模型的关联设计矩阵。然而,对林木田间试验而言,个体间存在竞争,且基因型与环境互作较为显著,即空间效应明显。如竞争效应和空间效应没有被考虑到混合线性模型中,遗传评估则会产生偏差,其中空间正相关易导致加性遗传方差偏大,中等竞争效应则抑制偏差[2]。因而,在林木遗传评估模型中,应纳入竞争效应和空间效应。
空间效应(或称环境异质性)可以是空间连续的,其表现为土壤和小气候效应的相似模式,也可以是不连续的,表现为育林措施和测量方法的不同效应,此外还可以是随机的,其表现为微环境的异质性[3]。针对空间效应的空间分析(Spatial analysis)已被广泛用于农作物品种田间试验研究。与传统试验设计相比,空间分析可减少试验误差并提高农作物品种评估的准确性。空间分析方法通常有趋势面分析法和空间自相关分析法等。其中,趋势面分析法是采用空间坐标的多项式函数来模拟环境的变异[4],而空间自相关分析法则通过行和列的空间自相关分析剔除地点内的空间相关误差来提高遗传试验分析的精确性[5]。
在林木研究上,Costa等[6]利用空间自相关分析法分析了辐射松Pinus radiata等3个树种的子代和无性系试验林,发现空间自相关分析可更有效拟合试验林的空间变异,显著降低了区组和小区的方差,且极大提高了遗传评估的准确度与遗传增益。另外,Dutkowski等[7-8]在更大数据量的树种与试验林上进一步采用空间分析,也发现空间分析可显著提高种源、家系、亲本和无性系遗传效应估算的准确性。但林木是多年生植物,除空间差异外,往往还存在个体竞争[9]。虽然Cappa等[10]采用遗传干预模型(Genotypic interference model, GIM)评估了林木遗传加性效应和竞争加性效应,但是没有考虑空间效应,因此模型比较受限。值得注意的是,Stringer等[11]提出空间效应与竞争效应的联合模型进行糖蔗品种试验的遗传分析,模型采用等根二阶自回归法(Equal-roots second-order autoregressive, EAR2)拟合环境的空间效应,而遗传的竞争效应则采用处理干预模型(Treatment interference model, TIM)。此外,Cappa等[12]提出拟合空间效应与竞争效应的林木单株混合模型,并采用贝叶斯法估计遗传参数。上述模型均在一定程度上提高了林木遗传评估的准确性,但模型涵盖的遗传竞争效应均为单一值,而实际上林木个体周边的遗传竞争效应可能不同,因此有必要提出新的单株混合模型,同时拟合不同的遗传竞争效应和异质的空间效应。基于此,本研究根据现代线性统计理论,构建了异质空间效应和不同竞争效应的联合模型,并通过数据模拟的方式,结合实际测定数据进行验证,同时运用REML 法估算遗传参数,对不同模型进行比较,进而探讨异质空间效应和不同竞争效应的分析模型在林木遗传测定中的有效性,以期为林木精确的遗传评估和后续可靠的良种选育提供理论参考和技术支撑。
1 方法 1.1 统计模型单株混合模型为:y=Xb+Zgug+e。式中: y为个体性状表型值构成的向量; b为固定效应构成的向量; u为随机遗传效应构成的向量,
遗传效应的方差–协方差矩阵为:
${\rm{var}}\left( {{{u}_{\rm{g}}}} \right) = \left[ \begin{array}{l}\begin{array}{*{20}{c}}{\delta _{{{\rm{g}}_{\rm{d}}}}^2} & {{\delta _{{{\rm{g}}_{{\rm{den}}}}}}} & {{\delta _{{{\rm{g}}_{{\rm{dwn}}}}}}}\\{{\delta _{{{\rm{g}}_{{\rm{den}}}}}}} & {\delta _{{{\rm{g}}_{{\rm{en}}}}}^2} & {{\delta _{{{\rm{g}}_{{\rm{ewn}}}}}}}\\{{\delta _{{{\rm{g}}_{{\rm{dwn}}}}}}} & {{\delta _{{{\rm{g}}_{{\rm{ewn}}}}}}} & {\delta _{{{\rm{g}}_{{\rm{wn}}}}}^2}\end{array}\begin{array}{*{20}{c}}{{\delta _{{{\rm{g}}_{{\rm{dsn}}}}}}} & {{\delta _{{{\rm{g}}_{{\rm{dnn}}}}}}}\\{{\delta _{{{\rm{g}}_{{\rm{esn}}}}}}} & {{\delta _{{{\rm{g}}_{{\rm{enn}}}}}}}\\{{\delta _{{{\rm{g}}_{{\rm{wsn}}}}}}} & {{\delta _{{{\rm{g}}_{{\rm{wnn}}}}}}}\end{array}\\\begin{array}{*{20}{c}}{{\delta _{{{\rm{g}}_{{\rm{dsn}}}}}}} & {{\delta _{{{\rm{g}}_{{\rm{esn}}}}}}} & {{\delta _{{{\rm{g}}_{{\rm{wsn}}}}}}}\\{{\delta _{{{\rm{g}}_{{\rm{dnn}}}}}}} & {{\delta _{{{\rm{g}}_{{\rm{enn}}}}}}} & {{\delta _{{{\rm{g}}_{{\rm{wnn}}}}}}}\end{array}\begin{array}{*{20}{c}}{\delta _{{{\rm{g}}_{{\rm{sn}}}}}^2} & {{\delta _{{{\rm{g}}_{{\rm{snn}}}}}}}\\{{\delta _{{{\rm{g}}_{{\rm{snn}}}}}}} & {\delta _{{{\rm{g}}_{{\rm{nn}}}}}^2}\end{array}\end{array} \right] \otimes {{I}_m},$ |
其中,对角线上依次为直接遗传效应、东向近邻效应、西向近邻效应、南向近邻效应和北向近邻效应的方差分量,非对角线为上述两两效应之间的协方差分量。var代表方差,δ2 和 δ 分别代表方差和协方差,Im代表单位矩阵。
为减少计算量和促进模型收敛,可通过因子分析法[1]将遗传效应的方差–协方差矩阵改写为如下矩阵:
${\rm{var}}\left( {u} \right) = \left( { \varGamma \times \varGamma ' + \varPsi } \right) \otimes {A},$ |
其中,Γ 是因子载荷矩阵,Ψ 是特殊方差矩阵,A是加性遗传相关矩阵。
在常规的空间分析[7]中,剩余误差(R)可分解为空间自相关误差方差(
${R} = \delta _\xi ^2\left( {\sum\nolimits_c {\rm{(}}{\rho _{_c}}{\rm{)}} \otimes \sum\nolimits_r {\rm{(}}{\rho _{_r}}{\rm{)}}} \right) + \delta _\eta ^2{I},$ |
其中,
$\sum {(\rho )} = \left| {\begin{array}{*{20}{c}}1 & \rho & {{\rho ^2}} & \cdots & {{\rho ^{n - 1}}}\\\rho & 1 & \rho & \cdots & \vdots \\{{\rho ^2}} & \rho & 1 & \cdots & \vdots \\ \vdots & \vdots & \vdots & \ddots & \vdots \\{{\rho ^{n - 1}}} & \cdots & \cdots & \cdots & 1\end{array}} \right|\text{。}$ |
采用R软件[1]和R程序包breedR[13]进行数据模拟。为简化数据模拟,参考下述松树数据,以该松树胸径均值(13)和估算的加性方差(
数据来源于文献[1]中的例6.4.4.2,该数据有46个全同胞家系,采用完全随机区组设计,设7个区组,5株行式小区,试验林共35行40列,测定性状是10年生胸径。
1.4 模型分析方法所有模型分析均采用ASReml软件V4.0[14],除总体均值作为固定因子外,所有其他试验因子均作为随机因子。空间相关的残差采用行列AR1自相关结构,遗传竞争效应的方差采用单因子XFA1结构。拟合模型的原始数据图和残差图采用R包AAfun[15]生成,空间变异图variogram由ASReml运行生成,所有估计的参数由ASReml运行获得。
单株狭义遗传力(
$h_i^2 = \delta _a^2/\left( {\delta _a^2 + \delta _e^2} \right),$ |
式中,
AIC值的计算公式如下:
${\rm AIC} = - 2 \times ({\rm ln}L - P)\text{,}$ |
式中,AIC值是赤池信息量准则,lnL是模型迭代收敛似然值,P是模型中随机参数的数量。
2 结果与分析 2.1 不同模型比较对模拟数据设置随机区组设计模型(RCBM)、空间模型(SM)、空间与测量误差模型(SUM)和空间与竞争模型(SCM)4种模型,模型的拟合结果(表1)显示:模型RCBM的参数最少(4个),收敛似然值(lnL)最小(–1 927),AIC值最大(3 863),剩余残差(7.56)也最大,所估算的单株狭义遗传力也最小(0.24±0.08);模型SM的参数为6个,收敛似然值(lnL)稍微增加,为–1 904,AIC值下降至3 820,空间自相关误差为7.04,加性方差为3.44,比RCBM提高了45%,但估计的行列相关值远低于模拟设置值,行自相关值为0.25,列自相关值为–0.07;模型SUM的参数中,随机测量误差方差为5.72,空间自相关误差方差为1.83,行列相关值与模拟设置值一致,加性方差为2.43,估计的单株狭义遗传力为0.30±0.10;模型SCM 参数最多(10个),似然值(lnL)最大(–1 853),但AIC值(3 734)和随机测量误差方差(3.13)均最小,估计的单株狭义遗传力最高,为0.43±0.13。此外,模型SCM估计的4个近邻效应方差在0.34~0.74,并且均显著。根据似然值lnL最大和AIC值最小原则,可知SCM为最佳模型。
残差图显示,RCBM的残差图(图1C)与原始数据图(图1A)非常相似,表明RCBM模型未能很好地拟合空间异质性。虽然SM模型比RCBM有所改进,但也只拟合了残差小的空间部分(图1D)。SUM和SCM模型,图形比较相似,有明显的行列空间趋势,但SCM模型更能消除部分残差高的区域(图1E、1F)。从空间变异图可以看出,SM模型的图形比较平坦(图2B),表明捕获的空间变异很小,而SUM和SCM模型的图形比较相似(图2C、2D),在行方向上存在明显的空间变异趋势,在列方向的突起尖峰表示列相关为负值,存在环境的竞争效应,表明SUM和SCM模型捕获空间变异的效果不错。此外,残差值与近邻表型均值的散点图(图2A)显示两者存在明显的负相关(r= –0.30),表明存在遗传竞争效应。上述结果与表1的结果相一致。
为了解不同初始参数值条件下,SCM模型拟合的效果,根据竞争相关值、行列自相关值设计了9组数据,结果(表2)表明,所有数据拟合的加性方差为2.25~2.52,随机测量误差方差为2.90~3.96,空间自相关误差方差为1.31~2.43,除了列相关值(pc= –0.3)拟合得一般外,其他均与初始值接近。各近邻效应方差多数在0.45~0.80,基本上都显著;此外,估计的单株狭义遗传力差异不大,均在0.40左右。可见,SCM 模型对于空间效应和竞争效应的拟合效果稳定性较好。
2.3 真实数据的模型拟合结果比较真实数据拟合的结果(表3)显示,常规的RCBM模型,估计的区组方差为0.19,加性方差为2.46,剩余残差为2.96;仅有空间行列自相关的模型SM,加性方差被高估为7.58;加上区组因子后的空间模型SM1,这种偏差被修正,加性方差为2.74,空间自相关误差方差为2.84,但行、列自相关值仅为0.11、0.13;“空间+测量”误差模型SUM,有较好的拟合结果,加性方差为2.21,随机测量误差方差为2.80,空间自相关误差方差为0.80;空间与竞争模型SCM为最优模型,加性方差为2.25,随机测量误差方差和空间自相关误差方差均比SUM有所降低,分别为2.41和0.72,所估算的近邻竞争方差在0.01~0.13,但没达到显著水平。从似然值lnL和AIC值可知,SCM模型仍是最优。
3 讨论与结论本研究利用R软件及其程序包breedR模拟数据,结合实测数据,采用XFA1结构拟合加性效应、近邻竞争效应和AR1结构拟合空间效应,建立了RCBM、SM、SUM和SCM 4种模型并进行比较,发现SCM是最优模型,该模型可大大降低随机误差方差,并估算到四周近邻竞争方差,设置参数的不同初始值,SCM估计的参数结果均较为稳定。
对林木研究而言,试验林面积一般较大,采用常规试验设计难以避免环境异质性,传统的分析模型因忽略环境的异质性而导致遗传参数和育种值估算的偏差,因此Dutkowski 等[8]提出林木遗传试验数据应尽可能采用空间分析。对于环境的空间变异,一般使用3次B样条函数(Cubic B spline)或行列自回归拟合。考虑到残差结构的灵活性,本文采用行列自回归法进行空间分析。从模拟数据的似然值lnL和AIC值结果可以看出,空间分析模型比常规分析模型拟合的结果更好,但是没有随机误差方差的空间模型,随机误差被纳入空间相关误差,从而导致空间相关误差方差被高估,因此,也造成空间变异趋势被抑制,使得行、列自相关值被低估,这从残差图和空间变异图也可得到佐证。然而,仅有行、列自相关的空间模型会导致加性方差被严重高估[3],本文也得到一致的结果。Magnussen[2] 指出,环境异质性一般存在于试验林早期,林木个体周围相似的微环境会导致个体间存在正相关,但随着林木生长,个体间的竞争日趋严重,将会导致个体间存在负相关。因此,对于多年生的林木试验林,个体间存在竞争效应是必然的,这是造成个体致死率的重要因素[9]。换言之,林木性状的表型中,会同时存在空间效应和竞争效应,忽略任何一个因素,都可能会造成遗传评估的偏差,并进一步影响后续的选择工作。对于遗传竞争效应,拟合的方法有GIM[10]、TIM[11]和CSM[12],但上述研究中遗传竞争效应均为单一值,忽略了林木个体周边的遗传竞争效应可能不同。本文提出的以加性效应及其四周近邻竞争效应的方差–协方差矩阵作为G结构,替代常规的加性G结构,从而达到拟合近邻竞争效应值不唯一的目的。考虑到降低运算量和模型收敛,采用因子分析法来拟合G结构。模拟数据的拟合结果表明,新的SCM模型,拟合方差虽然没有增加,但随机误差方差大大降低了,分别由7.56(RCBM)、5.72(SUM)降低到3.13(SCM),分别降低了58.6%、45.3%,这与同时考虑竞争效应和空间效应的研究模型结果[12]相一致。此外,拟合数据的SCM模型单株狭义遗传力保持在0.40左右,高于RCBM(0.24)和SUM模型(0.30)。而针对参数不同初始值的模拟数据,SCM模型的拟合结果都比较稳定,表明其可用于林木的遗传分析。最后,利用松树的实测数据验证SCM模型,估算结果与模拟结果比较一致,笔者认为,新单株混合模型即SCM模型可用于林木遗传分析,在提高遗传评估的准确性和可靠性方面具有优势。
[1] |
林元震, 陈晓阳. R与ASReml-R统计分析教程[M]. 北京: 中国林业出版社, 2014.
(0) |
[2] |
MAGNUSSEN S. A method to adjust simultaneously for spatial microsite and competition effects[J]. Can J For Res, 1994, 24(5): 985-995. DOI:10.1139/x94-129 (0) |
[3] |
DUTKOWSKI G. Improved models for the prediction of breeding values in trees [D]. Australia Tasmania: University of Tasmania, 2005.
(0) |
[4] |
YANG R C, YE T Z, BLADE S F, et al. Efficiency of spatial analyses of field pea variety trials[J]. Crop Sci, 2004, 44(1): 49-55. DOI:10.2135/cropsci2004.4900 (0) |
[5] |
GILMOUR A R, CULLIS B R, VERBYLA A P. Accounting for natural and extraneous variation in the analysis of field experiments[J]. J Agric Biol Environ Stat, 1997, 2(3): 269-293. DOI:10.2307/1400446 (0) |
[6] |
COSTA E SILVA J, DUTKOWSKI G W, GILMOUR A R. Analysis of early tree height in forest genetic trials is enhanced by including a spatially correlated residual[J]. Can J For Res, 2001, 31(11): 1887-1893. DOI:10.1139/x01-123 (0) |
[7] |
DUTKOWSKI G W, COSTA E S J, GILMOUR A R, et al. Spatial analysis methods for forest genetic trials[J]. Can J For Res, 2002, 32(12): 2201-2214. DOI:10.1139/x02-111 (0) |
[8] |
DUTKOWSKI G W, COSTA E S J, GILMOUR A R, et al. Spatial analysis enhances modeling of a wide variety of traits in forest genetic trials[J]. Can J For Res, 2006, 36(7): 1851-1870. DOI:10.1139/x06-059 (0) |
[9] |
LIN Y Z, YANG H X, IVKOVIĆ M, et al. Effect of genotype by spacing interaction on radiata pine genetic parameters for height and diameter growth[J]. For Ecol Manage, 2013, 304: 204-211. DOI:10.1016/j.foreco.2013.05.015 (0) |
[10] |
CAPPA E P, CANTET R J C. Direct and competition additive effects in tree breeding: Bayesian estimation from an individual tree mixed model[J]. Silvae Genet, 2008, 57(2): 45-56. (0) |
[11] |
STRINGER J K, CULLIS B R, THOMPSON R. Joint modeling of spatial variability and within-row interplot competition to increase the efficiency of plant improvement[J]. J Agric Biol Environ Stat, 2011, 16(2): 269-281. DOI:10.1007/s13253-010-0051-5 (0) |
[12] |
CAPPA E P, MUÑOZ F, SANCHEZ L, et al. A novel individual-tree mixed model to account for competition and environmental heterogeneity: A Bayesian approach [J]. Tree Genet Genomes, 2015, 11: 120. doi: 10.1007/s11295-015-0917-3.
(0) |
[13] |
MUÑOZ F, SANCHEZ L. R package version 0.7-16.2014: BreedR: Statistical methods for forest genetic resources analysts [CP/OL]. [2016-12-12]. https:// github.com/famuvie/breedR.
(0) |
[14] |
GILMOUR A R, GOGEL B J, CULLIS B R, et al. ASReml user guide release 4.0[M]. UKHemel Hempstead: Vsn International Ltd Hemel, 2016.
(0) |
[15] |
LIN Y Z. R package version 2.6.1.2016: AAfun: ASReml-R Added Functions [CP/OL]. [2016-12-12]. https://github.com/yzhlinscau/AAfun.
(0) |