Our population for GWA mapping consisted of 174 natural inbred lines (i.e., accessions) of A. thaliana (L.) Heyhn. selected from the 1001 Genomes Project (Weigel and Mott, 2009) based on marker information and seed availability (Supplementary Table S1). The laboratory strain Columbia-0 (Col-0) was chosen as the reference. The following lines (in the Col-0 background) were used to isolate T-DNA homozygous mutants of the studied genes: N572850, N586312, N616200, and N620707 (Supplementary Table S2). The ahp1 ahp2 ahp3 (Hutchison et al., 2006) mutants were also used. All lines used were obtained from the Nottingham Arabidopsis Stock Centre (NASC¹). Seeds were stored at 4°C for several weeks (>12) to break dormancy.

Seeds were surface-sterilized in 2% (w/v) NaClO and rinsed with sterile water before being transferred to 120 mm × 120 mm × 10 mm Petri dishes containing 75 mL of one-half-Murashige & Skoog (MS) medium with 2% sucrose, 8 g/L plant agar (Duchefa Biochemie, Netherlands) and 1× Gamborg B5 vitamin mixture (Duchefa Biochemie). After 4 days of stratification at 4°C in darkness, plates were wrapped in aluminum foil and were transferred (0 days after sowing) to an MLR-352-PE growth chamber (Panasonic, Japan) at 22 ± 1°C during 3 days in a nearly vertical position. Plates were unwrapped (3 days after sowing) and grew during another 3 days with continuous light (50 μmol⋅m⁻²⋅s⁻¹). For each accession, 12 seeds were sown per petri dish in triplicate (36 samples/line). Sixteen consecutive sowings including 11 accessions each were established. Additionally, Col-0 was also included in all the sowings to be used as the growth reference accession and for normalization purposes (Supplementary Figure S1A). Lines with germination percentage lower than 80% and with ambiguous marker information were discarded for further analysis (Supplementary Table S1).

To induce LR development during early seedling growth (Van Norman et al., 2014), we excised about 2 mm of the root tip using a sterile scalpel on a laminar flow hood at 6 days after sowing (Figure 1A). Next, samples were transferred back to the growth chamber and followed during 4 days for the analysis of several root traits as described below.

Petri dishes were daily imaged from 6 to 10 days after sowing using an Epson Perfection V330 Photo scanner (Seiko Epson Corporation, Nagano, Japan), at a resolution of 600 dpi, and saved as RGB color images in JPEG file format. Scanned images were processed using EZ-Rhizo (Armengaud et al., 2009) with available plug-in macros to convert them into binary images, remove noise, gap-filling, and skeletonize them prior to automated root detection (Supplementary Figures S1B–S1E). PR length was directly obtained from the EZ-Rhizo output files from scanned images of 6 days after sowing. A highly significant and positive correlation was found for PR length estimated by EZ-Rhizo and directly measured by ImageJ software² from hand-drawn roots (Supplementary Figure S1F). LR number was visually counted from the scanned images between 6 and 10 days after sowing. LR emergence onset corresponds to the day when the first newly emerged LR was visible. LR density was estimated at 8 and 10 days after sowing by the LR number/PR length ratio. Data values were normalized relative to Col-0 values in each sowing dividing each individual value by the Col-0 average (Supplementary Table S3).

Statistical analyses of the data was performed by using StatGraphics Centurion XV (StatPoint Technologies, United States) and SPSS 21.0.0 (SPSS Inc., United States) software packages. Data outliers were identified based on aberrant standard deviation values and were excluded for posterior analyses as described elsewhere (Aguinis et al., 2013). One-sample Kolmogorov–Smirnov tests were performed to analyze the goodness-of-fit between the distribution of the data and a theoretical normal distribution. Non-parametric tests and data transformation were applied when needed. To compare the data for a given variable, we performed multiple testing analyses with the ANOVA F-test or the Fisher’s LSD (Least Significant Differences) methods. Significant differences were collected with 5% level of significance (P-value < 0.05) unless otherwise indicated.

