Caracterización de poblaciones de miriápodos en zona costera (San Fernando de Sámara) y zona de bosque nuboso (Poás de Alajuela), Costa Rica

Nahomy Tamara Jirón-Rodríguez
Jeferson Cordero-Montiel
Juan Diego Rojas-Rodríguez
Melissa Torres-Sánchez
Larissa Guillén-Bonilla

ABSTRACT

Introduction: Myriapods include groups with important ecological niches, being that they are dead matter and small invertebrate controllers. In Costa Rica, the ecosystemic variability present induces a big diversity in animal groups, including myriapods, inside these ecosystems, dry forests in Sámara and humid tropical/cloudy forests like the ones in Poás can be found, and the differences between them can generate changes in the organisms that inhabit them, having to adapt to the conditions provided by each site.
Objective: Evaluate the ecosystemic effects that the forests of Sámara of Guanacaste and Poás of Alajuela have on the characteristics of myriapod populations.
Methods: For specimen detection of each group, a method of registry by 1m² aleatory quadrants was used, doing 30 repetitions each day for 6 days on each site, collecting data of soil types, taxonomic groups, phenotype, dimorphism in diplopods, length, width and weight, by which the myriapod families found were determined using field guides, furthermore, using these data for the corresponding statistical analysis with 95% confidence level, (α = 0.05).
Results: The initial sampling presented a bigger density of individuals in San Fernando of Sámara, in comparison, a higher variety of families were collected in San Rafael of Poás, in general, the data of 316 individuals, 3 soil types and 9 families were collected.
Using a fisher test, evidence of significant association between the amount of individuals by quadrant and the sampling site was found, acquiring as a result a bigger number of individuals in San Fernando of Sámara; utilizing one way chi-squared it was demonstrated that there was a significant difference in proportions of males and females in diplopod populations (X²= 13.337, p value= 0,00026). The two correlation tests done of the weight making use of the length and width variables were used to develop a multiple lineal regression model, with which it was demonstrated that both variables are significant predictors of body weight. Following with the previous point, the data was transformed into logarithmic form as a way to better adjust the model, providing evidence of an alometric relationship between body size and volume of individuals. Additionally, the model was able to perform a prediction for a myriapod of 13cm of length and 0.5cm of width. The Kruskal- Wallis non parametric variance analysis applied to the length of individuals showed significant differences for each taxonomic class of the same site, as well as for the different study sites. The second Kruskal-Wallis non parametric variance analysis for the weight of the individuals indicated that there was a significant differences in the weights of diplopods from San Fernando and chilopods from San Rafael of Poás. Alongside, it was evidenced that there were not significant differences between the weights of diplopods of both sites. The diversity index tests (Shannon-Weiner & Simpson) demonstrated that there was more uniformity in San Fernando, being dominated by the paradoxosomatidae family; meanwhile in Poás, the indexes have a higher variability tendency of registered families.
Conclusion: The study managed to demonstrate the differences that the morphological characteristics, population density and taxonomic diversity can have in myriapoda (diplopoda, chilopoda) populations depending of the ecosystem on which they live, however, it would be beneficial to the study to expand the methodology as a way to deepen the knowledge of how other characteristics present on each site can induce even more changes.

Keywords: Myriapoda, Chilopoda, Diplopoda, Sámara, Poás, Costa Rica, Diversity

RESUMEN

Introducción: Los miriápodos incluyen grupos con importantes nichos ecológicos, al ser controladores de materia muerta y poblaciones de invertebrados pequeños. En Costa Rica, la variabilidad de ecosistemas presentes propicia una gran diversidad en grupos animales, incluídos los miriápodos, dentro de estos ecosistemas se encuentran los bosques secos de Sámara y los bosques húmedos tropicales/nubosos de Poás, cuyas diferencias entre ellos pueden llegar a generar también cambios en los organismos que en ellos habitan, al adaptarse a las condiciones proporcionadas por cada sitio. Objetivo: Evaluar los efectos ecosistémicos que tienen los bosques de Sámara de Guanacaste y Poás de Alajuela en las características poblacionales de miriápodos.
Metodología: Para la detección de individuos de cada grupo se usó un método de registro por cuadrantes aleatorios de 1m², haciendo 30 repeticiones por 6 días en cada sitio, recolectando datos de tipo de suelo, grupo taxonómico, fenotipo, dimorfismo de diplópodos, largo, ancho y peso, gracias a esto se determinaron las familias de miriápodos presentes utilizando guías de campo, además de proporcionar la información para el análisis estadístico correspondiente con niveles de confianza al 95% (α = 0.05).
Resultados: El muestreo inicial presentó una mayor densidad de individuos en San Fernando de Sámara, mientras que se recolectaron una mayor variedad de familias en San Rafael de Poás, en general, se obtuvieron los datos de 316 individuos, 3 tipos de suelo y 9 familias. Por medio de la prueba de Fisher se logró evidenciar una asociación significativa entre la cantidad de individuos por cuadrante con el sitio de muestreo, consiguiendo como resultado un mayor número de individuos en San Fernando de Sámara, mediante la prueba de chi cuadrado de una sola vía se demostró una diferencia significativa entre la proporción de machos y hembras en la población de diplópodos (X²= 13.337 y un valor de p= 0,00026). Las dos pruebas de correlación con respecto al peso usando las variables de largo y ancho se usaron para desarrollar un modelo de regresión lineal múltiple la cual demostró que ambas variables son predictores significativos del peso corporal. Seguidamente se realizó una transformación logarítmica que mejoró el ajuste del modelo, evidenciando una relación alométrica entre el tamaño corporal y el volumen de los individuos. Además, el modelo permitió predecir el peso para un miriápodo de 13 cm de largo y 0.5 cm de ancho.
El análisis de varianza no paramétrico Kruskal-Wallis para el largo de los individuos presentó diferencias significativas para cada clase taxonómica del mismo sitio, al igual que para los diferentes sitios de estudio. El segundo análisis de varianza no paramétrico Kruskal-Wallis para el peso respecto a los individuos indicó diferencias significativas entre los pesos tomados para diplópodos de San Fernando y los quilópodos de San Rafael de Poás. A su vez se evidenció que no existen diferencias significativas entre los pesos de diplópodos en ambos sitios. Las pruebas con índices de diversidad (Shannon-Weiner y Simpson) demostraron más uniformidad en San Fernando, siendo dominado por la familia Paradoxosomatidae, mientras que en Poás los índices tienden a una mayor variabilidad en las familias registradas.
Conclusiones: El estudio logró demostrar las diferencias que pueden presentar las características morfológicas, densidad poblacional y diversidad taxonómica en las poblaciones de miriápodos (Diplopoda, Chilopoda) dependiendo del ecosistema en que viven, pero aún así se beneficiaría de una expansión metodológica para profundizar los posibles efectos que otras características presentes en cada sitio puedan tener.

Palabras clave: Miriápodos, Chilopoda, Diplopoda, Sámara, Poás, Costa Rica, Diversidad

INTRODUCCIÓN

La importancia de los miriápodos recae en el hecho que tienen un rol como bioindicadores terrestres. Estos organismos reciclan la materia orgánica presente en los suelos. En el caso de los mil pies estos son detritívoros por lo que facilitan la incorporación de materia orgánica al suelo, y en comparación a otros organismos como son los hongos y bacterias su proceso es más rápido. Al incorporar material orgánico permiten el aumento de los niveles de nitrógeno del suelo lo que activa a las bacterias fijadoras de nitrógeno en el suelo y plantas, además participan en el proceso de humidificación del suelo (Garcia, J. L. , 2017).

En Costa Rica existe una gran biodiversidad, en parte esto se debe a su posición geográfica, su clima neotropical y su variedad de hábitats. Entre los ecosistemas presentes en Costa Rica están las selvas tropicales, los bosques caducifolios, los bosques nubosos la costa del Atlántico y el Pacífico y los manglares. Con relación a los miriápodos estos proliferan en ecosistemas tropicales debido a la cantidad de recursos que pueden encontrar. Sin embargo, estos también se pueden encontrar en zonas templadas, solo que en un menor número de individuos (Garcia, J. L. , 2017).

El bosque seco tropical de Guanacaste, ubicado principalmente en el Área de Conservación Guanacaste (ACG), es el remanente más extenso de este ecosistema en Mesoamérica y se caracteriza por una marcada estacionalidad, con una estación seca de diciembre a mayo y una estación lluviosa durante el resto del año (Janzen, 1988). Como adaptación a la escasez de agua, muchas especies arbóreas son caducifolias y pierden sus hojas durante la época seca, generando una abundante capa de hojarasca que constituye un hábitat importante para diversos artrópodos, como hormigas, arañas y coleópteros (Tamarindo Estuary Tours Mi Rey, s. f.). A pesar de las condiciones de sequía y las altas temperaturas, este ecosistema mantiene una considerable riqueza biológica (Janzen, 1988).

Por su parte, el bosque nuboso del Volcán Poás se desarrolla entre 1.200 y 2.708 m s. n. m. y presenta una elevada precipitación, humedad constante y abundancia de epífitas, características que favorecen una gran diversidad de flora, fauna y artrópodos mediante la creación de numerosos microhábitats y una compleja estratificación vertical (Fernández-Landa et al., 2017). Además, la producción continua de materia orgánica sostiene comunidades estables de organismos descomponedores y detritívoros. La relevancia de estudiar estos ecosistemas es especialmente alta considerando que se estima que existen alrededor de 7 millones de especies de artrópodos terrestres y 5,5 millones de especies de insectos en el mundo, muchas de las cuales aún permanecen sin describir científicamente (Stork, 2018).

