Wetlands explain most in the genetic divergence pattern of Oncomelania hupensis

Infection, Genetics and Evolution 27 (2014) 436–444

Contents lists available at ScienceDirect

Infection, Genetics and Evolution journal homepage: www.elsevier.com/locate/meegid

Wetlands explain most in the genetic divergence pattern of Oncomelania hupensis Lu Liang a, Yang Liu b, Jishan Liao c, Peng Gong a,d,⇑ a

Department of Environmental Science, Policy and Management, University of California, Berkeley, CA 94720-3114, USA State Key Laboratory of Biocontrol and College of Ecology and Evolution, Sun Yat-sen University, Guangzhou 510275, China c Department of Biological Sciences, University of Notre Dame, Notre Dame, IN 46556, USA d Ministry of Education Key Laboratory for Earth System Modeling, Center for Earth System Science, Tsinghua University, Beijing 100084, China b

a r t i c l e

i n f o

Article history: Received 22 April 2014 Received in revised form 10 August 2014 Accepted 15 August 2014 Available online 23 August 2014 Keywords: Schistosomiasis japonicum Oncomelania hupensis Genetic divergence Landscape genetics Isolation by distance Spatial autocorrelation

a b s t r a c t Understanding the divergence patterns of hosts could shed lights on the prediction of their parasite transmission. No effort has been devoted to understand the drivers of genetic divergence pattern of Oncomelania hupensis, the only intermediate host of Schistosoma japonicum. Based on a compilation of two O. hupensis gene datasets covering a wide geographic range in China and an array of geographical distance and environmental dissimilarity metrics built from earth observation data and ecological niche modeling, we conducted causal modeling analysis via simple, partial Mantel test and local polynomial fitting to understand the interactions among isolation-by-distance, isolation-by-environment, and genetic divergence. We found that geography contributes more to genetic divergence than environmental isolation, and among all variables involved, wetland showed the strongest correlation with the genetic pairwise distances. These results suggested that in China, O. hupensis dispersal is strongly linked to the distribution of wetlands, and the current divergence pattern of both O. hupensis and schistosomiasis might be altered due to the changed wetland pattern with the accomplishment of the Three Gorges Dam and the South-to-North water transfer project. Ó 2014 Elsevier B.V. All rights reserved.

1. Introduction The study of coevolution between animal hosts and their parasites plays a fundamental role in understanding the mechanism of disease spread. Many emergence and reemergence of medically relevant diseases, such as malaria, schistosomiasis, are caused by parasites that associated with their animal hosts (Davis et al., 1995; Liu et al., 2011). The evolution of host and parasite is usually in a parallel direction, with host populations dispersed or diversified in one direction, the local parasite population should be in turn adapted in a same direction and strength, or become regionally extinct (Davis et al., 1999). Thus, it would be feasible to gain deeper insights into the prediction of parasite evolution in advance by understanding the genetic divergence pattern of the hosts (Pedersen and Babayan, 2011). This is especially necessary when the host species are under extensive medical control, which could potentially lead to higher mutation rate against those going through stable natural selection. ⇑ Corresponding author at: Department of Environmental Science, Policy and Management, University of California, Berkeley, CA 94720-3114, USA. E-mail address: [email protected] (P. Gong). http://dx.doi.org/10.1016/j.meegid.2014.08.012 1567-1348/Ó 2014 Elsevier B.V. All rights reserved.

Patterns of genetic variation in spatially separated host populations are believed to be strongly mediated by ecological processes (Manel et al., 2003), among which, isolation-by-distance (IBD) (Wright, 1943) and isolation-by-environment (IBE) (Wang and Summers, 2010) are considered as two most important ways (Wang et al., 2013). IBD is generally featured by reduced gene flow with enlarged geographic distance or enhanced landscape barriers, whereas IBE contributes to genetic differentiation through local adaption to distinct inhabiting environments (e.g. Liu et al., 2013). Although the idea that genetic divergence might be influenced by IBD and IBE has motivated a growing body of research (e.g. Coulson et al., 2011; Kappes and Haase, 2012), questions like what is the relative importance of IBD and IBE, and what are the key factors in the overall pattern, have not been fully addressed. One major challenge is the lack of the comprehensive dataset collections. Here, we examined those questions with Oncomelania hupensis, a species of tropical freshwater snail belonging to gastropod mollusk, with the aid of earth observation data and the stateof-art landscape genetics approaches. O. hupensis is the unique intermediate host of the human blood fluke Schistosoma japonicum (Gredler 1881). Schistosomiasis is a serious helminth infection that ranks the second severest zoonosis

L. Liang et al. / Infection, Genetics and Evolution 27 (2014) 436–444

