A Comparison of Methods to Estimate Additive–by–Additive–by–Additive of QTL×QTL×QTL Interaction Effects by Monte Carlo Simulation Studies

The goal of the breeding process is to obtain new genotypes with traits improved over the parental forms. Parameters related to the additive effect of genes as well as their interactions (such as epistasis of gene–by–gene interaction effect and additive–by–additive–by–additive of gene–by–gene–by–gene interaction effect) can influence decisions on the suitability of breeding material for this purpose. Understanding the genetic architecture of complex traits is a major challenge in the post-genomic era, especially for quantitative trait locus (QTL) effects, QTL–by–QTL interactions and QTL–by–QTL–by–QTL interactions. With regards to the comparing methods for estimating additive–by–additive–by–additive of QTL×QTL×QTL interaction effects by Monte Carlo simulation studies, there are no publications in the open literature. The parameter combinations assumed in the presented simulation studies represented 84 different experimental situations. The use of weighted regression may be the preferred method for estimating additive–by–additive–by–additive of QTL–QTL–QTL triples interaction effects, as it provides results closer to the true values of total additive–by–additive–by–additive interaction effects than using unweighted regression. This is also indicated by the obtained values of the determination coefficients of the proposed models.


Introduction
The goal of breeding is to obtain the most economically advantageous varieties. A very important aspect is to shorten the breeding process without losing the quality of the obtained genotypes. Knowledge of the genetic structure and genes actions can contribute to genetic advances resulting in a shorter breeding process. In traditional quantitative genetics, organisms are studied on the basis of only phenotypic observations, and conclusions about the genotype are drawn based on the ways in which quantitative traits are inherited, among other things. This is usually realized by describing the action of genes using genetic parameters, mean and variance functions of phenotypes [1]. Among these parameters is the additive effect of genes and interactions between them [2]. Additive effects are perpetuated in a population as it increases its homozygosity in subsequent generations. Thus, the significant additive effect of genes in the population means that selection starting in early generations offers hope of obtaining transgenic homozygous lines [3].
The introduction of different types of markers into genetic and breeding research has revolutionized these scientific disciplines. Markers, among other molecular markers, are used in selection and for genetic transformation [4]. Molecular markers show genetic differences between individuals with opposite traits in the form of nucleotide sequences and have been developed and used to speed up the breeding process in many crops [5]. Identifying markers significantly associated with plant traits can help develop new varieties with desirable traits. It is very important to localize the genes or groups of genes controlling a trait (that is, to give the linkage distances of the genes studied to the markers) and determine their effects [6]. The accuracy and correctness of the conclusions depend largely on the precision of determining the location of the gene on the chromosome and the estimation of parameters associated with gene effects [7].
Most important plant traits are quantitative in nature and are influenced by multiple quantitative trait locus (QTLs). The effects of these QTLs on the expression of quantitative trait has been considered in various species such as maize [8], barley [9], rice [10,11], oilseed rape [12][13][14], wheat [15] and sugar beet [16]. Many loci showing cumulative small effects are often rejected in standard QTL analyses [17,18]. However, by interacting with other loci, they can significantly affect the observed quantitative trait. The effect of epistasis (double interaction) on the observed trait has been considered and reported for, among others: wheat [19], Triticale [20], tomato [21], soybean [22], rice [23], pea [24], canola [25], maize [26] and barley [27]. However, higher-order interactions are omited, although it is difficult to imagine that a quantitative trait is not determined by the interaction of more than two genes. Based on the knowledge of physiological genetics and biochemistry, interactions between gene products are ubiquitous [28].
A previous report [29] described the results obtained for analytical and numerical comparisons of two methods for estimating the parameter associated with the additiveby-additive-by-additive (aaa) interaction effect, namely, the phenotypic method, based on extreme groups of homozygous lines, and the genotypic method, based on marker observations. One of the conclusions of these studies was that the estimate of the total additive-by-additive-by-additive interaction effect based on the marker data was in most cases smaller than that based on the phenotype. The explanation for this phenomenon may be simple. Phenotypic data can be used to estimate only the total three-way interaction effect of all hypothetical gene triplets determining a trait. Using marker data, which can be more or less precisely mapped into the genome, individual gene-by-gene-by-gene interaction effects can be estimated. For practical reasons, the number of QTL-QTL-QTL interactions is deliberately low. The sum of the effects obtained is lower than the phenotypic estimate, and this difference may be reinforced by the lack of markers in the regions where the genes are located.
In addition to the above explanation, other possible sources of differences between the calculated estimates should be considered. The above results concern QTL-QTL-QTL interaction effects obtained by the simplest possible methods, i.e., using multiple linear regression on the marker data. Cyplik et al. [30] showed that modifying this regression by using empirical weights provides a better agreement between phenotypic and genotypic estimates.
The purpose of this paper was to compare two methods for estimating the parameter associated with additive-by-additive-by-additive interaction effects: the phenotypic method, traditionally used in quantitative genetics, and the genotypic method, which is based on marker observation and is now routinely used in many species. The comparison was made using a Monte Carlo simulation study. Consideration was also given to modifying the estimation of the additive-by-additive-by-additive of QTL-QTL-QTL interaction effects by using weighted multiple linear regression (in two variants) to determine a bridge between the two compared methodologies.
In the approach described here, the estimation of the parameter associated with the additive-by-additive-by-additive interaction effect was based on extreme groups of homozygous lines and on data for genotypic markers. Weighted multiple linear regression was used to estimate the QTL-QTL-QTL interaction effects, using standard weighted regression with weights based on estimated variances for lines in two variants: using weighted regression to estimate only QTL-QTL-QTL triples interaction effects (unweighted regression for QTL selection and QTL-QTL epistasis) and using weighted regression on all three stages-for QTL selection, epistatic pairs and QTL-QTL-QTL triples. The effectiveness of the proposed method was tested using Monte Carlo simulations for different experimental variants, such as the number and location of chromosomes of QTL-QTL-QTL triples determining the quantitative trait (Table 1), the presence or absence of QTL-QTL epistasis ( Figure 1). Based on the knowledge of the literature, this is the first report of a comparison of methods for estimating a parameter related to the QTL-QTL-QTL triples interaction effect based on simulation studies.  (Table 1), the presence or absence of QTL-QTL epistasis ( Figure 1). Based on the knowledge of the literature, this is the first report of a comparison of methods for estimating a parameter related to the QTL-QTL-QTL triples interaction effect based on simulation studies.