La presente investigación tiene como objetivo comparar la biodiversidad de los miriápodos en dos zonas de Costa Rica que presentan ecosistemas diferentes, el bosque seco Tropical y bosque húmedo con la finalidad de estudiar la salud de los ecosistemas y cómo estos afectan las condiciones del desarrollo de los miriápodos, el peso, la longitud y la biodiversidad que hay en cada uno.

MATERIALES Y MÉTODOS

Se decidió realizar esta investigación para comparar la biodiversidad de miriápodos en dos zonas de Costa Rica, como una forma de establecer la salud y diversidad ecosistémica.

Muestreo de especies
En esta investigación se realizaron muestreos en San Fernando de Sámara, Guanacaste con una altura de 16 m.s.n.m en promedio (Lizano & Lizano, 2023) y en San Rafael de Poás, Alajuela cuya altura ronda los 1150 m.s.n.m (Municipalidad de Poás, 2026), al ser lugares cuyos ecosistemas son inversos en muchos aspectos, siendo San Fernando de Sámara caracterizado por su bosque seco (Fondo Nacional de Financiamiento Forestal, 2024) mientras que los bosques del cantón de Poás suelen ser de nubosos a tropicales húmedos (Sistema Nacional de Áreas de Conservación, 2025), es por esto que nace la duda de cómo influyen estas diferencias en el desarrollo y diversidad de los miriápodos.

Una vez en cada sitio de muestreo, se utilizó un cuadrante de 1m² construido de tubo PVC el cual se arrojó alrededor de la zona a manera aleatoria, contabilizando solo los individuos dentro del cuadrante; A estos mismos se les registraba el tipo de suelo donde se encontraron, hora de recolección las muestras, el ancho y largo utilizando una cinta métrica, peso mediante una balanza digital portátil y la cantidad de individuos por cuadrante, a estos se les clasificaba en diplópodos o quilópodos y finalmente se registraba su fenotipo para la posterior determinación de la familia de la especie utilizando coincidencias fenotípicas y de distribución geográfica comparando los registros fotográficos del estudio con los de registros públicos del Inaturalist y el sistema global de información sobre biodiversidad (GBIF) y la guía de campo “Spiders and Their Kin”; Adicionalmente, a los diplópodos únicamente, se les registró la presencia de dimorfismo sexual, esto porque a diferencia de los quilópodos, la clase diplópoda presenta un dimorfismo marcado al registrar la presencia o no de gonopodios. Para asegurar la veracidad de la presencia o ausencia de ejemplares en los cuadrantes, se removió el sedimento foliar, piedras y troncos pues estos animales se suelen encontrar debajo de estos. Se realizaron 30 repeticiones por día, durante 6 días en ambos sitios. Cabe destacar que los organismos fueron muestreados y devueltos al mismo sitio donde se encontraron al finalizar cada registro, no se registró muerte de individuos en el estudio.

Análisis Estadístico
Inicialmente para comprender cómo se distribuyen los datos obtenidos, se realizaron análisis de Estadística descriptiva para variables cuantitativas tales como la media, mediana, desviación estándar, varianza y máximos por sitio de muestreo. Todas estas se aplicaron para el largo dado en centímetros (cm), ancho en centímetros (cm) y peso en gramos (g). Seguido de esto, para evaluar el cumplimiento del supuesto de una distribución normal en las variables cuantitativas registradas en los muestreos, se aplicó la Prueba de Shapiro para cada grupo y zona de estudio previamente determinada con niveles de confianza al 95%. Como alternativa no paramétrica, se aplicó la Prueba de Wilcoxon para la evaluación de muestras pareadas, siendo utilizada para los datos biométricos en dos zonas diferentes (α = 0.05).
Para comparar el número de individuos en relación con las zonas de muestreo se realizó una prueba de chi cuadrado de dos vías (α = 0.05), con el fin de poder analizar la densidad poblacional en relación con los climas y la disponibilidad de alimentos en estos sitios.

Posteriormente se realizaron dos pruebas para analizar la correlación entre las variables obtenidas. A partir de las pruebas de correlación se analizaron las variables mediante la regresión lineal múltiple la cual se utilizó para predecir un dato. Se realizó una prueba de Chi cuadrado de una sola vía para evaluar la distribución de sexo de los Miriápodos.

Resultados

Muestreo:
Preliminarmente se vieron los 20 diferentes fenotipos constatados en la tabla 1; registrando datos de 316 individuos, esto de ambos sitios en conjunto, con 3 tipos de suelo, tierra seca con sedimento foliar, rocoso con vegetación y barrial húmedo de cafetal.

Fig. 1 Registro Completo de fenotipos encontrados en ambos sitios de estudio, el contenido muestra cómo se clasificaron los fenotipos vistos, de izquierda a derecha, PA: Patas Azules, NR: Negro Rojo, BN: Blanco Negro, NAR: Naranja, BLA: Blanco, CHOC: chocolate, LAD: Ladrillo, CPR: Café Patas Rayadas, PL: Patas Largas, REABEI: Real beige, REA: Real, REARAN: Real Rayas Naranja, REAGRI: Real Gris, REAXL: Real Extra Large, GPEQ: Gris Pequeño, GC: Gris Claro, G: Gris, NA: Negro Amarillo, GNA: Grandes Negro Amarillos, CA: Café Amarillo, fotografías propias.

Mediante las guías mencionadas en la metodología se pudo determinar que los fenotipos NA, GC, G, GPEQ son la misma especie, pasaron a ser todos NA, este también fue el caso para REABEI, REAGRI y REA, se les generalizó como REA, BN y CPR, registrados como BN y finalmente GNA y CA registrados como GNA, las diferencias fenotípicas atribuyeron a la etapa de vida del organismo así como el proceso de ecdisis, que en este subfilo tiene como consecuencia una pérdida de pigmentación del individuo. Unificadas estas especies, se analizó también factores como número de patas y segmentos. Gracias a esto se llegó a la distribución de familias taxonómicas de la tabla 1.

Tabla 1.
Resultado de la división de familias por fenotipo

Fenotipo Familia
NA Paradoxosomatidae
GNA Chelodesmidae
REA Spirobolidae
REARAN Spirobolidae
REAXL Spirostreptidae
NR Scolopendridae
NAR Cryptopidae
BLA Geophilidae
BN Scolopendridae
PL Scutigeridae
PA Scolopendridae
CHOC Lithobiidae
LAD Lithobiidae
Una vez teniendo la división de la tabla 1, se hizo el recuento de individuos por familia y por sitio por medio de las tablas 2 y 3 para su posterior procesamiento estadístico.

Tabla 2.
Cantidad de individuos de cada familia vista en San Fernando de Sámara

Familia Cantidad de Individuos
Paradoxosomatidae 226
Chelodesmidae 5
Spirobolidae 10

Tabla 3.
Cantidad de individuos de cada familia vista en San Rafael de Poás

Familia Cantidad de Individuos
Paradoxosomatidae 5
Chelodesmidae 2
Spirobolidae 30
Spirostreptidae 1
Lithobiidae 3
Scolopendridae 7
Scutigeridae 3
Geophilidae 2
Cryptopidae 22

Por otro lado, el total de diplópodos vistos en el estudio fue de 279 individuos, mientras que de chilopodos 37, cuya totalidad fue registrada en Poás, además de esto, 110 de los diplópodos fueron machos (presentaban gonopódios) y 169 fueron hembras.

Análisis estadístico: Seguidamente, se realizaron los análisis estadísticos, cabe destacar que todos se trabajaron con un nivel de confianza del 95%. Primeramente, como parte de la estadística descriptiva para variables cuantitativas tales como la media, mediana, desviación estándar, varianza y máximos. Se obtuvo que el largo presenta una desviación estándar mayor (2.26 cm) que la presente en ancho (0.17 cm) y peso (0.52 g) en la localidad de Poás, en San Fernando la desviación estándar del largo también es la mayor (1.82 cm), pero en menor medida en comparación con la de Poás. Este patrón se repite de igual manera con el máximo, Poás tiene una medición máxima mayor (12 cm) que en San Fernando (9 cm)

Se realizaron gráficos de densidad para comprobar la distribución de los datos, correspondientes a las figuras 2, 3 y 4. Se evidencia un sesgo de la curva positivo en las tres figuras, siendo también un indicador de la no existencia de normalidad para los datos seleccionados.


Fig. 2 Distribución por funcion de densidad del peso tomado de los Miriápodos en Samara y Poás
Fig. 3 Distribución por funcion de densidad del largo tomado de los Miriápodos en Samara y Poás


Fig.4 Distribución por funcion de densidad del ancho tomado de los Miriápodos en Samara y Poás 

shapiro.test(base$`Peso (g)`)
## 
##  Shapiro-Wilk normality test
## 
## data:  base$`Peso (g)`
## W = 0.58649, p-value < 2.2e-16
shapiro.test(base$`Largo (cm)`)
## 
##  Shapiro-Wilk normality test
## 
## data:  base$`Largo (cm)`
## W = 0.84509, p-value < 2.2e-16
shapiro.test(base$`Ancho (cm)`)
## 
##  Shapiro-Wilk normality test
## 
## data:  base$`Ancho (cm)`
## W = 0.81402, p-value < 2.2e-16

Mediante los datos obtenidos respecto al peso, largo y ancho de los individuos, se realizó la prueba de Shapiro para comprobar normalidad (ver Tabla 4).

Tabla 4.
Resultados de la prueba de normalidad