only after malaria in the world, with almost 200 million people infected worldwide (Gryseels et al., 2006), and is mainly prevalent in tropical and subtropical zones of East and Southeast Asia. The population divergence of O. hupensis is regarded to be strongly ecological mediated (Shrivastava et al., 2005). For instance, the subspecies O. h. robertsoni inhabits small irrigation channels in mountainous agricultural regions, whereas the nominate subspecies O. h. hupensis prefers lake, river and marshland environments in the downstream Yangtze River basin. There are some published works discussing about how physical conditions or social interactions influence the transmission of disease at a local scale. It is discovered that the seasonal cycling of O. hupensis was governed by temperature and episodic precipitation event with a population model (Remais et al., 2007), and the infection intensity was associated with crop types via statistical methods (Spear et al., 2004). Through a mathematical model, Gurarie and Seto (2009) found that schistosomiasis transmission was highly correlated with the hydrological systems and resident interactions between villages. However, till now, no study has been contributed to a consensus view about the relative roles of different environmental factors in shaping the divergence of O. hupensis at a nationwide scale. In this study, we first compiled a spectrum of geographical and environmental variables, as well as two sets of genetic markers, which were both included to avoid bias caused by various combination rates and effective population sizes (Davis et al., 1995; Zhou et al., 1995, 2008; Wilke et al., 2000). We specifically examined the impacts of geographical and environmental factors on genetic divergence. Further, a casual modeling framework with the exclusion of spatial autocorrelation was implemented to evaluate the impacts of IBD and IBE on gene flow. Finally, we discussed the potential underlying causes of spatial genetic divergence in Mainland China. 2. Materials and methods 2.1. Genetic dataset We searched in GenBank to obtain the most frequently sampled and deposited sequences of mitochondrial cytochrome C oxidase subunit I (CO1) gene and internal transcribed spacer fragments of the non-coding region of ribosomal DNA (ITS). The CO1 gene is praised for providing stable and reproducible identification of differentiation among rissooidean taxa from population to superfamily level (Davis et al., 1998), and the ITS gene is characterized by a relatively fast evolutionary rate (Stothard et al., 1996). To most sequences, the location information in their metadata was missing or only defined to country level. We thus referenced the source articles of all sequences, and retained those with exact geographic location description. Meanwhile, populations with only one sequence were also excluded. Finally, our CO1 dataset included 204 individuals from 29 localities with an average of 7 individuals per locality and 2.9 standard deviation, and ITS dataset included 165 individuals from a different set of 29 localities with 5.6 as the average and 1.4 standard deviation (GenBank accession number, location information, and population abbreviation can be found in Tables S1 and S2). Their spatial distribution is shown in Fig. 1. 2.2. Genetic analysis The DNA sequences of both gene datasets were aligned by CLUSTAL W algorithm (Thompson et al., 1997) in MEGA5.05 (Tamura et al., 2011) and followed by manual correction. The nucleotide substitution models were chosen according to the Akaike Information Criterion in Modeltest 3.7 separately (Posada

437

and Crandall, 1998). The HKY model (Hasegawa et al., 1985) assuming a rate variation across sites according to a gammashaped distribution was chosen for CO1 gene (gamma value: 0.78), and the General Time Reversible model (Tavaré, 1986) assuming a rate variation across sites according to a gammashaped distribution with an estimated proportion of invariant sites was selected for ITS genes (gamma value: 0.25). We constructed genetic diversification pattern using pairwise FST values and tested their significance by 10,000 permutations in Arlequin v3.5 (Excoffier and Lischer, 2010). Finally, we visualized the genetic distances in a spatially explicit way using Gephi (Bastian et al., 2009). 2.3. Building distance and dissimilarity metrics Five distance metrics for IBD model and two dissimilarity metrics for IBE model were constructed. 2.3.1. Distance metrics (i) Geographic separation, known as Euclidean distance, is the simplest IBD model, and serves as a baseline to assess the effect of alternative landscape elements on population differentiation structure. We used the straight line calculation with Hawth’s analysis tool (Beyer, 2004), and spherical distance calculation should be used instead if the scale is above country level (Liang et al., 2010). (ii) Stream path length between two localities was used as a surrogate of stream connectivity, which was considered as a crucial factor in facilitating the long distance dispersal of snails (Shi et al., 2002). The calculation of stream path length between localities was carried out with ArcGIS network analyst toolkit (ESRI, 2011). (iii) Habitat resistance describes the physical interconnection degree between populations. Resistance degree was represented by the least cost path (LCP as abbreviation hereafter) of the niche suitability layer between two localities. The niche suitability layer was built with the following steps. Firstly, we obtained the most comprehensive snail occurrence dataset from the list of 4532 highly endemic villages in China that was issued by the National Office for Schistosomiasis Control, State Council of China, which accounts for 39% of endemic villages, and 95% of the total infected people and animal hosts in China. Secondly, forty layers that are regarded important for the distribution of O. hupensis were rescaled to 1 km resolution for the construction of an ecological niche model (ENM). Nineteen WorldClim bioclim layers were used to reconstruct the bioclimatic environment of studied populations within a 3 km radius buffer. This data series are interpolations of observed monthly mean temperature and precipitation, which depict annual trends, seasonality, and extreme or limiting climatic factors, and are suggested to be more biologically meaningful than annual mean (Hijmans et al., 2005). The physicochemical properties of soil were derived from Harmonized World Soil Database, which was constructed based on the 1:1,000,000 scale soil map of China (FAO/IIASA/ISRIC/ISSCAS/JRC, 2012). Fourteen layers of properties were utilized, including topsoil gravel content, sand fraction, silt fraction, clay fraction, reference bulk density, organic carbon, pH, cation-exchange capacity, base saturation, total exchangeable bases, calcium carbonate, gypsum, sodicity and salinity. Seven topographic indices were calculated with Digital Elevation Model, namely elevation, slope, distance to nearest the wetland, Topography Position Index, Topographic Wetness Index, Stream Power Index, and Sediment Transport Capacity. Thirdly, we employed boosted regression trees, a novel machine learning method to build the ENM. In the parameter tuning step, we first set up an initial search grid for each parameter, and then used area under the Receiver Operating Characteristic curve (AUC) based on 10-fold cross validation as the criterion to select the best match (75% of random selected samples were used for

438

L. Liang et al. / Infection, Genetics and Evolution 27 (2014) 436–444

Fig. 1. Distribution of snail localities for CO1 and ITS genes overlaid on Digital Elevation Model.

model construction, and the rest were for assessment). We bootstrap tested our model for 100 times, and computed the mean and stand standard deviation. The model result with the highest AUC (0.99) was used in habitat resistance calculation (Fig. S1). Finally, based on niche suitability layer, the LCP for each pair locations was quantified with the minimum cumulated cost across all grids. In this study, the calculations of LCP were all conducted with PATHMATRIX (Ray, 2005). (iv) Landscape fragmentation indicates the qualities of certain type of land cover or land use that serves as obstacle or accelerator to dispersal (Sork et al., 1999). We extracted four landscape types (separately cropland, forest, grassland and wetland) from a 1 km resolution land cover product of China in year 2000, which was developed from 30 m Landsat Thematic Mapper scenes by visual interpretation at the scale of 1:100,000 and aggregated to 1 km raster (Liu et al., 2005). Each landscape type was represented by a raster layer, and each cell in the layer indicated the fraction percentage of particular land cover occupancy. And then, a resistance layer was designed for each class by subtracting the fraction percentage map with 100%. Lower value means that less cost is needed to move through this pixel, and vice versa. From all possible paths from one location to another, the route that can minimize the summation of resistance cost was selected, and a resistance matrix containing least cost distances between pairwise locations was derived for each landscape type (Coulon et al., 2004). (v) Slope is used as a representation of topographical barrier. We calculated slope based on the NASA Shuttle Radar Topographic Mission Digital Elevation Model V4 (USGS, 2004), and produced the topographic resistance layer via LCP. 2.3.2. Dissimilarity metrics (vi) Climate niche describes the climate requirements that allow a species to succeed, and most speciation events were proved

