Introducción

Los ecosistemas de agua dulce ocupan menos del 1 % de la superficie terrestre, sin embargo, albergan una alta diversidad de especies, entre las que se encuentra cerca del 43 % de todas las especies de peces descritas (^{Strayer & Dudgeon, 2010}). En los Andes tropicales se encuentra la mayor diversidad de peces de alta montaña (^{Miranda et al., 2021}) y un alto número de endemismos (^{Myers et al., 2000}); no obstante, en el norte de esta cordillera existe una mayor densidad poblacional humana y una mayor transformación de las coberturas naturales (^{Etter & van Wyngaarden, 2000}). A esto se le suman otras amenazas globales que prevalecen sobre los ecosistemas de agua dulce a nivel global, como la sobreexplotación de recursos, incluida la extracción y explotación de agua, la contaminación, el cambio climático y las especies invasoras (^{Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services, 2019}). Lo anterior amenaza la integridad ecológica, la disponibilidad de agua potable, la salud humana y los medios de subsistencia en los sistemas dulceacuícolas, poniendo en peligro su diversidad biológica (^{Bogardi et al., 2020}; ^{Dudgeon, 2019}), la cual ya ha sido catalogada como en alto riesgo de desaparecer (^{Maldonado-Ocampo et al., 2005}; ^{Myers et al., 2000}).

El estudio de la diversidad beta es útil para comprender los efectos de las modificaciones humanas en las comunidades biológicas (^{Dala-Corte et al., 2019}), al estimar el grado de cambio en la composición de la comunidad en relación con un gradiente complejo de entorno y evaluar los procesos que generan y mantienen la biodiversidad en los ecosistemas (^{Legendre & De Cáceres, 2013}; ^{Whittaker, 1960}). La diversidad beta se puede particionar en los fenómenos de reemplazamiento y anidamiento, esto permite evaluar si los cambios en la composición están dados por la sustitución de especies (reemplazamiento) y/o por la ganancia o pérdida de especies desagregadas de forma ordenada, siendo las biotas de sitios con menor riqueza, subconjuntos de las biotas de sitios con mayor riqueza (anidamiento) (^{Baselga, 2010}). Los fenómenos de reemplazamiento y anidamiento actúan de formas diferentes en las distintas regiones del mundo, según la latitud e historia climática (^{Leprieur et al., 2011}). En los sistemas andinos y en general en los ríos del neotrópico, el reemplazamiento suele ser el fenómeno dominante en la composición de las comunidades de peces (^{Carvalho et al., 2021}; ^{López-Delgado et al., 2020}). Este fenómeno ha sido relacionado con la clasificación ambiental (environmental sorting) (^{Baselga, 2010}; ^{Carvalho et al., 2021}; ^{López-Delgado et al., 2019}) la cual hace énfasis en la capacidad de los organismos para seleccionar y ocupar áreas con condiciones ambientales adecuadas (^{López-Delgado et al., 2020}), por lo que las diferencias ambientales entre parches promoverán una estructura heterogénea de las comunidades locales (^{Anderson & Maldonado-Ocampo, 2011}; ^{Jaramillo-Villa et al., 2010}). Adicionalmente, es posible identificar la contribución de los sitios individuales y la contribución de las especies individuales a la diversidad beta (^{Legendre & De Cáceres, 2013}). En este sentido, la primera representa el grado de singularidad de las unidades de muestreo en términos de su composición comunitaria y de sus particularidades ecológicas, mientras que la segunda indica la importancia relativa de cada especie sobre los patrones de diversidad beta, relacionándose con características como el grado de ocupación, la abundancia, la posición del nicho, la amplitud del nicho y los rasgos de cada especie (^{Heino & Grönroos, 2017}; ^{Legendre & De Cáceres, 2013}).

Por otra parte, gran parte de la información que se conoce sobre los ríos y su biota proviene de observaciones a escala espacial local, es decir, a partir de le evaluación de tramos cortos e inconexos que han generado una comprensión limitada de los procesos que allí suceden (^{Fausch et al., 2002}). Lo anterior subestima la complejidad de los ríos al no abordarlos como mosaicos complejos de múltiples tipos de hábitats y gradientes ambientales, incluyendo su carácter jerárquico (^{Fausch et al., 2002}), y al no reconocer su dependencia del paisaje que los circunda (^{Allan, 2004}; ^{Hynes, 1975}). Algunos estudios en ríos neotropicales han documentado que las comunidades de peces son afectadas principalmente por variables a escala de paisaje tales como los usos del suelo y la elevación (^{Herrera-Pérez et al., 2019}; ^{Pinto et al., 2006}). En otros, las características locales como la estructura física de la zona ribereña, las propiedades fisicoquímicas del agua, la morfología de la corriente, la composición del hábitat y la heterogeneidad mostraron tener una mayor influencia (^{Cortés-Hernández et al., 2023}; ^{Miranda et al., 2021}; ^{Roa-Fuentes et al., 2019}; ^{Roa-Fuentes & Casatti, 2017}). No obstante, las comunidades de peces pueden verse impactadas en igual medida por variables diversas como el uso del suelo, la hidrología, la geomorfología, y las características fisicoquímicas del hábitat (^{Cunico et al., 2012}). En los Andes tropicales se han realizado algunos trabajos sobre la afectación de variables a diferentes escalas sobre organismos acuáticos. Por ejemplo, en Colombia, Ecuador y Perú se evaluaron variables fisicoquímicas, hidromorfológicas, tipos de sustrato, índices de calidad ambiental y usos del suelo sobre la abundancia de las poblaciones de peces, donde las variables hidromorfológicas fueron las más importantes (Miranda et al., 2021). En ríos andinos de Colombia, cambios en la temperatura del agua y el oxígeno disuelto influyen sobre la composición de la ictiofauna (^{Carvajal-Quintero et al., 2015}).

