DOI: 10.5061/dryad.v6wwpzh1m
title: Genotype, phenotype and linkage data for Mimulus parishii x M. cardinalis hybrid incompatibility study

description: The evolution of genomic incompatibilities causing postzygotic barriers to
 hybridization is a key step in species divergence. Incompatibilities take
 two general forms – structural divergence between chromosomes leading to
 severe hybrid sterility in F1 hybrids and epistatic interactions between
 genes causing reduced fitness of hybrid gametes or zygotes
 (Dobzhansky-Muller incompatibilities). Despite substantial recent progress
 in understanding the molecular mechanisms and evolutionary origins of both
 types of incompatibility, how each behaves across multiple generations of
 hybridization remains relatively unexplored. Here, we use genetic mapping
 in F2 and RIL hybrid populations between the phenotypically divergent but
 naturally hybridizing monkeyflowers Mimulus cardinalis and M. parishii to
 characterize the genetic basis of hybrid incompatibility and examine its
 changing effects over multiple generations of experimental hybridization.
 In F2s, we found severe hybrid pollen inviability (&lt; 50% reduction
 vs. parental genotypes) and pseudolinkage caused by a reciprocal
 translocation between Chromosomes 6 and 7 in the parental species. RILs
 retained excess heterozygosity around the translocation breakpoints, which
 caused substantial pollen inviability when interstitial crossovers had not
 created compatible heterokaryotypic configurations. Strong transmission
 ratio distortion and inter-chromosomal linkage disequilibrium in both F2s
 and RILs identified a novel two-locus genic incompatibility causing
 sex-independent gametophytic (haploid) lethality. The latter interaction
 eliminated three of the expected nine F2 genotypic classes via F1 gamete
 loss without detectable effects on the pollen number or viability
 of F2 double heterozygotes. Along with the mapping of numerous
 milder incompatibilities, these key findings illuminate the complex
 genetics of plant hybrid breakdown and are an important step toward
 understanding the genomic consequences of natural hybridization in this
 model system.
