TGG 一篇关于G by E的综述

Li Y, Suontama M, Burdon R D, et al. 2017. Genotype by environment interactions in forest tree breeding: review of methodology and perspectives on research and application. Tree Genetics & Genomes, in press.

摘要

环境相互作用的基因型(G×E)是指不同环境之间不同基因型的比较性,代表基因型排名差异或环境遗传差异表达水平的差异。 G×E可以减少遗传力和整体遗传增益,除非育种计划的结构可以解决不同类别的环境。了解G×E的影响,环境在产生G×E中的作用以及问题和机遇对于有效的育种程序设计和部署遗传物质至关重要。我们回顾了目前用于识别G×E:因子分析模型,双标图分析和反应规范【反应规范:基因型对环境反应的幅度,即在一定的环境条件下,特定的基因型所产生的表型的变动范围。】的主要分析方法。我们还回顾了G×E对全球经济重要性的森林物种的生长,形态和木材性质的生物学和统计学证据,包括一些松树,桉树,道格拉斯冷杉,云杉和一些杨树。在这些物种中,高水平的G×E倾向于报告生长性状,对于形态特征和木材性质,G×E含量低。最后,我们讨论了利用G×E可以最大限度地发现森林树种育种中的遗传增益。表征环境在产生互动中的作用被视为基本平台,允许有效测试候选基因型。我们讨论相对于等级交互的表达式互动的重要性,比以前的许多报告更为重要,尤其是部署决策。我们研究了G×E对树木育种的影响,一些导致G×E的环境因素和G×E处理树种育种的策略,以及基因组学的未来作用。
关键词:环境相互作用的基因型。分析方法G×E司机。影响。探索策略。森林树种育种。基因型筛选。基因型部署

介绍

个体的表型受其基因型,环境和基因型与环境之间的任何相互作用(G×E相互作用或G×E)控制。据说当基因型的比较性能根据环境而变化时存在这种相互作用。在一种环境中优越的一种基因型的表现在另一种环境中可能较差(de Jong 1990; Falconer和Mackay 1996)。在林业中,G×E可导致某些基因型在某些环境中的表现不可预测。因此,表征和理解G×E以减少甚至消除这种不可预测性是一个总体目标。如果发生这种情况,G×E会使育种计划的设计变得复杂,但它也提供了将最合适的种植种群与目标部署条件相匹配的机会,以优化森林健康,生长和集约管理森林的木材质量。
G×E可以分为两大类(Lynch和Walsh 1998; White等,2007):( 1)等级变化相互作用,其中基因型在不同环境中按不同顺序排列; (2)表达水平相互作用[1],其中基因型差异的表达(例如育种值的扩散)在不同环境中变化,不一定是基因型排序顺序的任何变化。因此,育种者需要确定G×E的模式和大小,以便为森林工业获得最佳的遗传增益(Muir等1992; Raymond和Namkoong 1990)。由于育种者主要关注候选基因型的评估和选择,因此排名变化互动通常对他们更感兴趣。虽然表达水平相互作用对于育种通常不太重要,但是对于帮助确定在特定生长环境和特定终产物中选择用于部署的遗传物质可能是主要的兴趣。
在林木育种中,G×E通常可能是实质性的,并且在寻找始终如一的优良基因型方面产生问题,特别是对于广泛适应。育种者通常通过选择对环境变化不敏感的稳定基因型,或选择特定环境的基因型以最大化该位点的遗传增益来解决G×E(Raymond和Namkoong 1990)。在树木中,几乎所有商业上重要的物种,如辐射松(Pinus radiata D. Don; Raymond 2011; Wu和Matheson 2005),火炬松(Pinus taeda L .; McKeand等1997;保罗)都报道了显着的G×E等1997年),樟子松(Pinus sylvestris L .; Gullberg和Vegerfors 1987; Haapanen 1996),湿地松(Pinus elliottii Englem.var.eoliottii; Hodge and White 1992; Roth等2007),桉树(Eucalyptus spp.;例如Costa e Silva等2006; Hardner等2010),花旗松(Pseudotsuga menziesii(Mirb)Franco; Campbell 1992; Dungey等2012),云杉(Picea spp .; Bentzer等1988)和杨树(Populus spp .; Rae等人2008; Yu和Pulkkinen 2003)) 。
提高对G×E的规模和性质的认识有助于树木种植者和林农看到增加森林遗传增益的机会(Cullis等2014; Ivkovi等2013a,b)。随着这种意识的出现,从对G×E的思考转变为主要与育种群体的选择和结构相关,以及在部署方面更多地考虑G×E。对于育种而言,获取额外遗传增益的预期收益可能不会超过多环境育种计划的复杂性和额外成本(Carson 1991)。然而,对于部署,新西兰树木种植企业对某些地点实现遗传增益的担忧导致Radiata Pine Breeding Company(RPBC)重新评估了这种情况(Butcher,2015,个人通信)。多环境因素分析模型的应用也给出了新的统计效率水平和证据,通过在正确的位置部署正确的基因型,可以获得额外的遗传增益(Cullis等,2014)。
在这篇综述中,我们研究了用于鉴定G×E及其性质的可用分析方法,对于选定的商业重要物种,森林树木中G×E的重要性和模式的经验证据,驱动G×E的环境变量的信息在林业和树木育种计划中处理G×E的策略。基于此信息,讨论了探索G×E以优化林木育种中遗传增益的意义和挑战。还讨论了基因组选择在森林树育种中鉴定G×E的重要性。我们对审查统计方法的重点是分析现有试验数据以供选择和其他育种决策。确定G×E的环境驱动因素以供将来筛选环境选择的方面仅受到有限的关注,被视为单独综述的主题。

估算G×E相互作用的分析方法