Si bien en los Andes Tropicales se encuentran las cabeceras de las cuencas con mayor diversidad de peces neotropicales (Amazonas y Orinoco), la ictiofauna de las zonas montañosas y de piedemonte ha sido poco estudiada en comparación con la de tierras bajas (^{Anderson & Maldonado-Ocampo, 2011}). El Garagoa, por ejemplo, es un río andino de alta montaña que hace parte de la cuenca del río Orinoco y brinda servicios ecosistémicos a 33 municipios en los departamentos de Boyacá y Cundinamarca (Colombia), en los que se realizan diferentes actividades económicas, como minería, ganadería y agricultura (^{Corpochivor, 2018}). Las múltiples condiciones generadas por estas actividades, sumadas a sus características naturales, brindan una heterogeneidad de ambientes y condiciones que son ideales para evaluar cómo las variables ambientales a diferentes escalas espaciales influyen sobre las comunidades de peces de sistemas lóticos andinos. A pesar de que para esta cuenca existen estudios sobre la composición taxonómica de peces (^{Barrios-Alonso et al., 2025}; ^{Urbano-Bonilla et al., 2021}), se desconocen los factores que intervienen en la composición y distribución de estos organismos. En consecuencia, el objetivo de este trabajo es identificar la influencia de variables a escala local (oxígeno disuelto, pH, conductividad, temperatura, dureza, morfometría del canal, velocidad del agua, caudal y el índice de calidad de la vegetación de ribera QBR-And) y de paisaje (usos del suelo, características geomorfológicas y bioclimáticas) sobre la diversidad beta de la comunidad de peces de la cuenca del río Garagoa, Colombia, incluyendo los fenómenos de reemplazamiento y anidamiento, y el aporte individual de cada una de las especies y localidades sobre este componente de la diversidad. Teniendo en cuenta que la cuenca está siendo intervenida por múltiples actividades humanas y que el gradiente ambiental natural es heterogéneo, se plantean las siguientes hipótesis: 1) la diversidad beta de la comunidad de peces, así como sus fenómenos de reemplazamiento y anidamiento están afectados por variables en múltiples escalas (^{Diaz-Rojas et al., 2023}; ^{Miranda et al., 2021}) y 2) los gradientes ambientales propiciados por las condiciones naturales y de origen antrópico, actúan como filtros de la composición de especies en cada sitio, por lo que el reemplazamiento es el fenómeno dominante de la diversidad beta (^{Baselga, 2010}).

Materiales y métodos

Área de estudio: La cuenca del río Garagoa se ubica al oriente de la cordillera oriental de Colombia entre los 500 y los 3 450 m de elevación, nace en el páramo de Rabanal, desemboca en el río Guavio y forma parte de la macrocuenca del río Orinoco, gradiente en el cual se encuentran ecosistemas de piedemonte llanero, bosque andino, bosque altoandino y páramo, los cuales se encuentran transformadas en cerca del 80 % de su área total (^{Corpochivor, 2018}). Está conformada por 10 subcuencas pertenecientes a 33 municipios de los departamentos de Boyacá y Cundinamarca, extendiéndose a lo largo de 2 508 km² y presenta un régimen de lluvias monomodal, con picos máximos en mayo y junio y precipitaciones mínimas durante diciembre y enero (^{Corpochivor, 2018}).

Muestreo de ictiofauna: Se establecieron 36 sitios de muestreo en toda la extensión de la cuenca, buscando abarcar el gradiente ambiental. El sitio de mayor elevación (S1) se ubicó a 2 887 m y el de menor (S36) a 500 m (Material Suplementario - MSF1). En cada sitio se delimitó un transecto de 100 m de extensión para el muestreo de ictiofauna y el registro de variables ambientales locales. El muestreo se realizó entre febrero y julio de 2021, empleando dos métodos de pesca: i) electropesca con un equipo Smith Root LR-24, realizando dos recorridos en contracorriente en forma de zigzag, cada uno con un esfuerzo de 1 600 seg. aproximados y potencias eléctricas que variaron entre los 500 y 990 voltios dependiendo de la conductividad del agua; y ii) red de arrastre (2 x 1.25 m, abertura de malla 5 mm) durante 50 min en cada transecto. Los individuos capturados fueron anestesiados utilizando eugenol y posteriormente fijados con formaldehido al 5 % de concentración. Pasadas 72 horas se conservaron en alcohol al 70 %. Finalmente se determinaron taxonómicamente a partir de claves taxonómicas especializadas (^{Urbano-Bonilla et al., 2021}; ^{van der Sleen & Albert, 2017}) y se confirmaron las especies por medio de una visita a la colección de peces del Instituto Alexander von Humboldt. El material recolectado se depositó en la colección del Museo de Historia Natural Luis Gonzalo Andrade de la Universidad Pedagógica y Tecnológica de Colombia (UPTC-Pe-00167 al UPTC-Pe-00209).

