QTL×QTL×QTL Interaction Effects for Total Phenolic Content of Wheat Mapping Population of CSDH Lines under Drought Stress by Weighted Multiple Linear Regression

Abstract

This paper proposes the use of weighted multiple linear regression to estimate the triple3interaction (additive×additive×additive) of quantitative trait loci (QTLs) effects. The use of unweighted regression yielded an improvement (in absolute value) in the QTL×QTL×QTL interaction effects compared to assessment based on phenotypes alone in three cases (severe drought in 2010, control in 2012 and severe drought in 2012). In contrast, weighted regression yielded an improvement (in absolute value) in the evaluation of the aaa_gw parameter compared to aaa_p in five cases, with the exception of severe drought in 2012. The results show that by using weighted regression on marker observations, the obtained estimates are closer to the ones obtained by the phenotypic method. The coefficients of determination for the weighted regression model were significantly higher than for the unweighted regression and ranged from 46.2% (control in 2010) to 95.0% (control in 2011). Considering this, it is clear that a three-way interaction had a significant effect on the expression of quantitative traits.

1. Introduction

Approximately 80-95% of the fresh biomass of a plant is water, which plays an essential role in the plant growth, development and metabolism [1]. Drought is the most severe of all the environmental stresses affecting crop productivity and it belongs to the main factors restricting food production worldwide [2]. It is the single most significant threat to the future of global food security and has been a major famine trigger in the past. The unpredictable nature of drought depends on many factors: uneven and unreliable precipitation distribution, evapotranspiration and water-holding capacity around the rhizosphere [3,4]. Plants have evolved innate adaptations to stress conditions through a series of biochemical and physiological interventions involving the function of many stress-related genes.

Average growth is hampered by drought, which also alters the water cycle and makes plants less efficient water-users. This phenomenon is more complicated in plants since they exhibit a wide range of physiological and biochemical reactions at the cellular and organismal levels. Stomatal closure, membrane damage and disrupted activity of numerous enzymes, particularly those involved in adenosine triphosphate (ATP) generation, are the key factors reducing the rate of photosynthesis. Plants have a variety of defense mechanisms to tolerate drought. In response to drought stress, plants activate drought response mechanisms such as morphological and structural changes, the expression of drought-tolerance genes, synthesis of hormones and synthesis of osmoregulators to alleviate drought stress [5]. Those actions primarily aim to minimize further oxidative stress, water loss and other unfavorable effects of drought on plants. Secondary metabolites bind reactive oxygen species (ROS) to protect plants from oxidative damage under drought stress [6].

Phenols play multiple roles in plants, as they are structural components of cell walls, contribute significantly to growth and development, and improve plant tolerance to abiotic stresses [7]. These phenolic compounds are electron and proton sources of reactive oxygen species. Drought stress is associated with increased phenol accumulation and antioxidant activity. The accumulation of phenols under abiotic stress belongs to one of the defense mechanisms [8]. Salicylic acid (SA) is an essential phenolic compound currently recognized as a plant growth regulator that enhances plant responses to drought stress. SA is vital for transpiration, ion uptake, solute transport, photosynthesis, flower induction and protein synthesis [9].

Recent advances in genomics have enabled whole-genome sequencing of any genotype. These genomic references allow us to pinpoint each single nucleotide polymorphism’s (SNP’s) precise chromosomal location and the location is generated by sequence-based genotyping (GBS) [10]. This large number of SNPs is used in genome-wide association studies (GWAS) and quantitative trait locus (QTL) mapping to identify genomic regions or genes that may regulate target traits such drought tolerance, among others. By analyzing the genetics of complex traits, it is possible to better understand abiotic stress defense mechanisms and consequently create varieties with enhanced drought stress tolerance. The further inclusion of genetic interactions in such models can significantly increase the accuracy of such predictions [11].

