Variación espacial y temporal en la talla de Nerita scabricosta en áreas protegidas del Pacífico norte de Costa Rica

Bioestadística General I · Proyecto de investigación · Los Burgaditos

Profesora: Karol Ulate-Naranjo

Resumen

Los burgados (Nerita scabricosta) son caracoles del intermareal rocoso alto que se recolectan como alimento y carnada, por lo que su tamaño puede reflejar la presión de extracción. Este estudio evaluó si el largo y el ancho de la concha difieren entre cuatro áreas silvestres protegidas del Área de Conservación Tempisque (Costa Rica), entre años y épocas climáticas, y en relación con la temperatura, la precipitación y el semáforo del Programa Nacional de Monitoreo Ecológico (PRONAMEC). Se midieron 1 019 individuos en 36 muestreos entre 2022 y 2026. El largo medio fue de 2.08 cm. La Reserva Natural Absoluta Cabo Blanco, donde no se permite la extracción, presentó los burgados más grandes, aunque el tamaño del efecto fue pequeño (ε² = 0.01) y la diferencia no se mantuvo al usar el muestreo como unidad de análisis (p = 0.740). El tamaño varió entre años dentro de cada área, especialmente en la Reserva, donde se observó un efecto grande del año sobre la distribución de longitudes (ε² = 0.48). El ANOVA de dos vías para 2025 y 2026 mostró una interacción significativa entre área y año (η² parcial = 0.16), pero el año coincidió con la época de muestreo, por lo que ambos factores no pudieron separarse. El largo medio no se asoció con la temperatura ni con la precipitación de los 90 días previos, y el ancho siguió un patrón similar. Según el semáforo del PRONAMEC, 26 de 34 muestreos fueron «regular», 4 «bueno», 4 quedaron sin categoría y ninguno fue «malo». En conjunto, el tamaño de los burgados dependió más de las características del sitio y del momento del muestreo que del tipo de manejo, por lo que se recomienda un monitoreo balanceado entre sitios, años y épocas, y calibrar los intervalos del protocolo con datos locales.

Palabras clave: gasterópodos, áreas silvestres protegidas, talla de burgados, intermareal rocoso, PRONAMEC.

Abstract

Burgados (Nerita scabricosta) are snails found in the upper rocky intertidal zone that are collected for food and bait; therefore, their size may reflect harvesting pressure. This study evaluated whether shell length and width differ across four protected wilderness areas within the Tempisque Conservation Area (Costa Rica), across years and seasons, and in relation to temperature, precipitation, and the traffic-light system of the National Ecological Monitoring Program (PRONAMEC). A total of 1,019 individuals were measured during 36 sampling events between 2022 and 2026. The mean length was 2.08 cm. The Cabo Blanco Absolute Nature Reserve, where harvesting is prohibited, had the largest snails, although the effect size was small (ε² = 0.01) and the difference was not significant when using the sampling event as the unit of analysis (p = 0.740). Size varied more between years than between areas: within the Reserve, the year accounted for nearly half of the variation in length (ε² = 0.48). The two-way ANOVA for 2025 and 2026 showed a significant interaction between area and year (partial η² = 0.16), but the year coincided with the sampling season, so the two factors could not be separated. The average length was not associated with temperature or precipitation during the preceding 90 days, and the width followed a similar pattern. According to the PRONAMEC traffic-light system, 26 of 34 samples were rated “fair,” 4 were “good,” 4 were unrated, and none were “poor.” Overall, the size of burgados depended more on site characteristics and the timing of the sampling than on management type; therefore, it is recommended to conduct balanced monitoring across sites, years, and seasons, and to calibrate the protocol intervals using local data.

Keywords: gastropods, protected areas, burgado size, rocky intertidal zone, PRONAMEC.

Introducción

Cuando la marea baja, el estrato rocoso queda al descubierto, convirtiéndose en uno de los ambientes más desafiantes de la costa, puesto que se presenta una mayor exposición a la radiación solar, olas con mayor fuerza de impacto, así como el cambiante nivel del agua durante el día (Kunze et al., 2021; Leeuwis & Gamperl, 2022). Según Baliwe et al. (2022), a pesar de dichas condiciones, en estos hábitats se pueden encontrar caracoles, mejillones, algas, entre otros organismos, que más allá de sostener la vida del litoral, representan un recurso importante para las comunidades costeras (Faria et al., 2024; Gutiérrez-Rubio & Arzola-González, 2024).

Es por ello que el tamaño del cuerpo de los gasterópodos es utilizado como una señal del estado de las poblaciones explotadas. En un área protegida de Sudáfrica, las zonas donde está prohibida la extracción de caracoles de diversas especies de interés comercial, presentaron una conservación en cuanto a las densidades y las estructuras de sus tallas naturales (Baliwe et al., 2022). En las islas Azores, los moluscos de la especie con mayor interés fueron más grandes dentro de un área sin acceso humano, en contraparte, la especie menos valorada no mostró diferencias significativas en su tallaje (Faria et al., 2024). Ambos ejemplos sugieren que estas respuestas dependen de la especie y de las características del lugar de muestreo, a su vez, que el medir la forma de las conchas pueden ayudar a determinar los efectos de la exposición de estos organismos, a la intervención humana prolongada a lo largo de los años (Doyle et al., 2022).

El burgado Nerita scabricosta, es un caracol del intermareal alto, el cual posee una amplia distribución en el Pacífico oriental, donde recientemente se ha documentado una expansión de su rango geográfico asociada a procesos de tropicalización (Zarzyczny et al., 2026). Debido a su tamaño, peso y abundancia es considerado como un recurso pesquero potencial, extrayéndose como alimento y como carnada (Chacón-Monge et al., 2019; Gutiérrez-Rubio & Arzola-González, 2024).

Sin embargo, la talla de estos gasterópodos no depende únicamente de la presión por su recolección. Factores físicos como el oleaje influyen en la estructura de tallas en costas rocosas (Rodríguez-Grimón et al., 2024), al igual que la temperatura, dado que en estos ambientes las rocas pueden superar los 50 °C y elevar la temperatura corporal de los caracoles cerca de ese umbral (Marshall et al., 2023). Para evitar la mortalidad, estos organismos buscan microhábitats más frescos y ajustan su tolerancia térmica de acuerdo con las condiciones locales de cada sitio (Hui et al., 2022; Hui et al., 2025). En adición, las variaciones climáticas pueden alterar el reclutamiento de juveniles; por ejemplo, en Lottia gigantea, se ha observado que las anomalías cálidas incrementan la llegada de nuevos individuos, y con ello, la proporción de individuos pequeños (Nielsen et al., 2025).

Por otra parte, la época del año también puede influir en la talla en direcciones opuestas. Durante la época lluviosa, la mayor disponibilidad de agua reduce el estrés por desecación en el intermareal alto, condiciones que pueden favorecer la actividad de forrajeo y el crecimiento de los gasterópodos, ya que la desecación se asocia con pérdidas de peso y menor supervivencia en caracoles litorínidos (Maia & Troncoso, 2024). No obstante, debido a que N. scabricosta se reproduce exclusivamente en ese período (Collin et al., 2017), la llegada de juveniles podría disminuir la talla media de la población. Por estas razones, al evaluar las variaciones en la talla se deben considerar no solo el nivel de protección del área de estudio, sino también las características locales del sitio, la variabilidad interanual y las condiciones climáticas y estacionales.

En Costa Rica, el Sistema Nacional de Áreas de Conservación (SINAC) cuenta con el Programa Nacional de Monitoreo Ecológico (PRONAMEC), creado con el fin de generar y difundir información científica confiable que respalde las decisiones de manejo de la biodiversidad (Fallas-Flores et al., 2024). Uno de sus indicadores para el litoral rocoso es el largo medio de la concha de los burgados, que se clasifica mediante un sistema de semáforo como bueno, regular o malo según el grado de extracción que sugiere (SINAC, 2016). El Área de Conservación Tempisque (ACT), en el Pacífico Norte, alberga áreas silvestres protegidas con distinto régimen de manejo, desde una reserva natural absoluta en la que no se permite ninguna extracción hasta refugios de vida silvestre y un área marina de manejo con usos regulados, lo que permite comparar poblaciones bajo distintos niveles de protección. Con base en ello, el presente estudio tuvo como objetivo evaluar si el tamaño de la concha de N. scabricosta difiere entre cuatro áreas silvestres protegidas del ACT, entre años y entre épocas, si se asocia con la temperatura y la precipitación previas al muestreo, y cómo se clasifican sus poblaciones según el indicador del PRONAMEC. Se plantearon las hipótesis de que la Reserva Natural Absoluta Cabo Blanco, donde no se permite la extracción, presentaría individuos de mayor tamaño que en las áreas con uso regulado; que el tamaño sería mayor en periodos con mayor precipitación y menor temperatura, debido a un menor estrés por desecación, y que además diferirían entre épocas, sin una dirección definida.

Metodología