Results
Difficult experimental situations were examined in this large-scale simulation. Tables 2-9 show the results of the simulation performed to compare the estimates of additive-byadditive-by-additive interaction effects obtained by the four methods: the phenotypic methods and three genotypic methods (one unweighted and two weighted). Tables 2-5 show the  results obtained with an assumed error variance of 5, while Tables 6-9 show the results when the error variance was 10. Tables 2, 3, 6 and 7 contain results assuming that the total effect of epistasis was equal to 0, and Tables 4, 5, 8 and 9, assume that the effect of epistasis was significantly different from zero. The results of the simulation studies assuming that all the effects of the aaa triples were equal are shown in Tables 2, 4, 6 and 8; while Tables 3, 5, 7 and 9 contain the results obtained assuming different effects of the aaa triples. Table 2. Phenotypic and genotypic estimates of the total additive-additive-additive effect, mean squared errors and mean of coefficients of determination R 2 j obtained in the simulation study for: the error variance was equal to 5, no epistasis, equal effects for all triples.  Table 3. Phenotypic and genotypic estimates of the total additive-additive-additive effect, mean squared errors and mean of coefficients of determination R 2 j obtained in the simulation study for: the error variance was equal to 5, no epistasis, different effects for triples.  Table 4. Phenotypic and genotypic estimates of the total additive-additive-additive effect, mean squared errors and mean of coefficients of determination R 2 j obtained in the simulation study for: the error variance was equal to 5, presence of epistasis, equal effects for all triples.  Table 5. Phenotypic and genotypic estimates of the total additive-additive-additive effect, mean squared errors and mean of coefficients of determination R 2 j obtained in the simulation study for: the error variance was equal to 5, presence of epistasis, dif-ferent effects for triples.  Table 6. Phenotypic and genotypic estimates of the total additive-additive-additive effect, mean squared errors and mean of coefficients of determination R 2 j obtained in the simulation study for: the error variance was equal to 10, no epistasis, equal effects for all triples.  Table 7. Phenotypic and genotypic estimates of the total additive-additive-additive effect, mean squared errors and mean of coefficients of determination R 2 j obtained in the simulation study for: the error variance was equal to 10, no epistasis, different effects for triples.   Table 9. Phenotypic and genotypic estimates of the total additive-additive-additive effect, mean squared errors and mean of coefficients of determination R 2 j obtained in the simulation study for: the error variance was equal to 10, presence of epistasis, different effects for triples.