Most of the important plant traits are quantitative and are influenced by many genes of quantitative trait loci (QTLs). The effect of these genes on the expression of the quantitative trait has been considered in various species, such as barley [12], maize [13], oilseed rape [14,15,16], rice [17,18], sugar beet [19] and wheat [20]. Many loci showing cumulative minor effects are frequently being weeded out in routine QTL analyses. However, when interacting with other loci, they can significantly affect the observed quantitative trait. In recent years, the effect of epistasis (double interaction) on the observed trait has been considered and reported for barley [21], maize [22], oilseed rape [23], pea [24], rice [25], soybean [26], tomato [27], Triticale [28] and wheat [29], among others. Several methods have been developed for gene localization and estimation of the epistasis effect [30,31,32,33,34,35,36]. Higher-order interactions (additive×additive×additive, additive×additive×dominance, additive×dominance×additive, dominance×additive×additive, additive×dominance×dominance, dominance×additive×dominance, dominance×dominance×additive and dominance×dominance×dominance) are overlooked, although the concept of genetic interaction has been known for more than a century [37]. It is difficult to imagine that a quantitative trait is not determined by the interaction of more than two genes.

In the case of homozygous material, only additive×additive×additive interaction can be estimated. The research reported in this article aims to estimate parameters connected with the additive×additive×additive (QTL×QTL×QTL) interaction effects for total phenolic content of wheat mapping population from the cross Chinese Spring×SQ1 doubled haploid (CSDH) lines under drought stress by weighted multiple linear regression. The phenotypic and genotypic methods are presented and the selected methods offer a novel approach to determining abiotic stress responses in plants.

2. Materials and Methods

2.1. Plant Material

The plant material used was the mapping population consisting of 94 doubled haploid lines of wheat (CSDH). The genetic map of the population, with a total length of about 4040 cM, contains 920 molecular markers. A detailed description of the mapping population and map is presented in Czyczyło-Mysza et al. [38]. The 94 CSDH lines were vernalized for 7 weeks at 4 °C with an 8/16 h light/dark photoperiod (short day). Germinated seedlings were then planted in 3 L pots (diameter 15 cm and height 20 cm) in 3 replicates per each line and per treatment (564 plants per each year, 3 for control and 3 for drought). Plants were grown in pots in the vegetation tunnel and protected against rain in conditions similar to natural, from spring to autumn (May–September). Limited irrigation was maintained for four weeks starting at the stage of tillering. A 20–25% field water capacity (FWC) was adopted as a severe drought, and 65–70% FWC as well-watered (control). Soil water content was measured with HydroSense Soil Water System CS630 (Campbell Scientific, Logan, UT, USA) using 12 cm probes and weight methods. Every three days, the weight of randomly selected pots was checked to determine the water volume for watering. In the last day of drought treatment, flag leaves were sampled for evaluation of total phenolics contents. Sampled leaves were homogenized in ethanol; the total phenolic contents were determined using the Folin–Ciocalteau method of Singleton and Rossi [39]. The absorbance at 765 nm was read using a microplate reader (Synergy 2, BioTek, Winooski, VT, USA) and chlorogenic acid was used as a standard. The analyses were carried out separately for all years for severe drought and control groups. The experiment was repeated during three planting seasons between 2010 and 2012. The growth conditions and all experimental procedures were maintained.

2.2. Estimation of the QTL×QTL×QTL Interaction Effects

2.2.1. Estimation Based on the Phenotype

The total additive×additive×additive interaction of the homozygous loci (three-way epistasis) effect based on phenotypic ($aa{a}_p$) observations of total phenolic content was estimated using the formula [11]:

$$\widehat{aa{a}_p}=\frac{1}{2}\left(L_{max}+L_{min}\right)-\overline{L},$$

(1)

where $L_{min}$ and $L_{max}$ are the doubled haploid lines with minimal and maximal mean value, respectively; $\overline{L}$ is the mean of all DH lines.

The test statistic to verify the hypothesis about $aa{a}_P$ different than zero is given by [40]

$$F_{aa{a}_p}=\frac{M{S}_{aa{a}_p}}{M{S}_e},$$

(2)

where $M{S}_{aa{a}_p}$ and $M{S}_e$ are mean squares for epistasis $aa{a}_p$ and residual, respectively. The number of genes (number of effective factors, $\widehat{K}$) obtained, based on phenotypic observations only, was calculated using the formula presented by Kaczmarek et al. [41]:

