### Mapping population and genotyping

A recombinant inbred line (RIL) population was generated from the cross between the Mesoamerican bean cultivar Jamapa and the Andean cultivar Calima of *Phaseolus vulgaris* L. The RI family comprises 173 lines, which was propagated by single seed descent and in bulk afterwards to the F_{11:14} generation. The nuclear DNA was extracted from leaf of each RI line, and then *PstI* GBS libraries were prepared and submitted for genotyping by sequencing in the Illumina HiSeq platform. Sequence data were processed by multiple bioinformatics tools to obtain SNPs. Linkage map was constructed with 513 unique loci covering 943 cM, which included 442 SNP loci (DiM), 66 RFLP-based markers, three soybean-derived SNP markers, and two phenotypic marker loci for this RILs population. The two parents and RILs were genotyped for 513 molecular markers located on 11 linkage groups each covering a common bean chromosome [40].

### Experimental design and data collection

The mapping population (including both parents) was planted at two sites, Palmira and Popayan in southwestern Colombia during 2011–2012. These two sites have different temperature regimes [45]. In Palmira, the mean air temperature ranges from 19.5 to 28.8 °C; solar radiation, 14.7 MJ m^{− 2} d^{− 1}; day length, 15.6 h; and growing season, from 11 Nov 2011 to Jan 2012. In Popayan, the mean air temperature ranges from 13.7 to 25.5 °C; solar radiation, 11.8 MJ m^{− 2} d^{− 1}; day length, 12.1 h; and growing season, from 23 Mar 2012 to Jun 2012. A randomized complete block row-column design with three replicates (six for each parent) was employed at each site. Each RIL plot had between 30 and 50 plants [46]. One plant was harvested weekly from each RIL plot and from all three replicates. The first five leaves were harvested and measured independently for leaf area and mass (dry weight). The weekly leaf samplings started soon after the plants reached stage V0 (i.e. the most of leaves are fully-expanded) and ended when the plants reached stage R1 (i.e. the time at first anthesis). Data of the first trifoliate leaf at each time point were employed to obtain the mean of three replicates for each RIL for QTL mapping. The combined area of the leaflets of each leaf was measured using a Li-Cor® LI-3100C area meter after harvest. Leaf blades were dried at 65 °C in a drying oven for 3 days. Leaves were equilibrated to room temperature after removing them from the oven before weighing them on a balance with a 10^{− 3} g resolution. We obtained data on leaf area and leaf mass of 173 lines at five time points from each site.

### Statistical modeling

#### Mapping static allometry

Suppose there is a mapping population of *n* recombinant inbred lines (RILs) in which there are two alternative homozygous genotypes at each marker. We are interested in the allometric covariation of leaf area and leaf mass, which can be described by a power eq. (1a) for each time point (time 1, time 2, time 3, time 4, and time 5). Let *y*_{i} and *z*_{i} denote log-transformed values of leaf area and leaf mass measured at a time point for RIL *i* (*i* = 1, …, *n*), respectively. Consider a QTL with two genotypes *QQ* (coded as 1) and *qq* (coded as 2). Let *μ*_{1y} and *μ*_{1z} denote the genotypic values of leaf area and leaf mass for genotype *AA*, and *μ*_{2y} and *μ*_{2z} denote the genotypic value of leaf area and leaf mass for genotype *aa*, respectively. If this QTL affects leaf allometry, we can establish the following relationships from equation (1b):

$$ {\mu}_{1y}=b+\beta {\mu}_{1z} $$

(2a)

$$ {\mu}_{2y}=b+\beta {\mu}_{2z} $$

(2b)

by letting *b* = log*α*. The genetic effects of this QTL on leaf mass and leaf area are calculated as *a*_{z} = *μ*_{1z} − *μ*_{2z} and *a*_{y} = *μ*_{1y} − *μ*_{2y} = *βa*_{z}, respectively. The intercept of leaf area regressed on leaf mass is estimated as \( b=\frac{1}{2}\left[\left({\mu}_{1y}+{\mu}_{2y}\right)-\beta \left({\mu}_{1z}+{\mu}_{2z}\right)\right] \).

Using this QTL’s information, we formulate a bivariate likelihood for two leaf traits based on a mixture model, expressed as

$$ L\left(\varPhi |y,z\right)=\prod \limits_{i=1}^n\left[{\omega}_{1\mid i}{f}_1\left({y}_i,{z}_i\right)+{\omega}_{2\mid i}{f}_2\left({y}_i,{z}_i\right)\right] $$