Phenotypic Estimates
The results of the simulation showed that the phenotypic estimate was always greater than 15, the true-assumed value, except in four cases for five QTL-QTL-QTL triples: (1) the number of chromosomes containing QTLs was 1, 2, the error variance was equal to 5, the effect of epistasis was different from zero and the effects of aaa triples were different (Table 5), (2) the number of chromosomes containing QTLs was 1, 2, the error variance was equal to 10, the effect of epistasis was equal to zero and the effects of all aaa triples were equal (Table 6), (3) and (4) the number of chromosomes containing QTLs was 1, 2 (variant code-C10), and 2 (variant code-C11), respectively, the error variance was equal to 10, the effect of epistasis was different from zero, while the effects of aaa triples were different ( Table 9). The largest values were obtained for one QTL-QTL-QTL triple, also the smallest values for five QTL-QTL-QTL triples (taking into account the analogous distribution of triples in the linkage group). The largest values were obtained for the situation when the QTLs constituting the triple were located in different chromosomes, while the smallest values were for the situation when all three QTLs constituting the triple were located in one linkage group. The phenotypic estimates were larger when assuming no epistasis than when assuming the effect of epistasis on the values of the observed quantitative trait. When considering the effects of triples, equal or different, similar values of phenotypic estimates were obtained. The aaa p values were larger when the error variance was equal to 10 than for an error variance of 5, except when five triples were assumed.

Genotypic Estimates
The genotypic estimates of aaa effects were greater than true-assumed value 15 for 1, 2 and 3 QTL-QTL-QTL triples. For situations where five QTL-QTL-QTL triples were assumed in most cases (58.33%), genotypic estimates were less than 15. With an assumed error variance of 5, in 19 cases (51.78%) genotypic estimates were less than the true-assumed value of aaa equal to 15 (Tables 6-9). The largest values of genotypic estimates (unweighted and weighted) were obtained for one QTL-QTL-QTL triple, also the smallest values for five QTL-QTL-QTL triples. The largest values were obtained for situations when QTLs constituting triples were located in different chromosomes, except for the case when five QTL-QTL-QTL triples with different effects were assumed in the absence of epistasis and an error variance of 10 (Table 8). In contrast, the smallest values for cases when all three QTLs constituting the triple were placed in a single coupling group. In most cases, genotypic estimates were larger for different effects than for equal effects of triples. The genotypic estimates were larger with an assumed larger (10) error variance than with a smaller (5) variance in situations where the quantitative trait was determined by 1, 2 or 3 QTL-QTL-QTL triples. In contrast, when five QTL-QTL-QTL triples were assumed, genotypic estimates were larger for an error variance of 5 than for an error variance of 10.

Phenotypic vs. Genotypic Estimates
The differences between phenotypic estimates and the genotypic estimates for unweighted regression were always positive and were generally smaller when different effects of QTL-QTL-QTL triples were assumed. The relationships aaa p > aaa gw1 > aaa gw2 > aaa gu were generally observed, except in three cases: (1) when two QTL-QTL-QTL triples were located in three chromosomes for variance error equal to 5, lack of epistasis and different effects of triples ( aaa gw1 = aaa gw2 = 17.99, Table 3), (2) when three QTL-QTL-QTL triples were located in two linkage groups, variance error equal to 5, lack of epistasis and different effects of triples ( aaa gw1 > aaa p , Table 5), and (3) when three QTL-QTL-QTL triples with equal effects were located in two linkage groups and the effect of epistasis was different from zero: in which case aaa gw1 > aaa gw2 > aaa p ( Table 8). The weighted estimates were always closer to the phenotypic estimates than the unweighted estimates. The weighted method way (i) showed the best (closest to the true assumed value of 15) ratings of aaa triples.

