A Comparative Study of Genetic Responses to Short- and Long-Term Habitat Fragmentation in a Distylous Herb Hedyotis chyrsotricha (Rubiaceae)

The genetic effects of habitat fragmentation are complex and are influenced by both species traits and landscape features. For plants with strong seed or pollen dispersal capabilities, the question of whether the genetic erosion of an isolated population becomes stronger or is counterbalanced by sufficient gene flow across landscapes as the timescales of fragmentation increase has been less studied. In this study, we compared the population structure and genetic diversity of a distylous herb, Hedyotis chyrsotricha (Rubiaceae), in two contrasting island systems of southeast China. Based on RAD-Seq data, our results showed that populations from the artificially created Thousand-Island Lake (TIL) harbored significantly higher levels of genetic diversity than those from the Holocene-dated Zhoushan Archipelago (ZA) (π = 0.247 vs. 0.208, HO = 0.307 vs. 0.256, HE = 0.228 vs. 0.190), while genetic differences between island and mainland populations were significant in neither the TIL region nor the ZA region. A certain level of population substructure was found in TIL populations, and the level of gene flow among TIL populations was also lower than in ZA populations (m = 0.019 vs. 0.027). Overall, our comparative study revealed that genetic erosion has not become much stronger for the island populations of either the TIL or ZA regions. Our results emphasized that the matrix of water in the island system may facilitate the seed (fruit) dispersal of H. chrysotricha, thus maintaining population connectivity and providing ongoing resilience to the effects of habitat fragmentation over thousands of years.


Introduction
Over the past 200 years, human activities have altered more than 80% of the terrestrial world, leading to increasing habitat reduction and fragmentation of almost all ecosystems [1][2][3]. Theoretically, habitat fragmentation has severe negative impacts on community structures, population distributions and abundances, patterns of genetic variation, gene flow, etc. [4,5]. These changes, especially genetic alterations, can reduce the individual fitness and evolutionary potential of populations, further impeding the long-term persistence of species and increasing the risk of local extinction [6,7].
Detecting the genetic effects of habitat fragmentation on natural populations is challenging in practice, as there are many confounding factors involved [2]. For instance, Primula veris differed in fine-scale spatial genetic structure patterns at a small scale, which may be explained by partial self-compatibility of the pin morph combined with restricted seed dispersal and pollinator behavior. In our previous genetic study [13], we surveyed populations of H. chrysotricha from the recently fragmented TIL region only, using nuclear SSR markers. Although SSRs are highly polymorphic markers, a limited number of markers might result in low power for addressing landscape genetic questions at large spatial and temporal scales [30][31][32]. With the decreasing cost of next-generation sequencing, restriction site-associated DNA sequencing (RAD-Seq), which enables the fast identification of thousands of single nucleotide polymorphisms (SNPs) in non-model organisms without any prior information, has become a cost-effective approach in phylogeographic and population genetic studies [33]. In this study, we employed genome-wide SNPs derived from RAD-Seq data to examine and compare short-term vs. long-term effects of habitat fragmentation on H. chrysotricha in the TIL vs. ZA regions. Specifically, we aimed to: (1) investigate the effects of short-term vs. long-term habitat fragmentation on the species' population genetic variation (diversity, inbreeding) patterns; (2) examine the genetic differentiation and genetic structure of H. chrysotricha populations at small and large spatial scales; and (3) evaluate patterns of gene flow and demographic changes of H. chrysotricha populations in the two island systems. This comparative study will increase our understanding of the genetic responses of short-lived plants with strong seed dispersal capabilities to fragmentation and provide an insight into the management and conservation of subtropical forest dwellers in fragmented habitats at both small and large spatial-temporal scales.

RAD-Seq Data and SNP Filtering
Approximately 389.4 billion (389,411,194,200) raw reads were produced from 256 individuals of H. chrysotricha sampled across the TIL and ZA regions. After quality control, filtering, and trimming, an average of 1421.5 million and 1638.9 million high-quality reads were retained for the TIL and ZA samples, respectively (Table S1). A catalog containing 9,131,292 loci was constructed, and 9,067,359 loci were genotyped by GSTACKS using the data sets of all H. chrysotricha samples. The mean, minimum, and maximum values for effective per-sample coverage were 26.9×, 16.8×, and 53.0×, respectively. After filtering loci of low quality (minor allele frequency < 0.01; missing rate > 0.5), 185 and 188 unlinked SNPs were eventually identified in TIL and ZA populations, respectively, and used for all subsequent analyses (Tables S2 and S3).