$$\widehat{K}=\frac{{\left(L_{max}-L_{min}\right)}^2}{4V_L},$$

(3)

where $V_L$ is additive variance.

2.2.2. Estimation Based on the Genotype

The estimation of aaa_g was based on the assumption that genes responsible for the total phenolic content were completely linked to observed molecular markers. After the selection of p markers (among all observed), determining the total phenolic content, the phenotypic observation model for DH lines is given as

$$\mathit{y}=\mathbf{1}\mu +\mathit{A}\mathit{\beta }+\mathit{E}\mathit{\gamma }+\mathit{T}\mathit{\delta }+\mathit{e},$$

(4)

where $\mathit{y}$ are the mean values of total phenolic content, 1 is the n-vector of ones, $\mu$ is the general mean, A is the (n×p)-matrix of the form $\mathit{A}=\left[\begin{array}{cc}\begin{array}{cc}{\mathit{m}}_{l_1}& {\mathit{m}}_{l_2}\end{array}& \begin{array}{cc}\cdots & {\mathit{m}}_{l_p}\end{array}\end{array}\right]$, l₁, l₂, …, l_p ∈ {1, 2, …, q}, $\mathit{\beta }$ is the p-vector of unknown parameters of the form ${\mathit{\beta }}^{\prime }=\left[\begin{array}{cc}\begin{array}{cc}{a}_{l_1}& a_{l_2}\end{array}& \begin{array}{cc}\cdots & a_{l_p}\end{array}\end{array}\right]$, E is the matrix for which columns are products of some columns of matrix A, $\mathit{\gamma }$ is the vector of unknown parameters of the form ${\mathit{\gamma }}^{\prime }=\left[\begin{array}{cc}\begin{array}{cc}a{a}_{l_1{l}_2}& a{a}_{l_1{l}_3}\end{array}& \begin{array}{cc}\cdots & a{a}_{l_{p-1}{l}_p}\end{array}\end{array}\right]$, T is the matrix for which columns are three-way products of some columns of matrix A, $\mathit{\delta }$ is the vector of unknown parameters of the form ${\mathit{\delta }}^{\prime }=\left[\begin{array}{cc}\begin{array}{cc}aa{a}_{l_1{l}_2{l}_3}& aa{a}_{l_1{l}_2{l}_4}\end{array}& \begin{array}{cc}\cdots & aa{a}_{l_{p-2}{l}_{p-1}{l}_p}\end{array}\end{array}\right]$, and e is the n-vector of random variables such that E(e_i) = 0, Cov(e_i, e_j) = 0 for i≠j, i, j = 1, 2, …, n. The parameters $a_{l_1}$, $a_{l_2}$, …, $a_{l_p}$ are the additive effects of the genes controlling the total phenolic content, parameters $a{a}_{l_1{l}_2}$, $a{a}_{l_1{l}_3}$, …, $a{a}_{l_{p-1}{l}_p}$ are the additive×additive interaction effects and parameters $aa{a}_{l_1{l}_2{l}_3}$, $aa{a}_{l_1{l}_2{l}_4}$, …, $aa{a}_{l_{p-2}{l}_{p-1}{l}_p}$ are the additive×additive×additive interaction effects. It was assumed that epistatic and three-way epistatic interaction effects show only loci with significant additive gene action effects. Consequently, this assumption results in a decrease in the number of possible significant effects and makes the regression model more useful. Loci with significant additive effects of genes were located and estimated previously by Czyczyło-Mysza et al. [42].

Unweighted Regression

The selection of markers chosen for model (4) can be made, e.g., by a stepwise regression procedure [43]. Here, a three-stage algorithm was used: first, the selection was made by a backward stepwise search independently inside all linkage groups; then, markers chosen in this way were put in one group and subjected to the second backward selection [44]; finally, the third stage considered situations where chosen markers were located on the chromosome very close to each other (closer than 5 cM). Because these markers are probably linked to one QTL, only the marker with the largest test statistic value was retained in the set. A critical significance level of 0.001, derived from the Bonferroni correction [45], was used for the first and second stages.