已经提出了许多分析方法来测量树木育种中性状的G×E程度,包括相互作用与遗传方差的比率,稳定性分析,B型遗传相关双标图分析,因子分析模型和反应规范。所有这些方法都可以作为混合线性模型的特殊情况来实现。可以使用方差分析来估计相互作用与遗传方差的比率(例如Shelbourne 1972)。稳定性分析已用于鉴定多种环境中的稳定或敏感基因型(例如Eberhart和Russell 1966; Finlay和Wilkinson 1963; Huehn 1990)。 B型遗传相关性(Burdon,1977)和因子分析模型表征了多种环境中基因型排序变化的模式(Cullis等,2014)。似然比检验现在提供了与完美B型相关性偏离的相互作用的稳健测试。 Biplot分析结合方差分析(ANOVA)和主成分分析(PCA)并允许结果可视化(Gauch 1992; Mandel 1971; Yan等2007)。反应规范描述了在一系列环境中由单一基因型产生的一系列反应或表型(Lynch和Walsh 1998; Pierce 2005; Woltereck 1909)。前三种方法已经在作物和林木育种中的G×E书籍或评论文章中得到了很好的阐述(Freeman 1973; Kang and Gauch 1996; Zobel et al 1988)。本综述将涵盖双标图分析,因子分析模型和反应规范。

用于估计G×E相互作用的分析方法

Biplot分析

Biplot分析最常用于植物育种,目的是开发高产量和高质量的作物。由于G×E,从一种环境中选择的品种可能无法在另一种环境中保持其高性能。在多种环境中鉴定具有高性能和广泛适应性的栽培品种是最终目标双标图分析使用奇异值分解将数据分解为分量矩阵,同时显示列和行信息(Gabriel 1971)。主效可加互作可乘模型(AMMI)以及基因型主效应和G×E效应(GGE)是鉴定植物和林木育种中G×E模式的两种主要的双标图分析方法。 AMMI和GGE双标图分析测试G×E的显着性和G×E方差与遗传方差的相对大小,并可视化基因型和环境中基因型最佳表达的稳定性。
AMMI是结合ANOVA和PCA的统计方法(Gauch 1992; Mandel 1971; Mrode 2014)。它使用ANOVA分解变异来源,首先分解为基因型和环境的加性效应,然后使用PCA分解G×E的乘法效应(Zobel等,1988)。相互作用效应被分解为表示对G×E的实际响应的部分和由于随机变化引起的部分(Crossa等1990)。用于AMMI分析的线性模型如下(Gauch 1992):
y_{i j}=\mu+\alpha_{i}+\beta_{j}+\sum_{k=1}^{n} \lambda_{k} \xi_{i k} \eta_{j k}+\rho_{i j}+\varepsilon_{i j}

其中y_{i j}是环境j中基因型i的观察表型,\mu是平均值,\alpha_{i}是基因型主效应与\mu的偏差,\beta_{j}是环境主效应与\mu的偏差,\lambda_{k}是轴k相互作用主成分的奇异值(IPC) ),\xi_{i k}\eta_{j k}是轴k的基因型和环境IPC得分(即左和右奇异向量),\rho_{i j}是包含模型中未包括的所有乘法项的相互作用残差; n是模型保留的轴或主成分的数量,\varepsilon_{i j}是与环境j处的基因型i相关联的残差,假设独立且具有相同的分布。 AMMI模型(\mu\alpha_{i}\beta_{j})的可加部分是从ANOVA和来自PCA的乘法部分(\lambda_{k}\xi_{i k}\eta_{j k})估计的。
GGE方法也基于ANOVA和PCA,使用具有两个主要组成部分的站点回归模型(SREG)。 GGE分析中使用的线性模型如下(Yan and Hunt 2001):
y_{i j}-\beta_{j}=\sum_{k=1}^{2} \lambda_{k} \xi_{i k} \eta_{j k}+\varepsilon_{i j}

其中y_{i j}是环境j中基因型i的平均表现,\beta_{j}是环境j中种植的所有基因型的平均表现,\lambda_{k}是主成分k的奇异值,\xi_{i k}\eta_{j k}是基因型i和环境j的主成分k的得分。\varepsilon_{i j}是与基因型i和环境j相关的残差。
GGE方法通过ANOVA消除环境主效应,并在环境中心数据中保留基因型主效应(G)和相互作用效应(G×E)(Yan等2000)。它允许通过PCA直接显示跨多个环境的基因型的性能和稳定性。当环境是与基因型和G×E对总变异性的贡献相关的主要变异来源时,推荐该模型(Kandus等,2010)。在GGE双标图中,连接双标图原点和环境标记的线是环境向量。两个环境的矢量之间的角度与它们之间的相关系数有关。该角度的余弦近似于相关系数(Yan 2002),这相当于Ding等人(2008b)进行的研究中的B型遗传相关性。 GGE双标图的局限性在于,当基因型主效应远小于相互作用效应或G×E模式复杂时,它可以仅解释总GGE的一小部分(Ding等2008b)。 GGE双标图方法不适合严格的假设检验,因此它更适合作为假设 - 产生者而不是决策者(Ding et al 2008b)。
AMMI和GGE方法都使用线性模型,并将主要和交互效应视为固定效应(Crossa 2012)。残差方差具有正态分布,并且在整个环境中是均匀且独立的(Gauch 1992; Kang and Gauch 1996; Piepho 1995)。双标图程序用于通过绘制每种基因型和环境的相互作用效应的PCA评分来提供结果的图形解释(Crossa 1990; Crossa等1990; Kempton 1984; Yan等2000)。它们之间的区别在于GGE双标图分析基于以环境为中心的PCA,而AMMI分析是指以双重为中心的PCA(Kroonenberg 1995; Rad等人2013)。 GGE双标图有许多可视化解释,AMMI在展示“which-won-where”时则没有,特别是当可视化任何交叉G×E时(Ding et al 2008b)。
AMMI分析方法已广泛应用于作物育种中的多环境栽培试验的统计分析(Annicchiarico 1997; Crossa等1990,1999;表1; Gauch和Zobel 1989; Gauch和Zobel 1997; Hassanpanah 2010)。 AMMI也被用于表征森林树木中的G×E,例如在Eucalyptus(Baril等1997a; Karuntimi 2012; Lavoranti等2007),松树(Pinus spp .; Chambel等2008; Falkenhagen 1996; Kim等2008),杨树(Populus spp .; Rae et al 2008),白色云杉(Picea glauca(Moench)Voss; Rweyongeza 2011),桦树(Betula spp.; Zhao et al 2014)和lodgepole pine(Pinus contorta Douglas; Wu and Ying 2001)。
GGE方法已被用于研究许多农艺试验中的G×E(Yan等2000,2007)以及树育种(Correia等人2010; Ding等人2008b; Sixto等人2011)。 Correia等(2010)利用GGE双标图分析研究了来自葡萄牙,西班牙,法国和澳大利亚的30种海洋松(Pinus pinaster Aiton)种群的总高度,直径,茎形和存活率的G×E,多种环境来源在西班牙的审判。 GGE双标图分析用于分析杨树无性系的生物量产量和稳定性,根据欧洲西南部的平均表现和稳定性进行排序(Sixto等2011)。 Ding等(2008b)使用GGE双标图研究了澳大利亚五个环境中216个辐射松科的G×E。双标图分析允许根据它们的性能和展示的相互作用程度来鉴定不同的克隆组,从而为辐射松的选择过程提供有用的信息。