Variables predictoras: Las variables predictoras fueron agrupadas en cuatro categorías (MST1): i) variables locales (Lo), las cuales incluyen características evaluadas directamente en cada sitio de muestreo; ii) usos del suelo (Us); iii) factores naturales de la cuenca (Cu); y iv) variables bioclimáticas (Bio). La obtención de cada grupo de variables se describe a continuación.

Variables locales (Lo): En cada sitio se registraron características fisicoquímicas del agua con una sonda multiparamétrica Hanna HI 98 194 (porcentaje de oxígeno disuelto, pH, conductividad y temperatura) y un kit de análisis de aguas Aquamerck, con el que se estimó la dureza del agua. Se registró información morfométrica del canal dividiendo cada transecto en cinco segmentos de 20 m para crear un perfil transversal basado en la medición del ancho del río y la profundidad de la sección. Esta última se registró mediante vadeo en cinco ubicaciones distribuidas entre ambas orillas del río. En cada segmento se registró la velocidad superficial de la corriente midiendo el tiempo de recorrido de un objeto flotante durante una distancia determinada. La velocidad, oxígeno disuelto, pH, conductividad y temperatura, fueron tomadas con tres repeticiones en cada transecto, para los análisis posteriores se utilizaron los valores de la media y la desviación estándar (SD) de cada una de estas variables, esta última como medida de la heterogeneidad del hábitat (^{Yang et al., 2015}). Se caracterizó cada transecto por medio de estimación visual de la proporción de gramíneas, vegetación del ecotono, vegetación acuática, vegetación sumergida, la presencia de troncos y hojarasca y otros materiales dentro del canal, porcentaje de sustratos blandos (arena (diámetro 0.006-0.2 mm) y limo (diámetro < 0.006 mm) y duros (bloques (diámetro > 300 mm), cantos (diámetro 60-300 mm), guijarros (diámetro 20-60 mm) y grava (diámetro 0.2-20 mm)), y porcentaje de los mesohábitats como pozos, corriente y rápidos según ^{Motta-Díaz & Vimos-Lojano (2020)}. Adicionalmente se aplicó el protocolo de evaluación de la calidad ecológica de ríos andinos (CERA), el cual utiliza el índice QBR-And (índice de calidad de la vegetación de ribera Andina), en el que se tuvo en cuenta el grado de cubierta de la zona de ribera, la estructura de la cubierta, la calidad de la cubierta y el grado de naturalidad del canal fluvial (^{Acosta et al., 2009}).

Usos del suelo (Us): Las variables de uso del suelo se obtuvieron para cada área de drenaje calculada a partir de cada sitio de muestreo. Para delimitar cada área de drenaje primero se obtuvo un modelo digital de elevación (MDE) de la cuenca del río Garagoa a partir de imágenes satelitales ALOS PALSAR de 12.5 m de resolución (^{Alaska Satellite Facility, 2024}), posteriormente se corrigió mediante la función r.fill.dir con el fin de evitar sesgo en la dirección del flujo debido a depresiones pronunciadas. A continuación se generaron capas de dirección de flujo y drenaje mediante la función r.watershed con una escala mínima de 1 000 unidades, y por último, se delimitaron microcuencas o las áreas de drenaje para cada sitio de muestreo mediante la función r.water.outlet. Cada microcuenca fue aislada en una capa vectorial para finalmente obtener la proporción de cada uso del suelo a partir de una capa de coberturas para la cuenca del río Garagoa (^{Corpochivor, 2018}) (MSF2). Los usos del suelo obtenidos fueron la proporción de áreas destinadas a ganadería, cultivos de papa, áreas urbanas, embalses, pastizales, bosques, arbustales, otros usos antrópicos, otros cuerpos de agua y áreas degradadas.

Factores naturales de la cuenca (Cu): Se consideraron características geomorfológicas, el área, el perímetro y la elevación como factores naturales de la cuenca. La proporción de las características geomorfológicas (fondos de artesa, laderas de artesa, cimas y laderas, laderas estructurales, laderas erosionales, escarpes, planos estructurales, planos de glacis, planos de terraza y vegas) se obtuvo para cada área de drenaje a partir de una capa vectorial disponible para la cuenca del río Garagoa (Corpochivor, 2018). El área y el perímetro de cada área de drenaje se obtuvo a partir de la delimitación con base en el MDE y la elevación en cada sitio de muestreo a partir de un receptor GPS (Garmin Gps Etrex Venture Hc).

Variables bioclimáticas (Bio): Para cada sitio se obtuvo información de 19 parámetros bioclimáticos utilizando el conjunto de datos 2.1 de WordClim (1970-2000), que incluyen datos anuales, trimestrales y estacionales de las variaciones de precipitación, humedad y temperatura ambiental y que se encuentran disponibles de forma libre en la plataforma (^{Anderson-Teixeira et al., 2015}; ^{Fick & Hijmans, 2017}; ^{Hijmans et al., 2005}).

Análisis de datos: Con el fin de reducir el número de variables (69 en total) y extraer los gradientes ambientales principales se realizó un análisis de componentes principales (PCA) para cada una de las cuatro categorías de variables predictoras. Los componentes principales más importantes de cada PCA fueron seleccionados con base en el procedimiento de aleatorización de ter Braak (Rnd-Lambda) con el paquete PCDimension (^{Coombes & Wang, 2022}) y los puntajes de estos componentes seleccionados fueron extraídos para ser utilizados como variables predictoras en análisis posteriores (^{Machado-de Oliveira et al., 2023}). Las variables correlacionadas con los componentes seleccionados se identificaron con la función dimdesc del paquete FactoMineR (^{Lê et al., 2008}), con un valor p = 0.001 como criterio.

Para evaluar la variación en la diversidad beta entre los sitios de muestreo se empleó el índice de Sorensen particionado en sus fenómenos de reemplazamiento y anidamiento, evaluando de esta manera la sustitución y pérdida de especies a lo largo de la cuenca (^{Baselga, 2010}). Esto se calculó empleando la función beta.pair del paquete betapart (^{Baselga et al., 2023}). Para estimar la influencia de las variables predictoras (scores seleccionados en cada PCA) sobre la diversidad beta global y sus fenómenos de anidamiento y remplazamiento se utilizó el análisis de redundancia basada en distancia dbRDA (^{Legendre & Anderson, 1999}), el cual permite realizar el análisis de ordenación utilizando distancias diferentes a la euclidiana, en este caso se emplearon las matrices de disimilitud obtenidas con el índice de Sorensen. Se realizaron pruebas de permutación (n = 10 000, p = 0.05) para cada dbRDA con el fin de determinar la significancia de cada análisis, de cada eje y de cada variable predictora. Los análisis dbRDA se realizaron con la función dbdra utilizando la corrección de Cailliez para autovalores negativos y las pruebas de permutación con la función anova.cca, las dos funciones pertenecientes al paquete vegan (^{Oksanen et al., 2022}). Todos los análisis se realizaron en la plataforma R versión 4.2.3 (^{R Core Team, 2023}), mediante la interfaz RStudio 2023.03.1(^{Rstudio team, 2023}).

Finalmente, se evaluó el aporte de las especies a la diversidad beta a partir de la estimación de la variación de especies individuales en el área de estudio mediante el índice SCBD (Species Contribution to Beta Diversity) y la contribución de las comunidades locales, o singularidad de los sitios y su aporte a la diversidad beta, por medio del índice LCBD (Local Contribution to Beta Diversity), ambos a partir de datos de abundancia previamente tratados con la transformación de Hellinger (^{Legendre & De Cáceres, 2013}) y con la función beta.div del paquete vegan (^{Oksanen et al., 2022}).

Resultados

Diversidad de peces: Capturamos 3 472 individuos pertenecientes a 22 especies, 10 familias y cinco órdenes (MST2). El orden Siluriformes fue el más diverso con 10 especies y seis familias, seguido por Characiformes con ocho especies y dos familias. Los cuatro órdenes restantes presentaron una especie cada uno. La mayor riqueza de especies se registró en los sitios con elevaciones por debajo de 1 000 m, con una riqueza máxima de 13 especies a 523 m.

Registramos una diversidad beta global alta (βsor = 0.93), indicando que los sitios de muestreo presentan una elevada disimilitud en cuanto a su composición de especies, resultado dado principalmente por la contribución del reemplazamiento de especies (βsim = 0.81) (MSF3), y en menor medida por el anidamiento (βnes = 0.12) (MSF4). Los valores máximos de disimilitud se registraron en sitios de mayor elevación, entre los 2 013 y 2 887 m, lo que sugiere una composición de especies particular en la parte alta en relación con el resto de la cuenca (MSF5).

Variables predictoras: De cada conjunto de variables se seleccionaron los siguientes componentes. Para las variables locales, el PC1 (Lo1), PC2 (Lo2) y PC3 (Lo3) que acumularon el 48.7 % de varianza. El Lo1 (24.23 %) diferencia los sitios con aguas de alta variación de la velocidad y alto porcentaje de gramíneas de sitios con valores altos de dureza, ancho, temperatura del agua, alta calidad de la vegetación de ribera y alto porcentaje de rápidos y sustratos duros (MSF6, MST3). El Lo2 (14.12 %) diferencia sitios con altos niveles de conductividad, dureza y sustratos blandos, de sitios con alto porcentaje de vegetación del ecotono y sustratos duros (MSF7, MST3). El Lo3 (10.34 %) diferencia sitios con mayor ancho del cauce y valores altos para el caudal, de sitios con altos niveles de oxígeno disuelto (MSF7, MST3). En cuanto a los usos del suelo, el PC1 (Us1) y PC2 (Us2), que acumulan el 57 % de varianza. El Us1 (37.09 %) separa sitios con mayores porcentajes de áreas naturales de los sitios con mayor proporción de intervención y el Us2 (20.86 %) separa sitios con alta proporción de áreas ganaderas de sitios con uso minero, agrícola y cultivos de papa (MSF8, MST4).

Para las condiciones naturales de la cuenca, los PC1(Cu1), PC2 (Cu2), PC3 (Cu3) y PC4 (Cu4), con el 81 % de varianza acumulada. El Cu1 (37.09 %) separa sitios con área de drenaje y perímetro amplios y un alto porcentaje de vegas y planos de glacis, de áreas con una proporción alta de fondos de artesa y mayor elevación (MSF9, MST5); el Cu2(18.85 %) separa sitios con laderas estructurales y erosionales de áreas escarpadas (MSF9, MST5); el Cu3 (17.26 %) separa fondos de artesa, laderas de artesa y laderas erosionales de planos de terraza (MSF10, MST5); y el Cu4 (12.49 %), separa laderas de artesa y planos estructurales de cimas (MSF10, MST5). Por último y en relación a las variables bioclimáticas, el Bio1 con el 69.72 % de varianza, el cual diferencia sitios con valores altos de las variables bioclimáticas bio 7 y bio 2, que caracterizan el rango anual de temperatura y el rango diurno medio de temperatura, de sitios con valores mayores de bio 1, bio 4, bio 5, bio 6, bio 8, bio 9, bio 10, bio 11, que describen la temperatura media a diferentes escalas y bio 12, bio 13, bio 14, bio 16, bio 17, bio 19, que describen la precipitación a lo largo de los meses y distintas temporadas del año (MSF11, MST6).

Relación de las variables predictoras y la diversidad beta: La diversidad beta global está influenciada por una combinación de variables locales, usos del suelo, características naturales de la cuenca y variables bioclimáticas (Fig. 1A, Tabla 1), que en conjunto explican el 37 % de su variación total (R² ajustado = 0.37, p = 0.0004). Los dos primeros ejes acumularon una varianza del 77 %, el primer eje dbRDA1 del 49 % (p = 0.0001) y el segundo eje dbRDA2 del 27 % (p = 0.0039). En cuanto a las variables predictoras, las variables locales (Lo1, p = 0.00009; Lo3, p = 0.03580), de usos del suelo (Us1, p = 0.00020) y las naturales de la cuenca (Cu1, p = 0.01340; Cu2, p = 0.00310) influyeron a la diversidad beta total. Otros gradientes de las variables locales y bioclimáticas fueron marginalmente significativos (Lo2, p = 0.06599; Bio1, p = 0.07479). Estos resultados indican que las variables locales que influyeron a la diversidad beta fueron la temperatura, pH, dureza, caudal, la variabilidad de la velocidad de la corriente, concentración de oxígeno disuelto, presencia de gramíneas, presencia de materiales en el cauce, dureza del sustrato, ancho del canal y calidad de la vegetación de la ribera (índice QBR). Las de uso del suelo fueron la proporción de bosques, herbazales, actividades mineras, ganaderas, cultivos de papa, otras actividades agrícolas, procesos de urbanización y otras actividades antrópicas. Las características naturales de la cuenca más influyentes fueron la proporción de laderas estructurales, laderas erosionales, cimas y escarpes y en menor medida, la proporción de fondos de artesa, planos de glacis, vegas, áreas y perímetros. Y por último y de forma marginal, las características bioclimáticas importantes fueron el rango anual y medio de temperatura, es decir la variación entre temperaturas altas y bajas, y la temperatura y precipitación media evaluadas en diferentes escalas de tiempo a lo largo del año.

Fig. 1 Biplots del análisis de redundancias basado en distancias (dbRDA). A. La diversidad beta global. B. El reemplazamiento de especies de peces. En cada uno se incluyen los sitios de muestreo (S) y los componentes principales (PC) seleccionados para las cuatro categorías de las variables predictoras: variables a escala local (Lo1, Lo2, Lo3), variables de usos del suelo (Us1, Us2), variables naturales de la cuenca (Cu1, Cu2, Cu3 y Cu4) y variables bioclimáticas (Bio1). / Fig. 1. Biplots of the distance-based redundancy analysis (dbRDA). A. Overall beta diversity. B. Fish species turnover. Each plot includes the sampling sites (S) and the principal components (PC) selected for the four categories of predictor variables: local-scale variables (Lo1, Lo2, Lo3), land-use variables (Us1, Us2), natural watershed variables (Cu1, Cu2, Cu3, and Cu4), and bioclimatic variables (Bio1). 

De manera similar a la diversidad beta global, el reemplazamiento de especies fue influenciado por una combinación de condiciones locales, de usos del suelo, características naturales de la cuenca y marginalmente por las condiciones bioclimáticas (Fig. 1B, Tabla 1), estas variables explicaron el 7 % de la variación de los datos (R² ajustado = 0.07, p = 0.0001). El primer eje dbRDA1 (p = 0.0001) acumuló el 28 % de la varianza y fue el único significativo. Para la escala local, el Lo1 (p = 0.00009) mostró una influencia significativa sobre el reemplazamiento de especies y el Lo3 de forma marginal (p = 0.09349), con las variables pH, caudal, dureza del agua, dureza del sustrato, temperatura, oxígeno disuelto, variabilidad de la velocidad de la corriente, calidad de la vegetación de la ribera (índice QBR), presencia de gramíneas y materiales en el cauce. En la escala de paisaje fue significativo el Us1 (p = 0.00270) de los usos del suelo, que incluye áreas de bosque, áreas con herbazales, zonas mineras, ganaderas, sitios antropizados y procesos de urbanización. El Cu2 (p = 0.01660) de las variables naturales de la cuenca con las variables, proporción de laderas estructurales, laderas erosionales, cimas y escarpes. En cuanto a las variables bioclimáticas, ningún componente fue significativo sobre el fenómeno de reemplazamiento. Por otra parte, las variables predictoras no influenciaron el fenómeno de anidamiento (p = 0.1053), por lo tanto, el dbRDA correspondiente no fue interpretado.

En relación con el aporte de especies individuales a la diversidad beta, se obtuvieron valores de SCBD entre 0.0004 y 0.2638 (media = 0.0455; Tabla 2). Las especies que presentaron valores de importancia superiores a la media fueron Dolichancistrus fuesslii, Knodus deuterodonoides, Onchorhynchus mykiss, Chaetostoma dorsale, Astroblephus latidens y Grundulus bogotensis. O. mykiss y G.bogotensis son especies que se restringen a la parte alta de la cuenca, siendo en varios sitios las dos únicas especies presentes, por tanto son las que diferencian esta zona de la cuenca. Mientras que D. fuesslii, K. deuterodonoides, C. dorsale y A. latidens son las especies con distribución más amplia a lo largo de la cuenca. En cuanto al aporte de las localidades, se obtuvieron valores de LCBD entre 0.0115 y 0.0497 (media = 0.0278), donde los sitios de la parte de mayor elevación (S1, S2, S3, S4, S5 y S6), sitios de la parte baja (S30, S32, S33, S34, S35 y S36) y algunos de elevación intermedia (S12, S13, S24 y S25) son los conjuntos de sitios que presentaron una composición de especies singular (MST7).

Discusión

En este trabajo evaluamos la influencia de variables ambientales a escala local y de paisaje sobre la diversidad beta y los fenómenos de remplazamiento y anidamiento de las comunidades de peces en la cuenca del río Garagoa, Colombia. Registramos una alta diversidad beta en la cuenca, que alcanza sus valores más altos en la parte de mayor elevación, indicando altas diferencias entre microcuencas cercanas y una composición de especies particular en zonas de mayor altitud. En concordancia con nuestra primera hipótesis y con lo ya documentado en otros ríos neotropicales, inclusive andinos, la disimilitud responde principalmente al fenómeno de reemplazamiento, patrón que se atribuye principalmente a la clasificación ambiental generada por los diferentes gradientes (^{Carvalho et al., 2021}; ^{Diaz-Rojas et al., 2023}; ^{López-Delgado et al., 2020}; ^{Valencia-Rodríguez et al., 2023}). Además, el aumento del reemplazamiento a mayores altitudes en sistemas andinos ya ha sido documentado en otros estudios (^{Anderson & Maldonado-Ocampo, 2011}) y grupos taxonómicos como los macroinvertebrados (^{Diaz-Rojas et al., 2023}).

Tanto la diversidad beta global como el reemplazamiento de especies fueron influenciados por la acción conjunta de variables locales y de paisaje, apoyando nuestra segunda hipótesis según la cual, la diversidad beta es afectada por variables de ambas escalas. A escala local, las variables influyen en la diversidad beta siguiendo un gradiente altitudinal, en el que los sitios de mayores elevaciones (2 106 y 2 887 m), que presentan mayor disimilitud y reemplazamiento de especies, se caracterizan por presentar bajos caudales, alta variabilidad de la velocidad de la corriente, mayor concentración de oxígeno disuelto, alta presencia de gramíneas, alto contenido de materiales en el cauce y dominancia de sustratos blandos, valores que disminuyen con la altitud. Las zonas medias (861 y 1 960 m) presentan los mayores caudales y bajo porcentaje de oxígeno disuelto; mientras que las zonas bajas (500 a 711 m), bajos caudales, altos niveles de dureza, pH y temperatura, alta calidad de la vegetación de la ribera, dominancia de sustratos duros y tramos anchos. Estas variables afectan la distribución de la composición de especies por medio de la biología propia de cada una y sus niveles de tolerancia a diferentes condiciones del entorno. Por ejemplo, la temperatura y el oxígeno disuelto influyen directamente en la sobrevivencia de las especies, pues alteran sus procesos metabólicos, fisiológicos y reproductivos, afectando de manera importante la distribución altitudinal de las especies (^{Carvajal-Quintero et al., 2015}; ^{Jaramillo-Villa et al., 2010}). En el mismo sentido, la conductividad del agua condiciona la osmorregulación de los organismos, por lo que niveles altos de conductividad solo pueden ser soportados por especies tolerantes (^{Zhang et al., 2019}). Cambios en el pH del agua pueden optimizar el crecimiento y locomoción de los peces, en la tasa metabólica de los peces (^{Konstantinov et al., 2008}), sin embargo, valores extremos pueden ser tóxicos al interrumpir la regulación iónica de las branquias e incrementar la toxicidad de varios contaminantes (^{European Inland Fisheries Advisory Commission, 1969}; ^{Freda & Mcdonald, 1988}). Los caudales altos también intervienen al seleccionar especies según sus características morfológicas y comportamentales, ya que estas permiten que algunas especies mantengan la posición en corrientes rápidas, mientras que otras serán arrastradas aguas abajo (^{Winemiller et al., 2008}). La vegetación del ecotono, calidad de la cobertura vegetal y presencia de gramíneas filtran sedimentos, carga excesiva de nutrientes y pesticidas, contribuyendo a mantener la calidad del hábitat (^{Anbumozhi et al., 2005}) y, por lo tanto, impactando en el éxito reproductivo, la sobrevivencia de huevos y larvas, la evasión de depredadores, disponibilidad de hábitat y de alimento para la ictiofauna (^{Pusey & Arthington, 2003}).

En la escala de paisaje la diversidad beta global y el reemplazamiento son afectados por la distribución de las áreas naturales e intervenidas por el ser humano. En las zonas altas y bajas de la cuenca, por ejemplo, se presentan áreas con mayor proporción de bosques y herbazales; mientras que en la zona media hay una mayor proporción de áreas antropizadas, que incluyen actividades mineras, ganaderas, procesos de urbanización y otras actividades. Por lo tanto, las actividades humanas evaluadas desde una perspectiva de cuenca hidrográfica influyen claramente en la composición de la ictiofauna. Otros autores han registrado relaciones negativas entre las áreas de pastizales y cultivos sobre la presencia de especies que muestran baja tolerancia a estas perturbaciones, y relaciones positivas de especies tolerantes con los pastizales, la erosión del suelo y la huella humana (^{Miranda et al., 2021}). Asimismo, se ha documentado el impacto negativo de actividades mineras, usos agrícolas, urbanización y pérdida de bosques en la diversidad de peces (^{Azevedo-Santos et al., 2021}; ^{Diamond & Serveiss, 2001}; ^{Wright & Flecker, 2004}; ^{Zeni et al., 2019}). Estas actividades inducen cambios severos en la hidrología, geomorfología y conectividad fluvial, causando pérdidas significativas de hábitats ribereños y acuáticos, modificando los flujos de carbono, favoreciendo la presencia de especies invasoras y sumando contaminantes urbanos y agrícolas (^{Pelicice et al., 2021}).

Dentro de esta misma escala de paisaje, características naturales de la cuenca como la proporción de laderas estructurales, laderas erosionales, cimas y escarpes, y en menor medida la proporción de fondos de artesa, planos de glacis, vegas, área y perímetro, influyen en la diversidad beta global y reemplazamiento. Resultados semejantes a lo documentado en otros ríos de zonas templadas y subtropicales, donde se encontró que la diversidad de peces está influenciada por los gradientes geomorfológicos, sobre todo con la pendiente, la elevación y con el área de los arroyos (^{Rinne & Miller, 2006}; ^{Shields et al., 2021}; ^{Walters et al., 2003}), siendo en ocasiones variables más importantes que las características del hábitat, como la vegetación acuática o condiciones fisicoquímicas (^{Dauwalter et al., 2008}). Finalmente, las variables bioclimáticas mostraron una incidencia marginal sobre la disimilitud de especies, pero posiblemente se debe a que en las zonas altas hay un mayor contraste entre los cambios de temperatura tanto diarios como a lo largo del año, así como entre periodos secos y periodos lluviosos; mientras que en las zonas bajas hay mayores precipitaciones y mayor temperatura durante todo el año, generando condiciones más estables en el tiempo que se reflejan en la diversidad beta de la comunidad en esta cuenca andina.

En relación con el aporte de los sitios a la diversidad beta, se observa que hay sitios con composiciones particulares de especies en tres grupos de sitios, en la parte alta, la parte media y la parte baja del río, lo que concuerda con la separación de sitios según sus condiciones ambientales. Por su parte, las especies con mayores aportes a la diversidad beta fueron A. latidens, C. dorsale, D. fuesslii, G. bogotensis, K. deuterodonoides y O. mykiss. D. fuesslii presenta una amplia distribución altitudinal, sin embargo, sus mayores abundancias se concentran entre 861 y 1 732 m, K. deuterodonoides, se encuentra de la parte media hacia la parte baja, entre 500 y 1 432 m, siendo más abundante entre 1 275 y 1 340 m. A. latidens se encuentra entre 500 y 1 732 m, presentando algunos registros a 2 112 m, C. dorsale se concentra en la parte baja entre los 523 y los 1 022 m. La parte alta del río se caracteriza por la presencia de O. mykiss y G. bogotensis, la primera encontrándose entre los 2 046 y los 2 887 m y la segunda que únicamente se registra en tres localidades entre los 2 059 y 2 855 m. En este sentido, otros autores han documentado una menor riqueza y un aumento de endemismo a medida que incrementa la elevación (^{Anderson & Maldonado-Ocampo, 2011}), sin embargo, la dominancia de O. mykiss en las zonas altas de la cuenca es preocupante, considerando que es una especie invasora capaz de disminuir las poblaciones, incluso desplazar especies nativas y endémicas como G. bogotensis (^{Alexiades et al., 2022}).

En conclusión, los resultados indican que la diversidad beta de peces de sistemas lóticos andinos es alta y está dada principalmente por el aporte del reemplazamiento de especies. Además, se ve afectada por una combinación de variables de diferentes escalas, entre las que se encuentran los gradientes ambientales naturales, así como las alteraciones de origen antrópico. Dentro de estas variables se incluyen características fisicoquímicas del agua, morfología del canal, calidad de la vegetación de la ribera, características del relieve, así como la distribución de áreas con coberturas naturales y áreas intervenidas por procesos agrícolas, ganaderos, mineros o urbanos. Esto demuestra la importancia de incluir de manera simultánea información ambiental en múltiples escalas en los procesos de monitoreo y manejo de las cuencas andinas, complementando así las estrategias existentes para la conservación de estos ecosistemas acuáticos y su biodiversidad. La evaluación de variables a escala local y del paisaje es importante para predecir patrones en la composición de peces en paisajes influenciados por el ser humano (^{Sály et al., 2011}), así como la inclusión de aspectos geomorfológicos, que son el resultado de las características propias de la formación natural del suelo sumadas a las interacciones con diferentes organismos vivos y a las actividades humanas (^{Reinhardt et al., 2010}) y que proporcionan un hábitat para la biota y un marco físico para los procesos de los ecosistemas (^{Elosegi et al., 2010}).

Por ejemplo, en el plan de ordenación y manejo de la cuenca se realiza una evaluación de parámetros ambientales locales y de paisaje (^{Corpochivor, 2018}) pero, en nuestra opinión, no se analiza la implicación que tiene cada una de estas variables sobre los grupos taxonómicos. Además, en dicho plan se proponen estrategias para el manejo de descarga de aguas contaminadas, control de usos del suelo por medio de licencias ambientales y recuperación de ecosistemas, las cuales están enfocadas de forma prioritaria en áreas estratégicas como son los bosques, humedales, páramos y áreas protegidas (^{Corpochivor, 2018}). Si bien este enfoque puede mejorar la salud del ecosistema, consideramos que la gran variación en la composición de especies entre los sitios de muestreo de esta cuenca, observada también en otros grupos taxonómicos (^{Becerra-Infante et al., 2024}; ^{Diaz-Rojas et al., 2023}), podría generar que algunas especies y comunidades acuáticas no estén protegidas por los esquemas y/o áreas priorizadas. Por lo tanto, se requiere un abordaje en múltiples escalas (^{Erős et al., 2018}), que incluya el mantenimiento de la heterogeneidad del hábitat y la conectividad a escala de cuencas fluviales enteras (^{López-Delgado et al., 2020}). Por último, una porción importante de la variación en la diversidad beta no fue explicada por las variables analizadas, pudiendo responder a factores bióticos, como depredación y competencia (^{Jackson et al., 2001}; ^{Tesfay-Gebrekiros, 2016}), lo que invita a realizar más estudios en esta y otras cuencas andinas que permitan seguir comprendiendo la interacción de los peces con su entorno y cómo se refleja en los patrones comunitarios.

Declaración de ética: Los autores declaran que todos están de acuerdo con esta publicación y que han hecho aportes que justifican su autoría; que no hay conflicto de interés de ningún tipo; y que han cumplido con todos los requisitos y procedimientos éticos y legales pertinentes. Todas las fuentes de financiamiento se detallan plena y claramente en la sección de agradecimientos. El respectivo documento legal firmado se encuentra en los archivos de la revista.

Material suplementario

El material suplementario se encuentra en

disponible en el repositorio de acceso libre www.zenodo.org.

Ver material suplementario

a32v73n1-suppl1