to be associated with significant climate niche differentiation occurred (Warren et al., 2008). The 19 WorldClim bioclim layers in the habitat resistance were used again, but high correlations among some of those variables need to be removed from the dissimilarity calculation. Finally, twelve less correlated ones were selected based on pair-wise correlation analysis, separately annual mean temperature, mean diurnal range, isothermality, temperature seasonality, max temperatures of warmest month, min temperatures of coldest month, mean temperature of wettest quarter, mean temperature of driest quarter, annual precipitation, precipitation of wettest month, precipitation of driest month, and precipitation seasonality. Dissimilarity calculation was conducted with a niche overlap package in R (R Core Team, 2013). (vii) Environmental niche dissimilarity concerns about the environment distinction at different localities. Unlike habitat connectivity, it does not consider the environmental gradients between localities over space. Ecological niche dissimilarity was calculated with the same 40 environmental data layers as used in habitat resistance. Although both estimations share the same inputs, they are not inherently correlated. For each locality, the mean value for each variable was extracted and between-locality dissimilarity was calculated. 2.4. Landscape genetic analysis In the landscape genetic analysis, simple Mantel tests were first used to conduct regression between FST matrix and each distance/ dissimilarity matrix. Since IBD and IBE are not mutually exclusive, for instance, the likelihood of detecting environmental dissimilarity usually increases with greater geographic distances, we further implemented partial Mantel tests to exclude the eco-spatial autocorrelation effect from metrics (Smouse et al., 1986). We run both simple and partial Mantel tests in R with the significance test of association assessed by 1,000 permutations.

L. Liang et al. / Infection, Genetics and Evolution 27 (2014) 436–444

Finally, for each pair of comparison, we fitted a local polynomial regression to visualize the IBD and IBE patterns (Chambers and Hastie, 1991), which quantifies our assumption that populations close in geographic or environmental distance tend to be more genetically similar, and a maximum degree of genetic differentiation will be reached after a certain geographic or environmental distance. According to this assumption, the most typical theoretical polynomial curve shapes rise monotonically as a function of distance, and become stabilized afterwards. The curve shape can be characterized by sill, slope and range (Fig. 2). Specifically, the sill is the curve upper bound. The range denotes the distance at which the curve reaches the sill, and slope is the ratio of increment in genetic distance from the starting point to sill to the range. 3. Results 3.1. Genetic analysis of O. hupensis The alignment of the CO1 dataset resulted in a total of 645 characters including 9 sites with gaps and missing data that were excluded from further analyzes. Among the 168 polymorphic sites, 141 were parsimony informative. Pairwise FST values ranged from negative values to one, and 372 of the 406 pairwise CO1 gene comparisons were significantly larger than zero (Table S3). ITS dataset with 165 sequences contained a total number of 968 sites and 72 were parsimony informative among the 81 polymorphic sites. Its pairwise FST values showed a wider genetic differentiation range than CO1 genes, and 306 out of 406 pairwise comparisons results were significantly larger than zero (Table S4). The normalized genetic differentiation among O. hupensis populations with significant test values was spatially visualized in Fig. 3, from which, we can clearly observe that the genetic distances among CO1 populations were smaller than those among ITS populations, since a large portion of CO1 lines were in red or yellow, whereas most ITS lines were in green. One common pattern in both genes is that the genetic distances were larger among the populations in the upstream mountain region than those in the downstream plain region, despite of the smaller geographic distances. It is also evident that the genetic discrepancy between the mountain and floodplain populations was the highest in both cases. 3.2. Landscape genetic analysis In Mantel test, we observed that six significant IBD patterns in CO1 dataset, separately geographic separation, topographic barrier,

439

all landscape fragmentation metrics, and both IBE patterns, were significant. With ITS dataset, there are seven significant IBD correlations, namely geographic separation, stream connectivity, topographic barrier, and all landscape fragmentation metrics, as well as one significant IBE factor – the climate niche dissimilarity (Table 1). All significant correlations were positive. Euclidean distance explains 62% (CO1) and 65% (ITS) of total variance in the genetic differentiation matrix. Factors that achieved a more relevant relationship with genetic distance than purely Euclidean distance would suggest a stronger explanatory power. Thus, by comparing with the baseline scenario, we found that topographic barrier (70% for CO1 and 72% for ITS), cropland (70% for CO1 and 66% for ITS), water (72% for CO1 and ITS), and climate niche dissimilarity (68% for CO1) achieved high correlations. Meanwhile, forest (59% for CO1 and 64% for ITS), grassland (51% for CO1 and 34% for ITS) and environmental niche dissimilarity (17% for CO1) showed weaker support to explain genetic divergence variance. After controlling for spatial autocorrelation via partial Mantel tests, all the IBD correlations were weakened notably for CO1 genes, and some significant relationships even turned to insignificant. None of the partial correlation coefficient surpassed that of the Euclidean distance. When tested against ITS data, the spatial autocorrelation impact was even more drastic, with more than 50% relative difference declination in the correlation coefficient values for all the significant variables (Table 1). In contrast, the effect of spatial autocorrelation on IBE is less prominent with both genes. The correlation of environmental niche dissimilarity with the CO1 genetic distance became even stronger and the probability of partial Mantel r value got improved. Although all the significant IBD or IBE features were generally increasing with the increment of genetic distance for both genes, their polynomial curves can be divided into three groups according to their shapes (Fig. 4): (1) steep ladder, of which the curve is in ladder shape and the slope is large; (2) gradual ladder, of which the curve is in ladder shape but with gradually changed slope; (3) irregular shape curve. Both steep and gradual ladder curves follow an increase trend from the origin to the sill and then level off, but the former group has a faster speed of reaching stabilization than the later one due to the larger slope. We observed more IBD/IBE features in the form of gradual ladder shape for the CO1 sub-dataset, and a larger quantity of steep ladder shaped features were in the ITS sub-dataset. Moreover, metrics in the steep ladder shape for both sub-datasets include wetland, cropland fragmentation and topographical barrier, among which wetland presents the greatest slope, whereas the gradual ladder group solely contains geographic separation.

