Limitations of 16S rRNA gene as phylogenetic marker: a large-scale meta-omics analysis of plaque microbiota in periodontal diseases

INTERNATIONAL DOCTORAL SCHOOL OF THE USC Alba Regueira Iglesias PhD Thesis Limitations of 16S rRNA gene as phylogenetic marker: a large-scale meta-omics analysis of plaque microbiota in periodontal diseases Santiago de Compostela, 2022 Doctoral Programme in Dental Science  DOCTORAL THESIS Limitations of 16S rRNA gene as phylogenetic marker: a large-scale meta-omics analysis of plaque microbiota in periodontal diseases Alba Regueira Iglesias INTERNATIONAL PHD SCHOOL OF THE UNIVERSITY OF SANTIAGO DE COMPOSTELA PHD PROGRAMME IN DENTAL SCIENCE SANTIAGO DE COMPOSTELA 2022 DECLARACIÓNDELAUTOR/ADELATESIS        D./Dña.AlbaRegueiraIglesias  Títulodelatesis:Limitationsof16SrRNAgeneasphylogeneticmarker:alarge‐scalemeta‐omics analysisofplaquemicrobiotainperiodontaldiseases  Presentomitesis,siguiendoelprocedimientoadecuadoalReglamentoydeclaroque:  1) Latesisabarcalosresultadosdelaelaboracióndemitrabajo. 2) Deserelcaso,enlatesissehacereferenciaalascolaboracionesquetuvoestetrabajo. 3) Confirmoquelatesisnoincurreenningúntipodeplagiodeotrosautoresnidetrabajos presentadospormíparalaobtencióndeotrostítulos. 4) Latesiseslaversióndefinitivapresentadaparasudefensaycoincidelaversiónimpresaconla presentadaenformatoelectrónico. YmecomprometoapresentarelCompromisoDocumentaldeSupervisiónenelcasoqueeloriginalno estédepositadoenlaEscuela.  EnSantiagodeCompostela,02defebrerode2022.  Firmaelectrónica  AUTORIZACIÓN DE LOS DIRECTORES DE LA TESIS Limitations of 16S rRNA gene as phylogenetic marker: a large-scale meta-omic analysis of plaque microbiota in periodontal diseases Dª. Inmaculada Tomás Carmona D. Javier Tamames De la Huerta D. Víctor Manuel Arce Vázquez INFORMA/N: Que la presente tesis, se corresponde con el trabajo realizado por Dª. Alba Regueira Iglesias, bajo nuestra dirección, y a utorizamos su presentación , considerando que reúne l os r equisitos exigidos en el R eglamento de Estudios de Doctorado de la USC, y que como directores de esta no incurre en las causas de abstención establecidas en la Ley 40/2015. De acuerdo con lo indicado en el Reglamento de Estudios de Doctorado, declaramos también que la presente tesis doctoral es idónea para ser defendida en base a la modalidad Monográfica con reproducción de publicaciones, en los que la participación de la doctoranda fue decisiva para su elaboración y las publicaciones se ajustan al Plan de Investigación. En Santiago de Compostela, 2 de Febrero de 2022  “Nothing in life is to be feared, it is only to be understood. Now is the time to understand more, so that we may fear less.” Marie Curie Acknowledgments A mi tutora la Dra. Inmaculada Tomás, por la dedicación, minuciosidad, y pasión puestas en cada trabajo realizado, que han sido una fuente de inspiración y motivación constantes para querer dar lo mejor de mí. Poder contar con la capacidad y experiencia de una investigadora excepcional, en un ámbito de confianza y trabajo en equipo; ha sido un privilegio. Las enseñanzas académicas y humanas transmitidas a lo largo de estos años son incalculables. Gracias. A mis directores el Dr. Javier Tamames y el Dr. Víctor Manuel Arce, grandes investigadores, por su participación en el desarrollo de este trabajo. Al equipo de la Dra. María José Carrera, perteneciente al Centro de Investigación Singular en Tecnologías Inteligentes la Universidad de Santiago de Compostela (CiTIUS); especialmente a Carlos Balsa y a Lara Vázquez, por los conocimientos bioinformáticos y estadísticos, y por el tiempo aportados para que esta Tesis haya salido adelante. A la Dra. Marta Relvas por la aportación de muestras orales, y a la Dra. Manuela Alonso por llevar a cabo los análisis de laboratorio realizados en esta Tesis. A todas mis compañeras de la Unidad de Pacientes con Necesidades Especiales coordinada por la Dra. Tomás, por todos los momentos compartidos; y en especial, a Triana Blanco por su labor en la selección de los pacientes y en la recogida de muestras. Es un placer formar parte un equipo con personas tan maravillosas. A mi pilar fundamental: mi familia; a mis amigas, y a José, por haberme acompañado durante este largo y, a veces complicado, camino. Gracias por apoyarme cuando decidí seguir esta ruta y compartir conmigo la felicidad de cada pequeño paso hacia delante. Pero, sobre todo, gracias por estar ahí para levantarme cuando el sendero se acomplejaba y el entusiasmo desaparecía, y por siempre recordarme que puedo hacer lo que me proponga. 7 Resumo da Tese “Limitacións do xene ARNr 16S como marcador filoxenético: unha análise meta-ómica a gran escala da microbiota da placa nas enfermidades periodontais” INTRODUCIÓN O termo “enfermidades periodontais” refírese a unha serie de condición diferentes iniciadas pola biopelícula dental que afectan ós tecidos que rodean e soportan os dentes (1). A xenxivite é unha inflamación localizada que non se estende ó aparello de unión periodontal e é reversíbel reducindo os niveis de placa (2). Sen embargo, en individuos susceptíbeis, se non é tratada pode progresar a periodontite, a cal se caracteriza pola destrución gradual do aparello de unión do dente (1). Estas patoloxías xorden cando se perde o equilibrio entre a biopelícula microbiana e o sistema inmune debido ben a unha disbiosis ou a unha reacción esaxerada do hóspede ós microbios. Ademais, atópanse entre as enfermidades orais máis prevalecentes e con maiores consecuencias en todo o mundo (3), cunha prevalencia global da periodontite severa en 2015 do 7.4% (4). No mesmo ano, as estimacións para España revelaron que un 5% da poboación adulta de entre 34 e 44 anos e un 10% entre 65 e 74, tiveron bolsas periodontais profundas (≥6 mm) (5). A periodontite exerce un efecto negativo sobre a calidade de vida das persoas que a padecen, especialmente naquelas con periodontite severa, comprometendo aspectos relacionados tanto coa función como coa estética (6). Especificamente, se non se aplica ningún tratamento, os resultados poden levar á perda de dentes, un deterioro do rendemento mastigatorio, un peor estado nutricional, unha menor autoestima e, incluso, pode ter efectos negativos sobre a saúde xeral (7). En relación con isto, dende hai anos observouse que unha serie de enfermidades e afeccións sistémicas poden afectar o aparello de inserción periodontal. Na actualidade, unha gran cantidade de investigacións apoian esta relación bidireccional entre a periodontite e as enfermidades cardiovasculares (8), a diabetes mellitus (9), as afeccións ALBA REGUEIRA IGLESIAS 8 respiratorias (10), a artrite reumatoide (11), a enfermidade de Alzheimer (12) e os resultados adversos do embarazo (13). O inicio e progreso da periodontite están relacionados con múltiples factores etiolóxicos e factores de risco modificables e non modificables (7,14), sendo de especial importancia a interacción entre os microorganismos locais e a resposta inmune do hóspede. Como consecuencia, o desenvolvemento de enfoques terapéuticos efectivos require a identificación dos principais microbios e marcadores do hóspede asociados á periodontite. Poñendo o foco sobre o compoñente microbiano, o noso coñecemento da etioloxía e patoxénese das enfermidades periodontais cambiou ó longo do tempo grazas 1) á hipótese de que estas condicións están causadas por biopelículas e non por microbios nun estado planctónico, 2) ó emprego de conceptos ecolóxicos para estudar a microbiota oral (a poboación de microorganismos que coloniza unha parte do copo) e, especialmente, 3) ás melloras tecnolóxicas nos métodos empregados para estudar os microbios presentes nas mostras orais (15). XUSTIFICACIÓN E OBXECTIVOS Os enormes avances producidos no campo da microbioloxía nas últimas décadas debido ó desenvolvemento e implantación das tecnoloxías de secuenciación de próxima xeración son innegables (16). De xeito específico, a secuenciación do xene ARN ribosomal (ARNr) 16S, considerado por moitos como o marcador filoxenético definitivo grazas principalmente á súa presenza ubicua en bacterias e arqueas e á intercalación de zonas conservadas e variábeis (17); permitiu estudar as comunidades microbianas complexas como a oral a profundidades sen precedentes (18). Non obstante, os resultados de investigación poden verse afectados por múltiples fontes de posibles nesgos durante cada paso do fluxo de traballo da secuenciación do xene ARNr 16S (19-21). A selección do par de cebadores é un deses pasos (19-21). Os cebadores constrúense en base a secuencias de consenso, pero poden presentar discordancias con algúns taxons que poden levar á sobre- ou infrarepresentación dun grupo microbiano concreto (22). En consecuencia, o uso dun cebador non apropiado podería dar como resultado conclusións Resumo da tese 9 biolóxicas cuestionables sobre o nicho a estudar (22). Sen embargo, non hai ningunha análise sobre a cobertura dos cebadores empregados para detectar os microorganismos procariotas que habitan na boca humana, entendendo como cobertura a porcentaxe de coincidencias para un determinado grupo de secuencias/rango taxonómico. Por outro lado, algúns nesgos asociados co fluxo de traballo da secuenciación son limitacións inherentes do propio xene (23). De feito, nun primeiro exemplo, varios autores demostraron a existencia de múltiples copias do xene ARNr 16S nos xenomas procariotas (24- 29), o que afecta as estimanzas de abundancia baseadas en contas do xene de xeito que taxons cun menor número de xenes tenden a ser infraestimados e aqueles cun maior número son sobreestimados (24,27). En segundo lugar, as nove rexións variábeis do xene teñen diferentes graos de heteroxeneidade de secuencia (26,30). Ademais, algunhas especies diferentes poder ten amplicóns coincidentes (en inglés matching amplicons, MAs), definidos como aqueles cunha similitude do 100% e o mesmo número de nucleótidos. Outras especies, pola contra, comparten secuencias altamente similares, incluso por riba do comunmente empregado limiar do 97% para construír unidades operacionais taxonómicas (operational taxonomic units, OTUs) (27,31), o que significa que poden ser agrupadas de forma errónea no mesmo OTU. Este agrupamento afecta á construción das táboas de OTUs e, por extensión, ás asignacións taxonómicas e ós resultados de diversidade. Con todo, a pesar das cuestións mencionadas anteriormente, non se realizaron investigacións exhaustivas sobre a mellor maneira de avaliar o número de xenes intraxenómicos do ARNr 16S nas bacterias e arqueas que ocupan a cavidade oral, nin sobre o impacto que ten o cebador elixido para as diferentes rexións na detección de MAs ou de amplicóns moi similares de taxons distintos. Ademais, segue sendo necesario avaliar a cuestión de cantos taxons orais diferentes e que especies orais específicas poden agruparse erroneamente no mesmo OTU, en función do cebador utilizado. A comparación dos estudos do microbioma periodontal baseados na secuenciación é controvertida polas significativas diferenzas metodolóxicas en pasos relevantes dentro do fluxo de traballo típico. É amplamente coñecido que cada tecnoloxía de secuenciación funciona de xeito diferente na relación entre a lonxitude da lectura, o rendemento da secuencia e a taxa de error (21), sendo Illumina preferible a Roche 454 e Ion Torrent (19). Tamén, como mencionamos anteriormente, as diferentes rexións variábeis e, en consecuencia, os amplicóns ALBA REGUEIRA IGLESIAS 10 derivados delas teñen distintos graos de heteroxeneidade de secuencia (26,30). En consecuencia, parece cuestionable a comparación das secuencias de mostras orais obtidas mediante tecnoloxías de secuenciación e rexións xenéticas distintas. Con todo, ata a data ningunha investigación avaliou as secuencias obtidas en estudos sobre o microbioma periodontal presente en diferentes condicións de saúde, en particular as xeradas mediante a plataforma de alto rendemento Illumina e distinguidas pola rexión do xene ARNr 16S máis amplificada. Debido á falta de evidencia sobre as cuestións sinaladas anteriormente, esta Tese tiña os seguintes obxectivos: 1) Analizar in silico a cobertura dos pares de cebadores empregados nos estudos baseados na secuenciación da microbiota oral; para iso utilizaranse dúas bases de datos específicas para a boca que conteñen secuencias do xene ARNr 16S de especies bacterianas e de arqueas. 2) Analizar in silico o número de xenes ARNr 16S nos xenomas completos das especies bacterianas e de arqueas que habitan na boca humana. Ademais, avaliar como o uso de diferentes pares de cebadores dirixidos a rexións distintas afecta á detección de MAs de taxons diferentes, identificando así as especies orais que teñen MAs. 3) Analizar in silico o rendemento de diferentes pares de cebadores de distintas rexións para identificar distintas especies procariotas orais con valores de similitude do amplicón do xene ARNr 16S ≥97%, establecendo así as especies orais que poden agruparse erroneamente no mesmo OTU. 4) Analizar os perfís da comunidade microbiana na placa supraxinxival e subxinxival de 2045 doentes con diferentes condicións periodontais (sans, xenxivite, periodontite e periodontite tratada) en relación coa diversidade bacteriana, os patróns de redes de coocorrencia e os modelos preditivos; sendo as secuencias utilizadas da plataforma Illumina, cun enfoque na rexión 3-4, e tratadas co mesmo protocolo bioinformático. Resumo da tese 11 OBXECTIVO 1. Avaliación in silico e selección dos mellores cebadores do xene ARNr 16S para o seu uso na secuenciación de próxima xeración para detectar bacterias e arqueas orais 1.1 MATERIAL E MÉTODOS A través da base de datos PubMed e utilizando o software estatístico R (32) e o paquete RISmed (33), realizáronse procuras automáticas para elaborar unha listaxe de: 1) cebadores dos xenes ARNr 16S empregados para detectar e amplificar bacterias e arqueas en mostras orais antes da secuenciación masiva; e 2) especies de arqueas que habitan na boca humana. Tras aplicar técnicas avanzadas de análise de texto a tódolos resumos descargados utilizando o paquete tm de R (34), quedamos con 129 estudos sobre bacterias e 16 sobre arqueas que implicaban o uso de polo menos un cebador diferente do xene ARNr 16S, e con 53 artigos que contiñan información sobre especies orais de arqueas. Identificáronse un total de 444 cebadores do xene ARNr 16S: 204 directos (forward, F), 230 reversos (reverse, R) e 12 non identificados (unidentified, UI). Deles, 278 obtivéronse das procuras en PubMed e 166 extraéronse do artigo de Klindworth et al. (35). A tódolos cebadores asignóuselles un identificador único baseado na súa procedencia -"OP" para os cebadores orais e "KP" para os de Klindworth (35)- e a súa dirección (F, R ou UI), seguido dun número de tres díxitos. Trala comparación das secuencias 5'-3' de tódolos cebadores, identificáronse 75 coas mesmas secuencias; o que nos deixou con 369 cebadores do xene do ARNr 16S diferentes. Por outra banda, obtivemos 177 nomes diferentes de especies de arqueas orais. A base de datos de Escapa et al. (36), que inclúe 223.143 variantes da secuencia do amplicón (amplicon sequence variants, ASVs) das secuencias do xene ARNr 16S, contén erros de anotación que fan imposible calcular a posición correcta dos cebadores dentro de cada secuencia en caso de coincidencia. Co obxectivo de mellorala, desenvolvemos scripts en Python (37) e Bash (38). Primeiro, as secuencias do xene do ARNr 16S das ASV do mesmo nivel xerárquico separáronse en 769 arquivos fasta diferentes. Despois, inseriuse un identificador de especie a tódalas secuencias antes da xerarquía taxonómica. As secuencias da mesma xerarquía aliñáronse simultaneamente utilizando Clustal Omega (39) contra un conxunto de secuencias do xene ARNr 16S de Escherichia coli. Tódolos ocos creados por Clustal Omega (39) foron eliminados, salvo os inseridos dende o inicio ata o primeiro nucleótido de cada secuencia. Os ALBA REGUEIRA IGLESIAS 12 arquivos fasta aliñados combináronse nun único arquivo para crear unha base de datos de ASVs completamente aliñadas, sendo a posición un o primeiro nucleótido de E. coli J01859.1. Por último, recortáronse as secuencias aliñadas con bases nunha posición inferior á do primeiro nucleótido de J01859.1, e aquelas con nucleótidos por enriba da posición 2000. A continuación, na base de datos de nucleótidos non redundantes do Centro Nacional de Información Biotecnolóxica (National Centre for Biotechnology Information, NCBI) (40), procuramos os xenomas completos das 177 especies de arqueas orais. A través dun script en Python (37) e cos identificadores, puidemos descargar 193 xenomas de RefSeq (41) e oito de GenBank (42). O script completouse co módulo search_16S.py (43), baseado no algoritmo de Edgar (44), que nos permitiu: detectar e extraer as secuencias do xene ARNr 16S dos xenomas completos descargados, eliminar tódalas secuencias repetidas e almacenar as variantes identificadas nun arquivo fasta. O módulo e a integración da ferramenta "The Entrez Programming Utilities (E-utilities)" (45) en Biopython (46) permitiu obter e asignar fácil e automaticamente o rango taxonómico completo ós xenes. Ademais, buscáronse as secuencias do xene do ARNr 16S das especies sen identificadores xenómicos completos (47-49). Finalmente, tódalas secuencias do xene agrupáronse nun único arquivo fasta. As secuencias neste arquivo foron empregadas para facer BLASTN (50,51) contra a base de datos de nucleótidos non redundantes do NCBI (40). A continuación, utilizando unha cobertura de busca ≥98% e unha porcentaxe de identidade ≥99%, descargáronse as secuencias do xene ARNr 16S e as rexións aliñadas cos xenomas completos; e ambos tipos de secuencias foron tratadas como ASVs. A base de datos de arqueas orais creouse utilizando outro script, e contén 2842 ASVs. As secuencias da base de datos foron aliñadas e melloradas seguindo os mesmos pasos que na base de bacterias. Para levar a cabo a análise in silico dos cebadores, definíronse as coberturas a nivel de variante (variant coverage, VC): porcentaxe de coincidencias dun cebador concreto en relación co total de secuencias da base de datos; e a nivel de especie (species coverage, SC): porcentaxe de especies con coincidencias en polo menos unha das súas variantes de secuencia cando se utiliza un cebador concreto. As coincidencias entre cebadores e secuencias das bases de datos avaliáronse aplicando as expresións regulares do módulo regex (52) de Python (37). Resumo da tese 13 Os cebadores individuais cunha SC ≥75,00% foron escollidos e tódalas combinacións posibles entre F e R foron identificadas. Estimouse a lonxitude media entre estas dúas posicións para clasificar os pares de cebadores nunha das tres categorías de lonxitude media dos amplicóns: 1) curta (short, S)= 100 a 300 pares de bases (base pairs, bps); 2) media (medium, M)= 301 a 600; e 3) longa (long, L)= >600. Tódolos pares de cebadores obtidos foron avaliados con respecto ás bases de datos de bacterias e arqueas, o que nos permitiu determinar se eran específicos de dominio ou para ambos. 1.2 RESULTADOS Un total de 148 e 65 cebadores individuais tiveron valores de SC de bacterias e arqueas ≥75,00%, respectivamente. Tras aplicar os criterios de formación de pares de cebadores, 3993 combinacións bacterianas e 645 de arqueas foron posibles. Delas, 156 estiveron repetidas para ambos dominios, e o resto foron específicas de dominio. Pares de cebadores específicos de bacterias Na categoría de lonxitude S, 139 pares tiñan valores de SC bacteriana ≥95,00% (rango= 99,09% - 95,19%), mentres que 33 tamén tiñan unha SC de arqueas de 0,00%. Estes últimos amplificaban as rexións xénicas 3-4 ou 5-7 e tiveron valores de SC bacteriana que oscilaron entre o 97,92% e o 95,58%, o que significou que non se cubriron entre 16 e 34 especies de bacterias orais. Para a maioría deles, a lonxitude media dos seus amplicóns foi de ó redor de 186 (rango= 189 - 182). Destaca o par OP_F009-OP_R030 da rexión 5-7, cunha lonxitude media de 297 e un valor de SC bacteriano do 96,88%, polo que só 24 especies de bacterias orais non foron cubertas por este par. Na categoría de lonxitude M, 68 pares de cebadores tiñan valores de SC bacteriana ≥95,00% (rango= 98,83% - 95,06%), dos cales 45 tiñan unha SC de arqueas de 0,00%. Os seus valores de SC bacteriana tamén oscilaron no rango anterior o que significou que entre nove e 38 especies non foron cubertas. Ademais, estes pares dirixíanse ás rexións xénicas 3-5, 3-6 ou 4-7, e tiñan lonxitudes de lectura medias entre 566 e 454. Dos pares coas maiores lonxitudes medias de amplicón, os que proporcionaron a mellor cobertura foron, por orde: KP_F051- OP_R030; OP_F021-OP_R030; KP_F048-OP_R073; KP_F051-KP_R053; OP_F021- KP_R053; e OP_F050-OP_R073 (rango de SC bacteriana= 98,83% - 96,23%; rango de ALBA REGUEIRA IGLESIAS 14 lonxitude de lectura media= 566 - 546). Estes amplificaron as rexións 3-6 ou 4-7 e non cubriron entre nove e 29 especies de bacterias. Na categoría de lonxitude L, 20 pares de tiñan valores de SC bacteriana ≥95,00% (rango= 97,14% - 95,06%), e 17 tamén tiñan un valor de SC de arqueas de 0,00%. Estes últimos pares tiñan o mesmo rango de SC bacteriano e deixaban entre 22 e 38 especies sen cubrir. Todos eles dirixíanse á rexión xénica 3-7 e tiñan unha lonxitude media de lectura de entre 772 e 732. Os cebadores co mellor equilibrio entre a lonxitude media de lectura e a cobertura foron KP_F048- KP_R074 (SC bacteriana= 97,01%; lonxitude media de lectura= 767); e OP_F050-KP_R074 (96,36%; 766). Con todo, houbo opcións interesantes con lonxitudes >1000 bps e valores de SC bacteriana ≥90,00% (rango de SC bacteriana= 93,37% - 90,64%; rango de lonxitude media de lectura= 1066 - 1059). Neste sentido, KP_F048-KP_R060, KP_F048-KP_R076 e KP_F048- OP_R121 da rexión 3-9 tiñan lonxitudes de lectura medias de 1061, 1060 e 1060, respectivamente, e valores de SC bacteriana do 93,37%; deixando sen cubrir 51 especies de bacterias orais. Pares de cebadores específicos de arqueas Na categoría de lonxitude S, 12 cebadores tiñan valores de SC de arqueas ≥95,00% (rango= 98,45% - 95,36%). Deles, oito tiñan valores de SC bacterianas do 0,00%: OP_F066-KP_R013; KP_F059-KP_R013; KP_F016-KP_R002; KP_F018-KP_R003; OP_F066-KP_R006; KP_F018-OP_R102; KP_F059-KP_R006; e KP_F018-KP_R002. A súa SC de arqueas oscilou entre o 95,88% e o 95,36%, tiveron lonxitudes de lectura medias de 275 a 144, amplificaron as rexións xénicas 3 ou 5-6 e, por último, non cubriron entre oito e nove especies de arqueas orais. Dezanove pares de cebadores na categoría de lonxitude M tiñan valores de SC de arqueas ≥95,00% (rango= 97,42% - 95,36%). Entre eles, nove tiñan tamén un valor de SC bacteriana do 0,00%: KP_F018-KP_R031; KP_F018-KP_R032; KP_F018-KP_R035; KP_F018- OP_R020; KP_F018-OP_R070; KP_F020-KP_R006; KP_F020-KP_R013; KP_F016- KP_R032; e OP_F114-KP_R006. Estes amplificaban as rexións 3-5 ou 3-6 e tiveron unha lonxitude media de amplicón de 551 a 414. Os pares cubriron entre 95,88% e 95,36% das especies de arqueas, deixando entre oito e nove sen cubrir. Resumo da tese 15 Só un par de cebadores na categoría de >600 bps tiña un valor SC ≥95,00% na base de datos de arqueas: OP_F114-KP_R013; o cal tamén tiña un valor de SC bacteriana do 0,00%. Este par amplificaba a rexión 3-6, tivo unha lonxitude media de 679 e non detectou oito especies de arqueas. Vinte e sete pares de cebadores tiñan unha SC de arqueas ≥90,00%, unha SC bacteriana de 0,00% e unha lonxitude media >679, e 10 dos cales eran superiores a 1100 (rango de lonxitude media= 1131 - 681). Deles, o mellor equilibrio entre a cobertura e a lonxitude media do amplicón atopouse en: KP_F016-KP_R066; KP_F016-KP_R063; KP_F018- KP_R066; e KP_F018-KP_R063. A SC de arqueas foi do 92,78% para os dous primeiros pares e do 93,81% para os dous segundos, deixando 14 ou 12 especies, respectivamente, sen cubrir. Todos estes dirixíanse á rexión 3-9 e tiñan, por orde, lonxitudes medias de amplicón de 1129, 1128, 1119 e 1118. Pares de cebadores de bacterias e arqueas Dez pares da categoría S tiñan valores de SC bacteriana e de arqueas ≥95,00% (rango= 95,97% - 95,32%; e 99,48% - 97,94%, respectivamente). A súa lonxitude media oscilou entre 288 e 284 e todos amplificaron a rexión 4-5: KP_F020-KP_R031; KP_F020-OP_R070; KP_F020-KP_R032; KP_F020-KP_R035; KP_F020-OP_R020; KP_F020-KP_R038; KP_F020-OP_R010; KP_F020-OP_R014; KP_F020-OP_R036; e KP_F020-OP_R048. O número de especies bacterianas e de arqueas non cubertas por estes pares oscilou entre 31 e 36; e entre unha e catro, respectivamente. Na categoría M, dous cebadores tiñan valores de SC bacteriana e de arqueas ≥95,00%: OP_F114-OP_R070 (SC bacteriana= 95,58%; SC de arquea= 98,45%); e OP_F114-KP_R031 (95,71%; 98,45%). Ambos amplificaron a rexión 3-5 e tiveron lonxitudes medias de 460 e 457, respectivamente. Non cubriron 33 (OP_F114-KP_R031) ou 34 (OP_F114-OP_R070) especies bacterianas e tres de arqueas. Ó baixar o corte a SC ≥90,00%, atopamos seis pares cunha secuencia media máis longa. Entre eles destacou OP_F114-OP_R073, cunha lonxitude media de 549, dirixíase á rexión xénica 3-6 e presentaba uns valores de SC de bacterias e arqueas do 94,80% e o 93,30%, respectivamente. Non cubriu 40 especies de bacterias e 13 de arqueas. Ningún par de cebadores da categoría L tivo valores SC ≥95,00% en ningunha das bases de datos. Pola contra, 28 tiñan SC bacterianas e de arqueas ≥90,00% (rango= 94,54% - 90,64%; e 96,91% - 96,39%, respectivamente). Estes amplificaron as rexións 3-9, 4-9 ou 5-9, tiñan ALBA REGUEIRA IGLESIAS 16 lonxitudes medias entre 622 e 1063, e non cubrían de 42 a 72 especies bacterianas e de seis a sete de arqueas. A combinación de OP_F066 con KP_R060, KP_R076 e OP_R121 produciu os valores de cobertura máis altos, e todos eles dirixíronse á rexión 5-9. Estes cebadores non cubriron 42 especies de bacterias e 6 de arqueas. Con todo, as súas lonxitudes medias foron 623, 622 e 622, respectivamente. Os pares de cebadores formados por OP_F114 con KP_R060, KP_R076 ou OP_R121, da rexión 3-9, presentaban un mellor equilibrio entre os resultados de cobertura (SC bacteriana= 91,42%, SC de arqueas= 96,91%) e as lonxitudes medias das secuencias (1063, 1062 e 1062). Sesenta e seis bacterias e seis arqueas non foron detectadas. 1.3 CONCLUSIÓNS Tendo en conta as tres categorías de lonxitude media do amplicón, os pares de cebadores coa mellor cobertura estimada para detectar as bacterias orais dirixíronse ás rexións 3-4, 4-7 e 3-7, e foron: KP_F048-OP_R043 (posición do par de cebadores para E. coli J01859.1: 342- 529), KP_F051-OP_R030 (514-1079) e KP_F048-OP_R030 (342-1079). Para a detección de arqueas orais, os pares con mellor cobertura amplificaron as rexións 5-6, 3-6 e 3-6, e foron: OP_F066-KP_R013 (784-indefinido), KP_F020-KP_R013 (518-indefinido) e OP_F114- KP_R013 (340-indefinido). Os pares coa mellor cobertura dos dominios de bacterias e arqueas conxuntamente atopáronse nas rexións 4-5, 3-5 e 5-9, e foron: KP_F020-KP_R032 (518-801), OP_F114-KP_R031 (340-801) e OP_F066-OP_R121 (784-1405). Os pares de cebadores coa mellor cobertura identificados neste estudo non se atopan entre os máis empregados na literatura sobre o microbioma oral. Resumo da tese 23 Nos arquivos fasta, tódalas secuencias incluían un identificador de especie (SPn) e un identificador de variante (Vn) na súa cabeceira e, tamén, a cabeceira de cada secuencia incluía a xerarquía taxonómica ata o nivel da variante dentro de cada especie. Finalmente obtivéronse os amplicóns in silico de 186 especies de bacterias orais e 135 de arqueas. Por outro lado, desenvolveuse un script co wrapper NcbiblastnCommandline de Biopython (46) para manexar BLAST+ 2.11 (57) en modo local dende Biopython. Isto permitiu transferir facilmente os datos obtidos nos aliñamentos para a súa posterior análise en Python (37). Os parámetros de aliñación configuráronse para que fosen os mesmos que os predeterminados en MegaBLAST (58). Tódalas secuencias pertencentes ó mesmo arquivo fasta para o mesmo par de cebadores aliñáronse entre si; para iso, cada arquivo fasta inseriuse como suxeito e consulta en BLASTN (51) para obter a porcentaxe de similitude entre amplicóns in silico pertencentes a diferentes especies orais. Dos resultados obtidos, seleccionáronse os amplicóns in silico cunha cobertura de aliñación do 100% das secuencias de consulta e cun valor de similitude ≥97%. Dos aliñamentos obtidos, descartáronse os seguintes por non ser de interese: 1) amplicóns in silico co mesmo identificador único (SPn + Vn); 2) amplicóns in silico co mesmo identificador de especie; e 3) aliñacións duplicadas. Se dúas especies diferentes tiñan máis dun valor de similitude de amplicón in silico ≥97% (amplicon similarity value ≥97%, ASI97) entre elas, elixíase un ó azar. Os resultados dos pares de especies altamente similares almacenáronse empregando os módulos de Python (37) pandas (56) e xlsxwriter (59) . A continuación, creouse unha matriz de similitude para cada par de cebadores. Mediante un script en R (32) calculamos para cada par de cebadores 1) o número de especies con polo menos un ASI97 con outras especies; 2) o número total de ASI97 entre especies diferentes; 3) o número medio e máximo de ASI97 por especie. Tamén se estimou para cada par de cebadores a porcentaxe de especies detectadas (SC) e a porcentaxe de especies detectadas sen ASI97 (species coverage with no ASI97, SC-NASI97). ALBA REGUEIRA IGLESIAS 24 Este último parámetro empregouse como criterio para seleccionar aqueles cebadores asociados a un menor número de especies orais que puidesen estar agrupadas de forma errónea. Por último, describíronse os pares de especies que mostraban un ASI97 e avaliouse se pertencían a xéneros ou rangos taxonómicos superiores diferentes. 3.2 RESULTADOS Os pares de cebadores dirixidos ás bacterias obtiveron unha media de 91,88 (49,40%) especies bacterianas cun ASI97 e unha media de 153,46 ASI97 con especies distintas. No caso dos dirixidos a arqueas, estas cifras foron de 65,60 (48,59%) e 162,26, respectivamente. Se se exclúen os cebadores máis utilizados na literatura sobre o microbioma oral, os de lonxitudes de amplicón curtas (a diferenza das porcentaxes de SC) tiveron os valores máis baixos de SCNASI97 tanto para as bacterias (S= 39,54%) como para arqueas (S= 40,44%) en comparación cos cebadores de lonxitude media e longa (M= 45,82% e 46,35%, respectivamente; L= 48,39% e 44,32%, respectivamente). Polo que respecta ós pares de cebadores específicos para bacterias, o número de especies bacterianas cun ASI97 e o número total de ASI97 oscilou entre 37 e 32 cun dos cebadores máis utilizados, KP_F031-KP_R021 (M; SC-NASI97= 54,30%), e 120 e 277 con OP_F066- KP_R040 (S; SC-NASI97= 24,19%), respectivamente. Este último cebador tamén tivo o valor SC-NASI97 máis baixo, mentres que OP_F053-KP_R020 detectou o maior número de especies sen ASI97 (M; SC-NASI97= 65,05%). Ademais, excepto OP_F053-KP_R020, tódolos cebadores específicos de bacterias tiñan un número máximo de ASI97/especies superior a cinco (rango= 15 - 4 ASI97/especies). En canto ós pares de cebadores específicos para arqueas, o número de especies de arqueas cun ASI97 e o número total de ASI97 oscilou entre 24 e 96 co amplamente utilizado KP_F014- KP_R011 (L; SCNASI97= 12,59%) e 89 e 240 con OP_F066-KP_R013 (S; SC-NASI97= 29,63%), respectivamente. O primeiro cebador detectou o menor número de especies sen ASI97, e KP_F018-KP_R002 o maior (S; SCNASI97= 51,11%). Ademais, tódolos cebadores específicos para arqueas tiñan un número máximo de ASI97/especie ≥10 (rango= 13 - 10 ASI97/especie). Resumo da tese 25 Empregando os pares de cebadores para bacterias e arqueas, o número de especies bacterianas e de arqueas cun ASI97 e o número total de ASI97 oscilou entre 84 e 60 e 118 e 126, respectivamente, con OP_F114-KP_R002 (S; SC-NASI≥97= 47. 31% para bacterias e 54,81% para arqueas) a 124 e 95 e 239 e 286, respectivamente, con OP_F066-OP_R073 (S; SC-NASI≥97= 31,18% para bacterias e 22,96% para arqueas). Este último cebador tamén detectou o menor número de especies sen ASI97 e OP_F114-KP_R031 o maior (M; SCNASI97= 51,08% para bacterias e 53,33% para arqueas). A maioría das combinacións de cebadores de bacterias e arqueas tiñan un número máximo de ASI97/especies ≥10 (rango= 14 - 9 ASI/especies e 14 - 11 ASI/especies para ámbolos dous dominios, respectivamente). Doutra banda, 149 (80,11%) das especies de bacterias orais e 108 (80,00%) das especies de arqueas orais avaliadas tiñan un ASI97 con polo menos unha especie distinta. Entre elas destacan pola súa relevancia na cavidade oral as bacterias: Campylobacter concisus, Campylobacter curvus, Rothia dentocariosa, Streptococcus mitis, Streptococcus mutans, Streptococcus oralis, e Tannerella forsythia; e as arqueas: Halovivax ruber, Methanosalsum zhiliniae, Methanosarcina barkeri, Methanosarcina mazei, e Methanosarcina vacuolata; entre outras. Ademais, houbo 30 especies bacterianas e 27 de arqueas que puideron agruparse cun máximo de ≥10 especies diferentes cando se empregaron tódolos pares de cebadores analizados. A maioría destas pertencían ós xéneros bacterianos: Streptococcus and Staphylococcus; e de arqueas: Methanosarcina, Thermocococcus, and Pyrococcus. Pola contra, 37 (19,89%) especies bacterianas e 27 (20,00%) de arqueas non tiñan ASI≥97% con outros taxons orais. Tódolos primers dirixidos a bacterias permitíronnos detectar 4450 relacións dous a dous entre 408 pares de especies bacterianas diferentes cun ASI97. Dezaoito destes pares de especies foron obtidos cos 29 pares de cebadores avaliados (frecuencia= 29; número de veces que un par de taxons teñen un ASI97 nos diferentes pares de cebadores), e pertencían ós xéneros Actinomyces, Lactobacillus, Neisseria, Staphylococcus, e Streptococcus. Sen embargo, 50 pares de especies con ASI97 só se detectaron cun primer (frecuencia= 1). Aínda que as relacións dous a dous implicaban maioritariamente a especies dos mesmos xéneros (3641; 81,82%), 809 relacións (18,18%) estaban constituídas por taxons de xéneros diferentes. A combinación de especies de Klebsiella con outras de Cronobacter foi a máis frecuente (frecuencia= 99), seguida de Klebsiella-Serratia, Escherichia-Klebsiella, Cronobacter-Escherichia e Aggregatibacter- ALBA REGUEIRA IGLESIAS 26 Haemophilus (frecuencias= 67 - 28). En canto ós rangos taxonómicos superiores, 293 (6,58%) relacións déronse entre pares de especies cun ASI97 pertencentes a familias distintas e mesmo houbo 26 (0,58%) relacións entre pares de especies cun ASI97 de ordes distintas. Os cebadores dirixidos a arqueas permitiron detectar 3232 relacións dous a dous entre 340 pares de especies de arqueas diferentes cun ASI97. Tódolos pares de cebadores analizados identificaron sete pares de especies (frecuencia= 20), que pertencían ós xéneros Methanobrevibacter e Methanocaldococcus. Houbo 66 pares de especies detectadas unha soa vez por un só par de cebadores (frecuencia= 1). Unha vez máis, a maioría das relacións de dous a dous foron entre especies arqueas do mesmo xénero (2359, 72,99%), pero 873 (27,01%) relacións implicaban pares de taxons cun ASI97 de xéneros distintos. A combinación de especies de Pyrococcus e Thermococcus foi a máis frecuente (frecuencia= 428), seguida de Palaeococcus e Thermococcus (frecuencia= 109). Para os rangos taxonómicos superiores, 35 (1,08%) relacións eran pares de especies cun ASI97 de familias distintas, tres (0,09%) de ordes, e unha de clases (0,03%) distintas. 3.3 CONCLUSIÓNS Os pares de cebadores avaliados dirixidos a bacterias e/ou arquexas detectaron unha media de máis de 150 OTU potenciais que poderían conter especies diferentes, cando se utilizou o limiar de similitude do ≥97%. Segundo o parámetro SC-NASI97, os mellores pares de cebadores foron OP_F053-KP_R020 para bacterias (rexión 1-3; posición para E. coli J01859.1: 9-356); KP_F018-KP_R002 para arqueas (4; indefinido-532); e OP_F114-KP_R031 para ambas (3-5; 340-801). Ó redor do 80% das especies de bacterias e arqueas orais analizadas tiñan un ASI97 con polo menos outra especie distinta. Estas especies altamente similares desempeñan distintas funcións na microbiota oral e pertencen a xéneros bacterianos como Campylobacter, Rothia, Streptococcus e Tannerella, e a xéneros de arqueas como Halovivax, Methanosalsum e Methanosarcina. Ademais, o 20% e o 30% das relacións de similitude dous a dous establecéronse entre especies de diferentes xéneros bacterianos e de arqueas, respectivamente. Mesmo taxons de familias, ordes e clases distintas puideron agruparse no mesmo OTU potencial. Resumo da tese 27 Independentemente do par de cebadores empregado, a agrupación de secuencias cunha similitude ≥97% proporciona unha descrición inexacta das especies orais bacterianas e de arqueas, o que pode afectar en gran medida ós parámetros de diversidade microbiana. Como resultado, a agrupación de OTUs condiciona a credibilidade das asociacións entre algunhas especies orais e certas condicións de saúde e enfermidade. Isto limita de xeito significativo a comparativa dos resultados da diversidade microbiana reportados na literatura do microbioma oral. ALBA REGUEIRA IGLESIAS 28 OBXECTIVO 4. Unha análise meta-ómica a gran escala da microbiota da placa nas enfermidades periodontais 4.1 MATERIAL E MÉTODOS Un total de 120 participantes: 55 periodontalmente sans e 65 afectados por periodontite non tratada; que cumprían os criterios de inclusión preestablecidos foron recrutados. Os diagnósticos periodontais foron realizados por dous dentistas experimentados. A presenza de saúde periodontal ou de periodontite crónica xeneralizada de moderada a grave estableceuse segundo información clínica e radiográfica, aplicando criterios previamente publicados (60,61). Unha ou dúas semanas despois do exame inicial, recolléronse mostras de placa subxinxival de tódolos participantes. O ADN total das mostras foi extraído e illado para avaliar a súa calidade e concentración. Neste punto, excluíronse dúas mostras subxinxivales do grupo san por non cumprir os requisitos de calidade e concentración. A continuación, realizouse unha amplificación por reacción en cadea da polimerasa (PCR) da rexión 3-4 do xene ARNr 16S (62) A secuenciación levouse a cabo na plataforma Illumina MiSeq con lecturas de 2x300 bps. As secuencias obtidas depositáronse no arquivo de lecturas de secuencias (sequence read archive, SRA) (63) co número de acceso PRJNA773202. Doutra banda, tamén se incluíron na nosa investigación estudos previos sobre a diversidade microbiana na placa supraxinxival e subxinxival en individuos adultos con diferentes condicións periodontais. Incorporáronse tódolos estudos que empregaron cebadores da rexión 3-4, a tecnoloxía de secuenciación Illumina, e tiñan acceso ao repositorio das secuencias. O estándar de referencia para o diagnóstico dunha afección periodontal podía basearse unicamente en parámetros clínicos ou clínicos e radiográficos, independentemente dos criterios de referencia diagnósticos aplicados; pero tiña que ser reportado. Ademais, seleccionáronse aquelas investigacións nas que os metadatos de interese por mostra estaban correctamente asignados no repositorio; e para os que as secuencias almacenadas cumprían criterios adicionais de inclusión e exclusión. En xullo de 2021 realizáronse 120 buscas en cada unha das bases de datos electrónicas PubMed, Scopus e Embase para identificar estudos de secuenciación mediante Illumina sobre Resumo da tese 29 o microbioma periodontal utilizando dous conxuntos de termos relacionados con: 1) condicións de saúde periodontal, nichos orais e microbiota; e 2) a tecnoloxía de secuenciación do xene ARNr 16S. Realizáronse procuras adicionais na base de datos SRA (63) para garantir que examinaramos tódolos posibles bioproxectos de interese. A manipulación dos datos identificados nas procuras realizouse utilizando o software R (32), e os resultados de cada unha almacenáronse individualmente nun arquivo txt. (PubMed) ou csv. (Scopus e Embase). Os duplicados foron detectados e eliminados. A continuación, os resumos analizáronse de xeito computacional mediante sete conxuntos de termos positivos utilizando os paquetes tm e tratamento da linguaxe natural (natural language processing, NLP) (34,64) e cada un recibiu 100 ou un punto por cada termo presente. As publicacións con puntuacións predefinidas seleccionáronse para as posteriores avaliacións manuais dos seus resumos e texto completo. O proceso automatizado de extracción de datos validouse previamente nunha serie limitada de artigos que cumprían os criterios de inclusión. Os identificadores dos bioproxectos dos artigos seleccionados empregáronse para acceder á base de datos do SRA (63) e ó selector de execucións “SRA run selector” (https://www.ncbi.nlm.nih.gov/traces/study/). Entón, descargáronse e avaliáronse as táboas de metadatos alí depositadas. Os autores foron contactados nos casos necesarios para obter os metadatos ou con fins de aclaración. Neste punto, o noso bioproxecto, PRJNA773202, engadiuse ó total. Coa información da base de datos SRA (63) construíuse manualmente unha táboa de metadatos para cada un dos bioproxectos, que incluía información sobre aspectos relacionados tanto co bioproxecto como coas variables demográficas e clínicas dos doentes ós que pertencía cada mostra. En relación coas secuencias almacenadas, para cada bioproxecto descargouse a lista de identificadores (listas de acceso) correspondente ás mostras de interese en formato txt. Para descargar e almacenar as secuencias coas mencionadas listas de acceso, instalouse o software gratuíto SRA Toolkit (65) en modo local. A continuación desenvolveuse un script en Bash (38) en combinación cos comandos prefetch e fastq-dump do SRA Toolkit. As mostras de placas de ALBA REGUEIRA IGLESIAS 30 cada bioproxecto almacenáronse en arquivos fastq individuais para a súa posterior manipulación. O preprocesamento e a avaliación da calidade das secuencias de cada arquivo fastq realizáronse con USEARCH (66). As secuencias foron aliñadas e ensambladas, aceptándose un máximo de cinco desaxustes e unha porcentaxe mínima de similitude do 90% para os 2x250 bps, e 10 bps e o 80% para os 2x300 bps. Permitíronse un máximo de dous desaxustes na secuencia de cada cebador individual e catro nun par. Por último, descartáronse as secuencias cun erro máximo esperado >1 ou cunha lonxitude mínima <300 bps. Tódolos arquivos fasta dun determinado bioproxecto fusionáronse usando Bash (38). Isto, mais un script desenvolvido en R (32), permitiunos crear un arquivo de grupo para cada bioproxecto. Executáronse os comandos screen. seqs e unique. seqs de mothur (67) para obter o arquivo de nomes de cada bioproxecto. A continuación, tódolos arquivos fasta de tódolos bioproxectos fusionáronse nun único arquivo, xunto cos arquivos "grupos" e "nomes" para formar os arquivos “meta-omics”. A continuación, realizouse un filtrado de baixa abundancia, no que se eliminaron as secuencias das mostras globais con valores de abundancia <500 recontos. Aplicouse o pipeline de mothur (67) para ASVs con lixeiras modificacións, incluíndo a aplicación da base de datos específica oral para a clasificación taxonómica de ASVs descrita por Escapa et al. (36). Permitíronse as secuencias cunha lonxitude >400; e elimináronse as que tiñan máis de oito homopolímeros, as que se consideraron quimeras tras aplicar o algoritmo VSEARCH de mothur (67,68) e as clasificadas como taxons descoñecidos no nivel xerárquico máis alto. As secuencias non se agruparon en ningún nivel xa que o noso obxectivo era identificar e clasificar o maior número posible de secuencias no nivel ASV. Unha vez completado o pipeline de mothur (67), exportáronse os seguintes arquivos metaómicos a R-bioconductor (69) para a súa posterior análise: a táboa de reconto, a xerarquía taxonómica a nivel ASV, a árbore filoxenética e a táboa de metadatos. Resumo da tese 31 Ademais, dous autores avaliaron de xeito independente a calidade dos metadatos dos bioproxectos incluídos na nosa investigación empregando unha lista de comprobación que deseñamos para este fin, a cal contiña 19 variables relacionadas cos datos dispoñibles sobre os suxeitos da mostra. Segundo a puntuación obtida, os bioproxectos clasificáronse como: baixa calidade= 0,00 - 0,33; calidade media= 0,34 - 0,66; e alta calidade= 0,67 - 1,00. Tamén se avaliou o número de mostras por bioproxecto e o número medio de secuencias por mostra en cada proxecto como parámetro de calidade. O número medio de secuencias dividiuse por 10.000 e este parámetro denominouse puntuación media da secuencia (average sequence score, ASS). Os valores de ASS interpretáronse como representativos de 1) <0,25, secuencias de moi baixa cantidade; 2) 0,25 - 0,75, - secuencias de baixa cantidade; 3) 0,75 - 1,0 - secuencias de cantidade aceptable; 4) 1,0 - 2,0 - secuencias de alta cantidade; e 5) >2,0 - secuencias de moi alta cantidade. Por último, a análise estatística dos datos de secuenciación do ARNr 16S a nivel ASV realizouse segundo o protocolo proposto por McMurdie e Holmes (70), utilizando implementacións en R que incluían os paquetes phyloseq, DESeq2 e microbiome (71-73). Para eliminar as mostras cun baixo número de secuencias, excluíronse aquelas con menos de 2500 (n= 62), polo que quedaron 2124 mostras. A continuación, creáronse grupos segundo o tipo de placa dental e o estado de saúde periodontal dos participantes, obténdose un total de 11 grupos (ordenados alfabeticamente): 1) Placa supraxinxival, saúde periodontal, sitios sans (Sup_x0HHx; n= 210). 2) Placa supraxinxival, xenxivite, sitios enfermos (Sup_x0GDx; n= 79). 3) Placa supraxinxival, periodontite, sitios enfermos (Sup_x0PDx; n= 493). 4) Placa supraxinxival, periodontite, sitios enfermos tratados (Sup_x1PDx; n= 81). 5) Placa subxinxival, saúde periodontal, sitios sans (Sub_x0HHx; n= 155). 6) Placa subxinxival, xenxivite, sitios enfermos (Sub_x0GDx; n= 20). 7) Placa subxinxival, periodontite, sitios sans (Sub_x0PHx; n= 62). 8) Placa subxinxival, periodontite, zonas enfermas (Sub_x0PDx; n= 768). 9) Placa subxinxival, periodontite, sitios enfermos tratados (Sub_x1PDx; n= 197). ALBA REGUEIRA IGLESIAS 32 10) Placa submucosa, peri-implantite, sitios sans (Imp_x0IHx; n= 18). 11) Placa submucosa, peri-implantite, sitios enfermos ( Imp_x0IDx; n= 41). Os grupos Sub_x0GDx, Imp_x0IHx e Imp_x0IDx elimináronse debido ó seu baixo tamaño de mostra (n= <50), deixando un total de 2045 mostras para analizar. A relación entre as diferentes condicións de saúde periodontal e a microbiota da placa investigouse dende varias perspectivas. En primeiro lugar, utilizáronse os paquetes phyloseq e microbiome para obter os datos de diversidade alfa (71,73). Como indicadores da riqueza de taxons, calculáronse o reconto absoluto de ASVs e o índice de cobertura do 95%. Determináronse os índices de Shannon e Pielou como indicadores de diversidade e uniformidade (74,75). Utilizouse a U de Mann-Whitney para as análises comparativas. Segundo, empregouse unha análise de compoñentes principais (principal component analysis, PCA) para visualizar a agrupación das mostras de placas en relación co seu estado de saúde. O paquete mixOmics (76) empregouse para obter os gráficos de dispersión das dúas compoñentes principais baseadas na abundancia relativa das ASVs, amosando os centroides de cada grupo clínico e as elipses que representan o intervalo de confianza do 96%. Utilizouse unha análise non paramétrica multivariante permutada da varianza (permutational multivariate analysis of variance, PERMANOVA) (77) para medir as diferenzas a nivel de comunidade entre os grupos. Estas análises realizáronse co paquete vegan (78). Terceiro, utilizouse o paquete microbiome (73) para identificar as ASVs centrais ou do núcleo presentes cunha taxa de prevalencia de ≥75% en cada tipo de placa e cada condición periodontal. En cuarto lugar, utilizouse o paquete DESeq2 (72) para identificar as ASV cos cambios máis significativos na abundancia diferencial para as distintas condicións periodontais. As abundancias diferenciais medíronse co valor log2foldchange (log2 FC), e as diferentes condicións comparáronse utilizando a proba de Wald coa corrección de Benjamini-Hochberg. As medidas foron estatisticamente significativas se o p axustado era <0,01. Resumo da tese 39 En canto ás preditoras da enfermidade, centrándonos nos valores de abundancia relativa observados en Sub_x0PDx, as ASVs máis relevantes foron Peptostreptococcaceae [XI][G-5] saphenum ASV129, Dialister pneumosintes ASV194, Desulfobulbus HMT041 ASV149 e Mogibacterium timidum ASV640. A primeira ASV tamén predixo sitios tratados e sans na periodontite ó comparar ambos modelos coa saúde periodontal, mentres que a segundo e a terceira tamén predixeron a periodontite tratada. Dos resultados derivados dos modelos preditivos sobre a placa subxinxival, F. nucleatum subsp. vicentii foi a principal especie que mostrou un rendemento preditivo oposto. A ASV do núcleo e altamente abundante ASV10 foi unha forte preditora da periodontite detectado en varios modelos. Con todo, houbo outras ASVs menos abundantes que predixeron simultaneamente tanto a saúde como a periodontite. En canto ás ASVs preditoras da saúde na placa “supra” e “sub”, destacan pola súa abundancia as seguintes: R. dentocariosa ASV2, Haemophilus parainfluenzae ASV3, ASV78, ASV45, e ASV46, Kingella oralis ASV66, Streptococcus vestibularis ASV27 e Actinomyces HMT170 ASV119. Algunhas destas ASVs tamén se comportaron como fortes discriminantes de sitios sans na periodontite. En canto ás preditoras da periodontite en ambas placas, atribúese especial atención polos seus valores de abundancia a: T. forsythia ASV15, Filifactor alocis ASV19, Treponema denticola ASV38 e ASV150, Fretibacterium fastidiosum ASV97, Peptostreptococcaceae [XI][G-4] HMT369 ASV124, Streptococcus anginosus ASV142, e Peptostreptococcaceae [XI][G-6] nodatum ASV189. Tamén, F. fastidiosum ASV97 e S. anginosus ASV142 foron preditoras de xenxivite na placa supraxinxival. Todos estes taxons predixeron a periodontite tratada na placa subxinxival, e T. forsythia ASV15 mesmo na supraxinxival. As ASVs dos xéneros Alloprevotella, Fusobacterium, Gemella, Granulicatella, Lachnoanaerobaculum e Ruminococcaceae [G-1] discriminaron diferentes condicións clínicas entre ambas placas, emerxendo como preditoras da xenxivite na placa supraxinxival e da saúde na subxinxival. Exemplos disto foron: Fusobacterium periodonticum ASV11, Lachnoanaerobaculum umeaense ASV152 e Granulicatella elegans ASV207. ALBA REGUEIRA IGLESIAS 40 Centrándonos na correspondencia dos resultados dos modelos preditivos das principais ASVs con respecto ós seus respectivos resultados na abundancia diferencial, detectáronse numerosas inconsistencias biolóxicas nesta última. Na placa supraxinxival, 17 das 21 ASVs preditoras de saúde presentaban abundancias diferenciais, e o 88,23% (15/17) delas amosaban abundancias significativamente elevadas tanto en saúde como en periodontite. Das 25 ASVs preditoras de enfermidade, 22 eran diferencialmente abundantes e, delas, catro (18,18%) tiñan una tendencia a niveis elevados nas dúas condicións opostas. Na placa subxinxival, 29 das 31 ASVs preditoras de saúde foron diferencialmente abundantes, cun 31,03% (9/29) amosando abundancias significativamente elevadas tanto en saúde como en periodontite. Das 31 ASVs preditoras de enfermidade, 30 amosaron abundancias diferenciais e, delas, tres (10,00%) presentaron niveis significativamente elevados nas dúas condicións opostas. En ámbolos dous tipos de placa, 28 das 32 ASVs preditoras de saúde amosaron abundancia diferencial e, delas, o 71,42% (20/28) foron diferencialmente abundantes tanto para a saúde como para a periodontite. As 15 ASVs preditoras de enfermidade en ambas placas tamén foron diferencialmente abundantes e, delas, tres (20,0%) tiñan niveis significativamente elevados en ámbalas dúas condicións clínicas. 4.3 CONCLUSIÓNS A riqueza bacteriana asociada á periodontite é maior que na saúde periodontal na placa supraxinxival e menor na subxinxival; a uniformidade é maior na enfermidade que na saúde en ambos os nichos. A microbiota supraxinxival é máis rica e diversa que a súa homóloga subxinxival para o mesmo estado de saúde periodontal. A estrutura da comunidade bacteriana é diferente para as distintas condicións periodontais na placa supraxinxival e subxinxival, así como para o mesmo estado de saúde entre os dous nichos. O núcleo da microbiota da placa supraxinxival e subxinxival non permite caracterizar a saúde e a enfermidade periodontal, o que revela a gran heteroxeneidade da microbiota oral. A porcentaxe da comunidade bacteriana da placa dental que se organiza en redes de co-ocorrencia a nivel de ASV é moi pequena; a rede da periodontite non tratada da placa supraxinxival é máis Resumo da tese 41 extensa e contén máis nodos, interconexións e grupos bacterianos interconectados que a súa homóloga subxinxival. As principais ASV clave nas redes de saúde periodontal da placa supraxinxival son R. dentocariosa ASV2 e S. oralis subsp. dentisani clade 058 ASV1. O eixo principal nas redes de periodontite non tratadas da placa supraxinxival é S. sanguinis ASV228; e na placa subxinxival, é T. forsythia ASV15. Os principais taxons clave na rede de periodontite tratada do nicho subxinxival son T. forsythia ASV15, F. nucleatum subsp. vincentii ASV10 e S. oralis subsp. dentisani clade 058 ASV1. Unha pequena proporción dos taxons supra e subxinxivales teñen unha capacidade destacada para distinguir entre as condicións periodontais, e unha porcentaxe relevante son membros da microbiota central. Dende o punto de vista da metaxenómica clínica, a placa supraxinxival é un mellor biomarcador bacteriano que o seu homólogo subxinxival para diferenciar a saúde periodontal da periodontite non tratada e tratada. As principais ASVs preditoras da saúde periodontal na placa supraxinxival e subxinxival son R. dentocariosa ASV2; H. parainfluenzae ASV3, ASV78, ASV45 e ASV46; K. oralis ASV66; S. vestibularis ASV27; e A. HMT170 ASV119. Pola contra, as principais ASVs preditoras da periodontite en ámbolos tipos de placa son: T. forsythia ASV15; F. alocis ASV19; T. denticola ASV38 e ASV150; F. fastidiosum ASV97; P. HMT369 ASV124; S. anginosus ASV142; e P. nodatum ASV189. Destas, F. fastidiosum ASV97 e S. anginosus ASV142 tamén actuaron como preditoras de xenxivite no nicho supraxinxival. ALBA REGUEIRA IGLESIAS 42 REFERENCIAS (1) Kinane DF, Stathopoulou PG, Papapanou PN. Periodontal diseases. Nat Rev Dis Primers. 2017 Jun;3:17038. doi: 10.1038/nrdp.2017.38. (2) Chapple ILC, Mealey BL, Van Dyke TE, Bartold PM, Dommisch H, Eickholz P, et al. Periodontal health and gingival diseases and conditions on an intact and a reduced periodontium: consensus report of workgroup 1 of the 2017 World Workshop on the Classification of Periodontal and Peri-Implant Diseases and Conditions. J Periodontol. 2018 Jun;89 Suppl 1:S74-84. doi: 10.1002/JPER.17-0719. (3) Peres MA, Macpherson LMD, Weyant RJ, Daly B, Venturelli R, Mathur MR, et al. Oral diseases: a global public health challenge. Lancet. 2019 Jul;394(10194):249-60. (4) Kassebaum NJ, Smith AGC, Bernabé E, Fleming TD, Reynolds AE, Vos T, et al. Global, regional, and national prevalence, incidence, and disability-adjusted life years for oral conditions for 195 countries, 1990–2015: a systematic analysis for the global burden of diseases, injuries, and risk factors. J Dent Res. 2017 Apr;96(4):380-7. (5) Bravo-Pérez M, Almerich-Silla J, Ausina-Márquez V, Avilés-Gutiérrez P, Blanco-González J, Canorea-Díaz E, et al. Oral health survey in Spain 2015. RCOE. 2016 Jun;21(Suppl 1):8-48. Spanish. (6) Ferreira MC, Dias-Pereira A, Branco-de-Almeida LS, Martins CC, Paiva SM. Impact of periodontal disease on quality of life: a systematic review. J Periodont Res. 2017 Aug;52(4):651-65. (7) Chapple ILC, Bouchard P, Cagetti MG, Campus G, Carra M, Cocco F, et al. Interaction of lifestyle, behaviour or systemic diseases with dental caries and periodontal diseases: consensus report of group 2 of the joint EFP/ORCA workshop on the boundaries between caries and periodontal diseases. J Clin Periodontol. 2017 Mar;44 Suppl 18:S39-51. doi: 10.1111/jcpe.12685. Resumo da tese 43 (8) Carrizales-Sepúlveda EF, Ordaz-Farías A, Vera-Pineda R, Flores-Ramírez R. Periodontal disease, systemic inflammation and the risk of cardiovascular disease. Heart Lung Circ. 2018 Nov;27(11):1327-34. (9) Nascimento G, Leite F, Vestergaard P, Scheutz F, López R. Does diabetes increase the risk of periodontitis? A systematic review and meta-regression analysis of longitudinal prospective studies. Acta Diabetol. 2018 Jul;55(7):653-67. (10) Gomes-Filho I, Cruz SSD, Trindade SC, Passos-Soares J, Carvalho-Filho P, Figueiredo ACMG, et al. Periodontitis and respiratory diseases: a systematic review with meta-analysis. Oral Dis. 2020 Mar;26(2):439-46. (11) Fuggle NR, Smith TO, Kaul A, Sofat N. Hand to mouth: a systematic review and metaanalysis of the association between rheumatoid arthritis and periodontitis. Front Immunol. 2016 Mar;7:80. doi: 10.3389/fimmu.2016.00080. (12) Leira Y, Domínguez C, Seoane J, Seoane-Romero J, Pías-Peleteiro JM, Takkouche B, et al. Is periodontal disease associated with alzheimer's disease? A systematic review with metaanalysis. Neuroepidemiology. 2017;48(1-2):21-31. (13) Ide M, Papapanou PN. Epidemiology of association between maternal periodontal disease and adverse pregnancy outcomes – systematic review. J Clin Periodontol. 2013 Apr;40 Suppl 14:S181-94. doi: 10.1111/jcpe.12063. (14) AlJehani YA. Risk factors of periodontal disease: review of the literature. Int J Dent. 2021 Feb;2021:8735071. doi: 10.1155/2021/8735071. (15) Teles R, Teles F, Frias-Lopez J, Paster B, Haffajee A. Lessons learned and unlearned in periodontal microbiology. Periodontol 2000. 2013 Jun;62(1):95-162. (16) A genomic approach to microbiology. Nat Rev Genet. 2019 Jun;20(6):311. doi: 10.1038/s41576-019-0131-5. ALBA REGUEIRA IGLESIAS 44 (17) del Rosario-Rodicio M, del Carmen-Mendoza M. Bacterial identification by 16S rRNA sequencing: rationale, methodology and applications in clinical microbiology. Enferm Infecc Microbiol Clin. 2004 Apr;22(4):238-45. Spanish. (18) Zaura E. Next-generation sequencing approaches to understanding the oral microbiome. Adv Dent Res. 2012 Sep;24(2):81-5. (19) Nearing JT, Comeau AM, Langille MGI. Identifying biases and their potential solutions in human microbiome studies. Microbiome. 2021 May;9(1):113. doi: 10.1186/s40168-021- 01059-0. (20) Robinson CK, Brotman RM, Ravel J. Intricacies of assessing the human microbiome in epidemiologic studies. Ann Epidemiol. 2016 May;26(5):311-21. (21) de la Cuesta-Zuluaga J, Escobar JS. Considerations for optimizing microbiome analysis using a marker gene. Front Nutr. 2016 Aug;3:26. doi: 10.3389/fnut.2016.00026. (22) Hamady M, Knight R. Microbial community profiling for human microbiome projects: tools, techniques, and challenges. Genome Res. 2009 Jul;19(7):1141-52. (23) Rajendhran J, Gunasekaran P. Microbial phylogeny and diversity: small subunit ribosomal RNA sequence analysis and beyond. Microbiol Res. 2011 Feb;166(2):99-110. (24) Acinas SG, Marcelino LA, Klepac-Ceraj V, Polz MF. Divergence and redundancy of 16S rRNA sequences in genomes with multiple rrn operons. J Bacteriol. 2004 May;186(9):2629- 35. (25) Pei AY, Oberdorf WE, Nossa CW, Agarwal A, Chokshi P, Gerz EA, et al. Diversity of 16S rRNA genes within individual prokaryotic genomes. Appl Environ Microbiol. 2010 Jun;76(12):3886-97. Resumo da tese 45 (26) Sun D, Jiang X, Wu QL, Zhou N. Intragenomic heterogeneity of 16S rRNA genes causes overestimation of prokaryotic diversity. Appl Environ Microbiol. 2013 Oct;79(19):5962-9. (27) Větrovský T, Baldrian P. The variability of the 16S rRNA gene in bacterial genomes and its consequences for bacterial community analyses. PLoS One. 2013;8(2):e57923. doi: 10.1371/journal.pone.0057923. (28) Case RJ, Boucher Y, Dahllöf I, Holmström C, Doolittle WF, Kjelleberg S. Use of 16S rRNA and rpoB genes as molecular markers for microbial ecology studies. Appl Environ Microbiol. 2007 Jan;73(1):278-88. (29) Lee ZM, Bussema C 3rd, Schmidt TM. rrnDB: documenting the number of rRNA and tRNA genes in bacteria and archaea. Nucleic Acids Res. 2009 Jan;37(Database issue):D489- 93. doi: 10.1093/nar/gkn689. (30) Johnson JS, Spakowicz DJ, Hong BY, Petersen LM, Demkowicz P, Chen L, et al. Evaluation of 16S rRNA gene sequencing for species and strain-level microbiome analysis. Nat Commun. 2019 Nov;10(1):5029. doi: 10.1038/s41467-019-13036-1. (31) Schloss PD. Amplicon sequence variants artificially split bacterial genomes into separate clusters. mSphere. 2021 Aug;6(4):e0019121. doi: 10.1128/mSphere.00191-21. (32) R Core Team. R: a language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing; 2021; Available at: https://www.R-project.org/. (33) Kovalchik S. RISmed: download content from NCBI databases. 2017; Available at: http://www.CRAN.R-project.org/. (34) Feinerer, I., Hornik, K., Meyer, D. Text mining infrastructure in R. J Stat Softw. 2008 Mar;25(5):1-54. doi: 10.18637/jss.v025.i05. ALBA REGUEIRA IGLESIAS 46 (35) Klindworth A, Pruesse E, Schweer T, Peplies J, Quast C, Horn M, et al. Evaluation of general 16S ribosomal RNA gene PCR primers for classical and next-generation sequencingbased diversity studies. Nucleic Acids Res. 2013 Jan;41(1):e1. doi: 10.1093/nar/gks808. (36) F Escapa I, Huang Y, Chen T, Lin M, Kokaras A, Dewhirst FE, et al. Construction of habitat-specific training sets to achieve species-level assignment in 16S rRNA gene datasets. Microbiome. 2020 May;8(1):65. doi: 10.1186/s40168-020-00841-w. (37) Python Software Foundation. Python. 2020; Available at: http://www.python.org/. (38) GNU P. Free Software Foundation. Bash. 2020; Available at: http://www.gnu.org/. (39) Sievers F, Wilm A, Dineen D, Gibson TJ, Karplus K, Li W, et al. Fast, scalable generation of high-quality protein multiple sequence alignments using Clustal Omega. Mol Syst Biol. 2011 Oct;7:539. doi: 10.1038/msb.2011.75. (40) NCBI Resource Coordinators. Database resources of the National Center for Biotechnology Information. Nucleic Acids Res. 2016 Jan;44(D1):D7-19. doi: 10.1093/nar/gkv1290. (41) O'Leary NA, Wright MW, Brister JR, Ciufo S, Haddad D, McVeigh R, et al. Reference sequence (RefSeq) database at NCBI: current status, taxonomic expansion, and functional annotation. Nucleic Acids Res. 2016 Jan;44(D1): D733-45. doi: 10.1093/nar/gkv1189. (42) Clark K, Karsch-Mizrachi I, Lipman DJ, Ostell J, Sayers EW. GenBank. Nucleic Acids Res. 2016 Jan;44:D67-72. doi: 10.1093/nar/gkv1276. (43) Lyalina S. Search 16S py algorithm. 2019; Available at: https://github.com/slyalina/search_16S_py. (44) Edgar R. SEARCH_16S: a new algorithm for identifying 16S ribosomal RNA genes in contigs and chromosomes. Preprint at bioRxiv 2017:124131. doi: 10.1101/124131. Resumo da tese 47 (45) National Center for Biotechnology Information. Entrez programming utilities help. 2010; Available at: https://www.ncbi.nlm.nih.gov/books/NBK25501/. (46) Cock PJ, Antao T, Chang JT, Chapman BA, Cox CJ, Dalke A, et al. Biopython: freely available Python tools for computational molecular biology and bioinformatics. Bioinformatics. 2009 Jun;25(11):1422-23. (47) National Center for Biotechnology Information. NCBI RefSeq targeted loci project. Archaea FTP. 2008. ftp://ftp.ncbi.nlm.nih.gov/refseq/TargetedLoci/Archaea/. (48) Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, et al. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res. 2013 Jan;41(Database issue):D590-6. doi: 10.1093/nar/gks1219. (49) Parks DH, Chuvochina M, Chaumeil P, Rinke C, Mussig AJ, Hugenholtz P. A complete domain-to-species taxonomy for bacteria and archaea. Nat Biotechnol. 2020 Sep;38(9):1079- 86. (50) Morgulis A, Coulouris G, Raytselis Y, Madden TL, Agarwala R, Schäffer AA. Database indexing for production MegaBLAST searches. Bioinformatics. 2008 Aug;24(16):1757-64. (51) Chen Y, Ye W, Zhang Y, Xu Y. High speed BLASTN: an accelerated MegaBLAST search tool. Nucleic Acids Res. 2015 Sep;43(16):7762-8. (52) Barnett M. regex. 2020; Available at: https://pypi.org/. (53) F Escapa I, Chen T, Huang Y, Gajare P, Dewhirst FE, Lemon KP. New insights into human nostril microbiome from the expanded Human Oral Microbiome Database (eHOMD): a resource for the microbiome of the human aerodigestive tract. mSystems. 2018 Dec;3(6):e00187-18. doi: 10.1128/mSystems.00187-18. ALBA REGUEIRA IGLESIAS 48 (54) Schoch CL, Ciufo S, Domrachev M, Hotton CL, Kannan S, Khovanskaya R, et al. NCBI Taxonomy: a comprehensive update on curation, resources and tools. Database (Oxford). 2020 Jan;2020:baaa062. doi: 10.1093/database/baaa062. (55) Harris CR, Millman KJ, van der Walt, Stéfan J, Gommers R, Virtanen P, Cournapeau D, et al. Array programming with NumPy. Nature. 2020 Sep;585(7825):357-62. (56) McKinney W. Data Structures for statistical computing in Python. In: van der Walt S, Millman J, editors. Proceedings of the 9th Python in Science Conference; 2010; Austin. Texas: SciPy; 2010. doi: 10.25080/Majora-92bf1922-00a. (57) Camacho C, Coulouris G, Avagyan V, Ma N, Papadopoulos J, Bealer K, et al. BLAST+: architecture and applications. BMC Bioinformatics. 2009 Dec;10:421. doi: 10.1186/1471- 2105-10-421. (58) Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990 Oct;215(3):403-10. (59) McNamara J. xlsxwriter. 2013; Available at: https://xlsxwriter.readthedocs.io/. (60) Armitage GC. Development of a classification system for periodontal diseases and conditions. Ann Periodontol. 1999 Dec;4(1):1-6. (61) Page RC, Eke PI. Case definitions for use in population-based surveillance of periodontitis. J Periodontol. 2007 Jul;78(7 Suppl):1387-99. (62) Willis JR, González-Torres P, Pittis AA, Bejarano LA, Cozzuto L, Andreu-Somavilla N, et al. Citizen science charts two major “stomatotypes” in the oral microbiome of adolescents and reveals links with habits and drinking water composition. Microbiome. 2018 Dec;6(1):218. doi: 10.1186/s40168-018-0592-3. 55 Thesis summary “Limitations of 16S rRNA gene as phylogenetic marker: a largescale meta-omics analysis of plaque microbiota in periodontal diseases” The first objective of the present Thesis was to analyse in silico the coverage of the 16S ribosomal RNA (rRNA) gene primers used to study the composition of the oral microbiota, using two databases containing 16S rRNA sequences from oral bacteria and archaea; and to describe the best primer pairs for each domain. Searches were conducted in PubMed to create a list of 1) 16S rRNA gene primers used in sequencing-based studies of the oral microbiome, and 2) oral-archaea species inhabiting the human mouth. The individual primers found were evaluated against a previously reported database of 16S rRNA sequences from oral bacteria, which was modified by our group; and a self-created oral-archaea database, constructed based on the list of oral-archaeal species. Both databases contained the genomic variants detected for each included species. Primers were evaluated at the variant and species levels, and those with a species coverage (SC) ≥75.00% were selected for the pair analyses. All possible combinations of forward and reverse primers were identified and evaluated against the two databases. A total of 369 distinct individual primers were found in the literature. After applying the primer-pair formation criteria, 4638 primer pairs were identified. The best bacteria-specific pairs targeted the 3-4, 4-7, and 3-7 16S rRNA gene regions, with SC levels of 98.83% - 97.14%; meanwhile, the optimum archaea-specific primer pairs amplified regions 5-6, 3-6, and 3-6, with SC estimates of 95.88%. Finally, the best pairs for detecting both domains targeted regions 4- 5, 3-5, and 5-9, and produced SC values of 95.71% - 94.54% and 99.48% - 96.91% for bacteria and archaea, respectively. ALBA REGUEIRA IGLESIAS 56 Given the three amplicon length categories (100-300, 301-600, and >600 base pairs), the primer pairs with the best coverage values for detecting oral bacteria were: KP_F048-OP_R043 (region 3-4; primer pair position for Escherichia coli J01859.1: 342-529), KP_F051-OP_R030 (4-7; 514-1079), and KP_F048-OP_R030 (3-7; 342-1079). For detecting oral archaea, these were: OP_F066-KP_R013 (5-6; 784-undefined), KP_F020-KP_R013 (3-6; 518-undefined), and OP_F114-KP_R013 (3-6; 340-undefined). Lastly, for detecting both domains they were: KP_F020-KP_R032 (4-5; 518-801), OP_F114-KP_R031 (3-5; 340-801), and OP_F066- OP_R121 (5-9; 784-1405). The primer pairs with the best coverage identified herein are not among those described most widely in the oral microbiome literature. The purpose of the second study was to evaluate the number of 16S rRNA genes in the complete genomes of all the bacterial and archaeal species ever detected in the human oral cavity; and to assess how the use of different primer pairs would affect the detection and classification of redundant amplicons and matching amplicons (MAs) from different taxa. A total of 709 complete genomes (518 oral bacteria, 191 oral archaea) were downloaded from the NCBI database, and their complete 16S rRNA genes were extracted. The total number of genes and variants per genome were calculated. Next, 33 primer pairs selected from objective 1 and 6 commonly employed in the oral literature were used against all the genomes to obtain amplicons. For each primer pair, we calculated the number of 16S rRNA gene amplicons, variants, genomes, and species detected, as well as the percentage of coverage at the species level with no matching amplicons (SC-NMA). In total, 94.1% of oral bacteria and 52.59% of oral archaea had more than one 16S rRNA gene in their respective genomes. Between 46.70% - 1.29% of the bacterial species and between 38.89% - 4.65% of the archaeal species detected by the evaluated primer pairs had MAs, affecting relevant genera present in the oral environment such as Actinomyces, Fusobacterium, Lactobacillus, Methanosarcina, Staphylococcus, and Streptococcus. The best primer pairs were (SC-NMA; region; primer pair position for Escherichia coli J01859.1): KP_F048-OP_R030 for bacteria (93.55%; 3-7; 342-1079), KP_F018-KP_R063 for archaea (89.63%; 3-9; undefined- 1506), and OP_F114_OP_R121 for both bacteria and archaea (92.52%; 3-9; 340-1405). Thesis summary 57 In addition to the 16S rRNA gene redundancy, the considerable presence of MAs must be controlled to ensure the accurate interpretation of microbial diversity data. The SC-NMA is a more useful parameter than the conventional coverage percentage for selecting the best primer pairs. The performance of the primer pairs to detect non-MA species increases as the average length of the amplicons increases; none of these being the most widely used primer pairs in the oral literature. The choice of primer pair significantly affects diversity estimates and taxonomic classification, conditioning the comparability of oral microbiome studies using different primer pairs. The aims of the third study were to evaluate in silico the coverage of a set of previously selected primer pairs to detect oral species having 16S rRNA sequence segments with ≥97% similarity; and to describe oral species with highly similar sequence segments and determine whether they belong to distinct genera or other higher taxonomic ranks. Thirty-nine primer pairs were employed to obtain the in silico amplicons from the complete genomes of 186 bacterial and 135 archaeal species. Each fasta file for the same primer pair was inserted as subject and query in BLASTN for obtaining the similarity percentage between amplicons belonging to different oral species. Amplicons with 100% alignment coverage of the query sequences and with a similarity value ≥97% (ASI97) were selected. For each primer, the species coverage with no ASI97 (SC-NASI97) was calculated. Based on the SC-NASI97 parameter, the best primer pairs were OP_F053-KP_R020 for bacteria (region 1-3; primer pair position for Escherichia coli J01859.1: 9-356); KP_F018- KP_R002 for archaea (4; undefined-532); and OP_F114-KP_R031 for both (3-5; 340-801). Around 80% of the oral-bacteria and oral-archaea species analysed had an ASI97 with at least one other species. These very similar species play different roles in the oral microbiota and belong to bacterial genera such as Campylobacter, Rothia, Streptococcus, and Tannerella, and archaeal genera such as Halovivax, Methanosalsum, and Methanosarcina. Moreover, ~20% and ~30% of these two-by-two similarity relationships were established between species from different bacterial and archaeal genera, respectively. Even taxa from distinct families, orders, and classes could be grouped in the same possible operational taxonomic unit (OTU). ALBA REGUEIRA IGLESIAS 58 Regardless of the primer pair used, sequence clustering with a 97% similarity provides an inaccurate description of oral-bacterial and oral-archaeal species, which can greatly affect microbial diversity parameters. As a result, OTU clustering conditions the credibility of associations between some oral species and certain health and disease conditions. This significantly limits the comparability of the microbial diversity findings reported in oral microbiome literature. Lastly, the fourth objective was to analyse the supragingival and subgingival plaque microbiota at amplicon sequence variant (ASV) level of different periodontal conditions (periodontal health, gingivitis, and untreated and treated periodontitis) in terms of bacterial diversity, co-occurrence networks, and predictive models. A total of 120 patients (55 controls, 65 periodontitis) were selected for subgingival plaque collection. Sequencing of the 3-4 16S rRNA gene region was performed in Illumina MiSeq. The obtained sequences and metadata were uploaded to the sequence read archive (SRA). Searches were performed in PubMed, Scopus, Embase, and the SRA to identify previously published Illumina 3-4 sequencing studies on the supragingival and subgingival plaque microbiome in distinct periodontal conditions. Research that met the criteria for sequences and metadata were included in the meta-omics analysis, comprising a total of 2045 samples. Sequences were processed under the same bioinformatics protocol, which included the ASV- level classification and the use of an oral-specific database for taxonomic classification. The statistical analysis was conducted using the phyloseq, DESeq2, microbiome, mixOmics vegan, SpiecEasi, and igraph packages. Bacterial richness associated with periodontitis was higher than in health in supragingival plaque and lower in subgingival, but evenness was higher in disease in both niches. The supragingival microbiota was richer and more diverse than the subgingival for the same periodontal condition. The structure of the bacterial community differed among conditions in the supra- and subgingival plaque, as well as for the same health status between the two niches. In addition, the core microbiota of dental plaque did not allow the characterisation of periodontal health and disease; and the proportion of the bacterial community organised in cooccurrence networks at the ASV level was very small. However, a small proportion of supra- Thesis summary 59 and subgingival taxa had outstanding ability to distinguish between periodontal conditions, and a relevant percentage of them were core members. Supragingival plaque was a better bacterial biomarker than subgingival for discriminating periodontal health from untreated and treated periodontitis. The main health-predictor ASVs in supragingival and subgingival plaque were: Rothia dentocariosa ASV2, Haemophilus parainfluenzae ASV3, ASV78, ASV45, and ASV46, Kingella oralis ASV66, Streptococcus vestibularis ASV27, and Actinomyces HMT170 ASV119. The main predictor ASVs of periodontitis in dental plaque were: Tannerella forsythia ASV15, Filifactor alocis ASV19, Treponema denticola ASV38 and ASV150, Fretibacterium fastidiosum ASV97, Peptostreptococcaceae [XI][G-4] HMT369 ASV124, Streptococcus anginosus ASV142, and Peptostreptococcaceae [XI][G-6] nodatum ASV189. INTRODUCTION 63 Introduction I.1. PERIODONTITIS: EPIDEMIOLOGY, DIAGNOSIS AND CLASSIFICATION The term “periodontal diseases” refers to several different chronic inflammatory conditions that affect the tissue surrounding and supporting the teeth (1). Pathologies arise when the balance between the microbial biofilm and the immune system is lost, whether due to dysbiosis (biofilm imbalance) or an overreaction of the host to the microbes present (1). This is a matter of concern since the 2016 estimates of the World Health Organisation (WHO) report that oral conditions, including periodontal diseases, are the 10th cause of years of healthy life lost due to disability (YLDs) globally (2). Gingivitis is a localised inflammation of the gums that originates from the bacteria present in the dental plaque deposited on the teeth and the gingiva (1). The condition manifests clinically with swelling, redness, bleeding on probing (BOP), and discomfort during the probing process. Patients also usually have symptoms like bleeding and swollen red gums, pain, halitosis, and difficulties when eating (3). Gingivitis does not, however, extend to the periodontal attachment apparatus (cementum, periodontal ligament, and alveolar bone) and is reversible by reducing plaque levels (3). In susceptible individuals, untreated gingivitis can progress to periodontitis (1), which is characterised by the gradual destruction of the tooth-supporting apparatus. Considerable damage to both the connective tissue fibres and the apical extension of the junctional epithelium in response to the accumulation of plaque is evident in advanced periodontal lesions. The bone destruction produced at this advanced stage creates the periodontal pockets that are the hallmark of the disease (1). In a clinical exploration, periodontitis manifests with redness, a changed texture and swelling of the margin of the gums, BOP, the increased depth of the periodontal pockets, the destruction of the ligament and alveolar bone, the recession of the marginal gingiva, increased tooth mobility, and, eventually, tooth loss (4). ALBA REGUEIRA IGLESIAS 64 Figure 1. Schematic representation of healthy gingiva, gingivitis, early-to-moderate periodontitis, and advanced periodontitis. The image was taken from Kinane et al. (1) with the permission of Springer Nature. I.1.1. Epidemiology Periodontal diseases are a significant public health concern, being among the most prevalent and consequential oral conditions worldwide (5). In 2015, the global prevalence of severe periodontitis was estimated to be 7.4%, with 538 million cases (6). In the US, 42.2% of dentate adults ≥30 years old had some category of periodontitis, with 7.8% affected by the severe disease (7). In those ≥65 years, these figures increased to 68% and 11%, respectively (8). According to data from the 2015 Oral Health Survey, 5% of the adult population in Spain aged from 34 to 44 years old and 10% of those between 65 and 74 had deep periodontal pockets (≥6 mm) (9). The prevalences of moderate pockets (4-5 mm) in these groups were 18.5% and 27.0%, respectively. Although these outcomes had not changed substantially from the previous survey in 2010 (10), there was an increase in the presence of moderate pockets, which may signify a likely worsening periodontal status in these age groups in subsequent years (9). I.1.2. Impact on quality of life Periodontitis has been associated with a negative effect on the quality of life, especially in patients with severe periodontitis, compromising aspects related to both function and aesthetics (11). Specifically, if no treatment is provided, outcomes can include tooth loss, impaired masticatory performance, a poorer nutritional status, lower self-esteem and quality of life, and negative effects on general health (12). The disease may also be a source of social inequality (4). Introduction 71 I.2. PERIODONTITIS AND ITS IMPLICATIONS FOR GENERAL HEALTH Several systemic diseases and conditions, both inherent and acquired, can affect the periodontal attachment apparatus, causing the loss of periodontal tissue (32). This damage can either 1) influence the course of periodontitis, or 2) affect periodontal-supportive tissue, irrespective of dental plaque biofilm-induced inflammation (18). The first case includes both rare disorders, in which periodontitis is a manifestation of the systemic condition itself (e.g., genetic disorders), as well as more common diseases (e.g., diabetes mellitus). The damage in the second case arises from very rare conditions, many of which are neoplasms (32). In this section of the thesis, however, the focus is on the relationship between periodontitis and common systemic pathologies. As long ago as 2000, Williams and Offenbacher (37) used the term “Periodontal Medicine” to define a then rapidly emerging branch of periodontology based on data that established a strong relationship between periodontal and systemic health or disease. Today, there is a convincing body of scientific evidence that is supportive of this notion of a two-way relationship between periodontitis and cardiovascular diseases (38), diabetes mellitus (39), respiratory conditions (40), rheumatoid arthritis (41), Alzheimer’s disease (42) and adverse pregnancy outcomes (43). Furthermore, recent years have seen an increase in research linking periodontitis to other disorders such as metabolic syndrome (44), obesity (45), chronic kidney disease (46), and orodigestive cancers, including those of the oral cavity, gastrointestinal tract, and pancreas (47). ALBA REGUEIRA IGLESIAS 72 Figure 2. Representation of the systemic conditions associated with periodontitis. From left to right: cardiovascular diseases, respiratory disorders, diabetes, Alzheimer’s disease, rheumatoid arthritis, and adverse pregnancy outcomes. Three underlying mechanisms through which periodontitis may play a role in general health have been hypothesised (48,49):  Metastatic infection - an infectious disease caused by microorganisms from a distant part of the body (48). Microbes can spread in three different ways: 1) direct propagation through contiguous spaces or venous or lymphatic drainage, 2) aspiration, or 3) bacteraemia or the access of oral bacteria to the bloodstream.  Inflammation and inflammatory injury - the indirect damage caused to tissue and organs by microbes via the dissemination of bacterial exotoxins and endotoxins.  Adaptative immunity - the immune response of the host to oral microorganisms and their virulence factors. Soluble antigens can enter the bloodstream, bond to a specific circulating antibody, and form a macromolecular immunocomplex, with the latter potentially leading to multiple acute and chronic inflammatory reactions at the deposition sites (49). Introduction 73 I.3. AETIOLOGY OF PERIODONTITIS: MICROBIOTA AND HOST RESPONSE As explained previously, both the initiation and progression of periodontitis are related to multiple aetiologies and risk factors, with the interaction between local microorganisms and the host’s immune response being particularly relevant. Consequently, the development of effective therapeutic approaches requires the identification of the main periodontitis-associated microbes and host biomarkers. I.3.1. Microbiota The human body is an ecosystem formed not only by human eukaryotic cells but also by a tremendous diversity of bacteria, archaea, fungi, and viruses. Such populations of microbes colonise the gastrointestinal and genitourinary tracts, the oral cavity, the nasopharynx, the respiratory tract, and the skin (50). Indeed, highlighting their importance, microorganisms are responsible for more than 200 grams of the total body weight of an average human (70 kilograms), while approximately 3.8 x 103 of the body’s cells are contributed by bacteria and other microorganisms (51). The oral cavity in particular has a high abundance of microbes, exceeded only by the numbers present in the gastrointestinal tract (52). Two imprecise terms are usually employed to designate this group of microorganisms: the “microbiota” or the “microbiome”. However, there is a critical difference between them in that the former includes the population of microorganisms that colonises a body part, while the latter refers to the setting formed by microbes, their genes, and their metabolites in an ecological niche (53). As the focus of this Thesis is on the detection and identification of bacteria in different oral niches, the term microbiota is used throughout. Our knowledge of the aetiology and pathogenesis of periodontal diseases has changed over time for four different reasons: 1) technological advances in the methods used to study the microbes present; 2) the current assumption that these conditions are caused by biofilms, and not by bacteria in a planktonic state, and the adoption of ecological concepts for studying the oral microbiota; 3) the discovery of the impact of genetic and environmental factors on the initiation and progression of these diseases; and 4) our understanding of the role played by immune mechanisms (54). As a consequence, several microbial theories have been proposed ALBA REGUEIRA IGLESIAS 74 for the aetiology of periodontitis, ranging from the “specific plaque hypothesis” to the “polymicrobial synergy and dysbiosis model” (PSD model) (54,55). The study of the composition of the oral microbiota dates back to 1683, when Antony van Leeuwenhoek observed the microorganisms present in human dental plaque using the first prototype of a microscope. Almost two centuries later, Robert Koch, the father of modern microbiology, developed techniques for producing bacterial cultures, which allowed him to scrutinise any changes in the bacteria present over time (56). Since then, as in other microbiology disciplines, oral bacteria have been detected using culture-dependent methods (57). The first studies of the composition of dental plaque had used such techniques to identify important organisms like Fusobacterium spp., Neisseria spp., Streptococcus spp., and Veillonella spp. (58). Nonetheless, there are several issues with these traditional methods when it comes to cultivating species that require rigorous growth conditions. In fact, only 50% of oral bacteria are cultivable (59). Additionally, culture-dependent techniques are expensive and laborious, requiring experienced staff and sufficient time for their execution (57). The goal of overcoming the latter issue led to the development of several molecular deoxyribonucleic acid (DNA)-based technologies, including DNA microarrays and the polymerase chain reaction (PCR) test. These enabled researchers to analyse oral-microbiota communities more comprehensively and perform large-scale studies (57,58). In 1998, Socransky et al. (60) used checker-board DNA-DNA-hybridisation techniques to identify five different bacterial complexes that had distinct levels of association with health and the severity of periodontitis. This discovery was revolutionary because, until then, periodontitis, like other infectious diseases, was thought to be caused by a single pathogen rather than a series of organisms working with each other (58). Three species in particular - P. gingivalis, Tannerella forsythia, and Treponema denticola - were closely associated with the clinical parameters for periodontitis and together constitute the so-called “Red Complex”. This complex has also been related to other bacteria, including Campylobacter gracilis, Campylobacter rectus, Campylobacter showae, Fusobacterium nucleatum, Fusobacterium periodonticum, Peptostreptococcus micros, Prevotella intermedia, Prevotella nigrescens, and Streptococcus constellatus, which collectively form the “Orange Complex”. In contrast, members of the “Yellow Complex” (Streptococcus gordonii, Streptococcus intermedius, Streptococcus mitis, Introduction 75 Streptococcus oralis, and Streptococcus sanguis) and the “Purple Complex” (Actinomyces odontolyticus and Veillonella parvula) are associated with healthy states. Figure 3. Socransky’s microbial complexes in subgingival plaque. A PCR test is an in vitro molecular technique that amplifies a gene or DNA fragment directly or a ribonucleic acid (RNA) indirectly. A variant of the PCR, known as a quantitative PCR (qPCR) or real-time PCR (RT-PCR), makes it possible to amplify, and simultaneously quantify, the amplification product obtained from a sample. In comparison to traditional cultures, this technology has been found to have high diagnostic accuracy when it comes to detecting A. actinomycetemcomitans and P. gingivalis (61). Furthermore, a qPCR enabled our research group to obtain eight bacterial cluster-based models with good predictive accuracy at identifying a site with periodontal destruction in a periodontitis patient (62). All of the models used by our team had an area under the curve (AUC) of ≥0.760 and sensitivity and specificity scores of ≥75.0%, with the best values for the cluster formed by A. actinomycetemcomitans, F. nucleatum, Parvimonas micra, P. intermedia, T. forsythia and T. denticola (AUC= 0.789; sensitivity and specificity= 77.5%). Overall, we concluded that clusters formed by species that had different etiopathogenic roles, i.e., those belonging to distinct Socransky complexes, had good predictive accuracy for diagnosing periodontitis (62). ALBA REGUEIRA IGLESIAS 76 The human oral microbe identification microarray (HOMIM) is another microarray-based platform that has been utilised to detect and identify both cultivated and not-yet-cultivated oral bacteria (58). Employment of this tool has enabled species like Filifactor alocis and P. micra to be observed more frequently and in higher numbers in periodontitis samples than in healthy specimens (63). However, this technique is currently no longer available (http://homings.forsyth.org/index2.html). These DNA microarray methods (checkerboard-DNA, DNA-hybridisation, and HOMIM) are not without their limitations, since only a fixed number of species can be detected in a panel and a specific quantity of DNA is required to identify a microorganism (58). The sequence analysis of the 16S ribosomal RNA (rRNA) bacterial gene became the method of choice to overcome these shortcomings. The amplification of this gene, which is present in all prokaryotic organisms, was achieved through the use of universal primers followed by a subsequent sequencing step. This enabled the species present in a sample to be distinguished, even if they had not been identified previously (58). The characteristics of this gene will be explained in detail in Section I.4. Initial research on sequence analyses of the 16S rRNA gene was based on the Sanger method, which is one of the so-called “first-generation sequencing” technologies (64). This technique used universal primers for the 16S rRNA gene to amplify the DNA isolated from a specimen. The resulting amplicons were cloned into Escherichia coli, and the inserts obtained were subsequently sequenced to determine the identities of the species present (56). New periodontitis-associated genera were discovered with this technology, including Desulfobulbus spp., Eubacterium saphenum, F. alocis, Megasphaera spp., and Peptostreptococcus spp. (65). The years that followed saw the development of “second-generation sequencing” techniques, commonly known as “next-generation sequencing” (NGS). These enabled massive parallelisation and improved automation and speed, and were also less expensive (64). It thus became possible to complete large-scale sequencing projects in just a few days or sometimes even hours (57). The two NGS techniques employed the most - 454 pyrosequencing and Illumina - are described in depth in Section I.4. Introduction 77 Advances in molecular techniques and the subsequent development of NGS tools led to an entirely new field of research: omics. Whole-genome shotgun metagenomics, commonly known as just metagenomics, allowed researchers to sequence the complete DNA (genome) of a single microbial culture or a complex microbial population, enabling the generation of reference genomes (57). This provided information not only on phylogenetic compositions but also on the genetic potential of a community to carry out distinct functional activities (58). Nevertheless, this technique does not produce data on the fraction of the metagenome being expressed. which can instead be obtained using techniques like metatranscriptomics, metaproteomics, and metabolomics that, respectively, assess the synthesis of transcripts, proteins, or the metabolic products of a specific set of activities (58). The human microbiome project (HMP) (2008-2012) later emerged as an initiative of the United States (US) National Institutes of Health (NIH) and used 16S rRNA gene sequencing and metagenomics’ techniques to initially identify and characterise the microorganisms associated with human health in different body parts, including the oral cavity (66). The research found that microbial profiles varied significantly, even between healthy subjects, meaning that taxonomic characterisation alone is not enough to uncover the relationship between the microbiota and the healthy state (67). In addition, metabolic pathway reconstructions of metagenomic data found that several pathways were ubiquitous in subjects and body habitats (67). Consequently, the second phase of the HMP (2013-2016), known as the Integrative HMP, comprised studies on dynamic changes in the microbiota and host as a result of physiological (pregnancy and preterm birth) or pathological (inflammatory bowel disease and pre-diabetes) conditions (68). This research has recently completed its first phase (69-71). Although each of these integrative investigations revealed new biology within their respective areas of health and disease, a surprising range of immune and ecological features of the host microbiota were commonplace (72). The combination of shotgun metagenomics, untargeted metabolomics, and immuno-profiling measurements have efficiently captured the host and microbial properties linked to disease. As in most studies of the microbiota, changes that occurred within individuals, populations, or phenotypes were often much smaller than the baseline variations between them. Accordingly, it is clear from this research that healthassociated microbiota interactions in individuals can manifest in extremely diverse ways (72). ALBA REGUEIRA IGLESIAS 78 Each of the integrative HMP studies found that other aspects of these interactions are highly localised and subject-specific. Microbial changes and associated host responses in the three conditions (individuals, populations, or phenotypes) were strongest when captured at the time the changes occurred and, often, within the tissue of origin. It is evident from these and other investigations that host-microbiota interactions have both localised and systemic effects. The NIH’s HMP has now come to an end but has revealed multiple new avenues of research and technologies for studies in the future (72). Figure 4. The first and second phases of the Human Microbiome Project. The image was taken from The Integrative HMP Research Network Consortium (72), an open-access article distributed under a Creative Commons Attribution 4.0 International (CC BY 4.0) license (https://creativecommons.org/licenses/by/4.0/). According to the traditional Socransky viewpoint, a set of Red Complex bacteria is thought to be the causative agent behind periodontitis (60). Now, however, NGS and omics’ technologies and techniques have detected the presence of bacteria like P. gingivalis in the absence of disease (73,74). It has also been found that the periodontal microbiota is more heterogeneous and diverse than previously thought, with new component organisms identified (75,76). These results have confirmed the hypothesis that periodontitis is initiated by the PSD of the entire microbial community (55). The PSD model states that different members or specific gene combinations perform distinct roles to shape and stabilise disease-provoking microbiota. Of particular relevance are the so-called “keystone pathogens”, which impair the Introduction 79 host’s immune response and elevate the virulence of the entire community through interactive communication with accessory pathogens. This dysbiotic microbial community is mainly composed of anaerobic genera from the phyla Bacteroidetes, Firmicutes, Proteobacteria, Spirochaetes, and Synergistetes (77). Consequently, while P. gingivalis is regarded as a keystone pathogen, S. gordonii and T. forsythia are viewed as accessories (77). Figure 5. Polymicrobial synergy and dysbiosis in periodontitis. The image was taken from Hajishengallis (77) with the permission of Springer Nature. I.3.2. Host response As noted previously, periodontitis is considered to be an inflammatory disease initiated by bacteria. However, despite the advances made in both periodontal microbiology and pathobiology, the issue of which comes first - the inflammatory response or the change to a dysbiotic subgingival microbiota - is still a matter of debate (78). ALBA REGUEIRA IGLESIAS 80 In physiological conditions, there is a balance between the local immune response and the microbiota. Bacteria are undoubtedly the principal cause of gingivitis, but it is the uncontrolled host-inflammatory and host-immune responses that largely drive tissue destruction, i.e., the progression of the disease (79). It is important in studies of the pathogenesis of periodontitis to consider the temporal sequence of microbiota changes on the way to periodontal inflammation (78). The microbial shifts induced by inflammation go beyond the overgrowth of certain species. Indeed, growth conditions also provide an environment that changes the physiology, pathogenicity, and expression of the virulence factors of the polymicrobial biofilm community (79). Consequently, the inflammatory response and the resident microbiota are linked in a bi-directional balance state in health and an imbalance state in disease (78). If the constant interplay between microbes and the host inflammatory response is viewed as a continuum, it is thus evident that specific bacteria cannot be regarded as initial causal agents in the pathogenesis of periodontitis (78). In a recently published review, Van Dyke et al. (78) presented a new model for the pathogenesis of periodontitis. Scientific evidence led them to conclude that chronic inflammation enables the development of a periodontal pocket that changes the redox and nutrient environment, thereby increasing the diversity and species richness of the biofilm. This results in dysbiosis, which reinforces and exacerbates inflammation as a way to initiate bone resorption. In light of both how inflammation mediates dysbiosis and the associated exacerbation of periodontal damage, the “inflammation-mediated polymicrobial-emergence and dysbiotic-exacerbation” (IMPEDE) model was, therefore, proposed by Van Dyke’s team (78). This is designed to complement the current classification of periodontal diseases (CPD) approach (4) and suggests that inflammation can be present for each classification stage as a principal driver of the clinical condition. The IMPEDE model recognises five stages (0–4) through which health, gingivitis, and periodontitis may develop, be contained, or progress. These stages are as follows: 0) periodontal health; I) gingivitis; II) initiation of or early periodontitis; III) inflammation-mediated dysbiosis and opportunistic infection; and IV) late-stage periodontitis. As set out in figure 6.B, inflammation-mediated polymicrobial dysbiosis and tissue damage can be exacerbated if no Introduction 87 Although other molecular markers are available, there are several reasons why the 16S rRNA gene has been regarded as definitive (86). First, it is present in all bacteria. Moreover, its structure and function have remained constant over time, suggesting that sequence alterations reflect random changes. These changes occur slowly enough that researchers can obtain data on all the prokaryotes. Moreover, the variability is such that both distant and close organisms can be distinguished. In addition, the relatively large size (1500 base pairs -bps-) of the gene makes it suitable for informatic purposes, and the conservation in its secondary structures favours accurate alignment. Finally, the ease with which the gene can be sequenced means that extensive, and constantly expanding, databases are available. Nonetheless, the employment of this gene as a phylogenetic marker has its limitations. One of the most important is the presence of variations in the 16S rRNA operon copy numbers per bacterial genome, with values ranging from 1 to 15 (87,88) or 1 to 17 (89). A higher 16S rRNA operon copy number per genome for a specific taxon will overestimate its relative and absolute abundance values (89). Although the number of copies appears to be taxon-specific, there are also variations among strains of the same species (87). Furthermore, the 16S sequences obtained from the same species or within the same genome are often different. It was generally accepted that two gene sequences differing by 1 to 1.3% or more represented two distinct species (90). Nevertheless, Pei et al. (91) observed that 24 of the 586 species they analysed had an intragenomic diversity higher than 1-1.3%. Another investigation evaluated the complete genomes of 2013 bacteria and archaea and revealed that there was intragenomic heterogeneity in 952 of them (88). Even though the majority of the divergence was below 1%, 119 genomes presented with higher values (88). As intragenomic heterogeneity is believed to overestimate microbial diversity, these authors recommended using primers targeting regions 4 and 5, which had the fewest variations (88). I.4.2. DNA Sequencing techniques DNA sequencing is the process that allows the nucleotide sequence of a DNA sample to be determined. Although the double helix structure of DNA was discovered in 1953, it was not until the 1970s that the first sequence of the human genome was obtained using first-generation sequencing techniques (64). The first two methods for DNA sequencing were reported in 1997: “chemical cleavage sequencing”, developed by Maxam and Gilbert (92), and “chain terminator ALBA REGUEIRA IGLESIAS 88 or Sanger sequencing”, which was discovered by Sanger and colleagues (93). Due to its simplicity and reliability, the latter approach has been the gold standard over the last three decades (94,95). The Sanger technique belongs to the sequencing-by-synthesis methods, which means that it uses the DNA synthesised by the DNA polymerase to identify the nitrogenous bases present in a DNA sequence. In addition to the four deoxyribonucleotide-triphosphates (dNTPs: dATP, dCTP, dGTP and dTTP), the sequencing reaction employs specific chain-terminator dideoxynucleotides (ddNTPs), i.e., nucleotides that lack a 3’-OH group (64). The incorporation of a ddNTP into a growing DNA molecule hampers the integration of a new nucleotide, as no phosphodiester bond can be formed because of the absence of the 3’-OH group. This means that the DNA synthesis is interrupted in that position (64,95). After several repetitions, this technique involves the products of the reactions being loaded in an agarose gel and subjected to electrophoresis. The ordered banding pattern thus obtained enables the sequence of the DNA template to be determined. Over time, the Sanger method has incorporated a series of innovations involving the automation of the process using fluorescent terminator dyes linked to the ddNTPs and the development of software to interpret and analyse the sequences (95). Consequently, the method is still valuable when high-throughput is not required. The leader in the field of automated Sanger-sequencing is Applied Biosystems, whose current commercial sequencers can generate 600-1000 bps of a proper sequence (95). The contribution of the Sanger approach to scientific advances in diverse areas has been invaluable. In 2001, the field of periodontal microbiology saw Paster et al. (96) use the method in the first comprehensive characterisation of the subgingival microbiota. The authors found that the subgingival niche harboured 347 species, 215 of which were novel phylotypes. They also estimated the number of unseen species in the population and determined that there were 68 additional taxa, accounting for a total number of 415 subgingival species. Introduction 89 As is already known, the years that followed saw the development of second-generation technologies, revolutionising the study of microbial diversity. Described below are the two NGS techniques used the most to examine the human microbiota. I.4.2.1. 454 pyrosequencing The first NGS application, known as 454, was introduced in 2005 by the biotechnology firm 454 Life Sciences, which was bought by Roche two years later (58,64). The technology comprises an initial emulsion PCR step, followed by subsequent pyrosequencing, which is a method of DNA-sequencing-by-synthesis based on the generation of light after nucleotides are incorporated in a growing DNA chain (64,94). The process of 454 pyrosequencing is as follows (94):  The DNA is isolated, fragmented, ligated to special adapters, and separated into single strands.  The DNA is amplified using emulsion PCR. The emulsion contains: the PCR reagents, the DNA template to be sequenced, the capture beads with the primers attached to them (complementing one of the adaptors), and another primer for the PCR. Emulsification takes place after controlled and vigorous agitation of the oil-water system. Millions of aqueous droplets are formed, with the amplification occurring inside them. Optimisation of the process guarantees that there is only one template and one bead in each droplet, meaning that millions of copies of the template are generated on each bead.  The DNA is denatured and the beads containing single strands are transferred to the wells of a picotiter plate. Only one bead is deposited in each of the several hundred thousand wells.  The DNA is sequenced through synthesis. This step requires: a single-stranded DNA sample, the sequencing primer, and the enzymes DNA-polymerase, adenosine triphosphate (ATP) sulfurylase, luciferase, and apyrase. Two substrates are also included in the reaction: adenosine 5’ phosphosulfate (APS) and luciferin. Cycles are then performed where each well receives, sequentially, one dNTP at a time. If there is complementarity, the DNA polymerase catalyses the incorporation into the DNA strand. The polymerase introduces a dNTP into a DNA nascent molecule, releasing pyrophosphate (PPi) in an amount equivalent to the quantity of the incorporated ALBA REGUEIRA IGLESIAS 90 nucleotide. The ATP sulfurylase converts the PPi into ATP in the presence of APS. Such an ATP is utilised by the luciferase to turn luciferin into oxyluciferin. This step produces light at an intensity proportional to the amount of ATP used. The light is detected by a camera and registered as a peak in a pyrogram, with the height of the peak being proportional to the number of incorporated nucleotides.  The system is regenerated by the enzyme apyrase, which degrades the ATP and the unincorporated dNTPs. The next nucleotide is then added. As the process advances, the complementary DNA strand grows and the nucleotide sequence is determined according to the pyrogram values. Introduction 91 Figure 9. The 454-pyrosequencing approach to amplifying single-stranded DNA copies from a fragment library on agarose beads. The image was taken from Mardis et al. (97) with the permission of Annual Reviews, Inc. Over the years, 454 Life Sciences has increased the length of its sequence reads and the number of bases per run. Its initial tools yielded sequence reads of 100 bps and up to 60 million bases per run. Later on, these sequence-read and bases-per-run figures reached 400 bps and approximately 500 million, respectively, on the company’s well-known Genome Sequencer (GS) FLX Titanium platform (94); subsequently, read lengths up to 700 bps and an output of 700 mega bps (Mbps) were available with the GS FLX + series (98). ALBA REGUEIRA IGLESIAS 92 The capacity to obtain a high number of reads in a single run is one of the main advantages of the 454 technology over the Sanger technique, as it enables the acquisition of much more sequence data (94). As more reads are generated in a single run, the cost per base is much lower than with the Sanger method. Furthermore, as the cloning step is not required with the 454 technology, the bias inherent in that procedure are avoided (94). Pyrosequencing’s use of “barcodes” (sequences introduced into the PCR primers) that work as unique sample identifiers enables sequences from different samples to be mixed in the same run. This increases efficiency and the number of outputs and also reduces costs (94). The high-throughput of the 454 technology affects both the depth (number of sequences per sample) and breadth (number of samples evaluated) of the sampling (94). A greater sampling depth increases the opportunities available for detecting low-abundance or rare community members, while more breadth allows additional samples to be analysed. This means that the results are more robust for comparison purposes (94). If the goal of a study is to determine the composition of a community at distinct sites, such as the skin vs. the oral-cheek mucosa, more samples should be studied than more sequences per sample (99). In contrast, if specimens from the same site or one close-by (i.e., the tooth surface vs. the gingival crevice) are being compared, deeper sequencing is required to identify minor differences in the community’s composition (100). Pyrosequencing also has limitations. One of the most important is related to the detection of long homopolymers, which can lead to sequencing bias and an artificial increase in the richness estimators (94). The reagent cost is another notable drawback (64). In 2013, Roche made the decision to close 454 Life Sciences because its technology was no longer competitive. Reagents are, however, still available from several suppliers (95). I.4.2.2. Illumina The first Solexa sequencer, named Genome Analyzer, was launched in 2006. Illumina acquired the company a year later and, since then, has brought a number of different sequencers to the market, with data output rates more than doubling annually (101): the initial Genome Analyzer could sequence 1 giga bps (Gbps) of data in a single run; by 2014, this figure had Introduction 93 increased to 1.8 tera bps (Tbps) with the HiSeqX Ten sequencer (101) or 6 Tbps with the NovaSeq6000 (98). Currently, Illumina is undoubtedly the NGS technique used the most (58,95). Like the 454 technology, it is based on the sequencing-by-synthesis principle (64,94), and resembles the Sanger method in that it also relies on the incorporation of dye terminator nucleotides into the sequence (94). However, the Illumina terminators are reversible, meaning that the polymerisation can continue even after the fluorophore detection. The workflow of this technology is as follows (64,101):  Library preparation: the DNA or cDNA sample is fragmented randomly and then ligated to two distinct types of adaptor at the 5’ and 3’ ends. Alternatively, the fragmentation and ligation reactions can be combined in a single step (“tagmentation”), which increases the efficiency of the process. The next stages, amplification and sequencing, take place in a solid surface called “Flow cell”.  Cluster generation employed to amplify DNA using bridge PCR. The fragments of the DNA template are denatured and the single strands are deposited in the flow cell, the surface of which is covered with primers that complement the adaptors. This complementarity enables the creation of bridges by joining the adapters to the primers (Figure 10.B). The PCR reagents are then incorporated (i.e., the nucleotides and the DNA polymerase) and each fragment is amplified into distinct clonal clusters via bridge amplification. When the second strand is formed, the DNA is denatured again so that new amplification cycles can be performed. It has been calculated that clonal clusters comprising about 1000 copies of each DNA fragment can be obtained; meanwhile, each flow cell can support millions of parallel cluster reactions (95).  Reversible termination sequencing: when the cluster generation is complete, templates are then ready for sequencing using reversible terminator nucleotides. Each dNTP is marked with a different fluorescent molecule, allowing all four of them to be added at the same time. The complementary nucleotides are then incorporated and those that are not are eliminated. After laser excitation, the emitted fluorescence is recorded by a four-channel fluorescent channel and the first base is identified. A surface chemical treatment removes the fluorescent dye from the incorporated nucleotides. This unlocks the 3’ carbon, enabling a new version to be used to continue the sequencing reaction. ALBA REGUEIRA IGLESIAS 94 The reactions are repeated for 300 or more rounds (95). Stacking and overlying the images obtained from all the cycles enables software to reconstruct the sequence of the DNA template fragment. Figure 10. The Illumina sequencing process. The image was taken from Mardis et al. (97) with the permission of Annual Reviews, Inc. Illumina produces a variety of sequencing tools for different applications. These include genomic sequencing, targeted sequencing, metagenomics, chromatin immunoprecipitation Introduction 95 (ChiP) sequencing, RNA sequencing, and methylation sequencing (95,101). As stated above, the distinct platforms have varying throughput levels: MiniSeq -7.5 Gbps of data with 25 million reads/run and read lengths of 2x150 bps; MiSeq - 15 Gbps with 25 million reads/run and read lengths of 2x300 bps; and NextSeq - 120 Gbps with 400 million reads/run and read lengths of 2x150 bps (95). Unlike the 454 technology, which sequences in one direction, the Illumina tools rely on “paired-end sequencing”, where both ends of the DNA fragments are sequenced and the forward and reverse reads are aligned as read pairs. So, using the MiSeq platform, a final length of 600 bps would be obtained after joining two 300-bp strands, which is represented by 2x300 bps. As well as producing double the number of reads in the same time and for the same effort concerning the library preparation, sequences arranged as read pairs enable more accurate alignment and the ability to detect insertion-deletion variants (101). Finally, paired-end sequencing also facilitates the detection of genomic rearrangements and repetitive sequence elements like gene fusions and novel transcripts (https://emea.illumina.com/science/technology/next-generation-sequencing/planexperiments/paired-end-vs-single-read.html?langsel=/es/). Figure 11. Paired-end sequencing and alignment. The image was taken from Illumina Inc. (101) with the permission of Illumina, Inc. A further advantage of the Illumina technology is that the four dNTPs are present during each sequencing cycle, meaning that natural competition minimises incorporation bias and reduces raw error rates. As a result, it is possible to achieve highly accurate base-by-base sequencing that almost eliminates sequences’ context-specific errors, even in repetitive regions and homopolymers (101). Moreover, due to its direct rather than camera-based imaging, fluorescent recording improves the detection speed (95). ALBA REGUEIRA IGLESIAS 96 Nevertheless, issues may arise with Illumina sequencing. First, extreme base compositions, i.e., guanine-cytosine (GC)-poor or GC-rich sequences, lead to an uneven coverage or even no coverage of the reads across the genome (102). Furthermore, the amount of template-DNA employed has to be quantified accurately to prevent the “overclustering” of the system. An analysis of the sequencing errors generated during the process can, however, be performed to identify the real sequence variants and the protocol-induced artefacts (95). A common disadvantage of NGS techniques is that the length of the reads produced might not be as long as required, hampering the identification of bacteria (94). Indeed, genomes often contain numerous repeated sequences that are longer than the reads obtained with the technology, which may lead to misassemblies and gaps (103). However, although it is necessary to sequence the entire 16S rRNA gene to reliably identify some species and describe new ones, it is accepted that just sequencing the initial 500-bp region is sufficient for distinguishing a high number of bacteria (104). What is more, changes in community composition can be assessed using gene fragments as small as 100 bps (105). As most of the NGS platforms that are currently employed generate short reads, bacterial identification using these methods has focused on the hypervariable regions of the 16S gene which are, as is already known, extremely informative (94). Nonetheless, NGS tools like GS FLX Titanium or GS FLX + can achieve read lengths of more than 400 bps and up to 700 bps, respectively. A final issue with short reads relates to the detection and characterisation of large structural variations (SVs), which can be challenging. Smaller variants, like single-nucleotide variations (SNVs) and short indels, can however be identified with greater accuracy (103). In conclusion, each of the NGS techniques has its particularities, and getting the most out of them requires researchers to strike a balance between target size, read length, depth, sequence accuracy, usability, and cost (105). I.4.3. Third-generation sequencing Although NGS methods have revolutionised biology, there was nevertheless a need to develop new techniques capable of overcoming the drawbacks described above (103). The shift from “long read” (e.g., automated Sanger) to “short read” (e.g., Illumina) technologies has led Introduction 103 An alternative fourth-generation technique is a spatial transcriptomics. This enables the visualisation and quantitative analysis of the transcriptome, with spatial resolutions in individual tissue sections (113). Specifically, this method combines in-situ transcript mapping with ex-situ transcript identification by way of NGS (112). These fourth-generation technologies, although still in their infancy, are now being used in projects aiming to unravel the functions of the brain as well as in cancer research (112). ALBA REGUEIRA IGLESIAS 104 I.5. ANALYSIS OF SEQUENCING RESULTS The NGS platforms generate a huge amount of genomic data that is unmanageable with an ordinary computer, meaning that informatics has become vital for processing and analysis purposes (114). In parallel with the use of high-throughput 16S rRNA gene sequencing, bioinformatics has emerged as a discipline that conceptualises biology in terms of macromolecules and then applies informatic techniques (applied maths, computer science, and statistics) to understand and organise the huge amount of data associated with these molecules (115). In other words, it is the tool used to make sense of sequencing results: signals are converted to data, data to interpretable information, and information into actionable knowledge (116). The concept of a pipeline refers to a set of bioinformatic algorithms executed in a predefined sequence to process NGS data. Accordingly, the data flow is transformed into a process comprised of several sequential phases where the input of each is the output of the previous stage. Different local or web-based software packages have been developed to manage the amplicon sequence data from NGS, including the commonly used quantitative insights into microbial ecology (QIIME) (117), mothur (118), and the ribosomal database project (RDP) pipeline (119). QIIME 2 has recently emerged (120); the first version is therefore no longer supported, as efforts are now focused entirely on its successor (http://qiime.org/). This new pipeline has the potential to serve not only as a marker-gene analysis tool, but also as a multidimensional and powerful data-science platform that can be quickly adapted to analyse diverse microbiome features. It also makes use of many new interactive visualisation tools that facilitate exploratory analyses and the reporting of results. Furthermore, QIIME 2 provides a software-development kit that can be used both to integrate the technology with other systems and develop interfaces targeted towards users with different levels of computational knowledge and experience (120). Regardless of the software or platform employed, the complete bioinformatics processing of 16S rRNA gene amplicon data usually encompasses three steps: 1) the pre-treatment or quality filtering of raw sequence data; 2) the construction of operational taxonomic units (OTUs) or single-nucleotide resolution table and 3) advanced data analysis and visualisation (98,121) (Figure 15). Introduction 105 Figure 15. Example of a flow chart of a 16S rRNA gene amplicon data analysis pipeline using OTUs. The image was taken from Ju and Zhang (121) with the permission of Springer Nature. I.5.1. Pre-treatment of raw sequence data The sequencing process generates a sff binary file. This provides general information on a run (number of flows, order of the nucleotides in the flows…) and follows it up with a description of each run sequence, indicating the position where it was generated, the pipeline used, and the total number of bases. It also provides information on the bases incorporated and the quality assigned to each of them. The pipeline transforms the sff file into readable fasta and qual versions by applying a particular command (121). This enables the first step to begin, i.e., processing the raw barcoded sequence data, which is a process known as demultiplexing. Once this phase has been completed, all the raw sequence data can be now demultiplexed, which is a procedure that sets it into individual subsets belonging to different samples based on specific barcodes (121). Raw reads are then quality filtered by applying criteria such as the minimum average quality score allowed in a read, the maximum number of ambiguous bases, the minimum and maximum sequence lengths, and the maximum length of homopolymer and maximum mismatches in the primer or barcodes (121). Next, the sequence-alignment step determines ALBA REGUEIRA IGLESIAS 106 where each short DNA sequence aligns with the reference genome, i.e., the 16S rRNA gene (122). This makes it possible to discern the location of the sequenced fragment in the gene. The PCR amplification and sequencing can introduce bias into the process, including PCR single-base errors, PCR chimeras, and sequencing errors, which must be checked and removed (121). If this is not done, the true diversity of the bacterial community would be overestimated. Particular attention should be paid to the chimeras or chimeric amplicons. These are artificial DNA sequences generated during PCR amplification and consist of a combination of two (or more) true underlying sequences (84). They appear when the extension step for an amplicon is brought to an end, with the short product obtained functioning as a primer in the next PCR cycle. This amplicon anneals to the incorrect DNA template and continues the extension, synthesising a single sequence sourced from two different templates. These chimeric amplicons can be overamplified in the steps that follow those described, thus creating unreal taxa and distorting the results. AmpliconNoise (123) and Denoiser (implemented in QIIME) are two of the most widely used software applications to remove or correct the PCR and any sequencing errors; meanwhile, ChimeraSlayer (124) (QIIME default method) and UCHIME (125) (mothur default method) are some of the tools employed to filter the PCR chimeras (121). I.5.2. Sequence clustering The second step in processing the 16S rRNA gene amplicon data begins by clustering the clean sequences into OTUs. An OTU is a cluster of organisms that are similar at the sequence level beyond a particular threshold, and which are intended to correspond to taxonomic clades (84,126). Sequence differences in the selected variability radius are assumed to be due to the variation within the taxonomic group or to random sequencer noise (127), which avoids the problem of differentiating biological from technical sequence variations but at the cost of taxonomic resolution (128). Several identity cut-offs have been used for the different taxonomic ranks. Typically, sequences are clustered at the ≥97% similarity threshold, which has been conventionally regarded as the species-level correspondent (98,129). Conversely, the MEGAN pipeline recommends thresholds of ≥99% and ≥97% for the species and genus levels, respectively (121). Introduction 107 Nonetheless, the sequence-similarity levels used are imprecise measures of an imprecise concept of a “species”, and the sequence identity of a given region of the 16S rRNA gene does not reflect the precise identity of the entire gene (105). The assignment of sequences to OTUs is known as “binning” (84), and numerous OTU clustering algorithms have been integrated into the popular sequence-analysis pipelines, such as QIIME2 (120), mothur (118), and USEARCH (130). Overall, they use three different strategies (121):  De novo: sequences are clustered without a reference database.  Closed reference: sequences are matched against a reference database; those unmatched at the given identity cut-off are discarded.  Open reference: sequences are first picked for closed-reference OTUs and the unmatched reads are subsequently clustered for de novo OTU versions. Currently, there is mixed evidence on which strategy is best when attempting to define OTUs and reveal the observations closest to the true community (131). Although the de novo approach enables the exploration of uncharted territories in the microbiota (105) and has been shown to create higher quality OTU classifications (132), the reference-based method has several advantages. First, sequence data from different gene regions, or generated from distinct sequencing technologies, can be combined using reference databases (105). In these cases, de novo OTU-picking might wrongly assign the same organisms to different OTUs based solely on variations in the DNA region amplified or in the sequencing technique (105). Second, the reference-based approach is increasingly valuable as the scope of publicly available data is expands, enabling new research to be interpreted in the context of existing studies (105). Picking OTUs against a reference database can also diminish the impact of chimeras and noise data (105). However, a single OTU can contain groups of sequences that could be individually assigned to different, related taxon (98) and the three OTU clustering approaches produce different results in terms of obtaining OTUs, even when using the same dataset (132,133). Moreover, the same method can yield distinct results after only a minor parameter change (134). ALBA REGUEIRA IGLESIAS 108 More recently, distinct error-correction or denoising approaches have become available, which are based on algorithms that use a single-nucleotide resolution (i.e., 100% sequence similarity) by generating amplicon sequence variants (ASVs), thus improving the taxonomic determination (128). These methods attempt to model the error of the sequencer and to cluster reads in a way that their distribution within clusters is consistent with such error (127). Among the most widely-known ASV-based pipelines there are DADA2, Deblur, and UNOISE (135-137); and they differ in how the above-mentioned correction is done (128). For example, DADA2 generates a parametric error model that is trained on the entire sequencing run and then applies that model to correct and collapse the sequence errors into ASVs (135). For its part, Deblur aligns sequences together into “sub-OTUs” and, based on an upper error rate bound along with a constant probability of indels and the mean error rate, removes predicted error-derived reads from neighboring sequences (136). Also, the UNOISE3 pipeline uses a one-pass clustering strategy that depends on two parameters with pre-set values that were curated by its author to generate “zero-radius OTUs” (137). Lastly, other algorithms use the 100% sequence similarity to create oligotypes, or minimum entropy decomposition nodes (138). Despite the different nomenclatures indicated by the respective researchers to refer to the clusters, they all are commonly known as ASVs. Different investigations have compared the two sequence clustering approaches to discern which performs better (127,128,139-142). In these studies, authors contrasted one (128,140- 142), two (127) or three (139) OTU to one (140,142), two (141) or three (127,128,139) ASV clustering methods, using sequences derived from mock (127,128,139,141), human gut (128,139,141), shrimp gut (142) and soil (128) samples. Also, other researchers used 16S rRNA gene sequences from the rRNA operon copy number database (rrnDB) for its analysis (140,143). In general, the ASV pipelines had demonstrated superior sensitivity, specificity, and precision, and lower spurious sequence rates when compared to OTU algorithms (127,139). Moreover, they allow for easier inter-study integration of biological features as the ASVs have intrinsic meaning independent of the reference database used, contrary to the study-specific nature of OTUs (139,144). Still, ASV-level pipelines are not free of limitations and can fail to distinguish very closely related true biological sequences and clump them together into a single ASV (139). Also, when analysing 16S rRNA gene data, Schloss (140) has recently affirmed Introduction 109 that the risk of splitting a single genome into separate clusters when using ASVs is of higher importance than the risk of grouping together ASVs from distinct taxa into the same OTU. Lastly, there is no consensus regarding the influence of the method chosen on the diversity results obtained. Meanwhile, some authors obtained minor differences between pipelines using the two clustering methods, with comparable alpha- and beta-diversity profiles (concepts that will be further explained) (141,142); others evidenced distinct results even among those from the same approach (128,139). In fact, Nearing et al. (128) found that, despite the similar general community structures, the alpha-diversity metrics varied considerably among all pipelines evaluated, even within the ASV-based DADA2 (135) and UNOISE (137). So, they concluded that the clustering pipeline choice will largely impact the alpha-diversity results among samples. In table 3 there is a brief description of several tools available for sequence clustering into OTUs or single-nucleotide resolution. It should be noted that, as said above, the main pipelines QIIME (117), QIIME2 (120), mothur (118), and USEARCH (130) had their own approaches for OTU clustering. Table 3. Tools available for sequence clustering. The table was modified from Zaura et al. (98), an openaccess article distributed under a Creative Commons Attribution 4.0 International (CC BY 4.0) license (https://creativecommons.org/licenses/by/4.0/). Tool Description UPARSE (126) Implemented in USEARCH (130). Algorithm for OTU clustering QIIME BLAST (145) An approach that matches the reads to the closest sequence in the database and groups the reads based on the BLAST label CD-HIT-OTU-MiSeq (146) Approach for clustering and annotation of MiSeq-based 16S sequence data UNOISE (137) Implemented in USEARCH (130). Creates high-resolution OTUs referred to as zOTUs Minimum Entropy Decomposition (MED) (138) Information theory-based clustering algorithm for sensitive partitioning of sequences Provides single-nucleotide resolution (oligotypes or MED nodes) DADA2 (135) Corrects Illumina-sequenced amplicon errors, providing singlenucleotide resolution as ASVs Deblur (136) Produces sOTU with single-nucleotide resolution (putative errorfree sequences) ALBA REGUEIRA IGLESIAS 110 Once the sequences are grouped, a single sequence is selected as a representative of each cluster. This sequence can be random, the longest, the most abundant, or the first in a cluster (121). The fact that each cluster is now represented by a single sequence also speeds up the posterior analysis. I.5.3. Taxonomy assignment After clustering, a taxonomic identity has to be assigned to each of the representative sequences. Phylogenetic relationships among sequences can be inferred either de novo or by using a reference database with an associated phylogeny (105). In the latter case, the taxonomic assignment process can be performed via two different strategies, the best hit and the lowest common ancestor (LCA) (121). Both of these require an alignment tool to compare the sequences being analysed against a reference database containing related sequences whose taxonomy is known (121). The difference between the best hit and LCA is that the former assigns a sequence based on the alignment with the highest score, while in the latter this is achieved via multiple hits against a particular database (121). The de novo approach can be carried out using tools like NAST (for sequence alignment) (147) and FastTree (for making phylogeny inferences from aligned sequences) (148). RDP (119), Greengenes (149), and SILVA (150) are among the most widely used databases at the taxonomic assignment stage and are used in combination with pairwise alignment tools like BLAST (145) or USEARCH/UCLUST (130). Some databases, such as CORE (76) and the human oral microbiome database (HOMD), are specialised in oral microbiota (151). They emerged to provide a comprehensive and minimally redundant representation of the bacteria that usually reside in the human oral cavity, with computationally robust classifications at the genus and species levels. In fact, although larger public databases like GenBank (152) and RDP (119) return named matches for a slightly higher fraction of sequences identified in analyses of clinical samples, CORE and HOMD are much more likely to do so accurately (76). Nonetheless, the larger databases are still important supplements to the specialised versions when it comes to recognising rare species. Performing a diversity analysis first requires the generation of a phylogenetic tree of OTUs or ASVs. Initially, the representative taxa sequences have to be aligned using tools like Introduction 111 MUSCLE (153) or PyNAST (154). The tree can then be constructed, which allows the relationships between the sequences to be visualised in terms of their evolutionary distance from a common ancestor (Figure 16). Many different packages have been developed for inferring phylogenies and building trees for multiple sequence alignments, with MEGA (155) being the most popular and versatile (121). At the very least, this type of software generates a table detailing the number of times that an OTU/ASV is observed and which taxa it represents (156). Figure 16. Phylogenetic tree based on 1,642 HOMD (151) reference sequences. The image was taken from Edlund et al. (157), an open-access article distributed under a Creative Commons Attribution 2.0 Generic (CC BY 2.0) license (https://creativecommons.org/licenses/by/2.0/). I.5.4. Advanced data analysis and visualisation Understanding the compositional differences of microbial communities is essential in the field of microbial ecology (158). In this regard, an OTU or ASV table enables different taxonomic summaries to be obtained that show the bacteria present and their relative ALBA REGUEIRA IGLESIAS 112 abundances at all taxonomic levels (156). However, further analysis is required to understand the quality of the data, the diversity within and between samples, and, ultimately, which statistical comparisons are needed to determine whether the microbiota has experienced flux or dysbiosis (156). A metadata table is also required to perform advanced exploratory and inferential analyses. The term “metadata” refers to the information associated with the sequences, including the environmental conditions and the time and location of the sample collection (105). Metadata is very important, as it is required if the aim is to replicate a particular investigation (159). As stated above, it is also essential for performing meaningful comparisons between samples or with specimens from other studies. Consequently, genomic sequence data that lack an environmental context have no value (159), meaning that aspects like the host’s health, sex, age, and diet, as well as the method of sampling, the size of the sample and its preparation, should all be recorded (159). I.5.4.1. Analysis tools There are several methods for comparing sequencing data, including QIIME (117), mothur (118), MEGAN (160), metagenomic-rapid annotation using subsystem technology (MGRAST) (161), UniFrac (162), DOTUR (163) and Metastats (164). Most of these tools are used for conducting analyses of bacterial communities and can detect groups of related samples (165). However, they do not provide information on the phenotypes or environmental conditions associated with these communities, and do not usually identify the biological features responsible for bacterial group relationships (165). In fact, only Metastats (164) explicitly couples a statistical analysis (to assess if the metagenomes differ) with the identification of biomarkers (to detect features characterising the differences), based on repeated t statistics and Fisher's tests on random permutations (165). Nevertheless, none of the aforementioned approaches produce explanations on biological classes with which to establish statistical significance and biological consistency, or estimate the size of the effects of predicted biomarkers. In 2011, Segata et al. (165) developed the linear discriminant analysis (LDA) effect size (LEfSe) method, which is a logarithm for detecting and explaining high-dimensional biomarkers. This couples standard tests for statistical significance Introduction 119 In the equation, Sest is the estimated taxa richness, Sobs is the observed taxa richness, f1 is the number of singletons, and f2 is the number of doubletons. This index is especially useful for datasets skewed towards low-abundance classes, as is likely to be the position in the case of microbes (179,182). Related to Chao1 is the abundance-based coverage estimator (ACE) (183), which not only considers the ratio of singletons and doubletons but also of all the taxa observed up to an arbitrary count, usually set at 10 (178): In the ACE formula (Sace), fabund is the number of taxa above the abundance threshold and frare is the number below it (rare samples) (178). It should be noted that the sum of these two values equates to the total number of taxa observed (179). Additionally, Cace is a samplecoverage estimator and yace is the estimated coefficient of variation for rare taxa (178). In their respective equations, nrare refers to the total number of individuals in rare taxa and fi to the number of taxa observed “i” times (178). Essentially, this index uses the number of rare taxa (≤10) and the number of singletons (f1) to estimate how many more undiscovered taxa there might be. Nevertheless, both Chao1 and ACE underestimate true richness in small sample sizes (179). Moreover, there are richness measures that consider the phylogenetic diversity (PD) of populations. The faith phylogenetic diversity index, first described in 1992, is a qualitative divergence-based measure that calculates the total branch length in a phylogenetic tree that ALBA REGUEIRA IGLESIAS 120 includes all the taxa in a sample (175,184). Although it fulfils the requirement of being a measure of taxon richness in a community, this index is highly sensitive to the sampling effort because it assumes that the total diversity of the population has been sampled (175). Moreover, the PD depends on the method employed to infer branch lengths on a tree, making it sensitive to errors during the tree’s construction (175). B. Evenness Along with richness, it is also important to measure the evenness of a sample’s distribution (178). Take, for example, two specimens, A and B, with the following compositions: Table 4. Example explaining microbial evenness. Bacterial phyla Sample A (number) Sample B (number) Actinobacteria 690 330 Fusobacteria 200 330 Proteobacteria 110 340 TOTAL 1000 1000 Both samples have the same richness (three types of bacteria) and the same total number of bacteria. However, B is more even than A, because the total number of bacteria is evenly distributed between the three phyla. Conversely, most of the bacteria present in sample A are Actinobacteria, with a few representatives of Fusobacteria and Proteobacteria. Consequently, specimen A is less diverse than B. As can, therefore, be seen, evenness is a measure of the relative abundance of the different taxa in a sample (174). In general, when richness and evenness increase, so does diversity (174). Diversity can be viewed as a summary of a community’s structure since membership, abundance, and evenness are taken into account (177). Traditionally, the Shannon-Weaver (185) and Simpson (186) indices have been used to estimate diversity (174) and, taken together, are a measure of species richness and evenness. However, neither of them is free of bias and, while the former attaches greater weight to species richness, the latter accounts more for species evenness (174). The Shannon-Weaver (alternatively: Shannon entropy or, simply, Shannon diversity) index (185) was originally proposed as a measure of entropy within the text and quantifies the Introduction 121 uncertainty of predicting correctly what the next individual taken from a sample will be (177) (156). So, if a sample contains 1000 bacteria and 900 of them are Fusobacteria, the probability that the next one is also Fusobacteria is high and the Shannon index value would be low (close to 0). In contrast, if every 100 bacteria belong to 10 different species, the ability to estimate what the next one will be is low and the Shannon index score would be high. The value thus increases along with the number of species and as the distribution of individuals among the species becomes more even (174). The formula of the index is as follows, with S being the number of taxa and pi the proportion of the community represented by a taxa i: This estimate enables the employment of a further measure: the Pielou evenness index (187), which divides the observed value of the Shannon index by the highest possible value, i.e., the value if all the species in a sample are equally abundant (178). The Simpson Index (186), first described in 1949, estimates species dominance and reflects the probability that two individuals taken at random from a sample will belong to the same taxa (177,178). Its values range from 0 to 1, with 0 being “infinite diversity” and 1 “no diversity”. Consequently, the score produced by the index increases as diversity decreases (174). This is described mathematically as follows: Here,  represents the Simpson Index, S the total number of taxa in the community, and pi the proportional abundance of each taxa i (178). Since the value of the index increases as diversity decreases, it is usually represented as its inverse, which is known as the inverse Simpson index (1/). Accordingly, an increase in diversity is mirrored by an increased inverse Simpson value (177), which represents the probability that two individuals randomly selected from a specimen will belong to different taxa. In contrast, Theta (Ɵ) is an example of an alpha-diversity measure that accounts for both evenness and the divergence between taxa (175). Simply put, it calculates the average ALBA REGUEIRA IGLESIAS 122 difference between two randomly chosen sequences or individuals in a population. Nonetheless, Theta has not been widely used to measure microbial diversity (175). Meanwhile, over recent years, Cadotte et al. (188) developed three indices of PD that also consider the relative abundance of each taxon in a community. Conversely, other authors have extended metrics like the Shannon and Simpson indices, transforming them into phylogenetically weighted equivalents. These have been shown to outperform the standard measures when it comes to distinguishing healthy from disease-associated human microbiota communities (178). I.5.4.2.2. Beta-diversity Beta-diversity is the measure of diversity between multiple samples (156), and describes how many taxa are shared between communities, including the absolute or relative overlap (84). Thus, a beta-diversity measure estimates the similarity between populations (84). Consequently, aspects of microbial ecology that are unapparent when examining the composition of individual specimens can be revealed by assessing the differences between samples (189). There are many different approaches for evaluating the similarity between communities and some of those used the most are described below. Conceptually, these capture different aspects of diversity. Traditional measures like the Jaccard (190) or Bray-Curtis indices (191) focus on the taxa compositional overlap, which is quantified directly from the taxa-count data (192). Considered to be the earliest beta-diversity index, the Jaccard accounts for the relative taxa overlap between two samples, i.e., the ratio of shared taxa among all the organisms sampled (192). This is an incidence-based or unweighted (qualitative) index: it only considers the presence/absence of taxa. Over the years, different abundance-based or weighted variations of the original version have been proposed, including the Chao weighted Jaccard index (193) and the weighted Jaccard index (192). In addition, the widely employed Bray-Curtis similarity index describes the community overlap as the fractional minimum abundance of shared taxa between samples (191). The calculation is performed using the following formula: Introduction 123 S1 and S2 are two samples and S1i and S2i are the abundances of phylotype i in samples S1 and S2 (178). Although it is not very sensitive, this index is nonetheless appropriate for use with cero-inflated datasets (178). Unlike the traditional measures, the recently developed phylogenetically-informed indices do not treat taxa independently. Instead, these metrics consider the phylogenetic relationships between taxa and quantify the shared evolutionary history between communities (192). Among these new measures is the widely known unique fraction metric (UniFrac), of which there are different versions. The unweighted form was the first to be released and only considers species presence/absence and counts the fraction of the branch length unique to either community (162). Conversely, the weighted UniFrac uses species-abundance data and weights the branch length with the abundance difference (194). In other words, it detects changes in the number of sequences from each lineage, as well as changes in the types of taxa that are present (175). The unweighted version is most efficient for detecting abundance changes in rare taxa, while its weighted counterpart is most sensitive for identifying differences in abundant organisms (158). Nevertheless, neither is particularly powerful when it comes to recognising changes in moderately abundant lineages (158). The third version released was the variance-adjusted weighted (VAW) UniFrac, which moderates the branch proportion difference by its variance, increasing the index’s power over the weighted version for detecting the differences between two communities (195). The VAW-UniFrac was used by Chen et al. (158) to introduce generalised UniFrac distances that unify the weighted and unweighted UniFrac versions within a common framework. This combined metric adjusts the weight on the branches to cover a series of distances, ranging from weighted to unweighted, and is designed to apply to, and identify, a much wider range of biologically relevant changes in a microbiota’s composition (158). More recently, Schmidt et al. (192) proposed a novel family of beta-diversity indices that quantify community similarity in the context of taxa-interaction networks: the taxa interactionadjusted (TINA) and the phylogenetic interaction-adjusted (PINA). The authors argued that because the indices take into account interactions between taxa, they are capable of quantifying new aspects of diversity and can expand possible biological interpretations of diversity patterns in new ways (192). ALBA REGUEIRA IGLESIAS 124 The distinct approaches to community dissimilarity described, i.e., count-based vs. phylogenetic, can highlight different aspects of a population and how it functions. Consequently, combining these different analyses to gain a deeper insight into the system under study may be a valuable next step (178). A. Multivariate analysis Multivariate analyses are supplanting simple descriptive investigations of bacteria, and are widely used in microbial ecology, where complex, multidimensional datasets abound. However, the employment of OTUs or ASVs abundances makes it difficult to test the direct association between the composition of the microbiota and environmental factors, due to the high dimensionality, non-normality, and phylogenetic structure of the data (196). Consequently, multivariate analyses first require the researcher to select a methodology for measuring distance before conducting an analysis of estimated distances (196). A distance measure defined between any of two samples can be utilised. Among the numerous metrics that exist, the above-mentioned Bray-Curtis and UniFrac are two of those employed the most. Many types of multivariate statistical analyses have been used for the assessment of highthroughput datasets, and novel approaches for analysing large-scale datasets are also being developed (197). These methodologies can be categorised based on criteria such as the technique’s goal (e.g., interpret relationships, test statistical significance), the type of mathematical problem (regression, ordination, calibration, classification), or the variable response (e.g., linear, unimodal, mixture distribution) (197). These techniques can also be classified according to the primary research objectives, and three categories can be distinguished (197):  Exploratory methods: these are used to explore the relationships among objects (e.g., samples or sites) based on the values of the variables measured in those objects. These techniques provide a valuable visualisation of object similarities since similar objects are usually positioned close together on the visualisation plot, while dissimilar objects are wide apart.  Interpretive methods: these ‘constrained’ techniques use both the main set of measured variables and another of additional explanatory variables. Introduction 125  Discriminatory methods: these are an extension of the former techniques and are usually known as discriminant analyses (DAs). The goal of DAs is to define discriminant functions (synthetic variables) or hyperspace planes that maximise the separation of objects among different classes (groups). Of the exploratory approaches, the principal component analysis (PCA) is one of the most widely used and oldest (198). In the main, it is employed to calculate new synthetic variables (principal components), which are linear combinations of the original variables, and accounts for as much of the variance in the original data as possible (199). The first principal component (PC) represents the axis in the multidimensional data-space that would produce the largest dispersion of values. Other PCs are calculated as being orthogonal to their predecessors and are positioned along the largest remaining scatter plot of the values. Consequently, the PCA creates a rotation of the original system of coordinates, meaning that the PCs are orthogonal to one another and correspond to the directions of the greatest variance in the dataset (197). In this ordination plot, the first PC axis represents the largest variability gradient, PC2 is the second largest, and so on, until all the dataset’s variability has been assessed (197). Each object can be given a new set of coordinates in the PC space, and its distribution in such a space will correspond to the similarity of the variables’ scores for those objects (197). A conceptual extension of the PCA is the principal coordinate analysis (PCoA) (200). Similarly, this is used to order objects along PC axes in an attempt to explain the variance in a dataset (197). While a PCA organises objects by analysing a correlation or covariance matrix, the PCoA can be applied to any type of distance metric (197). The use of this method has recently increased in ecology since it can employ measures of phylogenetic distance and community composition to calculate the similarity among populations (197). As a distance matrix is an input file, the PCoA cannot directly relate any of the measured variables to individual coordinate axes (199). Instead, an indirect correlation or regression analysis of object values vs. object scores can be used to estimate the contribution of a variable to object dispersion along a particular PC axis (197). ALBA REGUEIRA IGLESIAS 126 Figure 22. Graphical representation of (A) a principal component analysis (PCA) and (B) a principal coordinate analysis (PCoA) plot. This image is adapted from Rao et al. (201), which is an open-access article distributed under a Creative Commons Attribution 4.0 International (CC BY 4.0) licence (https://creativecommons.org/licenses/by/4.0/). Finally, non-metric multidimensional scaling (NMDS) is another exploratory method in which a number of ordination axes are explicitly chosen in advance, after which data are fitted to those dimensions (202). As in the PCoA, a matrix of object dissimilarities is first calculated using a distance metric. The ranks of these distances for all the objects are calculated and then the algorithm identifies a configuration of objects in the N-dimensional ordination space that best matches the differences in ranks (197). In an NMDS ordination, the proximity between objects corresponds to their similarity, but the ordination distances do not correspond to the original distances between the objects (199). Introduction 127 Figure 23. Graphical representation of a non-metric multidimensional scaling (NMDS) plot. The image was taken from Drell et al. (203), an open-access article distributed under a Creative Commons Attribution 4.0 International (CC BY 4.0) license (https://creativecommons.org/licenses/by/4.0/). Interpretative methods for analysing large-scale datasets can be further subdivided into three types: symmetric, asymmetric, and statistical-significance testing. The first compares two datasets and does not distinguish between explanatory and response variables (197). Examples are the canonical correlation analysis (204) and the procrustes analysis (205). In contrast, the asymmetric approaches use two distinct sets of variables: one explanatory or independent and one response or dependent (197). The redundancy analysis (RDA) (206) and generalised linear models (207) are examples of asymmetric techniques. Specifically, the commonly used RDA is a type of constrained ordination that evaluates how much of the variation in one set of variables (response) can be explained by the variation in another set (explanatory) (206). This is a canonical version of the PCA and is based on similar principles, with the PCs constrained as linear combinations of the explanatory variables (197). The RDA provides a useful indication of, for example, how much the variation in species distribution is due to differences in the environmental factors between sites (197). ALBA REGUEIRA IGLESIAS 128 Figure 24. Graphical representation of a redundancy analysis (RDA) plot. The image was taken from Qi et al.(208) with the permission of Elsevier. The third interpretative method involves the statistical-significance testing of multivariate datasets (197). Several approaches are available for analysing among-group differences in microbiota data, such as the permutational multivariate analysis of variance (PERMANOVA) (209), the analysis of similarities (ANOSIM) (210), the multi-response permutation procedure (MRPP) (211), and the Mantel test (212). Of these, the PERMANOVA and ANOSIM are the most widely used in microbiota studies and are generally employed with a distance measure (196). These tests make it possible to evaluate elements like microbial divergence or similarity in populations or the factors affecting such communities. The significance of the results can also be confirmed through visualisation methods. Of the final examples of discriminatory methods, the discriminant function analysis (DFA) (213) and the random forest (214) should be highlighted. The DFA, better known as the LDA, is a method for evaluating how well a group of variables supports an a priori grouping of objects. Here, the measured variables are the predictor variables, while the variable defining the object classes is treated as the response variable (also called the grouping variable) (197). The LDA is closely related to other lineal methods like the PCA. However, unlike the PCA, it derives synthetic variables that specifically maximise the between-class group dispersion (197). As each discriminant function is a weighted linear combination of the measured predictor variables, the weights (called discriminant coefficients) can be used to define the contribution of each predictor variable to the observed discrimination between classes of objects (197). The Introduction 135 the network (229). Centralisation is high when a vertex has high centrality values and those of the other vertices are low. Conversely, if the centrality is distributed more evenly, the network centralisation is low. Additionally, groups of vertices may form a module or cluster, thus acting as a sub-network within the main network. A module is defined as a set of strongly related nodes that are, in turn, less related to those that do not belong to the group (223). A network is said to have high modularity if it presents dense connections within node clusters and sparse connections between different groups of vertices (223). All of the network measures described above can be calculated for the modules or clusters. Finally, as referred to earlier, the capacity to identify hubs or keystone taxa, which are highly connected OTUs or ASVs in the microbiota, is one of the most useful features of a cooccurrence network analysis (230). Different measures have been adopted to define these hubs in microbial communities. Banerjee et al. (222), for example, aimed to provide a quantifiable threshold for the consistent identification and validation of keystone taxa. Their findings led to a recommendation that the high mean degree, high CC, and low BC scores should be combined to this end. Nevertheless, it should be noted that the identification of highly connected OTUs or ASVs in a microbial network does not necessarily reveal their role as keystone taxa (231), which are very closely linked species that exert a considerable influence on the structure and functioning of the microbiota, irrespective of its abundance (222). Although further experimental evidence is required before network hubs can be defined as keystone taxa (231,232), identifying them is nonetheless a valuable step, since this will help researchers to target key community members (231). In conclusion, the findings from the co-occurrence network analyses described in the literature should be viewed with caution: they may be affected by methodological differences concerning, for example, the correlation values employed as cut-off points (233) or the use of different definitions of keystone taxa (222). ALBA REGUEIRA IGLESIAS 136 I.5.4.2.5. Predictive models Machine learning (ML) is a computer science discipline in which computers are programmed to learn patterns from the data in a multi-dimensional dataset and produce classifications or predictions based on statistical associations (234). The field has two main approaches, supervised and unsupervised, and the goals of the research determine which is the most appropriate (Table 5). The latter is employed to identify the underlying structures or relationships between variables (samples) in a dataset and is well suited to the visualisation of high-dimensional input data (234,235). Indeed, the PCA (198) and PCoA (200) mentioned earlier are examples of unsupervised ML algorithms. In contrast, supervised learning involves the classification of an observation into one or more categories or outcomes (235). As a consequence, this requires training data, with each training sample having values for a number of independent variables or features, as well as an associated classification label (236). Predictive modelling is a set of mathematical processes and computational techniques that enable the probability that an event will occur to be inferred from a set of previously obtained data. There is a close relationship between predictive analyses and ML, since predictive models typically include an ML algorithm. Unsupervised models do not require labelled data or the use of error-measurement metrics, with the algorithm instead searching for patterns among the input or the output variables (samples) during the modelling process. Conversely, a supervised model needs such data to generate a predictive model, as well as an error-measurement metric to improve it during the building process. Introduction 137 Table 5. Summary of the supervised and unsupervised machine-learning approaches. The table was taken from Reel et al. (237) with the permission of Elsevier. Learning approach Goal Description Unsupervised Identify clusters Unsupervised learning employs input variables without a target/output variable to find the underlying patterns in unlabelled data. It can be used for clustering, anomaly detection and dimensionality reduction. Supervised Predict new data Supervised learning involves fitting a model with labelled training data and then using it for predictive purposes. The problem it addresses can be classed in terms of either regression (the predicted variable is numeric) or classification (the predicted variable is categorical). The three steps of supervised learning are: 1) fitting a model from the sample’s input observations; 2) evaluating the model and then extensively tuning its hyper-parameters; and 3) setting up the model for the production stage and using it to make predictions. When a model is accurate in relation to both the training and test data, it is said to have learned properly. However, a particular ML output might predict the training data with a high degree of accuracy, but nonetheless fails to produce precise predictions with the test data (overfitting); it may even be unable to predict the training data correctly (underfitting) (234). In comparison to the very low number of samples and clinical conditions typically evaluated in this kind of work, the oral microbiota high-throughput data in this study is characterised by a large quantity of independent and predictor variables (OTUs or ASVs). This often indicates a high degree of multicollinearity and, as a result, produces very poorly conditioned problems (238). In a supervised framework, one solution is to reduce the dimensionality of the data, either by feature selection or introducing artificial variables that summarise most of the relevant information (238). To this end, several tools for supervised predictive modelling have been proposed, including the aforementioned random forest plots (214), support vector machine (SVM) (239), and regression models like the sparse partial leastsquares discriminant analysis (sPLS-DA) (238). The SVM device is a method for identifying a decision boundary to enable the classification of data. An SVM training algorithm is applied to a training dataset with information about the class to which each piece of data belongs, establishing a hyperplane that separates two classes. Next, the SVM seeks to optimise the width of the gaps between classes, i.e., the maximum-margin hyperplane. The resulting model can be used to determine whether a new data element is, or is not, a member of a particular class (240). Among the advantages of ALBA REGUEIRA IGLESIAS 138 this method are its efficiency at learning complex classification functions and its employment of powerful regularisation principles to prevent overfitting (241). However, in the case of highly dimensional datasets, the results obtained are often difficult to interpret given the large number of variables (238). Furthermore, multiclass classification problems require either their decomposition into several binary problems or the definition of multiclass objective functions (238). In contrast, the sPLS-DA (238), which is a natural extension of the PLS-DA, is based on the assumption that only a small number of features are responsible for driving a biological event or effect, enabling predictor variables to be selected and classified in a one-step procedure (238,242). The proper functioning of this method has been demonstrated previously (243), and multiple classes (e.g., clinical conditions) can be distinguished at the same time. Although there is some difficulty in construing the results obtained with this model compared to those that have only two classes, a graphical representation makes the interpretation process easier (238). Figure 27. Graphical representation of an sparse partial least-squares discriminant analysis (sPLS-DA). This image is adapted from Rohart et al. (244), which is an open-access article distributed under a Creative Commons Attribution 4.0 International (CC BY 4.0) licence (https://creativecommons.org/licenses/by/4.0/). The implementation of an sPLS-DA (238) in the mixOmics package (244) of R- Bioconductor (245) enables the following to be determined for each model (246): Introduction 139  The number of components or latent variables. There are as many dimensions of the sPLS-DA model as required.  A set of loading vectors, which are coefficients assigned to each independent variable (OTU or ASV) to define each component. In other words, they indicate the importance of each variable in the sPLS-DA and each loading vector is associated with a particular component. These vectors are obtained in such a way as to maximise the covariance between a linear combination of independent variables and the classes of interest.  A list of designated variables associated with each component. The procedure implemented in the mixOmics package (244) uses a k-fold cross-validation technique to automatically calculate the appropriate number of dimensions for each model. The rule of thumb is that the number of dimensions will be K-1, where K indicates the number of classes (e.g., clinical conditions) included in each model. Accordingly, a variable (e.g., OTU or ASV) that is present in multiple models will be identified as important. It becomes critical or sees its value increase if its presence is always associated with the same class.  Once the optimal parameters have been chosen (number of components and variables), the final model is run on the entire dataset. The model’s classification error rate is then estimated (244), and an additional accuracy evaluation using the receiver operating characteristic (ROC) and AUC can be performed (244). The use of predictive modelling to identify oral taxa that can distinguish between health conditions and are associated with specific disease states would be extremely valuable for determining the biomarkers of disease (236). However, to date, very few sequencing-based studies of the oral cavity have conducted predictivity analyses (247-250). ALBA REGUEIRA IGLESIAS 140 I.6. AN OWN STUDY ON THE RELATIONSHIP BETWEEN DENTAL AND PERIODONTAL HEALTH STATUS AND THE SALIVARY MICROBIOTA In a recently published 16S rRNA gene sequencing-based investigation of our investigation group (249), several of the above explained analysis tools were applied to examine the bacterial diversity and the co-occurrence network patterns of the salivary microbiota in patients who have been clinically classified by a self-designed and previously validated scale of overall oral health. Moreover, we evaluated the diagnostic potential of the salivary microbiota to discriminate between different clinical conditions. The own scale of overall oral health consists of three dental- and three periodontalassociated parameters (Table 6). Taking this into account, the participants’ dental and periodontal grades (DG and PG, respectively) corresponded to the grades assigned to at least two of the three variables analysed in each of these two categories. If there were differences between the grades allocated to each of the variables in a category, the parameters for “number of caries” and “number of periodontal pockets ≥4 mm” took precedence. If the same grade was allocated to two variables in a category, but the third variable’s grade was two levels higher, the value assigned to it was one grade higher than that of the matching variables. Lastly, the oral health grade (oral grade, OG) was determined by the category (dental or periodontal) with the highest ranking, enabling patients to be classified based on the score for their dental and periodontal health and the combination of both conditions. Table 6. The scale of overall oral health, involving grades of dental and periodontal health. The table was taken from Relvas et al. (249), an open-access article distributed under a Creative Commons Attribution 4.0 International (CC BY 4.0) license (https://creativecommons.org/licenses/by/4.0/). Caries severity, 1= affecting the enamel, 2= affecting the enamel and dentine, and 3= affecting the enamel, dentine, and pulp. Grade 0 Grade 1 Grade 2 Grade 3 Grades of dental health Supragingival plaque (O’Leary index) (251) 0 1–56 57–112 >112 Caries 0 1–4 5–8 ≥9 Severity of the caries (median) 0 1 2 3 Grades of periodontal health Gingival inflammation (Ainamo and Bay index) (252) 0 1–56 57–112 >112 Periodontal pockets ≥4 mm 0 1–56 57–112 >112 Severity of the pockets (mean) <4mm 4-4.9mm 5–5.9mm ≥6mm Introduction 141 Unstimulated saliva samples were collected from each participant. Sequencing of the 3-4 region was performed in an Illumina MiSeq platform with 2 × 300 bps reads, while the raw reads were processed according to the mothur pipeline (118). The statistical analysis of the 16S rRNA sequencing data at the species level was conducted using the phyloseq (166), DESeq2 (170), Microbiome (173), SpiecEasi (226), igraph (227), and mixOmics (238) packages. I.6.1. Results and discussion The overall oral health scale was used to produce a convenience sample of 81 patients that were given the following OGs: 0 for 17 of them; 1 for 25; 2 for 28; and 3 for 11. In relation to the subscales, 47 patients had a PG of 0 and different DGs (17 had a DG of 0, nine a DG of 1, 11 a DG of 2, and 10 a DG of 3), and 46 had a DG of 0 and different PGs (17 had a PG of 0, 14 a PG of 1; 14 a PG of 2; and one a PG of 3). The four patients excluded due to a low number of raw sequences obtained were: two of OG1 and two of OG2 (two of DG0, one of DG1, and one of DG2; two of PG0, one of PG1, and one of PG2). I.6.1.1. Impact of the dental and periodontal subscales and the scale of overall oral health on the salivary microbiota: alpha diversity indicators and the structure of the bacterial community Worsening dental or periodontal health revealed a trend of increasing alpha diversity in both the dental and periodontal subscales; although significant differences in the number of OTUs were only observed in the former: DG0 vs. DG123 (p= 0.009), and DG0 vs. DG23 (p= 0.006). The DG23 also showed a trend towards increased diversity (Shannon Index) and evenness (Pielou Index). On the contrary, other authors had observed that the saliva samples from the healthy and caries groups generally had similar levels of richness and diversity (253-256). This was true whether the diseased group was composed of subjects with active (254), inactive (255), or cavitated (256) caries. However, it should be noted that the dental health subscale used here incorporates variables that not only include the number of caries and their severity, but also the levels of supragingival plaque, all of which could affect the bacterial richness of the salivary community. In addition, there is some inconsistency between the alpha-diversity results of various studies in the literature and our findings when comparing periodontally healthy and periodontitis subjects. Several researchers described greater richness (248,257), diversity (248), ALBA REGUEIRA IGLESIAS 142 and evenness (257) in the saliva samples of patients with periodontitis than in those who were healthy; meanwhile, others, as observed here, only observed a trend of increased alpha diversity with worsening of periodontal health (258). In the overall oral health scale comparisons, grade increments were linked to progressive increases in bacterial richness and the Shannon Index values, especially in OG0 vs. OG23 (p= 0.013 and p= 0.026, respectively). The potential impact of the simultaneous presence of dental and periodontal disease on the richness and diversity of the salivary microbiota is in line with the results of Takeshita et al. (259). Like other studies in which the structure of the global salivary microbiota is similar in patients with good oral health (260-262), our PCoA revealed a grouping of the salivary samples taken from the participants with DGs, PGs, and OGs of 0. This contrasted with the picture for the other grades, whose compositional distributions were more diverse (Figure 28). The visual observation was confirmed by the PERMANOVA test, which produced significant results for the comparison of grades 0 and 123 (dental subscale, p= 0.0009; periodontal subscale, p= 0.0229; oral scale, p= 0.0008). These findings were mainly at the expense of the contrast between grades 0 and 23 (dental subscale, p= 0.0005; periodontal subscale, p= 0.0287; oral scale, p= 0.0008). Focusing on the group with the highest grade of oral pathology (OG23), PERMANOVA's test revealed that the structure of the salivary microbiota was different depending on the predominance of dental pathology (PG0_DG23) or periodontal pathology (DG0_PG23) (p= 0.027). Introduction 143 Figure 28. Principal Coordinate Analysis including PERMANOVA test values in the comparison between different grades of dental, periodontal, and oral health. The image was taken from Relvas et al. (249), an open-access article distributed under a Creative Commons Attribution 4.0 International (CC BY 4.0) license (https://creativecommons.org/licenses/by/4.0/). I.6.1.2. Impact of the dental and periodontal subscales and the scale of overall oral health on the salivary microbiota: composition of the core microbiota and testing differential abundance There are numerous 16S rRNA-based microbiome studies on salivary microbiota in the literature that have only analysed the differential abundance of the taxa associated with various oral conditions (256,257,263,264). In our study, we considered that it was essential to evaluate salivary microbiota from a dual perspective: the prevalence of the taxa determining the core microbiota; and their differential abundance in relation to the different DGs, PGs, and OGs. ALBA REGUEIRA IGLESIAS 144 I.6.1.2.1. Composition of the core microbiota The core microbiota associated with the participants’ dental and periodontal health contained 57 species, representing 14.14% of the total number of OTUs and 63.06% of the total abundance. There were only nine taxa in DG0 and eight in PG0 (the specific core of grade 0), exemplifying abundances of 7.80% and 1.34%, respectively. Of these specific core species, five were common to both the dental and periodontal health conditions: Neisseria macacae, Butyrivibrio sp. HMT 455, Campylobacter concisus, Porphyromonas catoniae, and Corynebacterium durum. There were 66 species in the core microbiota associated with the most severe dental disease (DG23) and 73 with the most serious periodontal disease (PG23), representing 16.37% and 18.11% of the total salivary microbiota, and 67.14% and 67.98% of the total abundance, respectively. There were only eight and 10 taxa present in DG23 and PG23 (the specific core of grade 23), exemplifying 2.54% and 2.46% of the abundance, respectively. Of these specific core species, only P. endodontalis was common to both pathological conditions. There were 35 taxa common to both the dental and periodontal subscales, regardless of the grade (non-specific core), representing abundances of 50.83% and 52.75%, respectively. Of these non-specific core species, 25 were common to both subscales, with the most abundant (abundance > 1%) being: Granulicatella adiacens, Haemophilus parainfluenzae, Leptotrichia sp., Porphyromonas pasteri, Prevotella sp., Prevotella melaninogenica, Prevotella salivae, Rothia mucilaginosa, Streptococcus sp., S. oralis subsp. dentisani clade 058, Streptococcus salivarius, Veillonella sp., and V. parvula. The core species in the literature on the salivary microbiome have various definitions and, as a consequence, any associated findings are difficult to compare (247,254,258,259). Despite this, our investigation demonstrated for the first time that the non-specific core of the salivary microbiota comprises a greater number of species in higher abundances than the specific-core associated with a particular dental or periodontal condition. Interestingly, more than half of the non-specific core species in the present series were the same as those previously identified by Takeshita in ≥ 75% of Japanese adults (259) as G. adiacens, H. parainfluenzae, or R. mucilaginosa, among others. These results confirm that several bacterial taxa in the salivary Objective 2 247 Table 1. Selected primer pairs with high in silico coverage percentages targeting oral bacteria and/or archaea and those most used in the sequencing-based studies of the oral microbiome. Primer pairs were selected based on the species coverage values (number of species detected /total species evaluated) in objective 1 (21). They were individually evaluated through regular expressions against Escherichia coli J01859 to define their positions. The U values represent a mismatch on the assessment and, therefore, the position cannot be confirmed with a guarantee. Gene regions were delimited as described by Baker et al. (30). *The primer pair gene region was 3 and **1-10 when considering the mode positions obtained in the analysis against the oral-bacteria and oral-archaea databases in objective 1 (21). ALC= amplicon length category; bps= base pairs; F= forward; KP= Klindworth primer; L= long mean amplicon length category, >600 base pairs; M= medium mean amplicon length category, 301-600 base pairs; OP= oral primer; Post= position; R= reverse; S= short mean amplicon length category, 100-300 base pairs; U= unaligned with Escherichia coli; V= variable. ALBA REGUEIRA IGLESIAS 248 The direct and reverse sequences of each primer pair selected were used in combination with Python’s regex module (31) to obtain, in silico, the amplicons of the 16S rRNA genes identified in all of the chosen genomes. For each primer pair, we determined: the mean size and number of the 16S rRNA gene amplicons; the number of gene variants; the number of genomes and species detected; and the percentage of coverage at the species level with no MAs (SCNMA). This coverage value was calculated as: SC-NMA (%)= [(Number of species detected - Number of species with MAs)/Total number of species evaluated] x 100 The overestimation of abundance at the species level (the overestimation factor -OF-) was also calculated. This represented for each species the combination of the number of copies of the 16S rRNA gene amplicons and the number of MAs. To remove the overestimation derived from the intragenomic gene redundancy, the OF of each species was divided by the number of gene copies, resulting in OF caused by the presence of MAs (OF-MA). Species with values equal to 1.00 did not have amplicons that matched other species for the corresponding primer pair, while those with estimates greater than 1.00 did. For each primer pair, both parameters were expressed cumulatively and as an average. The best primer pairs selected first were those with the highest SC-NMA value and of these, those with the lowest OF-MA value. The worst primer pairs were those with the lowest SC-NMA and the highest OF-MA. Objective 2 249 2.4. RESULTS 2.4.1. Number of intragenomic 16S rRNA genes in oral-bacteria and oral-archaea genomes Table 2 details the mean number of intragenomic 16S rRNA genes in the bacterial and archaeal phyla through seven taxonomic ranks. The 518 oral-bacteria genomes examined had a mean size of 2,933,660.68 bps and an average number of 4.55 intragenomic 16S rRNA genes, which in turn had a mean size of 1,501.32 bps and an average of 2.60 variants. Eleven of the 186 bacterial species (5.91%) had one gene/genome, 159 species (85.49%) showed a mean between two and six genes, and 16 species (8.60%), mean values of seven or more genes. The maximum mean number of intragenomic 16S rRNA genes observed was 10.83 in Bacillus anthracis, with five strains of this species having a total of 11 genes/genome. Concerning the average number of intragenomic gene variants, 63 bacterial species (33.87%) presented one variant/genome, 118 species (63.44%), between two and six, and five species (2.69%), seven or more. The 191 oral-archaea genomes had a mean size of 2,545,441.40 bps and an average of 1.95 intragenomic 16S rRNA genes, which in turn had a mean size of 1,471.25 bps and an average of 1.44 variants. Sixty-four out of the 135 archaeal species (47.41%) had a mean of one gene/genome, 67 species (49.63%) showed an average between two and three genes, and 4 species (2.96%) mean values above three (Methanobacterium formicicum, Methanococcus vannielii, Methanosphaera stadtmanae, and Methanospirillum hungatei). At the strain level, the maximum total number of genes/genome increased to five in Methanococcus maripaludis (unknown strain) and Sulfolobus acidocaldarius (unknown strain). Concerning the average number of intragenomic gene variants, 93 species (68.89%) had an average number of one variant/genome and 42 (31.11%) had between two and three. Appendices S3 and S4 contain the sizes of the bacterial and archaeal genomes and genes, the number of genes/genome, and the number of gene variants/genome across eight taxonomic ranks. ALBA REGUEIRA IGLESIAS 250 Table 2. Intragenomic 16S rRNA genes in the bacterial and archaeal phyla through seven taxonomy ranks. Mean number of intragenomic 16S rRNA genes Taxonomy level Phylum Phylum Class Order Family Genera Species Strain No. genomes Actinobacteria 3.12 3.19 – 2.00 3.41 – 1.33 4.55 – 1.10 4.55 – 1.00 5.00 – 1.00 5 - 1 91 Bacteroidetes 3.68 4.75 – 3.44 4.75 – 3.44 4.75 – 2.00 4.75 – 2.00 7.00 – 2.00 7 - 2 24 C. Saccharibacteria 1.00 1.00 1.00 1.00 1.00 1.00 1 1 Chlamydiae 1.00 1.00 1.00 1.00 1.00 1.00 1 - 1 5 Chlorobi 2.00 2.00 2.00 2.00 2.00 2.00 2 1 Chloroflexi 2.00 2.00 – 2.00 2.00 – 2.00 2.00 – 2.00 2.00 – 2.00 2.00 – 2.00 2 - 2 2 Firmicutes 5.43 5.52 – 3.25 6.61 – 3.25 9.85 – 2.00 10.18 – 2.00 10.83 – 2.00 11 - 2 177 Fusobacteria 4.35 4.35 4.35 4.40 – 4.33 4.75 – 3.00 5.00 – 3.00 5 - 2 21 Ignavibacteriae 1.00 1.00 1.00 1.00 – 1.00 1.00 – 1.00 1.00 – 1.00 1 - 1 2 Proteobacteria 5.21 6.13 – 2.17 6.98 – 2.17 7.14 – 2.00 8.00 – 2.00 8.00 – 2.00 8 - 2 170 Spirochaetes 2.00 2.00 2.00 2.00 2.00 2.00 - 2.00 2 - 2 11 Tenericutes 1.23 1.23 1.23 1.23 1.67 – 1.10 2.00 – 1.00 2 - 1 13 C. Thermoplasmatota 1.00 1.00 1.00 – 1.00 1.00 – 1.00 1.00 – 1.00 1.00 – 1.00 1 - 1 7 Crenarchaeota 1.10 1.10 1.25 – 1.00 1.25 – 1.00 1.29 – 1.00 2.00 – 1.00 5 - 1 43 Euryarchaeota 2.29 2.67 – 1.00 2.89 – 1.00 4.00 – 1.00 4.00 – 1.00 4.00 – 1.00 5 - 1 138 Thaumarchaeota 1.00 1.00 1.00 1.00 1.00 1.00 – 1.00 1 - 1 3 Mean number of intragenomic 16S rRNA gene variants Taxonomy level Phylum Phylum Class Order Family Genera Species Strain No. genomes Actinobacteria 1.54 1.56 – 1.20 2.00 – 1.00 3.00 – 1.00 4.00 – 1.00 4.00 – 1.00 4 - 1 91 Bacteroidetes 1.77 2.50 – 1.61 2.50 – 1.61 2.50 – 1.00 2.50 – 1.00 5.00 – 1.00 5 - 1 24 C. Saccharibacteria 1.00 1.00 1.00 1.00 1.00 1.00 1 1 Chlamydiae 1.00 1.00 1.00 1.00 1.00 1.00 1 - 1 5 Chlorobi 1.00 1.00 1.00 1.00 1.00 1.00 1 1 Chloroflexi 1.50 2.00 – 1.00 2.00 – 1.00 2.00 – 1.00 2.00 – 1.00 2.00 – 1.00 2 - 1 2 Firmicutes 3.18 3.50 – 1.75 4.85 – 1.75 8.00 – 1.00 9.00 – 1.00 9.00 – 1.00 10 - 1 177 Fusobacteria 3.55 3.55 3.55 3.73 – 3.00 3.73 – 3.00 5.00 – 1.00 5 - 1 21 Ignavibacteriae 1.00 1.00 1.00 1.00 – 1.00 1.00 – 1.00 1.00 – 1.00 1 - 1 2 Proteobacteria 2.87 3.55 – 1.00 4.93 – 1.00 5.45 – 1.00 6.17 – 1.00 8.00 – 1.00 8 - 1 170 Spirochaetes 1.36 1.36 1.36 1.36 1.36 2.00 – 1.22 2 - 1 11 Tenericutes 1.15 1.15 1.15 1.15 1.67 – 1.00 1.67 – 1.00 2 - 1 13 C. Thermoplasmatota 1.00 1.00 1.00 – 1.00 1.00 – 1.00 1.00 – 1.00 1.00 – 1.00 1 - 1 7 Crenarchaeota 1.00 1.00 1.00 – 1.00 1.00 – 1.00 1.00 – 1.00 1.00 – 1.00 1 - 1 43 Euryarchaeota 1.61 1.86 – 1.00 2.00 – 1.00 3.00 – 1.00 3.00 – 1.00 3.00 – 1.00 3 - 1 138 Thaumarchaeota 1.00 1.00 1.00 1.00 1.00 1.00 – 1.00 1 - 1 3 Ranges at the strain level are not mean values, they correspond to the maximum and minimum numbers of intragenomic genes in all strains from a given phylum. C. Saccharibacteria= Candidatus Saccharibacteria; C. Thermoplasmatota= Candidatus Thermoplasmatota; No.= number. Objective 2 251 2.4.2. Evaluation of the primer pairs taken from our previous research and those used most in oral microbiome studies Tables 3 and 4 detail the size and number of 16S rRNA gene amplicons detected by the primer pairs in the oral-bacteria and oral-archaea genomes. The mean number of 16S rRNA gene amplicons varied from 4.84 to 4.39 for bacteria (mean amplicon variants/genome= 2.69 to 1.09) and 2.43 to 1.58 for archaea (mean amplicon variants/genome= 1.34 to 1.08). All the primer combinations identified the maximum mean numbers of intragenomic genes for the bacterial and archaeal species examined (10.83 and 4.00, respectively). However, although most of the primer pairs were able to detect the highest mean value of the gene variants/genome for the archaeal species (i.e., 3.00), only one primer pair detected this maximum value for bacterial species (i.e., 9.00). ALBA REGUEIRA IGLESIAS 252 Table 3. Size and number of 16S rRNA gene amplicons detected by the primer pairs in the oral-bacteria genomes. Superkingdom level Species level ALC Bacteria-specific primer pairs Gene region Amplicon length (mean, bps) g/G (mean) gv/G (mean) Amplicon length (mean, range, bps) g/G (mean, range) gv/G (mean, range) S KP_F048-OP_R043 3 - 4 182.99 4.62 1.25 190.00 – 162.00 10.83 – 1.00 4.00 – 1.00 OP_F098-OP_R119 4 - 5 288.86 4.68 1.22 290.00 - 287.98 10.83 – 1.00 3.00 – 1.00 OP_F066-KP_R040 5 - 6 142.07 4.84 1.11 152.00 – 135.00 10.83 – 1.00 2.00 – 1.00 OP_F009-OP_R030 5 - 7 296.13 4.60 1.38 307.00 – 283.00 10.83 – 1.00 5.00 – 1.00 KP_F061-KP_R074 6 - 7 206.30 4.76 1.31 212.00 – 202.00 10.83 – 1.00 4.00 – 1.00 OP_F101-OP_R030 6 - 7 164.06 4.77 1.31 170.00 – 160.00 10.83 – 1.00 3.00 – 1.00 M OP_F053-KP_R020 1 - 3 351.74 4.61 1.97 547.00 – 315.00 10.83 – 1.00 8.00 – 1.00 KP_F048-KP_R031 3 – 5 454.84 4.61 1.42 462.00 – 433.00 10.83 – 1.00 5.00 – 1.00 KP_F048-OP_R073 3 - 6 546.81 4.67 1.49 554.20 – 520.00 10.83 – 1.00 5.00 – 1.00 KP_F051-KP_R041 4 - 6 410.90 4.84 1.32 421.00 – 404.00 10.83 – 1.00 3.00 – 1.00 KP_F051-OP_R030 4 - 7 566.17 4.62 1.55 577.00 – 552.00 10.83 – 1.00 5.00 – 1.00 OP_F116_KP_R060 7 - 9 308.84 4.54 1.35 1012.50 – 285.00 10.83 – 1.00 4.00 – 1.00 L KP_F048-OP_R030 3 - 7 733.00 4.61 1.71 742.56 – 707.00 10.83 – 1.00 5.00 – 1.00 KP_F048-KP_R060 3 – 9 1059.48 4.59 1.93 1070.00 – 1016.00 10.83 – 1.00 6.00 – 1.00 KP_F056-KP_R077 4 – 9 846.21 4.67 1.81 1551.50 – 821.00 10.83 – 1.00 6.00 – 1.00 Superkingdom level Species level ALC Bacterial and archaeal primer pairs Gene region Amplicon length (mean, bps) g/G (mean) gv/G (mean) Amplicon length (mean, range, bps) g/G (mean, range) gv/G (mean, range) S OP_F114-KP_R002 3 - 4 188.27 4.74 1.27 194.50 – 167.00 10.83 – 1.00 4.00 – 1.00 KP_F020-KP_R032 4 - 5 283.86 4.71 1.22 285.00 – 282.98 10.83 – 1.00 3.00 – 1.00 OP_F066-OP_R073 5 - 6 110.11 4.65 1.09 120.00 – 101.00 10.83 – 1.00 2.00 – 1.00 M OP_F114-KP_R031 3 - 5 456.84 4.60 1.42 464.00 – 435.00 10.83 – 1.00 5.00 – 1.00 OP_F114-OP_R073 3 - 6 548.82 4.67 1.49 556.20 – 522.00 10.83 – 1.00 5.00 – 1.00 KP_F020-OP_R073 4 - 6 375.93 4.72 1.29 386.00 – 366.00 10.83 – 1.00 3.00 – 1.00 L OP_F114-OP_R121 3 - 9 1060.47 4.59 1.93 1071.00 – 1017.00 10.83 – 1.00 6.00 – 1.00 KP_F020-OP_R121 4 - 9 889.30 4.70 1.82 1594.50 – 864.00 10.83 – 1.00 6.00 – 1.00 OP_F066-OP_R121 5 - 9 623.05 4.55 1.65 1328.50 – 598.00 10.83 – 1.00 6.00 – 1.00 Superkingdom level Species level ALC Most used primer pairs Gene region Amplicon length (mean, bps) g/G (mean) gv/G (mean) Amplicon length (mean, range, bps) g/G (mean, range) gv/G (mean, range) S KP_F078-OP_R010B+A 4 – 5 291.86 4.71 1.22 293.00 – 290.98 10.83 – 1.00 3.00 – 1.00 M KP_F031-KP_R021B 1 - 4 525.59 4.39 2.03 700.00 – 467.00 10.83 – 1.00 8.00 – 1.00 KP_F047-KP_R035B 3 - 5 460.25 4.61 1.42 467.00 – 438.00 10.83 – 1.00 5.00 – 1.00 OP_F009-OP_R029B 5 - 8 409.87 4.76 1.51 417.00 – 395.00 10.83 – 1.00 5.00 – 1.00 L KP_F034-KP_R065B 1 – 9* 1505.85 4.81 2.69 1677.00 – 1429.00 10.83 – 1.00 9.00 – 1.00 The amplicon length category and gene regions were determined in objective 1 according to the mean size of the amplicons generated by a given primer and to the mode first position of the forward primer and the mode last of the reverse primer, respectively. The most commonly used primer pairs in the literature were detected in objective 1 (21). *The primer pair gene region was 1-10 when considering the mode positions obtained in the analysis against the oral-bacteria and oral-archaea databases in objective 1 (21). A= archaea; ALC= amplicon length category; B= bacteria; bps= base pairs; F= forward; g/G= number of 16S rRNA gene amplicons per genome; gv/G= number of 16S rRNA gene variant amplicons per genome; KP= Klindworth primer; L= long mean amplicon length category, >600 base pairs; M= medium mean amplicon length category, 301-600 base pairs; OP= oral primer; R= reverse; S= short mean amplicon length category,100-300 base pairs. Objective 2 253 Table 4. Size and number of 16S rRNA gene amplicons detected by the primer pairs in the oral-archaea genomes. Superkingdom level Species level ALC Archaea-specific primer pairs Gene region Amplicon length (mean, bps) g/G (mean) gv/G (mean) Amplicon length (mean, range, bps) g/G (mean, range) gv/G (mean, range) S KP_F018-KP_R002 4* 154.68 1.96 1.09 860.00 – 138.00 4.00 – 1.00 3.00 – 1.00 OP_F066-KP_R013 5 - 6 274.32 1.99 1.11 277.00 – 274.00 4.00 – 1.00 2.00 – 1.00 M KP_F018-KP_R032 3 - 5 429.16 1.95 1.18 1914.00 – 407.00 4.00 – 1.00 3.00 – 1.00 KP_F018-OP_R073 3 - 5 526.11 1.98 1.22 2010.00 – 503.00 4.00 – 1.00 3.00 – 1.00 KP_F020-KP_R013 3 - 6 545.97 1.99 1.21 1325.00 – 539.00 4.00 – 1.00 3.00 – 1.00 KP_F022-KP_R063 5 - 9 586.51 2.00 1.22 670.00 – 530.00 4.00 – 1.00 3.00 – 1.00 L OP_F114-KP_R013 3 - 6 693.70 1.99 1.26 2178.00 – 671.00 4.00 – 1.00 3.00 – 1.00 KP_F018-KP_R063 3 - 9 1131.87 1.96 1.34 1828.00 – 1056.00 4.00 – 1.00 3.00 – 1.00 OP_F066-OP_R016 5 - 9 623.27 2.01 1.22 626.33 – 620.00 4.00 – 1.00 3.00 – 1.00 Superkingdom level Species level ALC Bacterial and archaeal primer pairs Gene region Amplicon length (mean, bps) g/G (mean) gv/G (mean) Amplicon length (mean, range, bps) g/G (mean, range) gv/G (mean, range) S OP_F114-KP_R002 3 - 4 162.29 1.95 1.10 868.00 – 146.00 4.00 – 1.00 3.00 – 1.00 KP_F020-KP_R032 4 - 5 289.43 1.95 1.14 1069.00 – 283.00 4.00 – 1.00 3.00 – 1.00 OP_F066-OP_R073 5 - 6 114.03 1.98 1.08 115.00 – 114.00 4.00 – 1.00 2.00 – 1.00 M OP_F114-KP_R031 3 - 5 436.72 1.95 1.19 1922.00 – 415.00 4.00 – 1.00 3.00 – 1.00 OP_F114-OP_R073 3 - 6 533.62 1.98 1.23 2018.00 – 511.00 4.00 – 1.00 3.00 – 1.00 KP_F020-OP_R073 4 - 6 385.76 1.99 1.18 1165.00 – 379.00 4.00 – 1.00 3.00 – 1.00 L OP_F114-OP_R121 3 - 9 1042.81 1.98 1.33 1741.00 – 1023.00 4.00 – 1.00 3.00 – 1.00 KP_F020-OP_R121 4 - 9 907.41 1.98 1.29 2279.00 – 891.00 4.00 – 1.00 3.00 - 1.00 OP_F066-OP_R121 5 - 9 635.78 1.98 1.22 1323.00 – 625.00 4.00 – 1.00 3.00 – 1.00 Superkingdom level Species level ALC Most used primer pairs Gene region Amplicon length (mean, bps) g/G (mean) gv/G (mean) Amplicon length (mean, range, bps) g/G (mean, range) gv/G (mean, range) S KP_F078-OP_R010B+A 4 – 5 292.79 2.43 1.21 294.00 – 291.50 4.00 – 1.00 3.00 – 1.00 L KP_F014-KP_R011A 3 - 6 606.18 1.58 1.14 603.00 – 608.00 4.00 – 1.00 3.00 – 1.00 The amplicon length category and gene regions were determined in objective 1 according to the mean size of the amplicons generated by a given primer and to the mode first position of the forward primer and the mode last of the reverse primer, respectively. The most commonly used primer pairs in the literature were detected in objective 1 (21). *The primer pair gene region was 3 when considering the mode positions obtained in the analysis against the oral-bacteria and oral-archaea databases in objective 1 (21). A= archaea; ALC= amplicon length category; B= bacteria; bps= base pairs; F= forward; g/G= number of 16S rRNA gene amplicons per genome; gv/G= number of 16S rRNA gene variant amplicons per genome; KP= Klindworth primer; L= long mean amplicon length category, >600 base pairs; M= medium mean amplicon length category, 301-600 base pairs; OP= oral primer; R= reverse; S= short mean amplicon length category,100-300 base pairs. Tables 5 and 6 show the percentages of detected taxa with and without MAs and overabundance estimators obtained by the primer pairs tested on the oral-bacteria and oralarchaea genomes. Our selected primer pairs detected 16S rRNA gene amplicons in a range from 99.46% to 88.71% for the bacterial species and 99.26% to 90.37% for the archaeal species; ALBA REGUEIRA IGLESIAS 254 these percentages were lower for the primer pairs used most in the oral microbiome literature (95.16% - 74.19% for the bacteria, and 63.70% and 30.37% for the archaea). Overall, excluding the most commonly used primer pairs in the literature, unlike the coverage values, the SC-NMA values increased as the mean length of the amplicons obtained by the primer pair increased. If we contrast the percentages of species detected with their respective SC-NMA, all the short primer pairs analysed showed the largest differences between both parameters (average difference= 21.34% for bacteria and 23.70% for archaea), followed by those of medium length (7.30% and 13.75%, respectively). The long primer pairs presented the smallest differences between the coverage and SC-NMA values (4.30% and 5.82%, respectively). According to the SC-NMA values, the best three bacteria-specific primer pairs were: KP_F048-OP_R030 and KP_F048-KP_R060 (L, SC-NMA= 93.55%, six MAs, OF-MA= 1.06 for both) and OP_F053-KP_R020 (M, 93.01%, six, 1.06). In contrast, the worst primer pair was OP_F066-KP_R040 (S, 47.31%, 77, 2.78). The most commonly used primer pairs in the literature did not stand out for their SC-NMA values among those in their category. Considering the three categories of amplicon lengths, the SC-NMA values for the archaeaspecific primer pairs ranged from 89.63% for the KP_F018-KP_R063 (L, six MAs, OF-MA= 1.11), 85.93% for KP_F022-KP_R063 (M, eight, 1.14) to 69.63% for the OP_F066-KP_R013 (S, 35, 1.99). Interestingly, the long primer pair KP_F014-KP_R011, which is the one used most in the literature to detect oral archaea, was only able to identify 30.37% of the species tested in this study, resulting in the lowest SC-NMA value (26.67%, five MAs, OF-MA= 1.14). In relation to the bacterial and archaeal primer pairs, the overall SC-NMA values ranged from 92.52% for OP_F114_OP_R121 (L, 12 MAs, OF-MA= 1.08), 88.79% for OP_F114- KP_R031 (M, 29, 1.26) to 54.21% for OP_F066-OP_R073 (S, 134, 3.45). In terms of overall SC-NMA, the second worst was KP_F078-OP_R010 (S, 66.67%, 48 MAs, OF-MA= 1.68), mainly due to its low capacity to detect archaea (63.70%), which directly affected the SC-NMA value for archaea (48.89%) (Table 6). However, this primer pair is the most widely used in the literature to detect bacteria and archaea. Objective 2 255 Table 5. Detected taxa with and without matching amplicons and overabundance estimators obtained by the primer pairs tested on the oral-bacteria genomes. ALC Bacteria-specific primer pairs Detected genomes ( % ) Detected species ( % ) Detected species with MAs ( % ) SC-NMA (%) OF OF-MA S KP_F048-OP_R043 508 (98.07) 180 (96.77) 22 (12.22) 158 (84.95) 5.82 1.30 OP_F098-OP_R119 493 (95.17) 177 (95.16) 28 (15.82) 149 (80.11) 8.28 1.73 OP_F066-KP_R040 455 (87.84) 165 (88.71) 77 (46.67) 88 (47.31) 13.16 2.78 OP_F009-OP_R030 504 (97.30) 181 (97.31) 29 (16.02) 152 (81.72) 6.01 1.32 KP_F061-KP_R074 468 (90.35) 169 (90.86) 39 (23.08) 130 (69.89) 7.48 1.61 OP_F101-OP_R030 460 (88.80) 167 (89.78) 39 (23.35) 128 (68.82) 7.44 1.61 M OP_R053-KP_R020 506 (97.68) 179 (96.24) 6 (3.35) 173 (93.01) 4.80 1.06 KP_F048-KP_R031 507 (97.88) 180 (96.77) 9 (5.00) 171 (91.94) 5.04 1.12 KP_F048-OP_R073 498 (96.14) 178 (95.70) 6 (3.37) 172 (92.47) 4.72 1.06 KP_F051-KP_R041 456 (88.03) 166 (89.25) 20 (12.05) 146 (78.50) 7.45 1.58 KP_F051-OP_R030 508 (98.07) 184 (98.92) 19 (10.33) 165 (88.71) 5.53 1.22 OP_F116_KP_R060 516 (99.61) 185 (99.46) 31 (16.76) 154 (82.80) 6.13 1.37 L KP_F048-OP_R030 507 (97.88) 180 (96.77) 6 (3.33) 174 (93.55) 4.72 1.06 KP_F048-KP_R060 507 (97.88) 180 (96.77) 6 (3.33) 174 (93.55) 4.72 1.06 KP_F056-KP_R077 495 (95.56) 180 (96.77) 10 (5.56) 170 (91.40) 4.89 1.10 ALC Bacterial and archaeal primer pairs Detected genomes (%) Detected species (%) Detected species with MAs (%) SC-NMA (%) OF OF-MA S OP_F114-KP_R002 485 (93.63) 172 (92.47) 22 (12.79) 150 (80.65) 5.82 1.30 KP_F020-KP_R032 488 (94.21) 176 (94.62) 28 (15.91) 148 (79.57) 8.28 1.73 OP_F066-OP_R073 502 (96.91) 182 (97.85) 85 (46.70) 97 (52.15) 15.90 3.31 M OP_F114-KP_R031 507 (97.88) 180 (96.77) 9 (5.00) 171 (91.94) 5.04 1.12 OP_F114-OP_R073 498 (96.14) 178 (95.70) 6 (3.77) 172 (92.47) 4.72 1.06 KP_F020-OP_R073 488 (94.21) 176 (94.62) 22 (12.50) 154 (82.80) 7.52 1.61 L OP_F114-OP_R121 507 (97.88) 180 (96.77) 6 (3.33) 174 (93.55) 4.72 1.06 KP_F020-OP_R121 489 (94.40) 177 (95.16) 10 (5.65) 167 (89.79) 4.89 1.10 OP_F066-OP_R121 516 (99.61) 185 (99.46) 16 (8.65) 169 (90.86) 5.20 1.16 ALC Most used primer pairs Detected genomes ( % ) Detected species ( % ) Detected species with MAs ( % ) SC-NMA (%) OF OF-MA S KP_F078-OP_R010B+A 488 (94.21) 176 (94.62) 28 (15.91) 148 (79.57) 8.28 1.73 M KP_F031-KP_R021B 347 (66.99) 138 (74.19) 2 (1.45) 136 (73.12) 4.50 1.02 KP_F047-KP_R035B 500 (96.53) 177 (95.16) 9 (5.09) 168 (90.32) 5.04 1.12 OP_F009-OP_R029B 469 (90.54) 164 (88.17) 24 (14.63) 140 (75.27) 5.62 1.24 L KP_F034-KP_R065B 440 (84.94) 155 (83.33) 2 (1.29) 153 (82.26) 4.50 1.02 A= archaea; ALC= amplicon length category; B= bacteria; F= forward; KP= Klindworth primer; L= long mean amplicon length category, >600 base pairs; M= medium mean amplicon length category, 301-600 base pairs; MAs= matching amplicons; OF= overestimation factor; OF-MA= overestimation factor associated with matching amplicons; OP= oral primer; R= reverse; S= short mean amplicon length category,100-300 base pairs; SC-NMA= species coverage with no matching amplicons. ALBA REGUEIRA IGLESIAS 256 Table 6. Detected taxa with and without matching amplicons and overabundance estimators obtained by the primer pairs tested on the oral-archaea genomes. ALC Archaea-specific primer pairs Detected genomes ( % ) Detected species ( % ) Detected species with MAs ( % ) SC-NMA (%) OF OF-MA S KP_F018-KP_R002 185 (96.86) 129 (95.56) 29 (22.48) 100 (74.07) 3.30 1.76 OP_F066-KP_R013 184 (96.34) 129 (95.56) 35 (27.13) 94 (69.63) 4.02 1.99 M KP_F018-KP_R032 186 (97.38) 130 (96.30) 20 (15.39) 110 (81.48) 2.68 1.49 KP_F018-OP_R073 177 (92.67) 122 (90.37) 18 (14.75) 104 (77.04) 2.65 1.46 KP_F020-KP_R013 183 (95.81) 128 (94.81) 20 (15.63) 108 (80.00) 2.61 1.35 KP_F022-KP_R063 180 (94.24) 124 (91.85) 8 (6.45) 116 (85.93) 2.26 1.14 L OP_F114-KP_R013 184 (96.34) 129 (95.56) 16 (12.40) 113 (83.70) 2.47 1.28 KP_F018-KP_R063 183 (95.81) 127 (94.07) 6 (4.72) 121 (89.63) 2.16 1.11 OP_F066-OP_R016 180 (94.24) 124 (91.85) 8 (6.45) 116 (85.93) 2.26 1.14 ALC Bacterial and archaeal primer pairs Detected genomes (%) Detected species (%) Detected species with MAs (%) SC-NMA (%) OF OF-MA S OP_F114-KP_R002 190 (99.48) 134 (99.26) 29 (21.64) 105 (77.78) 3.30 1.76 KP_F020-KP_R032 190 (99.48) 134 (99.26) 30 (22.39) 104 (77.04) 3.88 1.92 OP_F066-OP_R073 181 (94.76) 126 (93.33) 49 (38.89) 77 (57.04) 8.37 3.71 M OP_F114-KP_R031 190 (99.48) 134 (99.26) 20 (14.93) 114 (84.44) 2.68 1.49 OP_F114-OP_R073 181 (94.76) 126 (93.33) 18 (14.29) 108 (80.00) 2.65 1.46 KP_F020-OP_R073 180 (94.24) 125 (92.59) 26 (20.80) 99 (73.33) 3.74 1.85 L OP_F114-OP_R121 185 (96.86) 129 (95.56) 6 (4.65) 123 (91.11) 2.16 1.11 KP_F020-OP_R121 185 (96.86) 129 (95.56) 6 (4.65) 123 (91.11) 2.16 1.11 OP_F066-OP_R121 186 (97.38) 130 (96.30) 8 (6.15) 122 (90.37) 2.26 1.14 ALC Most used primer pairs Detected genomes ( % ) Detected species ( % ) Detected species with MAs ( % ) SC-NMA (%) OF OF-MA S KP_F078-OP_R010B+A 123 (66.40) 86 (63.70) 20 (23.26) 66 (48.89) 3.56 1.60 L KP_F014-KP_R011A 44 (23.04) 41 (30.37) 5 (12.20) 36 (26.67) 2.00 1.14 A= archaea; ALC= amplicon length category; B= bacteria; F= forward; KP= Klindworth primer; L= long mean amplicon length category, >600 base pairs; M= medium mean amplicon length category, 301-600 base pairs; MAs= matching amplicons; OF= overestimation factor; OF-MA= overestimation factor associated with matching amplicons; OP= oral primer; R= reverse; S= short mean amplicon length category,100-300 base pairs; SC-NMA= species coverage with no matching amplicons. Appendices S5-S11 contain more detailed information on the results of MA- and MA-free species coverage and the overabundance parameters (OF and OF-MA values) obtained by the primer pairs tested against the oral-bacteria and oral-archaea genomes. Appendices S5, S9, and S11 also include the results obtained by the bacterial and archaeal primer pairs for both domains. Objective 2 263 2.6. CONCLUSIONS In conclusion, nearly all oral bacteria and about half of the oral archaea have more than one 16S rRNA gene in their respective genomes. Depending on the primer pair used, up to almost half of the species present MAs, affecting relevant genera present in the oral environment such as Actinomyces, Fusobacterium, Lactobacillus, Methanosarcina, Staphylococcus, and Streptococcus. The performance of the primer pairs to detect non-MA species increases as the average length of the amplicons increases; none of these being the most widely used primer pairs in the oral microbiome literature. The best primer pairs were: KP_F048-OP_R030 (for bacteria; region 3-7; primer pair position for Escherichia coli J01859.1: 342-1079), KP_F018-KP_R063 (for archaea; 3-9; undefined-1506), and OP_F114_OP_R121 (for both bacteria and archaea; 3-9; 340-1405). In addition to the 16S rRNA gene redundancy, the considerable presence of MAs must be controlled to ensure the accurate interpretation of microbial diversity data. The SC-NMA is a more useful parameter than the conventional coverage percentage for selecting the best primer pairs. The choice of primer pair affects significantly diversity estimates and taxonomic classification, conditioning the comparability of oral microbiome studies using different primer pairs. ALBA REGUEIRA IGLESIAS 264 2.7. REFERENCES (1) Rajendhran J, Gunasekaran P. Microbial phylogeny and diversity: small subunit ribosomal RNA sequence analysis and beyond. Microbiol Res. 2011 Feb;166(2):99-110. (2) Woese CR. Bacterial evolution. Microbiol Rev. 1987 Jun;51(2):221-71. (3) del Rosario-Rodicio M, del Carmen-Mendoza M. Bacterial identification by 16S rRNA sequencing: rationale, methodology and applications in clinical microbiology. Enferm Infecc Microbiol Clin. 2004 Apr;22(4):238-45. Spanish. (4) Acinas SG, Marcelino LA, Klepac-Ceraj V, Polz MF. Divergence and redundancy of 16S rRNA sequences in genomes with multiple rrn operons. J Bacteriol. 2004 May;186(9):2629- 35. (5) Pei AY, Oberdorf WE, Nossa CW, Agarwal A, Chokshi P, Gerz EA, et al. Diversity of 16S rRNA genes within individual prokaryotic genomes. Appl Environ Microbiol. 2010 Jun;76(12):3886-97. (6) Sun D, Jiang X, Wu QL, Zhou N. Intragenomic heterogeneity of 16S rRNA genes causes overestimation of prokaryotic diversity. Appl Environ Microbiol. 2013 Oct;79(19):5962-9. (7) Větrovský T, Baldrian P. The variability of the 16S rRNA gene in bacterial genomes and its consequences for bacterial community analyses. PLoS One. 2013;8(2):e57923. doi: 10.1371/journal.pone.0057923. (8) Case RJ, Boucher Y, Dahllöf I, Holmström C, Doolittle WF, Kjelleberg S. Use of 16S rRNA and rpoB genes as molecular markers for microbial ecology studies. Appl Environ Microbiol. 2007 Jan;73(1):278-88. (9) Lee ZM, Bussema C 3rd, Schmidt TM. rrnDB: documenting the number of rRNA and tRNA genes in bacteria and archaea. Nucleic Acids Res. 2009 Jan;37(Database issue):D489-93. doi: 10.1093/nar/gkn689. Objective 2 265 (10) Johnson JS, Spakowicz DJ, Hong BY, Petersen LM, Demkowicz P, Chen L, et al. Evaluation of 16S rRNA gene sequencing for species and strain-level microbiome analysis. Nat Commun. 2019 Nov;10(1):5029. doi: 10.1038/s41467-019-13036-1. (11) Chen J, Miao X, Xu M, He J, Xie Y, Wu X, et al. Intra-genomic heterogeneity in 16S rRNA genes in strictly anaerobic clinical isolates from periodontal abscesses. PLoS One. 2015 Jun;10(6):e0130265. doi: 10.1371/journal.pone.0130265. (12) Durán-Pinedo AE, Frias-Lopez J. Beyond microbial community composition: functional activities of the oral microbiome in health and disease. Microbes Infect. 2015 Jul;17(7):505-16. (13) Verma D, Garg PK, Dubey AK. Insights into the human oral microbiome. Arch Microbiol. 2018 May;200(4):525-40. (14) F Escapa I, Huang Y, Chen T, Lin M, Kokaras A, Dewhirst FE, et al. Construction of habitat-specific training sets to achieve species-level assignment in 16S rRNA gene datasets. Microbiome. 2020 May;8(1):65. doi: 10.1186/s40168-020-00841-w. (15) Relvas M, Regueira-Iglesias A, Balsa-Castro C, Salazar F, Pacheco JJ, Cabral C, et al. Relationship between dental and periodontal health status and the salivary microbiome: bacterial diversity, co-occurrence networks and predictive models. Sci Rep. 2021 Jan;11(1):929. doi: 10.1038/s41598-020-79875-x. (16) Camelo-Castillo A, Novoa L, Balsa-Castro C, Blanco J, Mira A, Tomas I. Relationship between periodontitis-associated subgingival microbiota and clinical inflammation by 16S pyrosequencing. J Clin Periodontol. 2015 Dec;42(12):1074-82. (17) Camelo-Castillo AJ, Mira A, Pico A, Nibali L, Henderson B, Donos N, et al. Subgingival microbiota in health compared to periodontitis and the influence of smoking. Front Microbiol. 2015 Feb;6:119. doi: 10.3389/fmicb.2015.00119. ALBA REGUEIRA IGLESIAS 266 (18) F Escapa I, Chen T, Huang Y, Gajare P, Dewhirst FE, Lemon KP. New insights into human nostril microbiome from the expanded Human Oral Microbiome Database (eHOMD): a resource for the microbiome of the human aerodigestive tract. mSystems. 2018 Dec;3(6):e00187-18. doi: 10.1128/mSystems.00187-18. (19) Clark K, Karsch-Mizrachi I, Lipman DJ, Ostell J, Sayers EW. GenBank. Nucleic Acids Res. 2016 Jan;44:D67-72. doi: 10.1093/nar/gkv1276. (20) NCBI Resource Coordinators. Database resources of the National Center for Biotechnology Information. Nucleic Acids Res. 2016 Jan;44(D1):D7-19. doi: 10.1093/nar/gkv1290. (21) Regueira-Iglesias A, Vázquez-González L, Balsa-Castro C, Vila-Blanco N, Blanco-Pintos T, Tamames J, et al. In silico evaluation and selection of the best 16S rRNA gene primers for use in next-generation sequencing to detect oral bacteria and archaea. Accepted for publication in Microbiome. Preprint at Research Square. 2021. doi: 10.21203/rs.3.rs-516961/v1. (22) National Center for Biotechnology Information. Entrez programming utilities help. 2010; Available at: https://www.ncbi.nlm.nih.gov/books/NBK25501/. (23) Python Software Foundation. Python. Version 3.9.0. 2020. http://www.python.org/. (24) Schoch CL, Ciufo S, Domrachev M, Hotton CL, Kannan S, Khovanskaya R, et al. NCBI Taxonomy: a comprehensive update on curation, resources and tools. Database (Oxford). 2020 Jan;2020:baaa062. doi: 10.1093/database/baaa062. (25) O'Leary NA, Wright MW, Brister JR, Ciufo S, Haddad D, McVeigh R, et al. Reference sequence (RefSeq) database at NCBI: current status, taxonomic expansion, and functional annotation. Nucleic Acids Res. 2016 Jan;44(D1):D733-45. doi: 10.1093/nar/gkv1189. (26) Lyalina S. Search 16S py algorithm. 2019; Available at: https://github.com/slyalina/search_16S_py. Objective 2 267 (27) Edgar R. SEARCH_16S: A new algorithm for identifying 16S ribosomal RNA genes in contigs and chromosomes. Preprint at bioRxiv 2017:124131. doi: 10.1101/124131. (28) Harris CR, Millman KJ, van der Walt, Stéfan J, Gommers R, Virtanen P, Cournapeau D, et al. Array programming with NumPy. Nature. 2020 Sep;585(7825):357-62. (29) McKinney W. Data structures for statistical computing in Python. In: van der Walt S, Millman J, editors. Proceedings of the 9th Python in Science Conference; 2010; Austin. Texas: SciPy; 2010. doi: 10.25080/Majora-92bf1922-00a. (30) Baker GC, Smith JJ, Cowan DA. Review and re-analysis of domain-specific 16S primers. J Microbiol Methods. 2003 Dec;55(3):541-55. (31) Barnett M. regex. 2020; Available at: https://pypi.org/. (32) Coenye T, Vandamme P. Intragenomic heterogeneity between multiple 16S ribosomal RNA operons in sequenced bacterial genomes. FEMS Microbiol Lett. 2003 Nov;228(1):45-9. (33) Klappenbach JA, Saxman PR, Cole JR, Schmidt TM. rrnDB: the ribosomal RNA operon copy number database. Nucleic Acids Res. 2001 Jan;29(1):181-4. (34) Stoddard SF, Smith BJ, Hein R, Roller BR, Schmidt TM. rrnDB: improved tools for interpreting rRNA gene abundance in bacteria and archaea and a new foundation for future development. Nucleic Acids Res. 2015 Jan;43(Database issue):D593-8. doi: 10.1093/nar/gku1201. (35) Applied Maths NV. Kodon 2.0; Available at: https://www.applied-maths.com. (36) Lagesen K, Hallin P, Rødland EA, Stærfeldt H, Rognes T, Ussery DW. RNAmmer: consistent and rapid annotation of ribosomal RNA genes. Nucleic Acids Res. 2007;35(9):3100- 8. ALBA REGUEIRA IGLESIAS 268 (37) Teles R, Teles F, Frias-Lopez J, Paster B, Haffajee A. Lessons learned and unlearned in periodontal microbiology. Periodontol 2000. 2013 Jun;62(1):95-162. (38) Abranches J, Zeng L, Kajfasz JK, Palmer SR, Chakraborty B, Wen ZT, et al. Biology of oral streptococci. Microbiol Spectr. 2018 Oct;6(5):10.1128/microbiolspec.GPP3-0042-2018. (39) Jiang S, Gao X, Jin L, Lo EC. Salivary microbiome diversity in caries-free and cariesaffected children. Int J Mol Sci. 2016 Nov;17(12):1978. doi: 10.3390/ijms17121978. (40) Mitchell J. Streptococcus mitis: walking the line between commensalism and pathogenesis. Mol Oral Microbiol. 2011 Apr;26(2):89-98. (41) Thurnheer T, Belibasakis GN. Streptococcus oralis maintains homeostasis in oral biofilms by antagonizing the cariogenic pathogen Streptococcus mutans. Mol Oral Microbiol. 2018 Jun;33(3):234-9. (42) Deng ZL, Szafrański SP, Jarek M, Bhuju S, Wagner-Döbler I. Dysbiosis in chronic periodontitis: key microbial players and interactions with the human host. Sci Rep. 2017 Jun;7(1):3703. doi: 10.1038/s41598-017-03804-8. (43) Schloss PD. Amplicon sequence variants artificially split bacterial genomes into separate clusters. mSphere. 2021 Aug;6(4):e0019121. doi: 10.1128/mSphere.00191-21. (44) Zhang J, Ding X, Guan R, Zhu C, Xu C, Zhu B, et al. Evaluation of different 16S rRNA gene V regions for exploring bacterial diversity in a eutrophic freshwater lake. Sci Total Environ. 2018 Mar;618:1254-67. (45) Caporaso JG, Lauber CL, Walters WA, Berg-Lyons D, Lozupone CA, Turnbaugh PJ, et al. Global patterns of 16S rRNA diversity at a depth of millions of sequences per sample. Proc Natl Acad Sci U S A. 2011 Mar;108 Suppl 1(Suppl 1):4516-22. Objective 2 269 (46) Illumina, Inc. 16S Metagenomic Sequencing Library Preparation. 2013; Available at: https://support.illumina.com/content/dam/illuminasupport/documents/documentation/chemistry_documentation/16s/16s-metagenomic-library- prep-guide-15044223-b.pdf. (47) Zhang Y, Wang X, Li H, Ni C, Du Z, Yan F. Human oral microbiota and its modulation for oral health. Biomed Pharmacother. 2018 Mar;99:883-93. OBJECTIVE 3 Objective 3 279 Therefore, a total of 709 complete prokaryotic genomes from a total of 321 oral species were downloaded (more than one complete genome was analysed for several species). Finally, the complete taxonomic hierarchy (from superkingdom to strain) of all downloaded complete genomes was designated by the taxonomic identifier included in the annotated information in the NCBI database (26), all computationally performed with our script above. The taxonomy and NCBI identifiers of the oral-bacteria and archaea genomes are included in appendices S1 and S2, respectively. 3.3.2. Selecting the primer pairs and obtaining the in silico amplicons of the 16S rRNA gene Thirty-three primer pairs with the best in silico coverage, as identified in objective 1, were selected, along with the six primer pairs used the most in the oral microbiome literature (27). These primer pairs were classified according to the mean length of their amplicons into: short primer pairs (S, 100-300 base pairs), medium primer pairs (M, 301-600 bps), and long primer pairs (L, >600 bps); and the domain targeted (bacteria, archaea, or both) (Table 1). ALBA REGUEIRA IGLESIAS 280 Table 1. Selected primer pairs with high in silico coverage percentages targeting oral bacteria and/or archaea and those most used in the sequencing-based studies of the oral microbiome. Primer pairs were selected based on the species coverage values (number of species detected /total species evaluated) in objective 1 (27). They were individually evaluated through regular expressions against Escherichia coli J01859 to define their positions. The U values represent a mismatch on the assessment and, therefore, the position cannot be confirmed with a guarantee. Gene regions were delimited as described by Baker et al. (32). *The primer pair gene region was 3 and **1-10 when considering the mode positions obtained in the analysis against the oral-bacteria and oral-archaea databases in objective 1 (27). ALC= amplicon length category; bps= base pairs; F= forward; KP= Klindworth primer; L= long mean amplicon length category, >600 base pairs; M= medium mean amplicon length category, 301-600 base pairs; OP= oral primer; Post= position; R= reverse; S= short mean amplicon length category, 100-300 base pairs; U= unaligned with Escherichia coli; V= variable. Objective 3 281 Applying our script in combination with Python’s regex module (33), the direct and reverse sequences of each primer pair were used to obtain in silico sequence segments of the whole genomes analysed (hereafter referred to as in silico amplicons). An in silico amplicon was considered for subsequent analysis when the following conditions were present: 1) there is a zero mismatch of both primers (forward and reverse) of each pair; 2) the distance between the starting position of the forward primer and the ending position of the reverse primer is less than 2300 and higher than 100 nucleotides; 3) the in silico amplicon does not repeat within the same species. All in silico amplicons from the same species, even if from different strains, were considered for analysis. For each primer pair, a fasta file was created where all in silico amplicons found were stored. The stored sequences were identified with the same taxonomic hierarchy as the genomes from which they were detected. As many in silico-amplicons were detected within the same species (differing by at least 1 nucleotide), the sequence variants were identified with correlative numbering at a new hierarchical level below the species name. In the fasta files, all sequences included a species identifier (SPn) and a variant identifier (Vn) in their header, and then the header of each sequence also included the taxonomic hierarchy up to the variant level within each species. Finally, the in silico amplicons were obtained from 186 oral-bacterial and 135 oral-archaeal species. 3.3.3. Determination of the percentage of similarity between in silico amplicons of different oral species by MegaBLAST A script with the NcbiblastnCommandline wrapper from Biopython (34) was developed to manage BLAST+ 2.11 (35) in the local mode from Biopython. This enabled the data obtained in the alignments to be easily transferred for later analysis on Python (29). The alignment parameters were configured to be the same as the default settings in MegaBLAST (36) since these settings were appropriate for the alignment between sequences with a similarity ≥95%. All sequences belonging to the same fasta file for the same primer pair were aligned against themselves; in order to do this, each fasta file was inserted as subject and query in BLASTN ALBA REGUEIRA IGLESIAS 282 (37) for obtaining the percentage of similarity between in silico amplicons belonging to different oral species. From the results obtained, in silico amplicons with 100% alignment coverage of the query sequences and with a similarity value ≥97% were selected. That is, alignments with the following BLAST+ estimates (35): qcovs= 100%, qcovhsp= 100%, qcovus= 100%, and pident ≥97% were selected. Of the above alignments obtained, the following were discarded as they were not of interest: 1) in silico amplicons with the same unique identifier (SPn + Vn); 2) in silico amplicons with the same species identifier; and 3) duplicate alignments. If two different species had more than one in silico amplicon similarity value ≥97% (ASI97) among them, one of the alignments was chosen at random. The results of the highly similar species pairs, including taxonomic hierarchy data for both species, were then stored using the pandas (38) and xlsxwriter (39) Python modules. 3.3.4. Construction of a matrix with oral species showing in silico amplicon similarity values ≥97% and calculation of descriptive statistical estimators A similarity matrix was created for each primer pair, where rows and columns had the species identifiers, and cells indicated with a number 1 the presence of an ASI97 between two different species. We then developed a script in R (version 4.0.3) (40) through which we calculated the following estimates for each analysed primer pair: 1) the number of species with at least one ASI97 with other species; 2) the total number of ASI97 between different species; 3) the mean and maximum numbers of ASI97 per species. In addition, we estimated the percentage of detected species (species coverage, SC) and the percentage of detected species without ASI97 for each primer pair (species coverage no ASI97, SC-NASI97). This last parameter was then used as a criterion for selecting the primers associated with a smaller number of oral species that may be erroneously clustered. The SC-NASI97 parameter will be influenced not only by the number of species with ASI, but also by the coverage percentages of each primer pair. Finally, the bacterial and archaeal species pairs that showed an ASI97 were described and assessed whether they belonged to different genera or higher taxonomic ranks. Objective 3 283 3.4. RESULTS 3.4.1. Evaluation of the primer pairs for detecting oral species with in silico amplicon similarity values ≥97% The primer pairs that targeted bacteria had a mean of 91.88 (49.40%) bacterial species with an ASI97 and an average of 153.46 ASI97 containing distinct species. For those targeting archaea, these numbers were 65.60 (48.59%) and 162.26, respectively. If the primers used most in the oral microbiome literature were excluded, those with short amplicon lengths (unlike the SC percentages) had the lowest SC-NASI97 values for both bacteria (S= 39.54%) and archaea (S= 40.44%) compared to the medium length and long primers (M= 45.82% and 46.35%, respectively; L= 48.39% and 44.32%, respectively). Figures 2 and 3 show the number of species with ASI97 and the number of ASI97 with each primer pair evaluated against bacteria and archaea, respectively; while figures 4 and 5 detail the percentages of coverage and coverage considering the presence or absence of ASI97 for both domains. Concerning the bacteria-specific primer pairs, the number of bacterial species with an ASI97 and the total number of ASI97 ranged from 37 and 32 with the most widely used primer KP_F031-KP_R021 (M; SC-NASI97= 54.30%), to 120 and 277 with OP_F066- KP_R040 (S; SC-NASI97= 24.19%), respectively. This latter primer also had the lowest SCNASI97 value, while OP_F053-KP_R020 detected the highest number of species with no ASI97 (M; SC-NASI97= 65.05%). In addition, except for OP_F053-KP_R020, all the bacteriaspecific primers had a maximum number of ASI97/species above five (range= 4 - 15 ASI97/species). Concerning the archaea-specific primer pairs, the number of archaeal species with an ASI97 and the total number of ASI97 ranged from 24 and 96 with the widely used KP_F014- KP_R011 (L; SC-NASI97= 12.59%) to 89 and 240 with OP_F066-KP_R013 (S; SC-NASI97= 29.63%), respectively. The former primer detected the lowest number of species without an ASI97, and KP_F018-KP_R002 the highest (S; SC-NASI97= 51.11%). Moreover, all the archaea-specific primers had a maximum number of ASI97/species ≥10 (range= 10 - 13 ASI97/species). ALBA REGUEIRA IGLESIAS 284 Figure 2. Number of bacterial species with in silico amplicon similarity values ≥97% and number of in silico amplicon similarity values ≥97% with the primer pairs evaluated against the oral-bacteria genomes. (A) Estimates were obtained by the selected bacteria-specific primer pairs. (B) Estimates were obtained by the selected bacterial and archaeal primer pairs and the primer pairs used the most in the oral microbiome literature. Among the most commonly used primer pairs in the literature, those marked with an * are bacteria-specific and those with ** target both bacterial and archaea. ASI97= in silico amplicon similarity values ≥97%; F= forward; KP= Klindworth primer; No.= number; OP= oral primer; R= reverse. Finally, using both the bacterial and archaeal primer pairs, the number of bacterial and archaeal species with an ASI97 and the total number of ASI97 ranged from 84 and 60 and 118 and 126, respectively, with OP_F114-KP_R002 (S; SC-NASI≥97= 47.31% for bacteria and 54.81% for archaea) to 124 and 95 and 239 and 286, respectively, with OP_F066-OP_R073 (S; SC-NASI≥97= 31.18% for bacteria and 22.96% for archaea). The latter primer also detected Objective 3 285 the lowest number of species without an ASI97 and OP_F114-KP_R031 the highest (M; SCNASI97= 51.08% for bacteria and 53.33% for archaea) (Figures 2b-5b). Most bacterial and archaeal primer combinations had maximum numbers of ASI97/species ≥10 (range= 9 - 14 ASI/species and 11 - 14 ASI/ species for both domains, respectively). Figure 3. Number of archaeal species with in silico amplicon similarity values ≥97% and number of in silico amplicon similarity values ≥97% with the primer pairs evaluated against the oral-archaea genomes. (A) Estimates were obtained by the selected archaea-specific primer pairs. (B) Estimates were obtained by the selected bacterial and archaeal primer pairs and the primer pairs used the most in the oral microbiome literature. Among the most commonly used primer pairs in the literature, those marked with an * are archaea-specific and those with ** target both bacterial and archaea. ASI97= in silico amplicon similarity values ≥97%; F= forward; KP= Klindworth primer; No.= number; OP= oral primer; R= reverse. ALBA REGUEIRA IGLESIAS 286 Figure 4. Percentages of coverage and coverage considering the species with in silico amplicon similarity values ≥97% of the primer pairs evaluated against the oral-bacteria genomes. (A) Percentages were obtained by the selected bacteria-specific primer pairs. (B) Percentages were obtained by the selected bacterial and archaeal primer pairs and the primer pairs used the most in the oral microbiome literature. Among the most commonly used primer pairs in the literature, those marked with an * are bacteria-specific and those with ** target both bacterial and archaea. ASI97= in silico amplicon similarity values ≥97%; F= forward; KP= Klindworth primer; Non-SC= non-coverage of species; OP= oral primer; R= reverse; SC= species coverage; SCASI97= species coverage with in silico amplicon similarity values ≥97%; SC-NASI97= species coverage with no in silico amplicon similarity values ≥97%. Objective 3 287 Figure 5. Percentages of coverage and coverage considering the species with in silico amplicon similarity values ≥97% of the primer pairs evaluated against the oral-archaea genomes. (A) Percentages obtained by the selected archaea-specific primer pairs. (B) Percentages obtained by the selected bacterial and archaeal primer pairs and the primer pairs used the most in the oral microbiome literature. Among the most commonly used primer pairs in the literature, those marked with an * are archaea-specific and those with ** target both bacterial and archaea. ASI97= in silico amplicon similarity values ≥97%; F= forward; KP= Klindworth primer; Non-SC= non-coverage of species; OP= oral primer; R= reverse; SC= species coverage; SCASI97= species coverage with in silico amplicon similarity values ≥97%; SC-NASI97= species coverage with no in silico amplicon similarity values ≥97%. Figures 6 and 7 are networks showing the potential clusters (hereinafter referred to as potential OTUs) with a ≥97% similarity threshold obtained with the primer pairs that presented the lowest SC-NASI97 values (F066-KP_R040 for bacteria, OP_F066-KP_R013 for archaea, and OP_F066-OP_R073 for bacteria and archaea), as well as one of the most used primer pairs in the oral microbiome literature, KP_F078-OP_R010. Thus, for example, for the primer pair F066-KP_R040, focusing on the one indicated by a dashed dotted line, 24 bacteria formed a ALBA REGUEIRA IGLESIAS 288 potential OTU, in which 10 genera, five families, and two orders were involved. As can be seen, there were species such as Ligilactobacillus salivarius (spp. 162) that presented high similarity only with two others, Lacticaseibacillus paracasei (spp. 196) and Lacticaseibacillus rhamnosus (spp. 243); while Staphylococcus cohnii (spp. 297) presented high similarity with 11 species (among which, Enterococcus faecalis, spp. 155; Staphylococcus aureus, spp. 163; Levilactobacillus brevis, spp. 181; Lentilactobacillus buchneri, spp. 237) belonging to four genera, three families and two orders. For the primer pair OP_F066-KP_R013, the potential OTU indicated was formed by nine archaea, involving seven genera and two families. Thus, Desulfurococcus amylolyticus (spp. 37) showed high similarity only with Desulfurococcus mucosus (spp. 68), while Thermogladius caldera had high similarity with five species (Hyperthermus butylicus, spp. 24; Staphylothermus marinus, spp. 26; Staphylothermus hellenicus, spp. 57; Desulfurococcus mucosus, spp. 68; Pyrolobus fumarii; sp. 82) belonging to four genera and two families. Objective 3 295 3.5. DISCUSSION The high degree of similarity between full-length 16S rRNA sequences from distinct species, or even genera, has been reported in the literature (12,13), leading to questions about the reliability of diversity estimates derived from sequence clustering methods based on a given similarity threshold. Using full-length genes and a ≥97% similarity threshold, some authors have detected that around a quarter of constructed OTUs contain sequences from multiple species (12,13) and about a tenth from distinct genera (12). These estimates were obviously higher when gene regions were assessed instead of full sequences. Schloss et al. (13) found that, with a ≥97% similarity threshold and applying the OptiClust algorithm (42), 31.7%, 34.3% and 34.8% of the OTUs assessed had 16S rRNA amplicons from distinct species in the variable regions 3-4, 4, and 4-5, respectively (13). However, these investigations did not focus on taxa inhabiting a specific environment, despite the importance of conducting 16S rRNA gene-based research using habitat-specific databases (43). Consequently, we used primer pairs targeting several variable regions of the 16S rRNA gene (27) to determine the number of different oralbacterial and oral-archaeal species with in silico amplicon similarity values ≥97% (ASI97), as well as the potential OTUs that might contain distinct species. Moreover, for the first time in this kind of analysis, we described the specific taxa of the oral ecosystem with highly similar sequence segments, specifying if they belong to different genera or other higher taxonomic ranks. 3.5.1. Evaluation of the primer pairs for detecting oral species with in silico amplicon similarity values ≥97% In the present study, the primer pairs that targeted bacteria had a mean of 91.88 (49.40%) bacterial species with an ASI97 and an average of 153.46 potential OTUs containing distinct species. For those targeting archaea, these numbers were 65.60 (48.59%) and 162.26, respectively. Using the percentage species coverage with no in silico amplicons similarity ≥97% (SC-NASI97) as a selection criterion, the optimum primer pair for detecting oral bacteria was OP_F053-KP_R020. Although the primer used most in the oral microbiome studies, KP_F031-KP_R021 identified slightly fewer species with an ASI97 (37 vs. 58) and number of ASI97 (32 vs. 46); its SC-NASI97 was also lower than that of OP_F053-KP_R020 (54.30% vs. 65.05%). The primer pair producing the best estimates for detecting oral archaea was KP_F018- KP_R002. Again, the widely used primer KP_F014-KP_R011, although it only detected a few ALBA REGUEIRA IGLESIAS 296 species with an ASI97 (24 vs. 60) and number of ASI97 (96 vs. 125), however, also had a considerably lower SC-NASI97 than that of KP_F018-KP_R002 (12.59% vs. 51.11%). Lastly, we recommend the primer OP_F114-KP_R031 for detecting oral bacteria and archaea simultaneously. OP_F114-KP_R002, meanwhile, identified slightly fewer taxa with an ASI≥97% (for bacteria= 84 and for archaea= 60 vs. 85 and 62) and number of ASI97 (118 and 126 vs. 133 and 136) but had a lower SC-NASI≥97% (47.31% and 54.81% vs. 51.08% and 53.33%). In addition, as previously observed in objectives 1 and 2 (27,44), none of the primer combinations that are most commonly employed in sequencing-based studies of the oral microbiome were among the best. Specifically, the species coverage of KP_F078-OP_R010, a primer described by Caporaso (45), fell from 94.62% for bacteria and 63.70% for archaea as described in objective 2 (44) to 34.95% and 31.11%, respectively; when considering the species with an ASI97, possibly generating as many as 215 and 99 potential bacterial and archaeal OTUs, respectively; that contain different species. 3.5.2. Description of the distinct pairs of oral-bacteria species and oral-archaea species with in silico amplicon similarity values ≥97% Around 80% of the oral-bacteria and oral-archaea species analysed had an ASI97 with at least another species. The widely-known bacterial periodontopathogens F. nucleatum and T. denticola (46-48) had similar in silico amplicons to Fusobacterium hwasookii and Treponema putidum, respectively, which have also been detected in periodontal lesions (49,50). Interestingly, other bacteria with high in silico amplicon similarities had antagonistic roles in oral health and disease. Examples are: the health-associated C. concisus and the initially periodontitis-associated C. curvus (51); the health-related Rothia mucilaginosa (52) and the decay-abundant R. dentocariosa (53,54); the commensal S. mitis, oralis, and salivarius; the caries-associated S. mutans (46,55,56); and the periodontal health-related Tannerella sp. oral taxon HOT-286 (57,58) and the periodontitis-related T. forsythia (46-48). Furthermore, relevant oral-disease associated species, such as A. actinomycetemcomitans (46,59) and R. dentocariosa (53,54), were among those that had an ASI97 with taxa from distinct genera. Regarding the archaea, we found that four Methanosarcina species found in healthy and periodontitis pockets, namely barkeri, lacustris, mazeii, and vacuolata (60), were highly similar. Moreover, H. ruber, Methanotorris igneus, M. zhilinae, and N. occultus, which are reported to be among the 10 most Objective 3 297 abundant species in both healthy and periodontitis subjects (60), had an ASI97 with several taxa from distinct genera. Schloss (13) has recently stated that the risks of artificially splitting a genome into multiple amplicon sequence variants (ASVs) are greater than those of clustering ASVs from different species into the same OTU when using broad distance thresholds. However, considering the results obtained in the present study, our opinion is that the latter approach should be avoided in the analysis of the oral microbiota if the aim is to associate species with specific clinical conditions. In silico amplicons from species traditionally associated with contrary health conditions, like those described above, can be grouped with a ≥97% similarity threshold. This would result in both an overabundance of the single species representing the OTU and an underestimation of the diversity of the community, with other species within the OTU overlooked. Consequently, it would be better to use the lowest possible level of resolution, i.e., the variant level (23), and databases specifically designed for taxonomic identifications of taxa at this level (43). It has been demonstrated that distinct OTU clustering approaches, or even the same method, can yield uneven results for the same dataset (9-11). Therefore, we decided to analyse the 97% similarity relationships between oral species, without considering the influence of any clustering algorithm. Consequently, the results presented here are an approximation of the different oral species that could be grouped in potential OTUs. 3.5.3. Limitations of the present study The main limitation of our study is that we have only considered one, randomly selected, of all possible in silico amplicons with ASI97 between two different species to establish the existence of a close relationship between the two. Another consideration is that we were only able to evaluate 25% of the oral microorganism genomes listed on the eHOMD website, as the remainder were not fully sequenced. This absence of complete genomes reduced the number of species investigated to 35% of those set out on the site. Although the analysis could have been performed on annotations of the 16S rRNA gene sequences from oral microbes, we preferred to use complete genomes, thereby ensuring the high quality of the sequences reviewed. The reasons why we adopted this approach were: 1) Edgar (61) estimated that the taxonomy ALBA REGUEIRA IGLESIAS 298 annotation error rate of the ribosomal database project (RDP) database (62) is ∼10%; on the other hand, he found 249,490 identical sequences with conflicting annotations in SILVA v128 (63) and Greengenes v13.5 (64) at ranks up to phylum (7,804 conflicts), indicating that the annotation error rate in these databases is ∼17%; 2) we have verified in objective 1 that a very high percentage of 16S rRNA gene annotations present a loss of information of up to 60 - 70 nucleotides in regions 1 and 9 of the sequences, which invalidates their use (27); 3) most of the complete genomes evaluated here are isolates that were sequenced with Sanger technology or with second-generation technology (shorter sequences than Sanger). In both cases, contig scaffolding algorithms were used to construct the complete genomes from the sequences with a minimum coverage of 8x for Sanger sequences and 30x in the case of second-generation technologies (65). In these types of assemblies, positions within the genome that did not have high coverage included non-specific nucleotides. In the present study, we discarded genomes that included more than 20 consecutive unspecific positions; 4) in addition, many genomes were downloaded from the NCBI RefSef database (31), where the annotations of the complete genomes were manually curated or re-annotated concerning the information provided by the original author, including their taxonomic hierarchy. Thus, our results highlight only part of a much more extensive problem. Objective 3 299 3.6. CONCLUSIONS In conclusion, the tested primer pairs targeting bacteria and/or archaea detected an average of more than 150 potential OTUs that might contain different species, when ≥97% similarity threshold was used. According to the SC-NASI97 parameter, the best primer pairs were: OP_F053-KP_R020 for bacteria (region 1-3; primer pair position for Escherichia coli J01859.1: 9-356); KP_F018-KP_R002 for archaea (4 undefined-532); and OP_F114-KP_R031 for both (3-5; 340-801). Around 80% of the oral-bacteria and oral-archaea species analysed had an ASI97 with at least one other species. These very similar species play different roles in the oral microbiota and belong to bacterial genera such as Campylobacter, Rothia, Streptococcus, and Tannerella, and archaeal genera such as Halovivax, Methanosalsum,and Methanosarcina. Moreover, ~20% and ~30% of these two-by-two similarity relationships were established between species from different bacterial and archaeal genera, respectively. Even taxa from distinct families, orders, and classes could be grouped in the same potential OTU. Consequently, regardless of the primer pair used, sequence-clustering with ≥97% similarity provides an inaccurate description of oral-bacterial and oral-archaeal species, which can greatly affect microbial diversity parameters. As a result, OTU clustering conditions the credibility of associations between some oral species and certain health and disease conditions. This significantly limits the comparability of the microbial diversity findings reported in oral microbiome literature. ALBA REGUEIRA IGLESIAS 300 3.7. REFERENCES (1) Midha MK, Wu M, Chiu KP. Long-read sequencing in deciphering human genetics to a greater depth. Hum Genet. 2019 Dec;138(11-12):1201-15. (2) Davidson RM, Epperson LE. Microbiome sequencing methods for studying human diseases. Methods Mol Biol. 2018;1706:77-90. (3) Zaura E, Pappalardo VY, Buijs MJ, Volgenant CMC, Brandt BW. Optimizing the quality of clinical studies on oral microbiome: a practical guide for planning, performing, and reporting. Periodontol 2000. 2021 Feb;85(1):210-36. (4) Edgar R. UPARSE: highly accurate OTU sequences from microbial amplicon reads. Nat Methods. 2013 Oct;10(10):996-8. (5) Stackebrandt E, Goebel, BM. Taxonomic note: a place for DNA-DNA reassociation and 16s rRNA sequence analysis in the present species definition in bacteriology. Int J Syst Evol Microbiol. 1994 Oct;44(4):846-9. (6) Bolyen E, Rideout JR, Dillon MR, Bokulich NA, Abnet CC, Al-Ghalith G, et al. Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2. Nat Biotechnol. 2019 Aug;37(8):852-7. (7) Schloss PD, Westcott SL, Ryabin T, Hall JR, Hartmann M, Hollister EB, et al. Introducing mothur: open-source, platform-independent, community-supported software for describing and comparing microbial communities. Appl Environ Microbiol. 2009 Dec;75(23):7537-41. (8) Edgar R. Search and clustering orders of magnitude faster than BLAST. Bioinformatics. 2010 Oct;26(19):2460-1. (9) Wei Z, Zhang X, Cao M, Liu F, Qian Y, Zhang S. Comparison of methods for picking the operational taxonomic units from amplicon sequences. Front Microbiol. 2021 Mar; 12:644012. doi: 10.3389/fmicb.2021.644012. Objective 3 301 (10) He Y, Caporaso JG, Jiang XT, Sheng HF, Huse SM, Rideout JR, et al. Stability of operational taxonomic units: an important but neglected property for analyzing microbial diversity. Microbiome. 2015 May;3:20. doi: 10.1186/s40168-015-0081-x. (11) Westcott SL, Schloss PD. De novo clustering methods outperform reference-based methods for assigning 16S rRNA gene sequences to operational taxonomic units. PeerJ. 2015 Dec;3:e1487. doi: 10.7717/peerj.1487. (12) Větrovský T, Baldrian P. The variability of the 16S rRNA gene in bacterial genomes and its consequences for bacterial community analyses. PLoS One. 2013;8(2):e57923. doi: 10.1371/journal.pone.0057923. (13) Schloss PD. Amplicon sequence variants artificially split bacterial genomes into separate clusters. mSphere. 2021 Aug;6(4):e0019121. doi: 10.1128/mSphere.00191-21. (14) Edgar R. UNOISE2: improved error-correction for Illumina 16S and ITS amplicon sequencing. Preprint at bioRxiv. 2016. doi: 10.1101/081257. (15) Eren AM, Morrison HG, Lescault PJ, Reveillaud J, Vineis JH, Sogin ML. Minimum entropy decomposition: unsupervised oligotyping for sensitive partitioning of high-throughput marker gene sequences. ISME J. 2015 Mar;9(4):968-79. (16) Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJ, Holmes SP. DADA2: highresolution sample inference from Illumina amplicon data. Nat Methods. 2016 Jul;13(7):581-3. (17) Amir A, McDonald D, Navas-Molina JA, Kopylova E, Morton JT, Zech Xu Z, et al. Deblur rapidly resolves single-nucleotide community sequence patterns. mSystems. 2017 Mar;2(2):e00191-16. doi: 10.1128/mSystems.00191-16. (18) Caruso V, Song X, Asquith M, Karstens L. Performance of microbiome sequence inference methods in environments with varying biomass. mSystems. 2019 Feb;4(1):e00163-18. doi: 10.1128/mSystems.00163-18. ALBA REGUEIRA IGLESIAS 302 (19) Nearing JT, Douglas GM, Comeau AM, Langille MGI. Denoising the denoisers: an independent evaluation of microbiome sequence error-correction approaches. PeerJ. 2018 Aug;6:e5364. doi: 10.7717/peerj.5364. (20) Prodan A, Tremaroli V, Brolin H, Zwinderman AH, Nieuwdorp M, Levin E. Comparing bioinformatic pipelines for microbial 16S rRNA amplicon sequencing. PLoS One. 2020 Jan;15(1):e0227434. doi: 10.1371/journal.pone.0227434. (21) Abellan-Schneyder I, Matchado MS, Reitmeier S, Sommer A, Sewald Z, Jan Baumbach J, et al. Primer, pipelines, parameters: issues in 16S rRNA gene sequencing. mSphere. 2021 Feb;6(1):e01202-20. doi: 10.1128/mSphere.01202-20. (22) García-López R, Cornejo-Granados F, Lopez-Zavala A, Cota-Huízar A, Sotelo-Mundo R, Gómez-Gil B, et al. OTUs and ASVs produce comparable taxonomic and diversity from shrimp microbiota 16S profiles using tailored abundance filters. Genes (Basel). 2021 Apr;12(4):564. doi: 10.3390/genes12040564. (23) Callahan BJ, McMurdie PJ, Holmes SP. Exact sequence variants should replace operational taxonomic units in marker-gene data analysis. ISME J. 2017 Dec;11(12):2639-43. (24) F Escapa I, Chen T, Huang Y, Gajare P, Dewhirst FE, Lemon KP. New insights into human nostril microbiome from the expanded Human Oral Microbiome Database (eHOMD): a resource for the microbiome of the human aerodigestive tract. mSystems. 2018 Dec;3(6):e00187-18. doi: 10.1128/mSystems.00187-18. (25) Clark K, Karsch-Mizrachi I, Lipman DJ, Ostell J, Sayers EW. GenBank. Nucleic Acids Res. 2016 Jan;44:D67-72. doi: 10.1093/nar/gkv1276. (26) NCBI Resource Coordinators. Database resources of the National Center for Biotechnology Information. Nucleic Acids Res. 2016 Jan;44(D1):D7-19. doi: 10.1093/nar/gkv1290. Objective 3 303 (27) Regueira-Iglesias A, Vázquez-González L, Balsa-Castro C, Vila-Blanco N, Blanco-Pintos T, Tamames J, et al. In silico evaluation and selection of the best 16S rRNA gene primers for use in next-generation sequencing to detect oral bacteria and archaea. Accepted for publication in Microbiome. Preprint at Research Square. 2021. doi: 10.21203/rs.3.rs-516961/v1. (28) National Center for Biotechnology Information. Entrez programming utilities help. 2010; Available at: https://www.ncbi.nlm.nih.gov/books/NBK25501/. (29) Python Software Foundation. Python. Version 3.9.0. 2020; Available at: http://www.python.org/. (30) Schoch CL, Ciufo S, Domrachev M, Hotton CL, Kannan S, Khovanskaya R, et al. NCBI Taxonomy: a comprehensive update on curation, resources and tools. Database (Oxford). 2020 Jan;2020:baaa062. doi: 10.1093/database/baaa062. (31) O'Leary NA, Wright MW, Brister JR, Ciufo S, Haddad D, McVeigh R, et al. Reference sequence (RefSeq) database at NCBI: current status, taxonomic expansion, and functional annotation. Nucleic Acids Res. 2016 Jan;44(D1): D733-45. doi: 10.1093/nar/gkv1189. (32) Baker GC, Smith JJ, Cowan DA. Review and re-analysis of domain-specific 16S primers. J Microbiol Methods. 2003 Dec;55(3):541-55. (33) Barnett M. regex. 2020; Available at: https://pypi.org/. (34) Cock PJ, Antao T, Chang JT, Chapman BA, Cox CJ, Dalke A, et al. Biopython: freely available Python tools for computational molecular biology and bioinformatics. Bioinformatics. 2009 Jun;25(11):1422-23. (35) Camacho C, Coulouris G, Avagyan V, Ma N, Papadopoulos J, Bealer K, et al. BLAST+: architecture and applications. BMC Bioinformatics. 2009 Dec;10:421. doi: 10.1186/1471- 2105-10-421. ALBA REGUEIRA IGLESIAS 304 (36) Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990 Oct;215(3):403-10. (37) Chen Y, Ye W, Zhang Y, Xu Y. High speed BLASTN: an accelerated MegaBLAST search tool. Nucleic Acids Res. 2015 Sep;43(16):7762-8. (38) McKinney W. Data structures for statistical computing in Python. In: van der Walt S, Millman J, editors. Proceedings of the 9th Python in Science Conference; 2010; Austin. Texas: SciPy; 2010. doi: 10.25080/Majora-92bf1922-00a. (39) McNamara J. xlsxwriter. 2013; Available at: https://xlsxwriter.readthedocs.io/. (40) R Core Team. R: a language and environment for statistical computing. R package version 4.0.3. Vienna, Austria: R Foundation for Statistical Computing; 2020; Available at: https://www.R-project.org/. (41) Csardi G, Nepusz T. The Igraph software package for complex network research. InterJournal, Complex Systems. 2006; 1695; Available at: http://igraph.org. (42) Westcott SL, Schloss PD. OptiClust, an improved method for assigning amplicon-based sequence data to operational taxonomic units. mSphere. 2017 Mar;2(2):e00073-17. doi: 10.1128/mSphereDirect.00073-17. (43) F Escapa I, Huang Y, Chen T, Lin M, Kokaras A, Dewhirst FE, et al. Construction of habitat-specific training sets to achieve species-level assignment in 16S rRNA gene datasets. Microbiome. 2020 May;8(1):65. doi: 10.1186/s40168-020-00841-w. (44) Regueira-Iglesias A, Vázquez-González L, Balsa-Castro C, Blanco-Pintos T, Vila-Blanco N, Carreira MJ, et al. Impact of 16S rRNA gene redundancy and primer pair selection on the quantification and classification of oral microbiota in next-generation sequencing. Preprint at Research Square. 2021. doi: 10.21203/rs.3.rs-662236/v1. 311 Objective 4. A large-scale meta-omics analysis of plaque microbiota in periodontal diseases 4.1. ABSTRACT Aims: To analyse the supragingival and subgingival plaque microbiota at ASV level of different periodontal conditions (periodontal health, gingivitis, and untreated and treated periodontitis) in terms of bacterial diversity, co-occurrence networks, and predictive models. Material and methods: A total of 120 patients (55 controls, 65 periodontitis) were selected for subgingival plaque collection. Sequencing of the 3-4 16S rRNA gene region was performed in Illumina MiSeq. The obtained sequences and metadata were uploaded to the sequence read archive (SRA). Searches were performed in PubMed, Scopus, Embase, and the SRA to identify previously published Illumina 3-4 sequencing studies on the supragingival and subgingival plaque microbiome in distinct periodontal conditions. Research that met the criteria for sequences and metadata were included in the meta-omics analysis, comprising a total of 2045 samples. Sequences were processed under the same bioinformatics protocol, which included the ASV-level classification and the use of an oral-specific database for taxonomic classification. The statistical analysis was conducted using the phyloseq, DESeq2, microbiome, mixOmics, vegan, SpiecEasi, and igraph packages. Results and conclusions: Bacterial richness associated with periodontitis was higher than in health in supragingival plaque and lower in subgingival, but evenness was higher in disease in both niches. The supragingival microbiota was richer and more diverse than the subgingival for the same periodontal condition. The structure of the bacterial community differed among conditions in the supra- and subgingival plaque, as well as for the same health status between the two niches. In addition, the core microbiota of dental plaque did not allow the characterisation of periodontal health and disease; and the proportion of the bacterial community organised in co-occurrence networks at the ASV level was very small. However, a small proportion of supra- and subgingival taxa had outstanding ability to distinguish between ALBA REGUEIRA IGLESIAS 312 periodontal conditions, and a relevant percentage of them were core members. Supragingival plaque was a better bacterial biomarker than subgingival for discriminating periodontal health from untreated and treated periodontitis. The main health-predictor ASVs in supragingival and subgingival plaque were: R. dentocariosa ASV2, H. parainfluenzae ASV3, ASV78, ASV45, and ASV46, K. oralis ASV66, S. vestibularis ASV27, and A. HMT170 ASV119. The main predictor ASVs of periodontitis in dental plaque were: T. forsythia ASV15, F. alocis ASV19, T. denticola ASV38 and ASV150, F. fastidiosum ASV97, P. HMT369 ASV124, S. anginosus ASV142, and P. nodatum ASV189. 4.1.1. Keywords Meta-omics analysis; next-generation sequencing; 16S rRNA gene; dental plaque; supragingival; subgingival; microbiota; periodontal diseases. 4.1.2. Declaration of conflict of interest The doctoral candidate and the rest of the authors of the present study declare that they have no conflict of interest concerning the objectives proposed in this chapter. 4.1.3. Funding This investigation was supported by the Instituto de Salud Carlos III (General Division of Evaluation and Research Promotion, Madrid, Spain) and co-financed by the FEDER (European Regional Development Fund, ERDF) (“A way of making Europe”) under grant ISCIII/PI21/00588; the Consellería de Cultura, Educación e Ordenación Universitaria de la Xunta de Galicia (group with growth potential ED431B 2020-2022 GPC2020/27; A. Regueira- Iglesias support ED481A-2017/233) and the ERDF, which acknowledges the CiTIUS-Research Center in Intelligent Technologies of the Santiago de Compostela University as a Research Center of the Galician University System. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Objective 4 313 4.2. INTRODUCTION Hundreds of articles have been published in the last two decades on the use of nextgeneration sequencing (NGS) of the 16S ribosomal RNA (rRNA) gene as a way to study the oral microbiome. These have generally analysed different intra-oral niches like dental plaque (1), tongue coatings (2), the soft tissues (3), and saliva (4) to determine the microbial diversity associated with distinct conditions, including: periodontal health (5); periodontal and periimplant diseases (6); dental caries (7); and oral cancer (8). The large amount of such scientific production and the variation in the results obtained have caused researchers to conduct numerous narrative reviews in an attempt to achieve a consensus when defining the microbial profiles for distinct periodontal health statuses (9-12). However, these published studies on the periodontal microbiome vary in terms of the relevant steps undertaken within a typical 16S rRNA gene sequencing workflow. This has had significant effects on the diversity of the results obtained, making comparisons very difficult (13-15). It is well known that each sequencing technology performs differently in the trade-off between read length, sequence throughput, and error rate (13); being Illumina that with preferable performance over Roche 454 or Ion Torrent (15). On the other hand, we have recently demonstrated through the in silico analysis performed in objective 1 (16) that, even among primer pairs with coverage values ≥90%, the oral species detected by primers targeting a particular region tended to be not covered by others amplifying a different zone and vice versa. Consequently, it can be said that it is rather questionable to compare sequences and consequently microbial diversity data derived from distinct sequencing technologies and gene regions. On the other hand, more than 80% of recently published studies of the periodontal microbiome used the clustering of operational taxonomic units (OTUs) to perform their analyses. However, the 97% similarity threshold that is typically employed means that community descriptions based on this approach are wildly inaccurate, since 80% of oralbacterial and archaeal species have an amplicon sequence similarity ≥97% to at least one other oral species as described in objective 3 (17). It is therefore necessary to conduct periodontal microbiota analyses using techniques that are currently considered to be more reliable, for example by examining levels of amplicon sequence variants (ASVs) (18-20). Furthermore, ALBA REGUEIRA IGLESIAS 314 high-quality, oral-specific databases are required if accurate classifications of these ASVs are to be achieved (21). In an attempt to produce the strongest evidence to date on the periodontal microbiota, we conducted the large-scale meta-omics research described herein. This had the following objective: 1) to analyse the supragingival and subgingival plaque microbiota at ASV level of different periodontal conditions (healthy periodontal, gingivitis, periodontitis, and treated periodontitis) in terms of bacterial diversity, co-occurrence networks, and predictive models. To achieve our objective, we re-analysed sequences stored in public repositories from previously published Illumina 3-4 sequencing studies on the periodontal microbiome in supragingival and subgingival plaque. Our sample also included a bioproject with in-house sequences of the same region, which were taken from the subgingival plaque of periodontally healthy and periodontitis patients from our setting. The meta-omics analysis employed a unique bioinformatics protocol for high-quality filtering and sequence analysis. Objective 4 315 4.3. MATERIAL AND METHODS The complete analysis protocol applied in the present study is detailed in figure 1. Figure 1. The complete analysis protocol applied in the present meta-omics study. 4.3.1. Selection of study groups and collecting the subgingival plaque samples A convenience sample of 120 eligible participants, comprising 55 periodontally healthy controls (control group) and 65 subjects affected by untreated periodontitis (periodontitis group), were recruited from 350 consecutive patients in the general population who were referred to the School of Medicine and Dentistry (Universidade de Santiago de Compostela, Spain) and the Instituto Superior de Ciências da Saúde Norte, Cooperativa de Ensino Superior, Politécnico e Universitário (CESPU, Gandra, Paredes, Portugal) between 2018 and 2019 for an assessment of their oral health status. Patients were recruited if they fulfilled the following inclusion criteria: 1) age 24 to 75; 2) the presence of at least 15 natural teeth; 3) no previous periodontal treatment; 4) no medical ALBA REGUEIRA IGLESIAS 316 history of diabetes mellitus, hepatic or renal disease, or other severe medical conditions or transmittable diseases; 5) no intake of systemic antimicrobials during the previous six months; 6) no intake of anti-inflammatory medication in the previous four months; 7) no routine use of oral antiseptics; 8) no history of alcohol or drug abuse; 9) no pregnancy or breastfeeding; 10) no presence of implants or orthodontic appliances; 11) have smoked for at least one year; and 12) have never smoked or stopped more than three years ago. Two experienced dentists performed all the periodontal diagnoses. The bleeding on probing (BOP) and the bacterial plaque level (BPL) were recorded for the full mouth on a binary scale (presence/absence) at six sites per tooth. We also documented the probing pocket depth (PPD) and clinical attachment level (CAL) throughout the mouth, again at six sites per tooth, using a PCP-UNC 15 probe. Standardised radiographs of all the teeth were obtained to assess the alveolar bone status. The diagnosis of periodontitis was based on the clinical and radiographic information obtained. The control group included periodontally healthy patients who had: BOP ≤20%, no location with a PPD ≥4 mm, and no radiographic evidence of alveolar bone loss. The presence of periodontal health or moderate to severe generalised chronic periodontitis was established according to the clinical/radiographic information, applying previously published criteria (22,23). The "smoking habit" of the participants was evaluated using a questionnaire, with information collected on its extent, i.e., non-smoker, former smoker, current smoker, time spent as a former or current smoker, and the number of cigarettes consumed per day. The research was conducted following the principles of the Declaration of Helsinki (revised in 2000) on studies involving human experimentation (24), and its protocol was approved by the Galician Clinical Research Ethics Committee (registration number 2018/295) and the Instituto Superior de Ciências da Saúde-Norte, CESPU (registration number 35/CEIUCS/2019) (Appendix S1). All the participants provided their written informed consent to their involvement in the study. The plaque collection took place one or two weeks after the initial examination. Subgingival plaque samples from the controls and periodontal patients were collected and Objective 4 317 pooled from eight non-adjacent proximal sites using two paper strips inserted into the gingival sulcus or periodontal pocket for 30 seconds. In the first case, samples were taken from subgingival healthy sites in quadrants one and three, and in the second case from sites with the most in-depth PPD in each quadrant. The strips used with each recruit were inserted into labelled tubes with 300 ml of 0.01M phosphate-buffered saline (PBS) (pH=7.2) and frozen at - 80ºC until further genomic analysis. 4.3.2. 16S rRNA gene amplicon sequencing of subgingival samples Total DNA was extracted from the subgingival plaque samples using a commercial kit (MasterPure Complete DNA and RNA Purification Kit; Epicentre,Wisconsin, USA) according to the manufacturer’s instructions, albeit with minor modifications, including a mechanical disruption of bacteria (Pathogen Lysis Tube S; Qiagen, Hilden, Germany), and the addition of a lysozyme treatment (20 mg/ml at 37 ºC for 30 minutes). The isolated DNA was eluted in 50 µl of distilled and apyrogenic water, and its quality and concentration were assessed using a Nanodrop spectrophotometer (ND-2000 Spectrophotometer, Wilmington, USA). DNA samples with spectrophotometer ratios (Abs 260/280) between 1.5 and 2.0 were considered to be acceptable for inclusion in the study. Two subgingival samples from the control group were excluded due to non-compliance with this requirement. A polymerase chain reaction (PCR) amplification of the 16S rRNA gene was performed with the KAPA HiFi HotStart ReadyMixPCR Kit (Cat. No. KK2602, 7958935001; Kapa Biosystems, F. Hoffmann-La Roche Ltd, Basel, Switzerland). The 3-4 hypervariable region was amplified as previously described (25) using the following primers in a limited-cycle PCR: 3-4-Forward (5′-CCT ACG GGNGGC WGC AG-3). 3-4-Reverse (5′-GAC TAC HVGGG TAT CTA ATC C-3). A set of modified primers, 3-4-F and 3-4-R, were also used. This set contained a 1-3 base pair (bp) "heterogeneity spacer" that we designed to mitigate the issues caused by low-sequence diversity amplicons. ALBA REGUEIRA IGLESIAS 318 Each PCR amplification was carried out on a total volume of 10 μl, which comprised 4 μl of DNA, 0.2-Μm from each forward and reverse primer, and a Kapa ready mix (Kapa Biosystems). The PCR conditions were modified by conducting: 1) an initial denaturation at 95°C for 3 minutes; 2) 25 three-step cycles at 95°C for 30 seconds, 55°C for 30 seconds and 72°C for 30 seconds; and 3) a final 5-minute extension at 72ºC. Water, up to a total volume of 50 μl, was added after the first PCR step. The reactions were purified using AMPure XP beads (Beckman Coulter, Brea, CA, USA) with a 0.9X (3-4 amplicon) ratio, according to the manufacturer’s instructions. The PCR products were eluted from the magnetic beads with 32 μl of Buffer EB (Qiagen N.V, Hilden; Germany), with 30 μl of the eluate transferred to a fresh 96-well plate. The primers described above contain overhangs that enable the addition of full-length Nextera adapters. Barcodes are available for multiplex sequencing in a second PCR step, which produces sequencing-ready libraries. To this end, 5 μl of the first amplification was used as a template for the second PCR, with Nextera XT v2 adaptor primers added up to a final volume of 50 μl. The PCR mix and thermal profile employed for the first PCR were also used for the second, but only for eight cycles. After the second PCR, 25 μl of the final product was purified and normalised with the SequalPrep normalisation kit (Invitrogen, Carlsbad, CA, USA), according to the manufacturer’s protocol. Libraries were eluted in a 20 μl volume and pooled for sequencing. Final pools were quantified with a quantitative PCR (qPCR) using the Kapa library quantification kit for Illumina Platforms (Kapa Biosystems) on an ABI 7900HT real-time cycler (Applied Biosystems, Foster City, CA, USA). Sequencing using v3 chemistry with a loading concentration of 18 pM was performed in Illumina MiSeq (Illumina Inc., San Diego, CA, USA) with 2x300 bps reads. In all cases, 10% of the PhIX control libraries were spiked to increase the diversity of the sequenced samples. In parallel, negative control tests of the sample-collection buffer, DNA-extraction and PCR-amplification steps were conducted routinely under the same conditions and using reagents. One such non-template control was subjected to the library preparation and then Objective 4 319 sequenced. As expected, this yielded very few reads (1611 per sample). This was in contrast to an average of 249,747 reads/library in the sample-derived collections. The bacterial mock community as a positive control for the downstream procedures were taken from the ZymoBIOMICS Microbial Community DNA Standard (Catalog Number D6306, Zymo Research, Irvine, CA, USA), which is a mix of genomic DNA isolated from pure cultures of eight bacterial and two fungal strains. Mock DNAs were amplified and sequenced in the same way as all the other samples used in the experiment. The sequences obtained were deposited in the sequence read archive (SRA) database (26) under accession number PRJNA773202. 4.3.3. Characteristics of the studies for the meta-omics analysis: inclusion and exclusion criteria Studies (cross-sectional, longitudinal, or interventional) on the microbial diversity in both supragingival and subgingival plaque in adult individuals with different periodontal conditions were included in our research (periodontal health, gingivitis, periodontitis, treated periodontitis, periimplantitis, and treated periimplantitis). We incorporated all the studies in which the diversity of the periodontal microbiome was assessed using primers from the 3-4 region and the Illumina-sequencing technology. An associated bioproject number indicates the repository in which the sequences are stored. Studies were included in our analysis if the reference standard for diagnosing a periodontal condition was based on only clinical (PPD or CAL) or clinical and radiographic parameters (bone loss -BL-), irrespective of the diagnostic benchmarks applied. Consequently, in the absence of homogeneous criteria, any definition based on the author’s reported standards was accepted. Studies without a reference for diagnosing the periodontal condition were ineligible for inclusion, as were those that failed to assess the periodontal status of patients using at least one clinical parameter (either the PPD or CAL). ALBA REGUEIRA IGLESIAS 320 4.3.4. Characteristics of the metadata table and the stored sample sequences: inclusion and exclusion criteria After applying the criteria described above to the studies in the literature, we further selected those where the metadata of interest per sample was properly assigned in the repository. In relation to the characteristics of the stored sequences, the inclusion and exclusion criteria were as follows: 1) direct and reverse sequences were accepted, whether or not the primer-pair sequence was included; 2) contigs with or without primer pairs, whose minimum average length had to be ≥350 bps; 3) the primer sequences must have been aligned with the complete Escherichia coli J01859.1 16S rRNA gene sequence using BLAST (27) to determine their initial and end positions; if they corresponded to the region of interest, the bioproject was included in the analysis; 4) study samples without primers were accepted for the analysis if multiple sequences were selected and aligned with the full 16S rRNA gene of E. coli J01859.1 to confirm that they belonged to the region of interest; 5) bioprojects in which most of the samples had a very low number of stored sequences (≤7000 sequences) were rejected; and 6) bioprojects were excluded if the samples were multiplexed or had different barcodes in each file. 4.3.5. Search methods for the identification and selection of studies and bioprojects 4.3.5.1. Information sources and search strategy The searches were conducted in July 2021 using the electronic databases PubMed, Scopus, and Embase. The search strategy to identify Illumina sequencing-based studies of the periodontal microbiome encompassed two sets of terms relating to: 1) periodontal health conditions, oral niches and microbiota; and 2) the 16S rRNA gene sequencing technology (Appendix S2). All the searches in the three databases were filtered according to the publication year - 2000 to 2021 (inclusive). The searches of Scopus and Embase were also filtered by the type of: document/publication (Embase: article, article in the press; or review); source (journal); and language (English). Additional searches of the SRA database (26) were performed using the terms “periodontitis”, “periodontal health”, “periodontal disease”, “peri-implantitis”, “gingivitis”, Objective 4 327 The groups Sub_x0GDx, Imp_x0IHx, and Imp_x0IDx were removed due to their low sample sizes (n= <50), leaving a total of 2045 samples to be analysed. An independent filter had previously excluded from the statistical analysis the ASVs with an abundance of ≤10 counts and a presence in ≤2 samples (42), leaving a final total of 8379 ASVs. The relationship between the different periodontal health conditions and the plaque microbiota was investigated from several perspectives: 1) the alpha diversity indicators and the structure of the bacterial community; 2) the composition of the core microbiota and the testing of differential abundance; 3) the co‑occurrence network patterns; and 4) the predictive capacity of the plaque microbiota for discriminating the periodontal health condition. In general, where applicable, the comparative analyses were first performed between different clinical conditions within the same niche (supragingival plaque or subgingival plaque), and then in the same clinical condition between the two different plaques. 4.3.12.1. Alpha diversity indicators and the structure of the bacterial community The phyloseq and microbiome packages were used to obtain the alpha diversity data (39,41). As indicators of taxa richness, we calculated the absolute count data ("observed") and the coverage index, which defines how many of the more abundant ASVs are required to achieve a particular proportion of the occupied ecosystem (95%). The Shannon and Pielou indices were determined as indicators of diversity and the evenness of the ASVs present in the samples (43,44). The Mann-Whitney U test (two-tailed) was used to conduct different comparative analyses. A principal component analysis (PCA) was employed to visualise the clustering of the plaque samples in relation to their respective periodontal health condition. The mixOmics package (version 6.16.3) (45) was used to obtain the scatter plots of the first two principal components based on the relative abundance of the ASVs, showing the centroids of each clinical group and the ellipses representing the 95% confidence interval. A non-parametric permutational multivariate analysis of variance (PERMANOVA) (46) was used to measure the ALBA REGUEIRA IGLESIAS 328 multivariate community-level differences between the groups. These analyses were performed using the vegan package (version 2.5-7) (47). 4.3.12.2. Composition of the core plaque microbiota and testing the differential abundance The microbiome package (41) was used to identify the core ASVs present at a prevalence rate of ≥75% in each plaque type and each periodontal clinical condition. The DESeq2 package (40) was used to identify the ASVs with the most significant changes in differential abundance for the different periodontal conditions. Improvements to the stability and dispersion of the counts (variance) were required before it was possible to calculate the differential abundances. To this end, we used the estimate SizeFactors function in DESeq2 (40) to transform the stabilisation of the variance. The differential abundances were measured with the log2foldchange (log2FC) value and the different conditions were compared using the Wald test with the Benjamini–Hochberg correction (Q parameter= 0.1, false discovery rate (FDR) <10%). The differential-abundance measurements were statistically significant if the adjusted p-value was < 0.01 (−log10 adjusted p-value= 2). 4.3.12.3. Co‑occurrence networks in the plaque microbiota Co-occurrence network analyses were performed with the clinical groups with more than 100 samples, filtering out ASVs with an abundance of 0.01%. The SparCC method was used to generate the networks (48), as this allows researchers to detect with a high degree of accuracy the linear relationships in both a set of samples and a compositional dataset (49). The default parameters and the SpiecEasi package (version 1.1.1) (50) were used to run SparCC, and the correlation matrix obtained was filtered using an absolute correlation score greater than or equal to 0.5. The networks were then visualised with the igraph package (version 1.2.6) (51), where each node represents an ASV and each edge represents the correlations between the abundances of the ASVs. A set of measures was calculated to describe the topology of the resulting networks: 1) the network coverage, defined as the percentage of ASVs present in the co-occurrence network Objective 4 329 concerning the total number of ASVs detected in the corresponding group; 2) the number of nodes and edges; 3) the number of sub-networks; and 4) the number of modules (52). We calculated the betweenness centrality (BC) (53) to measure the relative importance of each ASV within the network (how influential a taxon is within a network). This determines the fraction of the shortest paths through one particular bacterial taxon to another. The BC of a taxon in a network reflects the importance of the control exerted by the taxon over the interactions of other taxa in the same network (53). In line with Banerjee et al. (54), a combined score based on a high degree value and a high BC value was used as a threshold to define the hub or keystone ASVs in the microbial communities. 4.3.12.4. Predictive capacity of the plaque microbiota for discriminating the clinical condition. We conducted a supervised classification in the form of a sparse partial least-squares discriminant analysis (sPLS-DA) (55) to facilitate the categorisation of the different clinical groups and identify the ASVs that best distinguished two groups within each plaque niche (supragingival and subgingival plaques); consequently, the discriminant models were calculated two-to-two. The sPLS-DA was performed using the mixOmics package (45), which is dedicated to the integrative examination of “omics” data. The ASVs with a relative abundance of less than 0.1% in the total samples were previously excluded from the development of predictive models. The number of components in each model was determined by applying the rule of thumb K-1, where K is the number of classes (in our case, two clinical groups). Consequently, all the predictive models were of one component. Receiver operating characteristic (ROC) curves were constructed with the true positivity rate (sensitivity) as a function of the false positivity rate (1-specificity), while area under the curve (AUC) values were used to distinguish between each clinical group in the supragingival and subgingival plaque. It should be noted that simulations with an AUC value equal to or higher than 0.70 are generally considered to be acceptable predictive models (56). ALBA REGUEIRA IGLESIAS 330 4.4. RESULTS 4.4.1. Clinical characteristics of the study groups from our setting Subjects in our own-recruited periodontitis group had a higher mean age and a higher number of smokers than those in the healthy group. Smokers with periodontitis consumed more cigarettes per day and had been smoking for more months. Regarding the clinical parameters associated with the periodontal status, patients in the periodontitis group had significantly higher BPL in the full-mouth, and BOP, PPD, and CAL values both in the full mouth and sampled sites records than subjects in the healthy group (p<0.001; Table 1). Table 1. Age, sex, smoking habit, and clinical characteristics associated with periodontal status in our ownrecruited healthy and periodontitis groups. Clinical parameters Study groups Control ( n= 53 ) * Periodontitis ( n= 65 ) p value Age (years) 46.13 (12.11) 52.31 (9.78) 0.002 Sex Female 27 38 NS Male 26 27 Smoking habit Non-smokers 44 30 <0.001 Smokers 9 35 Cigarettes/day (no.) 1.02 (2.78) 8.62 (9.89) <0.001 Months of smoking (no.) 39.55 (105.66) 175.75 (185.12) <0.001 No. of teeth 26.98 (2.38) 25.22 (4.01) 0.012 Full mouth BOP (%) 10.91 (6.27) 50.29 (20.43) <0.001 BPL (%) 22.62 (17.48) 55.02 (27.21) <0.001 PPD (mm) 2.03 (0.26) 3.61 (0.72) <0.001 CAL (mm) 2.20 (0.40) 4.40 (1.15) <0.001 Sampled sites BOP (%) 5.92 (7.23) 66.38 (24.88) <0.001 PPD (mm) 2.20 (0.25) 5.58 (0.76) <0.001 CAL (mm) 2.29 (0.32) 6.11 (1.05) <0.001 *Of the 55 initial control subjects, two were excluded due to non-compliance with the requirements for the amount of DNA extracted. Values indicate means (standard deviations) and the number of subjects. After applying the Shapiro-Wilks test and verifying the non-normal distribution of almost all the clinical variables, the Mann-Whitney U test (two-tailed) was used to compare the quantitative clinical variables between the control and periodontitis groups. The Fisher’s exact test (two-tailed) was used to assess the association of the qualitative variables between the two study groups. A significance level of p<0.05 was established. BOP= bleeding on probing; BPL= bacterial plaque level; CAL= clinical attachment level; mm= milimetres; n= sample size; No.= number; NS= No significant; PPD= probing pocket depth. 4.4.2. Studies and bioprojects obtained in the search process Figure 2 shows the flowchart of the search process, including the number of results obtained from each step. Objective 4 331 Figure 2. Flowchart of the search process. *The exclusion reasons of the rejected articles and bioprojects are indicated in appendices S4, S5, and S6. The abstracts of 30331 articles derived from the searches of the electronic databases were analysed computationally, using seven sets of positive terms to select the candidates for evaluation. A total of 1159 articles from these databases and 39 bioprojects from the SRA were ALBA REGUEIRA IGLESIAS 332 evaluated. These are listed in appendices S4 and S5, respectively; and the exclusion reason is indicated if applicable. Ultimately, 32 articles in which the sequence data had been deposited in 32 different bioprojects met the inclusion criteria (Appendix S6). The bioproject containing our sequences was added at this point. Ten authors were contacted to obtain or clarify data and three provided the information required. Five articles were excluded in the metadata assessment step, and a further four after the samples and sequences were evaluated. We were ultimately left with 23 articles (6,57-78) and 25 bioprojects for inclusion in the meta-analysis, involving a total of 2045 samples distributed in eight periodontal clinical groups (four groups in supragingival plaque and the remaining four in subgingival plaque). 4.4.3. Quality assessment of metadata and sequences 4.4.3.1. Quality of the metadata stored in repositories Three of the 25 included bioprojects had high-quality metadata (range= 0.95 - 0.77), four were medium quality (range= 0.45 - 0.34), and 18 low quality (range= 0.30 - 0.15). Overall, the those with medium and low-quality metadata did not include information about the periodontitis type and severity, ethnicity, or clinical periodontal parameters (total and sampling site). Moreover, most of the bioprojects with low-quality metadata did not provide information on the age or sex of the participants. 4.4.3.2. Sample size and number of sequences stored in repositories From a sample size point of view, 10 bioprojects had <50 samples (40%), 9 between 50 and 100 (36%), and six >100 (24%). No bioprojects had an ASS of <0.25, as these had been eliminated in previous steps because of their very low-quantity sequences. Four bioprojects (BP20, BP21, BP22, and BP25) involving a total of 167 samples, representing 8.17% of all of those analysed, had ASS values from 0.75 - 1.0. These were therefore of an acceptable quantity, with more than 7500 sequences per sample. Eight bioprojects (BP12, BP14, BP18, BP19, BP23, BP27, BP44, and BP45) and 422 samples, representing 20.63% of all of those processed, had ASS values from 1.0 - 2.0. These were thus deemed to be high quantity, with 10,000 - 20,000 sequences per sample. Finally, 13 bioprojects, with an overall total of 1456 samples, Objective 4 333 representing 71.20% of all of those processed, had an ASS >2.0, making them very high quantity, with more than 20,000 sequences per sample. Figure 3. The methodological quality of selected studies and bioprojects: (A) metadata; and (B) sample size and sequence quantity. Twenty-three articles and 24 bioprojects were included since one article referred to two distinct bioprojects (+ 1 bioproject associated with our samples). 4.4.4. Characteristics of the selected studies and bioprojects Appendix S7 contains a quantitative summary of the main descriptive characteristics of the sequencing-based studies of the periodontal microbiome that formed part of our meta-omics research. Less than a third (7/23 articles + 1 own unpublished bioproject; 29.17%) were able to establish the periodontal diagnosis with the new Classification of Periodontal and Peri-implant ALBA REGUEIRA IGLESIAS 334 Diseases and Conditions (79), with most using earlier classifications or the authors’ own criteria (17/24; 70.83%). In 13/24 investigations (54.17%), there was a comparison of the microbial profiles in relation to states of periodontal health and disease, while 10/24 (41.66%) only evaluated periodontitis. There was also one article (4.17%) where only five healthy samples were selected, as these were the only ones that could be assigned to a specific health condition. Subgingival plaque was used the most to study the periodontal microbiota (16/24; 66.66%), followed by supragingival plaque (4/24; 16.67%) or both types (4/24; 16.67%). Moreover, 4/24 studies (16.67%) assessed the changes produced in the microbiota after non-surgical periodontal therapy, including (in some cases) the adjuvant effect of antibiotics or toothpastes. 4.4.5. Alpha-diversity in supragingival and subgingival plaque microbiota 4.4.5.1. Supragingival plaque microbiota As shown in table 2, supragingival plaque richness decreased significantly from periodontal health to gingivitis and then increased strongly in the periodontitis condition (median number of ASVs observed= 610.50, 474.00, and 892.00, respectively; 95% coverage index= 220.00, 130.00, and 288.00, respectively). There was a significant decrease in both the number of ASVs and the 95% coverage index in the post-periodontal therapy samples compared to those collected before treatment (Sup_x1PDx vs. Sup_x0PDx: 781.00 and 263.00 vs. 892.00 and 288.00). However, these post-treatment estimates of richness did not reach the levels of the healthy group, with significant differences remaining between the two clinical conditions (Sup_x1PDx vs. Sup_x0HHx, 781.00 and 263.00 vs. 610.50 and 220.00). Conversely, the diversity and evenness indexes showed a continuous upwards trend from health to disease and even continued to improve after treatment (Shannon index range between 4.75 and 4.07; Pielou index range between 0.70 and 0.62). All the two-by-two group comparisons were significantly different, excepting Sup_x0HHx and Sup_x0GDx (Shannon index) and Sup_x0GDx and Sup_x0PDx (Pielou index). Objective 4 335 Table 2. Alpha diversity indicators in the different periodontal health conditions and dental plaque types. Groups (n) No. Observed ASVs Coverage index ( 95% ) Shannon index Pielou index Median values and IQR Sup_x0HHx (210) 610.50 (458.00) 220.00 (145.00) 4.07 (0.75) 0.62 (0.13) Sup_x0GDx (79) 474.00 (162.50) 130.00 (89.50) 4.19 (0.65) 0.68 (0.11) Sup_x0PDx (493) 892.00 (912.00) 288.00 (165.00) 4.51 (0.65) 0.67 (0.12) Sup_x1PDx (81) 781.00 (333.00) 263.00 (114.00) 4.75 (0.57) 0.70 (0.08) Sub_x0HHx (155) 478.00 (1142.50) 171.00 (169.00) 4.15 (1.18) 0.65 (0.12) Sub_x0PHx (62) 474.00 (320.50) 120.00 (93.75) 4.05 (0.94) 0.65 (0.13) Sub_x0PDx (768) 417.50 (455.25) 142.00 (130.00) 4.17 (0.98) 0.68 (0.10) Sub_x1PDx (197) 507.00 (461.00) 129.00 (105.00) 4.20 (0.83) 0.69 (0.10) Comparison of distinct periodontal health conditions in the supragingival plaque (p-value) Sup_x0HHx vs. Sup_x0GDx 1.0227E-06 8.0438E-12 NS 0.0002 Sup_x0HHx vs. Sup_x0PDx 7.0505E-14 7.7671E-13 8.5206E-24 2.7619E-08 Sup_x0HHx vs. Sup_x1PDx 0.0010 0.0008 9.9932E-20 2.9953E-16 Sup_x0GDx vs. Sup_x0PDx 9.6812E-21 1.2208E-29 4.3632E-08 NS Sup_x0GDx vs. Sup_x1PDx 5.4703E-15 1.7601E-17 4.5964E-11 0.0003 Sup_x0PDx vs. Sup_x1PDx 0.0030 0.0466 0.0003 1.2601E-07 Comparison of distinct periodontal health conditions in the subgingival plaque (p-value) Sub_x0HHx vs. Sub_x0PHx NS 0.0129 NS NS Sub_x0HHx vs. Sub_x0PDx 0.0005 0.0158 NS 0.0001 Sub_x0HHx vs. Sub_x1PDx 0.0139 0.0384 NS 2.4036E-05 Sub_x0PHx vs. Sub_x0PDx NS NS NS 0.0017 Sub_x0PHx vs. Sub_x1PDx NS NS 0.0475 0.0003 Sub_x0PDx vs. Sub_x1PDx NS NS NS NS Comparison of the same periodontal health condition between the supragingival and subgingival plaques ( p-value ) Sup_x0HHx vs. Sub_x0HHx 0.0363 9.2241E-06 NS 0.0050 Sup_x0PDx vs. Sub_x0PDx 8.0165E-72 5.6815E-97 2.2485E-23 8.4811E-05 Sup_x1PDx vs. Sub_x1PDx 6.2408E-13 3.8383E-20 1.8787E-12 0.0194 A significance level of p<0.05 was established. ASVs= amplicon sequence variants; IQR= interquartile range; n= sample size; No.= number; NS= No significant; Sub_x0HHx= subgingival plaque of periodontally healthy subjects, healthy sites; Sub_x0PHx= subgingival plaque of periodontitis subjects, healthy sites; Sub_x0PDx= subgingival plaque of periodontitis subjects, diseased sites; Sub_x1PDx= subgingival plaque of periodontitis subjects, diseased sites after therapy; Sup_x0GDx= supragingival plaque of gingivitis subjects, diseased sites; Sup_x0HHx= supragingival plaque of periodontally healthy subjects, healthy sites; Sup_x0PDx= supragingival plaque of periodontitis subjects, diseased sites; Sup_x1PDx= supragingival plaque of periodontitis subjects, diseased sites after therapy. 4.4.5.2. Subgingival plaque microbiota Regarding the subgingival microbiota, significantly lower values were detected in the number of ASVs and the 95% coverage index in diseased sites from periodontitis patients with respect to healthy patients (417.50 and 142.00 vs. 478.00 and 171.00, respectively). After periodontal treatment, there was a significant increase in bacterial richness, surpassing even the healthy levels (507.00 compared to 478.00), although the 95% coverage index remained lower than the healthy levels (129.00 compared to 171.00). The diversity and evenness indices tended to increase in the groups of diseased tooth sites, although only the Pielou index comparisons ALBA REGUEIRA IGLESIAS 336 were significant (0.65 in Sub_x0HHx vs. 0.68 and 0.69 in Sub_x0PDx and Sub_x1PDx, respectively). The number of observed ASVs and the diversity and evenness indices of the subgingival microbiota did not vary significantly between Sub_x0HH and Sub_x0PHx. Similarly, the alpha diversity indicators did not show significant variations between the different periodontal groups, except for the Pielou index in the comparisons Sub_x0PHx vs. Sub_x0PDx (0.65 vs. 0.68) and Sub_x0PHx vs. Sub_x1PDx (0.65 vs. 0.69), and the Shannon in Sub_x0PHx vs. Sub_x1PDx (4.05 vs. 4.20) (Table 2). 4.4.5.3. Supragingival and subgingival plaque When contrasting the “supra” and “sub” plaques of the subjects with the same periodontal health status, we observed that the number of ASVs, the 95% coverage index values, and the Shannon diversity scores were significantly higher in the supragingival niche than in the subgingival niche (“supra” vs. “sub”: ASV number range= 892.00 - 610.50 vs. 507.00 - 417.50; 95% coverage index range= 288.00 - 220.00 vs. 171.00 - 129.00; Shannon diversity range= 4.75 - 4.07 vs. 4.20 - 4.15); the exception was represented by the Shannon index of the healthy groups for both plaques. Conversely, the evenness values were significantly higher in the subgingival environment, except for the case of Sup_x1PDx (Table 2). 4.4.6. Structure of the bacterial community in supragingival and subgingival plaque microbiota The PCAs revealed a grouping of the supragingival and subgingival samples according to the periodontal health condition of the subject and the sampled site (the latter in the case of the Sub_x0PHx group) (Figures 4 and 5). The visual observations were confirmed by the PERMANOVA, which produced significant results for all the two-by-two group comparisons (Table 3). In the comparison between the different niches of the same periodontal health condition, the PCA revealed a clustering of the samples according to the type of plaque collected from the subject for the same periodontal health status (Figure 6). The visual observations were Objective 4 343 abundances (16.96% and 45.42% of the total detected by the two groups, respectively), while Sup_x0HHx vs. Sup_x0GDx, 945 ASVs and 198 species (15.09% and 39.52%, respectively). The comparison of Sup_x0HHx and Sup_x1PDx revealed 926 ASVs and 210 species both differentially abundant (13.12% and 41.02% of the detected taxa, respectively), and, again, Sup_x0HHx vs. Sup_x0PDx, a total of 918 (12.17%) and 272 (51.52%) ASVs and species, respectively. In contrast, the lowest relative numbers of ASVs and species with differential abundances were observed in the analysis of Sup_x0PDx vs. Sup_x1PDx (total= 660 ASVs, 8.95%; 145 species, 27.62%). The percentages of core ASVs and core species showing differential abundance ranged from 6.31% - 2.87% and 14.65% - 7.35%, respectively (Table 5). 4.4.8.2. Subgingival plaque microbiota The results for the subgingival plaque demonstrated that the highest relative numbers of ASVs and species with differential abundances were obtained when comparing the periodontal health group to both the non-treated and treated periodontitis groups. Accordingly, the comparison of Sub_x0HHx vs. Sub_x0PDx revealed a total of 1074 ASVs (12.64% of the total detected by the two groups) from 273 species with differential abundances (48.75%), while Sub_x0HHx vs. Sub_x1PDx had 1015 ASVs (14.45%) from 225 species (41.67%). Conversely, the lowest relative numbers of ASVs and species with differential abundances were observed in the analysis of Sub_x0PHx vs. both Sub_x0PDx (total= 364 ASVs; 4.28%; 156 species; 27.81%) and Sub_x1PDx (total= 339 ASVs; 5.64%; 117 species; 22.20%). The percentages of core ASVs and core species showing differential abundance ranged from 6.59% - 2.06% and 10.26% - 5.13%, respectively (Table 5). ALBA REGUEIRA IGLESIAS 344 Table 5. Number of total and core taxa that presented differential abundances in the different periodontal health conditions and dental plaque types, and the relative abundance values they represented. No. ASVs ( % detected ) No. Species ( % detected ) No. Core ASVs ( % detected ) * No. Core species ( % detected ) * Differential abundances of distinct periodontal health conditions in the supragingival plaque Sup_x0HHx vs. Sup_x0GDx 945 (15.09%) 198 (39.52%) 41 (4.34%) 29 (14.65%) Sup_x0HHx vs. Sup_x0PDx 918 (12.17%) 272 (51.52%) 33 (3.59%) 20 (7.35%) Sup_x0HHx vs. Sup_x1PDx 926 (13.12%) 210 (41.02%) 33 (3.56%) 21 (10.00%) Sup_x0GDx vs. Sup_x0PDx 1290 (16.96%) 243 (45.42%) 37 (2.87%) 25 (10.29%) Sup_x0GDx vs. Sup_x1PDx 507 (10.08%) 163 (33.61%) 32 (6.31%) 16 (9.82%) Sup_x0PDx vs. Sup_x1PDx 660 (8.95%) 145 (27.62%) 19 (2.88%) 14 (9.66%) Differential abundances of distinct periodontal health conditions in the subgingival plaque Sub_x0HHx vs. Sub_x0PHx 425 (6.62%) 160 (29.52%) 13 (3.06%) 12 (7.50%) Sub_x0HHx vs. Sub_x0PDx 1074 (12.64%) 273 (48.75%) 36 (3.35%) 25 (9.16%) Sub_x0HHx vs. Sub_x1PDx 1015 (14.45%) 225 (41.67%) 27 (2.66%) 18 (8.00%) Sub_x0PHx vs. Sub_x0PDx 364 (4.28%) 156 (27.81%) 24 (6.59%) 16 (10.26%) Sub_x0PHx vs. Sub_x1PDx 339 (5.64%) 117 (22.20%) 7 (2.06%) 6 (5.13%) Sub_x0PDx vs. Sub_x1PDx 604 (7.14%) 189 (33.87%) 14 (2.32%) 12 (6.35%) Differential abundances of the same periodontal health condition between supragingival and subgingival plaques Sup_x0HHx vs. Sub_x0HHx 802 (10.57%) 255 (46.88%) 36 (4.49%) 25 (9.80%) Sup_x0PDx vs. Sub_x0PDx 2367 (27.34%) 349 (62.21%) 48 (2.03%) 33 (9.46%) Sup_x1PDx vs. Sub_x1PDx 198 (3.85%) 72 (14.55%) 14 (7.07%) 4 (5.56%) The percentages of detected ASVs and species are calculated concerning the total number of different ASVs and species detected by at least one of the groups to be compared. *The percentages of core ASVs and species are calculated concerning the total number of different ASVs and species that showed differential abundances in the two groups compared. The taxa that could not be classified at the species level (“unclassified”) were counted once so the number of species detected is the minimum that could be obtained. ASVs= amplicon sequence variants; No.= number; Sub_x0HHx= subgingival plaque of periodontally healthy subjects, healthy sites; Sub_x0PHx= subgingival plaque of periodontitis subjects, healthy sites; Sub_x0PDx= subgingival plaque of periodontitis subjects, diseased sites; Sub_x1PDx= subgingival plaque of periodontitis subjects, diseased sites after therapy; Sup_x0GDx= supragingival plaque of gingivitis subjects, diseased sites; Sup_x0HHx= supragingival plaque of periodontally healthy subjects, healthy sites; Sup_x0PDx= supragingival plaque of periodontitis subjects, diseased sites; Sup_x1PDx= supragingival plaque of periodontitis subjects, diseased sites after therapy. 4.4.8.3. Supragingival and subgingival plaque microbiota As shown in table 5, the comparison of Sup_x0PDx vs. Sub_x0PDx revealed the highest relative numbers of ASVs and species with differential abundances (total= 2367 ASVs, 27.34% of the total detected by the two groups; 349 species, 62.21%). Conversely, the lowest relative estimates were observed in the analysis of Sup_x1PDx vs. Sub_x1PDx (total= 198 ASVs, 3.85%; 72 species, 14.55%). The percentages of core ASVs and core species showing differential abundance ranged from 7.07% - 2.03% and 9.80% - 5.56%, respectively (Table 5). Objective 4 345 4.4.9. Co-occurrence networks in supragingival and subgingival plaque microbiota Table 6 shows the topological parameters of the co-occurrence networks in the two supragingival and subgingival plaque groups that met the inclusion criteria for this analysis. Table 6. Topological parameters of the co-occurrence networks in the different periodontal health conditions and dental plaque types. Supragingival plaque Subgingival plaque Sup_x0HHx Sup_x0PDx Sub_x0HHx Sub_x0PDx Sub_x1PDx Network coverage* 2.26% 2.54% 2.75% 0.63% 1.54% Number of nodes 136 187 163 53 78 Number of edges 290 959 387 80 111 Number of positive correlations (%) 290 (100.0%) 958 (99.9%) 387 (100.0%) 80 (100.0%) 111 (100.0%) Number of negative correlations ( % ) 0 (0.0%) 1 (0.1%) 0 (0.0%) 0 (0.0%) 0 (0.0%) Ratio of positive correlations and nodes 2.13% 5.12% 2.37% 1.51% 1.42% Number of subnetworks 18 12 12 10 12 Number of modules 25 55 25 11 13 Number of modules with more than 3 nodes 12 12 11 4 6 *Percentage of ASVs present in the co-occurrence network with respect to the total number of ASVs detected in the correspondent group. Sub_x0HHx= subgingival plaque of periodontally healthy subjects, healthy sites; Sub_x0PDx= subgingival plaque of periodontitis subjects, diseased sites; Sub_x1PDx= subgingival plaque of periodontitis subjects, diseased sites after therapy; Sup_x0HHx= supragingival plaque of periodontally healthy subjects, healthy sites; Sup_x0PDx= supragingival plaque of periodontitis subjects, diseased sites. 4.4.9.1. Supragingival plaque microbiota The network coverage and the number of nodes in Sup_x0PDx were slightly higher than in Sup_x0HHx (2.54% and 187 vs. 2.26% and 136, respectively). Moreover, the number of edges was more than three times greater in the diseased than in the healthy group (959 vs. 290, respectively). Practically all these correlations were positive in both groups, except for that of Dialister invisus ASV68 and Streptococcus unclassified ASV4 in Sup_x0PDx (correlation value= -0.51). The diseased network had fewer subnetworks and a higher number of modules (12 and 55 vs. 18 and 25 in Sup_x0HHx), but both groups had the same number of modules with more than three nodes (Table 6). ALBA REGUEIRA IGLESIAS 346 In the Sup_x0HHx network, the three main hubs or keystone ASVs were: Streptococcus unclassified ASV90, Rothia dentocariosa ASV2, and Streptococcus oralis subsp. dentisani clade 058 ASV1. All were part of the Sup_x0HHx core microbiota but, despite having an abundance ≥0.5%, none of the three taxa were differentially abundant when compared to Sup_x0PDx. The principal keystone ASVs in the Sup_x0PDx network were: Streptococcus unclassified ASV85, Streptococcus sanguinis ASV228, and Streptococcus unclassified ASV121. Only the former had an abundance ≥0.5% in Sup_x0PDx, but none were core members or had differential abundance if compared to Sup_x0HHx. Figures 7 and 8 represent the main modules of the co-occurrence networks associated with Sup_x0HHx and Sup_x0PDx, respectively. Objective 4 347 Sup_x0HHx Samples= 210 ASVid Genus Species ASV Core Relative abundance AV00001 Streptococcus oralis_subsp.dentisani _clade_058 BTASV016027 Y 11.2800 AV00085 Streptococcus Unclassified unclassified Y 0.6454 AV00090 Streptococcus Unclassified unclassified Y 0.5700 AV00121 Streptococcus Unclassified unclassified Y 0.3356 AV00155 Streptococcus Unclassified unclassified Y 0.2059 AV00181 Streptococcus Unclassified unclassified N 0.1094 AV00220 Streptococcus Unclassified unclassified N 0.1252 AV00228 Streptococcus Sanguinis unclassified Y 0.2314 AV00282 Streptococcus Unclassified unclassified Y 0.1702 AV00320 Streptococcus oralis_subsp.dentisani _clade_058 unclassified Y 0.1281 AV00355 Streptococcus Unclassified unclassified N 0.0609 AV00392 Streptococcus Sanguinis unclassified Y 0.1019 AV00425 Streptococcus oralis_subsp.dentisani _clade_058 unclassified N 0.0682 AV00546 Streptococcus Unclassified unclassified N 0.0861 Figure 7. Main module of the co-occurrence network associated with the supragingival plaque of periodontally healthy subjects. In the graph, the most important taxa are highlighted in green, orange, and yellow according to the score obtained in the analysis group.The highest value is shown in green, the values belonging to the first quartile in orange, and those belonging to the second quartile, i.e. up to the median, in yellow. ASV= amplicon sequence variant; ASVid= amplicon sequence variant identifier; N= no; subsp.= subspecies; Sup_x0HHx= supragingival plaque of periodontally healthy subjects, healthy sites; Y= yes. ALBA REGUEIRA IGLESIAS 348 Sup_x0PDx Samples= 493 ASVid Genus Species ASV Core Relative abundance AV00004 Streptococcus unclassified unclassified Y 4.0960 AV00085 Streptococcus unclassified unclassified N 0.5505 AV00090 Streptococcus unclassified unclassified N 0.5082 AV00155 Streptococcus unclassified unclassified N 0.1948 AV00228 Streptococcus sanguinis unclassified N 0.1995 AV00390 Streptococcus unclassified unclassified N 0.0686 AV00392 Streptococcus sanguinis unclassified N 0.1129 AV00479 Streptococcus unclassified unclassified N 0.0901 AV00535 Streptococcus unclassified unclassified N 0.0837 AV00546 Streptococcus unclassified unclassified N 0.0796 AV00613 Streptococcus unclassified unclassified N 0.0663 AV00619 Streptococcus sanguinis unclassified N 0.0534 AV00624 Streptococcus unclassified unclassified N 0.0492 AV00741 Streptococcus unclassified unclassified N 0.0558 AV00990 Streptococcus unclassified unclassified N 0.0325 Figure 8. Main module of the co-occurrence network associated with the supragingival plaque of periodontitis subjects in the diseased sites. In the graph, the most important taxa are highlighted in green, orange, and yellow according to the score obtained in the analysis group.The highest value is shown in green, the values belonging to the first quartile in orange, and those belonging to the second quartile, i.e. up to the median, in yellow. ASV= amplicon sequence variant; ASVid= amplicon sequence variant identifier; N= no; subsp.= subspecies; Sup_x0PDx= supragingival plaque of periodontitis subjects, diseased sites; Y= yes. [Document text truncated for crawler view.]