The impact of urbanization on outcrossing rate and population genetic variation in the native wildflower, Impatiens capensis

Abstract

Introduction

Urbanization is one of the most prevalent forms of landscape transformations and poses a serious threat to biodiversity. Urbanization increases habitat fragmentation by replacing natural habitats with artificial surfaces, such as roads and buildings (Cadenasso et al. 2013). Urban populations of native species tend to be smaller and less connected than nonurban populations (Johnson, Thompson, and Saini 2015), which is predicted to reduce genetic diversity within populations and increase genetic differentiation among populations (Johnson and Munshi-South 2017; Miles et al. 2019; Schmidt et al. 2020). Because adaptive capacity relies on standing genetic variation within a population (Barrett and Schluter 2008), a loss of genetic diversity through increased genetic drift may hinder the ability of populations to adapt to ongoing urbanization, possibly leading to the extirpation of small urban populations. This effect may be exacerbated if urban development creates barriers to dispersal, leading to reduced gene flow between urban and rural populations (Schmidt et al. 2020). If urbanization reduces gene flow, then urban populations may become genetically differentiated from rural populations and/or contain a subset of the genetic variation found in rural populations. Establishing the effects of urbanization on population genetic diversity and differentiation, as well as individual mating strategies (e.g. outcrossing rates), is an important step to developing strategies for cities to conserve native species that provide important ecosystem services, such as resources for pollinators (Dearborn and Kark 2010).

The effects of urbanization on genetic variation may be exacerbated in plants when population sizes are reduced and/or interactions with pollinators are disrupted (Eckert et al. 2010). Small plant populations often exhibit increased biparental inbreeding due to fewer potential mates because individuals are more likely to be related to one another (Uyenoyama 1986; Barrett and Harder 2017). Urbanization also alters interactions with pollinators (Rivkin et al. 2020), where pollinators spend more time visiting multiple flowers on a single plant or within a single plant patch in urban habitats (Wenzel et al. 2020). This behaviour can alter the outcrossing rate of a population by increasing the incidence of geitonogamous self-fertilization (i.e. fertilization via the transfer of pollen from one flower to another on the same plant) of self-compatible species. As a result, urban plants may exhibit higher self-fertilization rates (Cheptou and Avendaño 2006), which is predicted to reduce genetic diversity within populations, lead to increased expression of deleterious mutations and may ultimately decrease the long-term survival of a population (Aguilar et al. 2019). Consequently, identifying the effects of urbanization on both outcrossing rate and genetic variation in self-compatible plant species is essential to understanding the ability of native plants to persist in cities.

The goal of our study was to investigate the effects of urbanization on outcrossing rate, population genetic diversity and population genetic structure of Impatiens capensis (Meerb). Impatiens capensis is ideally suited to test questions regarding mating system variation because it exhibits a mixed mating system that has the potential to be strongly influenced by the environmental changes associated with urbanization (Barker and Sargent 2020). Impatiens capensis is also a native plant in our study area that grows in habitats that are heavily impacted by urbanization and may experience altered patterns of genetic diversity and structure as urban development progresses. Lastly, I. capensis is an important nectar resource for native bumblebees and hummingbirds (Schemske 1978), thus documenting the effects of urbanization on outcrossing and population genetic variation is the first step to mitigating the impact of urban growth on this species.

Methods

Study system and sampling

Impatiens capensis is a shade-tolerant, annual herb native to eastern North America. Populations thrive in wet forest depressions, ditches and along the margins of wetlands, such as marshes, ponds, lakes and rivers (Gleason and Cronquist 1963). Populations grow in dense stands (Schmitt and Ehrhardt 1990), and flowering begins in mid-summer and lasts until the first frost (Schemske 1978). Impatiens capensis exhibits a mixed mating system, whereby each plant produces both open, outcrossing (chasmogamous) and closed, obligately self-fertilizing (cleistogamous) flowers. Chasmogamous flowers are pollinated by Bombus spp., Apis mellifera and Archilochus colubris (Ruby-Throated Hummingbird; Schemske 1978). Chasmogamous flowers are strongly protandrous and rarely autogamous; however, the outcrossing rate of the chasmogamous flowers is typically less than one due to geitonogamous self-fertilization (Waller and Knight 1989). Relative to chasmogamous flowers, cleistogamous flowers are drastically reduced in size and are only capable of autogamous self-fertilization. Fruits produced by cleistogamous flowers are typically smaller and contain fewer seeds than fruits produced by chasmogamous flowers (Waller 1984).