4. Discussion 4.1. The role of IBD and IBE in genetic divergence pattern

Fig. 2. A theoretical local polynomial regression curve and its characteristics.

Though selfing occasionally occurs in some gastropod species (Chase, 2007), it is rarely reported in the dioecious O. hupensis, and thus leads a quantitative effect on genetic divergence. In this study, we considered that our observed patterns were resulted from abiotic effects such as geography and environment. Both IBD and IBE have been considered as important reasons that can drive genetic divergence of populations, however, disentangling the effects of IBD and IBE has only received attention until recently, and some studies reached contradictory conclusions (Crispo et al., 2006; Surget-Groba et al., 2012; Wang et al., 2013). In our case, IBD plays a more prominent role than IBE in explaining the genetic structure of O. hupensis. On average, significant IBD factors explained about 60% of genetic divergence, and the percentage for IBE is around 40%. By excluding the eco-spatial autocorrelation effect from IBD/IBE via partial

440

L. Liang et al. / Infection, Genetics and Evolution 27 (2014) 436–444

Fig. 3. The spatially explicit visualization of normalized pairwise genetic distances for CO1 (up) and ITS (bottom) population. All distances are rescaled to the range [0, 1] for a standardized legend. Only distances with significant test values are shown. Colors of the lines symbolize the strength of genetic differentiation. Red stands for smallest difference (0), and green represents the largest scaled genetic distance (1).

Table 1 Mantel and partial Mantel test results for CO1 and ITS populations. Category

Isolation by distance Geography separation Stream connectivity Topographical barrier Habitat resistance Landscape fragmentation

Isolation by environment Environmental niche dissimilarity Climate niche dissimilarity

Mantel test

Cropland Forest Grassland Wetland

Partial mantel test

Relative difference

CO1

ITS

CO1

ITS

0.62 (0.001) ⁄⁄ 0.13 (0.925) 0.70 (0.001) ⁄⁄ 0.12 (0.123) 0.70 (0.001) ⁄⁄ 0.59 (0.001) ⁄⁄ 0.51 (0.001) ⁄⁄ 0.72 (0.001) ⁄⁄

0.65 (0.001) ⁄⁄ 0.51 (0.001) ⁄⁄ 0.72 (0.001) ⁄⁄ 0.08 (0.822) 0.66 (0.001) ⁄⁄ 0.64 (0.001) ⁄⁄ 0.34 (0.005) ⁄⁄ 0.72 (0.001) ⁄⁄

0.25 (0.995) 0.43 (0.001) ⁄⁄ 0.24 (0.991) 0.41 (0.001) ⁄⁄ 0.26 (0.004) ⁄⁄ 0.26 (0.007) ⁄⁄ 0.47 (0.001) ⁄⁄

0.23 (0.003) ⁄⁄ 0.34 (0.001) ⁄⁄ 0.21(0.991) 0.26 (0.003) ⁄⁄ 0.16 (0.129) 0.01 (0.470) 0.29 (0.002) ⁄⁄

⁄⁄

0.17 (0.03) 0.68 (0.001)

⁄⁄

0.05 (0.703) 0.49 (0.001) ⁄⁄

0.20 (0.008) 0.41 (0.001)

⁄⁄ ⁄⁄

CO1

ITS

39.29%

52.78%

41.43% 55.93% 49.02% 34.72%

60.61%

59.72%

0.01 (0.458) 0.03 (0.416)

Note: Correlation coefficient values for Mantel tests are given for each variable and the significant associated probabilities (in brackets, with p value less than 0.05) were marked by ⁄⁄. Relative difference was calculated for variables that shown significant relationship in both Mantel and partial Mantel test, and in the formula: (partial mantel r – mantel r)/(mantel r).

L. Liang et al. / Infection, Genetics and Evolution 27 (2014) 436–444

441

Fig. 4. Fitted local polynomial regression curves between genetic distance and various distance and dissimilarity metrics, for CO1 populations (i) and ITS populations (ii). (a) geography separation; (b) stream connectivity; (c) topographical barrier; (d) cropland fragmentation; (e) forest fragmentation; (f) grassland fragmentation; (g) wetland fragmentation; (h) habitat connectivity; (i) climate niche dissimilarity; (j) environmental niche dissimilarity.

Mantel tests, IBD effect is still more dominant than IBE. In the previous studies, IBD was also found to be significant in the O. h. robertsoni populations, and high degree of lineage endemicity was observed in respective phylogroups (Hauswald et al., 2011). One possible explanation is that species with poor dispersal ability cannot facilitate gene flow through self-moving or landscape vectors, like what the species with versatile movement ability can do (Liu et al., 2012), and thus larger geographic separation will probably allocate more barriers on gene flow. 4.2. What explains the most in shaping population divergence pattern? In the polynomial regression as part of the landscape genetic analysis, we observed the steep ladder and gradual ladder as two dominant curve shapes, and certain plausible explanations exist.

Under the genetic context, sill can be interpreted as the genetic distance maximum, and range represents the minimum least cost distance or minimum greatest dissimilarity that are needed by populations to reach the relatively highest degree of isolation among all pairs of localities. Slope is a signal of the strength of IBD or IBE effect, where a larger slope indicates a stronger effect. Thus, variables with a steep ladder shape – two landscape fragmentation factors and topographical barrier – should exert more prominent effects in the population dynamics and evolutionary process of O. hupensis, whereas the gradual ladder shaped feature – geographic separation in this case – exhibits only mild function. This is in congruent with the Mantel test results, which suggests that although geography isolation is significantly correlated with genetic distances, it explains a less portion of variation than topographic barrier, cropland and wetland fragmentation.