因子分析模型

因子分析(FA)模型可以为估计所有试验对之间的遗传相关性提供可靠,简约和整体的方法(Cullis等人2014; Smith等人2015),并为复杂多重模型中的G×E模式提供自然框架。 - 环境实验(Meyer 2009)。 FA模型对于选择育种群体和生产种群的部署决策最有用。
在多环境试验中使用FA模型是基于使用PCA的特征向量(Jolliffe 1986; Smith等2001)并扩展以适应加性和非加性效应(Oakey等2006a,b)。 Cullis等(2014)使用FA模型来适应大量环境以及动物模型减少的环境之间的连通性差。 FA模型旨在确定导致变量之间相关性的统计公共因子(Mrode 2014)。它们将在多种环境下评估的特征表示为几个潜在变量的线性组合(Cullis等人2014; Smith等人2001),称为共同因子(Hardner等人2010; Meyer 2009),从而减少了其中的维度 - 地点相对于G×E的变化。因子的数量称为模型的阶数,而阶数k的FA模型表示为FAk。假设线性混合模型:
\mathbf{y}=\boldsymbol{X} \boldsymbol{\beta}+\boldsymbol{Z} \boldsymbol{u}+\boldsymbol{e}
其中y是地点t的观测值向量; \beta是固定效应的向量; u\boldsymbol{u} \sim N(0, \boldsymbol{G} \otimes \boldsymbol{A})的随机加性遗传效应向量,其中G=[ \begin{array}{cccccc}{\sigma_{a_{1}}^{2}} & {\cdots} & {\sigma_{a_{1} a_{t}}} & {\vdots} & {\ddots} & {\sigma_{a_{t} a_{1}} \cdots \sigma_{a_{t}}^{2}}\end{array}] ,其中\sigma_{a_{i}}^{2}是地点i的加性遗传方差,\sigma_{a_{i} a_{j}}是站点i和站点j之间的加性遗传协方差,A是遗传关系矩阵,\otimes表示Kronecker积; ee \sim N(0,[\sigma_{e_{1}}^{2} \cdots 0 \vdots \ddots 0 \cdots \sigma_{e_{t}}^{2}])随机残差效应的向量,\sigma_{e_{i}}^{2}是地点i的残差方差,XZ是将表型与\betau关联的设计矩阵,t是试验次数。 在t试验中m个基因型的加性遗传效应的FAK模型可以建模为\boldsymbol{u}=\left(\boldsymbol{\Lambda} \otimes \boldsymbol{I}_{m}\right) \boldsymbol{f}+\boldsymbol{\delta}(Costa e Silva等2006; Cullis等2014),其中\Lambdat \times k的试验载荷矩阵, fmk \times 1的得分向量,\deltamt \times 1的遗传回归残差向量。 \operatorname{Var}(\boldsymbol{u})=\left(\Lambda \Lambda^{\prime}+\Psi\right) \otimes I_{m},假设\operatorname{var}(\boldsymbol{f})=\boldsymbol{I}_{m k}\operatorname{var}(\boldsymbol{\delta})=\boldsymbol{\psi} \otimes \boldsymbol{I}_{m},其中\psit \times t的对角矩阵,每个环境有一个方差(称为特定方差),随机效应向量 f\delta是相互独立的多元高斯分布,均值为零。 环境间遗传方差矩阵定义为G_{e}=\left(\Lambda \Lambda^{\prime}+\psi\right)。 可以使用REML算法估计Ge为G_{e}=\left(\Lambda \Lambda^{\prime}+\psi\right)可以转换为相关矩阵C_{e}=D_{e} G_{e} D_{e}D_{e}是对角矩阵,其元素是G_{e}的对角元素的平方根的倒数。 上面概述的FA模型等同于Meyer(2009)指定的扩展因子分析模型。
Latent 回归图用于显示对试验载荷的遗传响应,表明FA模型中多个环境中选择候选物的G×E(或稳定性)的大小(Chen等2017; Cullis等2014; Smith等2015;表格1)。具有较高斜率的选择候选者的潜在回归意味着候选者对环境更敏感具有零斜率的潜在回归表明候选者的表现在多个环境中是稳定的。该方法已用于研究大麦多种环境(Smith等2001),马铃薯(Burgue等2011),玉米(Burgue等2011)和小麦(Burgue等2011; Oakey等)的G×E。 al 2006a,b),以及畜牧业(Meyer 2009)。在林木育种中,这种方法已成功应用于Khasi松(Gus E)(Pinus kesiya Royle ex Gordon; Costa e Silva 2007),辐射松(Cullis等2014; Ivkovi等2015),火炬松( Zapata-Valenzuela 2012),Eucalyptus和Eucalyptus杂种无性系(Costa e Silva等2006; Hardner等2010)和挪威云杉(Picea abies(L.)H. Karst.; Chen等2017)的研究。尽管具有统计效力,FA方法尚未确定特定环境因素在推动G×E中的作用(B. Cullis,2016,个人通信)。