Área de estudio: El estudio se realizó en el litoral rocoso de la costa pacífica de la península de Nicoya (Guanacaste y Puntarenas, Costa Rica), dentro del Área de Conservación Tempisque (ACT). Se muestrearon 12 sitios distribuidos en cuatro áreas silvestres protegidas con distinto régimen de manejo: (a) la Reserva Natural Absoluta Cabo Blanco (en adelante, Reserva Cabo Blanco), el área con mayor protección, con los sitios San Miguel, Punta Flor y Playa Cabo Blanco, donde se excluye toda actividad extractiva (Decreto Ley N.º 3763, 1966; Alvarado et al., 2012); (b) el Área Marina de Manejo Cabo Blanco (en adelante, Área Marina Cabo Blanco), con los sitios Cabuya, Tambor y Playa Cuestas, donde la pesca y el turismo están regulados (Decreto Ejecutivo N.º 40442-MINAE, 2017); (c) el Refugio Nacional de Vida Silvestre Camaronal, con los sitios Boca Ora, Camaronal Este y Carrillo (Decreto Ejecutivo N.º 23640-MIRENEM, 1994); y (d) el Refugio Nacional de Vida Silvestre Ostional, con los sitios La India, Mojón 69 y San Juanillo (Ley N.º 9348, 2016). En la Reserva Cabo Blanco la pesca no está permitida (Arias et al., 2015) y su área marina está excluida del Área Marina de Manejo (Decreto Ejecutivo N.º 40442-MINAE, 2017, art. 2); en el Área Marina Cabo Blanco la pesca está permitida con licencia (Decreto Ejecutivo N.º 40442-MINAE, 2017, arts. 7 y 9), y en los refugios de Camaronal y Ostional también está permitida (Arias et al., 2015). En el Área Marina Cabo Blanco habitan comunidades de pesca artesanal (Cabuya, Tambor y Manzanillo), y en el primer kilómetro desde la pleamar ordinaria solo se permite la pesca con cuerda de mano (SINAC, 2017, como se citó en Solís Rivera & Fonseca Borrás, 2018). Esta descripción del régimen de uso no constituye una medición de la intensidad de vigilancia ni de la presión de extracción.

Colecta de datos biológicos: Los muestreos se realizaron en 2022, 2023, 2025 y 2026, durante las épocas seca (diciembre a abril) y lluviosa (mayo a noviembre), con una periodicidad irregular. La colecta siguió el Indicador 2 (tallas de burgados) del Protocolo para el Monitoreo Ecológico de las Playas Rocosas (PRONAMEC; SINAC, 2016). Cada muestreo se realizó durante bajamares de sicigia, entre una hora antes y una hora después del pico de marea baja. Según el protocolo, en cada sitio y fecha se seleccionan al azar 30 individuos de Nerita scabricosta en el intermareal rocoso alto, tomados de distintas agregaciones para no sobrerrepresentar un solo grupo. A cada individuo se le midió el largo y el ancho de la concha con un calibrador Vernier, con una precisión de 0.01 cm. La especie se identificó según el protocolo PRONAMEC, que define este indicador para N. scabricosta, la especie de Nerita dominante y sujeta a extracción en el intermareal alto del Pacífico costarricense (Azofeifa-Solano et al., 2018).

Depuración de la base de datos: La base inicial contenía 1 596 registros. Se excluyeron 328 registros copiados dos veces, 60 registros digitados en milímetros que duplicaban registros ya presentes en centímetros, 64 registros de otras especies o sin medición, cinco registros con largo menor a 0.5 cm, incoherentes con el resto de su muestreo, y 120 registros de cuatro muestreos de 2022 en los que el ancho superó al largo de forma sistemática (razón ancho/largo mediana > 1), lo que indica un criterio de medición distinto. La base final incluyó 1 019 individuos. Además, se anularon 34 valores de ancho y se conservó el largo de esos individuos: 30 de un muestreo en el que el ancho no correspondía a una medición (era exactamente la décima parte del largo) y cuatro en los que el ancho era menor que 0.4 veces el largo, una proporción imposible para la forma de la concha. También se mantuvieron seis individuos aislados con un ancho mayor que el largo. Estos 40 registros (30 + 4 + 6) se marcaron como dudosos para un análisis de sensibilidad. Como los datos fueron recolectados por terceros y no se contó con las hojas de campo, no fue posible verificarlos con los registros originales. La época de cada muestreo se asignó según el mes de colecta.

El número de individuos por muestreo varió entre 1 y 30: 32 de los 36 muestreos tuvieron 30 individuos, dos tuvieron 28 y 25, y dos tuvieron menos de 10 (Mojón 69 en marzo de 2025, con un individuo, y Tambor en febrero de 2026, con cinco). Por la misma razón, no fue posible determinar por qué esos muestreos quedaron incompletos. Todos los individuos se incluyeron en las comparaciones entre áreas, años y épocas, en las que el individuo es la unidad de análisis. Los dos muestreos con menos de 10 individuos se excluyeron solo de los análisis basados en el largo medio de cada muestreo (la clasificación del PRONAMEC y la correlación con el clima), porque su promedio no se consideró representativo.

Obtención y acoplamiento de datos ambientales: El Instituto Meteorológico Nacional suministró registros horarios del periodo 2022–2025 de las estaciones automáticas San José-Pinilla (72-1493), Paquera (72-1594), Santa Rita (72-1895) y ASADA Artola (74-0816) (Instituto Meteorológico Nacional (IMN), 2026). Los registros se resumieron en valores diarios y un día se consideró válido si contaba con al menos 20 horas con datos. La estación Santa Rita se excluyó porque su serie presentó un desfase de fechas creciente respecto a las demás estaciones, detectado mediante correlación cruzada. Cada área se asignó a la estación utilizable más cercana en línea recta: los sitios de Cabo Blanco a Paquera (aproximadamente 30 km) y los de Ostional a San José-Pinilla (25–31 km). Los sitios de Camaronal no se incluyeron en el análisis climático por encontrarse a más de 50 km de cualquier estación utilizable. Para cada muestreo se calcularon la temperatura media del aire y la precipitación acumulada de los 90 días previos (día −90 a día −1), con al menos 80 % de días válidos; la precipitación acumulada se estimó como el promedio diario de los días válidos multiplicado por 90. Se eligió una ventana larga porque el crecimiento de la concha es acumulativo (Frank, 1969) y no existen estimaciones publicadas de la tasa de crecimiento de N. scabricosta (Chacón-Monge et al., 2019). Los registros suministrados abarcaron hasta el 31 de diciembre de 2025, por lo que los muestreos de 2026 no contaron con datos climáticos y no se incluyeron en este análisis. Estas variables representan condiciones climáticas regionales y no el microclima de los sitios.

Análisis descriptivo: Para el largo y el ancho de la concha se calcularon el tamaño de muestra (n), la media, la desviación estándar (DE), el error estándar, la mediana, los cuartiles, los valores mínimo y máximo, el coeficiente de variación (CV), la asimetría y la curtosis. Los cálculos se hicieron por área, año y época, y para las combinaciones de área y año, y de área, año y época. La curtosis se reportó sin centrar, de modo que el valor de referencia de una distribución normal es 3. Para describir la forma de las distribuciones se construyeron tablas de frecuencias del largo. El número de clases se definió con la regla de Sturges (Daniel & Cross, 2013), lo que dio 11 clases de 0.4 cm. Estos mismos intervalos se aplicaron a todas las áreas para que las distribuciones fueran comparables entre sí, y se representaron como histogramas de frecuencia relativa. La variación del largo entre áreas y años se representó con diagramas de caja. Estas medidas se examinaron antes de las pruebas de hipótesis para evaluar la forma de las distribuciones y orientar la elección entre pruebas paramétricas y no paramétricas.

Cada muestreo se clasificó según el ámbito de variación permisible del Indicador 2 del PRONAMEC (SINAC, 2016), que evalúa el largo medio de los individuos medidos en cada sitio y fecha mediante un sistema de semáforo con tres categorías: verde o bueno (> 2.5–4 cm), que según el protocolo indica condiciones aptas y ausencia de extracción; amarillo o regular (1.5–2.5 cm), que sugiere un grado medio de extracción; y rojo o malo (< 1 cm), que sugiere un exceso de extracción. El protocolo no define una categoría para los promedios entre 1 y 1.5 cm, por lo que estos muestreos se registraron como “sin categoría”. Los dos muestreos con menos de 10 individuos no se clasificaron.

Análisis inferencial: El plan de análisis y los criterios para elegir cada prueba de comparación (entre áreas, años y épocas), de asociación con el clima y con el semáforo del PRONAMEC se registraron antes de realizar las pruebas de hipótesis. El ANOVA de dos vías y el análisis con el muestreo como unidad se agregaron después, como análisis complementarios. Todas las pruebas fueron de dos colas y se evaluaron con un nivel de significancia α = 0.05. Para comparar tres o más grupos se ajustó primero un ANOVA de una vía y se evaluaron la normalidad de los residuos (Shapiro-Wilk) y la homogeneidad de varianzas (Levene). Si se cumplían ambos supuestos se usó el ANOVA con comparaciones de Tukey; si los residuos eran normales pero las varianzas diferían, el ANOVA de Welch con pruebas t por pares sin varianza común; y si no había normalidad, la prueba de Kruskal-Wallis con comparaciones por pares de Wilcoxon. En ambos casos las comparaciones por pares se ajustaron con Bonferroni. Para comparar dos grupos se evaluaron la normalidad de cada grupo y la igualdad de varianzas (prueba F); se aplicó la prueba t de Student o de Welch si los datos eran normales, y la de Wilcoxon-Mann-Whitney para muestras independientes si no lo eran. Esta última supone observaciones independientes y aleatorias, una variable continua y distribuciones de forma similar entre grupos; cuando las varianzas difirieron, la diferencia se interpretó como una tendencia de un grupo a presentar valores mayores y no como un desplazamiento de la mediana. El tamaño del efecto se expresó como η² en el ANOVA y como ε² en la prueba de Kruskal-Wallis.

Con este criterio se comparó el largo de la concha entre las cuatro áreas para todo el periodo 2022–2026 y entre años dentro de cada área. Las épocas seca y lluviosa se compararon solo en las combinaciones de área y año muestreadas en ambas épocas: la Reserva Cabo Blanco y el Área Marina Cabo Blanco en 2022, y Ostional en 2025. En las comparaciones entre años (cuatro pruebas) y entre épocas (tres pruebas), los valores de p se ajustaron con Bonferroni. La asociación entre el largo medio de cada muestreo y la temperatura media y la precipitación acumulada de los 90 días previos se evaluó mediante correlación, con el muestreo como unidad de análisis, porque todos los individuos de un muestreo comparten los mismos valores climáticos (n = 18 muestreos con estación asignada y al menos 10 individuos); se aplicó la correlación de Pearson si ambas variables seguían una distribución normal y la de Spearman en caso contrario. La asociación entre el área y la categoría del PRONAMEC se evaluó con la prueba exacta de Fisher (n = 34 muestreos), porque varias frecuencias esperadas eran menores de cinco.

