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.

Study population

The study population of Impatiens capensis is in Glen Sutton, Quebec, Canada (45^o 02’ 37” N, 72^o 32’ 57” W). The plants occur in damp soil within an irregularly shaped area of ca. 150 m², beneath a canopy of a mixed, mature deciduous-evergreen (Acer saccharum-Tsuga canadensis) forest. I. capensis 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 (ca. 200-250 per m²) at the end of the season. Early season seedling density was at least twice as high. In I. capensis, the Pearson correlation between chasmogamous seed production and end of season biomass is r = 0.95, and between overall seed production and biomass is r = 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.

Impatiens capensis 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 Impatiens capensis (15 September), we recorded plant height and numbers of leaves from 97 plants collected at random in three, haphazardly placed 1 m² 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 (r² = 0.72, P < 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 process_radtags and clone_filter (STACKS 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 bwakit version 0.7.12 and converted to bam files using SAMtools version 1.13 (Lin et al., 2009). SNPs with quality scores > 30 were identified with gstacks (STACKS version 2.60) and further processed with populations (STACKS 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 < 0.01. The STACKS populations program was used to calculate population genetics statistics (nucleotide diversity and Wright’s Fst), 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 vcfR (version 1.12; Knaus & Grunwald, 2001) and SNPfiltR (version 1.01; DeRaad, 2022) R version 4.1.1 (R Core Team, 2021) to filter out SNPs with low read depths (< 7). To avoid spurious associations that might inflate relationship estimates of the GRM, SNPs with r² (squared correlation coefficient between the alleles at two loci) > 0.05 and within 1000 kb windows were filtered out with LDAK (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 Plink (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 LDAK. 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 KING 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 

LDAK v. 5.2 package of programs (https://dougspeed.com/ldak/).

STACKS version 2.60

bwakit v. 0.7.12 a

SAMtools v. 1.13

R version 4.1.1

vcfR v. 1.12 (an R package)

SNPfiltR v. 1.01 (an R package)

Plink v. 1.9