反应规范

反应规范,也称为反应规范,描述了在一系列环境中由单一基因型产生的一系列反应或表型(Lynch和Walsh 1998; Pierce 2005; Woltereck 1909)。 它适用于分析在环境梯度(例如温度)上逐渐且连续变化的性状数据(Kolmodin和Bijma 2004)。 单个性状的线性反应范数具有模型(Kolmodin和Bijma 2004; Strandberg等2000):
y_{i k}=b_{0}+b_{1} x_{k}+a_{0_{i}}+a_{1_{i}} x_{k}+e_{0_{i}}+e_{1_{i}} x_{k}

其中y_{i k}是环境k中基因型i的表型值; b_{0}b_{1}分别是基因型反应范数的截距和斜率的固定效应; a_{0_{i}}a_{1_{i}}分别是基因型i的反应范数的截距(对应于经典EBV表现潜力)和斜率(相当于环境敏感性的EBV)的随机加性遗传效应; e_{0_{i}}e_{1_{i}}是基因型i的反应范数的截距和斜率的随机残差效应; 和x_{k}是环境k对表型的影响。 假设随机加性遗传效应a_0a_1服从期望为零,方差分别为\sigma_{a_{0}}^{2}\sigma_{a_{1}}^{2}、协方差\boldsymbol{\sigma}_{\boldsymbol{a}_{0} \boldsymbol{a}_{1}}的正态分布。
反应规范的优点是选择响应不仅可以在任何环境中的表型表达中预测,而且可以通过线性反应范数的斜率(稳健性或对环境变化的响应性)来量化性状的环境敏感性(Kolmodin和Bijma 2004)。在进行反应规范分析时,疾病暴露、植株密度和营养质量被认为是影响牲畜生产的环境因素(Rauw和Gomez-Raya 2015)。同样,与林木育种相关的环境因素,如温度,降雨量和土壤营养,可以应用于反应规范分析,以确定林木育种中G×E的驱动因素反应规范的局限性在于分析中包含的任何环境变量需要清楚地识别,并且所识别的G×E模式仅在模拟环境条件的范围内有效(Gregorius和Kleinschmit 2001)。通过将观察到的表型值与一个或多个环境变量作图,可以显示反应范数模型中给定基因型的稳定性。反应范数的平坦斜率意味着在整个环境范围内产生的表型是恒定的。基因型的反应规范与整个群体的反应规范之间的任何差异构成了G×E的组成部分。例如,对于给定的基因型,特别陡峭的反应范数斜率意味着其表型对环境更敏感,并且基因型将对G×E有强烈贡献。 G×E的水平也可以表示为加性遗传斜率(\sigma_{a_{1}}^{2})与加性遗传截距(\sigma_{a_{0}}^{2})(Kolmodin和Bijma 2004)的方差比或加性遗传斜率与截距之间的遗传相关性(r_g),r_{\mathrm{g}}=\frac{\sigma_{a_{0}} a_{1}}{\sqrt{\sigma_{a_{0}}^{2} \sigma_{a_{1}}^{2}}}(Strandberg等2000)。
反应规范方法已应用于森林树木(Gregorius和Kleinschmit 2001),包括广泛分布在林地范围内的物种,特别是挪威云杉(Oleksyn等1998)和樟子松(Abraitiene等2002),以及小杆松(Rehfeldt等1999; Wang等2006)及其在加拿大的杂交种(Wu和Ying 2001)。 Oleksyn等(1998)研究了54种挪威云杉种群的植物生长,分配,净CO2交换率,组织化学和物候,以量化海拔梯度中种群间生长和相关植物性状的差异,并更好地了解它们之间的关系。这项研究表明,来自寒冷山区环境的挪威云杉种群可以通过几种潜在的适应性特征来表征,例如年平均温度和海拔高度。 Abraitiene等(2002)研究了樟子松花粉生活力和臭氧易感性的遗传变异。发现花粉对臭氧浓度增加的显着遗传变异,但只有5%的变异归因于G×E。
反演规范概念用于推导将气候变量纳入分析模型的响应函数,并大大提高了它们在小杆松种群中的可靠性(Wang等2006)和两种桉树之间特定结合能力的结构(Baril等1997b) )。响应函数预测,气候的微小变化极大地影响了森林种群的生长和存活,并且在全球变暖期间保持当代森林生产力需要在整个景观中批量重新分配基因型(Rehfeldt等1999)。通过在不列颠哥伦比亚省内部的57个种源试验地点测量的20年高度,发现从3个小杆松亚种(Pinus contorta ssp.controrta,ssp.latifolia和ssp.urrayana)取样的10个自然小杆松种群中的大量种群间的相互作用(Wu和Ying 2001)。

林木中G×E的证据