Denoting by ${\mathit{\alpha }}^{\prime }=\left[\begin{array}{cc}\begin{array}{cc}\mu & {\mathit{\beta }}^{\prime }\end{array}& \begin{array}{cc}{\mathit{\gamma }}^{\prime }& {\mathit{\delta }}^{\prime }\end{array}\end{array}\right]$ and $\mathit{G}=\left[\begin{array}{cc}\begin{array}{cc}\mathbf{1}& \mathit{A}\end{array}& \begin{array}{cc}\mathit{E}& \mathit{T}\end{array}\end{array}\right]$, the model was obtained:

$$\mathit{y}=\mathit{G}\mathit{\alpha }+\mathit{e}.$$

(5)

If G is of full rank, the estimate of ${\mathit{\alpha }}_{\mathit{u}}$ from traditional (unweighted) multiple linear regression model is given by [46]

$$\widehat{{\mathit{\alpha }}_{\mathit{u}}}={\left({\mathit{G}}^{\prime }\mathit{G}\right)}^{-1}{\mathit{G}}^{\prime }\mathit{y}.$$

(6)

The total three-way epistasis aaa_gu effect of genes influencing the total phenolic content from a traditional (unweighted) multiple linear regression model can be found as

$${\widehat{aaa}}_{gu}={{\displaystyle \sum }}_{k=1}^{p-2}{{\displaystyle \sum }}_{\begin{array}{c}{k}^{\prime }=2\\ k^{\prime }\ne k\end{array}}^{p-1}{{\displaystyle \sum }}_{\begin{array}{c}{k}^{″}=3\\ k^{″}\ne k^{\prime }\end{array}}^p{\widehat{aaa}}_{l_k{l}_{k^{\prime }}{l}_{k^{″}}}.$$

(7)

Weighted Regression

The modified version of total phenolic content regression on marker data is considered by taking a weighted multiple linear regression, that is, a regression with a diagonal matrix W of unknown variances of observations, which, however, may be empirically found by estimation. The selection of markers for the weighted regression is made by a similar method described for the unweighted case but concerns the weights for lines. In this model, the estimate of ${\mathit{\alpha }}_{\mathit{w}}$ is

$$\widehat{{\mathit{\alpha }}_{\mathit{w}}}={\left({\mathit{G}}^{\prime }{\mathit{W}}^{-1}\mathit{G}\right)}^{-1}{\mathit{G}}^{\prime }{\mathit{W}}^{-1}\mathit{y},$$

(8)

where $\mathit{W}=\left(w_{ii}\right)$ with $w_{ii}$ is the estimated variance for i-th DH ine, i = 1, 2, …, n. The selection of markers for the weighted regression is made by the same method as described for the unweighted case.

The total three-way epistasis aaa_gw effect of genes influencing the total phenolic content from the weighted multiple linear regression model can be found as

$${\widehat{aaa}}_{gw}={{\displaystyle \sum }}_{k=1}^{p-2}{{\displaystyle \sum }}_{\begin{array}{c}{k}^{\prime }=2\\ k^{\prime }\ne k\end{array}}^{p-1}{{\displaystyle \sum }}_{\begin{array}{c}{k}^{″}=3\\ k^{″}\ne k^{\prime }\end{array}}^p{\widehat{aaa}}_{l_k{l}_{k^{\prime }}{l}_{k^{″}}}.$$

(9)

The coefficients of determination were used to measure how both models (unweighted and weighted) fitted the data and were the amounts of the phenotypic variance explained by the total of three interactive models.

All analyses were conducted using the GenStat v. 22 statistics software.

3. Results

3.1. Estimation Based on the Phenotype

Phenotypic estimates of the total additive×additive×additive effect (aaa_p) are presented in

Table 1. Minimal and maximal values of average for lines, means for all doubled haploid lines for total phenolic content, coefficient of variation, phenotypic estimates of the total additive×additive×additive effect (aaa_p) and the number of genes (the number of effective factors).