Variables Valor de W Valor de p Normalidad
Peso (g) 0.58649 2.2e-16 No existe normalidad
Largo (cm) 0.84509 2.2e-16 No existe normalidad
Ancho (cm) 0.81402 2.2e-16 No existe normalidad
En base a los resultados de la prueba de normalidad, se optó por aplicar una prueba de muestras pareadas basadas en Wilcoxon, evaluando las posibles diferencias de las muestras tomadas en los dos sitios, cuyos resultados se encuentran en la Tabla 5.

Tabla 5.
Resultados obtenidos al aplicar la prueba par de Wilcoxon

library(knitr)
individuos  <- subset(base, `Peso (g)` > 0)
individuos3 <- subset(base, `Largo (cm)` > 0)
individuos4 <- subset(base, `Ancho (cm)` > 0)
w_peso  <- wilcox.test(`Peso (g)` ~ Sitio, data = individuos)
w_largo <- wilcox.test(`Largo (cm)` ~ Sitio, data = individuos3)
w_ancho <- wilcox.test(`Ancho (cm)` ~ Sitio, data = individuos4)
Variables Valor de W Valor de p
Peso (g) 12246 3.411e-06
Largo (cm) 8971 0.9237
Ancho (cm) 13759 3.606e-13

Posteriormente, se realizaron gráficos de cajas para evaluar las distribuciones de las variables entre los sitios de muestreo, como un apoyo visual a la prueba de Wilcoxon (ver Figuras 5, 6 y 7). Tal como se observa en la Tabla 5, las distribuciones del peso (g) y el ancho (g) presentan diferencias significativas, caso contrario a la distribución designada a los valores del largo (cm).

Fig.5 Distribución del peso,largo y ancho tomado para miriápodos según el sitio de muestreo

Distribución de los Miriápodos en los sitios estudiados: San Fernando, San Rafael de Poás utilizando chi cuadrado de dos vías. Se utilizó el criterio de fischer en lugar de la prueba de chi cuadrado para obtener resultados confiables. Se analizó la distribución de los organismos en relación con los sitios de muestreo y sus climas. 

## Número de registros con presencia de individuos:
## [1] 316
## Total de individuos registrados:
## [1] 790
## Distribución de individuos por sitio:
##   Sitio Individuo
## 1    SF       695
## 2   SPP        95
tabla_sitio_individuo <- table(
  Base_Sitios_Filtrada$Sitio,
  Base_Sitios_Filtrada$Individuo
)

tabla_sitio_individuo
##      
##        1  2  3  4  5  6  7  8  9 10 11 12 13 14
##   SF  88 55 33 21 13  9  7  5  4  2  1  1  1  1
##   SPP 57 16  2  0  0  0  0  0  0  0  0  0  0  0
## 
##  Fisher's Exact Test for Count Data with simulated p-value (based on
##  2000 replicates)
## 
## data:  tabla_sitio_individuo
## p-value = 0.0004998
## alternative hypothesis: two.sided
## [1] 316  13
## [1] 0

Se aplicó la prueba exacta de Fisher mediante la simulación de Monte Carlo (100.000 repeticiones) para evaluar la asociación entre el sitio de muestreo y la distribución de la abundancia de diplópodos por cuadrante, para evitar errores de conteo ya que SRP no obtuvo datos de diplópodos mientras que en SF si se obtuvieron . Se obtuvo un valor de significancia de p=0.00089, inferior al nivel de significancia establecido (=0.05). Por lo tanto, se rechaza la hipótesis nula de independencia, lo que indica que hay una diferencia altamente significativa entre la cantidad de individuos y el sitio de observación confirmando la hipótesis alternativa en la que se menciona que se encuentran más Diplópodos en las zonas cálidas. Además, se observó una mayor representación de cuadrantes con abundancias elevadas en el sitio SF en comparación con SRP.

##      
##        1  2  3  4  5  6  7  8  9 10 14
##   SF  33 22 12  8  4  2  2  1  2  1  1
##   SPP 41 14  2  0  0  0  0  0  0  0  0
set.seed(123)

resultado_fisher <- fisher.test(
  tabla_individuos_sitio,
  simulate.p.value = TRUE,
  B = 100000
)

resultado_fisher
## 
##  Fisher's Exact Test for Count Data with simulated p-value (based on
##  1e+05 replicates)
## 
## data:  tabla_individuos_sitio
## p-value = 0.00089
## alternative hypothesis: two.sided

Tabla 6.
Tabla de contingencia sitios y número de individuos por cuadrante

Sitio Presencia de un individuo Más de un individuo Total
SF 33 56 89
SRP 41 16 57
Total  74 72 146


Según la tabla 6 se aprecia que en SRP se observa que la población tiende a densidades críticamente bajas por unidad de muestreo, acumulando 123 cuadrantes examinados en total, con cero individuos en varios de los cuadrantes, sin altas densidades de individuos en la mayoría y un máximo absoluto de apenas 3 individuos por cuadrante. En cambio, San Fernando exhibió una distribución mucho más homogénea hacia densidades altas, alcanzando registros máximos de hasta 14 individuos en un solo cuadrante.

resultado_chi2 <- chisq.test(tabla_diplo)



print("Resultado_Chi_1V")
## [1] "Resultado_Chi_1V"
print(resultado_chi2)
## 
##  Chi-squared test for given probabilities
## 
## data:  tabla_diplo
## X-squared = 328.13, df = 4, p-value < 2.2e-16
print("Valores esperados")
## [1] "Valores esperados"
print(resultado_chi2$expected)
## Chilo Diplo    ND    No    Si 
## 126.4 126.4 126.4 126.4 126.4
print("Residuos de la prueba")
## [1] "Residuos de la prueba"
print(resultado_chi2$residuals)
##     Chilo     Diplo        ND        No        Si 
## -7.951773 13.573161 -7.951773  3.878046 -1.547661
tabla_diplo
## Chilo Diplo    ND    No    Si 
##    37   279    37   170   109
Base_Limpia$Dimorfismo <- factor(Base_Limpia$Dimorfismo,
                                     levels = c("No", "Si"),
                                     labels = c("Hembra", "Macho"))


tabla_diplo_final <- table(Base_Limpia$Dimorfismo)

print("--- TABLA CON NOMBRES CORREGIDOS ---")
## [1] "--- TABLA CON NOMBRES CORREGIDOS ---"
print(tabla_diplo_final)
## 
## Hembra  Macho 
##    170    109
print("--- PRUEBA DE CHI-CUADRADO FINAL ---")
## [1] "--- PRUEBA DE CHI-CUADRADO FINAL ---"
chisq.test(tabla_diplo_final)
## 
##  Chi-squared test for given probabilities
## 
## data:  tabla_diplo_final
## X-squared = 13.337, df = 1, p-value = 0.0002602

Prueba de una sola vía: Con la prueba de Fisher se busca demostrar la H0:La distribución de géneros en la muestra de diplópodos está equilibrada: 50% machos y 50% hembras mientras que la H1 indica que la distribución de sexo en los diplópodos no es proporcional machos ≠ hembras. Para analizar el sexo de los diplópodos se realizó una revisión de la presencia de los gonopodos en machos, los que no lograron ser identificados se marcaron con un “0” y al realizar la prueba de chi cuadrado fueron excluidos del análisis para evitar sesgos de error y conseguir una mayor fidelidad de los resultados obtenidos. Una vez con la base de datos preparada se realizó la prueba de chi cuadrado obteniendo los siguientes resultados: Chi-cuadrado = 13.337, gL = 1, valor p = 0.0002602 este valor de p rechaza la hipótesis nula aceptando la hipótesis alternativa donde el número de machos y hembras no es proporcional.

Se realizaron dos pruebas de correlación lineal para hallar la relación entre las variables (largo y ancho), con respecto al peso (lmtest)(carData). Antes de evaluar la relación entre las variables biométricas de los miriápodos, se verificó el supuesto de normalidad mediante la prueba de Shapiro-Wilk. En esta prueba, el estadístico W toma valores entre 0 y 1, donde valores cercanos a 1 indican una mayor aproximación a la normalidad. Para el peso se obtuvo W = 0.58649, para el largo W = 0.84509 y para el ancho W = 0.81402, mientras que en los tres casos el valor de significancia fue p < 2.2 × 10⁻¹⁶. Dado que los valores de p son mucho menores que 0.05, se rechaza la hipótesis nula de normalidad y se concluye que ninguna de las variables presenta una distribución normal. Debido a esta desviación de la normalidad, se empleó la correlación de Spearman (ρ) para evaluar relaciones monotónicas entre variables sin asumir normalidad.

shapiro_peso <- shapiro.test(base$`Peso (g)`)
shapiro_peso
## 
##  Shapiro-Wilk normality test
## 
## data:  base$`Peso (g)`
## W = 0.58649, p-value < 2.2e-16
shapiro_largo <- shapiro.test(base$`Largo (cm)`)
shapiro_largo
## 
##  Shapiro-Wilk normality test
## 
## data:  base$`Largo (cm)`
## W = 0.84509, p-value < 2.2e-16
shapiro_ancho <- shapiro.test(base$`Ancho (cm)`)
shapiro_ancho
## 
##  Shapiro-Wilk normality test
## 
## data:  base$`Ancho (cm)`
## W = 0.81402, p-value < 2.2e-16


Fig. 6 Correlación entre el peso y el largo


Fig. 7 Correlación entre el peso y el ancho 

cor_largo <- cor.test(base$`Peso (g)`, base$`Largo (cm)`, method = "spearman")
cor_largo
## 
##  Spearman's rank correlation rho
## 
## data:  base$`Peso (g)` and base$`Largo (cm)`
## S = 1092930, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.9562013
cor_ancho <- cor.test(base$`Peso (g)`, base$`Ancho (cm)`, method = "spearman")
cor_ancho
## 
##  Spearman's rank correlation rho
## 
## data:  base$`Peso (g)` and base$`Ancho (cm)`
## S = 931947, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.9626526