Como análisis de sensibilidad, las comparaciones entre áreas, años y épocas se repitieron sin los 40 registros marcados como dudosos, y la correlación con el clima se repitió sin el sitio Punta Flor, que presentó tallas pequeñas en todos sus muestreos.Como análisis complementario no incluido en el plan metodológico inicial, se ajustó un ANOVA de dos vías con el área, el año y su interacción como factores. Se usaron solo los años 2025 y 2026, los únicos en que se muestrearon las cuatro áreas, de modo que todas las combinaciones de área y año tuvieran datos (n = 691 individuos). Como el diseño fue desbalanceado, se usaron sumas de cuadrados tipo II, y el tamaño del efecto se expresó como η² parcial. Se evaluaron la normalidad de los residuos (Shapiro-Wilk) y la homogeneidad de varianzas entre las ocho combinaciones de área y año (Levene); como las varianzas difirieron, la tabla se repitió con errores estándar robustos a la heterocedasticidad (HC3). La interacción se representó con las medias de cada combinación y su intervalo de confianza de 95 %. Además, para evaluar el efecto de la pseudorreplicación, la comparación entre áreas se repitió con la prueba de Kruskal-Wallis usando el largo medio de cada muestreo como unidad (n = 34 muestreos).

Finalmente, el ancho se analizó como complemento con las mismas pruebas, incluidos la correlación con el clima y el ANOVA de dos vías. Los análisis se realizaron en R versión 4.6.1 (R Core Team, 2026), con los paquetes readxl (Wickham & Bryan, 2026) para leer los datos, writexl (Ooms & Denney, 2026) para exportar los resultados, moments (Komsta & Novomestky, 2022) para la asimetría y la curtosis, y car (Fox & Weisberg, 2019) para la prueba de Levene y el ANOVA de dos vías.

Resultados

En total se analizaron 1 019 individuos de Nerita scabricosta en 36 muestreos realizados entre 2022 y 2026. La Reserva Cabo Blanco presentó el mayor largo medio y la mayor variabilidad, y fue el área con más individuos mayores de 3 cm (29.4 %, frente a 10.9 % en el Área Marina Cabo Blanco y menos de 1 % en Camaronal y Ostional). Camaronal fue el área con menor variabilidad. En el ancho, la Reserva Cabo Blanco también presentó el valor medio y la variabilidad más altos (Tabla 1).

Tabla 1. Estadística descriptiva del largo y ancho de la concha (cm) de Nerita scabricosta en cuatro áreas silvestres protegidas del Área de Conservación Tempisque, Costa Rica, 2022–2026.

Variable	Área	n	Media ± DE	Mediana	Mín–Máx	CV (%)	Asimetría	Curtosis
Largo (cm)	RNA	238	2.25 ± 0.97	2.31	0.60–4.20	43.1	0.04	1.66
Largo (cm)	AMM	275	2.08 ± 0.64	2.00	0.67–4.78	30.9	0.88	4.14
Largo (cm)	Camaronal	270	2.02 ± 0.49	2.10	0.77–3.12	24.2	-0.59	2.74
Largo (cm)	Ostional	236	1.98 ± 0.55	2.00	0.83–3.13	27.7	-0.04	1.92
Largo (cm)	Total	1019	2.08 ± 0.69	2.09	0.60–4.78	33.1	0.40	2.98
Ancho (cm)	RNA	208	1.96 ± 0.73	2.10	0.59–3.41	37.0	-0.27	1.66
Ancho (cm)	AMM	275	1.59 ± 0.48	1.50	0.43–2.99	30.2	0.71	3.45
Ancho (cm)	Camaronal	269	1.64 ± 0.41	1.70	0.56–2.51	24.7	-0.72	2.78
Ancho (cm)	Ostional	233	1.57 ± 0.46	1.60	0.63–2.60	29.5	0.02	2.04
Ancho (cm)	Total	985	1.68 ± 0.54	1.66	0.43–3.41	32.3	0.31	2.58

Nota: RNA: Reserva Natural Absoluta Cabo Blanco; AMM: Área Marina de Manejo Cabo Blanco; Camaronal: Refugio Nacional de Vida Silvestre Camaronal; Ostional: Refugio Nacional de Vida Silvestre Ostional. n: número de individuos; DE: desviación estándar; CV: coeficiente de variación. Curtosis no centrada (valor de referencia de la distribución normal = 3). El n del ancho es menor porque se excluyeron anchos no medidos o con valores imposibles.

Las distribuciones del largo difirieron en forma entre áreas (Fig. 1; Tabla 1). En la Reserva Cabo Blanco la distribución fue aplanada, con un grupo de individuos entre 1 y 1.4 cm y otro entre 2.6 y 3.8 cm. La mayoría de los individuos pequeños de esta área provino del sitio Punta Flor (54 de 68). El Área Marina Cabo Blanco y Camaronal mostraron un solo pico, entre 1.8 y 2.6 cm, con asimetría positiva en el AMM y negativa en Camaronal. Ostional también presentó una distribución aplanada, con individuos repartidos de forma similar entre 1 y 3 cm. En todas las áreas fueron escasos los individuos menores de 1 cm (3.5 % del total).

Fig. 1. Distribución de frecuencias relativas del largo de la concha de Nerita scabricosta en cuatro áreas silvestres protegidas del Área de Conservación Tempisque, Costa Rica, 2022–2026. Intervalos de clase de 0.4 cm (regla de Sturges). RNA: Reserva Natural Absoluta Cabo Blanco; AMM: Área Marina de Manejo Cabo Blanco.

El largo de la concha difirió entre áreas (Kruskal-Wallis, H = 10.34, gl = 3, p = 0.016), aunque el tamaño del efecto fue pequeño (ε² = 0.01; Fig. 2). Se aplicó esta prueba, con un nivel de significancia α = 0.05nza de 95 %, porque los residuos no siguieron una distribución normal (Shapiro-Wilk, W = 0.994, p < 0.001) y las varianzas difirieron entre áreas (Levene, F = 79.57, gl = 3 y 1 015, p < 0.001). La Reserva Cabo Blanco difirió de Camaronal (p = 0.023) y de Ostional (p = 0.020), pero no del Área Marina Cabo Blanco (p = 0.470); esta última, Camaronal y Ostional no difirieron entre sí (p = 1.000). Al excluir los registros dudosos, la diferencia entre áreas se mantuvo y fue mayor (H = 40.36, gl = 3, p < 0.001, ε² = 0.04), y la Reserva Cabo Blanco difirió de las otras tres áreas (p < 0.001).

Ver comandos y salida de R

# Comparación del largo entre áreas (2022–2026)
m_area <- aov(Largo ~ Area, data = bio)
shapiro.test(residuals(m_area))                      # normalidad de residuos

#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  residuals(m_area)
#> W = 0.99361, p-value = 0.0002346

leveneTest(Largo ~ Area, data = bio, center = mean)  # homogeneidad de varianzas

#> Levene's Test for Homogeneity of Variance (center = mean)
#>         Df F value    Pr(>F)    
#> group    3  79.569 < 2.2e-16 ***
#>       1015                      
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

kw_area <- kruskal.test(Largo ~ Area, data = bio); kw_area

#> 
#>  Kruskal-Wallis rank sum test
#> 
#> data:  Largo by Area
#> Kruskal-Wallis chi-squared = 10.341, df = 3, p-value = 0.01588

unname(kw_area$statistic / (nrow(bio) - 1))          # épsilon cuadrado

#> [1] 0.01015839

pairwise.wilcox.test(bio$Largo, bio$Area, p.adjust.method = "bonferroni", exact = FALSE)

#> 
#>  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
#> 
#> data:  bio$Largo and bio$Area 
#> 
#>                 RNA Cabo Blanco AMM Cabo Blanco RNVS Camaronal
#> AMM Cabo Blanco 0.470           -               -             
#> RNVS Camaronal  0.023           1.000           -             
#> RNVS Ostional   0.020           1.000           1.000         
#> 
#> P value adjustment method: bonferroni

# Sensibilidad: sin los 40 registros marcados como dudosos
kruskal.test(Largo ~ Area, data = subset(bio, !Sensibilidad))

#> 
#>  Kruskal-Wallis rank sum test
#> 
#> data:  Largo by Area
#> Kruskal-Wallis chi-squared = 40.356, df = 3, p-value = 8.958e-09

Fig. 2. Largo de la concha de Nerita scabricosta por área silvestre protegida, Área de Conservación Tempisque, Costa Rica, 2022–2026. La línea central indica la mediana, la caja el rango intercuartílico, los bigotes 1.5 veces ese rango y los círculos los valores atípicos. n: número de individuos. RNA: Reserva Natural Absoluta Cabo Blanco; AMM: Área Marina de Manejo Cabo Blanco.

El largo varió entre años dentro de cada área, sobre todo en la Reserva y el Área Marina Cabo Blanco (Fig. 3). En la Reserva Cabo Blanco, la mediana pasó de 1.12 cm en 2023 a 3.15 cm en 2025. con un nivel de significancia α = 0.05nza de 95 %, las diferencias entre años fueron significativas en la Reserva Cabo Blanco, el Área Marina Cabo Blanco y Camaronal, pero no en Ostional. En la Reserva Cabo Blanco y en Camaronal, los residuos no siguieron una distribución normal (Shapiro-Wilk, W = 0.986, p = 0.024, y W = 0.970, p < 0.001) y las varianzas difirieron entre años (Levene, F = 38.88, gl = 3 y 234, p < 0.001, y F = 5.13, gl = 2 y 267, p = 0.007), por lo que se aplicó la prueba de Kruskal-Wallis (Reserva Cabo Blanco: H = 113.47, gl = 3, p < 0.001, ε² = 0.48; Camaronal: H = 12.99, gl = 2, p = 0.006, ε² = 0.05).

