Demographic and genetic characterization of harvested Corbicula fluminea populations

Commercial harvest and biomass–density relationship

First, commercial harvest data (in tons per year) for the period 1950 to 2017 were obtained from the Food and Agriculture Organization (FAO; http://www.fao.org) of the United Nations food balance sheets (2016). To the best of our knowledge this service constitutes the best available source of information for landings’ estimation. Data were aggregated by region (i.e., Asia) and country (i.e., China), fishing area (i.e., inland waters), animal group (i.e., freshwater mollusks) and species (i.e., Corbicula fluminea).

Second, we aimed to characterize C. fluminea fishery in two zones (lake and peripheral habitats) across six sites (two in lake habitats and four in peripheral habitats). Observations from multiple visits were used to estimate the number of harvest boats in each zone. Corbicula fluminea individuals were directly purchased from local fishers between 2018 and 2019. Fishers were identified by their on-site location and obtained clams were bagged, labelled, returned to the laboratory and preserved in 95% ethanol. Following, samples were counted and measured to the nearest 0.1 mm along the longest posterior-anterior axis using Vernier calipers.

In addition, we established a total 84 sampling points covering 28 sampling sections and eight sampling areas, which included (1) lake habitats representative of high harvest pressure, and (2) peripheral habitats representative of low harvest pressure (Fig. 1). Samples were collected at four sampling time points: December 2016 and April, July and October 2017. Two collection methods were used in this study according to habitat features. First, a modified Petersen grab (area of 1/16 m²) was used for sampling in Poyang Lake and the Yangtze River. Second, a quantitative Surber sampler for benthic macroinvertebrates (30 × 30 cm, 500 µm mesh, SN-312, China) was used in mid-size rivers. All samples were collected in triplicate at each sampling site. Retrieved C. fluminea samples were transported to the laboratory at the University of Nanchang in a low temperature incubator. Specimens were counted and weighed with an electronic balance (HANGPING FA1204B) to the nearest 0.1 g.

Environmental data for the 28 sampling sections were referenced from Li et al. (2019). Briefly, dissolved oxygen (DO), hydrogen ions (pH), turbidity (TURB), salinity (Sal) and water temperature (T) were measured with a multiparameter probe (YSI 650MDS, YSI Inc. Yellow Springs, OH, USA). Chlorophyll-a (Chl-a) was measured with a chlorophyll meter (HL - 168C06, China), and the water velocity (V) was measured with a velocity meter (FP111, Global Water Instrumentation, TX, USA).

Data were tested for homogeneity of variance using Levene’s Test and Kolmogorov–Smirnov’s Test for normality. Differences on the size of harvested C. fluminea among habitats were tested using Analysis of Variance (ANOVA) on ln (x) transformed data. No significant differences were found in density or biomass between sampling time points (ANOVA, P > 0.05; results not shown). Hence, the four sampling time points were pooled together as one sample in order to enrich the data used. In addition, the intercepts of biomass-density relationships were compared between habitats by a unifactorial Analysis of covariance (ANCOVA). The ANCOVA tested for the effect of a predictor (habitat; i.e., lake or peripheral) beyond the effects of one known covariate (density).

Sample collection and DNA extraction

Samples of C. fluminea (n = 216) were collected during 2018 and 2019 from the main area of Poyang Lake (PY; 17 specimens), the connected river channel of Poyang Lake (CR; 30 specimens), Ganjiang River (GJ; 30 specimens), Fuhe River (FH; 56 specimens); Xinjiang River (XJ; 44 specimens); and Raohe River (RH; 39 specimens) in Poyang Lake Basin. Tissues of C. fluminea samples were preserved in 95% ethanol until DNA extraction. The TINAamp Marine Animals DNA Kit was used to extract the genomic DNA. DNA was evaluated by spectroscopic methods (Nanodrop 2000, Thermo Scientific) and agarose gel electrophoresis.

We used the cytochrome c oxidase subunit-I (COI) primer LCO22me2 (5′-GGTCAACAAAYCATAARGATATTGG-3′) in combination with a reverse primer HCO700DY2 (5′-TCAGGGTGACCAAAAAAYCA-3′) to screen for reliable amplification using the 216 samples of C. fluminea. The PCR reaction was performed in a 25-µL volume containing 12.5 µL 2x Taq PCR MasterMix (TianGen), 9.5 µL ddH2O, 1.0 µL forward primer, 1.0 µL reverse primer and 1 µL genomic DNA. PCR amplifications were performed with an initial denaturation of 94 °C for 2 min, followed by 35 cycles of 94 °C for 1 min, 50 °C for 1 min and 72 °C for 1 min; a final extension temperature of 72 °C for 7 min. Gel electrophoresis was used to confirm successful amplification. An EZ-10 Spin Column PCR Product Purification Kit (Promega, Madison, WI) was used to purify the PCR products. An ABI 3730XL DNA Analyzer (Applied Biosystems, Carlsbad, CA) was used to sequence the purified DNA.

We used DNASP v5.10 (Librado & Rozas, 2009) to analyze the number of haplotypes (H), haplotype diversity (H d) and nucleotide diversity (π) for C. fluminea population from mtDNA COI sequences. MRBAYES v.3.2.2 (Ronquist et al., 2012) was used to construct a phylogenetic analysis using Bayesian inference methods. The fit of twenty-four models of evolution was selected using MRMODELTEST v.2.2 (Nylander, 2004). NETWORK 4.5 (Leigh & Bryant, 2015) was used to construct a network of COI haplotypes for C. fluminea population with a 95% connection limit. An analysis of molecular variance (AMOVA) was used to test the statistical significance of genetic divergence within and among populations and habitats (lake vs peripheral) based on ARLEQUIN 3.5 (Excoffier, Smouse & Quattro, 1992; Excoffier & Lischer, 2010). Neutrality tests and mismatch distribution was used to search for a signature of population expansion. ARLEQUIN 3.5 (Excoffier & Lischer, 2010) was used to analyze Tajima’s D and Fu’s Fs tests.

Mantel tests (Spearman’s method; Legendre & Legendre, 2012) with 9999 permutations were used to assess the correlations between pairwise water quality parameters matrices and the matrices of density, biomass and genetic diversity. Redundancy analysis (RDA) was used to evaluate variations in density, biomass and genetic diversity in relation to water quality parameters (Ter Braak & Verdonschot, 1995). Monte Carlo permutation tests were used with 499 permutations to assess the significance (P <0.05) of the RDA gradient, and the eigenvalues of the first 2 axes were used to measure their importance (Ter Braak & Verdonschot, 1995). All water quality parameters, density, biomass and genetic diversity were log10(X + 1) transformed to improve their normality before data analysis. CANOCO v 4.5 (Ter Braak & Verdonschot, 1995) was used to perform all the ordinations.

Commercial harvest and biomass–density relationship

Corbicula fluminea s.l. production in Asia increased up to 10,000 tons from 1950 to 1989. Further, it was stable from 1989 until a rapid increase occurred in 2003. Afterwards, production fluctuated between 40,000 and 20,000 tons, with most of this production (ca. 70%) occurring in China.

The total observation number of harvest boats was distributed between 120 and 160 in the Poyang Lake Basin, including lake and peripheral habitats. Lake habitats concentrate more harvest pressure (i.e., 50% of the harvest boats) than peripheral habitats; i.e., half of the harvest boats distributed across the four sites in peripheral habitats. In general, across the six sampling sites within Poyang Lake Basin, the mean harvestable size ± SD for C. fluminea was 19.05 ± 4.95 mm (n = 411). However, most of the captures (i.e., 74%) were ≥ 15 mm. There was not a significant effect of habitat on the size of harvested C. fluminea [F_(1,512) = 1.288, P = 0.257].

C. fluminea density and biomass was higher in lake habitats than in peripheral habitats (Table 1). A significant effect of habitat × density interaction (F_(1,73) = 30.753, p = 0.001) indicated different regression slopes. After accounting for the density effect, the individual of C. fluminea from peripheral habitats was heavier and larger than the individual from lake habitats (Fig. 2).

Population genetics analysis

Twelve mitochondrial COI haplotypes (GenBank accession numbers MN233672-MN233683) were found among 216 sequenced C. fluminea individuals around Poyang Lake Basin (Fig. 1). Lake habitats had the greatest variation, with 7 haplotypes (Table 2; Fig. 1). Haplotype diversity ranged from 0.355 to 0.765, and the collection population in Poyang Lake site (PY) had the greatest haplotype diversity. Nucleotide diversity ranged from 0.004 to 0.016, and the collection population Ganjiang River site (GJ) had the greatest nucleotide diversity (Table 2).

The twelve COI haplotype sequences formed two clades: four haplotypes from PY, GJ, CR, FH, XJ and RH (see Table 1 for code correspondences) formed one well supported clade, and eight haplotypes from PY, GJ, CR, FH, XJ and RH formed another clade (Fig. 3). The haplotype network of mitochondrial COI sequences showed that the most frequent haplotype of C. fluminea (Hap 4) contained 82 individuals and was shared among individuals from PY, GJ, CR, FH, XJ and RH. The analysis of haplotype network showed similar geographic structure with a general lack of geographic resolution as the phylogenetic analysis (Fig. 3).

Based on the analysis of molecular variance (AMOVA), 25.13% of the total genetic variance was found among populations, and the differentiation was significant among populations (F_ST = 0.25, p = 0.001, Table 3). Pairwise F_ST among C. fluminea populations ranged from 0 to 0.40 with generally significant values (Table 4). The total genetic variance found among peripheral populations was 34.33%, and the differentiation was significant among them (F_ST = 0.34, p < 0.0001). Among lake population, the genetic variance was 1%, but not significant (F_ST = 0.01, p = 0.25). Based on the AMOVA results, a total genetic variance of 0.002% was found between lake and peripheral populations, and the differentiation was not significant between them (F_ST = 0.0001, p = 0.47).There was a lack of significance for the mismatch distribution and neutrality tests based on mitochondrial COI sequences (p < 0.01), which indicated that C. fluminea population did not experience recently demographic expansion (Table 2).

Correlation between density, biomass, genetic and physicochemical parameters

Redundancy analysis showed that 93.6% of cumulative percentage variance of the species–environment relation occurred along the first axis, and along four axes comprised 100% (Fig. 4A). 97.0% of cumulative percentage variance of the species–environment relation occurred along the first axis, and along four axes comprised 100% (Fig. 4B). The density and biomass of C. fluminea were correlated with the dissolved oxygen, water velocity, turbidity and water temperature (Fig. 4A). The genetic diversity of C. fluminea was correlated with the dissolved oxygen, water velocity, turbidity, pH, salinity and water temperature (Fig. 4B). Therefore, variations in genetic parameters, density and biomass were highly correlated with water velocity (P < 0.05; Fig. 4). Moreover, water velocity had a significant effect on the density and genetic diversity of C. fluminea based on the Mantel test (P < 0.05; Table 5).

Commercial harvest and biomass–density relationship

Corbicula fluminea s.l. is an important economic resource in Asiatic markets. Among them, the Chinese shellfish fishery is the world’s biggest producer with about 22,500 tons per year in the period 2003–2017. It should be noted, however, that available information of total C. fluminea s.l. production in China is limited by data quality, availability and suitability. For instance, data on C. fluminea s.l. aquaculture production from China—which is of relevance for explaining our results—appear for the first time in the official statistics in 2003.

Here, we discuss demographic and genetic differences among lake and peripheral C. fluminea s.l. populations. From a demographic point of view, the harvest data confirm that C. fluminea s.l. fishery is focused on large individuals (>15 mm). These harvest data are necessary to verify our initial assumptions that a fishery such as that described here is phenotypically non-random. In this regard, we hypothesized that harvest can induce demographic and genetic changes in wild C. fluminea s.l. populations. Partially supporting this assumption is the observation that the stock of C. fluminea s.l. is in decline compared with 20 years ago (Maximum density: 156 ind. m⁻² in 1997–1999 to 30 ind. m⁻² in 2016–2017), which reflects the increasing fishing pressure (Wang et al., 2007). In adition, our results have shown that lake populations—exposed to intense harvest pressure—have lower biomass at all densities than peripheral populations—exposed to lower harvest pressure. Hence, it is reasonable to assume that higher harvest pressure in lake habitats biased populations towards younger ages and smaller sizes. It should be noted, however, that two quantitative sampling methods (modified Petersen grab and Surber sampler) were used in this study which could potentially bias the study results.

Effect of human activity on genetics diversity and structure of C. fluminea

From a genetic point of view, it is worth noting that Poyang Lake populations share common haplotypes which are also the most common in peripheral populations. This fact strengthens our first hypothesis that populations established in the lake are the source for the peripheral populations. An interesting result, is the presence of six different, some of them unique (refer to Fig. 3B), COI haplotypes in Poyang Lake populations, which could be related with repeated introductions from an alternative source; e.g., via aquaculture production. In general, higher genetic variation resulting from introduced populations allows the species to adapt and spread over large territories (Roman & Darling, 2007). In consequence, it could be expected a greater colonizing potential in C. fluminea s.l. population exposed to harvest management.

Even though we find that our results are robust with respect to changes at the population-level, other factors may also account for the population demographic and genetic characteristics. For instance, it should be noted that past drought events related to climate anomalies, water overexploitation or habitat alteration (e.g., Three Gorges Dam) have dramatically affected Poyang Lake biodiversity (Zhang et al., 2015; Wang et al., 2019). Environmental “bottlenecks” associated to drought conditions may create a new wave front where strong genetic drift may be assumed to be the main mechanism of genetic differentiation among populations (Peischl & Excoffier, 2015). In fact, the genetic structure among C. fluminea haplotypes was not well resolved based on the phylogenetic analysis of mitochondrial COI sequences. The higher level of genetic differentiation may be also attributed to geographical isolation as an important factor that affects distribution patterns and genetic structure of the species (Hayes et al., 2008; Lv et al., 2013; Liu et al., 2017).

Effects of environmental changes on density, biomass and genetic parameters

Anthropogenic alterations such as dam construction, water pollution and sand mining at Poyang Lake Basin are important factors responsible of environmental changes (Xiong, Ouyang & Wu, 2012; Jin et al., 2012). It has been noticed that these kinds of environmental changes may affect aquatic species richness and genetic diversity (Watson et al., 2016; Arroyo-Rodrıguez et al., 2013). In fact, our results indicated that the density, biomass and genetic diversity of C. fluminea were correlated with the dissolved oxygen, water velocity, turbidity and water temperature. In addition, these environmental factors have been also related with survival and growth of C. fluminea (McMahon, 1979; Karatayev, Burlakova & Padilla, 2005; Müller & Baur, 2011; Avelar, Neves & Lavrador, 2014; Ferreira-Rodríguez & Pardo, 2016). For example, increase in turbidity may cause siltation potentially increasing mortality rates (Avelar, Neves & Lavrador, 2014). Unusually low winter temperatures can result in massive die-offs of C. fluminea populations (Werner & Rothhaupt, 2008). Increasing temperature is also responsible of increasing metabolic rates but respiration is depleted at temperature up to 30 °C (McMahon, 1979). Müller & Baur (2011) showed that C. fluminea survival in water of 0 °C decreased from 100% to 17.5% with increasing exposure from 4 to 9 weeks. In addition, low river flow and low dissolved oxygen have been related with high mortalities in several macrozoobenthic species, such as C. fluminea (Ilarri et al., 2011; Sousa, Antunes & Guilhermino, 2008).

Relevance for species management in the invaded range

In addition to commercial harvest, harvest management has been also proposed as an effective strategy for the containment, and long-term control of invasive species (Pasko & Goldberg, 2014). Among freshwater animals, C. fluminea s.l. has become one of the most successful invasive species all over the world (Crespo et al., 2015). The species has measurable negative effects on native species and ecosystems (e.g., competition with native mollusks and habitat alteration; Ferreira-Rodríguez et al., 2018; Ferreira-Rodríguez, Iglesias & Pardo, 2019). On this regard, studying invaders in their native range can shed some light in the complex array of interactions that determine its successful establishment and management opportunities in other world regions (Hierro, Maron & Callaway, 2005). Nevertheless, prior the establishment of any management action, there must be a reasonable prospect that the risk of re-invasion is as low as possible. In this regard, haplotype differentiation may help to identify isolated populations where control or eradication measures will be implemented.