Los resultados mostraron una correlación positiva muy fuerte y altamente significativa entre el peso y el largo (ρ = 0.9562, p < 2.2 × 10⁻¹⁶) y entre el peso y el ancho (ρ = 0.9627, p < 2.2 × 10⁻¹⁶). El coeficiente de determinación aproximado (R² = ρ²) fue de 91.4 % para peso-largo y 92.7 % para peso-ancho. La tendencia a que aumente el peso si el largo o el ancho también aumentan se ve reflejado en las Figuras 6 y 7 respaldando estos resultados. Aunque existe cierta dispersión y algunos valores extremos en individuos de gran tamaño, la concentración general de los puntos sigue un patrón creciente constante, evidenciando que tanto el largo como el ancho constituyen predictores importantes del peso corporal en los miriápodos analizados.

A continuación, se realizó un análisis de regresión lineal múltiple con el objetivo de predecir el peso de los miriápodos a partir de las variables utilizadas anteriormente (largo y ancho). Se utilizaron los resultados y la interpretación de las correlaciones como fundamento para realizar esta prueba. Este análisis permitió evaluar la capacidad predictiva de las variables explicativas y determinar el cumplimiento de los supuestos estadísticos necesarios para la validez del modelo. Adicionalmente, se construyó un modelo de regresión múltiple logarítmica.

## 
## Call:
## lm(formula = `Peso (g)` ~ `Largo (cm)` + `Ancho (cm)`, data = base)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.0032 -0.1407  0.0414  0.1001  4.2702 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -0.10005    0.01679  -5.959 4.64e-09 ***
## `Largo (cm)`  0.10163    0.01210   8.398 4.20e-16 ***
## `Ancho (cm)`  0.78019    0.11768   6.630 8.33e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.266 on 528 degrees of freedom
## Multiple R-squared:  0.6569, Adjusted R-squared:  0.6556 
## F-statistic: 505.5 on 2 and 528 DF,  p-value: < 2.2e-16
library(car)
residuos_shapiro <- shapiro.test(residuals(modelo1))
residuos_shapiro
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(modelo1)
## W = 0.48527, p-value < 2.2e-16
ncvTest(modelo1)
## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 2127.735, Df = 1, p = < 2.22e-16
dwt(modelo1)
##  lag Autocorrelation D-W Statistic p-value
##    1      0.08754198      1.824622   0.032
##  Alternative hypothesis: rho != 0
vif(modelo1)
## `Largo (cm)` `Ancho (cm)` 
##     4.738001     4.738001
par(mfrow = c(2,2))
plot(modelo1)


Fig. 8 Diagnóstico del modelo de regresión en función al peso

Para evaluar la calidad del modelo se analizaron los gráficos diagnósticos y diversos estadísticos de validación. En el gráfico Residuals vs Fitted se observa una dispersión de los residuos alrededor de cero con algunos valores alejados de la tendencia general, lo que indica que el modelo explica gran parte de la variación del peso, aunque existen observaciones con errores de predicción elevados. El gráfico Normal Q-Q muestra una desviación evidente de los puntos respecto a la línea teórica, sugiriendo que los residuos no siguen una distribución normal. En el gráfico Scale-Location se aprecia un ligero incremento de la dispersión para determinados valores ajustados, aunque sin un patrón severo de heterocedasticidad. Finalmente, en Residuals vs Leverage destacan algunas observaciones influyentes, particularmente la observación 407, que ejerce un efecto importante sobre los coeficientes estimados.

Los resultados estadísticos del modelo mostraron un valor de F = 505.5, indicando que el conjunto de variables predictoras explica una proporción significativa de la variación observada en el peso. El error estándar residual (RSE = 0.226) indica que las predicciones del modelo se desvían, en promedio, aproximadamente ±0.226 g de los valores reales observados. La normalidad de los residuos fue evaluada mediante la prueba de Shapiro-Wilk, obteniéndose W = 0.485 y p < 2.2 × 10⁻¹⁶. El estadístico W mide el grado de ajuste a una distribución normal, donde valores cercanos a 1 indican normalidad; por lo tanto, el valor obtenido refleja una fuerte desviación de este supuesto. Asimismo, el valor de p inferior a 0.05 permite rechazar la hipótesis nula de normalidad. La prueba de heterocedasticidad arrojó un valor de χ² = 2127, mientras que el estadístico de Durbin-Watson (DW = 1.825) fue cercano al valor ideal de 2, indicando una independencia aceptable entre los residuos y ausencia de autocorrelación significativa. Por otra parte, el análisis de multicolinealidad produjo un VIF = 4.74, valor que indica una multicolinealidad moderada entre las variables explicativas, aunque insuficiente para comprometer seriamente el desempeño del modelo. En conjunto, estos resultados sugieren que el modelo posee una buena capacidad predictiva, aunque algunos supuestos fundamentales, especialmente la normalidad de los residuos, no se cumplen completamente.

datos_log <- subset(base, `Peso (g)` > 0 & `Largo (cm)` > 0 & `Ancho (cm)` > 0)
regresion_log <- lm(log(`Peso (g)`) ~ log(`Largo (cm)`) + log(`Ancho (cm)`), data = datos_log)
summary(datos_log)
##       Dia              Sitio       Cuadrante       Individuo        Dimorfismo 
##  Min.   :1.000   Length   :316   Min.   : 1.00   Min.   : 1.0   Length   :316  
##  1st Qu.:3.000   N.unique :  2   1st Qu.: 8.00   1st Qu.: 1.0   N.unique :  3  
##  Median :3.000   N.blank  :  0   Median :16.00   Median : 2.0   N.blank  :  0  
##  Mean   :3.522   Min.nchar:  2   Mean   :15.73   Mean   : 2.5   Min.nchar:  2  
##  3rd Qu.:5.000   Max.nchar:  3   3rd Qu.:24.00   3rd Qu.: 3.0   Max.nchar:  2  
##  Max.   :6.000                   Max.   :30.00   Max.   :14.0                  
##       Fenotipo     Largo (cm)       Ancho (cm)        Peso (g)    
##  Length   :316   Min.   : 0.500   Min.   :0.1000   Min.   :0.010  
##  N.unique : 13   1st Qu.: 3.000   1st Qu.:0.3000   1st Qu.:0.165  
##  N.blank  :  0   Median : 3.700   Median :0.4000   Median :0.470  
##  Min.nchar:  2   Mean   : 3.644   Mean   :0.3585   Mean   :0.482  
##  Max.nchar:  6   3rd Qu.: 4.200   3rd Qu.:0.4000   3rd Qu.:0.610  
##                  Max.   :12.000   Max.   :1.1000   Max.   :6.170  
##        Clase            Hora           Suelo            CAT     
##  Length   :316   Length   :316   Length   :316   Length   :316  
##  N.unique :  2   N.unique :  2   N.unique :  3   N.unique :  3  
##  N.blank  :  0   N.blank  :  0   N.blank  :  0   N.blank  :  0  
##  Min.nchar:  5   Min.nchar:  3   Min.nchar:  3   Min.nchar:  4  
##  Max.nchar:  5   Max.nchar:  5   Max.nchar:  6   Max.nchar:  4  
## 
par(mfrow = c(2,2))
plot(regresion_log)


Fig. 9 Diagnóstico del modelo de regresión en función al logaritmo del peso

Debido a las limitaciones observadas en el modelo lineal, se aplicó una transformación logarítmica utilizando la función:

\[ log(peso)=B_0 + B_1log(largo) + B_2log(ancho) \]

La transformación logarítmica permite convertir relaciones potencialmente exponenciales o alométricas en relaciones lineales, reduciendo la asimetría de los datos y la influencia de observaciones extremas. El análisis de los gráficos diagnósticos del modelo logarítmico, visualizados en la figura 9, mostraron una mejora en la distribución general de los residuos; sin embargo, los residuos aún presentan desviaciones respecto a la normalidad y una dispersión que no es completamente homogénea a lo largo de los valores ajustados. Esto indica que, aunque el modelo transformado describe mejor la relación biológica entre las variables, algunos supuestos de la regresión lineal continúan incumpliéndose parcialmente. Por lo tanto, las predicciones obtenidas deben interpretarse con cautela y como aproximaciones estadísticas más que como valores exactos.

A partir de este modelo logarítmico se realizó una predicción utilizando un individuo hipotético con 13 cm de largo y 0.5 cm de ancho. La sustitución de estos valores en la ecuación generó una predicción de:

\[ log(Peso)=(-1.29796)+((log(13) imes1.43756)+(log(0.5) imes1.55153)) \]

\[ log(Peso)=(-1.29796)+((2.56495 imes1.43756)+(-0.69315 imes1.55153)) \]

\[ log(Peso)=1.31 \]

Posteriormente, se aplicó la transformación inversa del logaritmo:

\[ Peso=e^{1.31} \]

obteniéndose un peso estimado de 3.72 g. Este resultado sugiere que un miriápodo con dichas dimensiones corporales tendría un peso cercano a 4 g según la tendencia observada en los datos, respaldando la utilidad del modelo logarítmico para realizar estimaciones biológicamente coherentes dentro del rango de valores analizados.