Year	2010		2011		2012
Stress	Control	Severe Drought	Control	Severe Drought	Control	Severe Drought
Minimal	7.104	6.499	7.164	6.231	7.645	4.626
Maximal	14.414	12.682	16.321	15.070	13.781	13.023
Mean	10.450	9.435	11.993	10.214	10.736	8.204
Coefficient of variation	18.06	18.72	15.39	15.51	14.34	23.63
aaap	0.309 *	0.156 *	−0.250 *	0.437 **	−0.023	0.620 ***
The number of genes (the number of effective factors)	4.617	3.824	5.050	6.021	4.164	4.846

* p < 0.05; ** p < 0.01; *** p < 0.001.

. In four of the six cases, the total aaa_p effect was positive. A negative effect was observed for the 2011 and 2012 control groups. The highest total additive×additive×additive effect was observed for the 2012 severe drought group and the lowest for the 2011 control group. Triples interaction effects were statistically significant for all groups except the control group in 2012.

Minimal and maximal values of averages for CSDH lines as well as means for all doubled haploid lines for total phenolic content were larger for the control groups than for the severe drought groups in all three years of study. The coefficients of variation were larger for severe drought groups than for the control groups in all cases. The coefficient of variation values ranged from 14.34 for the 2012 control group to 23.63 for the 2012 severe drought group, indicating low data spread relative to the mean values (Table 1). The number of genes (effective factors) varied between groups, ranging from 3.824 for the 2010 severe drought group to 6.021 for the 2011 severe drought group. The difference between the highest and the lowest number of genes in all groups was equal to 2.197; between the control and severe drought groups measured in the same year, the highest difference was observed in the 2011 group valued at 0.971 (Table 1).

3.2. Estimation Based on the Genotype

Genotypic estimates of the total additive×additive×additive effect calculated based on unweighted (aaa_gu) and weighted (aaa_gw) multiple linear regression are presented in

Table 2. Genotypic estimates of the total additive×additive×additive effects estimated based on unweighted (aaa_gu) and weighted (aaa_gw) multiple linear regression.

Year			2010		2011		2012
Stress			Control	Severe Drought	Control	Severe Drought	Control	Severe Drought
Unweighted	QTLs number		6	5	4	6	7	6
	Number of aaagu		8	1	0	2	14	8
	aaagu effects	Min.	−0.288	0.350		−0.489	−0.410	−0.138
		Max.	0.533	0.350		0.569	0.538	0.409
		Total	−0.043	0.350		0.110	0.170	0.688
	R2 [in %]		32.4	42.0		36.3	39.1	44.5
Weighted	QTLs number		16	15	20	14	26	14
	Number of aaagw		3	5	9	5	6	6
	aaagw effects	Min.	−0.659	−0.494	−0.607	−0.499	−0.526	−1.124
		Max.	0.598	0.612	0.548	0.368	0.557	1.182
		Total	0.443	−0.265	1.178	−1.302	−0.184	−0.228
	R2 [in %]		46.2	61.4	95.0	58.8	78.9	58.5

.

3.2.1. Unweighted Regression