On 17 September 2019, we sampled parent–offspring pairs (hereafter: family) from three urban and three rural populations in Toronto, Ontario, Canada (Fig. 1; Table 1). Impatiens capensis is native to the region where we were sampling, and the historical farming land-use of the region has likely led to I. capensis populations persisting in small, remnant populations located along the margins of watersheds and in urban parks. Urban populations were located within the City of Toronto borders and were found in parks that were surrounded by dwellings. In contrast, rural populations were located outside the city limits and were surrounded by farmland and few, if any, buildings (Fig. 1). We quantified the percentage of impervious surfaces within 1 km² of each site using the Global Man-made Impervious Surface Dataset (Brown de Colstoun et al. 2017). Urban populations were surrounded by approximately 30% impervious surfaces, whereas rural populations were surrounded by less than 4% impervious surfaces (Table 1). Because the mating system of I. capensis responds plastically to the environment (Waller and Knight 1989), we ensured that each population contained similar levels of light availability and were all located in similar environments on the border of streams or small ponds. Lastly, seed dispersal in I.capensis occurs primarily along waterways, so we ensured that the sites we selected were not part of the same watershed system to avoid sampling populations that were likely to be genetically related.

Figure 1:

Map of the study area in Toronto, Ontario, Canada. The complete study area with all populations is shown in the largest map. The smaller maps show each study site zoomed in. Urban populations are depicted with blue triangles, and rural populations are shown with red circles. Maps were produced with Google Satellite imagery (2018)

Open in new tabDownload slide

We haphazardly collected three newly opened leaves and 10–15 seeds from each of 10 individuals per population, spaced at least 2 m apart. Each individual represents a single family, with leaves from the parent and the seeds from the offspring. Our aim was to capture population-level averages of outcrossing rates; therefore, we avoided collecting the smallest and largest fruits on a plant to avoid bias in sampling selfed versus potentially outcrossed fruits. Our expectation was that the fruits sampled were potentially derived from both chasmogamous and cleistogamous flowers. This sampling design does not allow us to distinguish between autogamous self-fertilization of cleistogamous flowers and geitonogamous self-fertilization of chasmogamous self-fertilization. Consequently, any signal of selfing that we observe could be the result of cleistogamous fruit production or geitonogamy. We sampled 1–3 ripe fruits per individual, which produced 1–5 seeds each. Immediately after collection, we dried leaf samples on silica and placed the seeds in a cooler with ice. At the end of each day, we transferred the seeds to a −20°C freezer until DNA extraction.

During previous visits to each population (2016–8), we measured the population area and estimated census population size (Table 1). We measured the length and width of each population to determine the area of the population. We censused population size by placing three 0.25 m² quadrats in the population and counting all individuals growing within each quadrat. We divided the population into thirds and haphazardly tossed the quadrat into each section to select the location of each quadrat in as unbiased a manner as possible. We estimated population size as the area of the population multiplied by four times the mean number of individuals growing within the quadrats.

DNA extraction and sequencing

We extracted DNA from 84 individuals (Table 1). We selected the samples from eight families from each population except YCP_urban where we were only able to collect from three individuals due to construction fencing blocking access to most of the site. We selected this number as a preliminary investigation into population genetic variation in I.capensis to balance investigative power with cost-efficiency. We extracted DNA from a single leaf and a single seed per family using the DNeasy Plant Mini Kit in tube format from QIAGEN^® (Germantown, MD, USA). Prior to extraction, we ground 15- to 20-mg of leaf tissue for 1.5 min at 25 Hz in a TissueLyser (QIAGEN), and after removing the seed coat, we ground each seed for 3 min at 25 Hz in a TissueLyser. We followed a modified protocol to extract DNA from the samples with lengthened centrifugation to promote pellet formation (Toczydlowski and Waller 2019). We evaluated the DNA quantity from each extraction using a Qubit™ Fluorometer dsDNA High Sensitivity Assay (Thermo Fisher Scientific, Waltham, MA, USA). We also ran all of the samples on a 1.5% agarose gel to assess the quality of the DNA.