A continuación se realizaron dos modelos de análisis de varianza para evaluar cómo fluctúan las medias poblacionales de cada sitio y subfilo con respecto al largo y peso. Por último se evaluó la diversidad y uniformidad de las familias aplicando a cada uno de los sitios las pruebas de Shannon-Weiner y de Índice de diversidad de Simpson, según los métodos descritos por Pla, 2006 y Moreno, 2001.

Utilizando pruebas de análisis de varianza se compararon las fluctuaciones del largo de los miriápodos, al evaluar los supuestos del modelo, para la normalidad de los residuos se obtuvo un valor de p = 9.733e-15, lo que significa que estos no presentan normalidad, esto se ve reflejado también en la prueba de balance de datos, donde se denota que están altamente desbalanceados puesto que la masa de datos de diplópodos en San Fernando es de 241, mientras que los diplópodos de San Rafael son apenas 38, junto con 37 quilópodos; seguidamente, la evaluación de homocedasticidad Cochran para datos no paramétricos donde el valor de p fue de 4.158e-13, confirmando que en los sets de datos las varianzas se comportan diferente, visto esto se optó por el análisis de varianzas para datos no paramétricos de Kruskal-Wallis para el cual se obtuvo un p= 3.182e-12 con dos grados de libertad, indicando que efectivamente alguna de las medias es diferente a las otras, para comprobar cuál de los sets es responsable de las diferencias se hizo el test de anova post-hoc LSD (Least Significant Difference), para el cual se obtuvo que todos los sets de datos tenían diferencias significativas entre ellos, siendo mayormente significativa la diferencia entre diplópodos de ambos sitios (p= 9.6e-08), seguido por la diferencia entre diplópodos y quilópodos de San Rafael de poás (p= 2.3e-06) y por último la diferencia entre diplópodos de San Fernando y quilópodos de San Rafael de Poás (p= 1.4e-07), esto se refleja en el gráfico de la figura 10.


Fig. 10 Boxplot del ANDEVA de largo de miriápodos, SF_D: Diplópodos San Fernando, SR_C: Quilópodos San Rafael, SR_D: Diplópodos de San Rafael

Posteriormente se realizó un análisis de la varianza para la variable de peso de los miriápodos, siendo nuevamente verificados los supuestos. Para la normalidad según Shapiro, se obtuvo un valor de p = 2.2e-16 para los residuos utilizando una prueba Shapiro-Wilk, indicando así que los mismos no presentan una distribución normal. Como complemento se verificó que los datos se encuentran desbalanceados, como fue observado anteriormente para la variable del largo de los miriápodos. Al evaluar homocedasticidad con la prueba de Cochran para datos no paramétricos, se obtuvo un valor de p = 2.2e-16, el cual indica diferencias significativas entre el comportamiento de las varianzas. Visualizando este resultado, se realizó un segundo análisis de varianzas no paramétrico mediante Kruskal-Wallis, del cual se obtuvo un valor de p = 1.112e-13 con 2 grados de libertad, lo cual indica que el valor de alguna de las medias es diferente a las demás. Se procedió a verificar cual de los sets provoca las diferencias, esto de igual manera, mediante una prueba anova post-hoc LSD, donde se obtuvo que dos de los sets de datos tenían diferencias significativas entre sí, siendo esta entre diplópodos de San Fernando y quilópodos de San Rafael de Poás (9.4e-15), asimismo, se evidenció que no existen diferencias significativas entre los pesos de diplópodos de San Fernando y San Rafael de Poás (0.49). Lo anteriormente descrito se visualiza en la figura 11.
Fig. 11 Boxplot del ANDEVA de peso de miriápodos, SF_D: Diplópodos San Fernando, SR_C: Quilópodos San Rafael, SR_D: Diplópodos de San Rafael 

## Cargando paquete requerido: permute
ShannonSF <- diversity(IndicesSF$`Cantidad de individuos`, index = "shannon")
ShannonSF
## [1] 0.2727059
SimpsonSF <- diversity(IndicesSF$`Cantidad de individuos`, index ="simpson" )
SimpsonSF
## [1] 0.1184553
ShannonSPP <- diversity(IndicesSPP$`Cantidad de Individuos`, index = "shannon")
ShannonSPP
## [1] 1.636532
SimpsonSPP <- diversity(IndicesSPP$`Cantidad de Individuos`, index = "simpson")
SimpsonSPP
## [1] 0.736
sum(IndicesSPP$`Cantidad de Individuos`)
## [1] 75

Tabla 7.
Resultados obtenidos de aplicar índices de diversidad a ambos sitios de estudio

Sitio Prueba Resultado
SF Shannon-Weiner 0.273
SF Simpson (1-D) 0.118
SRP Shannon-Weiner 1.637
SRP Simpson (1-D) 0.736

Nota. SF: San Fernando, SRP: San Rafael de Poás.

Discusión

Las diferencias encontradas en la distribución de clases taxonómicas entre ambos sitios de estudio se puede llegar a una razón principal, los diplópodos han sido ampliamente constatados como una clase plástica en sus hábitats, si bien, al igual que los quilópodos, suelen buscar cierta humedad ambiental, hay registros de especies de diplópodos presentes hasta en desiertos, sin embargo, como regla general se les puede esperar con mayor frecuencia en lugares con amplia disponibilidad de hojarasca (Hopkin & Read, 1992), como lo son los bosques con especies caducifolias, esto es altamente compatible con el entorno estudiado en Sámara, mientras que la zona de estudio en Poás, si bien poseía parches boscosos, gran parte del territorio era cafetal, que debido al paso de trabajadores se limpia a menudo de la hojarasca para evitar accidentes con reptiles, por esto gran parte de los individuos de esta clase vistos en Poás se observaron directamente sobre las plantas; ahora, basado en las observaciones de la clase chilopoda, su distribución y preferencias de hábitat son más compatibles con el área templada, húmeda y con abundancia otros invertebrados pequeños para alimentarse, además de que varios de estos tienen la habilidad de hundirse en la tierra para enfriarse (Hickman at al, 2024), muchos de estos se vieron debajo de basura, troncos o piedras, incluso teniendo que escarbar para poder encontrarlos, con las especies de escolopendras siendo las que sí se vieron más activamente expuestos, por supuesto que no se descarta la presencia de quilópodos, pues hay especies de esta clase que incluso habitan zonas costeras (Korpella, 2017), no obstante es muy probable que el estado de sequía de la región hiciera que los individuos se movieran a zonas mucho más profundas del bosque a las que el estudio no llegó.

Respecto a la prueba de Shapiro aplicada para las variables de peso, largo y ancho, se evidenció que las tres no poseen una distribución normal. Lemba et al. 2026 establece que, en ecosistemas naturales, se encuentran patrones de distribución basados en parches para diversos organismos, enfocándose principalmente en Diplópodos, los cuales poseen preferencias marcadas en cuanto a la humedad, temperatura e intensidad de la luz. Este tipo de distribuciones, sumada a las diferencias ambientales entre ambos sitios de muestreo, influyen directamente en la falta de normalidad dentro de los datos. Al encontrarse agrupadas las muestras por localidad, se crea una combinación de “subpoblaciones” que están expuestas a condiciones de suelo y microclimas contrastantes, lo cual genera una asimetría considerable en las dimensiones corporales y masa de los individuos a los cuales se les efectuaron las respectivas mediciones.
La ausencia de diferencias significativas en la longitud corporal de los individuos entre San Fernando de Samara y San Rafael de Poás, en contraste con las diferencias significativas que fueron encontradas en el peso y el ancho mediante la prueba de Wilcoxon, concuerda con lo reportado por Milovanović et al. (2025). Dichos autores demostraron que las dimensiones lineales esqueléticas de los miriápodos suelen mantenerse conservadas entre distintas localidades, mientras que variables volumétricas y relacionadas a la forma, como lo son el peso y el ancho de un individuo, presentan una mayor plasticidad fenotípica como respuesta a los cambios en el ambiente de cada zona. Por lo tanto, la variabilidad vista en cuanto al peso y el ancho refleja ajustes condicionados por el microclima de cada sitio, los cuales no alteran la escala lineal de la longitud corporal de las especies presentes.

La marcada disparidad en las capturas entre Diplopoda (88,29%) y Chilopoda (11.71%) concuerda con los patrones reportados por Navas et al. (2025) para el Neotrópico, donde los milpiés representaron más del 70% de la abundancia total en los muestreos de miriápodos. Esta diferencia responde a las restricciones fisiológicas y climáticas asociadas al balance hídrico de ambos grupos, donde los diplópodos poseen adaptaciones tegumentarias y conductuales que les permiten conservar la humedad y mantener el forrajeo sobre la hojarasca (David, 2009), mientras que la clase Chilopoda muestra mayor susceptibilidad a la deshidratación en la superficie (Navas et al., 2025). Ante la estacionalidad seca o la variación del microclima local, los quilópodos evidencian una tendencia a migrar hacia estratos del suelo más profundos, o a buscar refugios en rocas o troncos (Bogyó et al., 2015). Tal conducta reduce las posibilidades de capturar organismos de la clase Chilopoda en relación a Diplopoda (Lemba et al., 2026; Stašiov et al., 2025)

Referente a la caracterización del dimorfismo sexual en la muestra de diplópodos, se observó que el 58.35% de los individuos no presenta dimorfismo sexual externo evidente, mientras que el 41.64% restante presentó una diferenciación clara. La mayor proporción de organismos catalogados con dimorfismo ausente se ve atribuida principalmente a la presencia de los estadios juveniles de muchos de los organismos en las capturas de campo, los cuales pueden llegar a tener una diferenciación progresiva a lo largo de sus mudas hasta que alcancen la madurez completa (Stojanović et al., 2023) A su vez, los individuos con dimorfismo presente corresponden a especímenes adultos que evidencian un marcado dimorfismo sexual de tamaño (SSD), en el cual las hembras constituyen el sexo de mayor tamaño y volumen corporal. Según Milovanović et al. (2025) e Illić et al. (2019), el patrón descrito es frecuente en los diplópodos, y se ve promovido por la selección por fecundidad, donde un cuerpo y ancho abdominal mayor otorga ventajas directas a las hembras, permitiendo que logren albergar mayor cantidad de huevos y a su vez, invertir más recursos energéticos en su descendencia.