The broad-sense heritability (H²) for the studied dataset was calculated as H² = σ_G²/(σ_G² + σ_GE²/e +σ²_e2/re), in which σ_G², σ_GE², σ_e², r, and e represent the estimated variances for the genetic effects, genotype-environment interactions, random errors, number of replications (12) and number of environments (three), respectively. The estimated variances for σ_G², σ_GE², and σ_e² were obtained by ANOVA using normalized data values as regards to the Col-0 reference accession from 106 of the studied lines.

The population structure of the selected accessions was estimated using the Bayesian model-based clustering algorithm (Porras-Hurtado et al., 2013) implemented in Structure v2.3.4 software³ (Falush et al., 2003). To this end, we used a collection of 319 randomly selected bi-allelic synonymous (likely evolutionary-neutral) SNP markers from available sequence data (Atwell et al., 2010; Cao et al., 2011; Seren et al., 2012). An in-house Matlab script (Supplementary Table S4) was used for single nucleotide polymorphism (SNP) data selection and file formatting. Structure analysis was performed for K = 1 to K = 10 clusters with 20 replicates and 50,000 burn-in period iterations, followed by 50,000 Markov chain Monte Carlo iterations while using a population admixture ancestry model. To determine the most likely number of subpopulations (K) we applied the ΔK method, as described elsewhere (Evanno et al., 2005).

Genome-wide association mapping was performed using the GWAPP web interface⁴, which contains genotypic information for up to 250,000 bi-allelic SNP markers (Seren et al., 2012). GWAS was conducted for the studied traits using the linear regression model (LM) to identify associations between the phenotype of 120 studied accessions and the 205,978 SNPs available in the database. Relative LR numbers were transformed using the y = x function to fit the theoretical normal distribution. Association mapping was obtained excluding from the analyses all SNPs with a frequency <0.12. SNPs with a −log₁₀(P-value) > 6.5 were considered significantly associated to the studied trait (Supplementary Table S5). Manhattan plots, representing the genomic position of each SNP and its association [−log₁₀(P-value)] with the studied trait, were downloaded from the GWAPP web interface. We analyzed the sampling bias on GWAS by systematically removing one or several geographically isolated accessions and found that it did not make any difference to the detected SNP associations (Supplementary Table S6). We selected non-synonymous SNPs with a −log₁₀(P-value) > 6.5 (P-value = 3.16 × 10⁻⁷) for further studies.

Seedlings with T-DNA homozygous insertions in the annotated genes were identified by PCR verification with T-DNA specific primers and a pair of gene-specific primers (Supplementary Table S2). Genomic DNA isolation and genotyping of T-DNA insertion loci PCR were performed as described elsewhere (Pérez-Pérez et al., 2004).

Primers amplified 81–178 bp of the cDNA sequences (Supplementary Table S2). To avoid amplifying genomic DNA, forward and reverse primers bound different exons and hybridized across consecutive exons.

RNA extraction and cDNA synthesis were performed as described elsewhere (Villacorta-Martín et al., 2015). For real-time quantitative PCR, 14 μl reactions were prepared with 7 μl of the SsoAdvanced Universal SYBR Green Supermix (Bio-Rad, United States), 4 μM of specific primer pairs, and 1 μl of cDNA- and DNase-free water (up to 14 μl of total volume reaction). PCR amplifications were carried out in 96-well optical reaction plates on a Step One Plus Real-Time PCR System (Applied Biosystems, United States). Three biological and two technical replicates were performed for each gene. The thermal cycling program started with a step of 10 s at 95°C, followed by 40 cycles (15 s at 95°C and 60 s at 60°C), and the melt curve (from 60 to 95°C, with increments of 0.3°C every 5 s). Dissociation kinetics of the amplified products confirmed their specificity.