En el Área Marina Cabo Blanco no se detectaron desviaciones de la normalidad de los residuos (W = 0.994, p = 0.303), pero las varianzas difirieron (F = 29.43, gl = 3 y 271, p < 0.001), por lo que se usó el ANOVA de Welch, que no requiere varianzas iguales (F = 22.41, gl = 3 y 120.6, p < 0.001, η² = 0.26). En Ostional, ninguno de los dos años siguió una distribución normal (2025: W = 0.955; 2026: W = 0.927; ambos p < 0.001), aunque las varianzas fueron similares (F = 1.32, gl = 145 y 89, p = 0.151); la prueba de Wilcoxon-Mann-Whitney no detectó diferencias (W = 7 512, p = 0.258). En la Reserva Cabo Blanco, todos los años difirieron entre sí excepto 2022 y 2026. En el Área Marina Cabo Blanco, el largo de 2025 fue mayor que el de los demás años (p < 0.001). En Camaronal, el largo de 2026 superó al de 2025 (p < 0.001).

Ver comandos y salida de R

# Largo entre años dentro de cada área (4 pruebas, Bonferroni x 4)
por_area <- split(bio, bio$Area)
for (a in names(por_area)) {
  d <- por_area[[a]]; d$Anio <- factor(d$Anio)
  cat("\n=====", a, "=====\n")
  if (nlevels(d$Anio) > 2) {
    print(shapiro.test(residuals(aov(Largo ~ Anio, data = d))))
    print(leveneTest(Largo ~ Anio, data = d, center = mean))
  } else {
    print(tapply(d$Largo, d$Anio, function(x) shapiro.test(x)$p.value))
    print(var.test(Largo ~ Anio, data = d))
  }
}

#> 
#> ===== RNA Cabo Blanco =====
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  residuals(aov(Largo ~ Anio, data = d))
#> W = 0.98647, p-value = 0.0237
#> 
#> Levene's Test for Homogeneity of Variance (center = mean)
#>        Df F value    Pr(>F)    
#> group   3  38.879 < 2.2e-16 ***
#>       234                      
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> ===== AMM Cabo Blanco =====
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  residuals(aov(Largo ~ Anio, data = d))
#> W = 0.99368, p-value = 0.3025
#> 
#> Levene's Test for Homogeneity of Variance (center = mean)
#>        Df F value    Pr(>F)    
#> group   3  29.426 < 2.2e-16 ***
#>       271                      
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> ===== RNVS Camaronal =====
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  residuals(aov(Largo ~ Anio, data = d))
#> W = 0.96974, p-value = 1.776e-05
#> 
#> Levene's Test for Homogeneity of Variance (center = mean)
#>        Df F value   Pr(>F)   
#> group   2   5.127 0.006533 **
#>       267                    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> ===== RNVS Ostional =====
#>         2025         2026 
#> 1.207055e-04 8.239215e-05 
#> 
#>  F test to compare two variances
#> 
#> data:  Largo by Anio
#> F = 1.3237, num df = 145, denom df = 89, p-value = 0.1511
#> alternative hypothesis: true ratio of variances is not equal to 1
#> 95 percent confidence interval:
#>  0.9017949 1.9100931
#> sample estimates:
#> ratio of variances 
#>           1.323681

oneway.test(Largo ~ factor(Anio), data = por_area[["AMM Cabo Blanco"]], var.equal = FALSE)

#> 
#>  One-way analysis of means (not assuming equal variances)
#> 
#> data:  Largo and factor(Anio)
#> F = 22.413, num df = 3.00, denom df = 120.61, p-value = 1.331e-11

p_anios <- c(
  RNA       = kruskal.test(Largo ~ factor(Anio), data = por_area[["RNA Cabo Blanco"]])$p.value,
  AMM       = oneway.test(Largo ~ factor(Anio), data = por_area[["AMM Cabo Blanco"]], var.equal = FALSE)$p.value,
  Camaronal = kruskal.test(Largo ~ factor(Anio), data = por_area[["RNVS Camaronal"]])$p.value,
  Ostional  = wilcox.test(Largo ~ factor(Anio), data = por_area[["RNVS Ostional"]], exact = FALSE)$p.value)
p.adjust(p_anios, method = "bonferroni")

#>          RNA          AMM    Camaronal     Ostional 
#> 7.865311e-24 5.324126e-11 6.050473e-03 2.582155e-01

Fig. 3. Largo de la concha de Nerita scabricosta por área y año, Área de Conservación Tempisque, Costa Rica, 2022–2026. La línea central indica la mediana, la caja el rango intercuartílico, los bigotes 1.5 veces ese rango y los círculos los valores atípicos. Las posiciones vacías corresponden a años sin muestreo en esa área. RNA: Reserva Natural Absoluta Cabo Blanco; AMM: Área Marina de Manejo Cabo Blanco.

La diferencia entre épocas no fue consistente entre áreas. En los tres casos evaluados, al menos una de las épocas se desvió de la normalidad (Shapiro-Wilk: Reserva Cabo Blanco 2022, época seca, W = 0.915, p < 0.001; Área Marina Cabo Blanco 2022, época seca, W = 0.925, p = 0.001; Ostional 2025, época lluviosa, W = 0.926, p < 0.001), por lo que se aplicó la prueba de Wilcoxon-Mann-Whitney para muestras independientes, con un nivel de significancia α = 0.05nza de 95 %. Las varianzas fueron similares entre épocas en el Área Marina Cabo Blanco (F = 0.58, gl = 59 y 29, p = 0.075) y en Ostional (F = 0.62, gl = 55 y 89, p = 0.060), pero difirieron en la Reserva Cabo Blanco (F = 3.04, gl = 57 y 29, p = 0.002).

En la Reserva Cabo Blanco en 2022, el largo fue mayor en la época lluviosa que en la seca (W = 187, p < 0.001; diferencia de 1.08 cm, IC 95 %: 0.69–1.43); como las varianzas difirieron, este resultado indica que los individuos de la época lluviosa tendieron a ser más grandes, y no estrictamente un desplazamiento de la mediana. En el Área Marina Cabo Blanco en 2022 ocurrió lo contrario (W = 1 450, p < 0.001; diferencia de 0.48 cm, IC 95 %: 0.32–0.63). En Ostional en 2025 no hubo diferencia entre épocas (W = 2 450, p = 1.000). Los resultados de años y épocas no cambiaron al excluir los registros dudosos.

Ver comandos y salida de R

# Largo entre épocas (3 combinaciones área-año, Bonferroni x 3)
casos <- list("RNA 2022"      = subset(bio, Area == "RNA Cabo Blanco" & Anio == 2022),
              "AMM 2022"      = subset(bio, Area == "AMM Cabo Blanco" & Anio == 2022),
              "Ostional 2025" = subset(bio, Area == "RNVS Ostional"   & Anio == 2025))
p_epocas <- sapply(names(casos), function(k) {
  d <- casos[[k]]; d$Epoca <- factor(d$Epoca, levels = c("Seca", "Lluviosa"))
  cat("\n=====", k, "=====\n")
  print(tapply(d$Largo, d$Epoca, function(x) shapiro.test(x)$p.value))
  print(var.test(Largo ~ Epoca, data = d))
  w <- wilcox.test(Largo ~ Epoca, data = d, conf.int = TRUE, exact = FALSE); print(w)
  w$p.value
})

#> 
#> ===== RNA 2022 =====
#>         Seca     Lluviosa 
#> 0.0005995446 0.3350223431 
#> 
#>  F test to compare two variances
#> 
#> data:  Largo by Epoca
#> F = 3.0428, num df = 57, denom df = 29, p-value = 0.001693
#> alternative hypothesis: true ratio of variances is not equal to 1
#> 95 percent confidence interval:
#>  1.547280 5.592256
#> sample estimates:
#> ratio of variances 
#>           3.042758 
#> 
#> 
#>  Wilcoxon rank sum test with continuity correction
#> 
#> data:  Largo by Epoca
#> W = 187, p-value = 1.872e-09
#> alternative hypothesis: true location shift is not equal to 0
#> 95 percent confidence interval:
#>  -1.4299859 -0.6900353
#> sample estimates:
#> difference in location 
#>              -1.080044 
#> 
#> 
#> ===== AMM 2022 =====
#>        Seca    Lluviosa 
#> 0.001247608 0.578786486 
#> 
#>  F test to compare two variances
#> 
#> data:  Largo by Epoca
#> F = 0.5777, num df = 59, denom df = 29, p-value = 0.07486
#> alternative hypothesis: true ratio of variances is not equal to 1
#> 95 percent confidence interval:
#>  0.2945203 1.0565591
#> sample estimates:
#> ratio of variances 
#>          0.5777037 
#> 
#> 
#>  Wilcoxon rank sum test with continuity correction
#> 
#> data:  Largo by Epoca
#> W = 1450, p-value = 2.547e-06
#> alternative hypothesis: true location shift is not equal to 0
#> 95 percent confidence interval:
#>  0.3200539 0.6300577
#> sample estimates:
#> difference in location 
#>               0.479974 
#> 
#> 
#> ===== Ostional 2025 =====
#>         Seca     Lluviosa 
#> 0.4289879323 0.0000733006 
#> 
#>  F test to compare two variances
#> 
#> data:  Largo by Epoca
#> F = 0.6236, num df = 55, denom df = 89, p-value = 0.06042
#> alternative hypothesis: true ratio of variances is not equal to 1
#> 95 percent confidence interval:
#>  0.3913667 1.0212781
#> sample estimates:
#> ratio of variances 
#>          0.6235991 
#> 
#> 
#>  Wilcoxon rank sum test with continuity correction
#> 
#> data:  Largo by Epoca
#> W = 2450.5, p-value = 0.7812
#> alternative hypothesis: true location shift is not equal to 0
#> 95 percent confidence interval:
#>  -0.2169777  0.1879537
#> sample estimates:
#> difference in location 
#>            -0.02382919