We dried the extracted DNA for 1 hour in a PCR Thermocycler heated to 65°C and shipped our samples to the Elshire Group (Palmerston North, New Zealand) for library preparation and sequencing using genotype-by-sequencing (Elshire et al. 2011). Drying was performed without a plate seal to allow the samples to be shipped in a stable form without needing to stay cold (Elshire, pers. comm.). Samples were rehydrated and then digested with the restriction enzyme, ApeKI, tagged with combinatorial barcodes, and multiplexed into a single library. The samples were sequenced in a single lane on an Illumina HiSeq XTen (Illumina, Inc., San Diego, CA, USA) with 150 bp paired end reads. In total, we generated 465 million sequence reads for all 84 individuals that we sampled.

De novo assembly and single nucleotide polymorphism calling

We demultiplexed the raw sequences using Axe v0.3.3 (Murray and Borevitz 2018), then trimmed adaptors using TrimGalore! v0.6.6, a wrapper program for Cutadapt v1.14 (Martin 2011). Our initial quality filtering step was to discard reads with a Phred score less than 10 in sliding widows that were 15% of the read length, discard reads with uncalled bases and trim reads to a uniform length of 110 bp using the process_radtags.pl pipeline from Stacks v2.55 (Rochette, Rivera-Colón, and Catchen 2019). Because I. capensis does not currently have a reference genome, we used the denovo_map.pl pipeline from Stacks to assemble reads de novo and call single nucleotide polymorphisms (SNPs). We followed the protocol outlined in Paris, Stevens, and Catchen (2017) to optimize our assembly parameters. We identified the optimal values for the minimum depth of each stack of reads called (-m 5), the distance between stacks (-M 4), and the distance between catalogue loci (-n 4). Using these values, 99.3% of reads were assembled into loci with average depth per sample of 26.7×, (SD = 8.7×). Prior to the population filtering step below, we assembled 133 million paired-end reads into 1.13 million loci, with an average of 154 sites identified per locus.

Once the reads were assembled, we reran the Stacks populations script to create filtered SNP libraries (Rochette and Catchen 2017). We retained SNPs only if they met the following criteria: (i) to reduce gaps within contigs, we retained only loci that were found in 75% of individuals sampled (-r 0.75); (ii) to reduce linkage between individual SNPs, we called only a single random SNP per locus (–write-single-snp); (iii) to reduce the likelihood of calling bp errors, we included only SNPs with minor allele frequency greater than or equal to 5% (–maf 0.05). Following these filtering steps, we retained 5627 SNPs for downstream analyses. We used this complete (parents + offspring) dataset to estimate outcrossing rates. We then performed two identical runs of populations using only samples that were generated using only leaf tissue (parents, 42 samples) or seed tissue (offspring, 42 samples). Following this filtering step, we identified 13 363 polymorphic loci from which 4357 SNPs were retained in the parent-only dataset and 37 096 polymorphic loci from which 13 510 SNPs were retained in the offspring-only dataset. We used these datasets in our analyses on population genetic diversity and structure.

Data analysis

Outcrossing rate estimation

We used BORICE v3 (Colicchio et al. 2020) to estimate outcrossing rate of our sampled populations. BORICE applies Bayesian methods to individual genotype data to estimate the probability that an individual offspring is selfed or outcrossed (Colicchio et al. 2020). BORICE also estimates the inbreeding history (IH) of each parent, where IH = 0 for fully outbred parents, IH = 1 for a parent produced via self-fertilization, and IH > 2 for a selfing event more than one generation proceeding the parental generation (Colicchio et al. 2020). Importantly, BORICE provides accurate estimates of outcrossing and inbreeding even when family sizes are small, as is the case in our study where families consisted of a single parent and offspring (Koelling, Monnahan, and Kelly 2012).

We generated a VCF file that contained calls from 82 individuals. BORICE requires data from both parent and offspring to correctly infer outcrossing, so we removed two individuals that did not have a matching parent or offspring prior to analysis. We used custom python scripts to format the VCF file containing all individuals into the input files for BORICE. We also removed SNPs with impossible genotype calls (i.e. an offspring in which both alleles are different from the parent; N = 3,361 SNPs). We ran BORICE with increased burn-in (500) and chain length (2000 steps) to maximize replication for each observation.

Population genetic diversity