Mean Squared Errors
In general, a decrease in the estimates was accompanied by an increase in their mean squared error (Tables 2-9). The relationships of mean squared error for aaa p < aaa gw1 < aaa gu < aaa gw2 were generally observed, except when two QTL-QTL-QTL triples with different effects were assumed to be present in two linkage groups in case of lack of epistatic effect and variance error was equal to 5, in which case mean squared error for aaa gw1 > aaa gw2 ( Table 3). The mean squared error was larger with an assumed error variance of 10 than with a smaller assumed error variance. The effect values of each triple (equal or different) did not affect the mean squared error values. On the other hand, assuming the presence or absence of the epistasis effect affected the mean squared error values for aaa: large mean squared error values were obtained when assuming no epistasis. The smallest mean squared errors were obtained when the trait was determined by one QTL-QTL-QTL triple, while the largest was obtained with five QTL-QTL-QTL triples. In general, the relationships of mean squared error values for assumed 1QQQ < 2QQQs < 3QQQs < 5QQQs were observed. The smallest mean squared error was obtained for QTL-QTL-QTL triples located in one linkage group, while the largest was obtained when QTL-QTL-QTL triples were located across three chromosomes.

Coefficients of Determination
Variance explained by the QTL-QTL-QTL triple interactions ranged from 83% (Tables 3 and 5) to 97% ( Table 8). Coefficients of determination for one QTL-QTL-QTL triple assumed ranged from 0.87 to 0.95, with a mean of 0.91. Coefficients of determination for two QTL-QTL-QTL triples assumed ranged from 0.83 to 0.97, with a mean of 0.89. Coefficients of determination for three and five QTL-QTL-QTL triples assumed were similar and ranged from 0.85 and 0.84 to 0.95 and 0.95, with mean of 0.90 and 0.90, for three and five QTL-QTL-QTL triples assumed, respectively.

Discussion
The goal of the breeding process is to obtain new genotypes with traits improved over the parental forms [31]. Parameters related to the additive effect of genes as well as their interactions (such as additive-by-additive epistasis of gene-by-gene interaction effect and additive-by-additive-by-additive of gene-by-gene-by-gene interaction effect) can influence decisions on the suitability of breeding material for this purpose [32][33][34][35][36][37][38]. Understanding the genetic architecture of complex traits is a major challenge in the postgenomic era, especially for QTL effects, QTL-by-QTL interactions [39] and QTL-by-QTLby-QTL interactions.
Monte Carlo simulations are a research tool very widely used to solve various problems by obtaining approximate results. Estimation of epistasis effects using simulation studies has been presented, for example, by Bocianowski [40], Ahsan et al. [41], de los Campos et al. [42], Wang et al. [43] and Sharma et al. [44]. With regards to the comparing methods for estimating additive-by-additive-by-additive of QTL×QTL×QTL interaction effects by Monte Carlo simulation studies, there are no publications in the open literature. The paper describes the first simulation study comparing methods for estimating the total additive-by-additive-by-additive effect of gene-by-gene-by-gene interactions. In addition, weighted multiple linear regression was used to estimate QTL-by-QTL-by-QTL interaction effects, with two variants proposed.
The Monte Carlo simulation studies conducted to compare estimation methods could not, of course, take into account all possible experimental situations. It is essential to note that these simulations may not fully capture the complexities and nuances of real-world genetic interactions. The findings should be interpreted with caution until further empirical validation is conducted. However, the parameter combinations used in the analyses presented here correspond to the cases most often encountered in actual QTL studies [45][46][47][48][49][50][51][52][53][54][55][56][57][58][59][60][61][62]. In simulation studies, the number of lines analyzed is most often assumed to be as high as 500 or 2000 [45,46]. Viana et al. [46] analyzed 400 plants, genotyped for 975 SNPs distributed in 10 chromosomes of 100 cm. One thousand SNP markers were considered in their study by Crawford et al. [47]. The ten chromosomes include Anthoxanthum odoratum, Brassica campestris, Sorghum sudanense, Zea maize. Practical experience shows that the number of genes located on individual chromosomes varies greatly [48][49][50][51]. The number of QTLs assumed in the present study was 14. The same number of QTLs was presented in their studies by Wu [52] and Balestre et al. [53]. The assumed number of five epistatic pairs was a result observed in many practical experiments [54][55][56][57]. Simulation studies based on 10,000 generations were also presented by Goutelle et al. However, it was found that the number of assumed QTL-QTL-QTL triplets, their assumed positions, in one, two or three chromosomes, and their effects (unequal or equal) affected the estimates and their comparison. Phenotype estimation decreased toward the true value as the number of presumed QTL-QTL-QTL triplets increased (Figure 2). For the largest number of QTL-QTL-QTL triplets, it may have been less than the true value, as genotypic combinations generating extreme lines were not adequately represented in the simulated sample. The phenotypic estimate also tended to the true value when the number of chromosomes with QTL-QTL-QTL triplets decreased, i.e., when there was an accumulation of coupled individual effects. Similarly, phenotypic estimation tended to have a true value when the assumed epistasis effect was different from zero. The phenotypic evaluation was overestimated for a larger assumed error variance. The genotypic estimate decreased its value in a very similar way to the phenotypic estimate. The difference between phenotypic and genotypic estimation was smallest when five QTL-QTL-QTL triplets were assumed, that is when a model close to truly polygenic was used. the proposed methods to estimate additive-by-additive-by-additive of QTL×QTL×QTL interaction effects can be used in association mapping, as a complement and supplement to the analysis of the effects of major genes and epistasis [65][66][67][68]. The total QTLs with additive QTL effects, additive-additive QTL-QTL epistatic effects and additive-additive-additive QTL-QTL-QTL triple interaction effects explained more than 83% of the phenotypic variation. The information obtained in this study will be useful for manipulating the QTLs for plant breeding by marker-assisted selection. Figure 2. Differences between estimated additive-by-additive-by-additive effect and trueassumed value (15) for particular methods of estimation and different number of QTL-QTL-QTL triples. p-phenotypic method, g u -unweighted regression, g w1 -weighted regression: weighted triple interaction effects only, g w2 -weighted regression: weighted all effects: QTL, epistatic and triple interaction.