森林树木中G×E的证据在森林树木中,G×E已经在各种经济上重要的物种中进行了研究,特别是在新西兰,澳大利亚,美国,欧洲,亚洲和非洲(表1)。大多数研究调查了G×E的生长和形态特征,并对木材密度特性和木材性状特征进行了一些G×E调查。表1列出了文献中各种林木研究中G×E的证据。两项研究提出了测量林木育种中G×E大小的标准。 Robertson(1959)建议将0.8或更高的遗传相关性解释为G×E,具有较低的生物学重要性。 Shelbourne(1972)提出,当相互作用方差达到遗传方差的50%或更多时,相互作用对选择和测试的遗传增益有严重影响。
通常发现高G×E用于树木生长(例如辐射松; Carson 1991; Johnson和Burdon 1990;表1; Wu和Matheson 2005)。在澳大利亚10个地点的双列实验中,估计的相互作用与胸径(DBH)的遗传方差之比高于0.50(Wu和Matheson 2005)。 DBH环境对之间的遗传相关性估计为0.39(Wu和Matheson 2005),0.34-0.38(Ding等人2008a)和-0.60至1.0(Raymond 2011)。 B型遗传相关系数为0.27-0.84(Dieters和Huber 2007; Li和Mckeand 1989; Roth等2007; Sierra-Lucero等2003),高度为0.18-0.95(Gwaze等2001; Li和Mckeand 1989; Owino等1977; Paul等1997)和0.27表示每公顷体积的平均年增量(Sierra-Lucero等2003)。 Baltunis和Brawner(2010)报道了澳大利亚遗址中克隆试验的生长性状高G×E,但新西兰遗址没有。在使用覆盖整个新西兰76个站点的数据的分析中,配对站点之间DBH的近三分之二的遗传相关性估计值低于0.6(McDonald和Apiolaza 2009)。 Ivkovi等(2015)报道,在20个对照授粉试验中,DBH的成对遗传相关性范围为-0.51至0.98,其中48%与1的完美遗传相关性显着不同,平均值为0.35。对试验之间B型遗传相关性的估计随着年龄的增长而增加,表明G×E的重要性似乎随着年龄的增长而下降,早期生长数据对于评估成熟时的G×E可能是不可靠的(Dieters等1995; Gwaze等2001); Roth等人2007; Zas等人2003)。
树形通常是树木育种计划中的重要特征。干直度和分枝等特征赋予轮伐年龄树木的价值(Cown等1984; Ivkovi等2006)。这些性状中G×E的程度在研究中变化很大。这与在分支大小,分叉数和分枝数(Wu和Matheson 2005)的形式特征以及澳大利亚和新西兰的一些地点(Baltunis和Brawner 2010)中的干直度和分枝质量的G×E的一些证据形成对比。同样,Gwaze等(2001)报道津巴布韦茎直线性的G×E含量高,Suontama等(2015)报道了新西兰分枝的高水平G×E。在新西兰的两项研究中发现了G×E的冲突:Dungey等人(2012)报道干直度的高水平的种源 - 地点相互作用,但花旗松中的家族内部物种相互作用水平较低; Kennedy等(2011)发现干直度和形式评分的高B型遗传相关性(超过0.80)和畸形的低B型遗传相关性(0.49)。
与木结构和质量相关的性状通常在针叶树中具有低G×E,例如,木材基本密度(Apiolaza 2012; Baltunis等人2010; Gapare等人2010,2012a; Johnson和Gartner 2006; Muneri和Raymond 2000),声速或弹性模量(Dungey等人2012; Gapare等人2012a; Jayawickrama等人2011; Johnson和Gartner 2006),木材化学性质(Sykes等2006),木材比重(Jett等1991)和树脂管性状(Westbrook等2014)。 Osorio等(2001)报道了哥伦比亚巨桉(Eucalyptus grandis)木材密度的最小G×E。然而,发现巴西的桉树杂种无性系(Eucalyptus grandisxEucalyptus urophylla)在四个地点具有显着的G×E木材基础密度(Lima等2000)。
统计学显着的B型遗传相关性可能不反映G×E的真实幅度。当估计遗传参数时,良好的实验设计在测试的所有阶段至关重要。幼苗种群缺乏随机性显然导致遗传变异分配问题,造成家庭间差异被夸大(Baltunis等,2007)。环境之间遗传相关性的差(即非常不精确)估计通常与环境之间共同的有限数量的父母有关(Apiolaza 2012; Raymond 2011)。传播效果可能与插条生根的季节有关。在冬季环境与任何其他试验的根系插条建立的试验之间观察到最差的遗传相关性,而最佳的遗传相关性来自田间试验,其中包括源自春季环境的根系插条(Baltunis等2005, 2007)

G×E对林木育种的影响

