Assessing the opportunity for selection to impact morphological traits in crosses between two Solanum species

Solanum cohort generation

Seeds for the original cohorts (S. lycopersicum (strain VF36), S. pennellii (strain 716), and F1 hybrids (strain 4135)) were obtained from the C.M. Rick Tomato Genetics Resource Center at the University of California, Davis. In September 2021, using seeds from the three strains from the resource center, we generated the crosses shown in Table 1. F1 seeds were used to generate the crosses but were unavailable to grow for data collection.

Seeds were treated with a 20% commercial bleach solution for 20 min and then placed in soil for germination. Seedlings were grown in a growth chamber (23 °C, 12-h day/night lighting, and 50% humidity) in four-inch pots. One month after germination, seedlings were transplanted into eight-inch pots and relocated to the experimental greenhouse at Texas A&M University (30.615 N and 96.339 W). The greenhouse temperature was maintained at 23 °C, while lighting and humidity were allowed to fluctuate. After seeds were generated for all cohorts, they were washed, allowed to dry completely, and stored until planted for data collection in August 2022.

Trait data collection

Seeds were planted in August of 2022. F1 individuals were used for generating the new cohorts, but F1 individuals were not available for data collection after the cohort set had been generated. Three tables in the greenhouse were used and cohorts of each type were distributed equally on all tables, ensuring that no cohorts experienced biased microclimates within the growing area. For the five cohorts being measured (P1, P2, F2, BC1, rBC2), mature branches and leaves were collected 77 days after planting. Specimens were pressed using a plant press for 27 days in a well-ventilated area and then fixed on blotting paper with pH-neutral PVA glue. The specimens were laid on white paper containing a black reference square (1 cm × 1 cm). The specimens were scanned as JPEG images using a Brother DCP L2540DW printer, producing an image like Fig. 1A (left). Scanned images were used for leaf area, leaf perimeter, leaf length, and leaf width measurements. The images were edited using ImageJ (Schneider, Rasband & Eliceiri, 2012). The images were cropped, and the leaves were isolated from their stems and other leaves using the “Paintbrush Tool”. Despeckling was applied to fill in the leaves to allow for accurate leaf measurements. The threshold was manually adjusted to convert the image to a silhouette, resulting in an image seen in Fig. 1A (right). The leaf midvein was not easily identified in the scanned images, so the samples were photographed using a Canon EOS 2000D camera to obtain images suitable for areal ratio measures (see “Areal Ratio” section for details). These methods were done for 85 leaves for which measurements were taken. With leaves collected from two plants per cohort, there were 23 leaves for P1, 14 leaves for P2, 13 leaves for F2, 14 leaves for BC1, and 21 leaves for rBC2. Given that multiple leaves were taken from each plant, in our statistical analyses the leaves from each plant were treated as non-independent individuals.

Leaf area: Leaf area was determined through ImageJ analysis and unit conversions. ImageJ measurements were set to calculate the area and display the results. The leaf was selected from the image with the “Wand (tracing) tool”, and the leaf area was found in pixels using the “Measure” command. The area of the reference square was calculated using the same method. The areas were converted to centimeters squared by dividing the pixel count for the leaf area by the pixel count for the 1 cm² reference square.

Leaf perimeter: For leaf perimeter, ImageJ measurements were set to calculate the perimeter and display the results. The leaf was selected from the image with the “Wand (tracing) tool”, and the leaf perimeter was found in pixels using the “Measure” command. The perimeter of the reference square was calculated in the same manner. The perimeters were converted to centimeters using the known perimeter of the reference square.

Leaf perimeter-area ratio: The ratio of perimeter to area was calculated using the previously determined calculations of leaf perimeter and leaf area as proposed by Yu et al. (2019). The leaf perimeter was divided by leaf area, resulting in units of cm⁻¹. We used the leaf perimeter-area ratio as a proxy for quantifying leaf shape. Larger perimeter-area ratios indicate a leaf with more or larger lobes, while smaller perimeter-area ratios indicate a leaf with less or smaller lobes.

Leaf length: In ImageJ, leaf length was measured as the longest perpendicular extension from a straight line at the protruding base of the leaf to a parallel line at the apex of the leaf. Leaf length was defined with parameters similar to Schrader et al. (2021) and shown in Fig. 1B. ImageJ measurements were set to calculate length and display the results. A line was drawn on the leaf using the “Straight Line tool” and “Angle tool” according to the specifications outlined above, and leaf length was calculated in pixels with the “Measure” command. The length of the reference square (in pixels) was also calculated using the “Straight Line tool” and “Measure” commands. The leaf length was converted to centimeters using the known length of the reference square.

Leaf width: In a similar manner to leaf length, leaf width was measured in ImageJ as the longest extension from any two points on the edges of the leaf perpendicular to the length axis (described in the previous section). The leaf width axis was measured according to parameters like Schrader et al. (2021) and shown in Fig. 1B. ImageJ measurements were set to calculate width and display the results. A line was drawn on the leaf using the “Straight Line tool” and “Angle tool” according to the specifications outlined above, and leaf width was calculated in pixels using the “Measure” command. The width of the reference square is the same as the length of the reference square found above. The leaf width was converted to centimeters using the known width of the reference square.