The number of QTLs detected using the traditional (unweighted) gene localization method (Table 2) was similar to the number of effective factors (Table 1) and ranged from four (control, 2011) to seven (control, 2012). The number of significant threes of QTL×QTL×QTL genes varied widely and ranged from zero (control, 2011) and one (severe drought, 2010) to fourteen (control, 2012). The range of effect values when the observed trait was determined by at least two threes was 0.547 (severe drought, 2012) to 1.058 (severe drought, 2011 (Table 2 and

Table 3. Genotypic estimates of the additive×additive×additive interaction effects for the particular QTL×QTL×QTL threes based on unweighted (aaa_gu) and weighted (aaa_gw) multiple linear regression.

Year			2010				2011				2012
Stress			Control		Severe Drought		Control		Severe Drought		Control		Severe Drought
QTL1 (LG)	QTL2 (LG)	QTL3 (LG)	aaagu	aaagw	aaagu	aaagw	aaagu	aaagw	aaagu	aaagw	aaagu	aaagw	aaagu	aaagw
dupw004b (1A)	wPt-3094 (3A)	wmc1 (3B)				0.471
dupw004b (1A)	wPt-668160 (1A)	gwm635b (7A)				0.612
wPt-3094 (3A)	gwm165.3 (4A)	wPt-0391 (4B)			0.350
wPt-3094 (3A)	gwm635b (7A)	wPt-671748 (7D)				−0.494
cfa2262 (3A)	dupw004a (4A)	blt101.t7 (4D)								−0.499
blt4.1 (3B)	wPt-4048 (3B)	gwm269.2 (4A)											−0.138
blt4.1 (3B)	gwm269.2 (4A)	gwm165.3 (4A)											0.049
wPt-6239 (3B)	blt101.t7 (4D)	wPt-2697 (5A)								−0.337
wPt-6239 (3B)	blt101.t7 (4D)	wmc83 (7A)								−0.348
rPt-8896 (3B)	wmc1 (3B)	gwm52 (3D)				−0.376
wPt-0021 (3B)	gwm52 (3D)	barc60 (4B)		−0.659
gwm269.2 (4A)	gwm165.3 (4A)	gwm205 (5D)	0.008
gwm269.2 (4A)	gwm165.3 (4A)	barc44 (5D)	−0.152
gwm269.2 (4A)	gwm165.3 (4A)	m69p78.1 (7A)	−0.100
gwm269.2 (4A)	gwm205 (5D)	barc44 (5D)	0.063
gwm269.2 (4A)	gwm205 (5D)	m69p78.1 (7A)	−0.263
gwm269.2 (4A)	barc44 (5D)	wPt-9834 (5A)	0.533
gwm165.3 (4A)	gwm165.3 (4A)	wPt-667091 (7D)									0.158
gwm165.3 (4A)	wPt-0391 (4B)	wPt-3883 (7A)									0.026
gwm165.3 (4A)	wPt-0391 (4B)	wPt-667091 (7D)									−0.179
gwm165.3 (4A)	wPt-0391 (4B)	wmc243b (7D)									0.247
gwm165.3 (4A)	gwm205 (5D)	m69p78.1 (7A)	0.156
gwm165.3 (4A)	barc44 (5D)	m69p78.1 (7A)	−0.288
gwm165.3 (4A)	gwm174 (5D)	wPt-3883 (7A)									0.026
gwm165.3 (4A)	gwm174 (5D)	wPt-667091 (7D)									0.031
gwm165.3 (4A)	gwm174 (5D)	wmc243b (7D)									−0.161
gwm165.3 (4A)	wPt-3883 (7A)	wmc243b (7D)									−0.008
wPt-0391 (4B)	wPt-0391 (4B)	wPt-667091 (7D)									−0.410
wPt-0391 (4B)	gwm174 (5D)	wPt-3883 (7A)									0.128
wPt-0391 (4B)	gwm174 (5D)	wmc243b (7D)									0.538
wPt-0391 (4B)	wPt-3883 (7A)	wPt-667091 (7D)									−0.151
barc152 (1B)	wPt-6239 (3B)	gwm191b (3D)										−0.400
barc152 (1B)	m65p64.8_4B (4B)	wPt-8149 (7A)										−0.526
barc152 (1B)	m92p78.10 (2A)	m60p64.13_3B (3B)										0.557
barc152 (1B)	wPt-8072 (2B)	gwm165.3 (4A)										−0.493
m65p64.8a (5B)	gwm271 (5B)	m69p78.1 (7A)		0.504
gwm174 (5D)	wPt-3883 (7A)	wPt-667091 (7D)									−0.130
gwm174 (5D)	wPt-3883 (7A)	wmc243b (7D)									0.055
psr648_1B (1B)	wmc181 (2A)	gwm285 (3B)		0.598
m17p65.1 (1B)	cfd73 (2D)	wmc468 (4A)				−0.478
wmc432 (1D)	wPt-3738 (1D)	psp2151.3 (2A)												−0.640
wPt-3738 (1D)	tPt-0202 (3A)	gwm269.2 (4A)												−1.124
wPt-3738 (1D)	tPt-0202 (3A)	dupw004a (4A)												1.182
wPt-3738 (1D)	gwm513 (4B)	wmc157 (7D)												0.387
wPt-3738 (1D)	wmc429 (1D)	tPt-0202 (3A)												0.354
wmc429 (1D)	cfd11 (2D)	wPt-10291 (3D)						−0.505
wmc429 (1D)	gwm30.1 (2D)	tPt-0202 (3A)						0.481
wmc429 (1D)	wPt-7466 (2D)	wPt-9749 (2D)						−0.607
wmc429 (1D)	gwm269.2 (4A)	dupw004a (4A)												−0.387
wmc429 (1D)	wPt-6316 (1D)	wPt-7466 (2D)						0.548
wmc429 (1D)	wPt-6316 (1D)	wPt-10291 (3D)						0.436
wmc429 (1D)	wPt-6316 (1D)	gwm161 (3D)						0.493
wmc429 (1D)	wPt-6316 (1D)	gwm165.3 (4A)						−0.511
wmc429 (1D)	wPt-6316 (1D)	gwm639.2 (5B)						0.419
wmc429 (1D)	wPt-732556 (1D)	gwm30.1 (2D)						0.424
m92p78.10 (2A)	wPt-6239 (3B)	gwm161 (3D)										0.360
m92p78.10 (2A)	wPt-6239 (3B)	wPt-8149 (7A)										0.317
wmc453a (2A)	wPt-3883 (7A)	wPt-8919 (7B)							0.569
wmc453a (2A)	wmc283.1 (7A)	wPt-8919 (7B)							−0.459
psp2151.3 (2A)	blt4.1 (3B)	wPt-4048 (3B)											0.201
psp2151.3 (2A)	blt4.1 (3B)	gwm269.2 (4A)											0.136
psp2151.3 (2A)	blt4.1 (3B)	gwm165.3 (4A)											−0.095
psp2151.3 (2A)	wPt-4048 (3B)	gwm269.2 (4A)											0.409
psp2151.3 (2A)	wPt-4048 (3B)	gwm165.3 (4A)											0.091
psp2151.3 (2A)	gwm269.2 (4A)	gwm165.3 (4A)											0.035
wPt-3949 (2B)	wmc257 (2B)	gwm165.3 (4A)								0.368
wmc257 (2B)	wPt-2697 (5A)	gwm292_5D (5D)								−0.485

)). The total aaa_gu effect ranged from –0.043 (control, 2010) to 0.688 (severe drought, 2012) (Table 2). The phenotypic variation of total phenolic content explained by the total triples interaction was similar in all five cases where the number of triples interactions was greater than zero and ranged from 32.4% (control, 2010) to 44.5% (severe drought, 2012) (Table 2).