p.adjust(p_epocas, method = "bonferroni")

#>      RNA 2022      AMM 2022 Ostional 2025 
#>  5.616677e-09  7.639524e-06  1.000000e+00

Con un nivel de significancia α = 0.05, el largo medio de los muestreos no se asoció con las condiciones climáticas de los 90 días previos (n = 18 muestreos; Fig. A1). Se aplicó la correlación de Spearman porque la temperatura (Shapiro-Wilk, W = 0.864, p = 0.014) y la precipitación (W = 0.848, p = 0.008) no siguieron una distribución normal, aunque el largo medio sí lo hizo (W = 0.956, p = 0.535). La correlación no fue significativa ni con la temperatura (ρ = 0.35, p = 0.150) ni con la precipitación (ρ = 0.34, p = 0.169), y la correlación de Pearson dio resultados equivalentes (temperatura: r = 0.36, IC 95 %: −0.13 a 0.71, p = 0.142; precipitación: r = 0.15, IC 95 %: −0.34 a 0.58, p = 0.555). El resultado se mantuvo al excluir el sitio Punta Flor (temperatura: ρ = 0.26, p = 0.324; precipitación: ρ = 0.20, p = 0.469). En estos muestreos, la temperatura media varió entre 25.9 y 28.3 °C, y la precipitación acumulada, entre 4 y 1 090 mm.

Ver comandos y salida de R

# Correlación entre el largo medio por muestreo y el clima de 90 días (n = 18)
p4 <- as.data.frame(read_excel("inferencial/Inferencial_Nerita.xlsx", sheet = "Datos_P4_muestreos"))
sapply(p4[c("Largo_medio", "T_90d", "P_90d")], function(x) shapiro.test(x)$p.value)

#> Largo_medio       T_90d       P_90d 
#> 0.535446017 0.014300224 0.007758464

cor.test(p4$Largo_medio, p4$T_90d, method = "spearman", exact = FALSE)

#> 
#>  Spearman's rank correlation rho
#> 
#> data:  p4$Largo_medio and p4$T_90d
#> S = 626.29, p-value = 0.1499
#> alternative hypothesis: true rho is not equal to 0
#> sample estimates:
#>       rho 
#> 0.3536719

cor.test(p4$Largo_medio, p4$P_90d, method = "spearman", exact = FALSE)

#> 
#>  Spearman's rank correlation rho
#> 
#> data:  p4$Largo_medio and p4$P_90d
#> S = 640.45, p-value = 0.1687
#> alternative hypothesis: true rho is not equal to 0
#> sample estimates:
#>       rho 
#> 0.3390639

cor.test(p4$Largo_medio, p4$T_90d, method = "pearson")

#> 
#>  Pearson's product-moment correlation
#> 
#> data:  p4$Largo_medio and p4$T_90d
#> t = 1.5466, df = 16, p-value = 0.1415
#> alternative hypothesis: true correlation is not equal to 0
#> 95 percent confidence interval:
#>  -0.1277385  0.7082585
#> sample estimates:
#>       cor 
#> 0.3606391

cor.test(p4$Largo_medio, p4$P_90d, method = "pearson")

#> 
#>  Pearson's product-moment correlation
#> 
#> data:  p4$Largo_medio and p4$P_90d
#> t = 0.60288, df = 16, p-value = 0.555
#> alternative hypothesis: true correlation is not equal to 0
#> 95 percent confidence interval:
#>  -0.3416027  0.5758392
#> sample estimates:
#>       cor 
#> 0.1490363

# Sensibilidad: sin Punta Flor
sin_pf <- subset(p4, Sitio != "Punta Flor")
cor.test(sin_pf$Largo_medio, sin_pf$T_90d, method = "spearman", exact = FALSE)

#> 
#>  Spearman's rank correlation rho
#> 
#> data:  sin_pf$Largo_medio and sin_pf$T_90d
#> S = 500.74, p-value = 0.3239
#> alternative hypothesis: true rho is not equal to 0
#> sample estimates:
#>       rho 
#> 0.2636233

cor.test(sin_pf$Largo_medio, sin_pf$P_90d, method = "spearman", exact = FALSE)

#> 
#>  Spearman's rank correlation rho
#> 
#> data:  sin_pf$Largo_medio and sin_pf$P_90d
#> S = 547.22, p-value = 0.4686
#> alternative hypothesis: true rho is not equal to 0
#> sample estimates:
#>       rho 
#> 0.1952697

La categoría del PRONAMEC se asoció con el área (prueba exacta de Fisher, p = 0.015; n = 34 muestreos). Se usó esta prueba en lugar de χ² porque varias frecuencias esperadas fueron menores de cinco. Tres de los cuatro muestreos en la categoría verde (bueno) correspondieron a la Reserva Cabo Blanco, y los muestreos sin categoría se presentaron en la Reserva Cabo Blanco (Punta Flor) y en Ostional (Mojón 69). En conjunto, la mayoría de los muestreos se ubicó en la categoría amarilla (regular) y ninguno en la roja (malo) (Tabla 2).

Ver comandos y salida de R

# Clasificación PRONAMEC por muestreo (n >= 10) y prueba exacta de Fisher
bio$Muestreo <- paste(bio$Sitio, as.Date(bio$Fecha))
mue <- do.call(data.frame, aggregate(Largo ~ Muestreo + Area, data = bio,
                                     FUN = function(x) c(n = length(x), media = mean(x))))
names(mue)[3:4] <- c("n", "Largo_medio")
mue <- subset(mue, n >= 10)
mue$Categoria <- cut(mue$Largo_medio, c(-Inf, 1, 1.5, 2.5, 4, Inf), right = FALSE,
                     labels = c("Malo", "Sin categoría", "Regular", "Bueno", "Sobre el ámbito"))
tab_p5 <- table(mue$Area, droplevels(mue$Categoria)); tab_p5

#>                  
#>                   Sin categoría Regular Bueno
#>   RNA Cabo Blanco             2       3     3
#>   AMM Cabo Blanco             0       8     1
#>   RNVS Camaronal              0       9     0
#>   RNVS Ostional               2       6     0

fisher.test(tab_p5)

#> 
#>  Fisher's Exact Test for Count Data
#> 
#> data:  tab_p5
#> p-value = 0.01494
#> alternative hypothesis: two.sided

Tabla 2. Número de muestreos de Nerita scabricosta por área silvestre protegida según el ámbito de variación permisible del largo de la concha del PRONAMEC, Área de Conservación Tempisque, Costa Rica, 2022–2026.

Área	Bueno	Regular	Sin categoría	Malo	Total
RNA	3	3	2	0	8
AMM	1	8	0	0	9
Camaronal	0	9	0	0	9
Ostional	0	6	2	0	8
Total	4	26	4	0	34

Nota. Cada muestreo corresponde a un sitio en una fecha y se clasificó según el largo medio de sus individuos (SINAC, 2016, Indicador 2): bueno, > 2.5–4 cm; regular, 1.5–2.5 cm; malo, < 1 cm. Sin categoría: largo medio entre 1 y 1.5 cm, intervalo no definido en el protocolo. Se excluyeron dos muestreos con menos de 10 individuos. RNA: Reserva Natural Absoluta Cabo Blanco; AMM: Área Marina de Manejo Cabo Blanco; Camaronal: Refugio Nacional de Vida Silvestre Camaronal; Ostional: Refugio Nacional de Vida Silvestre Ostional.

En el ANOVA de dos vías para 2025 y 2026 (n = 691), la interacción entre el área y el año fue significativa (F = 42.46, gl = 3 y 683, p < 0.001, η² parcial = 0.16), lo que indica que el cambio en el largo entre años no fue el mismo en todas las áreas (Fig. 4). De 2025 a 2026, el largo medio disminuyó en la Reserva Cabo Blanco (de 3.27 a 2.13 cm) y en el Área Marina Cabo Blanco (de 2.54 a 1.85 cm), aumentó en Camaronal (de 1.87 a 2.18 cm) y varió poco en Ostional (de 2.03 a 1.89 cm). El área (F = 47.27, gl = 3 y 683, p < 0.001, η² parcial = 0.17) y el año (F = 49.44, gl = 1 y 683, p < 0.001, η² parcial = 0.07) también fueron significativos, pero, por la interacción, sus efectos no deben interpretarse por separado. Los residuos se desviaron levemente de la normalidad (Shapiro-Wilk, W = 0.994, p = 0.010) y las varianzas difirieron entre combinaciones (Levene, F = 14.45, gl = 7 y 683, p < 0.001); con errores robustos a la heterocedasticidad, los tres términos se mantuvieron significativos (interacción: F = 41.11, gl = 3 y 683, p < 0.001). Por último, con el largo medio de cada muestreo como unidad, la diferencia entre áreas no fue significativa (Kruskal-Wallis, H = 1.26, gl = 3, p = 0.740; n = 34 muestreos).

Ver comandos y salida de R

# ANOVA de dos vías Área x Año (2025-2026)
d2 <- subset(bio, Anio %in% c(2025, 2026)); d2$Anio <- factor(d2$Anio)
table(d2$Area, d2$Anio)

#>                  
#>                   2025 2026
#>   RNA Cabo Blanco   60   60
#>   AMM Cabo Blanco   90   65
#>   RNVS Camaronal    90   90
#>   RNVS Ostional    146   90

modelo2 <- aov(Largo ~ Area * Anio, data = d2)
Anova(modelo2, type = 2)                              # sumas de cuadrados tipo II

#> Anova Table (Type II tests)
#> 
#> Response: Largo
#>            Sum Sq  Df F value    Pr(>F)    
#> Area       49.836   3  47.274 < 2.2e-16 ***
#> Anio       17.373   1  49.441 4.971e-12 ***
#> Area:Anio  44.759   3  42.458 < 2.2e-16 ***
#> Residuals 240.006 683                      
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