G×E对树木育种者的总体影响是使育种计划设计复杂化。 在另一个环境中追求遗传增益的一种环境中选择性状的效率与两种环境之间的遗传相关性以及两种环境中性状的总体遗传力成正比(Falconer和Mackay 1996)。 整个环境的整体遗传力计算公式为h^{2}=\frac{\sigma_{a}^{2}}{\sigma_{a}^{2}+\sigma_{g e}^{2}+\sigma_{e}^{2}},其中\sigma_{a}^{2}是加性遗传方差,\sigma_{g e}^{2}是相互作用方差,\sigma_{e}^{2}是残差方差。因此,高G×E可以在两个方面降低站点的整体遗传力。首先,分母中的大G×E方差直接降低了总遗传力。其次,G×E也损害了多个环境下遗传变异的估计,进一步减小了遗传力估计的大小。例如,在七个地点的分析中,平均年增量的遗传度为0.08,而在两个地区的每个地区,其年均增量为0.08(Sierra-Lucero等,2003)。如果在一个环境中估计,G×E也可以增加遗传力的估计。当G×E存在时,单个位置检验的估计值用于一般遗传预测时,遗传度估计将被扩大,因为G×E方差的一部分被划分为加性遗传方差。例如,发现添加剂G×E足够大以致可靠性估计和遗传增益预测高达60-100%的上升偏差(Owino等,1977)。与单站点估计相比,家庭平均遗传度估计高估了约15%,与现场分析相比,家庭平均遗传率估计高达15%,甚至网站之间的B型遗传相关性也是0.85。
G×E可以影响测试和选择情景中预测遗传增益的估计,因为它降低了整个环境的整体遗传力或准确性。 Sierra-Lucero等人(2003年)发现,选择1区的家庭在2区进行部署,导致每公顷的平均年增量减少4%。在忽略G×E的情况下,与考虑G×E的情景相比,遗传增益的损失达到10%以上,当地点之间的B型遗传相关性小于0.80时(Diaz Solar等2011)。 Leksono(2009)报道,直接选择产生的遗传增益明显高于间接选择产生的遗传增益,如果繁殖种群在繁殖区之间转移,遗传增益降低24%0%。澳大利亚昆士兰州最北部地区的直径与其他地点的直径差异很小(r_g = 0.39),如果个别选择基于此地点在任何其他地点种植,估计的遗传增益仅为29%~7%的效率与基于种植场地的选择程序一样(Woolaston等1991)。
Xie(2003)发现内部云杉,白云杉(Picea glauca(Moench)Voss)和Engelmann云杉(Picea engelmannii Parry ex Engelm)复合体(Xie和Yanchuk 2002)在五个种子规划区中的高度相当大G×E位于不列颠哥伦比亚省中北部的内陆地区,平均遗传相关系数为0.64。五个种子规划区聚为两个新种子区,新种子区内的遗传相关性分别为0.97和0.84,两个新区之间的遗传相关性为0.41。在每个区域内选择最佳25%的受试父母时,将整个区域视为一个区域时预期的遗传增益为19%,考虑五个原始区域时为24%,将五个原始区域合并为两个区域时为26%区域。尽管区域化育种的预期收益明显有所增强,但未能在性能与其他地方相关性较差的地点评估基因型,这将不可避免地牺牲这些地点潜在遗传增益的一半或更多。
显着的G×E并不总是导致遗传增益的显着损失。 Jett等(1991)发现了火炬松中特定木材重力的显着G×E。 18个家系中有4个被归类为性状不稳定,占观察到的G×E方差的一半。然而,这种G×E仅对潜在的遗传增益产生微不足道的影响。 Dieters等(1996)报道,对于梭形锈病抗性,B型基因遗传相关性超过0.67,而G×E似乎对美国湿地松的抗锈性没有重要意义。 Carson(1991)发现辐射松直径有显着的G×E,但几种区域化选择预测的遗传增益表明新西兰辐射松的G×E大小似乎太小而不能保证区域化繁殖种群。
潜在增益损失的度量(C)已被用作评估G×E对家庭选择育种计划(Matheson和Raymond 1984,1986),C=\left(1-\left[\left({V}_{{f}}+\frac{{V}_{{e}}}{N B S}\right) /\left({V}_{{f}}+\frac{{V}_{{i}}}{S}+\frac{{V}_{{e}}}{{N B S}}\right)\right]^{\frac{1}{2}}\right) \times 100
其中V_f是表型方差,V_i是相互作用的方差,V_e是误差均方; S是站点数,B是每个站点的重复数,N是每个站点的树数。这假设选择强度保持相同,并且残留和遗传组分在模型中保持相同,包括或不包括相互作用项。使用前面的关系,2%的潜在增益损失相当于遗传数值的近似减少5%。

树木育种中G×E的驱动因素

确定哪些环境因素是G×E的关键驱动因素对于育种和部署目的都很重要。它通知选择要测试和评估候选基因型的环境。同样重要的是要知道特定环境中基因型的表现可以告诉我们他们在其他环境中的预期表现,这些环境可能具有其属性,但不一定是经验G×E数据。已经进行了各种研究来表征环境在辐射松中产生G×E的作用。当环境因子对于至少一些基因型处于次优或非次优水平时,经常发现G×E反映基因型之间的差异应激反应(Kang 2002)。胁迫可以是生物(疾病或害虫)或非生物(例如温度,盐度和水或营养素的过量或不足)。对于给定的一组基因型,环境越多样化,可能的G×E的幅度越大(Li等人2015)。 G×E可能是由于子区域或个体基因型对夏季促进水分和光胁迫的环境条件的不同适应性,在某种程度上可能与观察到的生物因子的差异敏感性有关(Costa e Silva等2006)。
Johnson和Burdon(1990)在新西兰的两个地点类别之间获得了极好的区分,代表了Northland粘土(天然磷缺乏)和浮石遗址,这种模式与澳大利亚(Fielding和Brown 1961)和New新西兰(Burdon 1971; Burdon 1976)。 Wu和Matheson(2005)发现,先前的土地使用根据互动行为创建了一些地点。从那时起,Raymond(2011)对新南威尔士州辐射松直径增长的G×E研究发现,海拔高度是G×E的主要驱动因素,对先前的土地利用和地质母质的影响较小。然而,海拔高度与一系列重要的生物气候因素有关,包括降雨量,季节性和温度变量。高纬度种源和冬季凉爽,夏季干燥,高海拔种源以及降水量高,生长季节短的地方,对高度的G×E和白云杉的DBH贡献最大(Rweyongeza 2011)。
Wu和Matheson(2005)报道,辐射松的区域间相互作用的大基因型归因于两个较高海拔地区的大量积雪。在整个新西兰(McDonald和Apiolaza 2009)的76个辐射松试验和来源和后代水平的最低温度(Gapare等2015年)中收集的数据分析中,发现DBH的G×E可能受到极端最高温度的驱动。 5月和6月的平均日温低于3.2度解释了27.8%的G×E相互作用,与FA模型中的第一个因子中度相关,表明春季或秋季的霜冻天气条件可能是挪威云杉GxE的主要驱动因素(Chen et al 2017)。 Ivkovi等(2013a)报道,根据Cullis等人(2014)估计的DBH育种值,高降雨量和低温解释了G×E方差的25%,它们可能是新西兰G×E的驱动因子。
另一方面,与降雨量及其季节性密切相关的叶子病等生物因素,对于茎直径和体积增长来说是G×E的明显潜在驱动因素(参见Ades和Garnier-G閞?1997) )。由树木遗传学和基因组中差异暴露于瑞士针刺引起的G×E(2017)13:60 Page 11 of 18 60道格拉斯冷杉也是一个很好的例子,因为针头病在温暖的温度下更为普遍,而且可能有因此,寄主树种群中疾病负荷和差异生长反应完全不同(Dungey等,2012)。
Li等(2015)报道,生长性状的G×E水平与土壤养分中氮和总磷和年平均温度的位置差异显着相关,叶面钙含量和分子量的G×E水平显着相关土壤中镁和钾含量的地点差异。在最广泛的范围内(新西兰和澳大利亚),温度和降雨等气候变量是驱动观测到G×E的最重要因素;然而,在当地的区域范围内,土壤和地形因素更为重要(Ivkovi等,2013a)。