Population Variation
Based on genome-wide SNP markers, we detected significantly higher levels of genetic diversity in H. chrysotricha from the TIL region than from the ZA region (π = 0.247 vs. 0.208,  (Tables 1 and 2).  Pairwise F ST values between the TIL populations ranged from 0.001 to 0.393, with about 13% being significant (p < 0.05) after sequential Bonferroni correction (Table 3). Among the significant pairwise comparisons, the great majority (11 out of 16) involved population EM04. Pairwise F ST values between the ZA populations ranged from 0.001 to 0.436, with about 24% being significant (p < 0.05) after sequential Bonferroni correction (Table 4). Here, a large proportion of significantly high F ST values (6 out of 13) involved population ZI05. The average F ST value across populations was slightly higher in the TIL (F ST = 0.103) region than in the ZA (F ST = 0.092) region, yet this difference was only marginally significant (p = 0.046). No IBD pattern was found in either region (TIL: r = 0.181, p = 0.174; ZA: r = 0.072, p = 0.242) ( Figure S1).

Population Genetic Structure
Based on a STRUCTURE analysis, the most likely number of genetic groups was K = 3 for the TIL populations and K = 2 for the ZA populations. In the TIL region, Cluster I ('dark blue') was present at high frequency (80% of all local samples) in seven of eight island populations, whereas the great majority of individuals (89%) from the western mainland populations and the remaining island population (IP08) were assigned to Cluster II ('red') ( Figure 1 and Figure S2). In the eastern mainland populations, half of the individuals from EM04 were assigned to Cluster III ('green'); the remaining populations contained variable numbers of individuals belonging to both Cluster I and Cluster II ( Figure 1). In the ZA region, nearly all populations (excepting ZI04) consisted of variable numbers of individuals with a probability of ≥ 0.5 of belonging to both Cluster I ('blue') and Cluster II ('orange'), when K = 2 ( Figure 2). At K = 3-5, we also observed that the population ZM02 was assigned to a separate cluster ( Figure S3), suggesting potential genetic differentiation between this mainland population and the others.

Contemporary Gene Flow
In H. chrysotricha populations from the TIL region, the proportion of individuals that originated from within the same site ranged from 68% (IP04) to 81% (WM02), with an average value of 72% (see BAYESASS, Table S4). In contrast, rates of gene flow (or 'migration', m) between populations were low to moderate, and in the majority of populations no more than 5% of the individuals were exchanged with other populations (Table S4). Similarly, in the ZA region, populations largely consisted of individuals originating from the same site (73% on average) (see BAYESASS, Table S5). Nevertheless, the average rate of interpopulation gene flow in the ZA region (m ZA = 0.027) was significantly higher than in the TIL region (m TIL = 0.019; p = 0.003). The DIVMIGRATE method also suggested that the relative migration rate was significantly higher in the ZS region than in the TIL region (rm ZA = 0.349, rm TIL = 0.252, p < 0.01) (Tables S6 and S7). In addition, we observed significant asymmetric gene flow patterns in 10 out of the 16 populations in the TIL region, and 4 out of the 11 populations in the ZA region (Figures 3 and 4).

Demographic (Bottleneck) Analyses
Irrespective of the mutation model assumed (IAM, infinite allele model; SMM, stepwise mutation model; or TPM, two-phase mutation model), we found a significant excess of heterozygosity (p < 0.05) for one island (IP05) and one mainland population (EM04) in the TIL region, and one island population (ZI01) in the ZA region (Table 5), indicating that these populations may have experienced a recent bottleneck. However, a significant excess of heterozygosity was only found under the IAM model for one TIL population (IP08) and five ZA populations (ZI02, ZI07, ZM02, ZM03, ZM04; Table 5).

Discussion
Habitat fragmentation changes continuous habitats into small and isolated patches, and populations inhabiting these habitats are often considered to have low genetic diversity [8,34]. In our previous study [12], we found that island populations of H. chrysotricha in the TIL region had significantly lower mean genetic diversity than those from the mainland, based on nuclear SSR markers, suggesting that anthropogenic habitat fragmentation can lead to significant loss of genetic diversity over a few decades. In the present study, we further tested whether genetic erosion increases with an increase in spatio-temporal scales of fragmentation. Based on RAD-Seq-derived SNP markers, we observed that the average estimates of genetic diversity for ZA populations were significantly lower than for TIL populations (π = 0.247 vs. 0.208, p < 0.01; H O = 0.307 vs. 0.256, p < 0.01; H E = 0.228 vs. 0.190, p < 0.01; Tables 1 and 2), and we found no major differences in genetic diversity between island and mainland populations in either the TIL or the ZA region. Although it remains unknown how much of the initial genetic variation has been lost in each region since the formation of islands, these results may indicate that H. chrysotricha populations have retained a considerable proportion of their initial genetic diversity over either~60 years or thousands of years of habitat fragmentation.
Generally, the genetic effects of habitat fragmentation on plants vary, depending on their life-history traits [2]. In this study, the lack of a significant effect of habitat fragmentation on the genetic diversity of H. chrysotricha may be explained by demographic factors and/or strong gene dispersal capacity. As a distylous plant species, H. chrysotricha has two contrasting flower morphs that reciprocally differ in the spatial separation of stigmas and anthers, which in turn prevents selfing and intramorph mating. Morph bias is expected to reinforce genetic drift effects and inbreeding by reducing mating opportunities [35,36]. We also found that no biased morph ratio existed in either island or mainland populations, and overall levels of inbreeding in both island regions were negative and close to zero, suggesting a sufficient within-population heterozygosity (Tables 1 and 2). In addition, most of the populations showed no genetic signatures of recent bottlenecks (Table 5). Therefore, it is likely that population sizes of H. chrysotricha have remained sufficiently large to prevent loss of genetic diversity via inbreeding and genetic drift after fragmentation [37,38].
The gene dispersal ability via pollen or seed is another key characteristic that determines the potential of a plant species to counteract the negative effects of habitat fragmentation [8]. Hedyotis chrysotricha is an insect-pollinated plant; its flower visitors include solitary bees (Halictidae and Andrenidae), Diptera (Syrphidae and Bombyliidae), honey bees, small lepidopterans, and/or thrips [13,39]. Due to the limited foraging distances of these pollinators, they seem unlikely to be able to maintain significant levels of gene flow and population connectivity across fragmented habitats. In contrast, the indehiscent fruit of H. chrysotricha contains several very small seeds (1.5-2 × 2-2.5 mm), whose dispersal should be easily facilitated by wind [21]. In our previous study, we detected a moderate level of gene flow among 18 H. chrysotricha populations in the TIL region [12]. Based on high-quality SNP markers, our genetic analyses here have further verified the undisturbed population connectivity of H. chrysotricha populations in the TIL region, with an average of c. 28% individuals in each population being exchanged with other sites (Table S2). Given the relatively large spatial and temporal scales of habitat fragmentation in the ZA region ( Figure 5), one might expect high population differentiation and low levels of gene flow between island populations. However, both BAYESASS and DIVMIGRATE analyses suggested that the average level of gene flow among populations in the ZA region was higher than in the TIL region, although the results from these two methods may not be compared directly. Similarly, populations in the ZA region exhibited marginally lower levels of genetic differentiation than those in TIL region (F ST = 0.092 vs. 0.103). A previous study showed that water greatly facilitated the seed dispersal of a typical wind disperser across a fragmented landscape [40]. By referencing the population genetic study of the maritime Hedyotis strigulosa var. parviflora (Hook. et Arn.) Yamazaki [14], we speculated that water flow may also facilitate the fruit dispersal of H. chrysotricha. Although the hydration of H. chrysotricha could remain speculative due to the lack of detailed morphological and experimental evidence, the genetic evidence presented here further verified that the species has sufficient fruit (seed) dispersal capabilities to maintain moderate-to-high levels of ongoing gene flow and population connectivity across fragmented landscapes at large spatial and temporal scales.
Corresponding to the above results, little evidence of strong genetic structuring was found in H. chrysotricha within the two island systems. Based on SNP markers, the STRUC-TURE analysis grouped the TIL populations basically into three clusters ( Figure 1) and the ZA populations into two clusters ( Figure 2). In the TIL region, most of the island populations were genetically separated from western mainland populations, which is consistent with the results of our previous nSSR analysis, whereas the genetic composition of eastern mainland TIL populations was slightly different from the former study, consisting of a mixture of all three clusters [12]. In the ZA region, the genetic divergence between the mainland and island populations was not as high as expected, and most of the island populations still feature a high percentage of genetic admixtures, even after thousands of years of isolation ( Figure 2). The higher levels of population substructure in the TIL region (three clusters) than in the ZA region (two clusters) may be attributed to their contrasting landscape histories. For example, before water flooded the TIL region, mountains (e.g., the Baiji Mountains in the western part of the TIL region) may have acted as a stronger dispersal barrier than water for gene exchange between populations inhabiting the middle and western former hilltops. In the case of the ZA region, which separated from the nearby mainland 7000-9000 years ago due to a rise in sea level, water has been instrumental in facilitating gene exchange between populations for a much longer time than in the TIL region. When combined with the lack of isolation by distance in each island system, these results further suggested that population connectivity in H. chrysotricha has not been greatly modified in either island system, regardless of their significant differences in spatial-temporal characteristics.

Sample Collection and DNA Extraction
We collected a total of 160 H. chrysotricha individuals from eight islands and eight adjacent mainland populations in the TIL region, and 110 individuals from seven islands and four adjacent mainland populations in the ZA region (Tables 1 and 2; Figure 5). In each population, almost equal numbers of individuals with long-styled (L-morph) or short-styled (S-morph) flowers were collected separately. We extracted the total genomic DNA from the dry leaf material of each sample using DNA Plantzol Reagent (Invitrogen, Carlsbad, CA, USA), following the manufacturer's protocol. We then checked the quality of the DNA using 1% agarose gel electrophoresis and measured the DNA concentration using a NanoDrop 2000 spectrophotometer (Thermo Scientific, Wilmington, DE, USA). Finally, 256 of the 270 DNA samples met the minimum quality requirements and were retained for subsequent sequencing and genotyping.

RAD Sequencing and Data Processing
RAD libraries were prepared using the restriction enzyme EcoRI and sample-specific barcodes. Samples in the libraries were pooled and sequenced on an Illumina HiSeq 2500 to generate 150 bp pair-end reads, at the Novogene Bioinformatics Institute (Beijing, China). Standard quality control (QC) pipelines were used to process the raw sequencing data. Reads with > 5% of unknown nucleotides, low-quality reads (quality value ≤ 5), and reads that did not contain sample-specific barcodes were discarded. After filtering procedures, the quality-filtered reads data for each individual were then assembled into de novo loci, and SNP calling was performed using the following core programs: USTACKS, CSTACKS, SSTACKS, GSTACKS, and POPULATIONS, implemented in the STACKS pipeline v2.2 [41]. Briefly, all sequences were firstly processed in USTACKS to produce consensus sequences of RAD tags. The program USTACKS takes a set of short-read sequences as input and aligns them into exactly matching stacks. The minimum depth of coverage required for creating a stack was set at three (-m: 3) sequences, and the maximum distance allowed between stacks was set at four (-M: 4) nucleotides. The program CSTACKS was used to build a catalog of consensus loci containing all the loci from all the H. chrysotricha samples and merging all alleles together, with the number of mismatches allowed between sample loci set to four (-n: 4). Next, sets of stacks (i.e., putative loci) created in USTACKS were compared against the catalog using SSTACKS. Then, GSTACKS was used to align the reads to the locus and to call SNPs. Finally, the catalog of reads was filtered using the POPULATIONS program to produce a data set for downstream analyses. To identify highcredibility SNPs, polymorphic RAD loci that were present in at least 75% of the individuals in each population were retained (-r: 0.75). Potential homologs were excluded by removing loci showing heterozygosity of >0.5 within samples [42], and one SNP per RAD locus was kept, to reduce the impact of linkage disequilibrium.

Genetic Diversity and Population Structure
We calculated the observed heterozygosity (H O ), expected heterozygosity (H E ), nucleotide diversity (π), and inbreeding coefficient (F IS ) for each population using the POP-ULATIONS program in STACKS. Pairwise F-statistics (F ST ) (1000 permutations) among populations were calculated using ARLEQUIN v3.5.2 [43]. Population structure was analyzed using the program STRUCTURE v2.3.5 [44]. The number of genetic clusters (K) was set from 1 to 16 for TIL populations and from 1 to 11 for ZA populations, corresponding to the number of sampled populations in each region. Ten independent iterations were conducted for each K, with a burn-in of 10,000 and 100,000 Markov chain Monte Carlo replicates, by assuming the admixture model and independent allele frequencies. The optimal K value was determined from the ∆K values calculated by STRUCTURE HARVESTER [45]. A Mantel test for genetic differentiation [F ST /(1 − F ST )] against geographic distance (log10 transformed) was performed in GENALEX v6.1 [46] for each region, to test the pattern of isolation by distance (IBD). Statistical significance was determined with 1000 permutations.

Gene Flow and Demographic History Analyses
For each region, we estimated the levels of interpopulation gene flow using the program BAYESASS v3.03 [47]. First, we ran BAYESASS with the default delta values for allelic frequency, migration rates, and inbreeding coefficients. Then, subsequent runs were adjusted with different delta values to ensure that the acceptance rate ranged between 40 and 60% for each parameter [47]. We performed 10 independent runs (1 × 10 7 iterations with a burn-in of 10 6 generations), each with a different initial seed. Model convergence was assessed by the program TRACER v1.5 [48]. In addition, we also used the DIVMIGRATE function from the DIVERSITY package v1.9.90 in R [49] to calculate the relative direction of the gene flow between island and mainland populations, using Nei's G ST method. To test for asymmetric flow (significantly higher in one direction than the other), 95% confidence intervals were calculated from 1000 bootstrap replicates [50].
We used the program BOTTLENECK v1.2.02 [51] to determine whether H. chrysotricha populations underwent significant reductions in effective population size (N e ). We used Wilcoxon's signed rank test, which examines whether populations exhibit a greater level of heterozygosity than predicted in a population at mutation-drift equilibrium, to detect bottlenecks occurring over approximately the last 2-4 N e generations. For each H. chrysotricha population, we performed 10,000 simulations for each of three mutation models (IAM, infinite allele model; SMM, stepwise mutation model; and TPM, two-phase mutation model, with 95% single-step and 5% multi-step mutations). Statistical significance was set at the 0.05 level.

Conclusions
Based on high-resolution SNPs, we compared the genetic response of H. chrysotricha to short-term and long-term habitat fragmentation. The levels of genetic diversity and the population substructure of H. chrysotricha populations were both higher in the TIL region than in the ZA region, while gene flow was higher and less asymmetric among ZA populations. The differences in genetic diversity between island and mainland populations in both regions were not significant, suggesting that genetic erosion of island populations does not increase as the spatial and temporal scales of fragmentation increase. Population connectivity in H. chrysotricha has also not been greatly modified by either short-term or long-term habitat fragmentation. The strong fruit (seed) dispersal capacity of H. chrysotricha, facilitated by the matrix of water in the island system, could buffer against the negative effects of habitat fragmentation and maintain the population connectivity over thousands of years.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/plants11141800/s1. Table S1: Raw statistical data for Hedyotis chrysotricha samples from the Thousand-Island Lake (TIL) and Zhoushan Archipelago (ZA) regions, generated by Illumina sequencing; Table S2: Genotype information for the 185 filtered SNPs from the TIL dataset; Table S3: Genotype information for the 188 filtered SNPs from the ZA dataset; Table S4: Migration rates (the fraction of individuals in population i (pop i) from population j (pop j)) and confidence intervals between populations of H. chrysotricha from the TIL region. Bold values represent the proportion of non-migrants within a population, and italic values indicate asymmetrical migration rates between populations; Table S5: Migration rates (the fraction of individuals in population i (pop i) from population j (pop j)) and confidence intervals between populations of H. chrysotricha from the ZA region. Bold values represent the proportion of non-migrants within a population, and values in italics indicate asymmetrical migration rates between populations; Figure S1: Analysis of isolation by distance based on F ST /(1 − F ST ) and the geographic distance (log10 transformed) in (A) the Thousand-Island Lake (TIL) region and (B) the Zhoushan Archipelago (ZA) region; Figure S2: Histogram of the STRUCTURE analysis in the TIL region for the model with K = 2, 4, and 5. A vertical bar represents a single individual, and each color corresponds to a suggested cluster; Figure S3: Histogram of the STRUCTURE analysis in the ZA region for the model with K = 3-5. A vertical bar represents a single individual, and each color corresponds to a suggested cluster.
Author Contributions: N.Y. and Y.Q. conceived the study; N.Y., R.L. and S.L. collected the samples. N.Y. extracted the DNA samples. N.Y., S.W. and R.L. analyzed the genetic data and wrote the manuscript. H.P.C., R.L. and Y.Q. revised the manuscript. All authors have read and agreed to the published version of the manuscript.