https://orcid.org/0000-0001-8382-2533
Abstract
Habitat-specific thermal responses are well documented in various organisms and likely determine the vulnerability of populations to climate change. However, the underlying roles of genetics and plasticity that shape such habitat-specific patterns are rarely investigated together. Here we examined the thermal plasticity of the butterfly Bicyclus dorothea originating from rainforest and ecotone habitats in Cameroon under common garden conditions. We also sampled wild-caught butterflies from forest and ecotone sites and used RADseq to explore genome-wide population differentiation. We found differences in the level of phenotypic plasticity across habitats. Specifically, ecotone populations exhibited greater sensitivity in wing eyespot features with variable development temperatures relative to rainforest populations. Known adaptive roles of wing eyespots in Bicyclus species suggest that this morphological plasticity is likely under divergent selection across environmental gradients. However, we found no distinct population structure of genome-wide variation between habitats, suggesting high level of ongoing gene flow between habitats is homogenizing most parts of the genome.
Introduction
Phenotypic plasticity, particularly polyphenism (Mayr, 1963), refers to the phenomenon that multiple discrete phenotypes could arise from the exact same genotype in response to different environmental cues (Scheiner, 1993). Phenotypic plasticity is commonly observed in diverse organisms (Pigliucci & Pigliucci, 2001) and has been suggested to be one of the key mechanisms that allow organisms to quickly respond to short-term changes in the environment, especially within the span of an individual’s lifetime. The degree of phenotypic plasticity could be heritable, and if there is genetic variation in plasticity, it could also be the target of selection (Via et al., 1995), as shown in both laboratory research (Newman, 1994) and wild populations (Brommer et al., 2005; Nussey et al., 2005). For many species, populations inhabit heterogenous environments and distinct habitats to which they are subjected to differential selection pressures. This could lead to habitat-specific variation in trait plasticity for populations, with implications for responses to climate change impacts (Kellermann & Heerwaarden, 2019; Landry Yuan et al., 2018; Valladares et al., 2014).
Bicyclus butterflies represent one of the most studied model systems for the ecology and evolution of plasticity (Brakefield & Reitsma, 1991; Brakefield & Schemske, 2010; Brakefield et al., 2009; Windig et al., 1994). Species in the Bicyclus genus all share characteristic eyespots on their wings. They are endemic to Africa occupying different types of seasonal and aseasonal habitats such as savanna, rainforest, and mountains (Aduse-Poku et al., 2022; Condamin, 1973; Halali et al., 2021; Roskam & Brakefield, 1999). For species inhabiting seasonal habitats, the size and pattern of eyespots (among other traits) exhibit a striking discontinuous phenotypic plasticity—also known as seasonal polyphenism—over wet and dry seasons in their native habitats; such seasonal variations in morphology have been studied in B. anynana, B. cottrelli, B. vansoni, B. ena, and B. safitza in Malawi (Brakefield & Reitsma, 1991; van Bergen et al., 2017; Windig et al., 1994) and recently in B. dorothea in Cameroon (Dongmo et al., 2018). In general, these species exhibit wings with reduced eyespot size during the dry season and more prominent eyespots in the wet season in East and Southern Africa (Brakefield & Reitsma, 1991). Plasticity in the eyespots of wings of B. anynana is known to be adaptive. During the wet season, conspicuous eyespots deflect the attacks from predators, while dry season patterns aid in crypsis (Lyytinen et al., 2004; Prudic et al., 2015). Under natural conditions, seasonal variation in wing pattern is triggered by different climate cues including temperature and rainfall (Roskam & Brakefield, 1999). These cues themselves are translated into hormone signaling that will further regulate phenotypic expression depending on the intensity of the cues (Oostra et al., 2011). However, several studies have also shown that temperature variation alone, under laboratory conditions, can induce important variation in wing patterns—lower temperatures tend to yield individuals with reduced eyespots, while intermediate and high temperatures produce individuals with large eyespots (Nokelainen et al., 2018; Oostra et al., 2014; Roskam & Brakefield, 1996). Moreover, other life-history traits in Bicyclus can respond plastically to temperature gradients, indicating adaptation capacity to novel environments (van Bergen et al., 2017).
Bicyclus dorothea is a tropical butterfly widely distributed across diverse habitats in West and Central Africa. Like most Bicyclus species, phenotypic plasticity in the wing patterns of B. dorothea, especially eyespot size, has been observed in the wild between wet and dry seasons. Dongmo et al. (2018) found that the degree of plasticity is a function of habitat type under natural conditions. Specifically, populations from ecotone habitats are more variable over the wet and dry seasons relative to their forest counterparts. However, environmental cues for development are different between habitats. Whether the observed difference in the degree of phenotypic plasticity solely depends on the differences in environmental cues during ontogeny between forest and ecotone habitats or is intrinsic due to the difference in genomic makeups that give rise to different abilities of plasticity (Lafuente & Beldade, 2019; Suzuki & Nijhout, 2006) among these populations is still unclear.
To address this question, we collected B. dorothea individuals from both lowland rainforest and the rainforest–savanna ecotone and reared them in the laboratory under common garden conditions to obtain second-generation offspring, controlling for maternal effects. We used the second-generation offspring to quantify and compare levels of wing-pattern plasticity of individuals from different habitat types developing at the same set of temperatures. We also used Restriction site Associated DNA (RAD) sequencing to examine whether there is genome-wide differentiation between populations from forest and ecotone habitats and identify loci that are divergent between habitats that may explain the observed differences in phenotypic plasticity. While the environmental gradient between rainforest and ecotone has been the subject of long-term studies of divergence and speciation in vertebrates (Freedman et al., 2010, 2023; Smith et al., 2005; Zhen et al., 2017), this study is the first to examine plasticity and genomic differentiation together across rainforest and ecotone gradients. By exploring these mechanisms, the results have relevance to the origin of diversity in tropical rainforest (Smith et al., 1997) and climate change responses (Landry Yuan et al., 2018).
Methods
Study organism
Bicyclus dorothea is a widespread Bicyclus species whose distribution extends from Southern Cameroon in Central Africa across much of West Africa (Dongmo et al., 2017). Bicyclus dorothea inhabits diverse forest habitats across its range but appears to be especially attracted to disturbed open areas within forests (Dongmo et al., 2017). Grasses (Poaceae) serve as the larval host plants for B. dorothea, and it has been reared successfully on Pennisetum glaucum and Axonopus compressus (Dongmo et al., 2017).
Cameroon extends across 11 degrees of latitude (2°N–13°N), creating an ecological gradient with a range of diverse tropical habitats (Hansen et al., 2002). Southern Cameroon is dominated by rainforest habitats that give way to ecotone habitats, a mosaic of forest and savannah habitats, in the central part of the country (Figure 1). Bicyclus dorothea are readily found in abundance in both forest and ecotone sites throughout Cameroon (Dongmo et al., 2018). For the phenotypic assays and lab rearing, we focused on four populations, two from lowland rainforest (Mbalmayo and Somalomo) and two from ecotone (Ako and Ndikiniméki). To assess the variation of environmental temperatures in the four localities selected for the phenotypic data, we used microclimate simulations from Kearney et al., (2014) at 1 cm from the ground under two shade conditions (0% and 100%) following Dongmo et al. (2021). We downloaded long-term estimates of hourly minimum and maximum temperatures (~15 km resolution) of each locality (Hijmans et al., 2005; Kearney et al., 2014). Twenty-four layers (for each hour of the day) were downloaded, and temperature values were extracted for each lat/lon of the localities using the packages raster (Hijmans et al., 2005) and ncdf4 (Pierce, 2019). Mean monthly values of minimum and maximum temperatures under the two shade conditions were plotted (Supplementary Figure S1). To understand patterns of population structure across the range of B. dorothea in Cameroon, we sequenced samples collected from 18 sites, that is, 4 from ecotone and 14 from forest, employing RADseq (Figure 1; Supplementary Table S1).
Sampling locations for populations in Cameroon, Central Africa. White outlines of sites indicate those populations that were also used in phenotypic plasticity assay. Background forest cover was produced using a MODIS-derived data that estimates the percent of tree canopy cover across the study region (Hansen et al., 2002).
Sampling and lab rearing for phenotypic measurements
Adult B. dorothea individuals were collected during the wet season in the four localities for lab rearing (Figure 1). The number of field trips and adult individuals captured for laboratory cultures varied across localities: Ako (one field trip, N = 52 females), Ndikiniméki (three field trips, N = 8, 19, and 12 females, respectively), Mbalmayo (two field trips, N = 32 and 36 females, respectively), and Somalomo (two field trips, N = 19 and 29 females, respectively). All collected females were transported to the laboratory of the International Institute of Tropical Agriculture in Yaoundé, Cameroon. Once there, individuals were held in cages (24 × 24 × 24 cm; no more than 10 females per cage) made of white polyester mesh (to facilitate observation of insect activities and also to allow good ventilation) containing mashed banana, distilled water-soaked cotton (Brakefield & Frankino, 2009). The cages were kept in a room maintained at 26°C, 75%–80% relative humidity, and photoperiod D12:L12. In each cage, young potted millet (Pennisetum glaucum) plants were exposed to female butterflies as the egg-laying host. Eggs were then collected and kept in Petri dishes lined with a black moistened filter paper for eclosion. New hatchlings were transferred on the natural potted lawn of Axonopus compressus and reared through the first generation to avoid maternal effects (Jenkins & Hoffmann, 1994; Mousseau & Dingle, 1991; Mousseau & Fox, 1998). These hatchlings were maintained in the lab with respect to habitat, sampling localities, and field trip at each sampling locality. Adults obtained from the first generation were allowed to mate in cages, and the produced eggs (second generation) were allocated to unique thermal environments.
Second-generation eggs, collected from females of all populations from the first generation, were placed in incubators set at three different temperature regimes (22°C, 26°C, and 30°C), with relative humidity ranging between 70% and 80% and D12:L12 photoperiod. For each population, at least three batches of 100 eggs were lined on moist black filter papers in Petri dishes and exposed to a climate cabinet (Percival, model I-36VL) set at a specific temperature (i.e., a total of 300 eggs per temperature). Newly hatched larvae were counted and immediately transferred to potted lawn grass placed in a climate cabinet set at the given temperature. These hatchlings were allowed to develop to the adult stage in each climate cabinet. Host plants were monitored daily and watered or replaced as necessary.
Measurement and analysis of wing-pattern variation
Adult butterflies from the crossings at each temperature were frozen at −40°C, and their wings were separated from the body to assess morphological parameters. Bicyclus dorothea eyespots are similar in structure and position to most Bicyclus species (Dongmo et al., 2017); however, they are located on the ventral side of the wings and are not visible from the dorsal side as they are in many other Bicyclus species. There are nine eyespots in total, two on the forewings and seven on the hindwings (Dongmo et al., 2018; Supplementary Figure S2), though one can find three spots on the forewings and eight on the hindwings of some individuals. Eyespots are approximately circular in shape and made of numerous rings: a central white focus, a middle black disc, and a golden outer disc (Mateus et al., 2014; Monteiro et al., 1997). Following Dongmo et al. (2018), we measured nine wing-pattern characters and four proxies for the wing size (Zhen et al., 2023). As wing eyespots are not perfectly circular, all measurements were done parallel to the vein of the wing as in Dongmo et al. (2018) (see also Supplementary Figure S2). Moreover, measurements were taken only on the left fore and hind wings. To measure the diameters of the eyespots, the fore and hind wings were placed between two glass slides under a stereomicroscope fitted with a micrometer eyepiece at 6× magnification. Wing length, from the thorax to the apex, and the width of the fore and hind wings were measured using a caliper (0.1-mm accuracy).
All 13 wing-pattern measurements (Supplementary Table S2) were reduced using a principal component analysis (PCA). We performed the PCA using the “prcomp” function and the correlation matrix of the R base “stats” package in R version 3.5.1 (R Core Team, 2019) on the combined data from ecotone and forest habitats. To test the significance of the standardized loadings of each variable with the first two principal components (PC1 and PC2), we applied a nonparametric permutation test on the resulted PCA with 999 permutations implemented with the function “permut_pc_test” of the package “syndRomics” (Torres-Espín et al., 2021). Furthermore, we used a linear mixed model to analyze the effect of rearing temperature (considered here as a categorical variable), habitat, and sex (fixed factors), including two- and three-way interactions, on both PC1 and PC2 with field trip number as random effect. All models were computed using the function “lmer” in the lme4 package (Bates et al., 2015) in the R environment version 3.5.1 (R Core Team, 2019). Due to the nature of the development of butterflies (together in open cages), we were unable to track matrilines and instead use here “field trip” as random factor to account for batch effects. Post hoc pairwise comparison (Tukey–Kramer Honestly Significant Differences) tests were implemented using the package “emmeans” (Lenth et al., 2021) to compare PCs means between habitats, sites, sexes, and rearing temperatures. Additionally, we examined eyespot size relative to wing size across temperatures and treatments following the approach of Mateus & Beldade (2022) and compared the PCA results.
Sampling and DNA extraction for RADseq
For RAD sequencing, a total of 192 B. dorothea samples were collected from 18 geographically distant sites in Cameroon, including the four sampling sites in phenotypic assay (Figure 1; Supplementary Table S1). Sampling sites were classified into ecotone or forest habitat type by researchers in the field and had previously been confirmed using remote sensing data (Slabbekoorn & Smith, 2002; Smith et al., 2005, 2013), which included 4 ecotone populations and 14 forest populations. Each population had 3–28 samples (Supplementary Table S1) and low-quality samples were filtered in subsequent steps, resulting in a total of 156 samples with good quality. For analysis using population allele frequencies, we used only the 12 populations that had no less than eight samples. Samples were stored in Queens Lysis Buffer at -20ºC for approximately six months before DNA extraction (Smith et al., 1997, 2005). Genomic DNA was extracted using DNeasy kit (Qiagen, CA) using thorax and abdomen tissues.
RAD library preparation and sequencing
RAD libraries were prepared using a modified version of the protocol outlined in Ali et al. (2016). Briefly, genomic DNA was arrayed into 96-well plates, digested with PstI, and ligated to biotinylated adapters with a unique barcode for each well. The resulting samples were pooled using an equal volume (5 µL; approximately 1/3 of the total volume) from each well. The pooled sample was sonicated and biotinylated fragments (a.k.a. “RAD tags”) were purified using streptavidin-coated magnetic beads. An initial Illumina sequencing library was generated using the purified RAD tags as input to the NEBNext Ultra DNA Library Prep Kit. The resulting library was sequenced with paired-end 100 bp reads in one lane on an Illumina HiSeq 4000. After sequencing, individual samples were demultiplexed by requiring a perfect match to the barcode and partial restriction site, and the number of reads from each sample was determined (Ali et al., 2016). To standardize the amount of sequencing data generated from each sample, the remaining digested/ligated samples were pooled again using volumes based on the read counts from the initial sequencing (Prince et al., 2017). Additional RAD libraries were generated, sequenced, and demultiplexed as described earlier using the repooled samples. Sequences from the two batches were combined for each sample.
De novo RAD locus assembly and SNP discovery
RAD loci were discovered and assembled as previously described (Miller et al., 2012; Sağlam et al., 2016) using custom PERL scripts (Miller et al., 2012), the alignment program Novoalign and the genome assembler PRICE (Ruby et al., 2013). This pipeline resulted in 61,988 RAD contigs ranging from 80 to 953 bp (Zhen et al., 2023). Paired-end sequence data were aligned to the de novo RAD assembly using the MEM algorithm implemented in BWA (Li & Durbin, 2009) with default parameters. SAMtools (Li et al., 2009) was then used to sort, filter for proper pairs, remove PCR duplicates, and merge alignments for the same sample from the initial and additional rounds of sequencing, which generated a final BAM file set. SNP discovery and genotyping were performed by using the final BAM files as an input to ANGSD (Korneliussen et al., 2014) with a minimum mapping quality score (minMapQ) of 10, a minimum base quality score (minQ) of 10, the SAMtools genotype likelihood model (GL 1; Li, 2011), an SNP p-value cutoff (SNP_pval) of 1, inferring major and minor alleles (doMajorMinor 1; Skotte et al., 2012), estimating allele frequencies (doMaf 2; Kim et al., 2011), a minimum sample cutoff (minInd) of 0.5 or 0.8, and calling genotypes (doGeno) with a uniform prior (doPost 2) and posterior cutoff (postCutoff) of 0.95 or 0.99. To determine the total number of sites that passed these criteria, all sites were reported regardless of whether or not they contained a polymorphic (i.e., SNP_pval 1). However, subsequent analyses were restricted to sites with an SNP p-value less than or equal to 10–12. SNP genotypes were first discovered and called using all 192 samples. Samples with bad genotyping success, that is, from low sample/DNA/library quality, were identified from this preliminary run and removed from subsequent analysis. The remaining subset of 156 samples with good genotyping success was kept, and the analysis was rerun using this subset of 156 samples. The de novo locus discovery and assembly step has already removed highly repetitive RAD loci, and an additional filtering step was performed to further remove a small set of RAD loci that have a high frequency of polymorphic sites, which were probably results of sequencing and assembly errors.
Genetic diversity and population genomic differentiation
We used the resulting full SNP data set of 156 samples (Supplementary Table S1) to estimate several genetic diversity statistics for each population, including number of segregating sites (S), average pairwise differences (π), and Waterson’s θ (θw). To detect whether there was any population structure in genomic data, PCA was performed on SNP data of 156 samples using the bioconductor package SNPrelate (Zheng et al., 2012). A total of 255,473 SNPs with MAF ≥ 2% in 156 were used in the final PCA (Supplementary Table S3). The first three principal components were visually examined to identify any clustering pattern with habitat type or sampling location.
Populations with less than eight samples were subsequently removed, leaving 133 samples from 12 populations (Supplementary Table S1). Rare SNPs with MAF < 2% in these 133 samples were also removed using vcftools (https://vcftools.github.io/index.html), and the remaining SNPs were used for downstream analyses including Admixture, pairwise FST calculation and Bayescan outlier analysis. We also used Admixture (Alexander et al., 2009) to examine population structure, using only the first SNP of each RAD loci to reduce linkage disequilibrium. The analysis was run for K = 1–12 (Supplementary Figure S3). Phylogenetic analyses were performed using VCF2PopTree (Subramanian et al., 2019) with Jukes–Cantor distance and Neighbor-Joining method. The tree was drawn and visualized using iTOL version 6.4.1 (https://itol.embl.de/).
To quantify pairwise population differentiation, we calculated pairwise FST between populations using R package StAMPP (Pembleton et al., 2013). Statistical significance of FST was estimated using 1,000 bootstraps. Pairwise geographic distances between any two populations were calculated using geosphere package in R (Karney, 2013). The presence of isolation by distance was estimated by a simple Mantel test with 9,999 permutations using vegan package in R (Mantel, 1967; Oksanen et al., 2017). The correlation between genetic distance and whether a pair of population was from the same habitat was also tested while controlling for geographic distance using a partial Mantel test. Mantel tests and partial Mantel tests were performed using both raw FST and FST/(1 − FST), as well as both raw Euclidian geographic distances and log-transformed distances (Rousset, 1997; Slatkin, 1995).
Demographic inference using ∂a∂i
To compare alternative divergence with gene flow scenarios of B. dorothea populations and to estimate the strength of gene flow between forest and ecotone populations, we used the diffusion approximation method of software package ∂a∂i (Gutenkunst et al., 2009). We used an established analysis pipeline (Portik et al., 2017) and compared three two-population divergence models, that is, divergence with no migration, with symmetric migration, and with asymmetric migration (Supplementary Figure S4). Mbalmayo from forest and Ako from ecotone were used as representative populations because their relatively higher divergence. There are 384,612 SNPs detected in 22 individuals from Mbalmayo and 16 individuals from Ako. To avoid linage disequilibrium, first SNP of each RAD loci was used, including 41,327 unlinked SNPs. To account for missing data and to maximize the number of segregating sites used in our inference, sample sizes were projected down to 24 alleles for both Mbalmayo and Ako using Spectrum.from_data_dict from ∂a∂i. Folded 2D-SFS was used due to the lack of an outgroup to establish the ancestral state for each SNP.
Forty independent runs were performed for each model to ensure final optimization using ∂a∂i (Portik et al., 2017). For each run, there are five rounds with 50 replicates for each round. Each round was conducted by generating 50 sets of randomly perturbed parameters and optimizing each parameter set using Nelder–Mead method (optimize_log_fmin). Each optimized parameter set was used to simulate the 2D-SFS, and the log-likelihood of the 2D-SFS given the model was estimated using the multinomial approach. After each round, the parameters of the best-scoring replicate were used to generate perturbed starting parameters for the subsequent round. The generated perturbed starting values were threefold, threefold, threefold, twofold, and onefold for the first to the fifth round, respectively (Portik et al., 2017).
The best model was determined based on Akaike information criterion (AIC), and the replicate with the best likelihood of each model was used to calculate AIC scores, ΔAIC scores, and Akaike weights (ωi) (Burnham & Anderson, 2002). Four generations per year were used. Since there is no estimate of Bicyclus mutation rate, we used a mutation rate of 2.9 × 10–9 per nucleotide per generation from a different butterfly species Heliconius melpomene (Keightley et al., 2015). dadi.Godambe.FIM_uncert function was used to calculate confidence intervals for parameters from the highest likelihood of the best model (Gutenkunst et al., 2009).
Candidate SNPs under selection
To identify SNPs that are potential candidates under natural selection, we used BayeScan2.1 (Foll & Gaggiotti, 2008) to find highly differentiated SNPs between ecotone and forest habitats. This program takes a Bayesian approach to search for SNPs exhibiting extreme FST values that could be due to local adaptation. In this analysis, only the populations with a sample size no smaller than eight were used (i.e., 133 samples from 12 populations), and outliers were identified from SNPs with MAF > 2% across all good samples. We ran 5,000 iterations using prior odds of 10 and assessed the statistical significance of a locus being an outlier using a false discovery rate (FDR) of 5%.
To examine the possible role of candidate SNPs, we used RAD loci that contained these outlier SNPs to search the NCBI nucleotide collection (nr/nt) using blastn programs with default parameters. The best hit for each RAD locus with lowest e-value less than 1E-10 was identified. We also downloaded genome of a close relative species B. anynana from NCBI (GCF_900239965. 1, Bicyclus_anynana_v1.2).To view the chromosomal distribution of outlier SNPs, we map all RAD loci to B. anynana genome using blastn with default parameters. The best hit for each locus was identified using the same cutoff as above. We manually examined the annotation around the regions where RAD loci that contained these outlier SNPs mapped, and identified the closest gene. We note that the most updated B. anynana genome has 10,800 scaffolds, of which 2,233 scaffolds have RAD loci mapped. We then generated the manhattan plot using R package qqman (Turner, 2021) and q values from Bayescan analysis.
Results
Plasticity in morphological variation across habitat
In total, 1,052 second-generation individuals from four sampling localities were reared at three temperatures (22°C, 26°C, and 30°C; Supplementary Table S2). We measured 13 wing morphological characters (Zhen et al., 2023) for each individual and found that the percentage of total morphological variation explained by the first two principal components combined (PC1 and PC2) was 84.82. The diameters of eyespots loaded heavily on the first principal component (PC1: 57.93%), while the length and the width of both fore and hind wings were captured by the second principal component (PC2: 26.89%) (Supplementary Table S2). We found significant effects of sex and temperature on both PC1 and PC2, as well as the interactions between sex and temperature with habitat (forest vs. ecotone) (Table 1). Results from the linear mixed model revealed significant effects of sex and temperature on both PC1 and PC2 (p < .001), as well as the interactions between sex and habitat (p < .001), temperature and habitat (p < .001), while the three-way interaction between sex, habitat, and temperature was not significant (p = .743) (Table 1). However, no effect of habitat was detected for PC2, while the eyespots (PC1) had a strong response to habitat. The wing morphology of populations from different habitats responded differently to developmental temperature (Table 1), that is, different plasticity for populations from ecotone and rainforest habitats. We compared prominent morphological features representing PC1 (all eyespots; Supplementary Figure S2) and PC2 (wing size) across treatments (Figure 2, and ran statistical models with these features as done with the principal components (Supplementary Table S2). Eyespots grew increasingly large (almost linearly) for ecotone populations under increasing temperatures (p < .001), while forest populations eyespot size exhibited an increase from 22°C to 26°C but the reaction norm became discontinuous between 26°C and 30°C (p < .001) for both males and females (Figure 2). Results of the PCA were similar to patterns of eyespot size to wing length (Supplementary Figure S5).
Statistical summary results of the linear mixed models for the relationships between the first and second principal components (PC1 and PC2) and rearing temperatures, sex, and habitat (forest vs. ecotone) as fixed effects and sampling localities and field trip number as random effects.
| Response trait | Effects | df | Mean sq. | F value | p-value |
|---|---|---|---|---|---|
| PC1 (eyespots) | Sex | 1 | 68.63 | 20.30 | <.001 |
| Habitat | 1 | 117.03 | 34.62 | <.001 | |
| Temperature | 2 | 1,454.51 | 430.30 | <.001 | |
| Sex × Habitat | 1 | 275.40 | 81.47 | <.001 | |
| Sex × Temperature | 2 | 0.24 | 0.07 | .930 | |
| Habitat × Temperature | 2 | 219.10 | 64.81 | <.001 | |
| Sex × Habitat × Temperature | 2 | 1.00 | 0.30 | .743 | |
| Field trip (random) | 1 | 22.45 | <.001 | ||
| PC2 (wing size) | Sex | 1 | 1921.00 | 1511.88 | <.001 |
| Habitat | 1 | 2.45 | 1.92 | .165 | |
| Temperature | 2 | 89.53 | 70.46 | <.001 | |
| Sex × Habitat | 1 | 16.68 | 13.13 | .000 | |
| Sex × Temperature | 2 | 6.43 | 5.06 | .006 | |
| Habitat × Temperature | 2 | 22.84 | 18.00 | <.001 | |
| Sex × Habitat × Temperature | 2 | 0.48 | 0.38 | .685 | |
| Field trip (random) | 1 | 11.83 | <.001 |
| Response trait | Effects | df | Mean sq. | F value | p-value |
|---|---|---|---|---|---|
| PC1 (eyespots) | Sex | 1 | 68.63 | 20.30 | <.001 |
| Habitat | 1 | 117.03 | 34.62 | <.001 | |
| Temperature | 2 | 1,454.51 | 430.30 | <.001 | |
| Sex × Habitat | 1 | 275.40 | 81.47 | <.001 | |
| Sex × Temperature | 2 | 0.24 | 0.07 | .930 | |
| Habitat × Temperature | 2 | 219.10 | 64.81 | <.001 | |
| Sex × Habitat × Temperature | 2 | 1.00 | 0.30 | .743 | |
| Field trip (random) | 1 | 22.45 | <.001 | ||
| PC2 (wing size) | Sex | 1 | 1921.00 | 1511.88 | <.001 |
| Habitat | 1 | 2.45 | 1.92 | .165 | |
| Temperature | 2 | 89.53 | 70.46 | <.001 | |
| Sex × Habitat | 1 | 16.68 | 13.13 | .000 | |
| Sex × Temperature | 2 | 6.43 | 5.06 | .006 | |
| Habitat × Temperature | 2 | 22.84 | 18.00 | <.001 | |
| Sex × Habitat × Temperature | 2 | 0.48 | 0.38 | .685 | |
| Field trip (random) | 1 | 11.83 | <.001 |
p-value < 0.05: statistically significant.
Statistical summary results of the linear mixed models for the relationships between the first and second principal components (PC1 and PC2) and rearing temperatures, sex, and habitat (forest vs. ecotone) as fixed effects and sampling localities and field trip number as random effects.
| Response trait | Effects | df | Mean sq. | F value | p-value |
|---|---|---|---|---|---|
| PC1 (eyespots) | Sex | 1 | 68.63 | 20.30 | <.001 |
| Habitat | 1 | 117.03 | 34.62 | <.001 | |
| Temperature | 2 | 1,454.51 | 430.30 | <.001 | |
| Sex × Habitat | 1 | 275.40 | 81.47 | <.001 | |
| Sex × Temperature | 2 | 0.24 | 0.07 | .930 | |
| Habitat × Temperature | 2 | 219.10 | 64.81 | <.001 | |
| Sex × Habitat × Temperature | 2 | 1.00 | 0.30 | .743 | |
| Field trip (random) | 1 | 22.45 | <.001 | ||
| PC2 (wing size) | Sex | 1 | 1921.00 | 1511.88 | <.001 |
| Habitat | 1 | 2.45 | 1.92 | .165 | |
| Temperature | 2 | 89.53 | 70.46 | <.001 | |
| Sex × Habitat | 1 | 16.68 | 13.13 | .000 | |
| Sex × Temperature | 2 | 6.43 | 5.06 | .006 | |
| Habitat × Temperature | 2 | 22.84 | 18.00 | <.001 | |
| Sex × Habitat × Temperature | 2 | 0.48 | 0.38 | .685 | |
| Field trip (random) | 1 | 11.83 | <.001 |
| Response trait | Effects | df | Mean sq. | F value | p-value |
|---|---|---|---|---|---|
| PC1 (eyespots) | Sex | 1 | 68.63 | 20.30 | <.001 |
| Habitat | 1 | 117.03 | 34.62 | <.001 | |
| Temperature | 2 | 1,454.51 | 430.30 | <.001 | |
| Sex × Habitat | 1 | 275.40 | 81.47 | <.001 | |
| Sex × Temperature | 2 | 0.24 | 0.07 | .930 | |
| Habitat × Temperature | 2 | 219.10 | 64.81 | <.001 | |
| Sex × Habitat × Temperature | 2 | 1.00 | 0.30 | .743 | |
| Field trip (random) | 1 | 22.45 | <.001 | ||
| PC2 (wing size) | Sex | 1 | 1921.00 | 1511.88 | <.001 |
| Habitat | 1 | 2.45 | 1.92 | .165 | |
| Temperature | 2 | 89.53 | 70.46 | <.001 | |
| Sex × Habitat | 1 | 16.68 | 13.13 | .000 | |
| Sex × Temperature | 2 | 6.43 | 5.06 | .006 | |
| Habitat × Temperature | 2 | 22.84 | 18.00 | <.001 | |
| Sex × Habitat × Temperature | 2 | 0.48 | 0.38 | .685 | |
| Field trip (random) | 1 | 11.83 | <.001 |
p-value < 0.05: statistically significant.
Variation in the PC1 and PC2 from the mixed model (mean ± SE), representing the diameter of all eyespots and wing size, respectively. *Statistically significant (p < .001) effects of habitat type (i.e., forest vs. ecotone) on the PCs at each temperature, and ns indicates nonsignificative effects.
The model for the response of PC1 showed that the estimates of the slope of ecotone (slope = 1.65; t = 5.04; p = .025) was significantly higher than that of forest (slope = 1.35; t = 1.72; p = .084). Moreover, the Pearson’s correlation coefficient between the diameter of the eyespots and temperature was higher in ecotone (r = .736, p < .001) compared with forest (r = 0.418, p < .001) populations, suggesting that ecotone populations had stronger plasticity compared with their forest counterparts. The eyespots (PC1) varied consistently in females from both habitats at 30°C (p < .001) and in males at 22°C (p < .001), and 26°C (p < .001), while nonsignificant variation was detected in females at 22°C (p = .580) and 26°C (p = .131) and in males at 30°C (p = .656). The wing size, represented here by the second principal component (PC2), did not respond significantly to habitat (p = .165). However, within temperatures, some significant effects on the wing size of butterflies were detected; for instance, female butterflies from ecotone were larger in wing size than those from the forest at 26°C (p < .001) and males from the forest were also larger in wing size than those from ecotone at 22°C (p < .001) and 30°C (p < .001).
SNP discovery and low level of population differentiation
Bicyclus dorothea does not have a reference genome; thus, we used RAD sequencing to de novo assemble RAD loci and surveyed genome-wide diversity across the native range of B. dorothea from two different habitats, that is, rainforest and ecotone. A total of 784 million pairs of reads were obtained for 192 samples and 61,988 high-quality RAD loci was de novo assembled, ranging from 80 to 953 bp (Supplementary Figure S6). Using posterior cutoff of 0.99 and required that SNPs be genotyped in at least 80% of the final set of 156 good samples, we obtained 560,516 SNPs from 3,420,744 bp sites. The number of segregating sites per population ranges from 81,063 to 316,723. Waterson’s θ (θw) was estimated to be 0.0129 to 0.0213 per bp (mean = 0.0164), and π was estimated to be 0.0136–0.0180 per bp (mean = 0.0152), within the wide ranges of genetic diversity previous reported from other butterfly species in the subfamily Satyrinae (Mackintosh et al., 2019) (Supplementary Table S1). Across populations, forest population Mbalmayo has the highest estimates of genetic diversity.
We performed PCA and Admixture analysis and found no obvious population structure between habitats (Figure 3). Because there might be slight trend of clustering of some populations, we further focused on the four populations used in phenotypic assay, including two ecotone populations Ako and Ndikiniméki and two forest populations Mbalmayo and Somalomo (Supplementary Figure S7). We found no clustering of individuals by habitat, suggesting there was no significant genome-wide divergence between ecotone and forest habitats. By population, individuals from Ako are not only clustered but also intermixed with individuals from other populations. A range of minor allele frequency filters (from 0.5% to 20%), quality filters (posterior probability of genotypes and individual genotype missingness), and subsets of samples (removing specific population) had been tested, and the results were consistent with no clear population structure between populations from different habitats. Phylogenetic analyses of samples show a star shape relationship, with some population individual clustering but no major clustering pattern of habitats (Supplementary Figure S8).
Principal component analysis (A, B) and Admixture (C) analysis find no population structure across habitats. This analysis included all 156 good-quality samples. Samples in (A) are colored by population, while samples in (B) are colored by habitat.
For the 12 populations that had no less than eight samples, pairwise FST ranged from −0.0005 to 0.0345 (mean = 0.010; Figure 4; Supplementary Table S4), indicating no genomic differentiation across populations, consistent with PCA. In addition, there was no correlation between pairwise FST and geographic distance between populations (Mantel r = .051, mantel simulated p-value = .36; Figure 4), as well as between FST/(1 − FST) and log-transformed distances (Mantel r = .11, mantel simulated p-value = .23; Supplementary Figure S9), consistent with no population structure from PCA, suggesting that there was little isolation by distance and the levels of gene flow between populations from distant geographic locations were high to homogenize across populations. We also used partial Mantel tests to examine whether population pairs from different habitats have higher FST compared with population pairs from different habitats, controlling for geographic distance. We found no significant correlation between FST and whether a pair of populations comes from the same habitat (Mantel r = −.07, mantel simulated p-value = .65), as well as using FST /(1 − FST) and log-transformed distances (Mantel r = −.09, mantel simulated p-value = .70).
(A) Heatmap of pairwise FST between populations. Pairwise FST are low, and populations from different habitats are not more differentiated compared with populations from the same habitat. Colored lines on axis indicate the habitat type is ecotone (black line) or forest (red line). (B) No correlation between pairwise FST and geographic distance. Correlation and statistical significance were tested using Mantel and partial Mantel tests.
Demographic inference
We inferred the divergence time and strength of gene flow between populations from forest and ecotone habitats using ∂a∂i. Mbalmayo from forest and Ako from ecotone were used as representative populations because their relatively higher divergence. We compared three models of population divergence with no migration, with symmetric migration and with asymmetric migration, and found that the model with asymmetric gene flow has the best fit for the two butterfly populations (log-likelihood = −1,771.93, ωi = 1; Figure 5; Supplementary Figure S4; Supplementary Table S5). The divergence time between the forest and ecotone populations was estimated to 50,382 years ago (95% CI: 42,230–58,533). Effective population sizes after the split from common ancestor were estimated to be 171,193 (95% CI: 158,053–184,334) and 307,258 (95% CI: 280,518–333,998) for the Mbalmayo population and Ako population, respectively (Supplementary Table S6). Interestingly, gene flow was estimated to be higher from ecotone population Ako to forest population Mbalmayo than in the opposite direction. The immigrants from Ako to Mbalmayo and from Mbalmayo to Ako were approximately seven and one individuals per generation, respectively (Supplementary Table S5).
The best-fit model of divergence with asymmetric migration between forest population Mbalmayo and ecotone population Ako. Generation time of four generations per year and mutation rate of 2.9 × 10−9 per nucleotide per generation was used.
Candidate loci under selection
Although we did not detect genome-wide divergence between habitats or populations, our phenotypic assay identified significant differences in heritable wing pattern plasticity between ecotone and forest populations. This suggests there might be a small number of loci under selection for different plasticity in different habitats. To further explore potential candidate loci under selection, we identified SNPs with extreme allele frequency divergence between habitats, which could be enriched by targets of habitat adaptation. We identified 18 outlier SNPs between the meta-ecotone and meta-forest populations with an FDR of 5% using Bayescan (Supplementary Figure S10; Supplementary Table S7). These outlier SNPs were from 11 RAD loci. Nine of these RAD loci have hits in NCBI nr database using blastn with E-value less than 1 × 10−10. Two of these RAD loci have best hits on annotated genes, that is, predicted B. anynana odorant receptor Or1-like (LOC112057728) mRNA and predicted B. anynana B-cell receptor CD22 (LOC112052937) mRNA that may play a role in immune system. By mapping RAD loci to B. anynana genomes using blastn, we mapped 10 of the 11 loci to 10 different B. anynana scaffolds, including these two. We further identified four RAD loci mapped to the intronic regions of hemicentin-1 (a large extracellular member of the immunoglobulin superfamily), aminopeptidase N that involves in oxidative stress and redox pathway, as well as two RAD loci mapped to annotated but uncharacterized genes (Supplementary Table S8).
Discussion
Under a common garden environment, B. dorothea individuals from ecotone and forest habitats show contrasting patterns of phenotypic plasticity. Specifically, ecotone populations show greater plasticity in wing eyespot size than forest populations as eyespot size increased at high temperatures in ecotone populations but did not change in forest populations. However, no apparent genome-wide genetic differentiation was found between B. dorothea populations from forest and ecotone habitats. Demographic inferences support a model of divergence with asymmetric gene flow between forest and ecotone populations. Together the results suggest that habitat differences are important in structuring phenotypic responses to climate but that the lack of clear genomic differentiation between populations indicates that the phenotypic patterns are likely a consequence of a few loci or potentially epigenetic in nature.
The inferred demographic parameters suggest large effective population sizes and high level of genes flow between forest and ecotone populations, with more migrants from ecotone to forest than in the opposite direction. While we have no estimates of dispersal for B. dorothea in Cameroon, our observations suggest that individuals are capable of moving effectively through open habitats and forest edges, which could facilitate the gene flow between populations. In such a case, selection could promote rapid divergence in loci underlying divergent phenotypes, but a relatively high level of ongoing gene flow, frequent recombination, and low linkage disequilibrium could lead to relatively undifferentiated pattern in most parts of the genome, which has been reported in a number of studies, even for currently recognized species (Amaral et al., 2018; Branch et al., 2017; Campagna et al., 2012; Cooper & Uy, 2017; Mason & Taylor, 2015; Palacios et al., 2019; Semenov et al., 2018).
To identify potential candidate genes under divergent selection between forest and ecotone habitats, we detected 18 outlier SNPs located in 11 RAD loci using Bayescan. None of these outliers is located within or close to genes with known function in wing developmental plasticity. Two outlier RAD loci map to genes with functions in immune response. Insect pigmentation is also known to associate with immunity (Wittkopp & Beldade, 2009). Alternatively, immunity-related genes are among the classes of genes that are known to evolve most rapidly when organism adapt to different environmental conditions (Early & Clark, 2017). In addition, recent findings suggest that immune activity is associated with life span in a seasonally polyphonic butterfly (Freitak et al., 2019). Another outlier RAD loci maps to an odorant receptor associated with environmental cues sensing, aminopeptidase N in oxidative stress and redox pathway, and P450. Due to limited knowledge of Bicyclus genome annotation and the sparsity of RAD loci markers across the genome, most of the outlier RAD loci could not be linked to known functionality. Genomic differences that could lead to variation in the degree of wing eyespot plasticity may directly related to the trait or indirectly in the steps including sensing environmental cues and modulating these signals. Further investigation is needed to examine whether these candidate SNPs play roles in wing eye spot plasticity. Additionally, it is possible that there is only a small number of loci influencing wing plasticity between forest and ecotone populations, but they are not included in or linked to our RAD loci (0.27% genome coverage if using B. anynana genome size as an approximation). Further investigation with whole-genome resequencing and a higher quality reference genome would provide more power to identify candidate variants under selection between habitats. Lastly, it is also possible that the observed difference in phenotypic plasticity between habitats could be controlled by heritable epigenetic factors (Flatscher et al., 2012).
Two previous studies on the little greenbul (Andropadus virens), an African songbird, and an African rodent Praomys misonnei and a skink (Trachylepis affinis) found that differential natural selection along the ecotone–rainforest environmental gradient drives clear population differentiation at the genomic level (Freedman et al., 2023; Morgan et al., 2020; Zhen et al., 2017). While sampling populations along the same forest–ecotone gradient and using similar RADseq data, this pattern of genomic differentiation was not found in B. dorothea butterfly. This pattern of no differentiation is also found in another species of African sunbird (Arriens et al., unpublished). The distinct patterns of observed genomic differentiation across this environmental gradient among species suggest that species respond to environmental gradients differently. This implies that ecological speciation may happen with different probabilities for different species confronted by the same environmental challenges. However, the demographic histories of these species are not well understood, and we cannot rule out other confounding factors that could influence the neutral divergence patterns, including colonization timing or population size.
Bicyclus eyespots have been intensively studied across multiple species and the regulatory pathways (e.g., 20-hydroxyecdysone 20E) that determine plasticity in eyespot size are being increasingly characterized and understood (Beldade & Peralta, 2017; Monteiro, 2017). The reaction norm with a positive slope in wing eyespots (larger eyespots at higher temperatures) found here in B. dorothea is consistent with other studies in species of Bicyclus (Bhardwaj et al., 2020; Oostra et al., 2014). Furthermore, Roskam and Brakefield (1996) reared six species of Bicyclus, among which five originated from Malawi and one from Cameroon. They found that species originating from ecotone habitats are more sensitive to environmental fluctuations with respect to eyespot features. Our results further demonstrate that such reaction norms also vary between populations within the same species. In contrast to the eyespot results, the wing size reaction norms did not vary linearly across habitat (though differences were apparent), consistent with the supposition that this life-history trait is less likely to be under variable selection pressure relative to other morphological features like eyespots (De Jong et al., 2010; van Bergen et al., 2017). Furthermore, local adaptation in thermal tolerance has previously been demonstrated for B. dorothea in ecotone versus forest (Dongmo et al., 2021) and in other tropical butterflies across environmental gradients (e.g., Heliconius; Montejo-Kovacevich et al., 2020).
Plasticity in eyespot in wings of B. dorothea is likely to be adaptive, similar to B. anynana (Lyytinen et al., 2004; Prudic et al., 2015). Greater plasticity in ecotone is likely to be adaptive in ecotone habitat because of the greater environmental fluctuation across seasons (Supplementary Figure S1). Higher heterogeneity in vegetative cover in the ecotone will increase exposure to thermal extremes and could necessitate mechanisms for managing this variation. While we largely focused on wing traits here, Bicyclus species exhibit a wide range of polyphenic traits in response to temperature including morphological characteristics (Mateus et al., 2014; Monteiro et al., 2015), life-history traits such as body size and reproductive diapause (Halali et al., 2021; Pijpe et al., 2007), physiological traits such as relative fat content and resting metabolic rate (Oostra et al., 2010), and behavioral traits such as mate choice and predator avoidance behavior (Prudic et al., 2011; van Bergen & Beldade, 2019). Key interactions or the summation of these multiple traits across environmental gradients will ultimately structure population and demographic changes.
In summary, we find different levels of phenotypic plasticity between B. dorothea originally from different habitats across rainforest–ecotone gradients in common garden experiments. These habitat-specific thermal responses in B. dorothea are likely important for adaptation to different environments. However, we find no genome-wide differentiation among populations between habitats. Further investigation is necessary to fully understand the genomic or epigenomic underpinnings of such habitat-specific responses.
Data availability
RADseq data: NCBI SRA PRJNA890878. Morphological and genotype data: https://doi.org/10.5061/dryad.z08kprrj3. Scripts are available at https://github.com/ZhenLab-Westlake/bicyclus-RAD
Author contributions
T.B.S., T.C.B., R.H., and Y.Z. conceived of and supervised the project. K.R. supervised the laboratory work. Y.Z. analyzed genomic data. M.D. collected specimens and performed phenotypic assays. Q.F. performed demographic inference. Y.Z. and M.D. wrote the manuscript with input from all authors.
Conflict of interest: The authors declare no conflict of interest.
Acknowledgments
We thank Mike Miller for his assistance with the RAD library preparation and bioinformatic support and Kevin Njabo for helping in field sample collection. We acknowledge support from a QCB Collaboratory Postdoctoral Fellowship to Y.Z. and the QCB Collaboratory Community directed by Matteo Pellegrini. This research was supported by National Science Foundation grants IIA PIRE #1243524 to T.B.S., a General Research Fund award #17118317 from the Research Grants Council Hong Kong to T.C.B. and M.A.K.D., National Natural Science Foundation of China #31900315 and Zhejiang Provincial Natural Science Foundation of China grants #LR21C030001 to Y.Z.
References
Author notes
Y.Z. and M.A.K.D. contributed equally.