Se realizó una prueba de chi cuadrado de dos vías con el fin de analizar si la densidad de individuos se ve afectada por el sitio de muestreo de estos, esta prueba evidenció una asociación significativa entre el sitio de muestreo y la frecuencia de los individuos registrados por cuadrante con un valor de p correspondiente a 0.0004998, con este resultado se descartó la hipótesis nula y acepta la alternativa la cuál indica que la distribución de organismos difiere entre San Fernando (SF) y San Rafael de Poás (SRP). En San Fernando se obtuvieron hasta un máximo de 14 individuos por cuadrante con distribuciones mucho más homogéneas, mientras que en San Rafael de Poás predominó una menor cantidad de organismos por unidad de muestreo. Estas diferencias se ven relacionadas con las características ambientales de cada sitio, como la humedad, la temperatura y la disponibilidad de recursos y materia orgánica. En este sentido, (Lam et al,2024) señalan que la calidad de la hojarasca y su acondicionamiento dependiendo de estos factores influyen en el consumo de material orgánico por parte de los diplópodos, fortaleciendo la hipótesis alternativa donde los sitios con una mayor temperatura representarían una densidad mayor de individuos gracias a una mayor disponibilidad de alimento debido a una descomposición más acelerada de la materia orgánica por los procesos microbianos y fúngicos a los cuales benefician este ambiente . Así mismo, (Zeng et al, 2024) destaca la contribución de los invertebrados en la descomposición de la hojarasca, con una mayor tendencia en sistemas tropicales. Por lo tanto las diferencias observadas podrían estar relacionadas con la disponibilidad y calidad de los recursos en ambos sitios, sin embargo sería necesario realizar más mediciones ambientales, de pH, zonas de mayores densidades de materia orgánica el grado de descomposición y la calidad de la misma entre otras como el tipo de suelo para determinar con una mayor certeza cuáles factores influyen directamente en las distribuciones registradas.

El análisis de chi cuadrado de una sola vía evidenció diferencias significativas en la proporción sexual de los diplópodos en comparación con los resultados esperados de obtener una proporción de -1:1 entre machos y hembras (X²= 13.337, gl=1; p-Value= 0.0002602). Las hembras registraron un predominio mayor con un total de 170 individuos, lo cual, representa un porcentaje del 60.93%, mientras tras el número de machos registrados equivale a 109 individuos, representante tan solo un 39.07% de los individuos totales identificados (279), lo que representa aproximadamente 1,6 hembras por cada macho. Esta desviación entre las proporciones de los sexos en los individuos se podría ver relacionada con los comportamientos y horas de actividad de los individuos, la longevidad o los hábitos de cada individuo, según una investigación (Žabova et al, 2026). se obtuvieron resultados donde se aprecia que en la mayoría de insectos las hembras tienen una longevidad mayor que los machos lo cual podría ser una razón para encontrar más de estas, sin embargo no se han realizado estudios recientes de este factor en Miriápodos, el último registro de este hecho fue descrito por (Pinheiro et al, 2009), en el cual se estudiaron las especies Poratia Salvador Golovatch y Sierwald en Brasil (clima Tropical) donde la mayoría de hembras conseguidas en cautiverio mostraron una longevidad promedio de 147.2 días más que los machos, además de que al conseguir individuos juveniles mediante los huevos de las hembras capturadas las cuales se reprodujeron mediante partenogénesis, se obtuvo una proporción de 1 macho: 139 hembras lo que indica una proporción mucho mayor de hembras en esta especie. Sin embargo, los resultados obtenidos no permiten establecer una causa directa para este predominio en las especies analizadas. Es importante recalcar que para este análisis solamente se consideraron los individuos a los que se les logró identificar con un sexo a partir de la presencia de los gonopodios en machos, excluyendo aquellos cuya identificación no fue posible, esta decisión fue tomada con el fin se reducir los errores asociados a esta clasificación de género y trabajar con individuos cuya identificación si se sabía con certeza absoluta. Las pruebas realizadas para analizar la correlación entre el peso y otras dos variables se desarrollaron con el modelo de Spearman a causa del déficit de normalidad en la prueba. Además de la evaluación para determinar el modelo se analizó la desviación, siendo que el peso presentó la desviación más marcada (W = 0.586), seguido por el ancho (W = 0.814) y el largo (W = 0.845). Estos resultados sugieren la presencia de asimetría y posibles valores extremos, situación frecuente en organismos naturales donde coexisten individuos juveniles y adultos dentro de una misma muestra. Se reveló una asociación positiva muy fuerte entre el peso y el largo corporal (ρ = 0.956; p < 2.2 × 10⁻¹⁶). Se esperaba que existiera una correlación positiva entre estas variables respecto al peso ya que estas influyen en este. A nivel biológico se evidencia que los individuos que poseen un mayor tamaño corporal (longitud) y una mayor robustez la cual se relaciona al volumen de los individuos (ancho) van a presentar una mayor biomasa. Entre estas dos variables es más utilizada la longitud de los individuos ante su ancho debido a que usualmente se analiza el patrón de crecimiento (Sohlström et al., 2018).

El modelo de regresión lineal múltiple demostró que tanto el largo como el ancho contribuyen en la predicción del peso corporal en los miriápodos, esto demostrado por los coeficientes de las variables los cuales fueron altamente significativos. El coeficiente asociado al largo indica que, manteniendo constante el ancho, un incremento de un centímetro en la longitud corporal produciría un aumento promedio de aproximadamente 0.10 g en el peso. Así mismo, el coeficiente del ancho sugiere que un incremento de una unidad de ancho genera un efecto mucho mayor sobre la masa corporal. El coeficiente de determinación obtenido (R² = 0.657) indica que aproximadamente el 65.7 % de la variación observada en el peso puede explicarse mediante las variables largo y ancho. El valor del estadístico F (F = 505.5; p < 2.2 × 10⁻¹⁶) demuestra que el modelo es estadísticamente significativo y posee capacidad predictiva real. A pesar de estos resultados también se demostró algunas limitaciones del estudio con los análisis diagnósticos. La prueba de Shapiro-Wilk aplicada a los residuos mostró una fuerte desviación respecto a la normalidad (W = 0.485; p < 2.2 × 10⁻¹⁶). Esto indica que parte de la variabilidad presente en el peso no está siendo completamente explicada por el modelo. La prueba de heterocedasticidad mostró un resultado significativo (χ² = 2127.7; p < 2.2 × 10⁻¹⁶), evidenciando la existencia de heterocedasticidad. No obstante, no puede afirmarse que este supuesto se cumpla completamente. La prueba Durbin-Watson (DW = 1.825) se encuentra cercana al valor ideal de 2, lo que demuestra la presencia de independencia de los residuos y una autocorrelación baja.

A partir de estas pruebas y sus resultados correspondientes se establece que la relación entre el tamaño corporal de los miriápodos y su peso muestra una escala alométrica. Estas características son propias de los artrópodos y otros invertebrados. En cuanto a la predicción realizada para un individuo de 13 cm de largo y 0.5 cm de ancho, la estimación es que produzca un peso de 3.72 g. Según lo mencionado anteriormente esta predicción no va a ser del todo precisa, aunque el valor de R² obtenido indica que el modelo constituye una herramienta útil para estimar biomasa en individuos donde no sea posible obtener mediciones directas de peso. A nivel mundial existe evidencia que los principios alométricos son aplicables a los análisis en este estudio (Ruiz-Lupión et al., 2020; Versluijs et al., 2023).

Para el análisis de varianza Kruskal-Wallis del largo de los individuos, la diferencia más sencillamente explicable es la vista entre quilópodos y diplópodos de ambos sitios, está bien documentado en la literatura que la cantidad de segmentos suele ser mayor en los diplópodos, con un promedio de 75~100 segmentos, mientras que para los quilópodos la cantidad de segmentos no suele sobrepasar los 60, para poder acomodar esta diferencia en número de segmentos, es necesario que el plan morfológico de los milpiés sea, por lo general, un poco más largo (Mitchell & Zachariah, 2009). Además de esto, la diferencia podría estar relacionada a las necesidades de cada clase, el plan corporal de la clase chilopoda está hecha para ser veloz y escurridiza como ventaja para la carnivoría, mientras que para los milpiés, cuyo mecanismo de defensa principal es el enrollamiento defensivo y la secreción de sustancias deflectivas, les conviene que tengan una longitud que permita que sus órganos más importantes queden al centro del espiral, mientras que partes menos cruciales queden como plano principal para el depredador (Hickman et al, 2024).