442

L. Liang et al. / Infection, Genetics and Evolution 27 (2014) 436–444

Among all the geographic and environmental variables, wetland has the most substantial effect on the genetic diversification pattern of O. hupensis in China, as supported by three pieces of evidence. Firstly, it showed the strongest correlation with genetic distance matrix in the Mantel tests, either been controlled or not controlled for spatial-autocorrelation. Secondly, its polynomial regression fitting curve possessed the largest slope and shortest range, which indicates that the absence of wetlands would bring a higher degree of isolation between snail localities, if the distances or dissimilarity of all variables decrease by the same percentage. Thirdly, the maternally inherited mitochondrial CO1 and bi-parentally inherited nuclear ITS exhibited different levels of genetic divergence among populations, which can be a consequence led by several biological factor, such as mutation (Meirmans and Hedrick, 2011), genetic drift and selection (Ballard and Whitlock, 2004), and sex-specific gene flow (Liu et al., 2012). However, both marker sets show the concordant spatial pattern that the concentration of clusters of more genetically similar localities in the downstream floodplains, which reveals more frequent and stronger gene flow happening there (Fig. 3). In contrast, in the upstream mountain zones, although the geographic distances between localities are relatively smaller, the level of genetic distances is almost as high as that between mountain and plain populations. This is convincing given that the widely spread river–lake complex in the downstream Yangtze River plain could enhance the gene flow by facilitating their movement, whereas the mountain populations have a higher probability of isolation due to smaller area of wetland occupancy. Wetland is unique since its hydric soil supplies snails a combined function of dispersal vector (water) and shelter (mud land). Other landscape variables that can solely provide one function, show less to no significant impact on genetic differentiation pattern. A typical example is river. It seems intuitive that it should be important in facilitating the gene flow of snails because of its strong longitudinal transportation ability. Hauswald et al. (2011) observed clear IBD effects for both marker systems and concluded that local snail dispersal mainly occurred along waterways. However, in our case, the effect is relatively minor or even negative. This might be because the river distance from the Western mountainous zone to the central floodplain is much longer than that local study, and it can hardly lead to any survivals after such long distance transmission (Hang et al., 2011). As another illustration, rice paddies and irrigation ditches are ideal habitats of snails, but they are less efficient in the egg transportation without the water flow. Thus, their contributions to individual dispersal are not as remarkable as wetland. Some other geographic variables also contribute to the genetic structure of O. hupensis. The variation explained by slope ranks second in Mantel test. For species with low movement ability, differences in slope explained some of the genetic variation observed among localities, presumably due to reduced mobility as slope increases (Hauswald et al., 2011). Although IBE factors show less impact on the genetic differentiation, climate niche dissimilarity explained approximately half of the variance. Climate niche is proved to be one reason for micro-evolutionary processes (Sork et al., 1999), and certain critical temperature thresholds exist for the living and breeding of O. hupensis (Zheng et al., 1998). However, climate heterogeneity is not as great as the pattern of other variables across a national scale, and thus climate niche dissimilarity alone is insufficient in explaining the variance in genetic differentiation pattern. Nonetheless, environmental niche dissimilarity that employed the same bioclimate dataset in the climate niche calculation with additional soil and topography layers is far less significant. We are unclear about the reasons. But as revealed from Fig. 4, the overall environmental niche dissimilarity is obviously larger than the climate niche dissimilarity, where most points are

concentrated at the latter half of the x-axis, and thus the IBE pattern was not captured at the beginning part. It is possible that the variation in soil and topography variables among sampling localities is larger than that in climate variables, and resulted in the larger combined dissimilarity. 4.3. Implications for schistosomiasis control The estimated number of people infected with S. japonicum varied between 10.5 and 11.8 million in the mid 1950s, when the People’s Republic of China was just founded (Mao and Shao, 1982), and reduced to 1.52–1.64 million, after an era of China’s national schistosomiasis control program (Chen and Zheng, 1999). However, the re-emergency of this disease was recently found in seven endemic provinces (Zhou et al., 2005), and made the progress towards disease elimination slow because of high reinfection rate (Seto et al., 2011). In the same period, the area of natural wetlands in China shrank by about 33% from 1978 to 2008, among which, the artificial wetlands expanded by about 122% (Niu et al., 2012). We do not have direct evidences at current stage to prove the linkage between schistosomiasis reintroduction and artificial wetland area increase. However, since most artificial wetlands are located around human residential areas, its expansion could highly facilitate the gene exchange of local species with aliens, or the extension of snail habitats to unexplored regions. Moreover, besides urbanization and foreseeable climate change induced wetland changes, two giant hydrological projects in the endemic region, the Three Gorges Dam (TGD) and the South-to-North water transfer projects, should be noticed by their impacts on the hydrological cycle, landform structure, wildlife habitat and migration behavior, and relocation of human residences. It is predicted that the long term TGD will increase the area of permanent marshlands in downstream Yangtze River because of the regulation of flood, and the degeneration of rice paddies into marshlands as a result of underground water levels rising (Gray et al., 2012). These extensive environmental transformations may alter the current genetic diversification pattern of both parasite and host and bring more uncertainties in the prediction of their evolutionary path. The schistosomiasis disease persistence, establishment, and intervention optimization are highly dependent on the genetic diversification pattern of O. hupensis populations, which is demonstrated to be strongly environmentally associated. Landscape features have the potential to inform public health decisions such as where to focus surveillance efforts, or where to disrupt the connection to stop the gene exchange. There have been successful examples in proving the effectiveness of environmental management intervention in the schistosomiasis control (Utzinger et al., 2005; McManus et al., 2011). For the surveillance purpose in China, because of the significant role of wetlands in shaping the genetic diversification pattern, we suggest that first priority should be given to the management of hydrological projects, as well as the elimination of parasite contaminated water bodies, proper dispose of bovine feces in marshland regions to avoid them pollute the nearby rivers, and the disconnection of rivers with high risk of infection. 4.4. Uncertainties and limitations One biggest limitation in this study is the potential decreased power of landscape genetic inference due to the low number of samples. Although all the available georeferenced CO1 and ITS gene data of O. hupensis from GenBank was collected for the mainland China, the sample size is relatively small as compared to the vast geographic coverage of the study region. Thus, we need to acknowledge the power of our test may be affected, and some caution in our interpretation is warranted. However, it is impossible to