(3)

where *ω*_{j|i} is the conditional probability of QTL genotype *j* (*j* = 1 for *QQ* or 2 for *qq*), conditional on the marker interval that harbors the QTL, *f*_{j}(*y*_{i},*z*_{i};**Σ**) is the bivariate density function of leaf area and leaf mass for QTL genotype *j*, and *Φ* represents the unknown parameters that describe the location of the QTL and its genetic effects on leaf allometry and residual (co)variances. The conditional probability is expressed in terms of the recombination fraction between the marker interval and QTL [52]. *f*_{j}(*y*_{i},*z*_{i};**Σ**) is assumed to be a normal density function expressed as

$$ {f}_j\left({y}_i,{z}_i\right)=\frac{1}{2\pi {\sigma}_y{\sigma}_z\sqrt{1-{\rho}^2}}\exp \left[-\frac{1}{2\left(1-{\rho}^2\right)}\left[\frac{{\left({y}_i-{\mu}_{jy}\right)}^2}{\sigma_y^2}-2\rho \frac{\left({y}_i-{\mu}_{jy}\right)\left({z}_i-{\mu}_{jz}\right)}{\sigma_y{\sigma}_z}+\frac{{\left({z}_i-{\mu}_{jz}\right)}^2}{\sigma_z^2}\right]\right] $$

where *μ*_{jy} and *μ*_{jz} were defined as above and \( {\sigma}_y^2 \), \( {\sigma}_z^2 \) and *ρ* are the error variances of leaf area and leaf mass and their correlation coefficient, respectively.

The expectation-maximization (EM) algorithm was implemented to estimate the parameters *Φ* = (QTL position; *a*_{z}, *b*, *β*; \( {\sigma}_y^2 \), \( {\sigma}_z^2 \), *ρ*). The significant QTL was estimated by assuming a QTL at every 2 cM position over the linkage map. After the parameters are estimated, we need to perform hypothesis tests. First, whether there exists a significant QTL can be tested by formulating the following hypotheses:

$$ {\mathrm{H}}_0:\cdot {\mu}_{1y}={\mu}_{2y}\cdot \mathrm{and}\cdot {\mu}_{1z}\cdot =\cdot {\mu}_{2z} $$

(4)

H_{1}: At least one of the equalities in the H_{0} does not hold.

The log-likelihood ratio calculated from the H_{0} and H_{1} is compared against the genome-wide critical threshold determined from permutation tests. Second, we can test whether *β* is significantly different from a specific value, e.g. 1 or 3/4. The allometry theory can be used to interpret the biological meaning of these parameters [17].

#### Mapping ontogenetic allometry

Ontogenetic allometry states that leaf area (*A*) is scaled with leaf mass (*M*) over developmental time (*t*). Previous studies noted that as leaves grow, increases in surface area and mass are not synchronous [8]. To accommodate this phenomenon, we introduced an intercept *d* into the power equation [53], obtaining

$$ A(t)={\alpha M}^{\beta }(t)-\mathrm{d} $$

(5)

Let *Y*_{i}(*t*) and *Z*_{i}(*t*) denote the observed values of leaf area and leaf mass at time *t* (*t* = 1, …, *T*) for RIL *i* (*i* = 1, …, *n*), respectively. By statistical reasoning, we found a goodness-of-fit of eq. (5) to *Y*_{i}(*t*) and *Z*_{i}(*t*) data for individual RILs (Additional file 4).

To map how a QTL affects the ontogenetic allometry of leaf area vs. leaf mass, we implemented Zhao et al.’s [54] bivariate functional mapping. Consider a QTL with two genotypes *QQ* and *qq*. Let **Y**_{i} = (*Y*_{i}(1), …, *Y*_{i}(*T*)) and **Z**_{i} = (*Z*_{i}(1), …, *Z*_{i}(*T*)). The bivariate functional mapping is formulated as

$$ L\left(\varPhi |Y,Z\right)=\prod \limits_{i=1}^n\left[{\omega}_{1\mid i}{f}_1\left({Y}_i;{Z}_i\right)+{\omega}_{2\mid i}{f}_2\left({Y}_i;{Z}_i\right)\right] $$

(6)