Primer validation and gene expression analyses were performed by the absolute quantification method (Lu et al., 2012) by using a standard curve that comprised equal amounts from each cDNA sample. The housekeeping At4g26410 gene (RGS1-HXK1 INTERACTING PROTEIN 1, RHIP1) (Czechowski et al., 2005) was chosen as an internal control and to ensure reproducibility. In each gene, the mean of fold-change values relative to the Col-0 reference genotype was used for graphic representation. Relative expression values were analyzed using SPSS 21.0.0 (SPSS Inc., United States) by applying the Mann–Whitney U-test for statistical differences between cDNA samples (P-value < 0.05).

To validate our experimental approach (Figure 1A), we studied wound-induced LR formation in Columbia-0 (Col-0) during 7 days. PR length remained almost invariable after root tip excision during the whole experiment (Supplementary Figure S1G). The new LRs were already visible at 1 day after PR tip excision (1 dae) and reached 4.97 ± 2.36 (n = 467) LRs at 4 dae (Figure 1B). At 7 dae, the number of LRs slightly increased but it was not possible to measure it unambiguously due to overlap between the LRs of adjacent seedlings (Figure 1B,C). We did not observe a clear spatial pattern of LR emergence from the PRs except that, in all cases, the new LRs emerged from its convex side (Figure 1C, inset). We found a slight variation in PR length and LR number between the different sowings (Supplementary Figures S2A,B), which might be caused by subtle environmental differences at the growth chamber.

We studied wound-induced LR formation in a collection of 173 additional accessions selected from the 1001 Genomes Project (Weigel and Mott, 2009; Supplementary Table S1). In 34 of the studied accessions, the germination percentage at 6 days after sowing was lower than 80% and were discarded for further analysis; other 20 accessions were also discarded because of ambiguous genotypes at the GWAPP web interface (Seren et al., 2012) (Supplementary Table S1). We found variation in all the studied traits (Supplementary Figure S2B). Exceptionally, one or two LRs were observed before root tip excision (0 dae) in some samples, but these were not considered. The broad-sense heritability (H²) was calculated for each of the studied traits (see section “Materials and Methods”). Heritability estimates ranged between 0.90 (LR emergence onset and LR density) and 0.95 (LR number). Broad-sense heritability for PR length was 0.93. Interestingly, we found a positive and significant correlation between LR number and PR length at 4 dae (r = 0.83; Figure 2A), as well as a negative and significant correlation between LR number at 4 dae and LR emergence onset (r = −0.69; Figure 2B). LR number ranged from 1.38 ± 0.82 in Ru3.1-31 (PR length: 0.28 ± 0.07 cm; n = 24) and 9.42 ± 4.84 in Kidr-1 (PR length: 1.90 ± 0.50 cm; n = 24). Hence, reduced LR number in Ru3.1-31 compared to Kidr-1 was likely caused by its reduced PR length. Some accessions, such as Leo-1 and Voeran-1 displayed contrasting phenotypes as regards their LR number (7.07 ± 1.41 and 3.14 ± 1.48 LRs, respectively; n = 29) although they displayed similar PR lengths (Figure 2A,C). On the other hand, Leo-1 and Aitba-2 displayed similar LR number which were larger in Leo-1 likely due to its earlier LR emergence onset (0.55 ± 0.57 days in Leo-1 and 1.82 ± 0.50 days in Aitba-2; n = 29; Figure 2B,C). Ped-0 displayed very short PRs while their wound-induced LRs were much longer (Figure 2C). As regards LR density, Castelfed4.2 and Leo-1 displayed extreme phenotypes, with 2.98 ± 1.35 and 8.05 ± 2.37 roots/cm (n = 29), respectively.

