The daily extreme operative temperatures experienced in the field by limpets, gradually increased on both shores from March to June 2013 (Fig. 1a). The daily extreme operative temperatures were significantly different between the two shores (paired-samples t test, t = 4.20, df = 89, P < 0.01), revealing that the northwest-facing shore was a hotter habitat than the southeast-facing shore. During the survey period of 90 days, daily extreme operative temperatures beyond 34 ℃ occurred on 12 days (accounting for 13.3%) in the hotter habitat, but none was observed in the benign habitat (Fig. 1b). Temporal autocorrelation analysis showed that the daily extreme operative temperatures were more predictable (greater temporal autocorrelation) on the northwest-facing shore than that on the southeast-facing shore (Fig. 1c). Therefore, the in situ operative temperature data suggested that there were two thermal environments: the hotter and more predictable (HP) habitat on the northwest-facing shore and the relatively benign and unpredictable (BU) habitat on the southeast-facing shore.
Figure 1. Daily operative temperatures of Siphonaria japonica from March 2013 to June 2013 at study sites. a The 99th percentile of daily operative temperatures in the hot and predictable habitats (HP) (red line) and the relative benign and unpredictable habitat (BU) (blue line). Relative frequencies of daily extreme operative temperature in the HP habitat (b) and BU habitat (c). d Autocorrelation of daily extreme operative temperatures in the HP and BU habitats
Abundance of adult limpets gradually decreased (nested ANOVA, F = 4.86, df = 10, P < 0.01, Table 1) from March 24th 2013 (mean±SEM: HP habitat, 117 ±34/m2; BU habitat, 94±39/m2) to June 18th 2013 (HP habitat, 9±3/m2; BU habitat, 28 ±13/m2) (Fig. 2a). There was a significant difference in abundance between habitats within sampling occasions (nested ANOVA, F = 1.86, df = 11, P = 0.04, Table 1).
df SS MS F P Adult limpet density Among sampling occasions 10 1487.38 148.74 4.86 < 0.01 Among habitats within sampling occasions 11 625.03 56.82 1.86 0.04 Within habitat 308 9422.30 30.60 Total 330 28016.00 Egg ribbon density Among sampling occasions 10 4939.24 493.92 15.35 < 0.01 Among habitats within sampling occasions 11 1429.75 129.98 4.04 < 0.01 Within habitat 308 9907.95 32.17 Total 330 21152.00 Number of egg ribbons per limpet Among sampling occasions 10 197.07 19.71 13.97 < 0.01 Among habitats within sampling occasions 11 118.75 10.80 7.65 < 0.01 Within habitat 308 434.43 1.41 Total 330 851.71
Table 1. Two-factor nested ANOVA to test for variation in the abundance of adult limpets and egg ribbons and the number of egg ribbons per limpet among sampling occasions and among habitats within each sampling occasion
There were significant changes of egg ribbon density among sampling occasions (nested ANOVA, F = 15.35, df = 10, P < 0.01, Table 1). Abundance of egg ribbons initially rose in early April, reached maximum values in the middle of May (mean±SEM: HP habitat, 446±121/m2; BU habitat, 123 ± 55/m2), and disappeared in June. There was a significant difference in egg ribbon density between habitats within sampling occasions (nested ANOVA, F = 4.04, df = 11, P < 0.01, Table 1). The numbers of egg ribbons per limpet among sampling occasions (nested ANOVA, F = 13.97, df = 10, P < 0.01, Table 1; Fig. 2b) and between habitats within sampling occasions (nested ANOVA, F = 7.65, df = 11, P < 0.01, Table 1) were significantly different.
Figure 2. Population dynamics of Siphonaria japonica in the hot and predictable (HP) (circle) and benign and unpredictable (BU) (triangle) habitats at study sites in Dongtou, China. a Adult limpet density. b The number of egg ribbons per limpet. c Shell length of adult limpets. Values are reported as mean±SEM. There are significant differences in adult density, the number of egg ribbons per limpet, and shell length between habitats based on nested ANOVA. Asterisk indicates significant difference (P < 0.05) between HP and BU habitats at the same date
The mean body size was significantly smaller in the HP habitat than in the BU habitat in the most field surveys (Fig. 2c). When pooling all shell length data collected from all sampling occasions from each shore, the mean shell length (± SEM) of limpets in the HP habitat (8.08 ± 0.06 mm) was significantly lower than those in the BU habitat (9.13 ± 0.07 mm) (t test, t = 11.22, df = 1250, P < 0.05).
To evaluate effects of in situ real-time temperatures on population density, we conducted GLM analyses with a negative binomial distribution. It showed that population density significantly decreased with habitat temperature increase (GLMNB, P < 0.01, Table 2; Fig. 3). Population densities between two habitats were significantly different (P < 0.01, Table 2), with population density declining more rapidly in the HP habitat. No significant interaction of temperature and habitat on population density was found (P = 0.32).
Predictors df Deviance Residential df Residential deviance P Temperature 1 30.73 313 541.93 2.96e−08*** Habitat 1 24.47 312 517.46 7.57e−07*** Temperature:habitat 1 0.99 311 516.47 0.32 Negative binomial overdispersion parameter: φ = 25.37 (Siphonaria japonica)
Three asterisks indicate P < 0.01
Table 2. Negative binomial generalized linear model (NBGLM) analysis of deviance of the effects of temperatures and habitats on the numbers of individuals of Siphonaria japonica
There was no significant difference in expression of hsp70 mRNA between the two habitats (t test, t = − 0.04, df = 99, P = 0.94). Levels of hsp70 expression increased with onshore body temperature increase (GLM, χ2 = 219.92, P < 0.01, Table 3), and differed between two habitats (χ2 = 17.61, P < 0.01) (Fig. 4). There was no significant interaction between body temperature and habitat (χ2=1.56, P = 0.21). GLM analysis suggested that changes in hsp70 expression per unit change in body temperature were higher in BU than in HP habitats.
df χ2 P Body temperature 1 219.92 < 2.20e−16*** Habitat 1 17.61 2.72e−05*** Habitat×body temperature 1 1.56 0.21 Three asterisks indicate P < 0.01
Table 3. Analysis of deviance for generalized linear model (GLM) with gamma error distribution, to investigate the effects of habitats and on-shore body temperatures on the relative mRNA expressions of hsp70
Figure 4. On-shore relative mRNA expression of hsp70 of Siphonaria japonica at various body temperatures in the hot and predictable (HP) (circle) and benign and unpredictable (BU) (triangle) habitats. Lines represent best-fitting GLM models with a gamma error distribution, and shaded regions represent 95% confidence intervals. Red and blue lines represent data from HP and BU habitats, respectively
Daily extreme operative temperatures
Adult limpet and egg ribbon abundances
Expression of hsp70 mRNA
Biomimetic loggers (Robolimpets) were used to monitor operative temperatures of limpets, according to the method used by Lima and Wethey (2009), on Dongtou Isle, Zhejiang, China (121°10′ N, 27°51′ E; Fig. 5a). Six robolimpets were deployed on Marth 24th 2013 in the middle intertidal zone across heights that covered the main vertical distribution range of Siphonaria japonica (between 2.5 m and 5.5 m above chart datum), on two shores (northwest facing and southeast facing) in the same bay (three robolimpets on each shore). Loggers recorded temperature every 30 min, with a resolution of 0.06 ℃. During the experiment, damaged or broken Robolimpets were replaced with new ones at the same position, and temperature data were collected every seven days in case of loss of the robolimpets.
Figure 5. Study sites in the southeast-facing and northwest-facing shores, Dongtou Island, Zhejiang, China (a). Spawning of Siphonaria japonica with egg ribbons laid on the rocks (b) and investigation of the abundance of limpets and egg ribbons (c) in the field
The daily extreme operative temperature was obtained from the 99th percentile value for each day, taken from the pooled temperatures for all three loggers at a site (Helmuth et al. 2002). Daily extreme operative temperatures were transformed using the reciprocal transformation in case of heterogeneity, and then a paired-samples t test was used to compare the two groups of temperature using SPSS v22.00 with α = 0.05 (IBM Corp, NY, USA).
To examine the temporal patterns of stressful operative temperature on each shore, a temporal autocorrelation analysis was performed with 'acf' function in the stats package in R v3.5 (R Core Team 2014), according to the method by Helmuth et al. (2006). Temporal autocorrelation reflects the relationship between successive observations (lags) in a time series, i.e., the degree to which historical thermal conditions on one day affect future thermal conditions at different lag times. The autocorrelation in daily extreme temperatures on both shores was calculated for lag periods varying from one day to 14 days.
Abundances of adult limpets and egg ribbons were recorded during daytime low tides approximately weekly from March to June 2013 (11 sampling occasions on each shore). A 5-m-wide transect parallel to the sea was established on each shore, and investigations were carried out from 2.5 to 5.5 m above chart datum. A total of 15 100-point doublestrung quadrats (25 × 25 cm) were placed randomly along each shore. Within each quadrat, the numbers of limpets and egg ribbons of S. japonica were counted. The body size of all limpets was measured with a vernier caliper (0.05 mm), and rock surface temperature where the quadrat was placed was recorded as the real-time temperature.
Because our data were over-dispersed count data with a large number of zeros, which may impede accurate data analysis when using parametric variance analysis and nonparametric Friedman methods, the relationship between the abundance of limpets and environmental factors (real-time temperature and shore) was assessed using a generalized linear model (GLM) fitted with a negative binomial distribution (GLMNB) in R v3.5 (R Core Team 2014) with MASS package (generalized linear model "glm.nb") (Venables and Ripley 2002). The model was fitted using both factors (realtime temperature and shore) with two-way interaction and both factors without interaction respectively, and then the model with the lowest Akaike information criterion (AIC) value was selected.
Abundance data were square root transformed before statistical analyses to improve homogeneity of variances. A two-factor nested ANOVA (shores nested within sampling occasions) was used to analyze temporal variation in the abundance of adult limpets and egg ribbons and the number of egg ribbons per limpet using SPSS v22.00 with α = 0.05 (IBM Corp, NY, USA). Bonferroni-corrected alpha levels were used to assess the significance of ANOVA results.
Approximately 10 individuals (>8 mm in shell length) were collected on each shore during spring low tides (at two-week intervals) from March to June 2013. The sampling area was at least 10 m away from the area used for the abundance investigation. Collections were conducted within 40 min of low tide, to minimize possible physiological variations related to endogenous tidal rhythms (Tomanek and Sanford 2003). The body temperature of individuals was measured with a thermocouple inserted underneath the foot. After body temperature measurement, each limpet was immediately dissected, and the foot muscle was cut into small pieces (~ 0.1 cm2); these were then placed into a 1.5 ml microtube containing 0.5 ml microtube containing RNAlater® solution (Ambion, Applied Biosystems, Austin, USA). Samples were kept at 4 ℃ for ~ 12 h to allow the RNAlater to permeate tissues and stabilize the RNA (Hillyard and Clark 2012) and then were stored at − 80 ℃ until further processing.
Total RNA was extracted from~50 mg of foot muscle for each individual sample using Eastep Reagent Kit following the protocol provided by the manufacturer (Promega, Madison, WI, USA) and quantified using a NanoDrop ND-2000 photometer (Thermo Fisher Scientific, Waltham, MA, USA). Synthesis of cDNA was conducted using total RNA (~5 μg) by reverse transcriptase (RT) reactions with a PrimeScript RT Reagent Kit with gDNA Eraser (Takara, Shiga, Japan). We quantified the mRNA expression level of heat shock protein 70 (hsp70) gene using SYBR Green quantitative PCR following the manufacturer's protocol (Bio-Rad Laboratories, Inc., Hercules, CA, USA). For the normalization of hsp70 gene expression, three genes (18S ribosomal RNA, β-actin, and β-tubulin) which usually have relatively stable expression levels were selected as reference genes, and their expression levels were evaluated based on the expression stability measures (M values) as described by Etschmann et al. (2006). The partial sequences of these four genes were cloned (GenBank accessions: 18S, KX529884; β-actin, KX529885; β-tubulin, KX529886; hsp70, KX529887) and used to design the real-time PCR primers (Table 4) with Beacon Designer 7 software (Premier Biosoft International, Palo Alto, CA, USA). The amplification efficiency of each primer set was assessed by real-time PCR prior to sample analyses. All samples were analyzed in triplicate. Ct (dR) values were calculated with the Bio-Rad CFX Manager™ Software v3.0. Finally, the expression level of hsp70 mRNA was determined relative to the value of three reference genes from a reference individual.
Gene name Gene symbol Function Primers (5′–3′) Heat shock protein 70 hsp70 Molecular chaperone F: CGTTGCTCCTCTGTCTCTTG
Beta actin β-actin Reference gene F: ACCACCTACAACTCCATCAT
Tubulin beta chain β-tubulin Reference gene F: AGAACAAGAACTCATCCTACT
18S ribosomal RNA 18S Reference gene F: TTAGCCACACGAGATTGAG
Table 4. Functions and primers of selected genes
To examine the effects of on-shore body temperature and shore on the expression levels of hsp70 gene, a generalized linear model fitted with a gamma error distribution was employed, using R v3.5 with MASS package following Han et al. (2019). A logistic link function was used in the gamma error distribution. On-shore body temperature, shore, and their interactions were used as explanatory variables.
The study was support by National Natural Science Foundation of China (nos. 41776135, 41976142). Nature Science funds for Distinguished Young Scholars of Fujian Province, China (no. 2017J07003).
JW and YWD: conceived and designed the study; JW and XP: conducted the field surveys; JW: conducted heat shock protein quantification and other data analysis. All the authors contributed to manuscript writing and approved the final version for submission.