article: Vertical growth of plants is a dynamic process that is influenced by genetic and environmental factors and has a pronounced effect on overall plant architecture and biomass composition. We have performed six controlled growth trials of an interspecific Setaria italica x Setaria viridis recombinant inbred line population to assess how the genetic architecture of plant height is influenced by developmental queues, water availability and planting density. The non-destructive nature of plant height measurements has enabled us to monitor height throughout the plant life cycle in both field and controlled environments. We find that plant height is reduced under water limitation and high density planting and affected by growth environment (field vs. growth chamber). The results support a model where plant height is a heritable, polygenic trait and that the major genetic loci that influence plant height function independent of growth environment. The identity and contribution of loci that influence height changes dynamically throughout development and the reduction of growth observed in water limited environments is a consequence of delayed progression through the genetic program which establishes plant height in Setaria. In this population, alleles inherited from the weedy S. viridis parent act to increase plant height early, whereas a larger number of small effect alleles inherited from the domesticated S. italica parent collectively act to increase plant height later in development.
  Organismal growth is a complex process that is influenced by genetic and environmental factors in a time dependent manner. Understanding the dynamics of these interactions is a critical step in understanding how traits are determined. Height is one of the most influential components of plant architecture. Due to ease of measurement, high heritability and agronomic importance, plant height has been an attractive topic for scientific inquiry for over a century [1]. Plant investment in vertical growth has adaptive value because it helps a plant to out-compete neighbors for access to solar radiation. However, vertical growth requires a concomitant energetic allocation for construction of stem tissue at the expense of allocation to lateral growth and other processes. Greater height growth also amplifies hydraulic cost [2], and increases the risk of lodging [3]. The benefit to cost ratio of such adaptation is likely directly influenced by growth environment [4]. Adaptations introgressed into wheat and rice varieties during the Green Revolution that led to reduced height, increased yield, harvest uniformity, improved carbon partitioning and nutrient and water use efficiency. In biomass crops, height has been shown to be positively correlated with above ground biomass [5,6]. Thus, depending on the crop and breeding objective either increased or decreased height may be targeted. Plant height in the Poaceae, a family of grasses that includes the cereal crops and bioenergy grasses, is a function of internode length and number, the increase of which is terminated at reproductive maturity. Plant height within the grasses is a highly heritable, polygenic trait [7–14]. Genetic studies of height in maize, sorghum, sugarcane, wheat, barley, setaria and rice have identified well over 100 QTL [7–16], and enabled the cloning of multiple loci that contribute large effects [9,15–36]. Forward genetic screens of mutant populations has also been an effective approach to identify causative genes associated with height [24,37]. Although plant height exhibits a strong degree of genetic determinism, it is influenced by environmental factors such as water availability and planting density and exhibits dynamic behavior throughout the plant life cycle [10,38–42]. The strong influence of environment on plant height suggests that multiple signals impinge on the regulation of plant height throughout development. Surprisingly, relatively few studies have assessed how the temporal genetic architecture of this dynamic trait changes through developmental time [14,43]. A major challenge of performing such experiments is the difficulties associated with growing large populations in conditions with contrasting environmental variables while capturing high precision measurements of phenotypes at appropriate intervals. Recent use of modern high-throughput phenotyping (phenomics) technology is beginning to alleviate these technical obstacles [44,45]. Utilizing genetic model systems that possess desirable growth attributes and tractable genetics [46–48] in combination with phenomics technology improves our ability to obtain a more holistic systems view of growth through developmental time and in the face of environmental challenge. As a model system, plants in the genus Setaria, possesses many favorable experimental and life history attributes that make temporal genetic analysis of complex traits possible and the genetic dissection of these traits more feasible than in larger evolutionarily related crop plants like maize and sorghum [46,47]. In this study we conduct a large scale and multi-environment analysis of plant height. We identify thirty seven QTL associated with plant height throughout development within a Setaria recombinant inbred population (RIL) population [49–51] across twelve different experiments varying water availability and planting density. Analysis of ortholog positions suggests that known height affecting genes from the Poaceae account for only a fraction of the loci identified. The results of this study are discussed in light of development and environmental variation and suggest a strategy for the fine scale manipulation of plant height in the Poaceae.
  A Setaria F7 RIL population comprised of 217 individuals was used for genetic mapping. The RIL population was generated through an interspecific cross between the wild-type S. viridis accession, A10, and the domesticated S. italica accession, B100 [49–51]. Plant height measurements of individuals within this population were collected in six different trials. Four of the six trials were conducted at field sites located at University of Illinois-Champaign. Two field trials designed to assess plant phenotypic response to planting density and water availability were each conducted in both 2013 and 2014. The controlled environment trials were performed in a growth chamber at Carnegie Institute for Science, Stanford, CA and using the Bellweather Phenotyping Facility at Donald Danforth Plant Science Center [52]. Field experiments were conducted on the South Farms at the University of Illinois Urbana-Champaign in summers 2013 and 2014. The field site is rain-fed, tile-drained, has a deep and organically rich Flanagan/Drummer series type soil [53]. For each field experiment, seeds were first germinated in a greenhouse in plug trays (128sq, T. O. Plastics) containing a media mix composed of sphagnum, vermiculite, fine bark, bark ash, gypsum, slow-release nitrogen, dolomitic limestone, and a wetting agent (Metro-Mix 360, Sun Gro Horticulture). Plugs were provided with ample water daily by clear water or by fertigation every other day (EXCEL CAL MAG 15-5-5). Seedlings were hand transplanted into genotype specific subplots (dimensions 7 x 6) within a mechanically tilled field approximately 1 week after sewing. Control over water availability was achieved by growing plants under retractable awning rainout shelters at 5 cm spacing to block precipitation from experimental plots and re-applying differential drip irrigation regimes to wet and dry replicates. Each awning is referred to as a plot and the genotype specific blocks within each awning are called subplots. Planting density experiments were performed by growing each genotype specific subplot at 5 cm (dense) and 25 cm (sparse) spacing. Plants (F7 RIL population) were grown in deep pots (Stuewe & Sons, Oregon, D25L) with one plant per pot (one plant per RIL was phenotyped). Seeds were directly germinated in 14 inch Deepots filled with a soil mixture containing 75% Pro-Mix (PRO-MIX (r) PGX soil, Premier Tech, Canada) and 25% river sand, imbibed in water to pot capacity. Pot capacity, defined as the amount of water that soil in a pot can hold against the pull of gravity, was estimated to be 230 ml. Plants were grown in a growth chamber (12 hours light at 31°C and 12 hours dark at 23°C with constant relative humidity at 54–55%). Well-watered plants were watered back to field capacity once every third day whereas for water deficit experiments, seeds were sown in soil imbibed to pot capacity and no further water was added. To prevent water loss, the bottom of each pot was covered with a plastic bag. After a six week stratification in moist long fiber sphagnum moss (Luster Leaf Products Inc. , USA) at 4 C, Setaria seeds were planted in four inch diameter (10 cm) white pots pre-filled with ~470 cm3 of Metro-Mix 360 soil (Hummert, USA) and 0. 5 g of Osmocote Classic 14-14-14 fertilizer (Everris, USA). After planting, seeds were given 7 days to germinate in a Conviron growth chamber with long day photoperiod (16 h day/8 h night; light intensity 230 μmol/m2/s) at 31°C day/21°C night before being loaded onto the Bellweather Phenotyping System using a random block design. Plants were grown on the system for 25 days under long day photoperiod (16 h day/8 h night; light intensity 500 μmol/m2/s) cycling the temperature to match the photoperiod (31°C day/21°C night) and maintaining the relative humidity between 40–80%. Weighing and watering of plants was performed between 2–3 times per day to maintain soil volumetric water content at either 40% full-capacity (FC) (drought) or 100% FC (well-watered) as determined by [52]. Water limitation began 15 days after planting. Individual plants were imaged at 4 different angular rotations (0°, 90° 180°, 270°) every other day. The effect of pot size was evaluated by growing parental lines of Setaria viridis (A10) and Setaria italica (B100) in four different size pots: 160 cm3 (T. O. Plastics, 715348C), 473 cm3 (Pöppelmann, VCC 10 F US), 950 cm3 (Stuewe & Sons, MT38) and 6033 cm3 (Nursery Supplies, C600) filled with Metro-Mix 360 soil (Sun Gro Horticulture). Seeds were given 7 days to germinate in a greenhouse and watered daily by the Donald Danforth Plant Science Center greenhouse staff. Plants were grown for 50 days under long day photoperiod (14 h day/ 10 h night) cycling the temperature to match photoperiod (27°C day/ 23°C night) and maintaining the relative humidity between 40–100%. In the field, plant height was measured as the distance from the plant base to the uppermost leaf collar on the culm. In 2013 this was measured by judging the average height of the canopy in a subplot relative to a yard stick and recorded the measurement using a Janam XP20 hand held scanner [Janam Technologies; Woodbury, NY]. In 2014 height was directly measured repeatedly on three tagged plants in each subplot. In each field experiment, a direct culm height measurement was taken at the final destructive biomass harvest. In 2013 this was done on three representative plants from each subplot; in 2014 this was done on the same-tagged plants as the in-field repeated height measurements. In the 2013 Density experiment, height was measured at 19,25,31,38,45,52,59, and 67 days after seed sowing. Final biomass harvest began 87 days after seed sowing. In the 2013 Drought experiment, height was measured at 21,29,36,43,50,59, and 67 days after seed sowing. Final biomass harvest began 92 days after seed sowing. In the 2014 Drought experiment, height was measured at 25,33,40, and 47 days after seed sowing. Final biomass harvest began 59 days after seed sowing. In the 2014 Density experiment, height was measured at 25 and 46 days after seed sowing. Final biomass harvest began 67 days after seed sowing. Panicle emergence was measured as the number of days after sowing at which the panicle head was seen past the collar of the culm flag leaf in at least half of the individuals in a genotype specific subplot. At Carnegie, plant height was measured from the base of the plant to the tip of the panicle (at 48 days after seed sowing). Five different measurement functions encoded within the PlantCV software package were used to estimate plant height in the images collected at the Bellweather Phenotyping Facility [52] (S1 Fig). These include length of the plant object along the y-axis (extent_y), length of the plant object from the top of the pot to the maximal y-axis point identified as the plant object (height_above_bound), the distance from the top of the pot to the y-coordinate reported as the center of mass of the plant object (centroid_y), the distance from the top of the pot to the y-coordinate identified as the center of the ellipse that encompasses the plant object (ellipse_y) and an additional function designed to estimate the distance between the top of the pot and y-coordinate where average plant width is maximized (canopy_height). Canopy height is calculated by dividing the height_above_bound measure into 20 different equally sized scoring windows (divide the object into 5% intervals along the y-axis, scoring is based upon median width within each interval) and measuring the median plant width within each window. The function reports the average of y-coordinates that comprises the scoring window and median width within the scoring window where the object (plant) is the widest (canopy_width). Centroid_y and ellipse_center_y were the most and least heritable measurements of plant height respectively (S2a Fig and S2 Table). Height_above_bound and extent_y exhibited the largest proportion of variance attributed to water limitation (S2b Fig and S3 Table). Summarization of height measurements derived from four rotational side view images was performed by calculating the mean, median, maximum or minimum of each distribution. Summarization method does not appear to influence the heritability of this these traits (S3 Fig). All height metrics are reported as the average value estimated from all individual rotational images at each plant time point. Lengths of objects measured in pixels were converted to centimeters (cm) using the calibration reported by [52]. Plant height of parental lines grown in the greenhouse at the Donald Danforth Plant Science Center were measured as the distance from the base of the plant to most distant leaf collar on the culm using a tape measure. Measurements were recorded every other day during the workweek for a span of 50 days. A de novo genetic map of the A10 x B100 Setaria RIL population was constructed through genotyping by sequencing (GBS) using the methods described in [54,55]. Briefly, leaf tissue was collected from plants in the 2014 planting density experiment in addition to greenhouse trials. High molecular weight DNA was extracted using the CTAB method [56] and quantified using a Nanodrop 1000 Spectrophotometer [Thermo Scientific; Waltham, MA]. DNA was then double-digested using the restriction enzymes Pst1 and MspI and a unique barcode and common adaptors were ligated onto their respective fragment ends [57]. After size selection and cleanup using Agencourt AMPure XP beads [Beckman Coulter; Brea, Ca], barcoded samples were pooled and submitted for sequencing (100-bp single end reads) using two Illumina Hi-Seq lanes [Illumina; San Diego, CA]. The resulting sequence reads were processed using the Tassel 3. 0 GBS pipeline [58]. Raw sequence reads were filtered to remove any tags that were not present 20 times at minimum within in each lane of sequencing and anchored to the S. viridis genome version 1. 1 (Phytozome) using bowtie2 [59] with default parameters. Anchored SNP polymorphisms exhibiting an inbreeding coefficient less than 0. 9 were removed and identical taxa within and between library preps (2 sequencing lanes) were called as heterozygous if the ratio between alleles in and between all identical taxa and/or library preps was less than 0. 8. All polymorphisms categorized as heterozygous were called as missing values for the remainder of the analysis, yielding 17,612 SNPs. The genetic map was generated using the R/qtl software [60]. Genotype calls within the combined A10 taxa were compared to the S. viridis reference genome version 1. 1 (Phytozome) and all missing and non-congruent polymorphisms were removed. The remaining sites were filtered to remove SNPs not present in > 90% of the unified taxa yielding 217 RILs and 3,295 markers. All polymorphisms were then recoded to reflect their parent of origin. Dubious genotype calls were identified by examining the recombination frequency between adjacent markers. Based upon manual inspection of histograms generated of recombination fraction between adjacent markers, all markers exhibiting a recombination frequency between adjacent markers greater than 36% (80 recombinations / 217 individuals) were considered to be erroneous and removed from the map. Next imputation of missing genotypes was performed at all sites that did not exhibit recombination between missing makers (using the R/qtl fill. geno () function and specifying the ‘no_dbl_XO’ method) and duplicate markers were removed from the genetic map. Finally, the genetic distance between markers was calculated using the Kosambi mapping function assuming a genotype error rate of 0. 005. Prior to QTL mapping, individuals missing more than 10% of their genotype calls are removed yielding a genetic map containing 208 individuals genotyped at 1,595 sites with 99. 6% coverage. Based on pairwise similarity between identical individuals, the GBS map constructed here is highly consistent with other genetic maps of this population described in the literature [49,61]. Variance components corresponding to broad sense heritability and total variance explained were estimated using a mixed linear model using the R package lme4 [62]. Broad sense heritability was calculated using two methods. Within an individual experiment, broad sense heritability on a line-estimate basis was calculated using the following formula: H2experiment=σ2genotype/ (σ2genotype + (σ2genotype X treatment/ ntreatment) + (σ2residual/ nreplicates) ) in which ntreatment is the harmonic mean of the number of treatment blocks in which each line was observed and nreplicates is the harmonic mean of number of replicates of each genotype in the experiment. Heritability within treatment blocks was calculated by fitting a linear model with genotype as the only explanatory factor within each treatment block. The proportion of variance attributed to genotype divided by total variance within each treatment block is reported as broad sense heritability within treatment. Total variance explained was calculated by fitting a linear model including factors genotype, treatment, plot and genotype x treatment effects across all phenotypic values in all treatments. The proportion of variance that is incorporated into these factors divided by the total variance in the experiment is reported as total variance explained. Best linear unbiased predictors (BLUPs) of plant height for each genotype were predicted from a linear mixed effect model independently within all individual time points of field drought and density experiments in 2013 and 2014 using lme4 [62]. The following terms were included in the model: Yijk=μ+ genotypei+ treatmentj+ plot (treatment) k (j) + (genotype x treatment) ij+εijk In which Yijk is the individual observation of plant height, μ is the overall mean, genotypei is the effect of the ith genotype, treatmentj is the effect of the jth treatment, plot (treatment) k (j) is the effect of the kth plot nested within the jth treatment, and (genotype x treatment) ij is a term describing the interaction between the ith genotype and jth treatment and εijk is the random error term. As fixed effects we entered treatment and plot nested within treatment into the model. The terms genotype and genotype by treatment interactions were modeled as random effects. Logistic regression was used to model the line values of plant height in experiments performed at the Bellweather Phenotyping Facility. Three types of logistic regression models considered (three-parameter logistic, four-parameter logistic and Gompertz) and estimates from the model with the lowest Akaike information criterion (AIC) score were used. Modeling was performed using R functions described in [63]. QTL mapping was performed at each time point within treatment blocks and on the numerical difference, relative difference and trait ratio calculated between treatment blocks using functions encoded within the R/qtl and funqtl package [60,64]. The functions were called by a set of custom python and R scripts (https: //github. com/maxjfeldman/foxy_qtl_pipeline). Two complimentary analysis methods were utilized. First, a single QTL model genome scan using Haley-Knott regression was performed to identify QTL exhibiting LOD score peaks greater than a permutation based significance threshold (α = 0. 05, n = 1000). Next, a stepwise forward/backward selection procedure was used to identify an additive, multiple QTL model based upon maximization of penalized LOD score. Both procedures were performed at each time point, within treatment blocks and on the numerical difference relative difference and trait ratio calculated between phenotypic values measured in treatment blocks at each time point. QTL associated with difference or ratio composite traits may identify loci associated with genotype by environment interaction [65]. QTLs across different experiments were considered to be the same if they were in a 20cM window around the most commonly identified marker in that window. Tests for epistatic interaction between QTL were performed using the scantwo function in the R/qtl package [60]. Epistasis between QTL was evaluated by comparing the log10 likelihood ratio of a model describing the additive and interaction between two QTL with a model only containing the additive interaction between two QTL. The significance threshold was determined through permutation (α = 0. 05, n = 100). The function-valued approach described by Kwak et al. , 2016 [64] was used to identify QTL associated with the average (SLOD) and maximum (MLOD) score at each locus throughout the experiment [64]. The genotypic mean height within treatments was estimated using a logistic function and the QTL significance threshold was determined based upon permutation-based likelihood of observing the empirical SLOD or MLOD test statistic. Across all grow outs, separate; independent linkage mapping analysis performed at each time point identified a larger number of QTL locations relative to similar function valued analysis based on the SLOD and MLOD statistics calculated at each individual marker throughout the experimental time course. After refinement of QTL position estimates, the significance of fit for the full multiple QTL model was assessed using type III analysis of variance (ANOVA). The contribution of individual loci was assessed using drop-one-term, type III ANOVA. The proportion of variance explained by in addition to the allelic effect size of each locus were determined by comparing the fit of the full model to a sub-model with one of the terms removed. All putative protein coding genes (Setaria viridis genome version 1. 1) found within a 1. 5-logarithm of the odds (LOD) confidence interval were reported for each QTL. A search for orthologous genes associated with plant height (putative causal genes) was performed within the 37 unique QTL intervals identified throughout the genome. A curated list of plant height genes cloned from maize, sorghum and rice was generated through literature review and searching Maize GDB with search terms including “height” and “dwarf”. The coding sequence name and genomic coordinates of each of these genes was downloaded from their species-specific repository (Maize GDB, Rice Genome Annotation Project, Phytozome). Regions of genomic synteny between Setaria viridis and maize, sorghum or rice were identified using the synmap function encoded by the CoGe web tool (default parameters) [66]. The genomic coordinates of coding sequences of interest were extracted from the CoGe text output using a custom python script (https: //github. com/maxjfeldman/Feldman_Setaria_Height_2016). Reciprocal BLAST searches were used to identify the location of orthologous genes not found in syntenic blocks. The location of these sequences was plotted across the genome to examine their location relative to the coordinates of the QTL and median confidence interval at each location.
  Height in the field was measured using two different methods of approximation. An analysis of variance indicates that on average, measuring the exact height to the collar of the uppermost leaf on the culm of a subset of three individuals within a plot (as performed in 2014) was approximately 5. 8% more heritable than attempting to estimate the average height to the collar of the uppermost leaf on the culm of all plants within the plot (as performed in 2013) (S1 Table). An increase in heritability between years was not observed for flowering time between the 2013 and 2014 field seasons, suggesting that the increase in heritability was not due to differences in environment or management (S1 Table). Five different automated methods of measuring plant height (height_above_bound, extent_y, centroid_y, ellipse_center_y, canopy_height—see definitions in methods) from 2-D images captured at the Bellweather Phenotyping Facility were evaluated (S1 Fig). Among these measurements, height_above_bound, extent_y and centroid_y were highly correlated (R2 > 0. 97) (S2 Table) and all three were highly heritable on average (h2 > 0. 85) and showed clear effects of water limitation (S2 Fig, S3 and S4 Tables). Summary statistic does not greatly influence the heritability of height type measurements (S3 Fig). The mean value of height_above_bound based upon four rotational images was the primary metric used for the remainder of this study as it is intuitively simple to understand, highly heritable, strongly correlated with other height metrics, sensitive to treatment effects and previously validated against manual measurements [52]. Height of individuals within the S. viridis (A10) X S. italica (B100) RIL population was evaluated under very diverse growing conditions. Variation resulted from different water availability and planting density treatments. In addition, day length and temperature varied between plant-outs at three locations (Illinois field site, DDPSC growth chamber, Carnegie growth chamber) and two different planting dates each year at the Illinois field site. Although the mean and variance of final height differed significantly between grow outs (Fig 1a), the coefficient of variation remained relatively constant across all replicated experiments and treatments (Table 1 and S5 Table). The proportion of variance attributed to experimental factors (genotype versus treatment versus error) also varied between grow outs (Fig 1b) and throughout time in a given experiment (Fig 1c and S4 Table). Across all experiments the mapping population size averaged ~160 individuals with a minimum count of 146 (Table 1). Final height estimated near the end of each these experiments averaged 38 cm but every experiment exhibited at least a 20 cm difference between the 5th and 95th percentile (Table 1). Modest transgressive segregation was apparent across all experimental grow outs (Fig 1a). At the end of all eight field grow outs, the B100 parental line was taller than the A10 parent, whereas in three out of four grow outs conducted in artificial environments at Carnegie and Bellweather, A10 was taller (Fig 1a). Within each trial and treatment, the heritability of final plant height exceeded 0. 84 (S1 Table). The mean heritability of plant height in this population across all time points was 0. 82 and tended to be higher in controlled environments (average H2 = 0. 86 at Bellweather and Carnegie, average H2 = 0. 83 in field at Illinois). Generally, the proportion of height variance attributed to genotype increased as the plant approached maturity (Fig 1c and S6 Table). The effects of water limitation and planting density treatments did not dramatically influence the heritability of plant height in this study (S1 Table). The block structure of the field experimental design and the frequent repeated measurements on the Bellweather phenotyping system enabled modeling of the underlying trait phenotypes. Best linear unbiased estimated line values of height within individual experiment and time points were used to account for variation attributed to plot and treatment effects in field studies. Logistic regression was used estimate the expected value of repeated temporal measurements of height collected at the Bellweather Phenotyping Facility. Water limitation experiments were conducted in all three locations; each of which differed by genotypes planted, statistical power and the severity of treatment administered in addition to other latent variables inherent to growth location and experimental design. The effect of these factors was assessed at a single time point representative of final plant size using type III analysis of variance [67]. Within drought trials, genotype had the strongest effect on height, followed by trial location and treatment effect nested within location (S7 Table). Similar to the correlation of height and days to anthesis in maize [13], a significant positive correlation between final height and day of panicle emergence was observed in each field grow out (Table 1), suggesting that different planting dates, with their different day lengths and temperatures could explain some of the trial location variation. The water limitation treatment varied widely between experiments, simulating several possible scenarios plants could experience. Field trials at Illinois and the controlled growth trial at Carnegie were designed as progressive drought experiments, where after an initial growth period, water is withheld. In the Bellweather phenotyping experiment, a reduced water availability set point was maintained after an initial dry down. The controlled growth experiments all showed a significant effect of drought, as did the 2014 field experiment (despite an end of season flood). In retrospect, water was not withheld early enough in the 2013 field experiment to detect a significant drought effect on height although soil probes showed a substantial difference in water availability within 3 weeks (S8 Table). The effect of withholding water completely soon after germination in the Carnegie experiment made the observed height effect the most pronounced of any grow out observed (S5 Fig). Height appears less responsive to planting density than water limitation, with most of the variance attributed to genotypic and plot nested within treatment followed by the variance attributed to year and lastly planting density (S7 Table). Nonetheless, in both 2013 and 2014, planting density had a significant effect on height by the end of the experiment. Using the rank order of lines by phenotype as a proxy for how similar the experiments are to each other, the effect of experimental location is apparent, as well as the difference between controlled environment and the field. (Fig 2 and S4 Fig). The well-watered treatment at Bellweather, the drought treatment at Carnegie and the 2013 field density experiment are the most divergent environments encountered within this study (Fig 2). Two previous genetic maps have been constructed for this population [49,61] but neither included all of the lines used in this study. We constructed a de novo genetic map and genotyped all phenotyped RILs to verify the identity of genetic stocks [68]. DNA samples were prepared from 183 field grown and 34 greenhouse grown RIL lines (in addition to the two parental genotypes) and GBS libraries were prepared and sequenced [54,55]. After alignment to the recently completed S. viridis genome (Phytozome v1. 1) using bowtie2 [59] and filtering to select SNPs that match the sequenced parent samples and removal of duplicates, the final genetic map consists of 217 individuals, spanning 1,088 cM, calculated using 1,595 markers (S9 Table). After imputation, the average spacing between markers is 0. 7 cM with the greatest distance between markers being 14. 4 cM (S9 Table). The map size, coverage and identity of individuals in this RIL population are largely identical between map builds [49,61]. We implemented a standard stepwise forward/backward selection procedure to identify the optimal, additive, multiple QTL model based upon penalized LOD score [60]. Mapping was performed on all traits, at each time point within treatment blocks and on the numerical difference, relative difference and trait ratio between phenotypic values by treatment blocks. This procedure identified a total of 153 individual SNPs associated with height across all time points and treatments within the 12 experimental grow outs (Fig 3 and S10 Table). No significant epistatic interaction between QTL was detected. Many of the individually identified SNP positions group into clusters of linked loci that are likely representative of a single QTL location. Groups of closely linked QTL markers (10 cM radius) were collapsed into the most frequently observed SNP marker in each cluster, resulting in 37 unique QTLs (S11 Table and S6 Fig). A substantial proportion of these QTL were observed at least once in three (17/37) or four (12/37) treatment groups, suggesting that the genetic architecture of height is similar across all experimental treatments contrasts considered in this study (S11 Table). Confirming this observation, Mauro-Herrera et al, with the same population grown in the greenhouse and field and height measured at three different timepoints, identified 7 of the 12 QTLs [16]. A static QTL model including only the 12 QTL locations identified at least once in all treatment groups was fit at a time point representative of final plant size in each experiment. Across all experimental grow outs these 12 QTL explained between 25–64% of the additive genetic variance and 53% on average (Fig 4). Across all experimental time points in this study the three loci that explain the greatest proportion of variance in any single time point are included within this set of 12 QTL. The most frequently observed QTL (marker: S5_41999990) was identified in each of the 12 experiments and explained the largest percentage of additive genetic variance (mean of 19%). The B100 allele at this QTL position is associated with reduced plant height. The individual contribution of the other 11 QTL at this representative time point on average, did not explain more than 7% of the additive genetic variance alone, but cumulatively these 11 QTL account for greater than 33%. The B100 allele at nine of these twelve positions are associated with increased plant height. Across all field trials performed, the number of QTL shared between treatment blocks and across years is greater than the number of QTL unique to any treatment block across growth years or both treatment blocks within a single growth year (Fig 5 and S7 Fig). No overlap between QTL detected exclusively in well-watered or water-limited conditions was observed between any of the experiments (Fig 5a and S7 Fig). Within planting density experiments a single QTL (marker: S9_52254840) was identified uniquely within the high planting density treatment block in both years (Fig 5b and S7 Fig). A majority (13/16) of the QTL detected in controlled environments were also detected in field environments (Fig 5c). The results from QTL scans performed on the numerical difference, relative difference and the ratio of plant heights measured between treatment blocks at all time points largely overlapped with the results from the QTL analysis performed on measured value of plant height within treatments (S8 Fig and S12 Table). In total, 18 of the 29 unique QTL associated with these numerical metrics of difference in plant height between treatment blocks were also identified during our QTL analysis of plant height (S8 Fig and S12 Table). In contrast to the QTL locations identified for plant height, most of which were found in multiple experiments, of the 11 difference QTLs that do not overlap 10 were only found in one experiment (S8 Fig and S12 Table). While there are genetic loci that appear to be responsive to different environments, the majority of the environmentally responsive loci are due to relative differences in the timing of the loci growth affect between treatments. Height is a dynamic trait that changes at different rates throughout the plant life cycle. To determine which growth stage the identified loci influence, we examined the best fit multiple-QTL model produced at each separate independent time point and performed an analysis using the average LOD (SLOD) and maximum LOD (MLOD) function-valued approaches as described by [64,69]. Across all grow outs, linkage mapping analysis performed independently at each time point identified a larger number of QTL locations relative to the function-valued analysis. The positional locations of major QTL detected using each of these methods varied slightly, but in all cases it is likely that the same QTL are identified. To simplify discussion, we will refer to these QTL regions using notation that describes the approximate location of the QTL on each chromosome (QTL detected on chromosome 5 at position 99. 1 cM and at position 105. 2 will be reported as 5@100 for example). The relative significance of each QTL location changes substantially throughout the life cycle of plant (Fig 6 and S9 Fig). In all experimental grow outs, plant height at early stages of development is strongly associated with a single major QTL located at 5@100. In the dataset with the highest temporal resolution (Bellweather), the proportional contribution of this QTL diminishes over time while the contribution of other QTL located at 2@91 and 9@35 increase (Fig 7a). The proportion of variance explained by 5@100 diminished late in plant development in all field grow outs (S10 Fig). However, height measurements started later and continued longer through plant development under field conditions. This revealed that the proportion of variance explained by 5@100 peaked between 30–40 days after sowing, before starting to decline. This is in contrast to its declining role throughout the growth period (8–32 days) under controlled environment conditions. Three out of four field studies also detected 9@35, and again the proportion of variance explained by this QTL diminished at late developmental stages only assessed in the field (S10 Fig). The magnitude of allelic effects contributed by individual loci generally increase through time in all growth environments (Fig 7b and S11 Fig) highlighting the cumulative nature of developmental traits. Although the identity of the major loci that influence plant height in this population is relatively unresponsive to treatment, the proportional contribution of each locus throughout time appears responsive to water treatment (Fig 7). Water-limited plants exhibit slightly decreased average and maximal rates of growth (Fig 8). In addition, the day at which maximal growth rate is achieved is significantly delayed in water-limited plants (Fig 8). QTLs associated with maximal and daily growth rate are essentially identical to those detected for height, although the time at which these QTL exhibit the greatest influence appears shifted toward earlier time points in the experiment (S12 Fig). No significant QTL were identified for either the day at which maximal growth rate occurs or the difference of this trait between treatments. The SLOD based function-valued QTL analysis of our highest resolution dataset identified a simple, additive, QTL model comprised of loci on chromosomes 2,5 and 9 (2@91,5@100 and 9@35) which provides a framework to link many of the experiments. The same causative QTL loci were identified in the wet and dry treatment blocks of the grow outs at Bellweather. Early in plant development, individuals that inherit the B100 allele at position 5@100 will be shorter than those that inherit the A10 allele (Fig 9). The effects of this locus are then counterbalanced by the influence of the QTL at position 2@91 and later by 9@35 where the B100 allele increases plant height at maturity. If the RIL genotypes are divided into four genotypic classes based on these three loci: all A10 (AAA), all B100 (BBB) and combinations where the 5@100 loci is discordant (ABA and BAB), the effect of these loci is clear. The BAB class is the tallest and the ABA class is the shortest, throughout the growth cycle. The AAA and BBB classes, which include the parent genotypes, however, have an inflection point between the early growth, where AAA is initially taller and the later stages where the BBB class becomes and remains taller for the remainder of the experiment (Fig 9). In the well-watered treatment block of the Bellwether grow out, individuals that have inherited the AAA loci exhibit the largest initial growth rate of all genotypes but are over taken by both the BBB and BAB within 13 days after planting (Fig 9c). Individuals that inherit the BBB allele combination maintain the highest rate of vertical growth at the end of this experiment in the well-watered treatment block (Fig 9c). The 5@100 QTL is found in all of the experiments with the A10 contributing the tall allele. The 2@91 and 9@35 QTLs are present in many other experiments, but there are also several additional QTLs with large effect and B100 contributing the taller allele. In most of the experiments, a fairly simple model with the 5@100 QTL driving early growth and several B100 tall QTLs taking over later explains the observed variation (Fig 10). The exceptions are three field grow outs (2013 Dry, 2014 Dry/Wet) where A10 high QTLs drive early growth, but the B100 QTLs do not appear to be the major loci influencing growth at later developmental time points. Closer inspection of the variance explained by the QTLs in these experiments reveals that at later time points, the percent variance explained by the QTLs is low, despite the trait heritability remaining high (S13 Fig). Running the QTL analysis with a lower permutation threshold cutoff (α = 0. 25) reveals the presence of several additional small effect QTLs with B100 tall alleles. Given the limitations of the population size, it is reasonable to predict that there may be more small effect QTLs that could be contributing to make B100 taller in this experiment. This chronological genetic model of plant height may in part explain the differences in plant height observed when parental lines are grown in controlled environments versus field settings. Experiments performed in controlled environmental settings potentially end before the effects of the B100 alleles at positions 2@91 and 9@35 can fully materialize or may reflect limitations in growth due to pot size constraints. The results of a glass house experiment designed to test these factors, suggests that the pot size used during in the Bellweather trial does not dramatically limit vertical growth given the duration of this experiment (experiment ends prior to 35 days after planting; Fig 11). Clearly, pot size is a significantly influential factor for larger plants, later in development (S14 Fig). The 37 unique QTL locations and their corresponding median 1. 5 LOD confidence interval were compared to the genomic location of the Setaria orthologs of plant height genes previously cloned from monocot grasses (S13 Table). Analysis using CoGe [66] identified the genomic location of syntenic orthologs in S. viridis for 16 of the 20 previously cloned height genes examined. BLAST searches were utilized to identify the location of the remaining 4 reference genes. Comparison of these two lists identified several cases of genomic colocalization (Fig 12, S14 Table). Several promising candidate genes were identified for two of the three major QTL identified during this study. The QTL of greatest influence on plant height (5@100) colocalizes with the S. viridis ortholog of rice SD-1 (semi-dwarfing), a putative gibberellin 20-oxidase enzyme [36]. This gene lies within 3. 2 kB of the peak LOD score and exhibits non-synonymous sequence polymorphisms in parental species coding sequence [16]. An alternative candidate gene associated with this QTL was located 156 kB away from the peak LOD score (outside the 1. 5 LOD confidence interval) is the S. viridis orthologue of the maize Na1, a 3-oxo-5-alpha-steroid 4-dehydrogenase enzyme required for brassinosteroid biosynthesis [19]. Two potential candidate genes appear to colocalize with QTL on 9@35. The S. viridis orthologs of maize Kn1 (15 kB from peak LOD position) and maize D8 (836 kB from peak LOD position) are located within the 1. 5 LOD confidence interval of this QTL [23,24]. These candidate genes were also identified in another study using the same interspecific RIL population [16]. None of the S. viridis orthologs in our table of previously cloned height genes overlap with the QTL confidence interval on 2@91. This colocalization approach also identified reasonable candidates for two smaller effect QTL located on 2@70 and 5@48. The S. viridis orthologue of maize Dfl1, a basic leucine zipper protein involved in modulating flowering time and plant stature [29] is located on chromosome 2 approximately 252 kB from the position of maximal LOD score for the QTL on 2@70. The QTL on 5@48 colocalizes (103 kB from position of maximal LOD score) with the orthologue of maize D1, a gibberellic acid (GA) 3-oxidase enzyme, which is the first committed step of GA biosynthesis [20]. The candidate gene for the S. italic D2 dwarf locus (15) is located just outside the confidence interval for a single environment QTL on the distal arm of Chr. 3, suggesting that this locus has only a minor, if any, role in controlling height in this population. Although several other candidate genes were found within the confidence intervals of QTL on 1@49,3@69,7@59,9@1 and 9@54, the candidates were not in particularly close proximity to the position of maximal LOD score or the QTL were only found in a single environment (Fig 12).
  Although several end point analyses of plant height have been performed in grasses [7–14], we lack an understanding of the temporal interplay of genetics and environment over development. A major impediment is the challenge of performing replicated field trials that impose treatment effects and the difficulty of repeatedly measuring any trait manually. The use of model systems, robust field design strategies and high-throughput phenotyping platforms alleviate many of these obstacles. The proportion of energy that grasses invest into vertical growth influences vegetative biomass and grain yield, which are important feed and biofuel traits [2,70,71]. As demonstrated here, using platforms with high frequency and non-destructive measurements of plant height provides temporal resolution of QTL effects. The results of 12 coordinated grow outs of a Setaria RIL population in three different experimental locations suggest that plant height is specified as a functionally continuous developmental process that is under strong and robust genetic control but can be influenced by both water availability and planting density. As expected [40–42], more aggressive water limitation treatments led to a greater reduction in plant height. However, no significant difference in the maximal plant growth rate was observed, rather it was delayed by approximately 3 days in water limited plants. Another surprising finding was that RILs planted at lower planting densities were taller than their densely planted counterpart, counter to the expectation that the shade avoidance response would result in increased plant height in densely grown stands [42,72] and may reflect resource limitations that inhibit growth of densely planted stands. Trial location was the most influential variable affecting plant height after genotype which may be the result of differences in day length and temperature influencing flowering time. In all field trials, flowering time was positively correlated with plant height at maturity. Our results indicate that the identity of the genetic components that most strongly influence plant height are not dependent upon growth environment. Rather, delayed progression through this “core” genetic program that establishes mature plant height is a key element of G x E responses. Twelve genetic loci are found in each of the four different treatment regimens and are responsible for the strongest effects in every experiment, although many smaller effect QTL were identified across individual time points in multiple independent replicated trials. Overall, the identity of genetic loci identified in this study are highly congruent with an independent study of height performed on the same population [16]. The alleles that increase plant height are found in both parental lines, with the allele of largest positive effect (5@100) inherited from the shorter A10 parental line. The effects of this locus are realized early, as growth rate approaches its maximum value. Later in development, the number of loci that contribute to plant height increases. A vast majority of the loci that increase plant height late in development are inherited from the B100 parent. The dynamics of the late-appearing genetic loci are clearly dependent upon environmental factors as they are more variable between experiments. The balance of early alleles favoring height coming from the wild parent and later alleles coming from the domesticated parent may reflect the results of domestication on growth habit, with the wild parent favoring rapid growth and flowering while the domesticated favors longer growth cycles and larger plants. Given the number of genes that have been cloned for plant height in the Poaceae, it is not surprising to find several candidate orthologs located near our identified loci. The close linkage of the 5@100 peak LOD score marker with the ortholog of SD-1 strongly suggests that this gene is responsible for the phenotype, but confirmation will require further experiments. Several other promising candidates were found for other QTLs, but several have multiple possible candidates, while others are within linkage disequilibrium, but farther away from the peak LOD SNP. Several QTLs have no known candidate orthologs within their confidence intervals, suggesting that they result from novel height altering loci. Despite stark differences in rank order among genotypes and the relative height of the parents observed within field and controlled environments, the QTL detected largely overlap. The temporal QTL model provides the explanation for these conflicting observations across experiments conducted at different time windows. Plants grown at the Bellweather Phenotyping Facility are limited to the first 30–40 days of growth due to size constraints, whereas in the field data can be collected from more mature plants but cannot be reliably collected before 20 days after planting due to the need to establish seedlings after transplanting. These results demonstrate that adding temporal analysis to phenotyping can greatly improve our ability to detect and interpret the factors underlying complex traits. The dense marker map created for this study can be combined with diverse germplasm studies to fine map causal variants and identify tightly linked markers for use in breeding programs. These tools can be used to modulate plant height and vertical growth rate throughout plant development given a predicted precipitation regimen or planting density. These approaches should be applicable across diverse species and traits for improvement of both food and biofuel crops.
  By combining field and controlled environmental analysis with quantitative genetics we were able to dissect the genetic and environmental components of plant height in Setaria. The identity of genetic components that determine plant height within this RIL population exhibits a strong degree of genetic canalization, but the contribution of these loci is influenced by temporal and environmental factors. Temporal analysis revealed that alleles from the wild parent act early in development, while alleles from the domesticated parent act later. This approach can be used to dissect all of the traits that can be measured in a high throughput manner by emerging phenotyping systems to better understand plant adaptation.
summary: Growth is a dynamic process that responds to a changing environment. Most of the methods that we have for measuring are static and collecting information throughout an organisms lifecycle is labor and cost prohibitive. Advances in imaging and robotics technology have enabled novel approaches to understanding how plants adapt to the environment. Using the model grass Setaria and new methods for measuring parameters from images, we investigate the genetic architecture of plant height in response to water availability and planting density. Height is one of the most influential components of plant architecture, determining tradeoffs between competition and resource allocation and is an important trait for boosting yields. The non-destructive nature of plant height measurements has enabled us to monitor growth throughout the plant life cycle in both field and controlled environments. We identified several loci controlling height in a population derived from a wild strain of Setaria viridis and its domesticated relative Setaria italica, as well as the developmental time in which these loci act. In this population, alleles inherited from the wild parent act to increase plant height early, whereas a larger number of small effect alleles inherited from the domesticated parent collectively act to increase plant height later in development.
keywords: quantitative trait loci
  plant science
  model organisms
  experimental organism systems
  molecular biology techniques
  plants
  flowering plants
  research and analysis methods
  grasses
  gene mapping
  maize
  plant genetics
  setaria
  molecular biology
  genetic loci
  plant and algal models
  phenotypes
  genetics
  biology and life sciences
  organisms
title: Time dependent genetic analysis links field and controlled environment phenotypes in the model C4 grass Setaria