L. Liang et al. / Infection, Genetics and Evolution 27 (2014) 436–444

quantitatively assess the uncertainties at the current time since no random, linear or systematic sampling for all O. hupensis populations in this region has been conducted yet with high sample size, levels of polymorphism, and number of molecular marker, which was proved to be effective in identifying the landscape process of population structure (Oyler-McCance et al., 2013). To minimize bias in the study, we analyzed two different types of genetic markers, and their basically concordant results increased the reliability and validity of our conclusions. Meanwhile, although larger samples could decrease the variability of the correlation between landscape pattern and genetic structure, it was indicated that the power of landscape genetic inferences was highly dependent on the number and variability of loci, but was weakly related with the number of samples (Landguth et al., 2012). Another limitation is that we did not assess the uncertainties in the LCP calculation, which were determined by the accuracy of the generated landscape fraction. In our case, the errors can be from two sources: (1) the classification accuracy of the land cover map. Although the reported average accuracy was as high as 92.9% (Liu et al., 2005), the accuracy for cropland, grassland and wetland could be much lower than the mean since greater confusion may exist between these three land cover types and others (Gong et al., 2013; Liang and Gong, 2013). (2) Uncertainties caused by the scaling effect during the aggregation process from 30 m vector data to 1 km grid layer, which had been proved to be significant (Liu et al., 2001). We did not carry out the uncertainty assessment here due to data restriction, and it could be performed by varying the scale and scaling up methods. For instance, a series of raster data sets with 1 km, 2 km, 4 km and 10 km spatial resolution can be derived using nearest neighborhood, maximum area, or centric attribute decision law rules, and be applied in the LCP calculation, in order to observe how those two factors might affect the results of landscape genetic analysis. 5. Conclusions To the best of our knowledge, this is the first study that investigated the ecological reasons for driving the divergence pattern of O. hupensis populations in China. The integration of geospatial analyses and ecological niche modeling revealed that the distribution of wetlands predicted the genetic differentiation pattern of O. hupensis. Yet Isolation-by-distance, rather than environmental factors played a more important role in influencing gene flow of this zoonotic vector. Thus our results provide some insights into the surveillance of its parasite S. japonicum, and the analytical framework applied in this study can be useful to the management of environmentally mediated diseases. Acknowledgement This work was partially supported by a research grant from Tsinghua University (2012Z02287), and Guangdong Natural Science Foundation (S2013040016690). Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.meegid.2014.08. 012. References Ballard, J.W.O., Whitlock, M.C., 2004. The incomplete natural history of mitochondria. Mol. Ecol. 13, 729–744. Bastian, M., Heymann, S., Jacomy, M., 2009. Gephi: an open source software for exploring and manipulating networks. ICWSM, 361–362.

443