shapiro.test(residuals(modelo2))

#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  residuals(modelo2)
#> W = 0.99423, p-value = 0.009925

leveneTest(Largo ~ Area * Anio, data = d2)

#> Levene's Test for Homogeneity of Variance (center = median)
#>        Df F value    Pr(>F)    
#> group   7  14.454 < 2.2e-16 ***
#>       683                      
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova(lm(Largo ~ Area * Anio, data = d2), type = 2, white.adjust = "hc3")  # errores robustos

#> Analysis of Deviance Table (Type II tests)
#> 
#> Response: Largo
#>            Df       F    Pr(>F)    
#> Area        3 101.751 < 2.2e-16 ***
#> Anio        1  19.060 1.463e-05 ***
#> Area:Anio   3  41.106 < 2.2e-16 ***
#> Residuals 683                      
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

interaction.plot(d2$Anio, d2$Area, d2$Largo, xlab = "Año", ylab = "Largo medio (cm)",
                 trace.label = "Área")

# Kruskal-Wallis con el muestreo como unidad (pseudorreplicación)
kruskal.test(Largo_medio ~ Area, data = mue)

#> 
#>  Kruskal-Wallis rank sum test
#> 
#> data:  Largo_medio by Area
#> Kruskal-Wallis chi-squared = 1.2556, df = 3, p-value = 0.7397

Fig. 4. Largo medio de la concha de Nerita scabricosta (± IC 95 %) por área silvestre protegida y año, Área de Conservación Tempisque, Costa Rica, 2025–2026 (n = 691 individuos). Las líneas unen las medias de cada área para mostrar la interacción. RNA: Reserva Natural Absoluta Cabo Blanco; AMM: Área Marina de Manejo Cabo Blanco.

El ancho mostró un patrón similar al del largo. Difirió entre áreas (Kruskal-Wallis, H = 52.98, gl = 3, p < 0.001, ε² = 0.05; n = 985); a diferencia del largo, la Reserva Cabo Blanco difirió de las otras tres áreas (p < 0.001), y el Área Marina Cabo Blanco difirió de Camaronal (p = 0.007). Entre años, el ancho varió en la Reserva Cabo Blanco (Kruskal-Wallis, H = 76.61, gl = 2, p < 0.001, ε² = 0.37), el Área Marina Cabo Blanco (ANOVA de Welch, F = 27.71, gl = 3 y 116, p < 0.001, η² = 0.30) y Camaronal (Kruskal-Wallis, H = 10.89, gl = 2, p = 0.017, ε² = 0.04), pero no en Ostional (Wilcoxon-Mann-Whitney, W = 6 622, p = 1.000). Entre épocas, fue mayor en la época lluviosa en la Reserva Cabo Blanco en 2022 (W = 194, p < 0.001), mayor en la seca en el Área Marina Cabo Blanco en 2022 (t de Student, t = 5.27, gl = 88, p < 0.001) y no difirió en Ostional en 2025 (W = 2 450, p = 1.000); los valores de p entre años y épocas están ajustados con Bonferroni. El ancho medio de los muestreos no se correlacionó con la temperatura (Spearman, ρ = 0.26, p = 0.308) ni con la precipitación (ρ = 0.27, p = 0.300; n = 17 muestreos, porque uno no tuvo valores de ancho válidos).

En el ANOVA de dos vías para 2025 y 2026 (n = 687), la interacción entre el área y el año también fue significativa (F = 46.72, gl = 3 y 679, p < 0.001, η² parcial = 0.17), y se mantuvo con errores robustos a la heterocedasticidad (F = 44.57, p < 0.001); los residuos no se desviaron de la normalidad (W = 0.996, p = 0.051), pero las varianzas difirieron (Levene, F = 17.27, gl = 7 y 679, p < 0.001). Con el ancho medio de cada muestreo como unidad, tampoco hubo diferencia entre áreas (Kruskal-Wallis, H = 3.80, gl = 3, p = 0.284; n = 33 muestreos). En resumen, el ancho llevó a las mismas conclusiones que el largo en casi todas las comparaciones. La diferencia estuvo en las comparaciones entre áreas: la Reserva Cabo Blanco se distinguió de las otras tres áreas, incluida el Área Marina Cabo Blanco, y el Área Marina Cabo Blanco difirió de Camaronal.

Ver comandos y salida de R

# Resultados del ancho, generados con los scripts 02 y 03
res <- as.data.frame(read_excel("inferencial/Inferencial_Nerita.xlsx", sheet = "Resumen_pruebas"))
kable(subset(res, Variable == "Ancho"), caption = "Pruebas de hipótesis para el ancho (script 02).")

Pruebas de hipótesis para el ancho (script 02).

	Pregunta	Grupo	Variable	Factor	n	Niveles	Prueba	Estadistico	gl	p	Efecto	p_Bonferroni
3	P1 Área (ancho, complemento)	Todos	Ancho	Area	985	RNA, AMM, Camaronal, Ostional	Kruskal-Wallis	H = 52.98	3	0.000000	ε² = 0.054	NA
12	P2 Año (ancho, complemento)	RNA	Ancho	Anio	208	2022, 2025, 2026	Kruskal-Wallis	H = 76.61	2	0.000000	ε² = 0.370	0.00e+00
13	P2 Año (ancho, complemento)	AMM	Ancho	Anio	275	2022, 2023, 2025, 2026	ANOVA de Welch	F = 27.71	3, 116	0.000000	η² = 0.296	0.00e+00
14	P2 Año (ancho, complemento)	Camaronal	Ancho	Anio	269	2022, 2025, 2026	Kruskal-Wallis	H = 10.89	2	0.004310	ε² = 0.041	1.72e-02
15	P2 Año (ancho, complemento)	Ostional	Ancho	Anio	233	2025, 2026	Wilcoxon (Mann-Whitney)	W = 6622	NA	0.669000	Dif. de ubicación (Hodges-Lehmann) = 0.04 cm (IC 95 %: -0.12 a 0.18)	1.00e+00
22	P3 Época (ancho, complemento)	RNA 2022	Ancho	Epoca	88	Seca, Lluviosa	Wilcoxon (Mann-Whitney)	W = 194	NA	0.000000	Dif. de ubicación (Hodges-Lehmann) = -0.89 cm (IC 95 %: -1.10 a -0.61)	0.00e+00
23	P3 Época (ancho, complemento)	AMM 2022	Ancho	Epoca	90	Seca, Lluviosa	t de Student	t = 5.27	88	0.000001	Dif. de medias = 0.32 cm (IC 95 %: 0.20 a 0.44)	2.90e-06
24	P3 Época (ancho, complemento)	Ostional 2025	Ancho	Epoca	144	Seca, Lluviosa	Wilcoxon (Mann-Whitney)	W = 2450	NA	0.938000	Dif. de ubicación (Hodges-Lehmann) = 0.00 cm (IC 95 %: -0.14 a 0.17)	1.00e+00

np <- "inferencial/No_planificados_Nerita.xlsx"
kable(read_excel(np, sheet = "A3_ancho_clima"), caption = "Correlación entre el ancho medio y el clima (script 03).")

Correlación entre el ancho medio y el clima (script 03).

Variable	n	rho	p_spearman	r	IC95	p_pearson
T_90d	17	0.26	0.308	0.30	-0.21 a 0.68	0.235
P_90d	17	0.27	0.300	0.09	-0.41 a 0.55	0.727

kable(read_excel(np, sheet = "A3_ancho_ANOVA"), caption = "ANOVA de dos vías para el ancho (script 03).")

ANOVA de dos vías para el ancho (script 03).

Termino	gl	F	p	eta2_parcial	F_HC3	p_HC3
Area	3	47.66	0	0.174	144.08	0.0e+00
Anio	1	44.29	0	0.061	21.56	4.1e-06
Area:Anio	3	46.72	0	0.171	44.57	0.0e+00
Residuals	679	NA	NA	NA	NA	NA

Discusión y conclusiones

La protección absoluta de la Reserva Cabo Blanco se asocia con burgados más grandes, pero esa señal es débil frente a la variación natural entre sitios. La reserva, el área con mayor protección, presentó el mayor largo medio y la mayor proporción de individuos mayores de 3 cm, un patrón esperable donde no se permite la extracción: en zonas sin extracción de otras costas rocosas se han registrado individuos más grandes y abundantes de las especies explotadas, aunque ese efecto dependió del tipo de sustrato y de la región (Baliwe et al., 2022). Nerita scabricosta se extrae en la pesca de pequeña escala en Costa Rica (Azofeifa-Solano et al., 2018) y como alimento y carnada sin regulación en México (Gutiérrez-Rubio & Arzola-González, 2024), por lo que la extracción es una explicación plausible. Sin embargo, la evidencia de este estudio es débil: la reserva difirió de Camaronal y de Ostional, pero no del Área Marina Cabo Blanco, y el área explicó una proporción muy pequeña de la variación. En el ancho, la reserva sí difirió de las otras tres áreas, aunque el tamaño del efecto también fue pequeño. Además, la categoría de manejo es solo una aproximación a la presión de extracción, que no se midió (Nielsen et al., 2025), y las diferencias de talla atribuidas a la extracción pueden confundirse con diferencias ambientales entre sitios, como la exposición al oleaje (Doyle et al., 2022). Por otra parte, el área sin extracción estuvo representada únicamente por la Reserva Cabo Blanco, por lo que esta comparación no permite separar el efecto de la protección de las características particulares de esa área.

