Chromosomal evolution, environmental heterogeneity, and migration drive spatial patterns of species richness in Calochortus (Liliaceae)

We used nuclear genomic data and statistical models to evaluate the ecological and evolutionary processes shaping spatial variation in species richness in Calochortus (Liliaceae, 74 spp.). Calochortus occupies diverse habitats in the western United States and Mexico and has a center of diversity in the California Floristic Province, marked by multiple orogenies, winter rainfall, and highly divergent climates and substrates (including serpentine). We used sequences of 294 low-copy nuclear loci to produce a time-calibrated phylogeny, estimate historical biogeography, and test hypotheses regarding drivers of present-day spatial patterns in species number. Speciation and species coexistence require reproductive isolation and ecological divergence, so we examined the roles of chromosome number, environmental heterogeneity, and migration in shaping local species richness. Six major clades – inhabiting different geographic/climatic areas, and often marked by different base chromosome numbers (n=6-10) – began diverging from each other ~10.3 million years ago. As predicted, local species number increased significantly with local heterogeneity in chromosome number, elevation, soil characteristics, and serpentine presence. Species richness is greatest in the Transverse/Peninsular Ranges where clades with different chromosome numbers overlap, topographic complexity provides diverse conditions over short distances, and several physiographic provinces meet allowing immigration by several clades. Recently diverged sister-species pairs generally have peripatric distributions, and maximum geographic overlap between species increases over the first million years since divergence, suggesting that chromosomal evolution, genetic divergence leading to gametic isolation or hybrid inviability/sterility, and/or ecological divergence over small spatial scales may permit species co-occurrence.

Methods
We used custom-designed baits for anchored hybrid enrichment, preparing Illumina sequence libraries for 158 samples (156 Calochortus and 2 outgroups). Our final dataset consisted of 294 low-copy nuclear loci (provided here as alignments). Specific methods for generating these data are also follows:

Sampling. We included 1 to 2 samples per species and subspecies; herbarium vouchers for new samples were deposited in the herbaria noted in SI Appendix, Table S9. Total genomic DNAs were extracted from silica-dried leaf or floral tis- sue using DNeasy plant kits (Qiagen, Valencia CA) following the manufacturer’s instructions. We included all extant Calochortus species except extremely rare C. rustvoldii [known from only two sites/pixels in the western Transverse Ranges (76))]; including it would little affect our analyses.

Library Preparation. We prepared Illumina sequence libraries for 158 samples (156 Calochortus and 2 outgroups) following Lemmon et al. (77–79). We used a Covaris ultrasonicator (with reduced time for degraded samples) to fragment DNA to 140 to 400 bp; performed end-repair and A-tailing, ligated common Illumina adapters onto the template DNA ends using a Beckman Coulter FXp liquid-handling robot, and performed indexing PCR.

AHE Probe Design. Following Hamilton et al. (80) and Banker et al. (81), we developed hybrid enrichment probes for Lobelia and Lilium. We mapped sequences from two assembled transcriptomes—Lobelia siphilitica (from E. Carpenter) and Lilium superbum (from J. Leebens-Mack and C. dePamphi- lis)—to probe sequences from 27 references of the Angiosperm V1 AHE design (82, 83). Mapped transcriptome sequences were aligned to Angiosperm V1 reference sequences using MAFFT v7.023b (84). We used Geneious R9 (85) to visually inspect alignments, remove transcriptome sequences that were not clearly homologous, and trim transcriptome sequences to exons represented by the 27 reference sequences. Probes were tiled uniformly at 2.7× density. This probe set (Angiosperm V2 AHE) design contains 29 references (57,471 probes), including transcriptomes representing Lobelia and Lilium.

To improve phylogenetic resolution within Calochortus, we utilized data from two more species and expanded the targets into regions flanking the exons in Angiosperm V2 AHE. We first collected whole-genome sequence data (Illumina paired-end 200-bp protocol) for Calochortus albus (160M reads) and C. flexuo- sus (536M reads). We then used methods and scripts from Banker et al. (83) to identify loci and design probes. After merging overlapping reads following (86), we mapped merged reads to probe sequences from Angiosperm V2 AHE using Liliaceae as a reference. After extending consensus sequences into flanking regions using iterative mapping (80, 81), we aligned by locus the resulting con- sensus sequences to the reference, using Geneious to visually inspect alignments and trim poorly aligned regions from alignment ends. We removed 27 of the 517 alignments to ensure target loci did not overlap. After masking repetitive regions (80), final alignments contained 375,527 sites. Tiling 120 bp probes at 4.5× density for both references produced 16,600 probes. We used this probe set (AHE Cal1) to produce an Agilent Technologies Custom SureSelect XT kit for hybrid DNA enrichment.

Library Enrichment and Sequencing. We pooled indexed libraries in groups of 16 before enriching with the probe kit just described. Before sequencing, we pooled enriched libraries and assessed quality by Kapa qPCR (Roche). We sequenced samples on an Illumina NovaSeq6000 at Florida State with a PE-150 bp protocol and 8-bp (dual) indexing. After filtering poor-quality reads with the Illumina CASAVA v1.8 high-chastity filter and demultiplexing, we obtained an average of 5.4M read pairs per sample (~1.6 Gb).

Nuclear Assemblies and Alignment. We processed reads, corrected sequenc- ing errors, and trimmed adaptors following (80, 81, 86). Assembly used a quasi-de novo approach, mapping reads to probe-region sequences for both Calochortus species in the Cal1 design. We used the resulting consensus sequences with ≥98× coverage to determine orthology (80) and then formed orthologous clus- ters and removed clusters containing <50% of the individuals, resulting in 294 loci analyzed. For each locus, we aligned sequences in Mafft v7.023 and then trimmed/masked alignments with an automated procedure (80) using default parameters (i.e., MINGOODSITES=14, PROPGOOD=0.5). We verified alignment quality by inspecting the alignments in Geneious. Upon manuscript acceptance, accession codes will be provided for raw reads in the Sequence Read Archive and sequences and alignments in GenBank.