To identify an effect of urbanization on genetic diversity, we estimated population genetic diversity by generating individual-level estimates of observed heterozygosity (H_O), expected heterozygosity (H_E) under the assumption of Hardy–Weinberg equilibrium and the individual inbreeding coefficient (F), which quantifies the probability that an individual has two alleles at a particular locus that recently descended from a single common ancestor (i.e. via inbreeding). We generated PLINK v1.9 (Chang et al. 2015) compatible files, then used PLINK to filter out SNPs where >50% of genotype calls were missing, and calculate H_O, H_E and F. This step was performed separately on the parents and offspring datasets, and we retained 1365 SNPs from the parents, and 8417 SNPs from the offspring.

Habitat was a binary variable with two levels (urban and rural), and population size was a continuous variable that was log-transformed to reduce skew. We incorporated population as a random effect in the model to account for variation among populations due to differences unrelated to urbanization or demography. We ran separate models for parents and offspring, to better detect effects of biparental inbreeding in the offspring generation. We assessed the significance of fixed effects from both models using analysis of variance (ANOVA) implemented using the Anova function in the car v3.0-6 package (Fox and Weisburg 2011) to calculate Wald χ² test statistics with Type III sums-of-squares to test for significant interaction terms. Lastly, we used the DHARMa v0.4.1 package (Hartig 2021) to confirm that the residuals were normally distributed and demonstrated no autocorrelation or under-dispersion.

Population differentiation and genetic structure

To determine an effect of urbanization on population differentiation and structure, we generated estimates of genetic differentiation (F_ST) and population structure using the parental SNP dataset. We opted to use only the parents in these analyses to ensure that we were analysing only unrelated individuals. We estimated F_ST between populations with populations (–fstats) to calculate the pairwise fixation index analogue, φ_ST, using an analysis of molecular variance approach (Rochette, Rivera-Colón, and Catchen 2019).

We investigated the effects of urbanization on population genetic structure using a discriminant analysis of principal components (DAPC), an analysis that identifies and describes clusters of genetically related individuals (Jombart, Devillard, and Balloux 2010). We used the Plink-generated files generated for the population genetic diversity analyses for the structural analyses. We used the dapc function from adegenet v2.1.3 (Jombart and Ahmed 2011) to run a DAPC in R. We first identified the optimal number of clusters to use in our DAPC by examining the BIC scores of each cluster (K = 1–40). We found that K = 6 had the lowest BIC score, thus we proceeded to run the DAPC using the six sampled populations as the clustering variable in the analysis. Following optimization to avoid overfitting, we retained the first six PCs and all the discriminant functions in the DAPC. The final models explained 59% of the genetic variation between populations.

We complemented the results of the DAPC with an admixture analysis to further identify patterns of population clustering. We described population structure among all samples using ADMIXTURE v1.3.0 (Alexander, Novembre, and Lange 2009) implemented with the –cv flag. ADMIXTURE uses a likelihood model-based approach to estimate relatedness among individuals and allows the number of clusters (K) with the best predictive value to be determined using a cross-validation (CV) procedure (Alexander, Novembre, and Lange 2009). The value K with the lowest CV score indicates the optimal number of clusters in the dataset (Alexander, Novembre, and Lange 2009). We estimated ancestry iteratively for K = 1–6. We then used the R package ggplot2 v3.3.2 (Wickham 2011) to visualize the results from the ADMIXTURE analysis.

Results

Outcrossing rate

The majority of plants and populations exhibited outcrossing. Ninety-five per cent of offspring (39 out of 41 sampled) were outcrossed; two offspring sampled from KSR_rural were identified as being self-fertilized with 100% certainty (Fig. 1). Ninety-three per cent of parents (38 out of 41 sampled) originated from outcrossed matings; however, three parents were identified as being selfed. One parent was from the same family as one of the selfed offspring and was identified as being selfed with 60% certainty. The other parents were from ORT_rural and TCS_urban and were identified as being selfed with 100% certainty. A Kruskal–Wallis rank-sum test for unequal variance suggested that there was no effect of habitat on outcrossing (P = 0.193) or parental inbreeding history (P = 0.500).

Genetic diversity and differentiation

Population genetic diversity and differentiation differed among parents and offspring, and between populations. Mean H_O was 25% lower in parents than in offspring (0.16 and 0.20, respectively), H_E was the same in parents than offspring (0.33 for both), and F was 23% higher in parents than in offspring (0.52 and 0.41, respectively). However, the direction of difference in population genetic diversity varied among populations, with four of six populations having lower H_O and higher F in parents than offspring (Supplementary Table S1).