Leaf width-length ratio: Once the length and width of the leaf were determined, leaf width was divided by leaf length to calculate the ratio of width to length to indicate leaf roundness. The leaf width-length ratio was measured as a proxy for leaf roundness. Ratios close to zero indicate leaves elongated along the length axis, while ratios larger than one indicate leaves elongated along the width axis, and intermediate ratios close to one indicate more circular leaves.

Areal ratio: As a measure of leaf symmetry, the areal ratio of the leaf was determined by finding the ratio of the left half area to the right half area of the leaf, as done by Yu et al. (2019). Each leaf was laid next to a reference square (1 cm × 1 cm) on a flat white surface and photographed using a Canon EOS 2000D camera. The camera was mounted directly above the leaves at a consistent height (approximately 22.5 cm above the laid-out leaves). Direct lighting was applied to the leaves for photographing using an overhead light positioned 35.5 cm above the stage. Auto aperture settings were used for all photos. Canon EOS Utility software was used for remote shooting to minimize vibration. The highest contrast, sharpness, saturation, and color tone were used, resulting in images like Fig. 1C (left).

Using ImageJ, the images were edited with the “Paintbrush Tool” to remove the midvein and split the leaf into the left and right halves. The left and right halves were determined by orienting the image with the apex at the top of the screen and the base at the bottom. The leaves and square were converted to black and white by adjusting the threshold, and they were filled in, where necessary, to be solid black, as seen in Fig. 1C (right). Measurements were set to calculate the area. Each leaf half was selected from the image using the “Wand (tracing) tool”, and the area was measured in units of pixels with the “Measure” command. The area of the reference square was calculated using the same method. Leaf areas were converted to centimeters using the known area of the reference square. The area of the left half was divided by the area of the right half to determine the areal ratio of the left half to the right half as an indicator of symmetry.

Seed mass: Various seeds from several dates (09/13/2021–11/17/2021) during one seasonal harvest were collected, washed, and allowed to dry completely in a ventilated area on blotting paper. Seeds were individually weighed in grams. Fifteen seeds were weighed for each cohort, with the exception of rBC2, which had 13 seeds.

Line cross analysis

The cohorts included in this study define which genetic effects can be inferred. We include both parental cohorts (P1 and P2), an F2 cohort, and two backcrosses (BC1 and rBC2), allowing us to make inferences for ten different genetic effects. Effects include three types of additive genetic effects (autosomal additive, cytotype additive, and maternal effect additive), two types of dominance genetic effects (autosomal dominance and maternal effect dominance), and five types of epistatic interactions (autosomal additive-by-additive, autosomal dominance-by-dominance, autosomal additive-by-dominance, autosomal additive-by-cytotype additive, and autosomal dominance-by-cytotype additive) (Table 2). We recognize this is a relatively small cohort set, however previous studies have shown the ability of LCA to infer epistatic genetic effects with small cohort sets that align with inferences of epistasis within larger cohort sets. Burch et al. (2024) analyzed 40 datasets with a cohort set of 16, then reanalyzed the same datasets after they were reduced to a smaller cohort set of five. They found a slight, non-significant bias toward detecting less epistasis in datasets with fewer cohorts.

Line cross analysis was used to estimate the contributions of the composite genetic effects (additive, dominance, epistatic, etc.) by analyzing the mean and standard error of a trait for each cohort (P1, P2, F2, BC1, rBC2). Since each leaf was treated as an individual, and multiple leaves were collected and measured from each plant, we used an alternative approach to collecting means and standard errors for each dataset. This ensures that leaves measured from the same plant were treated as non-independent, allowing for more accurate quantification of the genetic architecture. Means were collected by taking the mean of all individual measurements from a single plant, repeating that for each plant, then a cohort mean was calculated as the mean of the plant means. For the standard error, we used the following Eq. (1):

(1) $$SE=\sqrt{{\sum }_{i=1}^n{\displaystyle \frac{{\sigma }_i^2}{s_i}}}$$where $\sigma$ is the standard deviation of the trait from samples of a plant, $n$ is the number of plants, and $s$ is the number of samples from that plant.

We used an information-theoretic LCA approach implemented in the R package SAGA 2.0 (Armstrong, Anderson & Blackmon, 2019; Blackmon & Demuth, 2016). This software automatically generates an appropriate matrix of coefficients (c-matrix), which describes the opportunity for each of the possible composite genetic effects to impact the trait in each cohort (Lynch & Walsh, 1998). After the c-matrix is generated, SAGA uses a weighted least squares regression to fit the measured phenotypes with models of the genetic architecture. From these models, we can evaluate the probability for each genetic effect to impact the phenotype of a cohort, thereby inferring the genetic architecture of the trait. The software uses a small sample size corrected version of the standard AIC (AICc) to account for small sample sizes that are typical in LCA studies. AICc was used to calculate model weights and construct a 95% confidence set of models that are used to generate model-averaged results, accounting for model selection uncertainty (Burnham & Anderson, 2003). Phenotypic data was not included for the F1 cohort, so the c-matrix was customized to exclude the F1. Under the analyzed cohort set, the contribution of ten composite genetic effects were evaluated: additive, dominance, cytotype, maternal effect additive and dominance, and epistatic interactions between additive, dominance, and cytotype additive effects (Table 2).