Beyer, H.L., 2004. Hawth’s analysis tools for ArcGIS. Chambers, J.M., Hastie, T.J., 1991. Statistical Models in S. CRC Press Inc. Chase, R., 2007. Gastropod reproductive behavior. Scholarpedia 2, 4125. Chen, M.G., Zheng, F., 1999. Schistosomiasis control in China. Parasitol. Int. 48, 11–19. Coulon, A., Cosson, J.F., Angibault, J.M., Cargnelutti, B., Galan, M., Morellet, N., Petit, E., Aulagnier, S., Hewison, A.J., 2004. Landscape connectivity influences gene flow in a roe deer population inhabiting a fragmented landscape: an individualbased approach. Mol. Ecol. 13, 2841–2850. Coulson, T., MacNulty, D.R., Stahler, D.R., Wayne, R.K., Smith, D.W., 2011. Modeling effects of environmental change on wolf population dynamics, trait evolution, and life history. Science 334, 1275–1278. Crispo, E., Bentzen, P., Reznick, D.N., Kinnison, M.T., Hendry, A.P., 2006. The relative influence of natural selection and geography on gene flow in guppies. Mol. Ecol. 15, 49–62. Davis, G.M., Yi, Z., Hua, G.Y., Spolsky, C., 1995. Population-genetics and systematic status of Oncomelania hupensis (Gastropoda, Pomatiopsidae) throughout China. Malacologia 37, 133–156. Davis, G.M., Wilke, T., Spolsky, C., Qiu, C.P., Qiu, D.C., Xia, M.Y., Zhang, Y., Rosenberg, G., 1998. Cytochrome oxidase I-based phylogenetic relationships among the Pomatiopsidae, Hydrobiidae, Rissoidae and Truncatellidae (Gastropoda: Caenogastropoda: Rissoacea). Malacologia 40, 251–266. Davis, G.M., Wilke, T., Zhang, Y., Xu, X.J., Qiu, C.P., Spolsky, C., Qiu, D.C., Li, Y., Xia, M.Y., Feng, Z., 1999. Snail-Schistosoma, Paragonimus interactions in China: population ecology, genetic diversity, coevolution and emerging diseases. Malacologia 41, 355–377. Excoffier, L., Lischer, H.E., 2010. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol. Ecol. Resour. 10, 564–567. ESRI, 2011. ArcGIS Desktop: Release 10. Environmental Systems Research Institute, Redlands, CA. FAO/IIASA/ISRIC/ISSCAS/JRC, 2012. Harmonized World Soil Database (version 1.2). FAO, Rome, Italy and IIASA, Laxenburg, Austria. Gong, P., Wang, J., Yu, L., Zhao, Y.C., Zhao, Y.Y., Liang, L., Niu, Z.G., Huang, X.M., Fu, H.H., Liu, S., Li, C.C., Li, X.Y., Fu, W., Liu, C.X., Xu, Y., Wang, X.Y., Cheng, Q., Hu, L.Y., Yao, W.B., Zhang, H., Zhu, P., Zhao, Z.Y., Zhang, H.Y., Zheng, Y.M., Ji, L.Y., Zhang, Y.W., Chen, H., Yan, A., Guo, J.H., Yu, L., Wang, L., Liu, X.J., Shi, T.T., Zhu, M.H., Chen, Y.L., Yang, G.W., Tang, P., Xu, B., Ciri, C., Clinton, N., Zhu, Z.L., Chen, J., Chen, J., 2013. Finer resolution observation and monitoring of global land cover: first mapping results with Landsat TM and ETM+ data. Int. J. Remote Sens. 34, 2607–2654. Gray, D.J., Thrift, A.P., Williams, G.M., Zheng, F., Li, Y.-S., Guo, J., Chen, H., Wang, T., Xu, X.J., Zhu, R., 2012. Five-year longitudinal assessment of the downstream impact on schistosomiasis transmission following closure of the Three Gorges Dam. PLoS Negl. Trop. Dis. 6, e1588. Gryseels, B., Polman, K., Clerinx, J., Kestens, L., 2006. Human schistosomiasis. Lancet 368, 1106–1118. Gurarie, D., Seto, E.Y., 2009. Connectivity sustains disease transmission in environments with low potential for endemicity: modelling schistosomiasis with hydrologic and social connectivities. J. R. Soc. Interface 6, 495–508. Hang, D.R., Tang, H.P., Huang, Y.X., She, G.S., Zhang, J.F., Huang, Y.J., Li, W., Zhu, X.G., Yang, K., 2011. Impact of simulation operation high water level on Oncomelania hupensis natural growth in water diversion rivers of east route of South-toNorth water Diversion Project. Chin. J. Schisto. Control 23, 664–667. Hasegawa, M., Kishino, H., Yano, T., 1985. Dating of the human-ape splitting by a molecular clock of mitochondrial DNA. J. Mol. Evol. 22, 160–174. Hauswald, A., Remais, J.V., Xiao, N., Davis, G.M., Lu, D., Bale, M.N., Wilke, T., 2011. Stirred, not shaken: genetic structure of the intermediate snail host Oncomelania hupensis robertsoni in an historically endemic schistosomiasis area. Parasit. Vectors 4, 206. http://dx.doi.org/10.1186/1756-3305-4-206. Hijmans, R.J., Cameron, S.E., Parra, J.L., Jones, P.G., Jarvis, A., 2005. Very high resolution interpolated climate surfaces for global land areas. Int. J. Climatol. 25, 1965–1978. Kappes, H., Haase, P., 2012. Slow, but steady: dispersal of freshwater molluscs. Aquat. Sci. 74, 1–14. Landguth, E.L., Fedy, B.C., Oyler-McCance, S., Garey, A.L., Emel, S.L., Mumma, M., Wagner, H.H., Fortin, M.J., Cushman, S.A., 2012. Effects of sample size, number of markers, and allelic richness on the detection of spatial genetic pattern. Mol. Ecol. Resour. 12, 276–284. Liang, L., Xu, B., Chen, Y.L., Liu, Y., Cao, W.C., Fang, L.Q., Feng, L.M., Goodchild, M.F., Gong, P., 2010. Combining spatial-temporal and phylogenetic analysis approaches for improved understanding on global H5N1 transmission. PLoS One 5, e13575. http://dx.doi.org/10.1371/journal.pone.0013575. Liang, L., Gong, P., 2013. Evaluation of global land cover maps for cropland area estimation in the conterminous United States. Int. J. Digital Earth., http:// dx.doi.org/10.1080/17538947.2013.854414. Liu, J.Y., Tian, H.Q., Liu, M.L., Zhuang, D.F., Melillo, J.M., Zhang, Z.X., 2005. China’s changing landscape during the 1990s: large-scale land transformations estimated with satellite data. Geophys. Res. Lett. 32, L02405. http:// dx.doi.org/10.1029/2004GL021649. Liu, M., Tang, X., Liu, J., Zhuang, D., 2001. Research on scaling effect based on 1 km grid cell data. J. Remote Sens. 5, 183–189. Liu, Y., Keller, I., Heckel, G., 2011. Range-wide genetic population structure of common pochard (Aythya ferina): a potentially important vector of highly pathogenic avian influenza viruses. Ecol. Evol. 1, 529–545. Liu, Y., Keller, I., Heckel, G., 2012. Breeding site fidelity and winter admixture in a long-distance migrant, the tufted duck (Aythya fuligula). Heredity 109, 108– 116.

444

L. Liang et al. / Infection, Genetics and Evolution 27 (2014) 436–444