Por lo tanto, las hipótesis planteadas recibieron un respaldo parcial: si bien la Reserva Cabo Blanco presentó individuos ligeramente mayores, este efecto fue débil y quedó superado por la marcada variación entre sitios y fechas de muestreo; además, la talla de N. scabricosta no siguió un patrón estacional consistente ni se asoció con la precipitación o la temperatura regional. Este patrón coincide con lo observado en el Golfo de California, donde la talla de la especie difirió entre islas cercanas según las condiciones de cada sitio (Gutiérrez-Rubio & Arzola-González, 2024), y en un plano distinto, el genómico, con lo reportado en Lottia gigantea, donde la variación climática dejó una huella selectiva mayor que la cosecha, aunque las áreas marinas protegidas sí mostraron una tendencia creciente en la abundancia de individuos grandes (Nielsen et al., 2025). En conjunto, estos resultados indican que, en el periodo evaluado, la identidad del sitio y el momento del muestreo influyen más en la estructura de tallas que la categoría de manejo del área.

La variación entre sitios de una misma área ayuda a explicar estos resultados. La distribución aplanada de la Reserva Cabo Blanco refleja dos grupos de individuos de diferente tamaño: la mayoría de los individuos pequeños provino de Punta Flor (54 de 68), mientras que San Miguel y Playa Cabo Blanco concentraron los individuos grandes. Punta Flor presentó tallas pequeñas tanto en 2022 como en 2023, lo que sugiere una característica estable del sitio y no un evento puntual. En el Golfo de California, el ancho medio de N. scabricosta varió de 12.2 a 26.7 mm entre islas cercanas, y los autores lo relacionaron con la pendiente, la disponibilidad de grietas y macroalgas, el oleaje y la temperatura (Gutiérrez-Rubio & Arzola-González, 2024). La exposición al oleaje también se ha asociado con la talla de otros gasterópodos (Rodríguez-Grimón et al., 2024) y con la distribución de N. scabricosta (Azofeifa-Solano et al., 2018). Como en este estudio no se midieron estas condiciones, no es posible separar su efecto del efecto del área.

De la misma manera, la composición de sitios explica parte de la variación entre años. En la Reserva Cabo Blanco, el muestreo de 2023 incluyó únicamente Punta Flor, mientras que en 2025 este sitio no se muestreó, de modo que el contraste entre esos años mezcla el efecto del año con el del sitio. La variación restante podría relacionarse con que N. scabricosta se reproduce solo durante la época lluviosa, sin sincronía con la marea ni con la fase lunar (Collin et al., 2017), y con que el tiempo de eclosión de sus cápsulas varía entre pozas de marea (Collin et al., 2016), lo que podría generar diferencias en el momento del reclutamiento entre sitios y años. En N. japonica, el reclutamiento modificó la densidad poblacional sin cambiar la talla media de los adultos (Takada, 2008), mientras que en Lottia gigantea una anomalía cálida aumentó el reclutamiento y la proporción de individuos pequeños (Nielsen et al., 2025). Estos procesos podrían contribuir a la variación interanual observada, aunque no pueden confirmarse con los datos disponibles.

La ausencia de un patrón estacional consistente apunta en la misma dirección. En la Reserva y el Área Marina Cabo Blanco en 2022, las dos épocas no incluyeron los mismos sitios: en ambos casos, la época seca incluyó un sitio que no se muestreó en la lluviosa. En cambio, en Ostional en 2025, el único caso en que se muestrearon los mismos tres sitios en ambas épocas, no se detectó diferencia. Esto sugiere que las diferencias entre épocas reflejan principalmente qué sitios se muestrearon y no un efecto de la estacionalidad sobre la talla.

El ANOVA de dos vías, aunque no estaba incluido en el plan inicial, agrega un matiz a esta interpretación. La interacción entre el área y el año indica que el largo cambió en direcciones opuestas entre 2025 y 2026 según el área. Sin embargo, en esos años el año coincidió con la época: en la Reserva y el Área Marina Cabo Blanco, los muestreos de 2025 se hicieron en la época lluviosa y los de 2026 en la seca, mientras que en Camaronal ocurrió lo contrario, y en las tres áreas el largo fue mayor en el año muestreado en la época lluviosa. Por ello, la interacción podría reflejar en parte un efecto de la época que no se observó de forma consistente en 2022, y este diseño no permite separar ambos factores.

La falta de asociación con el clima puede deberse tanto a las limitaciones de los datos como a la biología de la especie. Los 18 muestreos compartieron solo 11 valores distintos de temperatura, porque los sitios muestreados en una misma fecha recibieron el mismo valor, lo que reduce la capacidad de las pruebas para detectar una asociación. Además, la temperatura del aire registrada a 25–31 km no representa la que experimenta el caracol. En costas rocosas tropicales, la roca alcanza más de 50 °C y la temperatura corporal de especies de Nerita puede acercarse a los 50 °C (Marshall et al., 2023), muy por encima del intervalo de 25.9 a 28.3 °C registrado por las estaciones. Las especies de Nerita también seleccionan microhábitats entre 7 y 14 °C más frescos que su límite térmico (Hui et al., 2022), y su tolerancia térmica responde a la temperatura local de cada sitio (Hui et al., 2025). Por ello, el clima regional difícilmente captura las condiciones térmicas que influyen en el crecimiento.

Según el sistema de semáforo del PRONAMEC, la mayoría de los muestreos se ubicó en la categoría amarilla (regular) y ninguno en la roja (malo). De acuerdo con el protocolo, la categoría regular sugiere un grado medio de extracción y la roja, un exceso de extracción (SINAC, 2016), por lo que en ninguna de las áreas se detectaron señales de sobreexplotación severa. Sin embargo, la categoría regular también fue frecuente en la Reserva Cabo Blanco (tres de ocho muestreos), donde la extracción está prohibida, y abarcó todos los muestreos de Camaronal, lo que sugiere que el intervalo de 1.5 a 2.5 cm refleja en buena medida la talla habitual de la especie en estos sitios y no necesariamente extracción. Esto es coherente con que los umbrales del protocolo se basan en estudios de hace casi cuatro décadas (Ortega, 1986, 1987, como se citó en SINAC, 2016) y no se han calibrado para cada sitio. La asociación entre la categoría y el área (prueba exacta de Fisher, p = 0.015) se debió principalmente a que la Reserva Cabo Blanco concentró tres de los cuatro muestreos en la categoría verde (bueno), lo cual es consistente con la interpretación del protocolo, según la cual esta categoría indica ausencia de extracción (SINAC, 2016). Por otra parte, los muestreos sin categoría, dos de Punta Flor (Reserva Cabo Blanco) y dos de Mojón 69 (Ostional), se deben a que sus promedios se ubicaron entre 1.0 y 1.5 cm, intervalo para el cual el protocolo no define una categoría; estas observaciones no necesariamente representan un deterioro de las poblaciones, sino una limitación del protocolo para clasificar muestreos con predominio de individuos pequeños. Además, al evaluar solo el promedio, el indicador no distingue una población homogénea de individuos medianos de otra con una mezcla de tallas, como la de la Reserva Cabo Blanco. Definir este intervalo, calibrar los umbrales con datos locales y complementar el promedio con la distribución de frecuencias de tallas mejoraría la utilidad del indicador para el manejo costero.

Los individuos menores de 1 cm representaron solo el 3.5 % del total. Esta escasez probablemente refleja un sesgo de muestreo más que una ausencia real de juveniles: N. scabricosta se refugia en grietas y hendiduras para evitar la desecación y la depredación, y forma agregaciones con fines de protección (Gutiérrez-Rubio & Arzola-González, 2024). Si durante el muestreo no se revisan las grietas ni el interior de las agregaciones, los juveniles refugiados quedan subrepresentados en las mediciones.

Estos resultados deben interpretarse considerando varias limitaciones. Las pruebas entre áreas, años y épocas, y también el ANOVA de dos vías, trataron a cada individuo como una observación independiente, aunque los individuos de un mismo muestreo comparten las condiciones del sitio; esta pseudorreplicación tiende a subestimar los valores de p (Hurlbert, 1984), por lo que la interacción entre área y año debe considerarse exploratoria. En un análisis complementario, con el largo medio de cada muestreo como unidad (n = 34 muestreos), la diferencia entre áreas no fue significativa (p = 0.740). Además, el resultado entre áreas dependió de los sitios incluidos: al excluir los registros dudosos, que correspondían principalmente a Punta Flor en 2023, un sitio con ejemplares pequeños, la diferencia de la Reserva Cabo Blanco aumentó. Por último, los datos fueron recolectados por terceros y no se contó con las hojas de campo, lo que impidió verificar los registros dudosos, y el diseño fue desbalanceado, porque no todas las áreas se muestrearon en todos los años ni en ambas épocas; en particular, la época lluviosa de 2026 estaba en curso mientras se realizó este estudio, por lo que de ese periodo solo se contó con los muestreos de Camaronal de mayo de 2026.

Para futuras investigaciones sobre N. scabricosta en el litoral rocoso, se recomienda mantener los mismos sitios en todas las campañas y en ambas épocas, ya que la elección del diseño de monitoreo condiciona los parámetros que se pueden estimar (Wright et al., 2022). También se recomienda registrar cuantitativamente la presión de extracción, la afluencia turística y la exposición al oleaje; muestrear activamente dentro de grietas y agregaciones para estimar mejor la abundancia de juveniles menores de 1 cm; instalar, en la medida de lo posible, sensores de temperatura sobre las rocas intermareales para capturar el microclima que experimenta la especie, y analizar los datos con modelos lineales mixtos que incluyan el sitio de muestreo como efecto aleatorio.

En conclusión, la estructura de tallas de N. scabricosta en las áreas protegidas del Pacífico norte de Costa Rica depende principalmente de las características de cada sitio y del momento del muestreo, más que de la época o de las condiciones climáticas regionales. La Reserva Cabo Blanco presentó tallas mayores, pero esa diferencia es pequeña y no se sostiene al usar el muestreo como unidad, por lo que no es posible atribuirla con seguridad a la prohibición de extracción. Estos resultados destacan la necesidad de diseños de monitoreo balanceados y de calibrar los intervalos del protocolo PRONAMEC con datos locales.

Bibliografía