在树种育种计划中处理G×E的策略

已经提出了两个主要策略来处理强G×E的存在(Kang 2002; Raymond 2011):(1)选择跨站点稳定的人员;或(2)选择非常适合每个个体环境的个体以最大化遗传增益。当没有找到观察到的G×E的明显来源并且相互作用被认为基本上是“噪声”时,第一种策略是适用的。该策略旨在选择具有广泛适应性的基因型,预期该基因型可在多种环境中可靠地产生。选择稳定性是澳大利亚辐射松的G×E研究中推荐的方法(Ding 2008; Matheson和Raymond 1984),新西兰(Carson 1991)和西班牙(Codesido和Fern醤dez-L髉ez 2009) ;美国的火炬松(Owino 1977; Paul等1997);瑞典的挪威云杉(Bentzer等,1988)。该方法旨在消除与环境最相互作用的基因型。然而,为了获得良好的成功,育种者需要使用一组适当的测试环境来暴露交互式和稳定的基因型。
第二种策略是通过分析和解释遗传和环境差异来开发相互作用(Raymond 2011)。预计该策略将在每个环境中单独地最大化遗传力和遗传增益,因此,这需要创建单独的育种群体和种子生产,从而导致成本,管理,记录和血统控制的问题(Barnes等1984)。当环境数量大,跨越不同的气候和地理区域以及不同的土壤类型时,这种策略可能是不切实际的。这种策略的实际改编是将类似的环境分组到区域。该策略的应用依赖于识别哪些站点或环境因素导致互动的能力(Raymond 2011)。区域化通常用于在其自然范围内生长的北半球物种,其中通常与温度梯度相关的单一环境因素决定种子表现。例如,对于瑞典的苏格兰松,基于纬度温度梯度建立了育种和种子区,这些生长性状区内的G×E非常低(Hannrup等2008; Raymond 2011)。

研究结果表明,辐射松的圆形交配设计实验似乎有利于将DBH辐射松林育种区分为澳大利亚的两个主要地区:新南威尔士州的高海拔Tumut地区以及维多利亚州,南澳大利亚州和西澳大利亚地区Wu and Matheson 2005)。然而,在新西兰,Carson(1991)报道,辐射松的G×E的大小似乎太小,无法保证区域性的繁殖种群。 Matheson和Raymond(1984)提出,更好地解决G×E问题不是试图区分育种,而是忽略似乎特别容易受环境变化影响的家庭。需要充分评估权衡,因为这些方法对运营育种有严重的影响,可能需要进行成本/效益分析。具有单独繁殖种群的区域性育种计划在每个环境中可能具有较小的育种种群,较低的选择强度,因此比同一大小的一个国家计划降低遗传增益(Carson 1991)。虽然通过区域化可以获得额外的收益,但土地利用的进一步成本,多个种子果园的运营,附加记录和子代测试中的延伸将会成比例地增加(Barnes et al。1984; Carson 1991)。

应用和未来研究

多特征背景

我们已经回顾了研究G×E的统计方法,主要是关于单一特征的情况。然而,定量遗传学家和育种者在选择和部署的多个特征的背景下几乎总是必须考虑G×E。在方法学中,B型遗传相关性的研究(Burdon 1977)可以容易地扩展到多种环境之间的多种特征。具有多重特征,可能会出现G×E的几个更复杂的表现。不同的特征可能表现出不同的G×E模式。此外,遗传和表型性状(类型A)相关矩阵可以在环境之间不同,作为另一种形式的G×E。在不同环境中表达的不同性状之间的相关性(其中提出了类型-AB相关性的术语)在G×E的复杂表现中可能在成对环境之间不同,可能与健康和霜冻或耐旱性相关的性状相关。

多特征案例在选择和部署之间创建秩和变化级别G×E之间的相互作用。如上所述,表达水平的相互作用可以产生等级变化相互作用,一些典型的例子是疾病表达的影响,其通过影响树木生长或树形形式肯定会影响后者的基因型排名,具有明显的影响基因型表现的稳定性。一个这样的例子涉及新西兰辐射松的原产地试验,其中对针刺疾病耐受性较差的Cambria来源对易发生发病的部位的生长相对性能要差得多(Burdon等,1997)。另一个涉及新西兰沿海道格拉斯冷杉的原产地试验,其中比较增长表现受到瑞士针刺(由Phaeocrypopus gaeumannii引起)的强烈影响,易受影响的最低纬度原产地在最高纬度地区做得比较好( 46?S),其中疾病风险远低于低纬度地区(38?S)(Dungey et al 2012)。这种与疾病有关的影响对于育种和部署都是有意义的。

表达水平的相互作用

表达水平G×E不会导致选择候选者的等级变化,并且通常对育种不太重要(Muir等人,1992),但是对于跨多个环境进行部署是重要的(Burdon等人,2017)。如果没有排名变化的互动,在一个环境中的测试应该是足够的。如果部署选择涉及一系列环境和多特征选择,则级别的表达式相互作用可能需要将不同的基因型集合部署到不同的环境中,而不必进行排名变化的交互。在这种情况下,可能需要对几种环境进行评估,以便对所有感兴趣的特征提供良好的遗传差异分辨率。水平表达相互作用的典型例子可能与抗病性相关,其中基因型差异的分辨率可以很大程度上取决于疾病发病率(例如Dieters等人,1996; Sohn和Goddard,1979)。通常,分辨率往往最适合中度到重度的疾病,特别是在最耐药的基因型之间。然而,可能存在一种疾病可能以不具有直接实际重要性的水平存在的情况,但在该环境中解决遗传差异可能仍然很好。对这种疾病的抵抗不需要在这样一个环境中作为部署标准,而在那里表达的抗病性仍然可能是在其他环境中部署的有价值的选择标准,其中疾病是麻烦的,但发病率不高,因为高水平的“噪音”变化。