Population genetic diversity of offspring was significantly lower in urban habitats than in rural habitats (Supplementary Table S2). H_O was 16% lower (Habitat: χ²₁ = 6.23, P = 0.013) and F was 19% higher (Habitat: χ²₁ = 6.24, P = 0.013) in urban habitats relative to rural habitats. We also found a significant habitat × population size interaction on H_O (Habitat × Size: χ²₁ = 6.28, P = 0.012; Fig. 2) and F (Habitat × Size: χ²₁ = 6.32, P = 0.012; Fig. 2). Consistent with this trend, we found that urban census population sizes were 64% smaller than rural populations (approximately 7000 individuals vs 23 000 individuals, respectively; Table 1). There was no main effect of population size for any of our measures of genetic diversity. Lastly, we did not find an effect of habitat or population size on any measure of population genetic diversity in parents (Supplementary Table S2).

Figure 2:

The effects of habitat (urban: blue; rural: red) and population size (log transformed) on observed heterozygosity (H_O) and individual inbreeding coefficients (F) estimated for offspring (N = 42 individuals; 8417 SNPs) collected from six populations from urban and rural habitats in Greater Toronto Area. The x-axis has been log-transformed. Shading depicts the 95% confidence interval around each line

Open in new tabDownload slide

Population differentiation and structure

Population genetic differentiation was relatively high (mean F_ST = 0.18) and varied between populations (Supplementary Table S3), but did not differ between urban and rural habitats. We identified three genetic clusters of populations, and visualization of the DAPC results suggest that urban populations contain a subset of the genetic variation found in rural population. We found that WLN_rural was strongly separated from all other populations along the first axis (DPC1), which explained 76% of the variation between populations (Fig. 3A). We found that along the second axis (DPC2), ORT_rural and TCS_urban clustered together away from other populations, as did KSR_rural, YCP_urban and MUD_urban, which explained 18% of the variation between populations (Fig. 3A). The results of the DAPC were supported by the ADMIXTURE analysis, which identified K = 3 as the optimal number of clusters explaining ancestry among individuals (CV = 0.58; range: 0.58–0.82). Populations with shared ancestry were largely the same as those identified in the DAPC, although ORT_rural was admixed between the TCS_urban cluster and the KSR_rural, YCP_urban and MUD_urban (Fig. 3B).

Figure 3:

Population structure analyses on parents (N = 42 individuals; 1365 single nucleotide polymorphisms) sampled from urban and rural habitats in Greater Toronto Area. (A) Discriminant analysis of principal components highlighting population differentiation among three clusters of populations. Urban populations are represented by triangles and rural populations by circles. (B) Barplot of relatedness from ADMIXTURE analyses also identified three genetic clusters of populations. Each bar represents a single individual, and each colour represents a cluster of genetically related individuals that share ancestry

Open in new tabDownload slide

Discussion

Our study examined the effect of urbanization on the outcrossing rate and population genetic variation of I. capensis, a plant that is native to our study area. We found that individuals were primarily outcrossed and populations were genetically differentiated, but outcrossing rate and differentiation were largely independent of urbanization. However, we found that urban populations had lower genetic diversity and a stronger signature of past inbreeding relative to rural populations. Despite high levels of outcrossing, our results suggest that urban I. capensis populations have reduced genetic diversity and may be vulnerable to urbanization.

We did not find an effect of habitat on outcrossing rate. We found that most plants were outcrossed and those that were selfed were primarily located in rural habitats. These results are consistent with the only other study that examined reproductive output of urban I. capensis populations, which found no change in seed set or investment in self-fertilization (i.e. production of seeds produced from cleistogamous flowers) along an urbanization gradient (Barker and Sargent 2020). The degree of outcrossing we estimated is high but falls within the range reported from this species (Mitchell-Olds and Waller 1985; Waller and Knight 1989; Lu 2000; Steets, Knight, and Ashman 2007). Outcrossing rates are highly variable in I. capensis because the production of cleistogamous vs chasmogamous flowers is dependent on the surrounding environment and varies through the season. Plants always produce cleistogamous flowers; however, environmental factors such as herbivory (Steets and Ashman 2004; Steets, Knight, and Ashman 2007), intraspecific competition (Schmitt and Ehrhardt 1990), season (Schemske 1978) and light availability (Paoletti and Holsinger 1999) alter outcrossing rate by varying the number of chasmogamous flowers produced per individual. Consequently, the timing of sampling, demographics of the population or surrounding environment can influence the outcrossing rate within I. capensis populations.