Por otro lado, las diferencias significativas encontradas entre la misma clase en ambos sitios tiene implicaciones varias, actividades de carácter antropológico en la zona de Poás–visto que la zona de estudio era en gran parte un cafetal; como el uso de agroquímicos, cambios en la composición de la vegetación y el suelo pueden afectar el desarrollo de los individuos, en este caso podría traducirse como un cuerpo más largo en San Rafael, más para determinar esto se requieren estudios químicos los cuales se salen de las posibilidades del estudio, aún así, es un punto relevante a considerar en un futuro (de Albuquerque et al, 2026 ). Un punto más explicable con la metodología presentada es el análisis de la regla temperatura-tamaño, para sámara, en este año la temperatura promedio ha ido de 27~32°C, mientras que para Poás ha sido de 19~27°C (Weather Atlas, 2026), en este sentido, podría considerarse que la altitud junto con esta variación de temperatura estén influyendo sobre las características morfológicas de los diplópodos, donde debido a su carácter ectotermo, a menor temperatura mayor tamaño corporal para mayor almacén de energía, aún así, el área de estudio de las reglas de Bergmann y Allen varía mucho en sus resultados dependiendo del grupo de artrópodos en general, y en miriápodos hay una escasez de investigación en el área, por lo que es una rama que aún debe fortalecerse a futuro (Zuo et al, 2012; Shelomi & Zeuss, 2017). También se debe destacar que, hablando meramente del modelo estadístico utilizado en el análisis, los resultados podrían haber sido influenciados porque el individuo más largo fue encontrado en Poás, incurriendo en cierto nivel de sesgo en el cálculo como dato aislado.

Para el análisis del peso de los miriápodos mediante la prueba de Kruskal-Walli y la comparación post-hoc LSD (Least Significant Differences), la diferencia más marcada es la referente a los individuos de San Fernando y los quilópodos de San Rafael de Poás, con un valor de p de 9.4e-15. Tal diferencia responde a los planes corporales y roles en la cadena trófica de ambas clases: los diplópodos poseen un cuerpo cilíndrico con un exoesqueleto fuertemente calcificado y un tracto digestivo adaptado al proceso continuo de sedimentos y hojarasca (David, 2009; Stašiov et al., 2025), lo que incrementa sustancialmente la masa por individuo; mientras que los quilópodos son depredadores que requieren de agilidad, con un cuerpo menos calcificado y mayormente adaptado a la caza activa (Bogyó et al., 2015; Navas et al., 2025)

Por otro lado, la ausencia de diferencias significativas en el peso de diplópodos entre San Fernando y San Rafael de Poás, con un valor de p de 0.49, demuestra que la masa corporal individual se mantiene en rangos cercanos, aún con las variaciones de altitud, temperatura y del propio manejo del sustrato (Stašiov et al., 2025). Este resultado es concordante con lo reportado por Milovanović et al. (2025), quienes demostraron mediante el uso de técnicas morfométricas que las dimensiones estructurales y la masa corporal general en la clase Diplopoda tiende a mantenerse sin variaciones significativas entre poblaciones pertenecientes a un mismo gremio funcional. Aunque el microclima, la humedad o el hábitat de cada sitio puedan generar una plasticidad morfológica en otros caracteres de los individuos (David, 2009; Lemba et al., 2026; Stašiov et al., 2021), la masa se mantiene en una escala constante que no se encuentra alterada por los cambios de localidad.

Las pruebas de diversidad tuvieron resultados esperados conforme a los datos obtenidos, primeramente la prueba de diversidad de Simpson establece la probabilidad de que tomar dos individuos de una muestra, estos pertenezcan a la misma especie, en este caso a la misma familia, siendo este índice utilizado para evaluar dominancia de especies (Moreno, 2001), para la muestra poblacional de Sámara, se obtuvo un 88.2% de probabilidad de que dos individuos pertenezcan a la misma familia, mientras que para Poás fue de apenas un 26.4%, esto es un resultado sencillo de ver al comparar las densidades de los datos, puesto que de la muestra de San Fernando más de 200 de los individuos eran de la misma familia, con apenas 15 de otras familias, por lo que se evidencia que la población de miriápodos en este sitio está fuertemente dominada por la familia paradoxosomatidae, en cambio en la población de San Rafael las familias tuvieron una distribución más balanceada, siendo el índice más influenciado por las familias de spirobólidos y cryptopidos, ahora bien, se debe considerar que el tamaño de las muestras de cada sitio es muy desproporcionado, debido a que, como ya se comentó, los factores ambientales de cada hábitat propician a que la densidad poblacional se acumule más en la zona de Guanacaste (Bogyó et al, 2015), por esto no se descarta que si se llegara a hacer un estudio donde se asegure la misma cantidad mínima de individuos mínimos en cada sitio o se llevara a cabo un muestreo más extenso temporalmente, los resultados podrían variar. Continuando, el índice de Shannon-Weiner es más utilizado para evaluar la riqueza y equidad de la distribución de especies (abundancia relativa y especies presentes), en este caso el índice obtenido para San Fernando fue sumamente bajo, cercano a cero, nuevamente atribuyéndose a que solo se observaron 3 familias con una abundancia sumamente diferente para cada una, mientras que para Poás el índice fue medianamente cercano a dos, posicionándolo en un nivel de diversidad medio/bajo (Pla, 2006). Es sabido que este índice es sensible a especies raras por lo que es posible que el hecho de que en este último sitio varias familias tuvieran representaciones menores a 5, además de que la muestra fuese mucho más pequeña, influyeran en este resultado, por lo que nuevamente el estudio se fortalecería con la extensión de áreas y tiempo de muestreo. Un factor que pudo haber contribuido a los resultados de las pruebas de diversidad es la diferencia de temperatura entre sitios, pues en la literatura se ha visto por lo general una mayor cantidad de milpiés en sitios con temperaturas más altas y cobertura boscosa más densa (Bogyó et al, 2015), contrastantemente, debido al nivel de humedad que requieren los quilópodos, siendo este de hasta el 90%, es normal verlos más comúnmente en lugares más frescos que sean capaces de mantener este nivel de humedad ambiental por periodos prolongados de tiempo, oscilando temperaturas ideales de 20~25°C (Mitchell & Zachariah, 2009), vistas estos contrastes en necesidades intrínsecas de las clases, las diferencias en diversidad del subfilo de estudio cobran más sentido.

CONCLUSIONES

Los análisis realizados en las localidades de San Fernando de Sámara y San Rafael de Poás han permitido demostrar que características como el largo, peso y ancho pueden llegar a variar debido a las condiciones intrínsecas presentes en cada zona, donde se ve una tendencia aparente de ciertas características, como el largo, a verse más fuertemente influenciadas que otras como el peso. Similarmente, la desigualdad en densidad y diversidad poblacional en cada localidad puede deberse a las diferencias ambientales que propician un hábitat que cumple o no con las necesidades de cada grupo de estudio, como lo es la calidad de la hojarasca y la temperatura, que son factores beneficiosos para el crecimiento de la población y que según el grupo, las necesidades pueden ser diferentes. Así mismo, la alteración humana en ecosistemas similares a los estudiados es un factor determinante en la supervivencia y el desarrollo de las poblaciones de miriápodos, siendo factores como la remoción o eliminación de hojarasca, las alteraciones de temperatura y el posible uso de sustancias químicas en los alrededores, posibles detonantes de diferenciación.

En general, el estudio proporciona evidencia de la fuerte variabilidad poblacional que las características específicas de cada ecosistema pueden generar, incluso en grupos taxonómicos altamente relacionados, aún así, es necesario seguir reforzando los estudios para este subfilo, siendo recomendable aumentar el tiempo de estudio y el territorio analizado, además de tomar en cuenta variables que no se tomaron en consideración en el análisis estadístico como composición vegetal y del suelo, condiciones climáticas y relaciones de depredación presentes.

REFERENCIAS

Achim, Z. & Torsten, H. (2002). Diagnostic Checking in Regression Relationships. R News, 2(3), 7-10. https://CRAN.R-project.org/doc/Rnews/

Bogyó, D., Magura, T., Nagy, D. D. & Tothmeresz, B. (2015). Distribution of millipedes (Myriapoda, Diplopoda) along a forest interior – forest edge – grassland habitat complex. ZooKeys, 510, 181-195. https://doi.org/10.3897/zookeys.510.8657

Comunidad de Inaturalist. (2026). Observaciones de Miriápodos. Inaturalist. https://www.inaturalist.org/observations?verifiable=true&taxon_id=144128

David, J.-F. (2009). Ecology of millipedes (Diplopoda) in the context of global change. Soil Organisms, 81(3), 719. ISSN: 1864-6417

de Albuquerque Pereira, B., Ilha, C., Domingos Brescovit, A., Moretti Iniesta, L. F., Salvador Bouzan, R., Iuñes de Oliveira Filho, L. C., Baretta, D., Chagas Junior, A., de Souza Rezende, R. & Riviera Duarte Maluche Baretta, C. (2026). Patterns of taxonomic and functional diversity of Diplopoda and Chilopoda in land use systems. Pedobiologia, 114, 151115. https://doi.org/10.1016/j.pedobi.2026.151115

Fernández-Landa, A., Fernández-Moya, J., Algeet-Abarquero, N., Guillén-Climent, M. L., Pedroni, L., Marchamalo-Sacristán, M., Bonatti, J., Mencuccini, M., y Gilabert, M. A. (2017). High resolution biomass mapping in tropical forests with LiDAR-derived digital models: Poás Volcano National Park (Costa Rica). iForest - Biogeosciences and Forestry, 10(5), 899-905. https://doi.org/10.3832/ifor1744-009

Fondo Nacional de Financiamiento Forestal. (2024). Restauración y conservación de áreas de Bosque Seco en el Pacífico Norte de Costa Rica. Ministerio de Ambiente y Energía, Costa Rica. https://www.fonafifo.go.cr/media/4150/proyecto-bosque-seco.pdf

Fox, J. & Weisberg, S. (2019). An R Companion to Applied Regression. Sage, Thousand Oaks CA. https://www.john-fox.ca/Companion/