3.2.2. Weighted Regression

The use of weighted regression resulted in a definite increase in the number of QTLs detected compared to when classical unweighted regression was used. The increase ranged from 233% (from 6 to 14 in severe drought, 2011 and 2012) to 500% (from 4 to 20 in the control group in 2011) (Table 2). The number of threes using weighted regression was higher than using unweighted regression in three cases (severe drought, 2010; control 2011 and severe drought 2011). In the remaining cases (control, 2010; control, 2012 and severe drought 2012), it was smaller (Table 2). The spread of estimated effect values was larger in all cases except severe drought in 2011, where it decreased from 1.058 to 0.867 (Table 2 and Table 3). The coefficients of determination for the weighted regression model were significantly higher than for the unweighted regression and ranged from 46.2% (control in 2010) to 95.0% (control in 2011) (Table 2).

Genotypic estimates of the additive×additive×additive interaction effects for the particular QTL×QTL×QTL threes based on unweighted (aaa_gu), and weighted (aaa_gw) multiple linear regression are presented in Table 3. Sixty-seven statistically significant triples interactions were observed. None of the QTL×QTL×QTL threes were significant in at least two cases: threes significant using weighted regression were not significant using unweighted regression and vice versa (Table 3). The QTLs most frequently found in triples interaction were gwm165.3 (22 times), gwm269.2 (11 times) and wmc429 (11 times). Using unweighted regression, the genes most frequently found in triples interaction were gwm165.3 (19 times), gwm269.2 (11 times) and wPt-0391 (9 times) (Table 3). Using weighted regression, the genes most frequently observed in triples interaction were wmc429 (eleven times), wPt-3738 (five times), wPt-6239 (five times) and wPt-6316 (also five times) (Table 3).