The observed phenotypic distribution for the studied traits (PR length, LR emergence onset, LR number and LR density) suggested that these traits were controlled by multiple genes, that some of the causal alleles are pleiotropic (i.e., affect several of these traits), and that the studied population (n = 120 accessions) was polymorphic for those causal alleles. We determined the genetic relationship among the studied accessions using a Structure analysis with 319 genome-wide randomly selected and synonymous (likely evolutionary neutral) SNP markers already available (see section “Materials and Methods”). Structure analysis of these accessions identified two distinct genetic groups (Figure 3A) that closely correspond to their geographic regions of origin (Figure 3B): The so-called “West” subpopulation including 101 accessions, and the “East” subpopulation with the remaining 19 accessions. However, a detailed analysis of these results indicated a continuous genetic shift from “East” to “West” accessions that follow A. thaliana geographical distribution and that likely arose by local haplotypes, as it has been previously proposed (Platt et al., 2010).

We found a significant variation range for the studied traits between these two genetically distinct subpopulations (Figure 3C). Altogether, accessions belonging to the “East” subpopulation displayed longer PRs and an increasing number of LRs as regards the “West” subpopulation. However, some accessions of the “East” subpopulation, such as Shigu-1, displayed lower phenotypic values for the studied traits than most of their relatives (Figure 4). On the other hand, accessions of the “West” subpopulation and highly genetically divergent from those in the “East” subpopulation (i.e., Aitba-2, HKT2-4, Leo-1, Mrk-0, and Pra-6) displayed higher number of LRs compared with their closest relatives (Figure 4). Despite some population structure among the studied lines and due to high heritability estimates, there is potential for the identification of natural alleles affecting wound-induced LR formation responses through GWA mapping with our dataset.

To obtain insight into the genetic basis of the observed variation in wound-induced responses in young A. thaliana roots, we performed GWA mapping (see section “Materials and Methods”). The significance of association was first evaluated with three different statistical models (LM, KW and AMM; Supplementary Figures S3A–C) and no significant SNP associations were identified by randomization of the phenotypic values within the studied lines (Supplementary Figure S3D). Although LM, and KW usually include more false positives than AMM, they do not present any risk of overcorrection in P-value when applied to traits correlated with population structure (Filiault and Maloof, 2012). We used a conservative threshold of −log₁₀(P-value) > 6.5 and minor allele frequency (MAF) > 12% to select the SNPs being associated with a given trait. A total of 162 SNP associations were found with the LM method for the studied parameters (Supplementary Table S5). We found 32 SNPs associated with PR length with P-values ranging from 1.41 × 10⁻¹⁰ to 3.04 × 10⁻⁷. Thirty-two SNPs were significantly associated with LR emergence onset (P-values ranging from 1.22 × 10⁻⁸ to 3.09 × 10⁻⁷) and only one SNP was found associated with LR density. The larger number of significantly associated SNPs was found for LR number (n = 114), with P-values ranging from 8.27 × 10⁻¹¹ to 3.15 × 10⁻⁷. Consistently with our previous observation that PR length and LR number are significantly correlated, 11 SNPs were significantly associated with both traits; similarly, six significantly associated SNPs were shared between LR emergence onset and LR number (Supplementary Figure S4A).