where *ω*_{j|i} was defined as above, and *f*_{j}(**Y**_{i}; **z = Z**_{i}) is a bivariate longitudinal normal distribution function for QTL genotype *j* (*j* = 1 for *QQ* and 2 for *qq*). The genotype-dependent mean vector of *f*_{j}(**Y**_{i}; **Z**_{i}) is expressed as

$$ \left({\boldsymbol{\mu}}_{jy};{\boldsymbol{\mu}}_{jz}\right)=\left({\mu}_{jy}(1),\dots, {\mu}_{jy}(T);{\mu}_{jz}(1),\dots, {\mu}_{jz}(T)\right) $$

(7)

Based on the allometry law expressed by the power eq. (5), we model the genotypic value of leaf area by

$$ {\mu}_{jy}(t)={\alpha}_j{\upmu}_{jz}^{\beta_j}(t)-{d}_j $$

(8)

Based on the principle of functional mapping, we model the genotype-dependent growth of leaf mass *μ*_{jz}(*t*) by a growth equation, such as logistic equation [55]. Thus, we have

$$ {\mu}_{jz}(t)=\frac{A_j}{1+{B}_j{e}^{- Rt}} $$

(9)

where parameters (*A*_{j}, *B*_{j}, *R*_{j}) are the asymptotic growth, initial growth and relative growth rate of leaf mass for genotype *j* over time, respectively. Therefore, by substituting eqs. (8) and (9) into the mean vector (7), we can model the genotypic values of each QTL genotype for leaf area and leaf mass at different time points through parameters (*A*_{j}, *B*_{j}, *R*_{j}; *α*_{j}, *β*_{j}, *d*_{j}). Such modeling has incorporated the allometric scaling law and growth law, which makes QTL mapping of biological relevance and robustness.

The (co)variance matrix (**Σ**) of *f*_{j}(**Y**_{i}; **Z**_{i}) is a symmetric matrix containing three longitudinal (co)variance matrices within and between leaf area and leaf mass. Zhao et al. [54] derived a bivariate first-order structured antedependence (bi-SAD(1)) model to fit the structure of **Σ**. The bi-SAD (1) model requires the following parameters: antedependence parameters (*ρ*_{y} and *ρ*_{z}) and innovation variances (\( {\sigma}_y^2 \) and \( {\sigma}_z^2\Big) \) for leaf area and leaf mass, and correlation between the two traits (*ρ*_{yz}). Based on these parameters, Zhao et al. [54] derived the closed forms of the determinants and inverse of the bi-SAD(1)-structured longitudinal matrix. By implementing such closed forms into the likelihood (6), we can greatly increase the computational efficiency and precision of parameter estimation *Φ* = (QTL position; *A*_{j}, *B*_{j}, *R*_{j}; *α*_{j}, *β*_{j}, *d*_{j}; *ρ*_{y}, *ρ*_{z}, *ρ*_{yz}, \( {\sigma}_y^2, \)\( {\sigma}_z^2\Big) \).

Accordingly, the estimation of unknown parameters can be made by a hybrid between the EM and simplex algorithms. The QTL position can be estimated by a grid approach. The presence and location of an ontogenetic-allometry QTL can be tested through the null hypotheses:

$$ {\mathrm{H}}_0:\left({A}_j,{B}_j,{R}_j;{\alpha}_j,{\beta}_j,{d}_j\right)\equiv \left(A,B,R;\alpha, \beta, d\right)\ \mathrm{for}\ \mathrm{all}\ j=1,2 $$

(10)

These tests are carried out by calculating the log-likelihood ratio for the QTL effect at each position and comparing it to the genome-wide critical threshold determined from permutations tests. We can specifically test whether the QTL affects the growth trajectory of leaf mass (11), and/or its ontogenetic allometry with leaf area (12).

$$ {\mathrm{H}}_0:\left({A}_j,{B}_j,{R}_j\right)\equiv \left(A,B,R\right)\ \mathrm{for}\ \mathrm{all}\ j=1,2 $$

(11)

$$ {\mathrm{H}}_0:\left({\alpha}_j,{\beta}_j,{d}_j\right)\equiv \left(\alpha, \beta, d\right)\ \mathrm{for}\ \mathrm{all}\ j=1,2 $$

(12)

If both null hypotheses are rejected, then this means that the QTL pleiotropically affects leaf area and leaf mass growth. Possible functions of all the QTLs were annotated via BLAST in the “nr” database on website of National Center of Biotechnology Information (NCBI; http://blast.ncbi.nlm.nih.gov/).