Alvarado, J. J., Cortés, J., Esquivel, M. F., & Salas, E. (2012). Costa Rica’s marine protected areas: Status and perspectives. Revista de Biología Tropical, 60(1), 129–142. https://doi.org/10.15517/rbt.v60i1.2753

Arias, A., Cinner, J. E., Jones, R. E., & Pressey, R. L. (2015). Levels and drivers of fishers’ compliance with marine protected areas. Ecology and Society, 20(4), 19. https://doi.org/10.5751/ES-07999-200419

Azofeifa-Solano, J. C., Salas-Moya, C., Cubero-Campos, Y., & Sibaja-Cordero, J. A. (2018). Influencia de la vigilancia en una zona rocosa dentro de un área marina protegida, Pacífico Central de Costa Rica. Revista de Biología Tropical, 66(3), 984–995.https://doi.org/10.15517/rbt.v66i3.30487

Baliwe, N. G., Pfaff, M. C., & Branch, G. M. (2022). Assessing the effects of no-take zones in a marine protected area spanning two ecoregions and rock substrate types. Frontiers in Marine Science, 9, 893260.https://doi.org/10.3389/fmars.2022.893260

Chacón-Monge, J. L., Cubero-Campos, Y., & Sibaja-Cordero, J. A. (2019). Efecto de un área marina protegida sobre la talla de la concha de dos moluscos del intermareal en el Pacífico Central de Costa Rica. Intropica, 14(2), 138–147. https://doi.org/10.21676/23897864.3365

Collin, R., Kerr, K., Contolini, G., & Ochoa, I. (2017). Reproductive cycles in tropical intertidal gastropods are timed around tidal amplitude cycles. Ecology and Evolution, 7(15), 5977–5991. https://doi.org/10.1002/ece3.3166

Collin, R., Roof, K. E., & Spangler, A. (2016). Hatching plasticity in the tropical gastropod Nerita scabricosta. Invertebrate Biology, 135(2), 87–96. https://doi.org/10.1111/ivb.12119

Daniel, W. W., & Cross, C. L. (2013). Biostatistics: A foundation for analysis in the health sciences (10.ª ed.). John Wiley & Sons.

Decreto Ejecutivo N.º 23640-MIRENEM. (1994, 28 de septiembre). Declara Refugio Nacional de Vida Silvestre Camaronal. La Gaceta N.º 184. San José, Costa Rica.

Decreto Ejecutivo N.º 40442-MINAE. (2017, 9 de agosto). Creación del Área Marina de Manejo Cabo Blanco. La Gaceta N.º 151. San José, Costa Rica.

Decreto Ley N.º 3763. (1966, 19 de octubre). Creación de la Reserva Natural Absoluta Cabo Blanco. Asamblea Legislativa de la República de Costa Rica. San José, Costa Rica.

Doyle, D., Frias, J., Gammell, M. P., Lynch, M., & Nash, R. (2022). Assessing the morphological impacts of long-term harvesting in intertidal gastropods using historical data and morphometric tools. Journal of Molluscan Studies, 88(3), eyac019. https://doi.org/10.1093/mollus/eyac019

Fallas-Flores, S., Abarca-Guerrero, L., & Ulate-Naranjo, K. (2024). Propuesta de mejoras para la implementación del monitoreo de residuos sólidos mediante el PRONAMEC en las Áreas Marinas Protegidas de Costa Rica. Tecnología en Marcha, 37(4), 31–47. https://doi.org/10.18845/tm.v37i4.6882

Faria, J., Vale, M., Ribeiro, P., Hawkins, S. J., & Martins, G. M. (2024). The ecological value of fully enforced, no-entry, marine protected areas: A case study of harvested limpets. Aquatic Conservation: Marine and Freshwater Ecosystems, 34(2), e4086. https://doi.org/10.1002/aqc.4086

Fox, J., & Weisberg, S. (2019). An R companion to applied regression (3.ª ed.). Sage. https://www.john-fox.ca/Companion/

Frank, P. W. (1969). Growth rates and longevity of some gastropod mollusks on the coral reef at Heron Island. Oecologia, 2, 232–250. https://doi.org/10.1007/BF00379161

Gutiérrez-Rubio, Y., & Arzola-González, J. F. (2024). Densidad y estructura de tallas de Nerita scabricosta (Mollusca: Neritidae) en la zona intermareal de las islas del noroeste del Golfo de California, México. Revista Científica FCV-LUZ, 34, e34377. https://doi.org/10.52973/rcfcv-e34377

Hui, T. Y., Crickenberger, S., Chan, S. H. M., Todd, P. A., & Williams, G. A. (2025). How does local temperature shape thermal tolerance? A test using congeneric snails on tropical rocky shores. The Biological Bulletin, 249(1), 13–29. https://doi.org/10.1086/741537

Hui, T. Y., Crickenberger, S., Lau, J. W. T., & Williams, G. A. (2022). Why are ‘suboptimal’ temperatures preferred in a tropical intertidal ectotherm? Journal of Animal Ecology, 91(7), 1400–1415. https://doi.org/10.1111/1365-2656.13690

Hurlbert, S. H. (1984). Pseudoreplication and the design of ecological field experiments. Ecological Monographs, 54(2), 187–211. https://doi.org/10.2307/1942661

Instituto Meteorológico Nacional. (2026). Oficio IMN-DIM-CM-192-2026: Remisión de información meteorológica de las estaciones San José-Pinilla, Paquera, Santa Rita y ASADA Artola, periodo 2022–2025. Departamento de Información Meteorológica, San José, Costa Rica.

Komsta, L., & Novomestky, F. (2022). moments: Moments, cumulants, skewness, kurtosis and related tests (R package version 0.14.1). https://doi.org/10.32614/CRAN.package.moments

Kunze, C., Wölfelschneider, M., & Rölfer, L. (2021). Multiple driver impacts on rocky intertidal systems: The need for an integrated approach. Frontiers in Marine Science, 8, 667168. https://doi.org/10.3389/fmars.2021.667168

Leeuwis, R. H., & Gamperl, A. K. (2022). Adaptations and plastic phenotypic responses of marine animals to the environmental challenges of the high intertidal zone. En S. J. Hawkins et al. (Eds.), Oceanography and Marine Biology: An Annual Review (Vol. 60, pp. 625–680). CRC Press. https://doi.org/10.1201/9781003288602-13

Ley N.º 9348. (2016, 2 de marzo). Ley reguladora del Refugio Nacional de Vida Silvestre Ostional. La Gaceta N.º 43. San José, Costa Rica.

Maia, R. C., & Troncoso, J. S. (2024). The effect of extreme climatic events on littorinid snails in two estuarine environments, temperate (NW Spain) and tropical (NE Brazil). Ocean and Coastal Research, 72(suppl. 1), e24017. https://doi.org/10.1590/2675-2824072.23060

Marshall, D. J., Mustapha, N., & Monaco, C. J. (2023). Conservation of thermal physiology in tropical intertidal snails following an evolutionary transition to a cooler ecosystem: Climate change implications. Conservation Physiology, 11(1), coad056. https://doi.org/10.1093/conphys/coad056

Nielsen, E. S., Walkes, S., Sones, J. L., Fenberg, P. B., Paz-García, D. A., Grosberg, R. K., Sanford, E., & Bay, R. A. (2025). The effects of latitudinal gradients, climatic anomalies, and size-selective harvesting on the adaptive potential of an intertidal gastropod. Evolutionary Applications, 18(9), e70159. https://doi.org/10.1111/eva.70159

Ooms, J., & Denney, B. (2026). writexl: Export data frames to Excel 'xlsx' format (R package version 2.0.1). https://doi.org/10.32614/CRAN.package.writexl

R Core Team. (2026). R: A language and environment for statistical computing. R Foundation for Statistical Computing. https://www.R-project.org/

Rodríguez-Grimón, R. O., Campos, N. H., & Castro, Í. B. (2024). Efecto de la exposición al oleaje en el tamaño y el índice de condición de los neogastrópodos en Santa Marta, Colombia. Acta Biológica Colombiana, 29(1), 24–31. https://doi.org/10.15446/abc.v29n1.89120

Sistema Nacional de Áreas de Conservación. (2016). Protocolo PRONAMEC: Protocolo para el monitoreo ecológico de las playas rocosas. Proyecto Consolidación de las Áreas Marinas Protegidas, Programa de las Naciones Unidas para el Desarrollo / Fondo para el Medio Ambiente Mundial. San José, Costa Rica.

Solís Rivera, V., & Fonseca Borrás, M. (2018). El área marina de Cabo Blanco y su relación con las actividades de pesca de pequeña escala: Plan general de manejo. CoopeSoliDar R.L.

Takada, Y. (2008). Contrasting characteristics in increasing and decreasing phases of the Nerita japonica (Mollusca: Gastropoda) population on a boulder shore. Population Ecology, 50, 391–403. https://doi.org/10.1007/s10144-008-0091-7

Wickham, H., & Bryan, J. (2026). readxl: Read Excel files (R package version 1.5.0). https://doi.org/10.32614/CRAN.package.readxl

Wright, A. D., Grant, E. H. C., & Zipkin, E. F. (2022). A comparison of monitoring designs to assess wildlife community parameters across spatial scales. Ecological Applications, 32(6), e2621. https://doi.org/10.1002/eap.2621

Zarzyczny, K. M., Paz-García, D. A., Williams, S. T., Rius, M., MacLean, M., & Fenberg, P. B. (2026). Tropicalisation of rocky shore gastropod communities along the Baja California Peninsula. Hydrobiologia. https://doi.org/10.1007/s10750-026-06429-z

Anexos

Fig. A1. Relación entre el largo medio de la concha de Nerita scabricosta por muestreo y (A) la temperatura media y (B) la precipitación acumulada de los 90 días previos, Área de Conservación Tempisque, Costa Rica, 2022–2025 (n = 18 muestreos). Cada punto representa un muestreo. RNA: Reserva Natural Absoluta Cabo Blanco; AMM: Área Marina de Manejo Cabo Blanco.