Data Set
The difference between phenotypic and genotypic estimates of the total additive-byadditive-by-additive effect, observed in both numerical [29] and simulation comparisons, is partly due to the fact that genotypic methods do not find all gene-by-gene-by-gene triplets that determine the trait. However, the use of weighted regression shows that the difference between the phenotypic and genotypic estimation can be reduced, and the total genetic effects will be closer to the true value. Similar results were obtained by estimating QTL-QTL-QTL triple interaction effects for total phenolic content in the mapping population of wheat doubled haploid lines under drought stress conditions [30].
As a direct consequence of interactions, particularly the involvement of QTLs in epistatic [63] and higher-order interactions, the effects of single QTL loci are dependent on the genotypes of other loci. As the results presented here show, consideration of the presence of QTL-QTL epistasis effects had an impact on the evaluation of the effect of QTL-QTL-QTL triple interaction. The lack of epistasis assumed in the simulation studies can be considered in two ways: (1) as the actual absence of epistasis effect or (2) as the summed zero effect of QTL-QTL pairs with epistatic effects different from zero. Attempts to use QTLs in breeding programs must consider not only epistatic effects but also higherorder interactions. Determining the contribution of QTL-QTL-QTL triple interactions is important for understanding the genetic basis of complex traits. Therefore, the genetic models for QTL mapping assuming no QTL-QTL-QTL triple interaction can lead to a biased estimation of QTL parameters and QTL-QTL epistasis [63]. This is indicated by the results obtained by Cyplik et al. [64], where genetic parameters were analyzed for selected traits of maize inbred lines. Although the approach presented in this paper is based on, among other things, a genetic map, the proposed methods to estimate additiveby-additive-by-additive of QTL×QTL×QTL interaction effects can be used in association mapping, as a complement and supplement to the analysis of the effects of major genes and epistasis [65][66][67][68].
The total QTLs with additive QTL effects, additive-additive QTL-QTL epistatic effects and additive-additive-additive QTL-QTL-QTL triple interaction effects explained more than 83% of the phenotypic variation. The information obtained in this study will be useful for manipulating the QTLs for plant breeding by marker-assisted selection.