In 43 triples interactions, each QTL was located on a different linkage group (LG). In twenty-four cases, two QTLs (out of three) were located on one LG. There was no situation in which all three QTLs from a threesome were located on a single LG (Table 3). Using unweighted regression, 23 cases were observed with each of the QTLs of the 3 being at a different LG. In contrast, in ten cases, two QTLs were located on a single chromosome. Using weighted regression, 20 cases were observed where each QTL of the significant 3 was localized in a different linkage group. In contrast, in 14 cases, 2 QTLs were localized to a single LG (Table 3). The distribution of QTLs included in the triples interaction on each linkage group varied widely, ranging from 2 (on 7B) to 39 (on 4A) (Table 3).

4. Discussion

The aim of a breeding process is to obtain new genotypes with improved traits compared to the parental forms. The parameter connected with the additive×additive×additive interaction gene action can influence decisions about the usefulness of the breeding material for that purpose. A numerical comparison of the three methods of estimation of the additive×additive×additive interaction effects are presented in this paper. The comparison was conducted on 94 doubled haploid lines (DHLs) generated from the cross between hexaploid wheat (Triticum aestivum L.) genotypes Chinese Spring (CS) and SQ1 (a high abscisic acid breeding line) according to Quarrie et al. [47]. The total phenolic contents were measured separately for three years (2010–2012) for severe drought and control groups. The present results demonstrated the use of weighted regression to determine the triplets of QTLs and estimate the effects of their QTL×QTL×QTL interaction. To the best of our knowledge, there is no published evidence on the use of weighted regression to assess the QTL×QTL×QTL triple interaction. Most studies focus on the analysis of single genes as determinants of a quantitative trait. Weighted regression was used to estimate QTL×QTL epistasis interaction [48,49,50,51,52]. However, higher-order interactions are ignored. Therefore, the novelty of this proposal is in determining how to select triples of QTLs and estimate their effects using weighted regression.

Previously, in [11], we showed that QTL×QTL×QTL interaction effects assessments obtained with molecular marker observations (using traditional unweighted regression) had smaller absolute values than assessments based on phenotypic observations alone. The differences obtained were so large that it was decided to propose a new method with a chance to improve QTL×QTL×QTL assessments.

The present results obviously demonstrated that the detected QTL×QTL×QTL affected the total phenolic content. An important assumption to make is that aaa interaction effects show only loci connected to markers with significant effects. However, it should be emphasized that the selection of markers with significant effects was carried out using weighted regression. The use of unweighted regression yielded an improvement (in absolute value) in the effect of aaa_gu compared to assessment based on phenotypes alone in three cases (severe drought in 2010, control in 2012 and severe drought in 2012). In contrast, weighted regression yielded an improvement (in absolute value) in the evaluation of the aaa_gw parameter compared to aaa_p in five cases, with the exception of severe drought in 2012 (Table 1 and Table 2).

Considering the effect of three-way interaction on the expression of quantitative traits may also be important in the context of hypostasis, a situation in which the phenotype is altered by the expression of an allele at a separate locus, in an epistasis event [53]. Genetic interaction analysis may be a powerful genetic strategy for analyzing the fitness and phenotypes of genotypes to dissect the interactions between genes, reason genes into biological pathways, and characterize genes of unknown operation. Omitting to consider the effect of higher-order interactions in quantitative genetics can result in a lack of completeness in determining a quantitative trait. Cyplik et al. [39], in their study of maize inbred lines, observed the significance of a triple interaction in the absence of the epistasis effect. This phenomenon occurred in three situations: for cob diameter (twice) and for core length. The reverse situation (significant epistasis and non-significant triple interaction) did not occur in any case.

5. Conclusions

The results show that by using weighted regression on marker observations, the obtained estimates are closer to the ones obtained by the phenotypic method. Considering this, it is clear that a three-way interaction had a significant effect on the expression of quantitative traits. The proposed weighted regression method for estimation of the QTL×QTL×QTL interaction effect can be a bridge between the phenotypic method and the genotypic method.