基因组学将有助于解决G×E

正在越来越多地探索基因组选择作为森林树种育种选择的主要手段(Grattapaglia和Resende 2011; Isik et al。2011,2016; Lexer and St?ling 2012; Ratcliffe et al。2015; Resende et al。2011,2012a ,b)。基因组选择的主要优点是选择可以在正常的表型年龄之前进行 - 大约8岁的辐射松。 DNA可以提取,基因分型仅在几个针前进行,6个月之前。这意味着可以大大减少育种的产生间隔,并且可以增加每单位时间的预期遗传增益(Grattapaglia和Resende 2011; Isik 2014)。

已经对几种森林物种研究了基因组选择中的G×E。当G×E相互作用存在性状时,观察到的SNP对性状的影响在不同环境之间发生变化,而与性状的关系在一个环境中可能是显着的,而在其他环境中则不显着,如辐射松(Li et al。 )。与一个F2家族和三个F1家族的663个个体相比,发现了与DBH相关的QTL,基本木材密度,牛皮纸浆产量和纤维素,木质素木质素和提取物的木材化学成分以及木质素丁香醛对愈创木脂比的显着QTL在维多利亚州和西澳大利亚州塔斯马尼亚种植的桉树(Freeman等人,2011年,2013年)。然而,当使用一个群体的数据集开发模型来预测桉树中另一群体的表型时,基因组选择的准确性已显示显着下降(Resende等,2012a)。在美国的火龙果中,G×E严重影响了育种区域模型的转移能力(Resende et al。2012b)。在加拿大西海岸种植的室内云杉中,在预测不同地点的表型时,在多位点模型(其中拟合G×E项)的基因组选择准确度高于单站模型(El-Dien等人,2015)。

为了表征基因型和环境在驱动G×E在林木育种中的作用,我们可能需要在分子水平上收集基因型数据,评估更广泛的表型,并在所有季节中收集气候和地理数据以及土壤数据以及所测试的基因型增长的所有年份功能基因组学(Pevsner 2009),基因组序列注释(Ouzounis和Karp 2002; Wolf等人2001)和高分辨率表型(Crowell等,2016)可用于表征基因型和环境在驾驶G×E中的作用开发环境特异性基因组育种价值是在林木育种中使用基因组选择时,在多个环境中最大化遗传增益的下一个挑战。在上述分析方法中,因子分析模型(Cullis et al。,2014; Smith et al。,2015)是估计所有环境对之间的遗传相关性的简洁而整体的方法,最终可以发挥环境的重要作用特定的基因组育种值用于大量环境

概观

统计方法,AMMI,GGE二次分析,FA和反应规范在本文中进行了综述。表2总结了这些方法的优缺点。它们是确定森林树种繁殖中G×E的模式和幅度的工具。前三种方法推断环境和基因型分组,纯粹基于表型数据,并使用PCA来降低维度。 AMMI和GGE二分析分析测试G×E的显着性和G×E方差与遗传变异的相对大小,并且允许观察基因型最佳进行的稳定性和环境。这些方法最适合围绕跨多个环境部署基因型的决策。它们在森林数据的可视化中肯定会有用,以帮助您做出明智的部署决策。

FA模型能够在不使用过多数量的方差参数的情况下估计大量环境的非结构化遗传方差 - 协方差矩阵。这种类型的模型可以很容易地用于解释G×E的性质和程度。从FA模型估计的B型遗传相关性可以清楚地显示环境中是否存在基因型的等级变化。 FA模型在育种价值估计方面具有统计学的有效性,并提供可靠的整体方法来估计所有环境对之间的遗传相关性(Cullis等人,2014)。 FA模型对于选择育种种群和生产人口部署决策最有用。反应规范方法结合表型和环境数据,对G×E的环境驱动因素作了推论。适用于分析环境梯度逐渐不断变化的特征,需要明确界定的环境变量。稳定性分析对林木养殖具有直观的吸引力,有助于森林种植者识别稳定的基因型并减少长期风险。

所有这些统计方法为森林树种育种计划中高效稳定基因型的鉴定提供了解决方案。然而,追求性能的稳定性不能充分利用潜在的遗传在特定环境中获益。为了最大限度地发挥所有环境的潜在遗传增益,选择最适合特定环境的正确个体可能是最好的方法。大多数G×E林木育种研究调查了G×E的生长特征的模式和幅度。报道了这些特征的高水平G×E,特别是在辐射松,火炬松,桉树和杨树种。一些研究还研究了G×E的形态特征和木材性状特征的模式和幅度。本评价中的一半研究报告了G×E的高水平的形态性状,而木材性状没有报告G×E或最小G×E。总之,G×E可以通过多种方式进行量化,采取几种不同的方法可能会提供补充的见解,以帮助了解所涉及的相互作用,并确保最佳结果。对于新西兰和一般的林业,使用因子分析模型等新方法,育种价值估计显然是最有效的。为了理解和量化环境的详细作用,仍需要更多的工作。我们认为基因组技术(如Elshire等,2011; Neves et al。2013)和基因组序列注释(Ouzounis和Karp 2002; Wolf et al。2001)将提供更多信息,帮助解释遗传水平的原因和结果。为了最大限度地利用森林种植园的遗传收益和经济效益,应采用选择具体和已知环境的个人的策略。最终结果将是森林生产力和森林管理体系的进步,这将确保森林工业的长期可持续性和增加的利润。


  1. 我们区分表达水平与尺度效应的相互作用,为后者通过数据转换消除交互的子集的子集。

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 212,029评论 6 492
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,395评论 3 385
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 157,570评论 0 348
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,535评论 1 284
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 65,650评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,850评论 1 290
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,006评论 3 408
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,747评论 0 268
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,207评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,536评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,683评论 1 341
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,342评论 4 330
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,964评论 3 315
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,772评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,004评论 1 266
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,401评论 2 360
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,566评论 2 349

推荐阅读更多精彩内容