We observed a strong degree of genetic differentiation resulting in three genetically distinct population clusters. Given the high levels of outcrossing, it is likely that mating occurs between individuals in the same population rather than between individuals of different populations. This pattern is consistent with pollination occurring via bees rather than hummingbirds, which disperse longer distances than bees (Thomas et al. 2008; Courter et al. 2013), and could lead to a greater degree of admixture than we observed in our study. We found that urban populations were genetically clustered with rural populations, indicating shared genetic variation among these populations. Because of the linear sampling pattern that we employed (Fig. 1), we would expect urban and rural populations to cluster separately if urbanization was influencing genetic differentiation, or if isolation-by-distance (IBD) was operating. However, we did not detect a clear signal of urbanization or of IBD, suggesting that other landscape features, such as rivers and watersheds (Toczydlowski and Waller 2019), are contributing to connectivity patterns in our I. capensis, potentially in combination with isolation-by-environment (e.g. urbanization) or IBD.

We found that urban populations had lower genetic diversity than rural populations in offspring. This result is consistent with other studies that identify urbanization as a mechanism that reduces the genetic diversity within populations (Miles et al. 2019; Schmidt et al. 2020). We also found an effect of habitat on inbreeding in offspring, but not in parents. These results may have arisen due to biparental inbreeding in the parental generation because of smaller urban populations, leading to lower genetic diversity (Ritland and Jain 1981). If reproduction occurred between close relatives rather than through selfing, this may also explain why we did not detect a signature of urbanization on outcrossing rates, and why there was no effect of habitat on parents. Alternatively, I. capensis seeds have a survival rate of 50% (Mitchell-Olds and Waller 1985), and it is possible that we sampled seeds that would not have survived to adulthood due to selection against progeny with high inbreeding loads. However, in such a situation we would expect to find inbreeding coefficients higher in offspring than in parents, which was not the case. Further tests on the viability of these seeds would be required to completely identify the realized influence of inbreeding load in these populations. Lastly, it is possible that we were able to detect an effect of habitat in offspring but not parents because we identified four times as many SNPs in the offspring, increasing the sensitivity of our analyses to population-specific variants.

We also found that habitat interacted with population size to affect genetic diversity, where heterozygosity increased, and inbreeding decreased, with increasing population sizes in urban habitats but not in rural habitats. The effect of population size on genetic diversity in urban habitats is consistent with predictions of how genetic drift operates in small populations (Kimura 1954); however, it is difficult to explain why we did not observe a similar effect in rural habitats. With only three populations in each habitat, we are cautious to over-interpret this trend and we suggest that studies aim to include more populations per habitat to clarify this effect. Lower genetic diversity may increase the vulnerability of small, native plant populations to urbanization and pose a challenge for conservation efforts in cities by reducing the amount of genetic material available for adaptation.

Although we included a small number of populations, our results inform conservation strategies by assessing the vulnerability of a native wildflower to urbanization. Our results on population genetic diversity indicate that urbanization may increase biparental inbreeding and strengthen the effects of genetic drift acting on populations, potentially reducing the adaptive capacity of urban populations. Despite the limitations of our study design, our findings suggest that cities would benefit from conservation strategies that aim to maximize the genetic diversity of native plant populations.

Supplementary data

Supplementary data are available at JUECOL online.

Acknowledgements

We thank John Kelly, who generously provided custom python scripts and valuable advice on analyses with BORICE, and Lindsay Miles for helpful comments on the manuscript. We also thank Zahra Sina, Samreen Munim, and Vanessa Nhan for help with the population censuses. LRR was funded by NSERC CGS-D and Queen Elizabeth II scholarships, and MTTJ was funded by an NSERC Discovery Grant, an NSERC E.W.R. Steacie Fellowship, and a Canada Research Chair. LRR collected the data, extracted DNA, analysed the data, and wrote the manuscript with input from MTJJ at all stages. The Elshire Group prepared the genomic libraries and sequenced the DNA.

Data availability

All data, sequences and code will be archived in NCBI repository upon acceptance of this manuscript.

Conflict of interest statement. None declared.

References