Data Set
Consider a population of n significantly differentiated homozygous (recombinant inbred-RI or doubled haploid-DH) plant lines from a cross between two homozygous parents, as is currently practiced in the breeding of selfing species. For this population, an n-vector of phenotypic mean observations y = y 1 y 2 · · · y n and q n-vectors of marker genotype observations m l , where l = 1, 2, · · · , q, are obtained. The i-th element (i = 1, 2, · · · , n) of vector m l is equal to -1 or 1, depending on the parental genotype exhibited by the i-th line.

Estimation Based on the Phenotype
The total additive × additive × additive interaction of homozygous loci (three-way interaction) effect on the basis of phenotypic (aaa p ) observations of total phenolic content can be estimated by the formula [29]: where L min and L max are the doubled haploid lines with minimal and maximal mean values, respectively; L is the mean of all homozygous lines. The number of genes (number of effective factors, K) obtained, based on phenotypic observations only, was calculated using the formula presented by Kaczmarek et al. [69]: where V L is additive variance.

Estimation Based on the Genotype
A multiple linear regression model [70,71] of trait values on observations of markers acting as explanatory variables [72] was used to map QTL. It allows the elimination of markers that reveal an association with QTL when tested individually, but in multiple analyses are not characterized by any effect on the phenotypic trait [73]. The estimation of aaa g was based on the assumption that genes responsible for the quantitative trait were completely linked to observed molecular markers. After deciding which p markers out of all observed sufficiently well explain the variability of the quantitative trait, we can model phenotypic observation for homozygous lines as where y-mean values of the observed trait, 1-the n-vector of ones, µ-the general mean, A-(n × p)-matrix of the form A = m l 1 m l 2 · · · m l p , l 1 , l 2 , . . . , l p ∈ {1, 2, . . . , q}, β-the p-vector of unknown parameters of the form β = a l 1 a l 2 · · · a l p , E-matrix which columns are products of some columns of matrix A, γ-the vector of unknown parameters of the form γ = aa l 1 l 2 aa l 1 l 3 · · · aa l p−1 l p , T-matrix which columns are three-way products of some columns of matrix A, δ-the vector of unknown parameters of the form δ = aaa l 1 l 2 l 3 aaa l 1 l 2 l 4 · · · aaa l p−2 l p−1 l p , e-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 quantitative trait, parameters aa l 1 l 2 , aa l 1 l 3 , . . . , aa l p−1 l p , are the additive×additive interaction effects and parameters aaa l 1 l 2 l 3 , aaa l 1 l 2 l 4 , . . . , aaa l p−2 l p−1 l p are the additive × additive × additive interaction effects. We assume that the epistatic and three-way interaction effects show only loci with significant additive gene action effects. This assumption significantly decreases the number of potentially significant effects and causes the regression model more useful.

Unweighted Regression
The selection of markers chosen for model (3) was made by a stepwise trait selection by Akaike Information Criteria [74]. The procedure consisted of two steps: first we divided markers into groups based on the chromosomes they were located on and performed stepwise trait selection by AIC; after that, we combined the remaining markers into one group and repeated selection as above. All of the remaining markers were combined into the final group and the last feature selection was performed on a model with additive × additive × additive interaction effect included. To counteract the multiple comparison problems, we used Bonferroni correction [75].
Denoting by α = µ β γ δ and G = 1 A E T we obtain the model: If G is of full rank, the estimate of α u from traditional (unweighted) multiple linear regression model is given by [76]: The total three-way interaction aaa gu effect of genes influencing the total phenolic content from traditional (unweighted) multiple linear regression model can be found as Weighted Regression The modified version of trait regression on marker data in this paper is considered by taking a weighted multiple linear regression, that is, regression with a diagonal matrix W of unknown variances of observations, which, however, may be empirically found by estimation. In this model, the estimate of α w is: where W = (w ii ) with w ii being the estimated variance for i-th homozygous line, i = 1, 2, . . . , n. Selection of markers for the weighted regression is made in two ways: (i) by the same method as described for the unweighted case (weighted triple interaction effects only) and (ii) using weighted regression to select markers for model (3) (weighted all effects: QTL, epistatic and triple interaction). The total three-way interaction {aaa gw1 -for (i) way and aaa gw2 -for (ii) way} effect of genes influencing the quantitative trait from weighted multiple linear regression model can be found as where x is equal to 1 {for (i) way} or 2 {for (ii) way}.

