DOI: 10.5061/dryad.xwdbrv1jg
title: The heritability of size in a wild annual plant population with hierarchical size structure

description: The relative magnitude of additive genetic versus residual variation for
 fitness traits is important in models for predicting the rate of evolution
 and population persistence in response to changes in the environment. In
 many annual plants, lifetime reproductive fitness is correlated with
 end-of-season plant biomass, which can vary significantly from plant to
 plant in the same population. We measured end-of-season plant biomasses
 and obtained SNP genotypes of plants in a dense, natural population of the
 annual plant species Impatiens capensis with hierarchical size structure.
 These data were used to estimate the amount of heritable variation for
 position in the size hierarchy and for plant biomass. Additive genetic
 variance for position in the size hierarchy and plant biomass were both
 significantly different from zero. These results are discussed in
 relationship to theory for the heritability of fitness in natural
 populations and ecological factors that potentially influence heritable
 variation for fitness in this species.
description: Study population The study population of
 <em>Impatiens capensis</em> is in Glen Sutton, Quebec, Canada
 (45<sup>o</sup> 02’ 37” N, 72<sup>o</sup> 32’ 57”
 W). The plants occur in damp soil within an irregularly shaped area of
 <em>ca.</em> 150 m<sup>2</sup>, beneath a canopy
 of a mixed, mature deciduous-evergreen (<em>Acer saccharum-Tsuga
 canadensis</em>) forest. <em>I. capensis</em> plants in
 this population form nearly pure stands that emerge as a near continuous
 carpet of seedlings on the forest floor. The density of individuals
 remains high (<em>ca.</em> 200-250 per
 m<sup>2</sup>) at the end of the season. Early season seedling
 density was at least twice as high. In <em>I. capensis</em>,
 the Pearson correlation between chasmogamous seed production and end of
 season biomass is <em>r</em> = 0.95, and between overall seed
 production and biomass is <em>r</em> = 0.92 (Waller, 1979). 
 Small plants have been shown to produce no chasmogamous flowers or fruit
 at all (Waller, 1979), thus making position in the size hierarchy an
 interesting fitness component for study.
 <em>Impatiens capensis</em> reference
 genome From a single plant collected in the study
 population, 10 g of young leaves were harvested and frozen in liquid
 nitrogen. Genomic DNA was extracted from this tissue and sequenced by
 Dovetail Genomics/Cantata Bio LLC (Scotts Valley, California). For each
 Dovetail Omni-C library, chromatin was fixed in place with formaldehyde in
 the nucleus and then extracted. Fixed chromatin was digested with DNAse I,
 chromatin ends were repaired and ligated to a biotinylated bridge adapter
 followed by proximity ligation of adapter containing ends (Jordan
 Zhang—Dovetail Genomics, pers comm.). After proximity ligation, crosslinks
 were reversed, and the DNA purified. Purified DNA was treated to remove
 biotin that was not internal to ligated fragments, and sequencing
 libraries were generated using NEBNext Ultra enzymes and
 Illumina-compatible adapters (Jordan Zhang—Dovetail Genomics, pers comm.).
 Biotin-containing fragments were isolated using streptavidin beads before
 PCR enrichment of each library. The library was sequenced on an Illumina
 HiSeqX platform to produce approximately 30x sequence coverage.
 The input de novo assembly and Dovetail OmniC library reads were
 used as input data for HiRise, a software pipeline designed specifically
 for using proximity ligation data to scaffold genome assemblies (Putnam et
 al., 2016). Dovetail OmniC library sequences were aligned to the draft
 input assembly using bwa (https://github.com/lh3/bwa). The separations of Dovetail OmniC read pairs mapped within draft scaffolds were analyzed by HiRise to produce a likelihood model for genomic distance between read pairs, and the model was used to identify and break putative mis-joins, to score prospective joins, and make joins (Jordan Zhang—Dovetail Genomics, pers comm.). Plant biomass distribution in a random sample At the end of the 2020 growing season for <em>Impatiens capensis</em> (15 September), we recorded plant height and numbers of leaves from 97 plants collected at random in three, haphazardly placed 1 m<sup>2 </sup>qu­adrats within the study population. All 97 plants were weighed fresh to the nearest 0.01 g, and linear regression was used to establish a predictive relationship for plant biomass based on plant height and leaf number per plant (<em>r</em><sup>2</sup> = 0.72, P &lt; 0.001). Size (biomass) inequality was examined by ranking plants from lightest weight to heaviest and graphing cumulative biomass against rank order and calculating the Gini coefficient, a measure of inequality of resource distribution (Dorfman, 1979). K-means clustering (Pedregosa et al., 2011) of the biomasses of 97 randomly sampled plants was used to determine the level of support for distinct plant size clusters, as well as the biomass cutoff that separates the clusters. Genotyping-by-sequencing and population genetics analyses Leaf tissue from the small and large sampled plants was collected and preserved for DNA extraction using silica gel. DNA from samples of the preserved leaves was extracted and processed for genotyping-by-sequencing (GBS) (Elshire et al., 2011) at the University of Wisconsin Biotechnology Center. Genomic DNA was digested with the restriction enzyme ApeK1 and ligated to adapters and barcodes to create the GBS libraries. A NovaSeq6000 sequencer was used to obtain paired end (150 bp) sequence reads from the libraries. The depth of coverage was approximately 180x. Raw reads were demultiplexed and filtered for Illumina adapter sequences and PCR duplicates with <em>process_radtags</em> and <em>clone_filter</em> (<em>STACKS</em> version 2.60; Rochette, Rivera-Colon and Catchen, 2019). The demultiplexed and filtered reads were then aligned to the reference genome using the Burrows-Wheeler Alignment (Li and Durbin, 2009) tool as implemented in <em>bwakit</em> version 0.7.12 and converted to bam files using <em>SAMtools</em> version 1.13 (Lin et al., 2009). SNPs with quality scores &gt; 30 were identified with <em>gstacks </em>(STACKS version 2.60) and further processed with <em>populations</em> (<em>STACKS </em>version 2.60), which was used to filter out loci not found in at least 95% of individuals in each of the five quadrats sampled and where the minimum allele frequency for the less common SNP allele was &lt; 0.01. The <em>STACKS</em> <em>populations</em> program was used to calculate population genetics statistics (nucleotide diversity and Wright’s <em>Fst</em>), for reporting the results of SNP filtering, and to produce a vcf file of SNP genotypes for each sampled individual. Prior to the estimation of the genomic relationship matrix (GRM; see below) we applied a <em>vcfR </em>(version 1.12; Knaus & Grunwald, 2001) and <em>SNPfiltR </em>(version 1.01; DeRaad, 2022) R version 4.1.1 (R Core Team, 2021) to filter out SNPs with low read depths (&lt; 7). To avoid spurious associations that might inflate relationship estimates of the GRM, SNPs with <em>r</em><sup>2</sup> (squared correlation coefficient between the alleles at two loci) &gt; 0.05 and within 1000 kb windows were filtered out with <em>LDAK </em>(version 5.2; https://dougspeed.com/ldak; Speed et al., 2020), which used bed genotype files as input, created from the original vcf file using <em>Plink </em>(version 1.9, Chang et al., 2015). Analyses of genetic variance and heritable variation for position in the size hierarchy and plant biomass We treated size either as a threshold trait and analyzed it on a liability scale (small plants versus large plants) as is done in some agricultural genetic studies for traits such as secondary compound content, and disease resistance (Merrick et al., 2023), or we directly assayed biomass. We used four different methods to estimate heritability of plant size. These are implemented in the <em>LDAK. </em>The first two methods, Phenotype Correlation–Genotype Correlation (PCGC) and TetraHer, use a liability threshold model, such that the binary outcome (small versus large)  indicate whether the unobserved liability is above or below a threshold. Both PCGC and TetraHer estimate heritability on the liability scale by measuring the extent that the estimated relatedness between pairs of individuals correlates with their estimated liabilities. PCGC considers all pairs of individuals, and measures pairwise relatedness based on allelic correlations (inferred from the GRM). By contrast, TetraHer considered only the 3,606 pairs of individuals identified as having at least 17.5% (IBD) by descent using the kinship analysis <em>KING</em> software (Manichaikul et al., 2010). Both methods adjust for ascertainment (i.e., the fact that our sample was enriched for large plants, relative to the natural population). Prevalence of large plants was determined as the proportion of large plants observed in the random sample of 97 plants, as classified by K-means clustering—see above). We additionally obtained estimates of heritability on the observed scale (a continuous variable). For this we used restricted maximum likelihood (REML) and QuantHer, which are similar to PCGC and TetraHer. We note that sampling weights from the smallest and largest plants for REML and PCGC constitutes “extreme sampling”, but this has been shown previously to have a minimal biasing effect on heritability estimation (Golan et al., 2014). Nevertheless, we also conducted our own simulation analyses to gauge the possible effects of sampling protocol on our own results (Supplementary Methods for details). References  Anderson, J. T. (2016). Plant fitness in a rapidly changing world. <em>New Phytologist</em> 210:81-87. Bérénos, C., Ellis, P. A., Pilkington, J. G., & Pemberton, J. M. (2014). Estimating quantitative genetic parameters in wild populations: A comparison of pedigree and genomic approaches. <em>Molecular Ecolog</em>y, 23(14), 3434-3451. Bontemps, A., Lefèvre, F., Davi, H., & Oddou‐Muratorio, S. (2016). In situ marker‐based assessment of leaf trait evolutionary potential in a marginal European beech population. <em>Journal of Evolutionary Biology</em>, 29(3), 514-527. Burt, A. (1995). The evolution of fitness. <em>Evolution</em>, 49(1), 1-8. Castellanos, M. C., González‐Martínez, S. C., & Pausas, J. G. (2015). Field heritability of a plant adaptation to fire in heterogeneous landscapes. <em>Molecular Ecology</em>, 24(22), 5633-5642. Chang, C. C., Chow, C. C., Tellier, L. C., Vattikuti, S., Purcell, S. M., & Lee, J. J. (2015). Second-generation PLINK: rising to the challenge of larger and richer datasets. <em>Gigascience</em>, <em>4</em>(1), s13742-015. Chen, J., Glémin, S., & Lascoux, M. (2017). Genetic diversity and the efficacy of purifying selection across plant and animal species. <em>Molecular biology and evolution</em>, 34(6), 1417-1428. Day, P. D., Pellicer, J., & Kynast, R. G. (2012). Orange balsam (Impatiens capensis Meerb., Balsaminaceae): A re-evaluation by chromosome number and genome size. <em>The Journal of the Torrey Botanical Society</em>, 139(1), 26-33. DeRaad, D. A. (2022). SNPfiltR: An R package for interactive and reproducible SNP filtering. <em>Molecular Ecology Resources</em>, 22(6), 2443-2453. Donohue, K. (2002). Germination timing influences natural selection on life‐history characters in <em>Arabidopsis thaliana</em>. <em>Ecology</em>, 83(4), 1006-1016. Dorfman, R. (1979). A formula for the Gini coefficient. <em>The Review of Economics and Statistics</em>, 61, 146-149. Elshire, R. J., Glaubitz, J. C., Sun, Q., Poland, J. A., Kawamoto, K., Buckler, E. S., & Mitchell, S. E. (2011). A robust, simple genotyping-by-sequencing (GBS) approach for high diversity species. <em>PloS One</em>, 6(5), e19379. Fisher R.A. (1930). The genetical theory of natural selection. Clarendon Press, UK. Gienapp, P., Fior, S., Guillaume, F., Lasky, J. R., Sork, V. L., & Csilléry, K. (2017). Genomic quantitative genetics to study evolution in the wild. <em>Trends in Ecology & Evolution</em>, 32(12), 897-908. Golan, D., Lander, E. S., & Rosset, S. (2014). Measuring missing heritability: Inferring the contribution of common variants. <em>Proceedings of the National Academy of Sciences</em>, 111(49), E5272-E5281. Gomulkiewicz, R., & Holt, R. D. (1995). When does evolution by natural selection prevent extinction? <em>Evolution</em>, 49, 201-207. Harper, J. L. (1977). <em>Population biology of plants</em>. New York, NY: Academic Press. Hill, W. G., & Robertson, A. (1968). Linkage disequilibrium in finite populations. <em>Theoretical and Applied Genetics</em>, 38, 226-231. Hill, W. G., & Zhang, X. S. (2009). Maintaining genetic variation in fitness. In J. van der Werf, H. U. Graser, R. Frankham, & C. Gondro (Eds.), <em>Adaptation and Fitness in Animal Populations</em> (pp. 67-85). Dordrecht: Springer. https://doi.org/10.1007/978-1-4020-9005-9_5 Hohenlohe, P. A., Bassham, S., Etter, P. D., Stiffler, N., Johnson, E. A., & Cresko, W. A. (2010). Population genomics of parallel adaptation in threespine stickleback using sequenced RAD tags. <em>PLoS Genetics</em>, 6(2), e1000862. Johnston, S. E., Chen, N., & Josephs, E. B. (2022). Taking quantitative genomics into the wild. <em>Proceedings of the Royal Society B</em>, 289, 20221930. Knaus, B. J., & Grünwald, N. J. (2017). vcfr: A package to manipulate and visualize variant call format data in R. <em>Molecular Ecology Resources</em>, 17(1), 44-53. Kruuk, L. E., Clutton-Brock, T. H., Slate, J., Pemberton, J. M., Brotherstone, S., & Guinness, F. E. (2000). Heritability of fitness in a wild mammal population. <em>Proceedings of the National Academy of Sciences</em>, 97(2), 698-703. Kulbaba, M. W., Sheth, S. N., Pain, R. E., Eckhart, V. M., & Shaw, R. G. (2019). Additive genetic variance for lifetime fitness and the capacity for adaptation in an annual plant. <em>Evolution</em>, 73(9), 1746-1758. Li, H., & Durbin, R. (2009). Fast and accurate short read alignment with Burrows–Wheeler transform. <em>Bioinformatics</em>, 25(14), 1754-1760. Mackay, T. F., & Lyman, R. F. (2005). Drosophila bristles and the nature of quantitative genetic variation. <em>Philosophical Transactions of the Royal Society B: Biological Sciences</em>, <em>360</em>(1459), 1513-1527. Merilä, J., & Sheldon, B. C. (1999). Genetic architecture of fitness and nonfitness traits: Empirical patterns and development of ideas. <em>Heredity</em>, 83, 103-109. Mitchell-Olds, T. (1986). Quantitative genetics of survival and growth in <em>Impatiens capensis</em>. <em>Evolution</em>, 40(1), 107-116. Mitchell-Olds, T., & Bergelson, J. (1990). Statistical genetics of an annual plant, Impatiens capensis. II. Genetic basis of quantitative variation. <em>Genetics</em>, 124(2), 407-415. Mousseau, T. A., & Roff, D. A. (1987). Natural selection and the heritability of fitness components. <em>Heredity</em>, <em>59</em>(2), 181-197. Peschel, A. R., Boehm, E. L., & Shaw, R. G. (2021). Estimating the capacity of <em>Chamaecrista fasciculata</em> for adaptation to change in precipitation. <em>Evolution</em>, <em>75</em>(1), 73-85. Peschel, A. R., & Shaw, R. G. (2024). Comparing the predicted versus realized rate of adaptation of <em>Chamaecrista fasciculata</em> to climate change. <em>The American Naturalist</em>, <em>203</em>(1), 14-27. Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, et al. (2011). Scikit-learn: Machine learning in Python. <em>the Journal of machine Learning research</em>, <em>12</em>, 2825-2830. Putnam, N. H., O'Connell, B. L., Stites, J. C., Rice, B. J., Blanchette, M., Calef, R., & Green, R. E. (2016). Chromosome-scale shotgun assembly using an in vitro method for long-range linkage. <em>Genome Research</em>, 26(3), 342-350. Price, T., & Schluter, D. (1991). On the low heritability of life‐history traits. <em>Evolution</em>, 45(4), 853-861. R Core Team (2021). RL A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/. . Rochette, N. C., Rivera‐Colón, A. G., & Catchen, J. M. (2019). Stacks 2: Analytical methods for paired‐end sequencing improve RADseq‐based population genomics. <em>Molecular ecology</em>, 28(21), 4737-4754. Schmitt, J., Eccleston, J., & Ehrhardt, D. W. (1987). Dominance and suppression, size-dependent growth and self-thinning in a natural <em>Impatiens capensis</em> population. <em>The Journal of Ecology</em>, 75, 651-665. Schwaegerle, K. E., & Levin, D. A. (1991). Quantitative genetics of fitness traits in a wild population of Phlox. <em>Evolution</em>, 45(1), 169-177. Schwinning, S., & Weiner, J. (1998). Mechanisms determining the degree of size asymmetry in competition among plants. <em>Oecologia</em>, <em>113</em>, 447-455. Shaw, R. G., Platenkamp, G. A., Shaw, F. H., & Podolsky, R. H. (1995). Quantitative genetics of response to competitors in <em>Nemophila menziesii</em>: a field experiment. <em>Genetics</em>, <em>139</em>(1), 397-406. Shaw, R. G., & Etterson, J. R. (2012). Rapid climate change and the rate of adaptation: insight from experimental quantitative genetics. New Phytologist, 195(4), 752-765. Shaw, R. G., & Shaw, F. H. (2014). Quantitative genetic study of the adaptive process. <em>Heredity</em>, 112(1), 13-20. Sheth, S. N., Kulbaba, M. W., Pain, R. E., & Shaw, R. G. (2018). Expression of additive genetic variance for fitness in a population of partridge pea in two field sites. <em>Evolution</em>, <em>72</em>(11), 2537-2545. Seppey, M., Manni, M., & Zdobnov, E. M. (2019). BUSCO: assessing genome assembly and annotation completeness. <em>Gene prediction: methods and protocols</em>, 227-245. Silvertown, J., & Charlesworth, D. (2009). <em>Introduction to plant population biology</em>. John Wiley & Sons. Speed, D., Holmes, J., & Balding, D. J. (2020). Evaluating and improving heritability models using summary statistics. <em>Nature Genetics</em>, 52(4), 458-462. Stanton‐Geddes, J., Yoder, J. B., Briskine, R., Young, N. D., & Tiffin, P. (2013). Estimating heritability using genomic data. <em>Methods in Ecology and Evolution</em>, 4(12), 1151-1158. Stevens, L., Goodnight, C. J., & Kalisz, S. (1995). Multilevel selection in natural populations of <em>Impatiens capensis</em>. <em>The American Naturalist</em>, 1<strong>45</strong>(4), 513-526. Thomas, S. C., & Bazzaz, F. A. (1993). The genetic component in plant size hier-archies: norms of reaction to density in a Polygonum species. <em>Ecological Monographs</em>, 63(3), 231-249. Turner, M. D., & Rabinowitz, D. (1983). Factors affecting frequency distributions of plant mass: the absence of dominance and suppression in competing monocultures <em>of Festuca paradoxa. Ecology</em>, 64(3), 469-475. Waller, D. M. (1985). The genesis of size hierarchies in seedling populations of <em>Impatiens capensis</em> Meerb. <em>New Phytologist</em>, 100(2), 243-260. Weiner, J. (1985). Size hierarchies in experimental populations of annual plants. <em>Ecology</em>, 66(3), 743-752. Weiner, J. A. (1988) The influence of competition on plant reproduction” In J. L. Doust & J. D. Doust (Eds.), Plant reproductive ecology: patterns and strategies (pp. 228-245). Oxford Univ. Press, NY. Weiner, J. (1990). Asymmetric competition in plant populations. Trends in Ecology & Evolution, 5(11), 360-364. Yang, J., Lee, S. H., Goddard, M. E., & Visscher, P. M. (2011). GCTA: a tool for genome-wide complex trait analysis. <em>The American Journal of Human Genetics</em>, 88(1), 76-82.
description: <em>LDAK</em>
 v. 5.2 package of programs
 (https://dougspeed.com/ldak/). STACKS version 2.60 bwakit v. 0.7.12 a <em>SAMtools</em> v. 1.13 R version 4.1.1 <em>vcfR</em> v. 1.12 (an R package)<br> <em>SNPfiltR</em> v. 1.01 (an R package)<br> Plink v. 1.9
description: \**­­­­Title of Data Set: ** Heritability of fitness in a wild annual
 plant population with hierarchical size structure (*Note: This data set
 was updated on 19 June 2024, after discovering several quality control
 issues with the genomic data. We apologize for any inconvenience*). This
 data set includes measures of plant weight (fresh biomass in g) from a
 sample of 97 plants collected in three quadrats in the study population at
 random (random sample), and from 367 plants collected in five additional
 quadrats where large and small plants were sampled (extreme sample). The
 367 plants in the extreme sample were genotyped using
 genotyping-by-sequencing, and the vcf files for these plants and their SNP
 genotypes are also included in this DRYAD repository. The BAM files used
 for constructing this vcf file along with the *Impatiens capensis*
 reference genome are deposited in GenBank under SRA accession number
 PRJNA945897. **Description of the data and file structure** File:
 *Impatiens_capensis_geneome_assembly_MCG3086_report.pdf   *Hi-rise
 scaffolding report for reference genome. File:
 *Phenotypes_Random_sample_N=97.xlsx*    Plant height (cm), Leaf number per
 plant, and plant weight (fresh biomass in g) from a sample of 97 plants of
 *Impatiens capensis* (random sample) in an Excel file. File:
 *Plant_weights_extreme_sample.txt * Plant weight (fresh biomass in g) from
 a sample of 367 plants (186 small and 181 large plants: see text of
 article) of *Impatiens capensis* (extreme sample) in a text file with 3
 columns. The first two columns combined are plant identifiers and
 correspond to plant identifiers in the file
 “Impatiens_367_genotyped_individuals.vcf” File:
 *Impatiens_367_genotyped_individuals.vcf* Genotypic data used along with
 data in the file “Plant_weights_extreme_sample.txt” to estimate the
 heritability of plant weight and position in the size hierarchy using the
 scripts in the file ‘LDAK.code File: *Position_map.txt* Positions
 (quadrats) of each genotyped plant. **Code** File:
 *Simulation_Extreme_Sampling.R*  R program to simulate extreme sampling to
 examine the potential bias in estimating heritability and its significance
 level.  Uses data from File: *LDAK_code  *Code used with the LDAK software
 version 5.2 ([https://dougspeed.com/ldak/](https://dougspeed.com/ldak/))
 to estimate the heritability of plant weight and position in the size
 hierarchy. File: *K_means.py *Code for K means analysis of plant weight
 data in random sample. File: *VCFR.py *Code for quality control. See
 article. File: *Snakefile* Code for processing raw sequence reads and
 producing bam files for input into STACKS. Used to automate the processing
 of the sequence data, including trimming adapters, creating sorted .bam
 files, and removing duplicates. File: *Linearmodel_wt_f(ht+lf_no).py.*
 Used to analyze the relationship between individual plant biomass versus
 leaf number per plant via linear model analysis. File:
 *Impatiens_capensis_hirise.fasta*: This file contains all assembled
 scaffolds for the *Impatiens capensis* genome assembly.  Scaffolds 1
 through 9, and scaffold 11 account for all but 27 of them 36,044 genotyped
 regions that map to the hirise assembly (with scaffolds 10 and 12
 accounting for only 1 and 26, respectively).