Next, we classified the selected SNPs based on its molecular effect (Supplementary Figure S4B). About 36% of the significantly associated SNPs were located in intergenic regions and 17.2% of the SNPs laid in the coding region of the annotated gene causing amino acid changes in the protein (Supplementary Figure S4B). Previous reports have shown that, due to linkage disequilibrium, multiple significantly associated SNPs should be found within a small chromosome region for true associations (Rajarammohan et al., 2018). To reduce the number of selected loci for further studies, we focused on 19 candidate genomic regions based on the following criteria (Supplementary Figure S4C): (1) P-value of associated SNPs < 3.16 × 10⁻⁷ (which corresponded to a LOD score > 6.5), (2) presence of multiple significantly associated SNPs within an average of 10 Kpb genomic window, and (3) presence of, at least, one non-synonymous SNP within the selected region. We found five genomic regions putatively contributing to PR length variation in the studied population (Figure 5A), with one (At1g04260), two (At1g04470), three (At2g22660), one (At4g01090), and two (At4g22920, At4g22940) non-synonymous SNPs each (Supplementary Table S5). Three genomic regions were identified as regards their effect on variation in LR emergence onset (Figure 5B), each with one non-synonymous SNP (At1g72250, At1g72300 and At5g20980, respectively; Supplementary Table S5). We found 14 genomic regions putatively involved in the observed variation in LR number among the studied accessions (Figure 5C). Interestingly, the non-synonymous SNPs of three of these regions (dubbed as 2′, 4′, and 5′) overlapped with three genomic regions also selected as being involved in PR length variation (Supplementary Table S5). Hence, the affected genes in these cases, At1g04470, At4g01090, and At4g22940, might indirectly contribute to the LR number differences likely due to their direct effect on PR length before root tip excision. The remaining regions identified in the GWA analysis for LR number might include genes that directly contribute to wound-induced LR formation (Supplementary Table S5) and therefore deserve further studies. On the contrary, no other genomic region fulfilled our selection criteria as regards LR density and this trait was not considered (Figure 5D).

To identify allelic variation in genes contributing to the observed differences in LR number after root tip excision, we selected 20 non-synonymous SNPs for further studies (Supplementary Table S6). Although SNPs in intergenic regions could also be causative, we decided to focus on non-synonymous polymorphisms as their later characterization can be performed on an easier way using reverse genetics tools. Based on both geographic distribution and genotype (Figure 3), we hypothesized that Nemrut-1 and Yeg-1 might represent atypical accessions due to ancient migration and genetic isolation. To discard the false-positive SNPs of this spurious association, we repeated the GWA mapping by removing either one or both accessions, which allowed us to reduce to 11 the number of significantly associated SNPs contributing to LR number (Supplementary Table S6). Due to the genetic structure of the studied accessions, we performed ANOVA analyses as regards the “East” and “West” subpopulations independently. Polymorphisms at eight SNPs affecting six genes (At1g17700, At3g58220, At4g01090, At4g33530, At5g16220, and At5g19710; Table 1) were found significantly associated with LR number variation (Supplementary Table S6). Four haplotypes were detected for selected SNPs within the At4g01090 gene (CAA, CAT, CTA, and TAT). The accessions containing the TAT haplotype displayed a significant increase in LR number, irrespectively of their subpopulation of origin (Supplementary Figure S5A). Indeed, the T to A polymorphism at Chr4:472726 accounted for the quantitative differences in LR number by its own. In addition, we observed a haplotype-dependent relationship between SNPs at At5g16220, and At5g19710, which are separated by 1.4 Mb and were previously assigned to two different candidate genomic regions (Figure 5C). The GA haplotype for these two genes corresponded to the higher phenotypic values for LR number (Supplementary Figure S5B). To our knowledge, this is the first example of two-linked quantitative trait nucleotides (QTNs) detected through GWA mapping and further experiments will account for the functional relationship between the two genes and wound-induced LR number.

To validate the identification of novel genes involved in wound-induced LR formation in A. thaliana seedlings, we chose At4g01090, At4g33530, and At5g19710 for further studies. We searched for available T-DNA insertions in those three genes and identified homozygous mutants by means of PCR and sequencing (Supplementary Table S1). Based on haplotype studies, we found that the Chr4:472726 C/T polymorphism at the At4g01090 gene (Supplementary Figure S6A) was significantly associated with LR number variation, even in those accessions belonging to the same subpopulation such as Fei-0 (5.19 ± 1.88; n = 36) and Star-8 (8.36 ± 2.33; n = 36; Supplementary Figures S6B,C). Additionally, T-DNA homozygotes from the Salk_086312 segregating line displayed a reduced number of wound-induced LR at 4 dae (2.82 ± 1.27; n = 34) in comparison to their wild-type siblings (6.68 ± 1.87; n = 63; Supplementary Figures 6D,E). The homozygous seedlings were also characterized by their longer hypocotyl and shorter PRs (Supplementary Figure 6D). Our results confirmed that the hypothetical protein encoded by the At4g01090 gene participates in wound-induced LR development and that the observed natural variation in their protein sequence might affect their biochemical activity.