Coefficient of Determination
The coefficients of determination were used to measure how the model (3) fitted the data and, in this study, was the amount of the phenotypic variance explained by total QTLs with additive effects, QTL-pairs with epistatic effects and QTL-triples with additive-byadditive-by-additive effects (R 2 ). The R 2 was calculated for unweighted regression and in both cases for weighted regression.

Simulation Studies
In the Monte Carlo simulation studies comparing the "phenotypic" and "genotypic" (unweighted and two ways weighted) estimates of the additive-by-additive-by-additive interaction of QTL effects the following variants of the selected parameters were adopted. The true value of the parameter connected with additive × additive × additive interaction effect was assumed to be equal to 15 (aaa = 15) and the total mean value of the quantitative trait to 100. In all 500 homozygous lines were analyzed. Each homozygous line was represented by ten plants. We assumed genotyping for 1000 molecular markers distributed in ten chromosomes of 100 cM. The number of molecular markers per chromosome was equal to 100. The distance between markers was equal to 1 cM. The number of QTLs affecting the trait was 14 (each with an additive effect of two). A QTL was marked perfectly so that the marker and QTLs were in complete linkage disequilibrium with no recombination. We allocated one QTL in chromosome 1 (Q1), two QTLs in chromosomes 3 (Q2 and Q3), 5 (Q4 and Q5), and three QTLs in chromosomes 7 (Q6, Q7 and Q8), 9 (Q9, Q10 and Q11), and 10 (Q12, Q13 and Q14) ( Table 1). We assumed the action of five epistatic pairs: Q1-Q2, Q4-Q5, Q6-Q7, Q8-Q11, and Q1-Q14 ( Figure 1). Thus, we defined two epistatic pairs in the same chromosome: Q4 and Q5 (in chromosome 5) and Q6 and Q7 (in chromosome 7), and three epistatic pairs in distinct chromosomes. Equal epistatic effects were assumed for all five pairs in two variants: two or zero (no epistasis). The number of QTL-QTL-QTL threes with additive-by-additive-by-additive interaction effects that affected the trait was assumed to be one, two, three or five. The QTLs with additive-byadditive-by-additive interaction effects were located in one chromosome (three triples on one chromosome and two triples on two chromosomes), two chromosomes or three chromosomes (each QTL was in a different chromosome) ( Table 1). Effects of particular triples of genes were assumed to be: (i) equal for all triples, or (ii) one QTL-QTL-QTL triple effect was much large than the other (for two triples: 10 and 5; for three triples: 9, 3 and 3; for five triples: 7, 2, 2, 2 and 2). The error variance was equal to 5 of 10. A total of 10,000 data sets containing the vector of phenotypic observations and vectors of marker genotype observations were generated for each combination of the parameters. For each data set the line variances were estimated and the additive-by-additive-by-additive interaction effect was estimated by the phenotypic method aaa jp , and by unweighted aaa jgu and weighted aaa jgw1 , aaa jgw2 (j = 1, 2, . . . , 10, 000) versions of the genotypic method presented above. Additionally, the coefficients of determination R 2 j were estimated. Then, the mean values of parameter estimates aaa p , aaa gu , aaa gw1 and aaa gw2 for each series were calculated, together with the mean squared errors. Mean values of R 2 were calculated. Statistical analysis was performed using GenStat 22nd edition [77].

Conclusions
In practical research, one should expect estimates of the total effect of additive-byadditive-by-additive for triplets of gene-by-gene-by-gene interactions quantitative traits to be smaller than the total phenotype estimate. If a different situation is found, it should be carefully analyzed, checking whether the genetic assumptions of correct segregation and lack of close linkage between markers are met. Weighted multiple linear regression is useful for estimating the parameter associated with additive-by-additive-by-additive of the gene-by-gene-by-gene interaction effect. Regardless of the variant of the approach in the weighted regression, an improvement in the estimation of the parameter under consideration was obtained compared to that obtained using unweighted regression. This indicates the applicability of this estimation method in various experimental situations.
The results indicate that the estimation of additive-by-additive-by-additive interaction effects by the weighted regression method can be applied to different plant species.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.