Contrasting effects of habitat discontinuity on three closely related fungivorous beetle species with diverging host‐use patterns and dispersal ability

Abstract Understanding how landscape structure influences biodiversity patterns and ecological processes are essential in ecological research and conservation practices. Forest discontinuity is a primary driver affecting the population persistence and genetic structure of forest‐dwelling species. However, the actual impacts on populations are highly species‐specific. In this study, we tested whether dispersal capability and host specialization are associated with susceptibility to forest discontinuity using three closely related, sympatric fungivorous ciid beetle species (two host specialists, Octotemnus assimilis and O. crassus; one host generalist, O. kawanabei). Landscape genetic analyses and the estimation of effective migration surfaces (EEMS) method consistently demonstrated contrasting differences in the relationships between genetic structure and configuration of forest land cover. Octotemnus assimilis, one of the specialists with a presumably higher dispersal capability due to lower wing loading, lacked a definite spatial genetic structure in our study landscape. The remaining two species showed clear spatial genetic structure, but the results of landscape genetic analyses differed between the two species: while landscape resistance appeared to describe the spatial genetic structure of the specialist O. crassus, genetic differentiation of the generalist O. kawanabei was explained by geographic distance alone. This finding is consistent with the prediction that nonforest areas act more strongly as barriers between specialist populations. Our results suggest that differences in host range can influence the species‐specific resistance to habitat discontinuity among closely related species inhabiting the same landscape.