Hickman, C. P. Jr., Keen, S. L., Eisenhour, D. J., Larson, A. & l’Anson, H. (2024). Integrated Principles of Zoology. McGraw Hill.

Hopkin, P. S. & Read, J. H. (1992). The biology of Millipedes. Oxford University Press. ISBN: 0-19-857699-4. https://www.scribd.com/document/1010337842/The-Biology-of-Millipedes

Ilić, B. S., Vujić, V. D., Jovanović, Z. S., Pavković-Lučić, S. B., Dudić, B. D., Lučić, L. R. & Makarov, S. E. (2019). Sexual dimorphism in some morphological traits of three European millipedes (Diplopoda, Julida, Julidae). Animal Biology, 69(4), 483-496. https://doi.org/10.1163/15707563-20191113

Janzen, D. H. (1988). Tropical dry forests: The most endangered major tropical ecosystem. En E. O. Wilson & F. M. Peter (Eds.), Biodiversity (pp. 130-137). National Academies Press. https://www.ncbi.nlm.nih.gov/books/NBK219281/

Komsta, L. (2022). Outliers: Tests for Outliers. https://doi.org/10.32614/CRAN.package.outliers

Korpella, R. (2017). Habitats Of Centipedes. Sciencing, Estados Unidos de América. https://www.sciencing.com/habitat-centipede-8518236/

Lam, W. N., Slade, E. M. & Wardle, D. A. (2024). Effects of leaf litter traits on terrestrial isopod and millipede consumption, assimilation and growth. Functional Ecology, 38(5), 1018–1031. https://doi.org/10.1111/1365-2435.14520

Lemba, R., Zainal, S., Rafiqa, R., Dhafir, F., Masrianih, M. & Trianto, M. (2026). Abundance and Distribution of Millipedes (Diplopoda) in Wera Nature Recreation Park. Biology, Medicine, & Natural Product Chemistry, 15(1), 597-603. https://doi.org/10.14421/biomedich.2026.151.597-603

Levi, H. W. & Levi, R. L. (2002). Spiders and their kin. St. Martin’s Press. ISBN: 1-58238-156-9.

Lizano Araya, M. A. & Lizano Rodríguez, O. G. (2023). Escenarios de inundación ante el aumento del nivel del mar por Cambio Climático, para las Playas del Coco, Tamarindo y Sámara, Costa Rica. Entorno Geográfico, 25, 1–24. https://doi.org/10.25100/eg.v0i25.11903

Milovanović, J., Vujić, V., Jovanović, Z., Dudić, B., Lučić, L., Ilić, B. & Tomić, V. (2025). Morphological Variability in Two Populations of Millipede Glomeris hexasticha Brandt, 1833 (Diplopoda: Glomerida: Glomeridae): Traditional and Geometric Morphometric Approaches. Acta Zoologica, azo.70011. https://doi.org/10.1111/azo.70011

Moreno, E. C. (2001). Métodos para medir la biodiversidad. Sociedad Entomológica Aragonesa, 49(3-4), 1-86. ISSN: 1576 – 9526.

Municipalidad de Poás. (2026). Datos de Cantón. Poás Gobierno Local, Poás, Alajuela, Costa Rica. https://munipoas.go.cr/index.php/conozcanos/mi-canton/datos-de-canton.html

Navas-S., G. R., Carrillo-Pallares, M., Zapata-M., W., Martínez-Torres, D. & Bermúdez-Tobón, A. (2025). Diversity of centipedes and millipedes (Arthropoda: Myriapoda) in tropical dry forest fragments in Bolívar, Colombia. Revista de Biología Tropical, 73, e2025546. https://doi.org/10.15517/9w79j526

Oksanen, J., Simpson, G., Blanchet, F., Kindt, R., Legendre, P., Minchin, P., O’Hara, R., Solymos, P., Stevens, M., Szoecs, E., Wagner, H., Bedward, M., Bolker, B., Borcard, D., Carvalho, G., De Caceres, M., Durand, S., Evangelista, H., Hannigan, G., Hill, M., Lahti, L., Martino, C., Ouellette, M., Ribeiro Cunha, E., Smith, T., Stier, A., Ter Braak, C. & Weedon, J. (2026). Vegan: Community Ecology Package. https://doi.org/10.32614/CRAN.package.vegan

Pinheiro, T. G., Marques, M. I. & Battirola, L. D. (2009). Life cycle of Poratia salvator (Diplopoda: Polydesmida: Pyrgodesmidae). Zoologia, 26(4), 658–662. https://doi.org/10.1590/S1984-46702009000400010

Pla, L. (2006). Biodiversidad: Inferencia basada en el índice de Shannon y la riqueza. Interciencia, 31(8), 583-590. ISSN: 0378-1844.

Ruiz-Lupión, D., Gómez, J. M. & Moya-Laraño, J. (2020). Mass–length allometry covaries with ecosystem productivity at a global scale. Global Ecology and Biogeography, 29(1), 87-101. https://doi.org/10.1111/geb.13010

R Core Team. (2026). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. https://doi.org/10.32614/R.manuals

Shelomi, M. & Zeuss, D. (2017). Bergmann’s and Allen’s Rules in Native European and Mediterranean Phasmatodea. Frontiers in Ecology and Evolution, 5(25). https://doi.org/10.3389/fevo.2017.00025

Sistema Nacional de Áreas de Conservación. (2025). Parque Nacional Volcán Poás. Alajuela, Costa Rica. https://www.sinac.go.cr/es/ac/accvc/pnvp/paginas/default.aspx

Sohlström, E. H., Marian, L., Barnes, A. D., Haneda, N. F., Scheu, S., Rall, B. C., Brose, U. & Jochum, M. (2018). Applying generalized allometric regressions to predict live body mass of tropical and temperate arthropods. Ecology and Evolution, 8(24), 12737-12749. https://doi.org/10.1002/ece3.4702

Stašiov, S., Diviaková, A., Pätoprstý, V. & Svitok, M. (2025). The impact of submontane meadows management on the structure of millipede (Diplopoda) communities. Basic and Applied Ecology, 88, 94–103. https://doi.org/10.1016/j.baae.2025.09.002

Stojanović, D. Z., Vujić, V. D., Jovanović, Z. S., Milovanović, J. Z., Dudić, B. D., Ilić, B. S. & Makarov, S. E. (2023). Morphological variation during post-embryonic development in the centipede Lithobius melanops: Traditional and geometric morphometrics approaches. Contributions to Zoology, 92(4), 316-348. https://doi.org/10.1163/18759866-bja10044

Stork, N. E. (2018). How many species of insects and other terrestrial arthropods are there on Earth? Annual Review of Entomology, 63, 31-45. https://doi.org/10.1146/annurev-ento-020117-043348

Tamarindo Estuary Tours Mi Rey. (s. f.). Bosque seco tropical de Guanacaste: Árboles y fauna. https://www.tamarindoestuary.com/guides/tropical-dry-forest-guanacaste/?lang=es

Universidad Estatal a Distancia (UNED). (s. f.). Poás: El volcán que respira vida. Un laboratorio natural de biodiversidad en las alturas. https://www.uned.ac.cr/ecologiaurbana/notas-informativas/783-poas-el-volcan-que-respira-vida-un-laboratorio-natural-de-biodiversidad-en-las-alturas

Vargas, M. (2016). Myriapoda of Costa Rica (INBio), Version 1.4. Instituto Nacional de Biodiversidad (INBio), Costa Rica. Occurrence dataset. https://doi.org/10.15468/pftzpd

Versluijs, T. S. L., Zhemchuzhnikov, M. K., Kutcherov, D., Roslin, T., Martin Schmidt, N., van Gils, J. A. & Reneerkens, J. (2023). Site-specific length-biomass relationships of arctic arthropod families are critical for accurate ecological inferences. PeerJ, 11, e15943. https://doi.org/10.7717/peerj.15943

Weather Atlas. (2026). Clima y Previsión Meteorológica Mensual, Poás, Costa Rica. Emiratos Árabes Unidos. https://www.weather-atlas.com/es/costa-rica/poas-clima

Weather Atlas. (2026). Clima y Previsión Meteorológica Mensual, Sámara, Costa Rica. Emiratos Árabes Unidos. https://www.weather-atlas.com/es/costa-rica/samara-clima

Wickham, H., François, R., Henry, L., Müller, K. & Vaughan, D. (2026). Dplyr: A Grammar of Data Manipulation. https://doi.org/10.32614/CRAN.package.dplyr

Wickham, H. (2016). ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. ISBN: 978-3-319-24277-4. https://ggplot2.tidyverse.org

Žabova, B., Knapp, M., Kaasik, A. & Teder, T. (2026). Shared ancestry, sexual size dimorphism, and reproductive ecology jointly shape sex differences in adult lifespan in insects [Data set]. Dryad. https://doi.org/10.5061/dryad.xksn02vxm

Zachariah, T. & Mitchell, M. A. (2009). INVERTEBRATES. En Manual of Exotic Pet Practice. Elsevier. https://doi.org/10.1016/B978-141600119-5.50006-8

Zeng, X., Gao, H., Wang, R., Majcher, B. M., Woon, J. S., Wenda, C., Eggleton, P., Griffiths, H. M. & Ashton, L. A. (2024). Global contribution of invertebrates to forest litter decomposition. Ecology Letters, 27(4), e14423. https://doi.org/10.1111/ele.14423

Zuo, W., Moses, M. E., West, G. B., Hou, C. & Brown, J. H. (2012). A general model for effects of temperature on ectotherm ontogenetic growth and development. Proceedings of the Royal Society B: Biological Sciences, 279(1734), 1840–1846.