At4g33530 (Supplementary Figure S7A) encodes a potassium (K⁺) uptake transporter which is highly expressed in root hairs (Ahn et al., 2004). We found a statistically significant and non-synonymous SNP (Chr4:16128906) correlated with wound-induced LR phenotype variation (Supplementary Figure S7B). The accessions Pu2-7 and Aitba-2 differed in their LR number (6.44 ± 1.95; n = 32 and 9.53 ± 1.94; n = 34, respectively) and carried alternative alleles of the Chr4:16128906 marker (Supplementary Figure S7C). We identified T-DNA homozygotes from two Salk insertion lines interrupting the coding region of this gene (Supplementary Figure S7A). None of the studied homozygous mutants from Salk_120707 and Salk_072850 lines displayed significant differences in wound-induced LR number as regards their wild-type siblings (Supplementary Figures S7D,E).

The At5g19710 gene (Figure 6A) encodes a histidine phosphotransfer protein (AHP) whose function on the CK transduction pathway has not yet been elucidated. There are six other known AHPs involved in CK responses (Hutchison et al., 2006). Phylogenetic tree reconstruction of AHPs including At5g19710 (Figure 6A), suggested that the annotated AHP protein encoded by this gene was incorrectly predicted due to an exon skipping, and clustered together with the AHP4 negative regulator of CK signaling (Figure 6B; Moreira et al., 2013). We found that the Chr5:6665363 G/A polymorphism at the third exon of this gene (Figure 6A) was significantly associated with LR number variation (Figure 6C), even in those accessions belonging to the same subpopulation such as Ll-0 (3.91 ± 2.20; n = 32) and Pra-6 (7.91 ± 1.91; n = 33; Figure 6D). We identified a homozygous T-DNA insertion line for the At5g19710 gene whose seedlings showed a reduced number of wound-induced LRs (Figure 6E,F) due to a significant miss-regulation of At5g19710 gene expression (Figure 6G). We also confirmed that the triple ahp1 ahp2 ahp3 mutants, which was defective in CK root responses (Hutchison et al., 2006), displayed a significant reduction in wound-induced LRs at 4 dae compared to their wild-type background (Figure 6E,F). Taken together, our results seem to indicate that altered homeostasis of AHP proteins required for CK signaling interferes with wound-induced LR formation, a statement that requires further investigation.

Finally, we wondered whether there was an epistatic interaction between the allelic variants for some of the studied non-synonymous SNP markers (Chr4:472726, Chr4:16128906 and Chr5:6665363, respectively) that contributed to the observed variation in wound-induced LR formation in the studied population. Accessions sharing the CCG haplotype for these three markers displayed the smallest number of wound-induced LRs (Supplementary Figure S8A), while only the accessions with the haplotype containing a single polymorphism in the At4g01090 gene showed a significant increase (LSD; P-value < 0.01) in LR numbers (Supplementary Figure S8A). The individual contribution of the SNP polymorphisms in the other two genes considered, At4g33530 and At5g19710, hardly increased wound-induced LR numbers alone but in combination (CGA haplotype) their effects on wound-induced LR formation were enhanced (Supplementary Figure S8A). Similar interactions were found between the other SNP pairs (TGG and TCA haplotypes). Interestingly, we identified several accessions in the “West” subpopulation of all the haplotypes correlated with an increase in wound-induced LR numbers, such as Mrk-0 (8.17 ± 2.12; n = 36) and Vie-0 (6.33 ± 1.80; n = 33) (Supplementary Figure S8B). However, all accessions with the TGA haplotype combining the allelic variants that contribute to an increase of wound-induced LR formation belonged to the “East” subpopulation, which suggested that this trait might be ancestral.