and reduce local population size. However, the actual impacts on populations can be highly species-specific. While theoretical and empirical studies have reported negative effects of forest discontinuity on population persistence (Fahrig, 2002;Gibson et al., 2013), some species are less sensitive to forest discontinuity (Didham, Hammond, Lawton, Eggleton, & Stork, 1998;Lampila, Monkkonen, & Desrochers, 2005). This variation in sensitivity may be related to dispersal capacity and several ecological characteristics of the species (Henle, Davies, Kleyer, Margules, & Settele, 2004). Among ecological predictors of species sensitivity, specialization in habitat use and diet resources have been hypothesized to be key determinants (Keinath et al., 2017;Khimoun et al., 2016). Specialist species are less likely to disperse through areas where habitat patches are sparsely distributed, because, compared to generalists, they fulfill their resource requirements in smaller subsets of habitat patches and are more susceptible to local fluctuations of resources. Thus, nonforest areas act more strongly as barriers between specialist populations. In addition, specialist species tend to be patchily distributed, which increases differentiation among populations (Janz, Nylin, & Wahlberg, 2006) relative to generalist species. This pattern is expected to be more conspicuous in landscapes with discontinuous habitat. Correlations between ecological specialization and numerical responses of populations and communities to habitat fragmentation have been demonstrated in several taxa (e.g., birds: Devictor, Julliard, & Jiguet, 2008;butterflies: Steffan-Dewenter & Tscharntke, 2000). However, such changes in population and community structure can be driven by several factors (e.g., environmental change accompanied by fragmentation, correlation between specialization, and movement behavior). Therefore, it is important to quantify the dispersal patterns of organisms in discontinuous habitats to improve our understanding of the effects of habitat discontinuity on population structure.
While the direct observation and quantification of movement behavior are costly and nearly impossible to conduct, the spatial genetic structure of a population enables us to infer the extent and routes of effective dispersal. Reduced dispersal between habitat patches will decrease gene flow among populations and thus increase genetic differentiation. Recent developments of landscape genetic methods allow researchers to test the effects of environmental change and habitat connectivity on gene flow between populations (Balkenhol, Cushman, Storfer, & Waits, 2015). In particular, a pairwise F ST approach has been employed to test the effects of landscape quality on gene flow rates under different scenarios based on a null hypothesis of the absence of geographic structure (Balkenhol, Waits, & Dezzani, 2009). In this approach, an isolation-by-distance (IBD) scenario assumes that genetic differences increase with geographic distance due to limited dispersal across space, whereas an isolationby-resistance (IBR) scenario predicts a relationship between genetic differentiation and resistance distance, indicating the differential effects of landscape features on dispersal (McRae, 2006). The IBR concept aims to characterize how genetic differentiation is shaped in heterogeneous landscapes, and "resistance" represents the cost to an organism to cross a particular environment, whereby a low resistance denotes ease of movement and a high resistance denotes restricted movement (Zeller, McGarigal, & Whiteley, 2012). When applying these scenarios to population responses to forest discontinuity, the IBD model indicates limited dispersal but the absence of impacts of habitat isolation, and the IBR model indicates significant effects of the loss of habitat continuity on population structure.
Recently, a number of empirical landscape genetics studies have been conducted for a variety of taxa (Balbi et al., 2018;Beninde et al., 2016;Cleary, Waits, & Finegan, 2017;Crawford, Peterman, Kuhns, & Eggert, 2016;Frantz et al., 2012;Goldberg & Waits, 2010;Reid, Mladenoff, & Peery, 2017). However, most studies have focused on a single species or multiple species that largely differ in several characteristics (but see Engler, Balkenhol, Filz, Habel, & Rodder, 2014;Kelley, Farrell, & Mitton, 2000). Comparisons of closely related species that differ in their extent of ecological specialization on the same landscape would facilitate the examination of the effects of specialization on sensitivity to forest discontinuity.
Here, we perform a comparative population genetic study among closely related, sympatric ciid beetle (Coleoptera: Ciidae) species to test whether host specialization is associated with susceptibility to forest discontinuity. Ciid beetles are fungivorous and inhabit and feed on the basidiomes (fruiting-bodies) of bracket fungi (Basidiomycetes). Most species of Ciidae feed on a relatively restricted number of fungal taxa (Fossli & Andersen, 1998;Lawrence, 1973;Økland, 1995;Orledge & Reynolds, 2005;Paviour-Smith, 1960). Because their hosts, wood-rotting bracket fungi, depend on the existence of dead woods, forests are considered potentially suitable and resource-rich habitats for ciid beetles. The basidiomes of fungi are a relatively ephemeral and highly fluctuating resource, and they can occasionally disappear from small, isolated habitats.
Fungus-feeding species that can use multiple fungal species are expected to have a greater likelihood of fulfilling their resource requirements in such patches. Ciid beetles provide an ideal system for the study of spatial ecology in forest ecosystems, because they are abundant in number and depend on the basidiomes of bracket fungi at all stages of their life cycle. Several colonization experiments of insects on deadwood, including Ciidae, have suggested that the ability of insects to colonize isolated patches is highly species-specific (Jonsell, Nordlander, & Jonsson, 1999;Komonen, 2008 In this study, we compare population genetic structure among the above three Octotemnus species inhabiting the same landscape to test the prediction that ecologically specialized species (specialists; species with narrower host range) are more sensitive to forest discontinuity than generalist species (species with a broader host range). We hypothesized that nonforest areas will act more strongly as a barrier for specialist species than for generalist species and that compared to the simple IBD model, IBR scenarios will better explain the population structure of the specialist species when they do not differ in their dispersal abilities. We used microsatellite data and performed resistance surface optimization and applied the estimation of effective migration surfaces (EEMS) model to landscape population genetic structure of individual species. We evaluated the dispersal ability of focal species using morphological data. We found that different species showed varying levels of response to forest discontinuity, which can be explained by differences in dispersal ability and host specialization. has remained almost unchanged for 100 years. Therefore, the discontinuity among the study sites was considered to be longstanding. Insects were collected from the basidiomes of Trametes and Lenzites (and unidentified) species ( Figure 1). All beetle specimens were preserved in 99% ethanol until DNA extraction. See Supporting Information Tables S1 and S2 for detailed information of specimens used in this study. Overlapping paired reads were merged using PANDAseq (Masella, Bartram, Truszkowski, Brown, & Neufeld, 2012). The selection of merged reads containing microsatellites and the design of primers were conducted using QDD 3.1.2 (Meglecz et al., 2014). The universal tail sequence for fluorescent labeling of PCR fragments (Blacket, Robin, Good, Lee, & Miller, 2012) was added to forward primers. Loci were screened for PCR amplification success and polymorphism.
F I G U R E 1 Maps of land cover and sampling sites in Kyoto (Japan). Left, land cover types in the study area: deciduous broad-leaved forest (yellow), evergreen broad-leaved forest (green), conifer plantation (blue), arable land (brown), city (white), and others (gray). Right, sampling points of host fungi. Symbols represent host-fungal species from which beetles were collected. Areas of forest (gray) and nonforest (white; mainly city) are also shown TA B L E 1 Characteristics of the newly developed microsatellite markers Note. NA indicates that the locus was not amplified.

| Genotyping and summary statistics
We amplified 9-10 microsatellite loci for each species with two multiplexes of five to six loci. Multiplex PCR was performed in 4.5 µl reaction volumes containing 1X Type-it Multiplex PCR Master Mix In the following analysis, individuals collected from sites close to one another (typically <1 km), as well as those collected from the same fungal bodies, were treated as belonging to the same population. On average, four to six individuals per population were genotyped for each species. Populations with fewer than four individuals were excluded from the calculation of G″ ST (Meirmans & Hedrick, 2011). We checked null alleles using the Micro-checker software

| Landscape genetics analyses
We conducted landscape resistance analyses to test our hypothesis that the species differ in their responses to landscape type (forest/ nonforest). We obtained our land cover data from the database of the Biodiversity Center of Japan (http://www.biodic.go.jp/trialSystem/top_en.html). The land cover data of our study site are based on vegetation surveys conducted since 1999. The original vector format data were rasterized at 100-m resolution (the smallest census unit of vegetation data) to perform subsequent landscape analyses.
In addition, the original vegetation types were reclassified into two (forest and nonforest) or six (deciduous broad-leaved forest, evergreen broad-leaved forest, conifer plantation, arable land, city, and others) categories (Figure 1). Land cover types occupying <5% of the study area were reclassified as "other." We followed the framework of optimization and selection of resistance surfaces using the "ResistanceGA" package (Peterman, 2018) in R. This method uses a genetic algorithm (GA; Scrucca, 2013) to optimize resistance surfaces to the pairwise genetic distances and conducts model selection to determine the best-supported resistance surface. A linear mixed-effects model with MLPE is fit to the data in model selection. We used pairwise G″ ST values between sampling sites as input data and assessed model fits using the Akaike information criterion (AIC). We assessed the relative support of three competing models: the IBD model, which proposes that gene flow is a function of the Euclidian distance among populations; the IBR model, which proposes that gene flow is a function of the resistance distance; and a null model (absence of geographic structure). Bootstrap analyses were conducted using the resist.boot function to evaluate the relative support of competing distance models. In each bootstrap replication, pairwise response and distance matrices are subsampled and fitted to the MLPE model to the data to obtain statistics. The percentage of instances of the IBD or IBR model being the best-fit model was used as the support level.

| Estimated effective migration surfaces
We visualized how the IBD relationship varies across geographic space using Estimated Effective Migration Surfaces software (EEMS; Petkova, Novembre, & Stephens, 2016). This method estimates effective migration rates based on genetic distances and then creates a visual representation of effective migration rates by interpolation. EEMS estimates the effective migration across space without the need to observe environmental variables and thus provides an exploratory tool for spatial population structure. This exploratory approach is complementary to the hypothesis-driven resistance surface approach described above. We set the number of demes to 200 and ran three independent analyses with 1,000,000 burn-in Markov chain Monte Carlo steps and 2,000,000 iterations. The results of three runs were combined using the rEEMSplots R package (Petkova et al., 2016).

| Estimation of potential flight capability
It is believed that the study beetle species usually disperse by flight, because they have well-developed hind wings and are frequently collected by flight-intercept traps. We compared flight morphology of three species to evaluate relative dispersal ability. Specimens were collected from host fungi within the study site of landscape genetic analyses below (see Supporting Information Table S1 in detail). Beetles were killed and preserved in 100% ethanol for at least 48 hr and dried at room temperature for 24 hr. Body mass was measured using a digital balance (Sartorius BP 210D, Göttingen, Germany) to the nearest 0.01 mg. Subsequently measured beetles were digested in Nuclei Lysis Solution (Promega, Madison, WI, USA) with proteinase K (×mg/ml) at 55°C overnight, to easily dissect the hind wings. The left wing was removed and mounted in drops of mounting medium (Euparal). The length and width of the pronotum and elytra and the length, width, and area of the hind wings were measured using a VW-9000 microscope with a VW-600C camera and VH-Z 100R zoom lens (Keyence, Osaka, Japan). In total, 48 individuals (eight males and eight females of each species) were measured. Wing loading (body mass divided by wing area) and wing aspect ratio (wing length divided by wing width) for each individual were calculated. Body mass was highly variable among individuals ( Figure 2), likely because of differences in sexual development and gut contents; therefore, body length (sum of pronotal and elytral length) was used as a proxy of body mass to avoid such confounding influences. Pairwise differences between sex and species were examined using t tests, and Bonferroni adjustments were applied to p-values.

| Genetic diversity and population differentiation
In total, 21 microsatellite loci were used for the three Octotemnus species in this study (Table 1). In all, 9, 10, and 10 loci were poly-

| Landscape resistance analyses
The model selection results differed among the three Octotemnus species, as did the optimized circuit resistance distance in ResistanceGA (Table 2). In O. crassus, the IBR model with six land cover categories was the best-fit model, followed by the IBR model with two land cover categories. In the 6-land cover IBR model, deciduous broad-leaved forests and conifer plantations had lower resistance values (1 and 58, respectively), and evergreen broad-leaved forests, arable land, and city had higher resistance values (1,212, 2,415, and 1,159, respectively). In the 2-land cover IBR model, forests had a lower resistance value than nonforest land cover (1.0 vs. 13.6). The two IBR models were selected with a higher bootstrap percentage (65.4, 33.8%) than the distance model (1.3%), indicating effects of forest cover on population genetic structure. For O. kawanabei, the IBD model was supported, suggesting relatively limited dispersal; however, the estimated resistance values of forest and nonforest areas did not significantly differ (1.2 vs. 1.0 in the 2-land cover model). We found no significant effects of IBD or IBR on genetic variation for O. assimilis.

| Flight morphology
Octotemnus assimilis was smaller and lighter and had significantly lower wing loadings compared to O. crassus and O. kawanabei ( Figure 6). On average, the wing loadings of both O. crassus and O. kawanabei were 1.41 times higher than that of O. assimilis. No significant differences were detected between sexes in any of the three species. Wing aspect ratios did not differ between species or sexes.

| D ISCUSS I ON
The closely related fungus-feeding Octotemnus beetles exhibited differences in response to forest discontinuity. Analyses using both resistance surface and EEMS methods yielded similar results for each species. The interspecific differences could be associated with species' differences in dispersal capability and ecological specialization. Among the three beetle species examined, O. assimilis presumably has higher dispersal capability than the other two species because its wing loading is much lower (Figure 6).
The lack of spatial genetic structure on our study landscape for O. assimilis (Figures 3 and 4) (Dudley & Srygley, 1994). For example, wing loading and flight distance of monarch butterflies are negatively correlated in flight-mill experiments (Bradley & Altizer, 2005). However, it is uncertain whether wing loading is actually a reliable predictor of fight and dispersal capability in ciid beetles. Further laboratory and field studies are needed on the flight behavior of ciid species to clarify this matter.
Host use plays an essential role in the evolution and diversification of various organisms (Forbes et al., 2017;Hoberg & Klassen, 2002;Poulin & Morand, 2000). Specialization in host use can likely facilitate population differentiation and promote species diversification, because suitable habitats are generally more patchily distributed for specialists than for generalists and hence gene flow is more limited in specialist compared to generalist populations (Janz et al., 2006). This long-standing hypothesis has been tested in a variety of taxa, and the results of many of these studies have agreed with the prediction (Brouat, Chevallier, Meusnier, Noblecourt, & Rasplus, 2004;Engler et al., 2014;Kelley et al., 2000;Zayed et al., 2005; but see e.g., Peterson & Denno, 1998). Our results from fungus-feeding organisms are consistent with the prediction that specialization promotes population genetic subdivision.
A number of studies have explored the relationship between ecological specialization and sensitivity to habitat discontinuity; however, few studies have explicitly incorporated spatial genetic structure into the analyses. Such studies not only provide guidelines for conservation practices but also offer insight into the mechanisms of species diversification and biogeography. Recent developments of high-throughput sequencers enable us to analyze many species and individuals in a single study at a low cost. By conducting additional comparative studies of multiple sets of closely related species, more generalized patterns can be explored.

ACK N OWLED G M ENTS
This study was supported by a Grant-in-Aid for Japan Society for the Promotion of Science (JSPS) Fellows (Grant no. 17J07928 to T.K.) and a JSPS KAKENHI (Grant no. 15H02637 to T.S.).
F I G U R E 6 Measurements of (a) wing loading (cube of body length divided by wing area) and (b) wing aspect ratio (wing length divided by wing width) of the three Octotemnus species. Means with different letters are significantly different from each other (t test and Bonferroni adjustments). Error bars represent standard deviation