description: Study system and
 plant lines The plants in this
 study were all derived from two highly (&gt;10 generations) inbred
 lines of Sierran <em>M. cardinalis</em> (CE10) and
 <em>M. parishii</em> (PAR), which were also used in previous
 investigations of species barriers (Bradshaw <em>et al.</em>
 1998; Schemske and Bradshaw 1999; Ramsey <em>et al.</em> 2003;
 Bradshaw and Schemske 2003; Fishman <em>et al.</em> 2013,
 2015; Nelson <em>et al.</em> 2021a). We generated
 PAR x CE10 F<sub>1</sub> hybrids by hand-pollination (with
 prior emasculation of the PAR seed parent in the bud) and
 F<sub>2</sub> hybrids by self-pollination of
 F<sub>1</sub> hybrids. The F<sub>2</sub> hybrids
 were grown in two separate greenhouse common gardens at the University of
 Montana (UM-F<sub>2</sub>; total <em>N</em> = 524)
 and the University of Connecticut (UC-F<sub>2</sub>
 <em>N</em> = 253), along with parental control lines, and were
 phenotyped for numerous floral and vegetative traits including the pollen
 fertility traits presented here. Recombinant inbred lines (RILs) were
 generated by single-seed-descent from additional F<sub>2</sub>
 individuals grown at the University of Georgia and California State
 Polytechnic University, Pomona; a total of 167 RILs were formed through
 3–6 generations of self-fertilization. DNA extraction and
 sequencing Genomic DNA was
 extracted from bud and leaf tissue of the greenhouse-grown
 F<sub>2</sub> and RIL mapping populations using a CTAB
 chloroform protocol modified for 96-well plates (dx.doi.org/10.17504/protocols.io.bgv6jw9e). We used a double-digest restriction-site associated DNA sequencing (ddRADSeq) protocol to generate genome-wide sequence clusters (tags), following the BestRAD library preparation protocol (dx.doi.org/10.17504/protocols.io.6awhafe), using restriction enzymes PstI and BfaI (New England Biolabs, Ipswich, MA). Post-digestion, half plates of individual DNAs were labeled by ligation of 48 unique in-line barcoded adapters, then pooled for size selection. Libraries were prepared using NEBNext Ultra II library preparation kits for Illumina (New England BioLabs, Ipswich, MA). Each pool was indexed with a unique NEBNext i7 adapter and an i5 adapter containing a degenerate barcode and PCR amplified with 12 cycles. The F<sub>2</sub> libraries were size-selected to 200-700bp using BluePippin 2% agarose cassettes (Sage Science, Beverly, MA) and sequenced (150-bp paired-end reads) in a partial lane of an Illumina HiSeq4000 sequencer at GC3F, the University of Oregon Genomics Core Facility. The RIL library was sequenced (150-bp paired-ends) without size-selection on an Illumina HiSeq4000 at Genewiz (South Plainfield, NJ). Sequence processing and linkage mapping After sequencing, two separate ddRAD datasets were analyzed: one with samples from both F<sub>2</sub> populations (<em>N</em> = 283 UM-F<sub>2</sub> hybrids with 3 <em>M.</em> <em>parishii</em> and 2 <em>M. cardinalis </em>controls, and 253 UC-F<sub>2</sub> hybrids, with 3 each F<sub>1</sub>, <em>M. parishii </em>and <em>M. cardinalis </em>controls) and one with samples from the RIL population (<em>N</em> = 167). Samples from both datasets were demultiplexed using a custom Python script (dx.doi.org/10.17504/protocols.io.bjnbkman), trimmed using Trimmomatic (Bolger <em>et al.</em> 2014), mapped to the <em>M. cardinalis</em> CE10 v2.0 reference genome (http://mimubase.org/FTP/Genomes/CE10g_v2.0) using BWA MEM, and indexed using SAMtools (Li <em>et al.</em> 2009). The RIL dataset was also filtered in SAMtools using a mapping quality ≥ 29. We called SNPs in both datasets using HaplotypeCaller in GATK v3.3 in F<sub>2</sub>s, v4.1.8.1 in RILs (McKenna <em>et al.</em> 2010). Next, we performed a series of filtering steps to generate sets of high-quality SNPs. In the F<sub>2</sub> dataset, we filtered using vcftools (Danecek <em>et al.</em> 2011), retaining sites with read depth ≥ 5, mapping quality ≥ 10, and &lt; 40% missing data. We also filtered out loci deviating from Hardy-Weinberg Equilibrium at a p-value &lt; 0.00005. In the RIL dataset, we filtered a combined GVCF file using GATK, retaining sites with read depth ≥ 4*<em>N</em> (with <em>N </em>= number of RIL samples) and &lt; the mean + 2 * the standard deviation, QD score &lt; 2.0, FS score &gt; 60, MQ &lt; 40, MQ rank sum &lt; -12.5, ReadPosRankSum &lt; -8.0, and &lt; 10% missing genotypes. For both datasets, we used custom scripts to remove sites that were not polymorphic in the parents and heterozygous in the F<sub>1</sub> hybrids (F2: https://github.com/bergcolette/F2_genotype_processing). We excluded individuals from the F<sub>2</sub> dataset with &gt; 10% missing data and from the RIL dataset with low coverage, high missingness, or excessive heterozygosity (&gt; 50%, indicating line contamination). These filtering steps produced an F<sub>2</sub> dataset with 18,119 SNPs (<em>N</em> = 252 UM-F<sub>2</sub> and 253 UC-F<sub>2</sub>) and a RIL dataset with 47,851 SNPs (<em>N</em> = 145). To produce sets of high-quality marker genotypes for mapping, we binned each dataset into 18-SNP windows using custom Python and R scripts (provided at GitHub links above), requiring ≥ 8 sites to have SNP genotype calls to assign a windowed genotype. In the F<sub>2</sub> binning script, <em>M.</em> <em>cardinalis </em>homozygotes were coded as 2, <em>M.</em> <em>parishii </em>homozygotes as 0, and heterozygotes as 1. We called windows with mean values &lt; 0.2 as <em>parishii </em>homozygotes, &gt; 1.8 as <em>cardinalis </em>homozygotes, and between 0.8 and 1.2 as heterozygotes. Windows with means outside of these ranges were coded as missing genotypes. For the RILs, we required ≥ 88% of SNP calls to match each other to assign each parental homozygous genotype (<em>e.g</em>., 16/18 sites = homozygous for <em>M. cardinalis</em> alleles; https://github.com/vasotola/GenomicsScripts). We generated linkage maps for each dataset using Lep-MAP3 (Rastas 2017). First, we used the <em>SeparateChromosomes2</em> module to assign markers to linkage groups (F<sub>2</sub>: LodLimit = 25, theta = 3, RIL: LodLimit = 28, theta = 0.2). In the RIL dataset, 10 markers were assigned to linkage groups inconsistent with the reference genome assembly; we manually re-assigned these markers to linkage groups corresponding to their reference assembly chromosomes. Next, we performed iterative ordering using the <em>OrderMarkers2</em> module (Kosambi mapping function; 6 iterations/per linkage group in the F<sub>2</sub>s, 10 in the RILs); the order with the highest likelihood for each linkage group was chosen. This resulted in an F<sub>2</sub> map with 997 markers in seven linkage groups and a RIL map with 2,535 markers in eight linkage groups. In the RIL dataset, the genotype matrix output by Lep-MAP3 differed in two important respects from the input file. First, due to stringent thresholds for calling windowed genotypes, our input file includes a high percentage of missing data (23% of genotypes are coded as ‘no call’), whereas the output file contains no missing data (Lep-MAP3 converts each ‘no call’ genotype to a called genotype). Second, the Lep-MAP3 output file contains more heterozygous genotype calls than the input file. The reason for this increase in heterozygosity is that Lep-MAP3 disproportionately converts ‘no call’ genotypes to heterozygotes: relative to the input file, the output genotype matrix includes 115% more heterozygotes, compared to only 18% more <em>M. cardinalis</em> homozygotes and 20% more <em>M. parishii</em> homozygotes. Notably, Lep-MAP3 frequently converted ‘no call’ genotypes to heterozygotes when they occur at single markers between recombination breakpoints. Because most recombinational switches in this RIL population are between alternative homozygotes, any window that contains an actual breakpoint will carry a mixture of <em>M. cardinalis</em> and <em>M. parishii</em> homozygotes at the 18 SNPs (and thus be coded as ‘no call’ in our windowed genotype matrix). To circumvent these problems, for all downstream analyses, we used a modified version of the genotype matrix output from Lep-MAP3 in which genotypes were recoded as ‘no call’ as in the input file.    QTL mapping of pollen traits In the UM-F<sub>2</sub> and RIL populations, we directly assessed male fertility by collecting all four anthers of the first flower from each plant into 50 ml of lactophenol-aniline blue dye. We counted viable (darkly stained) and inviable (unstained) pollen grains using a hemocytometer (≥ 100 grains/flower). We estimated total pollen grains per flower (count per ml x 50) and pollen viability as viable grains/total counted. For a handful of individuals with &lt; 100 pollen grains counted, pollen viability was scored as missing data. Phenotyping in the RIL population was performed on siblings or selfed progeny of the individuals genotyped. In a few cases in the RILs (<em>N </em>= 12), the first flower had &lt; 100 pollen grains, so both traits were instead collected from the second flower. We mapped pollen QTLs in Windows QTL Cartographer 2.5 (Wang et al. 2012) using composite interval mapping (CIM; Zeng 1993, 1994), with forward-backward stepwise regression, a window size of 10 cM, five background markers, and a 1-cM walk speed. We used permutations (<em>N</em> = 1000) to set genome-wide significance thresholds for QTL peaks and calculated 1.5-LOD-drops to determine confidence intervals for QTL locations. Because pollen viability QTLs exhibited complex interactions in F<sub>2</sub>s (see Results), we directly estimated QTL effects (at each peak marker) and interactions using the Generalized Linear Model module in JMP16 (SAS Institute, Cary, NC), starting with a full factorial model of all QTLs and 2-way interactions and removing nonsignificant (P &gt; 0.05) interactions. To test for pollen viability signatures of a two-locus gametic incompatibility (LG4-LG8) detected from TRD and LD patterns, we contrasted Least Squared Means for each extant F<sub>2</sub> class from an ANOVA also including the large-effect translocation breakpoint marker at 60.55Mb on Chr 6. QTL mapping of underdominant effects in RILs has low power due to (relatively) few heterozygotes, plus most RIL QTL-mapping algorithms (including those in WinQTLCart) exclude non-homozygous genotypes. To directly test for persistent underdominant effects of the Chr 6-7 translocation in RILs, we screened for genotype-pollen viability associations at uniquely positioned markers with a joint sample size &gt;70 across the breakpoint region (57.07-61.08 Mb of Chr6 and 7.51- 13.05 of Chr 7; n = 42), using t-tests in the response screening module in JMP16 (SAS Institute, Cary, NC). However, we caution that because RIL hybrid male sterility was measured on siblings or descendants of the genotyped individuals, underdominant effects will likely be underestimated. That is, regions genotyped as heterozygous might actually be homozygous in phenotyped individuals, potentially explaining why a few of them are highly fertile. We conservatively controlled for multiple tests by using an FDR-corrected <em>P</em>-value of 0.05. 
description: The data are supplied as a multi-sheet Excel workbook. They
 should be openable by text handlers such as BBedit or similar.