Results were pooled in various ways to illustrate trends for genetic effects in each trait (e.g., additive, dominance, epistatic, maternal, etc.). Using criteria proposed by Blackmon & Demuth (2016), genetic effects were excluded if they did not have a variable importance greater than 0.5 and a standard error that excluded zero. Variable importance is the sum of the model weights in the confidence set of models that include that particular composite genetic effect. To allow for comparison among traits, we also scaled the contributions of genetic effects to sum to one. Briefly, we took the absolute values of the contributions, summed all significant genetic effects, and divided the estimated effects by this sum.

Epistatic effects

Estimating the role of epistatic interactions within the genetic architectures underlying trait divergence can help to understand how accessible a trait is to selection (Roff & Emerson, 2006). For instance, traits more closely aligned with fitness typically show more epistasis than traits less aligned with fitness (Burch et al., 2024; Crnokrak & Roff, 1995; Roff & Emerson, 2006).

Traits closely associated with fitness (i.e., reproductive traits) are usually the target of selective breeding. We often see more epistatic genetic effects and fewer additive genetic effects within the genetic architectures of these traits. Through selection on these traits, additive genetic effects are exhausted as beneficial alleles are fixed, leaving behind a genetic architecture composed primarily of epistatic genetic effects (Burch et al., 2024; Walsh & Lynch, 2018). By determining the contribution of epistatic effects to the trait divergence of these phenotypes, we were able to estimate how accessible these traits could be to selection given the absence of additive genetic effects. In Table 3, we show the observed traits in the order of their epistatic interactions, with higher numbers indicating more epistatic contribution to trait divergence. Four traits (leaf perimeter-area ratio, leaf perimeter, leaf width, and leaf area) had an epistatic contribution greater than 0.5, meaning that their genetic architectures were dominated by epistatic effects. One trait (seed mass) had 7% epistatic contribution, while the remaining traits had zero.

Maternal effects

In our inference of genetic architectures, we show significant maternal effects for two of eight traits: leaf width and seed mass (Table 3). Leaf width shows additive maternal genetic effects dominating 37% of the trait’s genetic architecture. Dominance maternal genetic effects account for 15% of the seed mass genetic architecture. Maternal effects are often inferred in analyses of traits in animals (McAdam et al., 2002; Willham, 1972). While Cockerham (1963) suggested that maternal effects in plants were minimal and did not require consideration, we have shown that two of the eight observed traits display maternal genetic effects as a significant proportion of their genetic architecture and that quantifying maternal effects is key to understanding the dynamics of trait evolution.

Ours is not the first study to identify significant maternal genetic effects in plants. In Lolium multiflorum, maternal effects were found to underlie the variation in third leaf size, as leaf length was due to additive maternal effects, and leaf width was due to interactions between maternal effects (Edwards & Emara, 1970). Lolium species show maternal effects in the size of leaves early in development, but maternal effects diminish during later stages of development (Roach & Wulff, 1987). A study in Cannabis sativa showed that maternal additive effects influence cannabichromene expression (Campbell, Dufresne & Sabatinos, 2020). Taken together with previous work, our study demonstrates that maternal effects are a key component of the genetic architecture of traits in plants. By quantifying the genetic architecture of traits in model plant species, we are able to improve our understanding of how accessible a trait is to selection, thereby enhancing breeding efforts across the field of agriculture.

Compound traits

Another challenge to our understanding of genetic architecture involves our definition of traits and how it impacts inference. For instance, in this study, we analyzed eight datasets, six of which (area, perimeter, length, width, areal ratio, seed mass) can be reasonably thought of as elemental traits. We also include two compound traits (perimeter-area ratio, width-length ratio) that are measured as the quotient of two elemental traits. However, the expectation for the genetic architectures of compound traits (i.e., a trait that is a function of more than one elemental trait) has not, to our knowledge, been explored. For instance, what genetic architectures will be inferred if two traits that contain only additive variation are measured and converted to a ratio and then analyzed in an LCA framework? Our study provides a single data point towards understanding what we might expect. Specifically, we found that divergence in the elemental trait of leaf width was dominated by epistatic effects (63% of total composite genetic effects) and to a lesser extent exhibited maternal effects (37% of total composite genetic effects). In contrast, the elemental trait of leaf length was composed of only autosomal additive effects. Finally, when we compare this to the compound trait of leaf width-length ratio we see that it, like leaf length, is dominated by additive effects, but both autosomal additive and cytotype additive effects are observed. However, whether this concordance between compound traits and underlying elemental traits is typical is, to our knowledge, unknown. A systematic study should evaluate the expectations for compound traits across a range of underlying elemental trait genetic architectures.