Liu, Y., Webber, S., Bowgen, K., Schmaltz, L., Bradley, K., Halvarsson, P., Abdelgadir, M., Griesser, M., 2013. Environmental factors influence both abundance and genetic diversity in a widespread bird species. Ecol. Evol. 3, 4683–4695. Manel, S., Schwartz, M.K., Luikart, G., Taberlet, P., 2003. Landscape genetics: combining landscape ecology and population genetics. Trends Ecol. Evol. 18, 189–197. Mao, S.P., Shao, B.R., 1982. Schistosomiasis control in the People’s Republic of China. Am. J. Trop. Med. Hyg. 31, 92–99. Meirmans, P.G., Hedrick, P.W., 2011. Assessing population structure: Fst and related measures. Mol. Ecol. Res. 11, 5–18. McManus, D.P., Gray, D.J., Ross, A.G., Williams, G.M., He, H.-B., Li, Y.-S., 2011. Schistosomiasis research in the Dongting lake region and its impact on local and national treatment and control in China. PLoS Negl. Trop. Dis. 5, e1053. http:// dx.doi.org/10.1371/journal.pntd.0001053. Niu, Z.G., Zhang, H.Y., Wang, X.W., Yao, W.B., Zhou, D.M., Zhao, K.Y., Zhao, H., Li, N.N., Huang, H.B., Li, C.C., 2012. Mapping wetland changes in China between 1978 and 2008. Chin. Sci. Bull. 57, 2813–2823. Oyler-McCance, S.J., Fedy, B.C., Landguth, E.L., 2013. Sample design effects in landscape genetics. Conserv. Genet. 14, 275–285. Pedersen, A.B., Babayan, S.A., 2011. Wild immunology. Mol. Ecol. 20, 872–880. Posada, D., Crandall, K.A., 1998. Modeltest: testing the model of DNA substitution. Bioinformatics 14, 817–818. R Core Team, 2013. R: A language and environment for statistical computing. R foundation for Statistical Computing. Ray, N., 2005. PATHMATRIX: a geographical information system tool to compute effective distances among samples. Mol. Ecol. Notes 5, 177–180. Remais, J., Hubbard, A., Wu, Z.S., Spear, R.C., 2007. Weather-driven dynamics of an intermediate host: mechanistic and statistical population modelling of Oncomelania hupensis. J. Appl. Ecol. 44, 781–791. Seto, E.Y.W., Remais, J.V., Carlton, E.J., Wang, S., Liang, S., Brindley, P.J., Qiu, D.C., Spear, R.C., Wang, L.D., Wang, T.P., 2011. Toward sustainable and comprehensive control of schistosomiasis in China: lessons from Sichuan. PLoS Negl. Trop. Dis. 5, e1372. http://dx.doi.org/10.1371/journal.pntd.0001372. Shi, C.H., Wilke, T., Davis, G.M., Xia, M.Y., Qiu, C.P., 2002. Population genetics, microphylogeography, ecology, and susceptibility to schistosome infection of Chinese Oncomelania hupensis hupensis (Gastropoda: Rissooidea: Pomatiopsidae) in the Miao River system. Malacologia 44, 333–347. Shrivastava, J., Qian, B.Z., Mcvean, G., Webster, J.P., 2005. An insight into the genetic variation of Schistosoma japonicum in mainland China using DNA microsatellite markers. Mol. Ecol. 14, 839–849. Smouse, P.E., Long, J.C., Sokal, R.R., 1986. Multiple regression and correlation extensions of the Mantel test of matrix correspondence. Syst. Zool. 35, 627–632. Sork, V.L., Nason, J., Campbell, D.R., Fernandez, J.F., 1999. Landscape approaches to historical and contemporary gene flow in plants. Trends Ecol. Evol. 14, 219–224. Spear, R.C., Seto, E., Liang, S., Birkner, M., Hubbard, A., Greenwood, B., Boelaert, M., Rijal, S., Regmi, S., Singh, R., 2004. Factors influencing the transmission of

Schistosoma japonicum in the mountains of Sichuan province of China. Am. J. Trop. Med. Hyg. 70, 48–58. Stothard, J.R., Hughes, S., Rollinson, D., 1996. Variation within the internal transcribed spacer (ITS) of ribosomal DNA genes of intermediate snail hosts within the genus Bulinus (Gastropoda: Planorbidae). Acta Trop. 61, 19–29. Surget-Groba, Y., Johansson, H., Thorpe, R.S., 2012. Synergy between allopatry and ecology in population differentiation and speciation. Int. J. Ecol. 2012. http:// dx.doi.org/10.1155/2012/273413. Tamura, K., Peterson, D., Peterson, N., Stecher, G., Nei, M., Kumar, S., 2011. MEGA5: molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Mol. Biol. Evol. 28, 2731–2739. Tavaré, S., 1986. Some probabilistic and statistical problems in the analysis of DNA sequences. Lect. Math. Life Sci. 17, 57–86. Thompson, J.D., Gibson, T.J., Plewniak, F., Jeanmougin, F., Higgins, D.G., 1997. The CLUSTAL_X windows interface: flexible strategies for multiple sequence alignment aided by quality analysis tools. Nucleic Acids Res. 25, 4876–4882. USGS, 2004. Shuttle Radar Topography Mission, 1 Arc Second scene SRTM_u03_n008e004, Unfilled Unfinished 2.0, Global Land Cover Facility, University of Maryland, College Park, Maryland, February 2000. Utzinger, J., Zhou, X.N., Chen, M.G., Bergquist, R., 2005. Conquering schistosomiasis in China: the long march. Acta Trop. 96, 69–96. Wang, I.J., Glor, R.E., Losos, J.B., 2013. Quantifying the roles of ecology and geography in spatial genetic divergence. Ecol. Lett. 16, 175–182. Wang, I.J., Summers, K., 2010. Genetic structure is correlated with phenotypic divergence rather than geographic isolation in the highly polymorphic strawberry poison – dart frog. Mol. Ecol. 19, 447–458. Warren, D.L., Glor, R.E., Turelli, M., 2008. Environmental niche equivalency versus conservatism: quantitative approaches to niche evolution. Evolution 62, 2868– 2883. Wilke, T., Davis, G.M., Chen, C.E., Zhou, X.N., Zeng, X.P., Zhang, Y., Spolsky, C.M., 2000. Oncomelania hupensis (Gastropoda: Rissooidea) in eastern China: molecular phylogeny, population structure, and ecology. Acta Trop. 77, 215– 227. Wright, S., 1943. Isolation by distance. Genetics 28, 114. Zheng, Y., Qang, Q., Zhao, G., Zhong, J., Zhang, S., 1998. The function of the overlaying climate data in analysis of Oncomelania snail distribution. Chin. Public Health 14, 724–725. Zhou, X.N., Hong, Q.B., Sun, L.P., Xu, Q., Lu, A.S., Wu, Z., Kristensen, T.K., 1995. Population genetics of Oncomelania spp. in mainland China. I. Genetic variation among populations of Oncomelania spp. Chin. J. Schisto. Control 7, 65–71. Zhou, X.N., Wang, L.Y., Chen, M.G., Wu, X.H., Jiang, Q.W., Chen, X.Y., Zheng, J., Jürg, U., 2005. The public health significance and control of schistosomiasis in China – then and now. Acta Trop. 96, 97–105. Zhou, Y.B., Zhao, G.M., Jiang, Q.W., 2008. Genetic Variability of Schistosoma Japonicum (Katsorada, 1904) Intermediate Hosts Oncomelania hupensis (Gredler, 1881) (Gastropoda: Rissooidea). Ann. Zool. 58, 881–889.