UNIVERSIDAD DE COSTA RICA SISTEMA DE ESTUDIOS DE POSGRADO ANÁLISIS METAGENÓMICO DE FUENTES TERMALES DE COSTA RICA: POTENCIAL FUENTE DE ENZIMAS HIDROLÍTICAS TERMOESTABLES. Tesis sometida a la consideración de la Comisión del Programa de Estudios de Posgrado en Ciencias Biomédicas para optar al grado y título de Maestría Académica en Bioinformática y Biología de Sistemas LAURA BRENES GUILLÉN Ciudad Universitaria Rodrigo Facio, Costa Rica 2021 ii Dedicatoria A Enrique, mi papá, mi mamá, mi hermana y a los perris que me acompañaron todo este tiempo. Agradecimientos A Dios, mi familia, compañeros del CIBCM, a mi tutora Lorena, y lectores César y Rebeca por sus comentarios, sugerencias y recomendaciones. A Daniela, Kaylen y Daniel por su ayuda. A todos los profesores que de una u otra forma estuvieron en el camino y contribuyeron a mi formación profesional. "Esta Tesis fue aceptada por la Comisi6n del Programa de Estudios de Posgrado en Ciencias Biomedicas de la Universidad de Costa Rica, come requisite parcial para optar al grado y tftulo de Maestrfa Academica en Bioinformatica y Biologfa de Sistemas" M.S.c. Marielos Mora Lopez Representante del Decana Sistema de Estudios de Posgrado M.Sc. Lorena Uribe Lorfo Profesora Guia PhD. Rebeca Campos Sanchez PhD. Cesar Rodrfguez Sanchez Lectora ~ er Vaglio Cedeiio Representante Programa de Posgrado en Ciencias Biomedicas Laura Brenes Guillen Sustentante ... lll iv Tabla de contenido Dedicatoria .............................................................................................................. ii Agradecimientos ...................................................................................................... ii Hoja de aprobación ................................................................................................ iii Resumen ................................................................................................................ vi Lista de cuadros .................................................................................................... vii Lista de figuras ..................................................................................................... viii Cuadro de abreviaciones.........................................................................................x Introducción ............................................................................................................ 1 Antecedentes ......................................................................................................... 6 Caracterización de comunidades de bacterias y arqueas en fuentes termales... 6 Metagenómica shotgun para el estudio de comunidades microbianas ............... 7 Ensamblaje y anotación de (meta)genomas ..................................................... 11 Enzimas hidrolíticas termoestables ................................................................... 16 Peptidasas ..................................................................................................... 18 Lipasas .......................................................................................................... 20 Amilasas ........................................................................................................ 22 Enzimas involucradas en la remoción de metales pesados y agroquímicos . 23 Justificación .......................................................................................................... 26 Objetivos .............................................................................................................. 27 Objetivo General ............................................................................................... 27 Objetivos Específicos ........................................................................................ 27 Materiales y Métodos ........................................................................................... 28 Características de las muestras ........................................................................ 28 Análisis Bioinformático ...................................................................................... 35 Análisis de calidad de los datos de secuenciación ........................................ 37 Ensamblaje, alineamiento y anotación de metagenomas .............................. 37 Perfil taxonómico y funcional de las comunidades microbianas .................... 38 Genomas obtenidos a partir de (meta)genomas (MAGs) .............................. 41 Creación de base de datos ............................................................................ 41 Mapeo y abundancia de variantes alélicas .................................................... 42 Mapeo de genes de resistencia a metales pesados y agroquímicos en tapetes microbianos de fuentes termales. .................................................................. 43 Resultados y discusión ......................................................................................... 44 v Análisis de calidad de las secuencias ............................................................... 44 Ensamblaje de los metagenomas ..................................................................... 46 Genomas obtenidos a partir de (meta)genomas ............................................... 49 Comparación del perfil taxonómico de las diferentes fuentes termales ............ 55 Fuentes termales de Costa Rica .................................................................... 58 Porcelana, Chile ............................................................................................ 59 Mound y Bijah, PNY ....................................................................................... 60 Obsidian Pool, PNY ....................................................................................... 61 Liberty Cap, Mammoth springs, PNY. ............................................................ 62 Isla Raoult, Nueva Zelanda ............................................................................ 62 Perfil funcional de las comunidades de las fuentes termales de Costa Rica, PNY, Chile y Nueva Zelanda............................................................................. 63 Identificación y clasificación de las enzimas hidrolíticas ................................... 65 Peptidasas ..................................................................................................... 73 Proteasas de serina termofílicas................................................................ 73 Metalopeptidasa (Metaloproteínas) ........................................................... 78 Lipasas .......................................................................................................... 92 Esterasas ....................................................................................................... 97 Amilasas ...................................................................................................... 102 Genes involucrados en la remoción de metales y agroquímicos encontrados en los tapetes microbianos .................................................................................. 105 Arsénico ....................................................................................................... 107 Cobre ........................................................................................................... 111 Plomo........................................................................................................... 117 Hierro ........................................................................................................... 118 Aluminio ....................................................................................................... 120 Paraquat ...................................................................................................... 121 Atrazina ........................................................................................................ 125 Análisis de correlación entre los genes relacionados con la remoción de agroquímicos y los principales filos encontrados en los tapetes microbianos . 126 Conclusiones ...................................................................................................... 130 Recomendaciones .............................................................................................. 132 Referencias ........................................................................................................ 133 Anexos ............................................................................................................... 190 vi Resumen En el presente estudio independiente de cultivo, se investigó a nivel taxonómico y funcional las comunidades microbianas de seis metagenomas de fuentes termales cercanas a la Cordillera de Tilarán, volcán Miravalles, volcán Rincón de la Vieja y zona sur, y datos de 13 metagenomas de fuentes termales del Parque Nacional Yellowstone, Nueva Zelanda y Chile. Se obtuvieron en total 229 Gbp de secuencias cortas obtenidas con la plataforma Illumina. Los resultados de la comparación del perfil taxonómico muestran que los tapetes microbianos de Chile y Costa Rica tienen 60% similitud, junto a las muestras de Nueva Zelanda que provienen de agua filtrada, siendo Proteobacteria, Firmicutes y en algunas muestras Cyanobacteria y Chloroflexi los filos más abundantes. Se creó una base de datos de potenciales enzimas hidrolíticas tales como peptidasas, amilasas, lipasas y esterasas, así como enzimas involucradas en la remoción de metales pesados y agroquímicos. La base de datos tiene más de 48 000 secuencias, muchas de las cuales podrían ser novedosas y objeto de estudios posteriores. La mayoría de las secuencias tienen un porcentaje de identidad entre el 70% y el 90% con las de microorganismos que toleran altas temperaturas, lo que sugiere que podrían tener funciones similares. Muchas de las enzimas encontradas dentro de estas categorías tienen una función housekeeping, por lo que son muy importantes para el mantenimiento de las funciones básicas de las células, sin embargo, podrían participar en los procesos de degradación e hidrólisis de herbicidas y diversos xenobióticos. En los MAGs obtenidos se encontraron genes que participan en el transporte, acumulación y oxidación-reducción de arsénico, cobre, plomo, hierro, aluminio y los herbicidas paraquat y atrazina, estos resultados sugieren la presencia de diversos mecanismos celulares, que les confieren a estos microorganismos resistencia/tolerancia a alta concentraciones de estos metales o por el contrario tienen mecanismos para poder obtenerlos cuando las concentraciones son muy bajas. Este estudio proporcionó un primer acercamiento diversidad funcional de los microorganismos que habitan estos sitios, las secuencias encontradas en este estudio podrían utilizarse para futuras investigaciones en el área de la biorremediación y la geomicrobiología. vii Lista de cuadros Cuadro 1. Localización geográfica, temperatura, pH de muestras utilizadas en este proyecto. ....................................................................................................... 28 Cuadro 2. Estadísticas de los datos crudos de secuenciación de las muestras de las fuentes termales de Costa Rica ...................................................................... 30 Cuadro 3. Características generales de cada muestra utilizada en el análisis. .... 31 Cuadro 4. Categorías funcionales según COG. ................................................... 40 Cuadro 5. Esquema de clasificación numérica para enzimas, basado en las reacciones químicas que catalizan ....................................................................... 42 Cuadro 6. Resultados de análisis de calidad de las muestras antes y después del trimming................................................................................................................ 46 Cuadro 7. Resultados de MetaQUASTt obtenido a partir de los ensamblajes con MEGAHIT ............................................................................................................. 47 Cuadro 8. Resultados de MetaQUAST obtenido a partir de los ensamblajes con MetaSPAdes ........................................................................................................ 48 Cuadro 9. Resumen de los MAGs encontrados en los tapetes microbianos. ....... 51 Cuadro 10. Resultados generales de la anotación de los metagenomas utilizando Prokka. ................................................................................................................. 66 Cuadro 11. Número total de variantes/alelos en cada metagenoma según la clasificación de las enzimas (bajo la categoría EC3). .......................................... 68 Cuadro 12. Resultados de UniProt-Blast para las variantes de metaloproteínas más abundantes de las variantes ftsH.................................................................. 81 Cuadro 13. Resultados de Uniprot-Blast para las variantes de metaloproteínas más abundantes de las variantes de los genes mmpA y rip................................. 85 Cuadro 14. Resultados de UniProt-Blast de las variantes de metaloproteínas más abundantes de los genes pmA. ............................................................................ 87 Cuadro 15. Resultados de UniProt-Blast para las variantes de metaloproteínas . 89 Cuadro 16. Resultados de UniProt-Blast para las variantes de lipasas ............... 96 Cuadro 17. Resultados de Uniprot-Blast para las variantes de esterasas ......... 102 Cuadro 18. Resultados de UniProt-Blast para las variantes de amilasas........... 104 Cuadro 19. Genes que participan en la degradación de atrazina encontrados en las fuentes termales. .......................................................................................... 126 viii Lista de figuras Figura 1. Descripción general del análisis bioinformático para el perfil microbiano del gen 16S ARNr y la metagenómica shotgun _________________________ 10 Figura 2. Clasificación de peptidasas según su especificidad y tipo de catálisis que realizan ________________________________________________________ 19 Figura 3. Distribución de peptidasas encontradas en bacterias y/o arqueas en clanes y familias según MEROPS ___________________________________ 20 Figura 4. Árbol filogenético de las lipasas y esterasas clasificadas en 19 subfamilias ____________________________________________________ 221 Figura 5. Fotografías de fuentes termales costarricenses _________________ 33 Figura 6. Fotografías de fuentes termales de Yellowstone_________________ 34 Figura 7. Diagrama general de la metodología de análisis bioinformático______38 Figura 8. Cantidad de MAGs aislados de tapetes microbianos de seis fuentes termales según el filo al que pertenecen ______________________________ 50 Figura 9. Diagrama de calor (Primer 7) que muestra las relaciones globales entre grupos de muestras ______________________________________________ 58 Figura 10. Ordenación de escala multidimensional no métrica (NMDS) de las comunidades bacterianas a nivel taxonómico de filo de las fuentes termales __ 58 Figura 11. Mapa de calor de la abundancia estandarizada de cada una de las variantes encontradas en los tapetes microbianos _______________________ 71 Figura 12. Abundancia de las enzimas clasificadas como serina proteasa en las fuentes termales. ________________________________________________ 75 Figura 13. Alineamiento de secuencias clasificadas como serina proteasa encontradas en los tapetes microbianos de las fuentes termales costarricenses___________________________________________________ 77 Figura 14. Abundancia de las diferentes variantes probablemente asociadas a metalopeptidasas en los tapetes microbianos de las fuentes termales costarricenses __________________________________________________ 80 Figura 15. Alineamiento de secuencias clasificadas como FtsH encontradas en los tapetes microbianos de las fuentes termales costarricenses _______________ 84 Figura 16. Árbol filogenético de las secuencias de lipasas y esterasas encontradas en los tapetes microbianos. ______________________________ 93 ix Figura 17.G ráfico de sombras de la abundancia relativa de cada uno de los genes involucrados en la tolerancia a metales encontrados en los tapetes microbianos. ___________________________________________________ 107 Figura 18. Representación de los operones arsRDABC y aioXSRBA en los MAGs encontrados en los diferentes tapetes microbianos de las fuentes termales de estudio. _______________________________________________________ 109 Figura 19. Representación de los diferentes genes encontrados en los tapetes microbianos de fuentes termales costarricenses que le confieren a las comunidades microbianas de tapetes microbianos de ambientes termales resistencia al arsénico ___________________________________________ 110 Figura 20. Representación de los diferentes genes encontrados en los tapetes microbianos de fuentes termales costarricenses que le confieren a las comunidades microbianas de tapetes microbianos de ambientes termales resistencia al cobre _____________________________________________ 112 Figura 21. Representación de los operones copBAC y copRMS en los MAGs encontrados en los diferentes tapetes microbianos de las fuentes termales de estudio _______________________________________________________ 115 Figura 22. Representación de los diferentes genes encontrados en los tapetes microbianos de fuentes termales costarricenses que le confieren a las comunidades microbianas de tapetes microbianos de ambientes termales resistencia al plomo _____________________________________________ 118 Figura 23. Representación de los operones chrBACF y adeABC en los MAGs encontrados en los diferentes tapetes microbianos de las fuentes termales de estudio. _______________________________________________________ 124 Figura 24. Resultados de correlación de Spearman entre los genes relacionados con la remoción/degradación de paraquat y atrazina ____________________ 127 x Cuadro de abreviaciones ARNr Ácido Ribonucleico arCOG Archaeal clusters of orthologous genes ATP Adenosín Trifosfato ADN Ácido Desoxirribonucleico CAzy The Carbohydrate-Active EnZymes database COG Clúster de Grupos Ortólogos CDS Región codificante CDD Conserved Domains Database dBg Bruijn graph dNTPs Desoxirribonucleótidos trifosfatos FISH Hibridación in situ de fluorescencia GO Gene Ontology HMM Modelos ocultos de Markov ITS Regiones Espaciadoras Transcritas KEGG Enciclopedia de Kyoto de Genes y Genomas KO KEGG Orthology LCA Last Common Ancestor MAG Genomas obtenidos a partir de metagenomas MDR MultiDrug Resitance MATE Multidrug And Toxin Extrusion SMR Small Multidrug Resistance MPs Metalopeptidasas MGRast Metagenomics Rapid Annotation using Subsystem Technology MFS The major facilitator superfamily NMDS Non-metric multidimensional scaling NGS Secuenciación de Nueva Generación Pb Pares de bases PCR Reacción en cadena de la polimerasa PNY Parque Nacional Yellowstone Pfam protein families database PRMER7 Plymouth Routines in Multivariate Ecological Research RAM Random Access Memory RND Resistance-nodulation-division SBS Sequencing by synthesis SRA Sequence Read Archive SRH Second Region of Homology SIMPROF The similarity profile routine SIMPER Similarity percentage análisis UNIREF UniProt Reference Clusters UniProt Universal Protein database https://www.ncbi.nlm.nih.gov/Structure/cdd/cdd.shtml Autorización para digitalización y comunicación pública de Trabajos Finales de Graduación del Sistema de Estudios de Posgrado en el Repositorio Institucional de la Universidad de Costa Rica. Yo, _______________________________________, con cédula de identidad _____________________, en mi condición de autor del TFG titulado ___________________________________________________ _____________________________________________________________________________________________ _____________________________________________________________________________________________ Autorizo a la Universidad de Costa Rica para digitalizar y hacer divulgación pública de forma gratuita de dicho TFG a través del Repositorio Institucional u otro medio electrónico, para ser puesto a disposición del público según lo que establezca el Sistema de Estudios de Posgrado. SI NO * *En caso de la negativa favor indicar el tiempo de restricción: ________________ año (s). Este Trabajo Final de Graduación será publicado en formato PDF, o en el formato que en el momento se establezca, de tal forma que el acceso al mismo sea libre, con el fin de permitir la consulta e impresión, pero no su modificación. Manifiesto que mi Trabajo Final de Graduación fue debidamente subido al sistema digital Kerwá y su contenido corresponde al documento original que sirvió para la obtención de mi título, y que su información no infringe ni violenta ningún derecho a terceros. El TFG además cuenta con el visto bueno de mi Director (a) de Tesis o Tutor (a) y cumplió con lo establecido en la revisión del Formato por parte del Sistema de Estudios de Posgrado. FIRMA ESTUDIANTE Nota: El presente documento constituye una declaración jurada, cuyos alcances aseguran a la Universidad, que su contenido sea tomado como cierto. Su importancia radica en que permite abreviar procedimientos administrativos, y al mismo tiempo genera una responsabilidad legal para que quien declare contrario a la verdad de lo que manifiesta, puede como consecuencia, enfrentar un proceso penal por delito de perjurio, tipificado en el artículo 318 de nuestro Código Penal. Lo anterior implica que el estudiante se vea forzado a realizar su mayor esfuerzo para que no sólo incluya información veraz en la Licencia de Publicación, sino que también realice diligentemente la gestión de subir el documento correcto en la plataforma digital Kerwá. https://es.wikipedia.org/wiki/Responsabilidad https://es.wikipedia.org/wiki/Perjurio 1 Introducción En el pasado y el presente, la búsqueda de vida requiere de una sólida comprensión del origen y evolución de los organismos en el planeta Tierra. Los procariotas como las arqueas y bacterias son posiblemente los únicos organismos que fueron capaces de tolerar las condiciones ambientales y atmosféricas de la tierra de hace 3.8 mil millones de años, tales como falta de oxígeno, alta concentración de gases (metano, dióxido de carbono, nitrógeno, hidrogeno y monóxido de carbono), alta radiación y elevadas temperaturas (Javaux 2006), y persistir hasta el día de hoy. Los ambientes extremos actuales son una representación de algunas de las condiciones de la Tierra hace millones de años, por lo que estudiar los microorganismos que allí se encuentran resulta de gran interés. Un organismo que es capaz de sobrevivir en un ambiente extremo se considera extremófilo. Los extremófilos pueden tolerar condiciones extremas físicas como la temperatura, la radiación, la presión y concentración de nutrientes, y extremos geoquímicos como por ejemplo la desecación, salinidad, pH, especies de oxígeno o potencial redox (Rothschild y Mancinelli, 2001). Según Kay (1980) un organismo extremófilo debe ser capaz de desarrollar un mecanismo para limitar, reducir o eliminar los efectos ambientales en su metabolismo o simplemente desarrollar estrategias funcionales para “vivir” con ese factor. Los estudios de ambientes extremos, como lo son las fuentes termales, han generado mucho interés en la comunidad científica desde 1960, con las investigaciones en el Parque Nacional Yellowstone realizadas por el microbiólogo Thomas Brock (Brock 2001), y los análisis elementales en este sitio realizados por Allen and Day en 1935 y Schoen en 1969 (Brock y Mosser, 1975), los cuales se consideran pioneros. La caracterización microbiológica, geológica y química de las fuentes termales ha sido objeto de estudio alrededor del mundo. La temperatura es uno de los principales factores ambientales que controlan las funciones de los organismos, reproducción y el crecimiento (Kay, 1980). Las fuentes termales se caracterizan por tener temperaturas elevadas (45ºC-100°C), pueden ser alcalinas 2 o ácidas, y generalmente tienen una alta concentración de iones (Rothschild y Mancinelli, 2001, Djokic et al., 2017). La metagenómica ha permitido realizar estudios de comunidades de ambientes extremos a detalle, no solo a nivel taxonómico sino también a nivel metabólico y funcional (López-López et al., 2013, Cowan et al., 2015, van der Walt et al., 2017). Los estudios para la identificación taxonómica son un paso fundamental para describir la distribución de los miembros de la comunidad, ya que brindan información acerca de la abundancia de los diferentes microorganismos (Jiang y Takacs-Vesbach, 2017, Thiel et al., 2016). La secuenciación de alto rendimiento o también llamada Secuenciación de Nueva Generación (NGS) ha permitido ampliar los estudios de comunidades microbianas mediante la secuenciación de marcadores filogenéticos como los genes para el ARNr 16S y las regiones ITS (regiones espaciadoras transcritas), lo cual permite realizar estudios de estructura. Por otro lado, la metagenómica permite realizar una caracterización de sus potencialidades a nivel funcional (Boteva y Kambourova 2018). Esta herramienta ha generado grandes avances en el estudio de las características metabólicas, fisiológicas y ecológicas de las comunidades termotolerantes (Hiraoka et al., 2016 y Saxena et al., 2017). Las secuencias de los metagenomas ensamblados de las comunidades microbianas han revelado la existencia de genes que le confieren a diferentes organismos la capacidad adaptarse o tolerar condiciones extremas. Más específicamente para esta investigación, se han realizado estudios para la búsqueda y caracterización de genes que codifican para enzimas hidrolíticas termoestables (Saxena et al., 2017, Sahay et al., 2017), las cuales juegan un rol importante dentro del metabolismo de los microorganismos. Además, se han realizado investigaciones sobre la dinámica y vías metabólicas, procesos celulares y biogeoquímicos de los microorganismos que conforman las comunidades bacterianas y de arqueas (Mehetre et al., 2016a, Mehetre et al., 2018). Todas estas investigaciones ofrecen una amplia información sobre posibles mecanismos o herramientas que pueden utilizarse en procesos industriales o biotecnológicos, y que pueden servir para solucionar problemas ambientales y de salud pública como la contaminación (Saxena et al.,2017, Bao et al., 2017). 3 La contaminación ambiental por compuestos tóxicos es una gran amenaza para el medio ambiente y la salud humana a nivel mundial (Singh y Naidu, 2012), por lo tanto, es de gran importancia conocer las vías metabólicas que tienen los microorganismos relacionados con la eliminación de sustancias contaminantes. El estudio e identificación de las vías metabólicas involucradas en la degradación de xenobióticos por parte de microorganismos podría tener un gran impacto en la biorremediación de compuestos contaminantes (Bao et al., 2017). En este sentido, los microorganismos aislados de ambientes extremos como las fuentes termales son buenos candidatos para la remoción de ciertos compuestos contaminantes, ya que la mayoría poseen diferentes estrategias metabólicas para remover y/o tolerar sustancias tóxicas, que resultan importantes para la biorremediación (Segretin et al.,2016). Cabe destacar las investigaciones que involucran bacterias termófilas en la remoción de metales pesados como hierro, manganeso, cromo y plomo (Vargas et al.,1998, Kashefi y Lovley 2000, Masaki et al.,2015). La búsqueda de enzimas con uso potencial en procesos de biorremediación en estos ambientes ha sido discutida previamente en la literatura, debido a la estabilidad de estas, inherente aún en condiciones de temperatura y pH extremos (Demirjian et al., 2001, de Miguel Bouzas et al., 2006), y en presencia de algunos solventes y detergentes (Bhalla et al., 2013). Los avances en el estudio de estas enzimas han sido posibles gracias al aislamiento de microorganismos termófilos de diferentes ambientes y su posterior extracción y purificación (Burrows, 1973, Antranikian et al., 1987, Groboillot, 1994, Bhushan y Hoondal 1998, Bauer et al., 1999, Kohilu et al., 2001). En Costa Rica se han realizado estudios de las comunidades microbianas y de organismos eucariotas de fuentes termales que se caracterizan por tener temperaturas desde los 35 ºC, rangos de pH de 3 a 9, y en su mayoría poseen una alta concentración de iones (Sittenfeld et al., 2002, Sittenfeld et al., 2004, Fisinger et al., 2008, Morales et al., 2008, Caldwell et al., 2010, Hernández-Ascencio, 2012, Hynek et al., 2018, Arce-Rodríguez et al., 2019, Uribe-Lorío et al., 2019). El objetivo principal de estos proyectos ha sido el aislamiento, identificación molecular y análisis de estructura de las comunidades de bacterias, arqueas y eucariotas presentes. Muchas de las fuentes termales estudiadas poseen orígenes 4 geotectónicos distintos debido a su asociación con volcanes y/o fallas tectónicas (Alvarado y Vargas, 2017), aunque también difieren en el grado de mineralización, y en las concentraciones iones como hidrocarbonato, sulfuro, sulfato, cloruro, nitrato, amonio, entre otros (Alvarado y Vargas, 2007, Alvarado y Vargas, 2017, Uribe-Lorío et al., 2019). Según la información consultada en artículos científicos y libros disponibles en internet en el presente año tanto en inglés como en español (utilizando palabras clave como “metagenómica”, “Costa Rica”, “enzimas”, entre otras) no se ha llevado a cabo hasta el momento investigaciones sobre la búsqueda de enzimas hidrolíticas en estos ambientes de nuestro país y en la región centroamericana utilizando secuenciación masiva paralela, por lo que resulta de un estudio novedoso. La identificación de secuencias de enzimas encontradas en ambientes con altas temperaturas es muy importante para generar información en ciencia básica y biotecnología, ya que son recursos que todavía no han sido descritos, y que por ende aumentan nuestro conocimiento de las comunidades microbianas y las características metabólicas que poseen estos organismos. Es por lo que en este proyecto se plantea la utilización de metagenómica para i) caracterizar parcialmente a nivel metabólico y funcional las comunidades de microorganismos de tres fuentes termales de Costa Rica, y ii) buscar en ellas secuencias de enzimas hidrolíticas que están posiblemente relacionadas con el metabolismo, degradación y remoción de xenobióticos como metales pesados y agroquímicos. Dentro de los principales resultados encontrados destacan que las comunidades microbianas de los seis metagenomas de tapetes microbianos de fuentes termales de Costa Rica son similares a otras alrededor del mundo, sin embargo incluyen algunos microorganismos únicos. Se identificaron secuencias completas de proteasas, lipasas, amilasas, esterasas e hidrolasas, categorías que según diversos estudios incluyen enzimas capaces de hidrolizar compuestos como los agroquímicos. La mayoría de las variantes/alelos encontrados son similares (80% -95% identidad) a las de microorganismos que toleran altas temperaturas, por lo que estas enzimas muy probablemente sean termotolerantes, sin embargo, en la mayoría de los casos, el porcentaje de identidad es menor al 95%, lo que sugiere 5 que estas comunidades constituyen un novedoso reservorio de potenciales enzimas de posible aplicación en el largo plazo en la degradación o remoción de xenobióticos como agroquímicos. La búsqueda de genes que posiblemente participan en la remoción de metales pesados en los 84 MAGs (Genomas obtenidos a partir de metagenomas) reconstruídos permitió asociar los mecanismos de remoción tanto de xenobióticos como de metales pesados con los diferentes grupos de microorganismos presentes. 6 Antecedentes Caracterización de comunidades de bacterias y arqueas en fuentes termales Los microorganismos son parte fundamental de la dinámica de los ecosistemas del planeta, sin embargo, a pesar de su importancia se desconocen muchas de sus funciones ecológicas y papeles dentro de una comunidad (Konopka et al., 2015, Saini et al., 2017). Las comunidades microbianas están típicamente compuestas por especies dominantes y algunos grupos taxonómicos poco comunes o desconocidos (Sogin et al., 2006). Estos organismos poco abundantes pueden tener importantes roles dentro de los ciclos biogeoquímicos y los flujos metabólicos del carbono, nitrógeno, azufre y fósforo (Musat et al., 2008, Wrighton et al., 2012). Las comunidades microbianas de ambientes extremos, específicamente los consorcios termofílicos, han sido ampliamente estudiadas y caracterizadas en algunas regiones del mundo (Amin et al., 2017, Saini et al., 2017, Oliverio et al., 2018), como Canadá, Nueva Zelanda, Estados Unidos, Japón, India y Malasia (Saini et al.,2010). Según Dick y Shock (2013) la estructura microbiana de las fuentes termales depende de gradientes de temperatura y pH, así como de interacciones ecológicas, la química de las rocas y las concentraciones de compuestos inorgánicos (Dick y Shock, 2013). Se ha encontrado que las comunidades microbianas de fuentes termales a lo largo de Costa Rica con temperaturas entre los 37 °C y 63 °C se rigen por las variables fisicoquímicas, tales como pH, temperatura, concentración de gases, cationes y aniones (Uribe-Lorío et al., 2019), por lo que la estructura microbiana al igual que en otras fuentes alrededor del mundo, va a depender principalmente de las características ambientales. Dentro de los grupos de microorganismos más abundantes en estos ambientes se encuentran las bacterias y arqueas (Bowen De León., 2013, Huang et al., 2013). Las fuentes termales poseen condiciones muy específicas que pueden favorecer el crecimiento de bacterias de los filos Acidobacteria, Actinobacteria y Proteobacteria, así como de otros filos que se encuentran en baja abundancia. En fuentes con temperaturas menores a los 70 °C, filos como Cianobacteria y Chloroflexi están entre los principales grupos de bacterias autótrofas que se pueden encontrar (Miller 7 y Castenholz, 2000), siendo Chloroflexi, un conjunto de microorganismos fotoheterótrofos capaces de tolerar un rango de temperaturas muy amplio (Ward et al.,2018). Entre las técnicas clásicas que se utilizan para describir y estudiar las comunidades de microorganismos se encuentran los microarreglos y análisis de Hibridación in situ de fluorescencia (FISH), reacción en cadena de la polimerasa (PCR), así como técnicas caídas en desuso como clonaje de ADN en fósmidos y cósmidos (Kublanov et al., 2009, Kumar et al., 2015), y que han sido progresivamente sustituidas por la secuenciación masiva de un fragmento de ADN como lo es el gen que codifica el ARNr 16S, y la secuenciación shotgun. Las técnicas de NGS permiten generar información de la estructura de la comunidad microbiana mediante la amplificación masiva de un marcador filogenético, así como características metabólicas y funcionales de las comunidades microbianas (Sharpton et al., 2014). Según Lozupone y Knight (2006), en un estudio en donde se comparan 202 muestras de diferentes ambientes existe una amplia variación filogenética y de diversidad de especies, siendo los sedimentos y los suelos los más diversos filogenéticamente. En el caso específico de las fuentes termales, la mayoría de las investigaciones señalan que el aumento de la temperatura afecta la diversidad de microorganismos (Tank et al., 2017, Pagaling et al., 2012, Ross et al., 2012), siendo las fuentes con temperaturas mayores a los 70 °C las menos diversas. No obstante, se ha demostrado que sitios calientes como las ventilas térmicas y las aguas calientes terrestres son fuente de alelos nuevos (Zhang et al., 2008, Auguet et al., 2010, Kumar et al., 2017, Campbell et al., 2017). Metagenómica shotgun para el estudio de comunidades microbianas El término metagenómica fue utilizado por primera vez por Jo Handelsman (Handelsman et al., 1998) y se refiere al análisis genómico de ADN extraído directamente a partir de una muestra ambiental. Los enfoques de metagenómica históricos consisten en el clonaje y la secuenciación por el método de Sanger de librerías de plásmidos y fósmidos (Ghai et al., 2010). No obstante, la 8 implementación de estos métodos consume más tiempo y requiere una concentración mayor de ADN, además de que son técnicas económicamente más costosas. La secuenciación de alto rendimiento permite examinar miles de organismos en paralelo y obtener de forma exhaustiva la mayoría de los genes, proporcionando así información de la biodiversidad y las características funcionales de las comunidades microbianas (Ranjan et al., 2016, Tessler et al., 2017). Con los avances NGS se desarrolló una visión mucho más detallada de las comunidades microbianas, utilizando secuenciación de amplicones para la identificación taxonómica de las comunidades y la secuenciación shotgun, la cual permite caracterizar a nivel funcional y metabólico los microorganismos (Figura 1). El término secuenciación shotgun fue acuñado en 1979 por Rodger Staden quien propuso que la secuenciación shotgun, era en la que se utilizaban vectores bacterianos para clonar fragmentos aleatorios de una molécula de ADN larga, Luego, los fragmentos se secuencian en paralelo y las lecturas se ensamblan utilizando sus superposiciones, lo que permite secuenciar genomas más grandes (Staden 1979), posteriormente, propone en 1980 el término contig, para describir los datos obtenidos después del ensamblaje de las lecturas de secuenciación shotgun. En el contexto de la publicación “un contig es un conjunto de lecturas de gel que se relacionan entre sí por superposición de sus secuencias” (Staden, 1980). Estos conceptos aplican desde entonces para diferentes áreas de la biología molecular, tipos de secuenciación y análisis bioinformáticos (Giani et al., 2019). 9 Figura 1. Descripción general del análisis bioinformático para el perfil microbiano del gen 16S ARNr y la metagenómica shotgun. Modificado de Milani et al., 2017. Dentro de las principales plataformas para secuenciar metagenomas se encuentra la de la casa comercial Illumina Inc, la cual ofrece algunas opciones para la preparación de librerías, la secuenciación y el análisis e interpretación de los datos. Esta genera una gran cantidad de información y alta precisión (tasa de error de 0.1- 1%) en comparación con otras plataformas, como Ion Torrent (~ 1 Gb, 1,5%) y Pacific Biosciences (100 Mb,10-15% tasa de error) (Quail et al., 2012, Buermans y den Dunnen, 2014, Laehnemann et al., 2016, Roumpeka et al., 2017, Driscoll et al., 2017). Es importante considerar la química que se utiliza y el equipo, ya que estas características cambian según las combinaciones. La tecnología de Illumina, que es la más utilizada, es capaz de generar secuencias cortas de entre 100 pb, 300 pb, o más, y miles de millones de secuencias para una única muestra (Breitwieser et al., 2017, van der Walt et al., 2017). Esta plataforma utiliza el método de Secuenciación por Síntesis (SBS). En este tipo de técnica se le adiciona adaptadores en los extremos a los fragmentos de ADN por secuenciar, 10 que consisten en secuencias cortas de nucleótidos que permiten la unión por complementariedad a la celda de flujo, y así poder ser enriquecidos mediante la amplificación tipo puente. Los adaptadores contienen: i) secuencias que permiten que la biblioteca se una y genere cluster en el flow cell (secuencias P5 y P7), ii) secuencias de sitios de unión del primer para iniciar la secuenciación (Rd1 SP y Rd2 SP) y iii) index (index 1 y, cuando corresponda,index 2), que son identificadores de muestra que permiten la multiplexación / agrupación de múltiples muestras en una única secuenciación. La formación de estos grupos o clusters permite amplificar la señal generada por nucleótidos fluorescentes incorporados en la cadena creciente, que es detectada por una cámara durante un proceso cíclico (Roumpeka et al., 2017). Los 2'-desoxinucleósido-trifosfatos o dNTPs se incorporan a la cadena de ADN en crecimiento, donde al menos uno de los cuatro nucleótidos (A, C, T, G) tienen una etiqueta fluorescente diferente que sirve para identificar la base y actuar como terminador reversible, esto con el fin de evitar múltiples eventos de extensión. En circunstancias ideales, todas las bases dentro de un grupo se amplificarán en la misma fase, por lo que los nucleótidos deben agregarse en las cadenas de cada clúster al mismo tiempo, de lo contrario la calidad de la señal irá disminuyendo hacia los extremos de las lecturas (Buermans y Den Dunnen. 2014). La precisión en cada fase ayudará a disminuir la cantidad de falsos positivos y falsos negativos. La señal de fluorescencia es detectada por una cámara, lo cual permite identificar los nucleótidos que conforman cada una de las secuencias (Roumpeka et al., 2017). Como resultado final se pueden obtener secuencias single end (un único extremo, hacia una única dirección) o paired end (extremos pareados, ambas direcciones) según el kit utilizado. Con la secuenciación single end el inserto de la biblioteca se lee solo desde el extremo P5. Con la secuenciación paired end, se genera una segunda lectura desde el otro extremo del inserto (extremo P7). La secuenciación paired end incrementa la cobertura de secuenciación, este tipo de estrategia se utiliza generalmente para la secuenciación del genoma o del exoma completo, y en algunos casos de secuenciación del transcriptoma o metagenómica. La secuenciación de un solo extremo es el método de elección para la secuenciación 11 de ARN, por ejemplo, para el análisis diferencial de la expresión génica. Los extremos pareados ofrecen datos de secuencias de alta calidad debido a que se conoce la distancia entre cada lectura, mejora la alineación y el ensamblaje, además de que permiten la detección precisa de variantes estructurales, fusiones de genes e isoformas. Por otro lado, la utilización de secuencias hacia una única dirección permite generar grandes volúmenes de datos de alta calidad, de forma rápida y económica. Por otro lado, la minería de alelos/genes es uno de estos enfoques novedosos que permite identificar alelos nuevos de organismos diferentes, algunas estrategias que se utilizaron en este estudio son el mapeo y el alineamiento. El mapeo consiste en alinear secuencias de nucleótidos cortas (100 pb – 300 pb) contra secuencias de nucleótidos largas (por ejemplo, contigs), esto ayuda a determinar el origen de la secuencia. En el alineamiento se tienen secuencias de nucleótidos o proteínas que se quiere alinear entre sí, lo que permite identificar similitudes y/o diferencias entre estas. Ensamblaje y anotación de (meta)genomas El ensamblaje de los (meta)genomas es un paso importante para caracterizar la ecología y la fisiología de los microorganismos ambientales. Se han desarrollado una serie de algoritmos y herramientas que permiten generar secuencias de mayor tamaño llamados contigs apartir de secuencias cortas (100 pb, 300 pb o más). Estos contigs pueden unirse en secuencias más largas y generar scaffolds, que se espera contengan al menos secuencias completas de genes. La mayoría de los algoritmos trabajan con k-mers, que son secuencias de oligonucleótidos de tamaño k, también llamados l-tuples o n-grams (Dubinkina et al., 2016). Según estos autores, este método es más simple y rápido que el análisis basado en referencias y alineamientos, ya que los k-mers son la representación comprimida de las secuencias. Existen al menos tres enfoques para el ensamblaje de secuencias: el paradigma de Greedy, la superposición (OLC) y los grafos de Bruijn, siendo este último el más utilizado al considerarse mejor para ensamblar secuencias cortas (< 100 pb) (Li et 12 al., 2012, Nagarajan y Pop 2013). Los algoritmos OLC y grafos de Bruijn utilizan grafos, que son un conjunto de objetos llamados vértices (o nodos) unidos por aristas (edges) que pueden representar relaciones binarias entre elementos de un conjunto. Típicamente, un grafo se representa mediante una serie de puntos (los vértices) conectados por líneas (las aristas). Estos permiten generar posibles soluciones al ensamblaje mediante conexiones de tal manera que las secuencias se vinculan para formar secuencias contiguas (contigs). El enfoque OLC fue muy exitoso para ensamblar secuencias generadas a partir de la tecnología de Sanger (> 200 pb). Algunos de los ensambladores que utilizan este algoritmo son CAP3 (Huang y Madan, 1999), Celera (Myers et al., 2000), Arachne (Batzoglou et al., 2002), PCAP (Huang y Yang, 2005) y Newbler (Margulies et al., 2005). La mayoría de los ensambladores mencionados están en desuso. Los grafos de Bruijn es un algoritmo propuesto por Ramana y Waterman en 1995 (Idury y Waterman, 1995), el cual utiliza k-mers para inferir la secuencia de ADN ensamblada. Existen dos tipos de métodos que utilizan los grafos de Bruijn, los sencillos en los que se utiliza un solo tamaño de k-mer como solución del ensamblaje y los múltiples, los cuales combinan los resultados de diferente tamaño de k-mers para generar el ensamblaje, lo que mejora el resultado de este (Pell et al., 2012). El primer ensamblador que usó este algoritmo fue EULER, en el año 2001 (Pevzner y Waterman, 2001), y actualmente la mayoría de los ensambladores utilizados en análisis de metagenomas utilizan modificaciones de este, tales como la eliminación de lecturas erróneas para múltiples valores de k-mers o alineamiento de otras secuencias, corrección de llamado de bases erróneas basado en la frecuencia de k-mer o el valor de calidad de cada base, partición del grafo, detección de contigs mal ensamblados, entre otros (Miller et al., 2010, Vollmers et al., 2017, Sczyrba et al., 2017, van der Walt et al., 2017, Forouzan et al., 2018). En el Cuadro I se mencionan algunas de las características de los ensambladores Metavelvet (Namiki et al., 2012), MetaSPAdes (Nurk et al., 2017) y MEGAHIT (Li et al., 2015), todos estos ensambladores de novo, o sea que no utilizan un genoma de referencia como molde para ensamblar los fragmentos de ADN y generar secuencias más largas. Debido a que no se utiliza un genoma de referencia con el cual se pueda comparar, el rendimiento del ensamblaje de novo debe considerar 13 qué tanto del total de secuencias se incluye dentro en el ensamblaje y qué tanta información del genoma se puede encontrar (Baker 2015). MEGAHIT es actualmente el ensamblador de novo más eficiente en términos de memoria, ya que para ensamblar metagenomas se necesita menor uso de memoria RAM (Memoria de Acceso Aleatorio) (Olson et al., 2017). La memoria se reduce al eliminar k-mers por debajo de un umbral de frecuencia definido. Este enfoque minimiza el impacto negativo de los errores de secuencia en el ensamblaje. Para conservar los k-mers de organismos de baja abundancia, distinguiéndolos de los errores, MEGAHIT reconsidera a los k-mers eliminados en las regiones de baja cobertura en el ensamblaje. MetaSPAdes es una versión específica para metagenomas del ensamblador SPAdes (Nurk et al., 2017). Una de las principales innovaciones de este ensamblador es el uso de información de extremos pareados durante el proceso de ensamblaje en lugar de hacerlo posteriormente. Esta información se incorpora en el gráfico utilizando un par de k-mers separados por una distancia estimada (Olson et al., 2017). En MetaSPAdes las micro-variaciones entre posibles cepas muy similares se combinan para formar secuencias de consenso de alta calidad, lo que genera una mejor representación de cada “especie”. Existen algunos factores que afectan directamente el desempeño de los ensambladores que utilizan grafos de Bruijn, como los errores de secuenciación, las regiones repetidas, la cobertura de secuenciación y la secuenciación de cepas muy relacionadas genéticamente. Las regiones repetidas de ADN pueden encontrarse en el mismo organismo o estar compartidas entre distintos organismos. Las repeticiones crean un mayor número de conexiones (edges) incrementando el número de posibles soluciones, lo que ocasiona ambigüedad en el ensamblaje. Por otro lado, el proceso de anotación consiste en la delimitación de características relevantes en una secuencia del genoma o metagenoma (Richardson y Watson, 2012). Dentro de las herramientas que existen para la predicción de genes bacterianos y de arqueas se encuentran los paquetes Glimmer (Delcher et al., 2007), GenemarkHMM (Lukashin y Borodovsky, 1998), Easygene (Larsen y Krogh, 2003) y Prodigal (PROkaryotic DYnamic programming Gene-finding ALgorithm, Hyatt et al., 2010). Este último algoritmo ha sido ampliamente utilizado en estudios 14 de metagenomas (Jungbluth et al., 2017, Tamames y Puente-Sánchez 2018, Tully et al., 2018). En 2018, se desarrolló METAWRAP el cual incluye numerosos módulos para el análisis de MAGs, incluida la asignación de taxonomía, estimación de abundancia, anotación funcional y visualización de datos de secuenciación tipo shotgun (Uritskiy et al., 2018). Existen otros algoritmos de predicción de genes ab initio los cuales utilizan propiedades de las regiones codificantes como características y una amplia variedad de técnicas matemáticas como redes neuronales para encontrar las regiones codificantes (Do y Choi 2006). Por otro lado, algunos de los algoritmos de anotación más implementados utilizan un conjunto de secuencias de referencia para entrenar un modelo, y luego utiliza ese modelo para predecir las regiones de codificación del genoma de interés. Algunos de estos algoritmos se han implementado en programas de predicción de genes de código abierto, como Prodigal, el cual identifica CDS (región codificante), sin asignar alguna función. Para la predicción de los marcos abiertos de lectura (ORF, por sus siglas en inglés) y de ARNs, la herramienta Prodigal es entrenada previamente con secuencias. En su modo básico este programa toma las secuencias que se le proporcionan, las estudia, aprende sus propiedades y luego es capaz de predecir los genes basándose en esas propiedades. Dentro de las propiedades necesarias se encuentran los codones de iniciación (ATG, GTG, TTG), el sitio de unión ribosomal (RBS), contenido de GC, estadísticas de codificación de hexámeros y otra información necesaria para construir un perfil de entrenamiento completo (Hyatt et al., 2010). El algoritmo tiene que determinar automáticamente un conjunto de genes "reales" putativos sobre los cuales entrenarse. Por otro lado, existe el modo anónimo, en donde se utiliza un archivo de entrenamiento precalculado, y predice los genes basándose en los mejores resultados. También existe el modo de entrenamiento en donde se aplica el modo normal, pero Prodigal guarda la información de entrenamiento para usarla en futuros análisis. Todas esas variantes lo que buscan es mejorar la predicción, el reconocimiento del sitio del inicio de la transducción y la reducción de falsos positivos. Las proteínas están generalmente compuestas por una o más regiones funcionales, comúnmente denominadas dominios. Las diferentes combinaciones de dominios 15 dan lugar a la diversidad de proteínas que se encuentran en la naturaleza. Por lo tanto, la identificación de los dominios dentro de las proteínas puede proporcionar información sobre su función. Posterior a la identificación de genes, se debe realizar la búsqueda de homología entre las secuencias y las bases de datos que sirven para los análisis de funciones como Pfam (colección de secuencias de familias de proteínas, Finn et al., 2016), COG (Clúster de Grupos Ortólogos, Tatusov et al., 2000) y KEGG (Enciclopedia de Kyoto de Genes y Genomas). Esta búsqueda se puede realizar utilizando la herramienta DIAMOND (Buchfink et al., 2015) o BLASTX (Altschul et al., 1990), lo que permite posteriormente asignar una clasificación funcional a cada gen y buscar variantes alélicas, que resulta de utilidad para la descripción de enzimas y vías metabólicas. La base de datos Pfam es una gran colección de dominios y familias de proteínas, cada una representada por múltiples alineaciones de secuencias y modelos ocultos de Markov (HMM) (El-Gebali et al.,2019). La versión 34.0 (marzo 2021) que es la más reciente contiene un total de 19179 familias proteicas las cuales se basan en datos de secuencias de UniProtKB, NCBI y en secuencias de proyectos seleccionados de metagenómica (Mistry et al., 2020). Esta base tiene 420933 secuencias de 199 arqueas, 20190161 secuencias de 5337 bacterias, 20073593 secuencias de 1009 eucariotas y 370423 secuencias de 7459 virus. Por otro lado, la base COG, es un intento de clasificación filogenética de las proteínas codificadas en genomas completos de bacterias, arqueas y eucariotas (Tatusov et al.,1997, Galperin et al., 2015, Makarova et al., 2015). Se basa en una comparación exhaustiva de todas las secuencias de proteínas, lo que ha generado una asignación de genes ortólogos (misma función en especies diferentes) y paralógos (homólogos en una misma especie, función distinta). Además, incluye una cuidadosa curación manual de los COGs, lo que evita errores de anotación y sobrepredicción (Galperin et al., 2015). La base de datos COG se creó en 1997 y la actualización más reciente es de enero del 2021. La última versión, incluye genomas completos de 1187 bacterias y 122 arqueas, y 4877 COG, para un total de 3236575 genes únicos mapeados en 3455867 loci genómicos, lo que representa el 73,5% de las 4401819 proteínas. Los planes futuros incluyen una mayor expansión de la colección COG agregando COG de arqueas (arCOG), dividiendo 16 los COG que contienen múltiples parálogos y el refinamiento continuo de las anotaciones COG (Galperin et al., 2021). Por otra parte, la base de datos KEGG es utilizada como referencia para la interpretación biológica de secuencias genómicas y otros datos generados a partir de NGS. Esta base se ha desarrollado para comprender la conservación y la variación de genes y genomas (Kanehisa et al., 2019). En particular, el sistema KO (KEGG Orthology) para genes ortólogos funcionales se ha desarrollado para representar características conservadas de genes y proteínas. La base de conocimientos de referencia de los mapas de ruta de KEGG dibujados como redes de nodos KO representa las características conservadas de procesos celulares, lo que permite la reconstrucción automática de las vías metabólicas y los genes involucrados. KEGG se ha convertido, según Kanehisa y colaboradores (2019), en un recurso ampliamente utilizado para la interpretación biológica de diferentes tipos de datos como genomas, transcriptomas, metabolomas y metagenomas. La implementación de todas las herramientas de análisis de metagenomas arriba mencionadas han permitido generar mucha información acerca de la diversidad y metabolismo de los microorganismos en diferentes ambientes. Dentro de los principales estudios de metagenómica funcional se encuentra el análisis de 20 metagenomas de sitios calientes y ácidos del Parque Nacional Yellowstone (Inskeep et al., 2013), 7 metagenomas provenientes de fuentes termales de la India (Saxena et al., 2017), además de otros estudios realizados en Taiwán (Lin et al., 2015), India (Mehetre et al., 2016b, Mangrola et al., 2018), Japón (Sato et al., 2017) y España (López-López et al., 2015). Todos estos estudios se han enfocado en la caracterización funcional de las comunidades y la asociación de éstas con variables ambientales en cada sitio, así como a la elucidación de relaciones ecológicas en los consorcios microbianos. Enzimas hidrolíticas termoestables Una enzima o proteína se considera termoestable cuando tiene una temperatura de desnaturalización mayor a 55 °C, o una vida media larga a alta temperatura (Turner et al., 2007). Las enzimas hidrolíticas aceleran las reacciones en las que un compuesto se rompe en componentes más simples. La adaptación de 17 biomoléculas a condiciones extremas está relacionada con la estabilidad y la flexibilidad de optimizar los estados funcionales de las proteínas (Cuadro II). La adaptación de proteínas de organismos extremófilos a ambientes extremos ha sido ampliamente estudiada (Ueno et al., 2016). Las enzimas termoestables se aíslan principalmente de organismos termófilos y poseen un gran número de aplicaciones debido a su estabilidad general inherente (Demirjian et al., 2001). Los avances en el estudio de estas enzimas han sido posibles gracias al aislamiento de microorganismos termófilos de diferentes ambientes y su posterior extracción y purificación (Burrows, 1973, Antranikian et al., 1987, Groboillot, 1994, Bhushan y Hoondal 1998, Bauer et al., 1999, Kohilu et al., 2001). Los organismos de ambientes termales han demostrado ser fuente de enzimas con actividad hidrolítica como proteasas, lipasas, glicosidasas, amilasas e hidrolasas principalmente (Aanniz et al., 2015, Dalmaso et al., 2015, De Castro et al., 2016, Mohammad et al., 2017, Sahay et al., 2017, Baltazi et al., 2017). Muchas de estas enzimas juegan un rol importante en las comunidades, facilitando la degradación de nutrientes y permitiendo el flujo de energía entre autótrofos y heterotrófos. Es importante mencionar que las enzimas más representativas involucradas en la biorremediación incluyen hidrolasas, esterasas, proteasas y lipasas, las cuales han mostrado un potencial para la degradación de diversos compuestos como polímeros, hidrocarburos aromáticos, compuestos halogenados, tintes, detergentes, compuestos agroquímicos (Bhandari, et al., 2021). 18 Peptidasas En la literatura existe ambigüedad con respecto a los términos relacionados con las enzimas que degradan péptidos o proteínas. El concepto de "enzima proteolítica", "proteasa", "proteinasa" y "peptidasa" a menudo son utilizados como sinónimos. El término "proteolítico" fue utilizado por primera vez en 1877 por Michael Foster (Rawlings, 2013) y aplicado a la enzima pepsina. La "proteólisis" fue documentada por primera vez en 1880 (Roberts 1880), cuando se realizaron estudios de la digestión en el estómago. Vines acuñó el término "proteasa" en 1903 (Vines 1903), posteriormente "proteinasa" se usó como sinónimo de "proteasa" (Fischer 1907). Finalmente, "peptidasa" se utilizó por primera vez en 1918 (Petersen y Short 1918), y en 1923 se distinguió una peptidasa de una verdadera proteasa: una proteasa degradadas proteínas, mientras que una peptidasa hidroliza solo péptidos o "peptonas". Estas distinciones no son aceptables en la actualidad, ya que se consideran términos similares (Cerdá-Costa y Gomis-Rüth, 2014). En este trabajo utilizaremos el término peptidasa para referirnos a las proteasas, proteinasas y peptidasas. Las peptidasas se pueden clasificar según su especificidad en endopeptidasas y exopeptidasas o por el tipo de catálisis que realizan (Fig. 2) También se pueden clasificar en peptidasas extracelulares o intracelulares. Existen clasificaciones por homología y evolución. Una clasificación ampliamente aceptada es la propuesta en MEROPS (https://www.ebi.ac.uk/merops/cgi- bin/clan_index?type=P), la cual clasifica conjuntos homólogos de peptidasas e inhibidores de proteínas jerárquicamente en especies, familias y clanes de proteínas en función de la similitud de secuencia, las distancias evolutivas y el tipo de residuo que utilizan para la catálisis de los sustratos peptídicos (versión más reciente v12.3). Ésta es similar a la encontrada en UniProt (https://www.uniprot.org/). Muchas de estas familias de peptidasas se encuentran ampliamente distribuidas en los diferentes filos de bacterias, sin embargo, algunas se pueden encontrar solamente en arqueas o bacterias (Fig. 3). 19 Figura 2. Clasificación de peptidasas según su especificidad y tipo de catálisis que realizan Fuente: Elaboración propia Figura 3. Distribución de peptidasas encontradas en bacterias y/o arqueas en clanes y familias según MEROPS. Modificado de Nguyen et al., 2019. 20 Algunos estudios señalan que las metalopeptidasas y las peptidasas de serina son las más abundantes en bacterias y arqueas (Nguyen et al., 2019). Las peptidasas aspárticas son en su mayoría codificadas por hongos, las metalopeptidasas son comunes en las bacterias, y las peptidasas de cisteína y serina son universales en los microorganismos (Caldwell, 2005). Las proteinasas tienen diversos roles bioquímicos, fisiológicos y reguladores en las bacterias y arqueas, sin embargo, su abundancia difiere entre los dos grupos de microorganismos y existen algunas familias únicas en cada uno de ellos (Nguyen et al., 2019). La distribución de las peptidasas también varía en función de los microhábitats ecológicos ocupados por diferentes taxones microbianos y las condiciones ambientales (Rawlings et al., 2019). Lipasas Las enzimas hidrolíticas de éster carboxílico constituyen un gran grupo de enzimas que pueden catalizar la hidrólisis, síntesis o transesterificación de un enlace éster. Se pueden encontrar en bacterias y arqueas hipertermófilas (Levisson et al., 2009). Hay dos grupos bien conocidos dentro de la familia de las hidrolasas de éster carboxílico: lipasas y esterasas. Las esterasas difieren de las lipasas al mostrar preferencia por los ésteres de acilo de cadena corta (más cortos que 10 átomos de carbono) y que no son activos en sustratos que forman micelas (Chahinian et al.,2002). Otros grupos incluyen, por ejemplo, arilesterasas y fosfolipasas. Según Kovacic y colaboradores (2018), la mayoría de las secuencias de lipasas disponibles han sido identificadas en bacterias. A pesar de la baja similitud entre las secuencias de aminoácidos entre las enzimas lipolíticas de bacterias, actualmente se clasifican en 19 familias según criterios filogenéticos, secuencias conservadas y funciones biológicas (Figura 4) (Kovacic et al., 2018). Las lipasas participan en diferentes mecanismos celulares tales como el crecimiento de las células, la adhesión a diferentes estructuras, tienen funciones específicas junto a otras enzimas, además tienen un rol muy importante en las infecciones microbianas (Stehr et al, 2003, Hausmann y Jaeger, 2010). El papel más destacado de las lipasas extracelulares para un microorganismo es la digestión de lípidos para la adquisición de nutrientes (Stehr et al, 2003). http://www.frontiersin.org/people/u/658073 21 Figura 4. Árbol filogenético de las lipasas y esterasas clasificadas en 19 subfamilias. Las familias están ordenadas por números romanos. Las familias I-VIII descritos por Hausmann y Jaeger 2010 están resaltados por recuadros grises. Las familias restantes IX al XIX están coloreadas de azul y rojo, las lipasas no clasificadas se encuentran en cajas sin colorear. Tomado de Kovacic et al., 2018. 22 Amilasas Las amilasas se pueden clasificar en dos grandes grupos, las α-amilasas (EC 3.2.1.1), las cuales son enzimas extracelulares que hidrolizan los enlaces glicosídicos α-1,4 del almidón al azar, liberando dextrinas como subproductos, y las β-amilasas (EC3.2.1.2), exoenzimas que hidrolizan los enlaces α-1,4 del extremo no reductor, lo que provoca la inversión de la configuración anomérica de la maltosa liberada a su forma β (Mehta y Satyanarayana, 2016). Las α-amilasas muestran una amplia gama de degradación del sustrato, pueden degradar amilosa, amilopectina, ciclodextrinas, glucógeno y dextrinas, pero poseen la mayor especificidad hacia el almidón (Antranikian, 1992). Varios iones metálicos pueden influir en la actividad de estas enzimas, tal como lo es el Ca2+, ya que puede mejorar la actividad de la enzima y le proporciona estabilidad cuando la temperatura aumenta (Khajeh et al., 2001). Actualmente, las α-amilasas se clasifican en las familias GH13, GH57 y GH119, siendo que las α-amilasas de las familias GH57 y GH119 se encuentran únicamente de procariotas (Janeček et al., 2012). La familia GH13 es la principal familia de α-amilasas, consta de más de 30 diferentes especificidades enzimáticas, y junto con GH70 y GH77 forma el clan GH-H (Cantarel et al., 2009). Los miembros del clan GH-H comparten el dominio α (β / α) y puede ser reconocido porque posee de 4 a 7 regiones conservadas de aminoácidos que contienen tres residuos catalíticos (Chen et al., 2012). Existen bases de datos específicas como CAZy (The Carbohydrate-Active EnZymes database), la cual describe las familias que degradan, modifican o crean enlaces glicosídicos (Lombard et al., 2014). La mayoría de amilasas son extracelulares, sin embargo existen algunas que se encuentran dentro la célula, como lo son algunas encontradas en la bacteria hipertermófila Thermotoga maritima MSB8 (Ballschmiter et al., 2006) y la encontrada en Streptococcus mutans (Simpson y Russell, 1998). https://www.ncbi.nlm.nih.gov/pubmed/?term=Ballschmiter%2520M%255BAuthor%255D&cauthor=true&cauthor_uid=16517673 https://www.ncbi.nlm.nih.gov/pubmed/?term=Simpson%2520CL%255BAuthor%255D&cauthor=true&cauthor_uid=9721315 https://www.ncbi.nlm.nih.gov/pubmed/?term=Russell%2520RR%255BAuthor%255D&cauthor=true&cauthor_uid=9721315 23 Enzimas involucradas en la remoción de metales pesados y agroquímicos Un xenobiótico es una sustancia contaminante con características químicas y estructurales complejas para las cuales se conocen pocas enzimas con la capacidad para degradarlo, y que han sido liberados al medio ambiente debido a actividades industriales y agrícolas (Gianfreda y Rao 2017). Algunas de las vías metabólicas relacionadas con el metabolismo, transporte intracelular y/o degradación de metales pesados, halocarbonos, herbicidas y pesticidas han sido descritas previamente en la base de datos como BacMet (http://bacmet.biomedicine.gu.se/) y KEGG (https://www.genome.jp/kegg/). Los estudios basados en NGS permiten la identificación y caracterización de los genes que codifican por enzimas en los microorganismos independientemente de sistemas de cultivo a nivel de laboratorio. Esto amplía el panorama de estudio y podría revelar la presencia de enzimas con características únicas, entre las cuales están las involucradas en el metabolismo y degradación de xenobióticos (Saxena et al.,2017, Bao et al., 2017). La hidrólisis de plaguicidas por enzimas hidrolíticas de microorganismos puede servir como un mecanismo de desintoxicación o activación para la selectividad o resistencia de plaguicidas, ya que la especificidad de sustrato puede variar drásticamente entre microorganismos (Hoagland y Zablotowicz, 2001). Los estudios de metagenómica han permitido identificar y caracterizar enzimas hidrolíticas relacionadas con la degradación y metabolismo de contaminantes (Pushpanathan et al., 2014, Ufarté et al., 2015, Saxena et al., 2017, Jeffries et al., 2018), lo que ha generado nuevas herramientas e información relevante en aspectos como la biorremediación y otras aplicaciones industriales. Es de gran importancia mencionar que existen muchos agroquímicos que utilizan en sus formulaciones metales tales como cobre y arsénico. Existen en el mercado herbicidas de cobre quelatado, herbicidas para plantas acuáticas a base de cobre (Mastin y Rodgers, 2000, Wagner et al., 2016), así como compuestos antimicrobianos a base de este metal (Lamichhane et al., 2018). En el caso del arsénico, existen insecticidas y fungicidas utilizados en la agricultura como arsenito de sodio (Na2HAsO3), arsenito de calcio (Ca (AsO2)2), cobre-acetoarsenito- ParisGreen (Cu (CH3COO)2– 3 Cu (AsO2)2), piroarsenito de cobre (Cu2As2O5), https://www.researchgate.net/scientific-contributions/Robert-M-Zablotowicz-39740142 https://www.researchgate.net/scientific-contributions/Robert-M-Zablotowicz-39740142 24 arsenato de calcio (Ca3(AsO4)2), y arsenato de sodio (Na2HAsO4. 12 H2O) y arseniato de plomo (Bencko and Foong, 2017). El uso de arsénico en combinación con Pb, Ca, Mg, Mn, Fe y arsenicales inorgánicos está prohibido en Costa Rica desde 1999 (Servicio Fitosanitario del Estado, Ministerio de Agricultura y Ganadería, 2021). No obstante, existen herbicidas como glifosato, fungicidas como folpam y teldor e insecticidas como Pyrinex que pueden tener metales pesados como As, Br, Co, Pb y Ni en su formulación (Defarge et al., 2018). De estos compuestos se tiene reporte de importación a Costa Rica desde el año 2005 (Servicio Fitosanitario del Estado, MAG, 2021). La mayoría de los estudios de remoción de metales se basan en procesos de adsorción a nivel de la membrana celular, en el que interfieren fuerzas electrostáticas. No obstante, existen otros mecanismos como la bioacumulación, el cual es un proceso metabólicamente activo en el que los microorganismos absorben metales en su espacio intracelular utilizando complejos importadores que crean una ruta de translocación a través de la bicapa lipídica (sistema de importación) (Diep et al., 2018). Una vez dentro del espacio intracelular, los metales pueden ser secuestrados por proteínas y ligandos peptídicos (Malik, 2004, Mishra y Malik, 2013). Las bacterias se protegen de los compuestos tóxicos, para esto utilizan diversos mecanismos, entre ellos las bombas de eflujo. Algunos transportadores permiten la salida de un compuesto en específico o una clase de compuestos, y existen también transportadores de múltiples compuestos estructuralmente no relacionados (Ikeda y Yoshimura, 2002). Cinco superfamilias de bombas de flujo están asociadas a la resistencia/tolerancia de los microorganismos a diferentes drogas: 1) MDR: extrusión de múltiples fármacos y toxinas (MATE), 2) pequeña resistencia a múltiples fármacos (SMR), 3) superfamilia de facilitadores principales (MFS), 4) casete de unión a ATP (ABC) y 5) Resistencia-nodulación y división (RND). Según Alav y colaboradores (2018), las bombas de flujo de RND solo se han encontrado en bacterias Gram negativas y están organizadas como sistemas tripartitos que constan de una bomba de membrana citoplásmica, una proteína adaptadora periplásmica y un canal proteico de membrana externa. Todas las superfamilias de bombas de flujo utilizan energía de la fuerza motriz protón / sodio, https://www.researchgate.net/profile/Ilyas_Alav 25 excepto la superfamilia ABC, que son transportadores primarios que utilizan energía de la hidrólisis de ATP. Las bombas son un componente clave para el flujo de fármacos, que es uno de los principales mecanismos de resistencia a los antibióticos en las bacterias. En las bacterias Gram positivas, la superfamilia de bombas de flujo de MFS es la más estudiada. En las células bacterianas existen diferentes enzimas que regulan la concentración de los metales en su interior. Entre estas se encuentran las metalotioneínas, fitoquelatinas y metalochaperonas, las cuales pueden utilizar los metales como cofactores enzimáticos y transportadores de membrana, y además participan en procesos de homeostasis y desintoxicación de iones metálicos. El control de la expresión de genes que codifican por transportadores de metales y proteínas de almacenamiento generalmente se da mediante sensores de metales, que incluyen sistemas de dos componentes y más de siete familias de reguladores transcripcionales de unión a ADN y de unión a metales (Waldron y Robinson, 2009). En las bacterias, la limitación de metales activa las vías que participan en la importación y movilización, mientras que el exceso induce la salida y el almacenamiento (Chandrangsu, et al., 2017). En estudios previos de cepas de cianobacterias aisladas de fuentes termales de nuestro país se ha demostrado la capacidad de tolerar y remover herbicidas como paraquat y bromacil (Brenes-Guillén et al., 2017), y metales pesados tales como aluminio, cobre, plomo y arsénico y hierro (Núñez, 2012). Estos resultados sugieren que algunos de los microorganismos aislados de ambientes con altas temperaturas son capaces de tolerar altas concentraciones de estos compuestos, y además remover sustancias tóxicas, lo que resulta de gran utilidad para el desarrollo de nuevas estrategias en biorremediación. Además del hecho que la mayoría de las aguas termales tienen alto contenido de algunos metales, por lo que las estrategias de resistencia a estos ya están activas en los microorganismos que ahí se encuentran. https://www.ncbi.nlm.nih.gov/pubmed/?term=Chandrangsu%20P%5BAuthor%5D&cauthor=true&cauthor_uid=28344348 26 Justificación La metagenómica proporciona una serie de estrategias para evaluar el acervo genético completo de todos los microbios en un entorno particular. Con esta técnica se ha revela do una diversidad sin precedentes en la composición de la comunidad microbiana y la diversidad de funciones codificada en los genomas que revela genes y alelos nuevos. Las fuentes termales poseen nichos con condiciones extremas en los cuales se hipotetiza que se asemeja al ambiente en que se originó la forma de vida primitiva, lo que los convierte en recursos muy particulares para la identificación y caracterización de enzimas, así como sus variantes. Estos ambientes, con condiciones extremas de temperatura, pH y altas concentraciones de sales e iones, propician que los microorganismos que los habitan tengan genes que codifiquen para enzimas relacionadas con la sobrevivencia a dichas condiciones. Es por eso por lo que, son candidatos atractivos para la identificación y búsqueda de secuencias de ADN que codifican para moléculas como las enzimas con actividad hidrolítica, involucradas en la degradación de compuestos contaminantes como xenobióticos. La presencia de genes con capacidad de codificar enzimas que degradan compuestos complejos ha sido ya descrita en algunos estudios, sin embargo, no se ha realizado en nuestro país, cuyas fuentes termales tienen comunidades microbianas, diferentes tanto molecular como fenotípicamente a los reportados en otras fuentes termales. La identificación y clasificación de las enzimas hidrolíticas y sus variantes alélicas provenientes de microorganismos termotolerantes ofrecerá una oportunidad para caracterizar estas secuencias y generar nuevas investigaciones. Las bases de datos actuales son un reservorio importante para comparar los datos generados en fuentes termales costarricenses con los de otros sitios, lo que permitiría realizar análisis de la micro-diversidad de ambientes extremos, así como comparaciones funcionales y metabólicas entre las diferentes muestras. Además, la información generada permitirá incrementar el conocimiento a nivel ecológico y metabólico de 27 estos ambientes, lo que justificaría futuros esfuerzos de investigación para la exploración y conservación de los recursos biológicos disponibles en las fuentes termales. Objetivos Objetivo General Analizar y comparar datos de secuenciación de seis metagenomas de fuentes termales de Costa Rica y otros metagenomas públicos en búsqueda de secuencias de potenciales enzimas hidrolíticas. Objetivos Específicos 1- Comparar la composición microbiana y anotación funcional de las fuentes termales costarricenses con 13 metagenomas de fuentes termales del Parque Nacional Yellowstone, Nueva Zelanda y Chile según las características fisicoquímicas como temperatura, pH y concentración de iones. 2- Identificar y clasificar secuencias de enzimas hidrolíticas tipo proteasas, lipasas, amilasas y esterasas encontradas en seis fuentes termales costarricenses. 3- Clasificar y describir las enzimas posiblemente involucradas en la remoción de metales y agroquímicos en seis fuentes termales costarricenses 28 Materiales y Métodos Características de las muestras Las muestras utilizadas en este análisis pertenecen a la base de datos de metagenomas y microbiomas del Área de Microbiología Ambiental del CIBCM. Esta base está conformada por 14 metagenomas obtenidos a partir de secuenciación tipo shotgun de tapetes microbianos de diferentes partes del país y 30 muestras, de las cuales se tienen datos de secuenciación de amplicones 16S ARNr y algunas de estas, secuenciación del gen 18S ARNr. En la figura 5, se muestra una fotografía de cada una de las fuentes termales. En el cuadro 1, se muestra la ubicación, temperatura y pH de cada una de las fuentes termales de donde fueron colectados los tapetes microbianos utilizados en este estudio. Cuadro 1. Localización geográfica, temperatura, pH de muestras utilizadas en este proyecto. Muestra Localización geográfica- asociación a volcán Temperatur a (°C) pH Rocas Calientes (RC1) Zona Sur 63.0 6.0 Las Lilas Volcán Rincón de la Vieja 74.0 6.0 Miravalles Volcán Miravalles 51.0 7.03 La Luz Cordillera de Tilarán 49.9 6.24 Salitral Cordillera de Tilarán 34.9 6.1 Río Naranjo Cordillera de Tilarán 59.8 6.53 29 Se extrajo el ADN de tapetes microbianos y agua de seis fuentes termales que han sido previamente caracterizados mediante secuenciación masiva de la región V3- V4 del gen 16S ARNr, y que poseen análisis fisicoquímicos asociados (Uribe-Lorío et al., 2019, Brenes-Guillén et al, datos no publicados). Para cada una de las muestras se tienen datos ambientales como pH y temperatura in situ, los cuales fueron obtenidos mediante una sonda multiparamétrica (Combo, Hanna Instruments). La extracción de ADN ambiental se realizó utilizando el kit de extracción Nucleospin Soil (Mackerey-Nagel) de acuerdo con el protocolo de la casa comercial. La calidad del ADN fue visualizada mediante electroforesis en un gel de agarosa al 1%. La concentración de ADN, así como la pureza relación 260/280 y 260/230 se cuantificó utilizando un espectrofotómetro Nanodrop 2000/2000c. El ADN extraído fue almacenado en un congelador a -20 °C. Su usó el TruSeq Nano DNA Kit para la preparación de librerías y se secuenciaron en la plataforma HiSeq (IIlumina Inc., USA) por servicio de la empresa Macrogen Inc. (Corea). Las estadísticas generales de las secuencias se muestran en el Cuadro 2. Los archivos SRA (Sequence Read Archive) de las muestras de Nueva Zelanda, Parque Nacional Yellowstone (PNY) y Chile fueron obtenidas de GenBank, sus características generales se muestran en el cuadro 3. Para algunas de las muestras se obtuvieron fotografías de diferentes fuentes (Fig. 6). 30 Cuadro 2. Estadísticas de los datos crudos de secuenciación de las muestras de las fuentes termales de Costa Rica. Las columnas Q20 y Q30 corresponden al porcentaje de bases con un puntaje de calidad de 20 y 30 respectivamente. ID Total de bases Total de secuencias (paired-end) GC (%) AT (%) Q20 (%) Q30 (%) La Luz 16,368,673,878 162,066,078 55.25 44.75 97.49 92.21 Miravalles 9,027,602,402 89,382,202 55.03 44.97 96.55 89.60 Las Lilas 8,939,704,930 88,511,930 45.47 54.53 96.70 90.60 Salitral 14,987,773,194 148,393,794 55.32 44.68 97.18 91.49 Río Naranjo 13,632,200,280 134,972,280 51.14 48.86 97.26 91.60 Rocas Calientes 9,627,099,416 95,317,816 60.12 39.88 96.60 91.30 31 Cuadro 3. Características generales de cada muestra utilizada en el análisis. ID Bioproyecto Secuenciador Sitio Temp. °C pH SRR5248290 PRJNA366268 Illumina HiSeq 2500 Obsidian Pool, Parque Nacional Yellowstone, Estados Unidos 42-90 ∼ 6.7 SRR5248299 PRJNA366343 Illumina HiSeq 2000 Obsidian Pool, Parque Nacional Yellowstone, Estados Unidos 42-90 ∼ 6.7 SRR5248365 PRJNA366346 Illumina HiSeq 2000 Obsidian Pool, Parque Nacional Yellowstone, Estados Unidos 42-90 ∼ 6.7 SRR5248302 PRJNA366345 Illumina HiSeq 2000 Obsidian Pool, Parque Nacional Yellowstone, Estados Unidos 42-90 ∼ 6.7 SRR5451033 PRJNA382437 Illumina HiSeq 2000 Tapete microbiano, Porcelana, Los Lagos, Chile 48 ∼ 6.5 SRR5451032 PRJNA382437 Illumina HiSeq 2000 Tapete microbiano, Porcelana, Los Lagos, Chile 58 ∼ 6.5 SRR5451031 PRJNA382437 Illumina HiSeq 2000 Tapete microbiano, Porcelana, Los Lagos, Chile 66 ∼ 6.5 SRR10063240 PRJNA412936 Illumina HiSeq 1500 Raoul Island (Green Lake), Nueva Zelanda 34.5 6.96 SRR10063241 PRJNA412936 Illumina HiSeq 1500 Raoul Island (Marker Bay Pool), Nueva Zelanda 29.9 5.9 SRR10063242 PRJNA412936 Illumina HiSeq 1500 Raoul Island (Eastern Pool), Nueva Zelanda 27.2 5.43 SRR5650826 PRJNA378813 Illumina NextSeq 500 Mound Spring, Parque Nacional Yellowstone, Estados Unidos 60-100 9 SRR5650808 PRJNA378813 Illumina NextSeq 500 Bijah Spring, Parque Nacional Yellowstone, Estados Unidos 60-100 7 SRR4030101 PRJNA336659 Illumina HiSeq 2000 Mammoth Hot Spring (Liberty Cap Streamers), Parque Nacional Yellowstone, Estados Unidos 72 6.5-7 Río Naranjo Illumina HiSeq 4000 Tapete microbiano Salitral, Costa Rica 34.9 6.53 https://www.ncbi.nlm.nih.gov/bioproject/PRJNA366268 https://www.ncbi.nlm.nih.gov/bioproject/PRJNA366343 https://www.ncbi.nlm.nih.gov/bioproject/PRJNA366346 https://www.ncbi.nlm.nih.gov/bioproject/PRJNA366345 https://www.ncbi.nlm.nih.gov/bioproject/PRJNA382437 https://www.ncbi.nlm.nih.gov/bioproject/PRJNA382437 https://www.ncbi.nlm.nih.gov/bioproject/PRJNA382437 https://www.ncbi.nlm.nih.gov/bioproject/PRJNA412936 https://www.ncbi.nlm.nih.gov/bioproject/PRJNA412936 https://www.ncbi.nlm.nih.gov/bioproject/PRJNA412936 https://www.ncbi.nlm.nih.gov/bioproject/PRJNA378813 https://www.ncbi.nlm.nih.gov/bioproject/PRJNA378813 32 Salitral Illumina HiSeq 4000 Tapete microbiano Río Naranjo, Costa Rica 59.8 6.1 La Luz Illumina HiSeq 4000 Tapete microbiano La Luz, Costa Rica 49.9 6.24 Miravalles Illumina HiSeq 4000 Tapete microbiano Miravalles, Costa Rica 51 7.03 Las Lilas Illumina HiSeq 4000 Tapete microbiano Las Lilas, Costa Rica 74 6 Rocas Calientes Illumina HiSeq 4000 Tapete microbiano Rocas Calientes, Costa Rica 63 6 33 Figura 5. Fotografías de fuentes termales costarricenses. A. Rocas Calientes. B. Río Naranjo. C. Miravalles. D. Las Lilas. E. Salitral. F. La Luz. A: Fotografía facilitada por Lorena Uribe. B, E y F: Fotografías facilitadas por Guillermo Alvarado Induni. C y D: tomadas de Morales 2008 y Brenes Guillén et al., 2019. 34 Figura 6. Fotografías de fuentes termales de Yellowstone, Nueva Zelanda y Chile A. Liberty Cap. B. Opsidian Pool. C. Mound Spring, D. Bijah Spring, E. Fuentes termales en Nueva Zelanda. F. Fuente termal Porcelana, Chile. En los recuadros se señalan las diferentes fuentes bibliográficas de donde se tomaron las fotografías. 35 Análisis Bioinformático Se realizó un análisis de los datos utilizando como referencia diferentes programas que han sido utilizados previamente para el análisis de datos de metagenómica shotgun. En la figura 7, se muestra un resumen de la metodología utilizada y a continuación se detallan las características de cada sección del análisis. 36 Figura 7. Diagrama general de la metodología de análisis utilizada en este trabajo. 37 Análisis de calidad de los datos de secuenciación La mayoría de los análisis se realizaron utilizando el clúster Kabre (CENAT, CONARE). Para el análisis de calidad de las secuencias se utilizó el módulo FastQC (v 0.11.7). Este análisis permitió obtener datos relacionados con la calidad de la secuencia, cantidad de secuencias y permitió visualizar los valores de calidad de las secuencias. Posteriormente, se utilizó el programa Trimmomatic v0.32 (USADELLAB.org) con los parámetros SLIDINGWINDOW:4:20 y MINLEN: 100, para eliminar las secuencias con baja calidad (phred score menor a 20), eliminar adaptadores, primers y, secuencias cortas (menor a 100 pb). Ensamblaje, alineamiento y anotación de metagenomas La estandarización del análisis bioinformático de metagenomas de las fuentes termales se realizó utilizando como base el conjunto de programas ANVIO (Universidad de Chicago) desarrollado por Eren y colaboradores (2015). Para el ensamblaje de los metagenomas se utilizaron dos ensambladores: MEGAHIT (Li et al., 2015) y MetaSPAdes (Nurk et al., 2017). Debido a que se utilizó un ensamblaje de novo, se evaluaron ciertas métricas basadas en el número de contigs y scaffolds generados, la proporción de secuencias que fueron ensambladas, y el tamaño en pares de bases de esas secuencias (porcentaje de secuencias mapeadas al ensamblaje mayor al 60%), para escoger aquellas con un tamaño mayor a 1000 pb y un tamaño de k-mer menor o igual a 255 (parámetro por default). La evaluación de las estadísticas que se generen en cada ensamblaje se realizó con el programa MetaQUAST (Mikheenko et al., 2016). Para el alineamiento de las secuencias y obtención de estadísticas como porcentaje de misassemblies se utilizó la herramienta Bowtie2 (buscadores locales y globales) (Langmead y Salzberg, 2012). La predicción de genes se realizó utilizando el programa PRODIGAL (Hyatt et al., 2010) y la herramienta HMMER (Finn et al.,2015) para mapear los resultados de las anotaciones proteicas obtenidas entre las bases de datos GO (Gene Ontology) y Pfam 32.0 (https://pfam.xfam.org). Se utilizaron los parámetros por defecto de Prokka que incluyen un valor de cut-off de 1 x 10 -6 (probabilidad que el alineamiento ocurra al azar y no por similitud, es significativo si es menor a 0.001) (Seemann, 2014). La herramienta HMMER se basa en modelos probabilísticos ocultos de Markov (HMM) e incorpora la clasificación taxonómica a las muestras utilizando las bases de datos 38 propuestas por Rinke y colaboradores (2013) y Campbell y colaboradores (2013). Para la identificación de proteínas ya reportadas en la literatura se utilizó el programa DIAMOND (Buchfink et al., 2014) y las bases de datos COG (versión 2014) y Pfam (versión 32.0). Todo este proceso se realizó tanto para las muestras costarricenses como para los datos de 13 fuentes termales de Yellowstone, Nueva Zelanda y Chile. Perfil taxonómico y funcional de las comunidades microbianas Para determinar el perfil taxonómico de las comunidades se utilizaron dos clasificadores:, Kraken2 dentro de OmicsBox (https://www.biobam.com/omicsbox/), MetaPhlAn2 (https://huttenhower.sph.harvard.edu/metaphlan) y MG-RAST (https://www.mg-rast.org/). El análisis en MG-RAST solamente se realizó para las muestras de Costa Rica. MG-RAST es una aplicación web en la cual se pueden elegir una serie de parámetros para el control de calidad. Luego, los datos crudos pasan automáticamente por varios pasos y, produce perfiles de abundancia a diferentes niveles taxonómicos, características funcionales y visualizaciones. Por otro lado, Kraken2, la última versión de Kraken, es un clasificador de secuencias que asigna etiquetas taxonómicas a lecturas cortas de ADN. Lo hace examinando los k-mers dentro de una lectura y consultando una base de datos con esos k- mers. Esta base de datos contiene un mapeo de cada k-mer en la biblioteca genómica de Kraken al ancestro común más bajo (LCA) (Wood et al., 2019). MetaPhlAn (Análisis filogenético metagenómico) es una herramienta computacional para perfilar la composición de comunidades microbiana, se basa en genes marcadores específicos de clado únicos identificados a partir de ~ 17.000 genomas de referencia (~ 13.500 bacterianos y arqueales, ~ 3500 virales y ~ 110 eucariotas) (Segata et al., 2012). Una vez obtenida la abundancia de cada filo por muestra, se realizó un análisis de estadística multivariada para determinar la similitud entre las muestras, para ello se utilizó como insumo solamente los resultados obtenidos en Kraken. Para ello los datos fueron normalizados utilizando la raíz cuadrada y se realizó una prueba del perfil de similitud (SIMPROF), recomendada para muestras no divididas a priori en grupos (muestreo, los valores obtenidos en el cálculo del índice de similitud de Bray-Curtis y diseño experimental). SIMPROF utiliza como insumo permutaciones para estimar el patrón de agrupamiento entre las muestras. Posteriormente se realizó la prueba SIMPER (Similarity percentage analysis) y BEST para determinar la contribución de cada uno de los filos a las diferencias entre las muestras. BEST permite determinar las variables que "explican mejor" el patrón multivariado en las variables de respuesta. La prueba SIMPROF se utiliza para decidir qué subgrupos en https://www.biobam.com/omicsbox/ https://huttenhower.sph.harvard.edu/metaphlan https://www.mg-rast.org/ 39 el análisis jerárquico son estadísticamente significativos. Posteriormente, se graficó el shadeplot, basado en los resultados de SIMPROF, en el cual se puede observar la abundancia de cada filo por muestra agrupadas, las líneas punteadas indican que el patrón de agrupamiento no tiene un soporte estadístico y por lo tanto no pertenecen al mismo grupo. El análisis de NMDS muestra el patrón de agrupamiento entre las muestras y el porcentaje de similitud entre las mismas según los análisis preliminares. Estos análisis se realizaron en el programa Primer 7 (Plymouth Routines In Multivariate Ecological Research, Clarke y Goyle 2015). Para realizar las comparaciones del perfil funcional de las comunidades, se utilizaron los datos generados por Prokka (Pfam y COG). Para ello se determinó el número de secuencias de proteínas dentro de las 17 categorías funcionales COG (Cuadro 4). Posteriormente se construyó una base de datos para comparar las diferentes fuentes termales, para esto se utilizaron los 2448 COGs encontrados y se graficó la abundancia de cada COG con el programa estadístico Primer7. Es importante mencionar que los pares COG se distribuyen en una amplia gama de identidad de secuencia y valor esperado. Ni el porcentaje de identidad de secuencia ni el valor esperado del número de aciertos esperados de calidad similar (puntuación) que podrían encontrarse por casualidad (e-value), pueden dar una idea completa de la relación entre las dos proteínas. Si dos proteínas tienen una identidad de secuencia entre 60% - 70%, tienen cerca del 90% de probabilidad (quizás más) de compartir el mismo proceso biológico para los niveles del índice de GO. 40 Cuadro 4. Categorías funcionales según COG. Clasificación Categoría Letra Almacenamiento y procesamiento de información Modificación y procesamiento de ARN A Traducción, estructura ribosómica y biogénesis J Transcripción K Replicación, recombinación y reparación L Metabolismo Conversión y producción de energía C Metabolismo y transporte de aminoácidos E Metabolismo y transporte de nucleótidos F Metabolismo y transporte de carbohidratos G Metabolismo y transporte de coenzimas H Metabolismo y transporte de lípidos I Metabolismo y transporte de sustancias inorgánicas P Biosíntesis, transporte y catabolismo de metabolitos secundarios Q Señalización y procesamiento celular Control del ciclo celular, división celular, partición cromosómica D Biogénesis de la pared celular / membrana / envoltura M Motilidad celular N Modificación postraduccional, proteínas y chaperonas O Mecanismos de traducción de señal T Tráfico intracelular, secreción y transporte vesicular U Mecanismos de defensa V Estructuras extracelulares W Citoesqueleto Z Poco caracterizado Función general R Función desconocida S Sin clasificación Mobiloma: profagos y transposones X 41 Genomas obtenidos a partir de (meta)genomas (MAGs) El ensamblaje fue utilizado para la obtención de MAGs. Los MAGs de los tapetes microbianos fueron obtenidos mediante la interfaz interactiva de Anvio versión 5.0 (Eren et al., 2015). Posteriormente se filtraron aquellos que tenían un porcentaje de completitud mayor al 60% y de redundancia menor al 10%. Posteriormente se corroboró esta información utilizando CheckM (Parks et al., 2015) y su afiliación taxonómica con GTDB-Tk (Chaumeil et al., 2020) en la plataforma KBase (Arkin et al., 2018). Una vez obtenidos los MAGs se pasó a hacer el mapeo de los diferentes genes relacionados con la remoción de metales con la metodología antes mencionada. Creación de base de datos Para la creación de las bases de datos de secuencias proteicas, se utilizaron los resultados de la anotación de los 19 metagenomas de fuentes termales. De cada archivo producto de la anotación de las proteínas, se construyó un archivo que incluía todas las secuencias en formato fasta. Se utilizó la herramienta MMseqs2 (Steinegger y Soeding et al.,2017) para seleccionar secuencias de proteínas representativas de cada una de las bases de datos creadas, utilizando los parámetros --min-seq-id 1 --cluster-mode 1 --kmer-per-seq 1 -- alignment-mode 3. Estos parámetros fueron seleccionados con el fin de agrupar secuencias con alta identidad (=100%), para ello se realiza un alineamiento que cubra al menos el 80% de las secuencias para crear el clúster. Se consideró una variante o alelo, aquella secuencia que tenga un porcentaje de identidad menor a 100% con otra secuencia en la base de datos y al menos un mismatch al realizar el mapeo. Posteriormente, se ordenó la clasificación de las enzimas generadas. Para ello se utilizó el esquema de clasificación numérica, basado en las reacciones químicas que catalizan (Comisión de Enzimas: número EC) (Cuadro 5). Se utilizaron diferentes términos de búsqueda para cuantificar y encontrar los diferentes grupos de proteínas dentro de la categoría de enzimas hidrolasas (EC3). Para la categoría de peptidasas se utilizaron los términos de búsqueda “peptidasa”, “proteinasa” y “proteas” en inglés. De cada metagenoma se extrajeron las secuencias clasificadas dentro de cada categoría. Esta información se utilizó como referencia para la construcción de la base de datos empleada en el mapeo de genes. 42 Cuadro 5. Esquema de clasificación numérica para enzimas, basado en las reacciones químicas que catalizan (Comisión de Enzimas: número CE Clasificación Descripción EC1 Oxidoreductasas Reacciones de oxidación/reducción y de transferencia de átomos de H, O o electrones desde un sustrato a otro. EC2 Transferasas Transferencia de un grupo funcional desde un sustrato a otro. El grupo puede ser metil-, acil-, amino- o fosfato. EC3 Hidrolasas Formación dos productos a partir de un sustrato. Facilitan la ruptura de los enlaces de C-C, C-O, C-N, S-S, S-N, S-P, C-P. EC4 Liasas Adición o eliminación no hidrolítica de grupos de los sustratos. Pueden romper los enlaces C-C, C-N, C-O o C-S. EC5 Isomerasas Isomerización de una molécula. EC6 Ligasas Unión de dos moléculas por síntesis de nuevos enlaces C-O, C-S, C-N o C-C con la rotura simultánea de ATP. EC7 Translocasas Catalizar el movimiento de iones o moléculas a través de las membranas o su separación dentro de las membranas. Mapeo y abundancia de variantes alélicas Posterior al establecimiento de la base de datos, se realizó el mapeo. Para el mapeo de genes se utilizó el protocolo establecido por Castro y colaboradores (2019). Se utilizaron las secuencias de nucleótidos correspondientes a las secuencias proteicas, para ellos se utilizó Bowtie2 con los parámetros -N 1 y alineamiento local, donde N corresponde al número de mismatches permitidos. Si se utiliza un valor más alto, el alineamiento será más lento pero aumenta la sensibilidad, ya que el alineamiento local no requiere que las lecturas se alineen de un extremo a otro. Las alineaciones locales se pueden "recortar" en uno o ambos extremos de una manera que optimice la puntuación de alineación, además, genera alineamientos en regiones con alta similitud o regiones conservadas de secuencias con diferente longitud. En el caso de las bases de datos que contenían muchas secuencias, se https://es.wikipedia.org/wiki/Oxidaci%25C3%25B3n https://es.wikipedia.org/wiki/Reducci%25C3%25B3n-oxidaci%25C3%25B3n https://es.wikipedia.org/wiki/%25C3%2581tomo https://es.wikipedia.org/wiki/Hidr%25C3%25B3geno https://es.wikipedia.org/wiki/Ox%25C3%25ADgeno https://es.wikipedia.org/wiki/Electr%25C3%25B3n https://es.wikipedia.org/wiki/Enlace_qu%25C3%25ADmico https://es.wikipedia.org/wiki/Isomer%25C3%25ADa https://es.wikipedia.org/wiki/Mol%25C3%25A9cula https://es.wikipedia.org/wiki/S%25C3%25ADntesis_qu%25C3%25ADmica https://es.wikipedia.org/wiki/Adenos%25C3%25ADn_trifosfato 43 procedió a realizar un mapeo preliminar por muestra para seleccionar las variantes de los diferentes grupos de enzimas con un porcentaje de cobertura mayor a 0.02 % del total de secuencias para la graficación de la abundancia de cada una por muestra, donde este 0.02 % abarca el 80 % de las secuencias mapeadas. Para efectos de este estudio se va a considerar abundancia como el número de secuencias mapeadas. Para realizar comparaciones más detalladas de las secuencias de proteínas, se utilizó la herramienta UniProt-Blast (https://www.uniprot.org/blast/) (The UniProt Consortium), utilizando la base de datos UniRef90 (UniProt Reference Clusters 90% idéntidad). Mapeo de genes de resistencia a metales pesados y agroquímicos en tapetes microbianos de fuentes termales. Para el mapeo de genes de resistencia a metales (capacidad de tolerar altas concentraciones de metales) y agroquímicos se utilizó como referencia las bases de datos BacMet y KEGG. BacMet (v2.0, marzo 2018) es una base de datos curada manualmente' de genes bacterianos, que se han confirmado experimentalmente que confieren a los microorganismos resistencia a metales y/o biocidas antibacterianos, con referencias completas. Además, incluye una base de datos de genes de resistencia putativos (Pal et al., 2014). Esta base contiene 7