Especiação mediada por características e extinções induzidas por humanos em proboscídeos reveladas por redes neurais bayesianas não supervisionadas
Abstrato
Traços da história de vida das espécies, paleoambiente e interações bióticas provavelmente influenciam as taxas de especiação e extinção, afetando a riqueza de espécies ao longo do tempo. Modelos de nascimento-morte que inferem o impacto desses fatores geralmente assumem relações monotônicas entre preditores e taxas individuais, limitando nossa capacidade de avaliar efeitos mais complexos e sua importância e interação relativas. Apresentamos um modelo bayesiano de nascimento-morte usando redes neurais não supervisionadas para explorar efeitos multifatoriais e não lineares nas taxas de especiação e extinção usando dados fósseis. Ele infere taxas específicas de linhagem e tempo e desembaraça os efeitos e a importância dos preditores por meio de técnicas de inteligência artificial explicáveis. A análise do registro fóssil proboscídeo revelou taxas de especiação moldadas pela flexibilidade alimentar e eventos biogeográficos. O surgimento dos humanos modernos aumentou as taxas de extinção, causando declínio recente da diversidade, enquanto o clima regional teve um impacto menor. Nosso modelo abre caminho para uma melhor compreensão da dinâmica intrincada que molda a diversificação de clados.
INTRODUÇÃO
A especiação e a extinção remodelam constantemente a biosfera em escalas temporais evolutivas. O equilíbrio mutável entre esses dois processos dita os padrões de biodiversidade e sua variação entre clados e ao longo do tempo. Isso inclui eventos de rotatividade que revendem o funcionamento do ecossistema regional ( 1 , 2 ) , o surgimento de padrões biogeográficos globais, como o gradiente de diversidade latitudinal ( 3–5 ), a ascensão e queda de linhagens inteiras sob o domínio da perturbação ambiental ou competição crescente ( 6–8 ) e grandes transições biológicas provocadas por extinções em massa e subsequentes recuperações de diversidade que virtualmente reinicializam a vida na Terra ( 9 ). Crucialmente, os processos de especiação e extinção também desempenham um papel de liderança na formação da disparidade biológica. Quando a poda e a brotação da árvore da vida são seletivas e favorecem certas propriedades das espécies em detrimento de outras [isto é, classificação ou seleção de espécies; ( 10 )], a especiação e a extinção tornam-se construtores eficazes de tendências evolutivas ativas que estimulam linhagens em regiões inexploradas do morfoespaço [( 11 – 13 ) ver “nota adicionada na prova”]. Assim, não é de surpreender que a investigação evolutiva tenha procurado durante muito tempo ligar padrões de diversidade a processos evolutivos ( 14 – 16 ). Ao fazê-lo, um objetivo fundamental da investigação macroevolutiva tem sido quantificar os efeitos das propriedades das espécies na diversificação ( 17 – 19 ), os efeitos da diversificação seletiva nas tendências evolutivas ( 20 ) e o papel da força ambiental nesta interação ( 21 , 22 ).
Como esses processos evolutivos deixam sua pegada no registro fóssil e nas árvores filogenéticas moleculares ( 23 ), estudar ambas as fontes de informação se tornou um foco central de pesquisa para biólogos evolucionistas. Avaliações quantitativas de padrões de diversificação ao longo do tempo e entre clados têm uma longa tradição em paleobiologia ( 24 , 25 ) e levaram ao desenvolvimento de uma ampla gama de métodos e modelos para inferir taxas de especiação e extinção ( 26 – 30 ). Da mesma forma, reconhecer que as filogenias moleculares encapsulam até certo ponto o sinal de eventos de diversificação passados ( 31 – 33 ) alimentou o desenvolvimento de modelos de diversificação adaptados para analisar árvores filogenéticas datadas de táxons existentes ( 34 – 37 ).
O reino da inferência de especiação e extinção testemunhou o florescimento de vários modelos que levam em conta a dependência de características e séries temporais no processo de diversificação [por exemplo, ( 8 , 18 , 28 , 38 – 44 )]. No entanto, a maioria dos modelos de diversificação permite apenas a estimativa dos efeitos de preditores únicos (por exemplo, características ou mudanças paleoambientais) nas taxas de especiação e extinção. Esta é uma limitação crítica, pois é improvável que apenas um único fator tenha uma alavancagem constante na especiação ou extinção em escalas de tempo macroevolutivas. Além disso, há muito tempo se reconhece que a relação causal entre características fenotípicas e diversificação deve variar com as mudanças nas configurações ambientais ( 21 ). Por exemplo, espécies com tamanho corporal maior podem ter maior probabilidade de desaparecer durante eventos severos de extinção ( 45 ). No entanto, à medida que a produtividade e o funcionamento do ecossistema se recuperam, nichos relacionados a tamanhos corporais maiores podem ter maior potencial de diversificação quando guildas de menor tamanho se tornam saturadas ( 13 , 46 ). Embora alguns dos modelos permitam uma análise conjunta de múltiplos fatores combinados (por exemplo, múltiplas características ou variáveis de tempo), eles normalmente assumem que seus efeitos são aditivos, independentes e limitados por funções monotônicas (frequentemente lineares) ( 47 – 51 ). Assim, quaisquer potenciais efeitos não lineares e a interação mutável entre atributos de espécies e mudanças paleoambientais na regulação das taxas de especiação e extinção têm se mostrado difíceis de quantificar.
Aqui, propomos um modelo flexível para analisar dados de ocorrência de fósseis que permite variação de taxa dependente de tempo e característica, integrando assim aspectos de modelos exploratórios e baseados em hipóteses em uma única estrutura coesa. Nosso modelo é probabilístico e baseado em um processo estocástico de nascimento-morte no qual as taxas são moduladas ao longo do tempo e entre linhagens, usando uma função flexível de tempo, variáveis paleoambientais e múltiplas características categóricas e/ou contínuas, juntamente com sua interação. Nosso modelo não assume uma função particular ligando características e tempo a taxas variáveis de especiação e extinção a priori. Em vez disso, o modelo estima essa função a partir dos dados por meio de uma rede neural não supervisionada. Aproveitando metodologias de inteligência artificial explicável (xAI), podemos estimar os efeitos de fatores individuais e classificar sua importância na formação de especiação e extinção. Aplicamos nosso modelo para analisar o rico registro fóssil de proboscídeos, examinando os papéis do paleoclima, ecomorfologia, biogeografia e sobreposição espaço-temporal com os humanos na formação dos padrões de diversificação dessa ordem de mamíferos.
RESULTADOS
Um novo modelo de diversificação dependente de características e tempo
Desenvolvemos um modelo bayesiano para inferir a dinâmica de diversificação de ocorrências fósseis com base em um processo de nascimento-morte acoplado a um processo de preservação e amostragem. Nossa abordagem se baseia no modelo PyRate ( 30 ) que usa um algoritmo bayesiano para amostrar conjuntamente os tempos de origem e extinção de todos os táxons em um conjunto de dados, bem como as taxas de especiação, extinção e preservação. Ao implementar processos de nascimento-morte variáveis no tempo, o PyRate pode estimar variações nas taxas de especiação e extinção ao longo do tempo. Nosso novo modelo estende essa abordagem ao permitir que as taxas de especiação e extinção variem tanto ao longo do tempo quanto em um nível específico de linhagem. Essa variação é inferida como uma função do próprio tempo, da relação filogenética (por exemplo, com base na taxonomia ou em uma árvore filogenética) e virtualmente qualquer número de características fenotípicas contínuas e categóricas e variáveis ambientais dependentes do tempo. Em vez de assumir uma função de resposta predefinida ligando esses fatores à especiação e extinção, usamos uma rede neural não supervisionada para modelar sua correlação, permitindo assim respostas não lineares e interações entre os fatores. Usamos um algoritmo de Monte Carlo de cadeia de Markov (MCMC) para amostrar os parâmetros da rede neural, a partir dos quais as taxas de especiação e extinção específicas de linhagem e tempo são derivadas. Essas taxas então alimentam a função de probabilidade de um processo de nascimento-morte que, juntamente com os priores, determina a probabilidade de aceitação no algoritmo MCMC. Como em outras implementações do PyRate, o processo de nascimento-morte é acoplado a um modelo de preservação, permitindo a variação das taxas de amostragem ao longo do tempo e entre as linhagens [ver ( 52 )].
Por meio da implementação de diferentes algoritmos adaptados do xAI, avaliamos o impacto dos fatores na modelagem da dinâmica de diversificação. Especificamente, visualizamos a magnitude e a forma do efeito inferido de características e variáveis dependentes do tempo nas taxas usando gráficos de dependência parcial (PD) ( 53 ), que mostram o efeito de um preditor enquanto marginaliza os outros. Classificamos os preditores por importância usando três abordagens xAI com base em permutações de preditores ( 54 ), probabilidades marginais e valores de SHapley Additive exPlanation (SHAP) ( 55 ). Essas métricas capturam diferentes aspectos do impacto do preditor: alavancagem na probabilidade de nascimento-morte ( Eqs. 8 e 9 ), consistência da direção de seu efeito e tamanho do efeito.
Desempenho do modelo de rede neural de nascimento-morte
Usamos simulações para avaliar a capacidade de nossa estrutura de inferir com precisão taxas específicas de linhagem-tempo e identificar corretamente os fatores que causam variação de taxa. Especificamente, geramos 100 conjuntos de dados de ocorrência de fósseis para nove cenários de diversificação. Esses cenários acomodam taxas de especiação e extinção que são constantes ao longo do tempo e iguais entre linhagens (cenário 1), bem como taxas que dependem do próprio tempo (por exemplo, por meio de mudanças de taxa), características categóricas ou contínuas específicas de linhagem ou uma série temporal de paleotemperatura e interações entre esses fatores (ver fig. S1). Cenários com taxas dependentes do tempo são interpretados como casos em que um fator que varia com o tempo impulsiona a mudança de taxa. Por exemplo, casos com mudanças de taxa (cenários 2, 3 e 8) podem refletir mudanças na dinâmica de especiação e extinção determinadas por uma grande transição ambiental. As simulações que retratam taxas crescentes de extinção que finalmente atingem a especiação (cenário 5) capturam a dinâmica antecipada sob um modelo de competição dependente da diversidade ( 6 , 36 ).
Todos os conjuntos de dados incluíam múltiplas características, variáveis dependentes do tempo e parentesco filogenético, mesmo quando não afetavam as taxas. Em simulações com taxas mudando em tempos predefinidos de mudança (ou seja, cenários 2, 3 e 8), os tempos de mudanças de taxa foram considerados desconhecidos e não fornecidos como parte dos conjuntos de dados. Da mesma forma, no cenário 9, as taxas eram dependentes de uma característica, que, no entanto, omitimos do conjunto de dados para avaliar o desempenho do modelo de rede neural de nascimento-morte (BDNN) quando faltavam os verdadeiros preditores. Nesses casos, esperamos que o sinal de variação de taxa seja capturado pelo próprio tempo ou pelo parentesco filogenético [aqui quantificado como autovetores filogenéticos ( 56 ); veja Materiais e Métodos] como proxies para os tempos de mudança e características faltantes.
As taxas de especiação e extinção específicas do tempo de linhagem inferidas foram precisas em todos os cenários, com os erros relativos absolutos medianos sendo mais baixos em simulações de taxa constante e variando entre 0,23 e 0,32 na maioria dos outros cenários ( Fig. 1 e tabela S1). No entanto, em simulações com taxas constantes por partes (cenário 3), o erro aumentou para ∼0,38 a 0,66, pois o modelo BDNN tendeu a suavizar as mudanças de taxa instantâneas simuladas (fig. S4D). Cenários com taxas dependentes de características categóricas ou contínuas (simulações 5 e 6) produziram erros baixos em torno de 0,23 a 0,25. Simulações mais complexas envolvendo os efeitos interativos de características e variáveis temporais (cenários 4, 7 e 8) ainda resultaram em precisão ligeiramente menor com erros em torno de 0,26 a 0,48, e precisão semelhante foi obtida mesmo quando os verdadeiros preditores foram omitidos (cenário 9). Com a implementação bayesiana do modelo BDNN, fomos capazes de calcular intervalos de credibilidade (ICs) de 95% ou taxas específicas de linhagem-tempo para avaliar com que frequência eles incluem os valores verdadeiros usados para simular sua origem e extinção (ou seja, cobertura). O BDNN produziu uma cobertura variando de 0,85 a 0,99 na maioria das simulações, com valores mais baixos no cenário 3 (0,63 a 0,79; tabela S1). A cobertura foi em torno de 0,85 em cenários de simulação complexos com interações entre características ou com dependência temporal, enquanto excedeu 0,90 em cenários mais simples.
A precisão das taxas de especiação e extinção inferidas do modelo BDNN excedeu substancialmente a de modelos alternativos baseados em modelos de nascimento-morte variáveis no tempo ( 52 ) ou no método de cruzamento de limites ( 26 ) em todas as simulações envolvendo taxas dependentes de características. Por exemplo, em simulações apresentando uma característica categórica afetando a especiação e a extinção (cenário 5), o erro BDNN foi de 0,25, comparado a 0,57 com base em modelos alternativos (tabela S1). Enquanto a precisão das taxas BDNN foi semelhante à de modelos alternativos em simulações de taxa constante (cenário 1), o modelo BDNN foi superado por modelos de variáveis no tempo em simulações com taxas sigmoidais variáveis no tempo ou taxas constantes por partes (cenários 2 e 3). A cobertura das estimativas BDNN foi comparável ou substancialmente maior do que aquela baseada em modelos alternativos em todas as simulações, exceto no cenário 3, onde um modelo de nascimento-morte com mudanças de taxa teve melhor desempenho.
Avaliação de variação significativa da taxa
Usamos o coeficiente de variação nas taxas específicas estimadas de linhagem-tempo (ver Materiais e Métodos; Fig. 2 ) para determinar se ele excedeu o que seria esperado sob um processo constante de nascimento-morte. Para esse objetivo, usamos as simulações geradas sob um processo constante de nascimento-morte (cenário 1) para estabelecer um coeficiente de variação de linha de base, acima do qual rejeitamos um modelo constante de nascimento-morte em favor de um modelo BDNN. Escolhemos uma linha de base correspondente a uma especificidade de 95%, ou seja, levando à rejeição errônea de um modelo constante de nascimento-morte em 5% das simulações. Com base nesses limites, a sensibilidade (porcentagem de rejeição correta do modelo de taxa constante em favor de um modelo BDNN) teve uma média de 86,3% para especiação e 90,8% para extinção. A sensibilidade variou entre os cenários, excedendo 80% na maioria das simulações (fig. S2). Foi mais baixa no cenário 9, onde a variação da taxa de condução do traço foi omitida da análise. A variação da taxa inferida sobrepôs-se à simulada, com uma ligeira tendência a subestimar o verdadeiro coeficiente de variação nas taxas específicas do tempo de linhagem (fig. S3).
Identificação dos preditores da variação da taxa
Realizamos várias etapas de pós-processamento para estimar a capacidade do modelo BDNN de identificar os fatores que determinam a variação da taxa. Primeiro, visualizamos a magnitude e a forma dos efeitos de características e variáveis dependentes do tempo nas taxas usando gráficos PD. Eles mostraram uma boa correspondência entre a relação simulada e inferida de taxas com características e variáveis dependentes do tempo ( Fig. 3 e fig. S4), mesmo no caso de efeitos não monotônicos, como ligações em forma de sino entre características contínuas e taxas, e interações entre características e variáveis dependentes do tempo. Embora a forma dos efeitos tenha sido inferida corretamente, embora com variação entre simulações, sua magnitude foi geralmente subestimada. Por exemplo, a diferença de 5 vezes na taxa entre os dois estados de uma característica categórica (cenário 5; fig. S4E) foi inferida com precisão para extinção [média em 100 réplicas: 4,1; Intervalo de confiança de 95%: 2,2 a 6,5], enquanto encontramos apenas uma diferença de 1,9 vezes (1,3 a 2,8) no caso da especiação. A magnitude da variação da taxa atribuída a preditores individuais foi, portanto, menor do que a heterogeneidade detectada entre as taxas específicas de tempo de linhagem marginal ( Fig. 1 ). Isso ocorre porque os gráficos de PD isolam o efeito de preditores individuais, removendo assim parte da variação atribuída pelo modelo a diferentes fatores.
Em seguida, classificamos a importância relativa de características e variáveis dependentes do tempo de acordo com três métricas baseadas em permutações, probabilidades marginais e valores SHAP. Usar um método de classificação de consenso ( 57 ) entre essas métricas provou ser mais robusto do que qualquer métrica individual sozinha (fig. S5). Em cenários onde as taxas variaram ao longo do tempo com base em tempos de mudança de taxa não explicitamente incluídos nas análises (ou seja, cenários 2 e 3), o próprio tempo foi identificado como o preditor mais importante em 76 a 100% das simulações onde um modelo de taxa constante foi rejeitado (fig. S5). Em até 20% das simulações, o fator mais importante foi identificado como um dos autovetores filogenéticos, que refletem inerentemente tanto a relação filogenética quanto o tempo. Quando a variação da taxa foi conduzida por uma característica omitida da análise (cenário 9), o sinal foi principalmente atribuído a autovetores temporais ou filogenéticos (75 a 85% das simulações), embora o modelo de taxa constante tenha sido rejeitado neste caso apenas em cerca de 50% dos casos (fig. S2). No cenário 4, onde as taxas dependem tanto da paleotemperatura quanto de uma característica discreta, a primeira ficou entre os dois principais preditores em 96 a 97% dos casos, enquanto a característica discreta foi corretamente identificada apenas em 26 a 46% das simulações. Em 40 a 57% das simulações, o sinal da característica foi atribuído a autovetores filogenéticos, destacando a dificuldade de distinguir entre o efeito de uma característica evoluindo ao longo da filogenia e a própria filogenia. Em casos onde uma ou mais características conduzem as mudanças nas taxas de especiação e extinção (cenários 6 e 7), a(s) característica(s) correta(s) foram classificadas como os principais preditores em 21 a 90% das simulações, com tempo ou autovetores filogenéticos sendo favorecidos na maioria dos outros casos. Em outros cenários onde uma característica afeta as taxas com efeitos variáveis ao longo do tempo (cenários 5 e 8), a característica correta foi classificada entre os dois principais preditores em 49 a 100% das simulações, enquanto tempo ou autovetores filogenéticos são encontrados entre os dois principais em quase todos os casos.
No geral, descobrimos que o tempo e os autovetores filogenéticos podem capturar efetivamente o sinal de variação de taxa, mesmo quando os verdadeiros preditores são incluídos, mas mais frequentemente quando não são. Isso indica que eles fornecem hipóteses nulas valiosas que permitem a variação de taxa sem implicar um fator específico. Como resultado, a probabilidade de identificar erroneamente uma característica como o preditor de classificação mais alta quando ela não determinou a variação de taxa (ou seja, um falso positivo) foi de 3,1%.
Os motores da diversificação proboscídica
Usamos o modelo BDNN para analisar a dinâmica e os fatores de especiação e extinção em proboscídeos com base em seu rico registro fóssil. Usamos um conjunto de dados publicado consistindo de 2118 ocorrências fósseis e informações ecomorfológicas detalhadas para 175 espécies de proboscídeos ( 58 ). Ele inclui 17 características (por exemplo, tamanho do corpo, morfologia das presas e características dentárias), que foram resumidas por meio de escala multidimensional não métrica (NMDS) em um ecomorfoespaço bidimensional que capturou 93% da variação ecomorfológica total. Além disso, incorporamos informações sobre distribuição geográfica (ou seja, presença na África, Eurásia, Américas ou ilhas), paleotemperatura regional e paleovegetação (aproximada como uma variável binária que reflete o tempo de surgimento de pastagens de habitat aberto em cada área geográfica). Além disso, usamos a flexibilidade do modelo BDNN para incorporar informações filogenéticas (resumidas em dois autovetores) para explicar a não independência de características e capturar os possíveis efeitos de preditores não observados. Por último, incluímos a sobreposição espaço-temporal com os primeiros humanos [codificada para espécies africanas e eurasianas entre o Pleistoceno Inferior e Superior, 1,8 milhões de anos (Ma) a 129 mil anos (ka)] e com os humanos modernos (após 129 ka na África e Eurásia, e no Holoceno 11,7 ka nas Américas) como um potencial preditor de extinção. Como esses limites de tempo foram determinados com base em estágios geológicos (seguindo a classificação usada para as outras variáveis paleoambientais), eles podem apenas aproximar o tempo da evolução humana e expansão para fora da África.
Os coeficientes de variação entre as taxas específicas de tempo de espécies estimadas foram 0,80 para especiação e 1,14 para extinção. Eles excederam fortemente os limites de 0,20 e 0,25, respectivamente, que foram inferidos de simulações baseadas em modelos de taxa constante (ver Materiais e Métodos), indicando assim evidências convincentes de variação significativa de taxa. Resumimos as taxas específicas de tempo de espécies para obter trajetórias gerais de taxas ao longo do tempo e descobrimos que as taxas de especiação declinaram quatro vezes entre o Eoceno Superior (40 Ma) e o fim do Mioceno (5,3 Ma), seguido por um aumento de taxa de magnitude semelhante entre o Plioceno e o Pleistoceno (fig. S6A). As taxas de extinção foram globalmente estáveis até o Pleistoceno Inferior (1,8 Ma), quando aumentaram quatro vezes até o fim do Pleistoceno, seguido por um aumento adicional de seis vezes no Holoceno (11.700 anos atrás; fig. S6B). Também encontramos evidências de forte heterogeneidade nas taxas de especiação e extinção entre as espécies. Por exemplo, no Pleistoceno, inferimos uma taxa de especiação de 2,91 (IC: 1,06 a 4,92) para o elefante anão siciliano Palaeoloxodon falconeri ( Fig. 4A ), que é 8,3 vezes maior do que a taxa do Stegodon huananensis eurasiano (média: 0,35; IC: 0,10 a 0,65). Da mesma forma, as taxas de extinção estimadas variaram 20 vezes entre o mamute lanoso recentemente extinto ( Mammuthus primigenius ; 3,91, IC: 0,42 a 7,35) e a espécie neogênica Gomphotherium angustidens (0,20, IC: 0,02 a 0,46).
Descobrimos que os três fatores mais importantes que afetam a taxa de especiação foram uma característica ecomorfológica, tempo e distribuição geográfica ( Fig. 4 , tabela S2 e dados S1). O gráfico dependente parcial mostra que o surgimento de fenótipos com durabilidade mastigatória dentária aprimorada (NMDS1) levou a um aumento na taxa de especiação ( Fig. 4B ). O modelo BDNN evidenciou ainda que a correlação entre NMDS1 e taxas de especiação se tornou mais forte ao longo do tempo. Embora a amplitude dessa característica ecomorfológica seja semelhante no Oligoceno Inferior e no Mioceno Superior (∼30 e ∼7 Ma, respectivamente), a mudança na taxa de especiação atribuída a essa característica aumentou de 7 a 16 vezes durante esse período. No entanto, nosso proxy para o início e expansão de pastagens de habitat aberto ∼20 a 16 Ma não foi identificado como um fator vital na diversificação proboscídica. Por fim, as espécies que ocorrem em ilhas apresentaram, em média, uma taxa de especiação ∼2,5 vezes maior do que aquelas nas Américas e na Eurásia e um aumento de ∼4,2 vezes em relação às espécies africanas.
Estimamos que as taxas de extinção dos proboscídeos foram mais fortemente afetadas pela sobreposição com os humanos, seguidas por efeitos mais limitados de distribuição geográfica e ecomorfologia ligadas às formas de presas e mandíbulas ( Fig. 4C e tabela S2). Com base em gráficos dependentes parciais, a sobreposição espaço-temporal com os primeiros humanos começando em torno de 1,8 Ma foi associada a um aumento de 5 vezes na taxa de extinção, enquanto o efeito do Homo sapiens moderno no Pleistoceno Tardio e no Holoceno foi associado a um aumento de 17 vezes. A geografia foi conectada a uma variação de 2,7 vezes, com espécies insulares sujeitas a maiores taxas de extinção. As principais morfologias craniodentais (NMDS2) foram identificadas como o terceiro preditor e associadas a uma variação de 3,3 vezes na taxa de extinção, com a maior taxa associada a mandíbulas longas, molares cúspides adaptados à navegação e presas inferiores em forma de pá ( Fig. 4C ). Notavelmente, a paleotemperatura regional ficou em quarto lugar entre os preditores (fig. S7), com as maiores taxas de extinção associadas às baixas temperaturas, enquanto o tempo em si não foi considerado um preditor importante.
Para avaliar a robustez dos nossos resultados empíricos em relação à seleção de preditores, repetimos as análises com base em um subconjunto deles. Especificamente, omitimos pastagens de habitat aberto, substituímos os traços ecomorfológicos por massa corporal discretizada e simplificamos a biogeografia para um traço de insularidade binária. Apesar dessas modificações, obtivemos previsões semelhantes, com taxas consistentes específicas de espécie-tempo inferidas em 98% das espécies ( R 2 = 0,81 para especiação e R 2 = 0,86 para extinção; fig. S8). A insularidade e os humanos ainda permaneceram entre os mais importantes impulsionadores da variação da taxa (tabela S3 e dados S2). No entanto, descobrimos que a massa corporal sozinha não capturou o sinal codificado nos eixos NMDS e não foi identificada como um preditor importante, enfatizando que as adaptações cranianas e mastigatórias tiveram um impacto maior na diversificação proboscídica do que o tamanho sozinho. Em vez disso, esta análise identificou ambos os autovetores filogenéticos entre os principais preditores de mudanças de taxa, refletindo o sinal filogenético relativamente forte na distribuição de ecomorfologias dentro do clado ( 58 ).
DISCUSSÃO
Um modelo de rede neural bayesiana não supervisionada de diversificação
Nossa compreensão da dinâmica macroevolutiva depende de nossa capacidade de estimar com precisão as taxas de especiação e extinção e as variáveis ambientais e características do histórico de vida que podem conduzi-las. Embora seja provável que múltiplos fatores tenham contribuído conjuntamente para moldar a diversidade de clados ao longo do tempo, a maioria das análises macroevolutivas de diversificação se concentrou em preditores individuais, por exemplo, com especiação dependente da diversidade ( 6 ), taxas dependentes do tempo ( 35 ) ou efeitos de paleotemperatura na diversificação ( 41 ). Embora alguns modelos possam incorporar múltiplos preditores de especiação e extinção, eles normalmente fazem suposições simplistas devido a restrições metodológicas. Por exemplo, eles geralmente assumem efeitos monotônicos e aditivos sem interações e permitem variação de taxa ao longo do tempo ( 48 ) ou entre linhagens ( 49 , 59 ). Apresentamos aqui o modelo BDNN, que relaxa muitas das suposições subjacentes a modelos alternativos de diversificação, permitindo simultaneamente variação de taxa entre linhagens e ao longo do tempo. A abordagem BDNN nos permite inferir conjuntamente os efeitos de múltiplos preditores sobre especiação e extinção a partir de dados de ocorrência fóssil, ao mesmo tempo em que considera vieses de preservação. Nossas simulações mostraram que o modelo tem bom desempenho em uma ampla gama de cenários evolutivos, incluindo efeitos não monotônicos de variáveis de tempo e características, e sua interação. Em uma ampla gama de simulações, notamos que o modelo BDNN superou consistentemente abordagens alternativas baseadas em processos de nascimento-morte com mudanças de taxa ( 53 ) e o método de cruzamento de limites ( 27 ), ambos em termos de maior precisão e melhor cobertura. Quando o processo de geração foi um nascimento-morte com mudanças de taxa instantâneas (cenário 3), então o modelo de nascimento-morte com mudanças, como esperado, e o método de cruzamento de limites tiveram melhor desempenho do que o BDNN. Por outro lado, em cenários onde as taxas foram moduladas por uma característica contínua omitida das análises (cenário 9), descobrimos que o modelo BDNN alcançou sensibilidades mais baixas, ou seja, favorecendo conservadoramente um modelo de nascimento-morte constante em uma fração maior de simulações. Mesmo neste caso, no entanto, ele alcançou melhor precisão em comparação aos modelos alternativos de nascimento-morte e de cruzamento de fronteiras. Esses resultados, juntamente com estudos de simulação anteriores explorando o desempenho de diferentes métodos para inferir taxas de especiação e extinção a partir de dados fósseis ( 52 , 60 – 62 ), sugerem que o modelo BDNN pode levar a uma melhoria substancial das taxas estimadas em uma gama de cenários evolutivos.
O modelo BDNN combina preditores baseados em hipóteses, como taxas de especiação mediadas por características ou extinção relacionada ao clima, com variação de taxa que é agnóstica sobre expectativas biológicas, como taxas variáveis no tempo e variação baseada em parentesco filogenético. Essa integração é importante porque a suposição irrealista de um processo de taxa constante como expectativa nula, típica de muitos modelos de nascimento-morte ( 8 , 18 , 48 ), pode levar à identificação espúria de características ou preditores ambientais da dinâmica de diversificação ( 42 , 63 , 64 ). Embora a parametrização complexa do modelo BDNN torne inviável modelar explicitamente o teste e avaliar a significância de cada preditor, por exemplo, por meio de fatores de Bayes, mostramos que ferramentas de IA explicáveis podem ser reaproveitadas para identificar os fatores mais importantes que afetam a diversificação. Isso reforça descobertas recentes que evidenciam que os valores de SHAP derivados de modelos de árvore de regressão podem identificar características vinculadas à sobrevivência em eventos de extinção em massa ( 65 ). Embora esses métodos não nos permitam atualmente descartar quaisquer preditores, descobrimos que eles podem quantificar com precisão sua importância relativa e com desempenho comparável para especiação e extinção. Da mesma forma, altos níveis de precisão foram recuperados em diferentes configurações de simulação, incluindo cenários com características interativas e variáveis dependentes do tempo. Assim como com modelos de regressão múltipla, podemos esperar que o BDNN tenha poder limitado em desembaraçar os efeitos de características ou séries temporais altamente correlacionadas. Em algumas de nossas simulações, autovetores filogenéticos, que são inerentemente correlacionados com características evolutivas preditivas, foram classificados entre os principais fatores. Portanto, é importante selecionar cuidadosamente variáveis e características biologicamente significativas ao configurar a análise ( 64 ) e levar em conta o tempo e a relação filogenética, seja por meio de autovetores filogenéticos ou taxonomia, como uma forma de acomodar preditores não observados ( 42 , 66 ).
Redes neurais foram usadas recentemente para prever taxas de especiação e extinção a partir de dados fósseis ou filogenéticos sob modelos simples dependentes de idade ou característica ( 67 – 69 ). Esses métodos foram aplicados dentro de uma estrutura de aprendizado supervisionado, na qual (i) conjuntos de dados de treinamento rotulados (ocorrências fósseis ou filogenias) são simulados sob taxas conhecidas de especiação e extinção; (ii) redes neurais são treinadas para prever as taxas de geração a partir de características que descrevem os dados simulados; e (iii) os modelos treinados são então aplicados ao conjunto de dados empíricos não rotulados, onde as taxas são desconhecidas. Embora essas abordagens tenham produzido resultados precisos, sua escalabilidade para cenários mais complexos envolvendo variação de taxa e dependência de múltiplas características e variáveis de tempo, conforme testado em nosso modelo BDNN, provavelmente exigiria a geração de dados de treinamento sob uma gama inviável de cenários simulados. Em contraste, o modelo apresentado aqui não é supervisionado, ou seja, não depende de um conjunto de dados de treinamento rotulado. Em vez disso, ele infere as taxas de especiação e extinção com base na probabilidade dos dados fósseis empíricos (não rotulados) sob um modelo de nascimento-morte e preservação. Essa estrutura nos permite inferir a dinâmica de diversificação e como elas são moldadas por diferentes preditores sem a necessidade de definir a priori o intervalo de possível variação de taxa e correlações.
Embora o uso de uma rede neural ofereça alta flexibilidade em termos de quais e quantos preditores podem ser analisados e como eles afetam as taxas, é provável que represente um modelo superparametrizado, que inclui vários parâmetros (pesos) que não estão contribuindo significativamente para explicar os dados. Isso é particularmente importante em um uso frequentista de redes neurais geralmente aplicado a tarefas de aprendizado supervisionado, onde a otimização dos pesos deve ser contrabalançada por técnicas de regularização, por exemplo, usando dropout, e por regras para evitar overfitting, por exemplo, interrompendo a otimização com base no desempenho do modelo em um conjunto de validação ( 70 ). No entanto, em redes neurais bayesianas, os priores nos pesos têm um efeito regularizador, mitigando efetivamente o risco de overfitting no contexto do aprendizado supervisionado ( 71 ). Isso foi ainda mais facilitado em nossa implementação pela introdução de uma camada de saída regularizadora (ver Materiais e Métodos, Eqs. 6 e 7 ). Como resultado da regularização, nosso modelo foi capaz de fornecer taxas altamente precisas e exatas mesmo quando o processo de geração era constante (tabela S1 e Fig. 1 ). Outra vantagem das redes neurais bayesianas sobre sua contraparte frequentista é sua capacidade de fornecer estimativas diretas da incerteza em torno dos parâmetros. Ao amostrar pesos de sua distribuição posterior conjunta em vez de depender de estimativas pontuais, somos capazes de obter estimativas posteriores de taxas de especiação e extinção junto com seus ICs de 95%. Assim, apesar da complexidade da parametrização do BDNN, a superparametrização da rede neural bayesiana não leva a resultados espúrios e é refletida principalmente em ICs inflados em torno das taxas, em oposição a desvios significativos dos valores verdadeiros ( Fig. 3 e fig. S4). Dado que os ICs bayesianos em modelos de nascimento-morte tendem a diminuir com conjuntos de dados maiores ( 30 ), o uso de um modelo complexo como o BDNN, especialmente em conjunto com múltiplos preditores, é, portanto, mais adequado para clados com um rico registro fóssil.
A macroevolução multifacetada dos proboscídeos
Nossa análise do registro fóssil proboscídeo revelou que a diversificação deste clado foi moldada por múltiplos fatores e suas interações. Embora os fatores identificados estejam amplamente de acordo com descobertas anteriores baseadas em análises independentes de alguns deles separadamente ( 58 ), nosso modelo BDNN nos permitiu classificá-los e avaliar suas interações. Especificamente, a especiação dentro do clado foi predominantemente impulsionada por características alimentares [conforme estimado anteriormente; ( 58 )], mas nosso modelo também revelou um efeito temporal não identificado anteriormente que amplificou essa variação impulsionada por características ao longo do tempo. A adaptação alimentar é provavelmente um impulsionador relevante da especiação entre clados animais, desbloqueando o acesso a novos recursos e permitindo a diferenciação ( 47 , 72 – 74 ). Embora não tenhamos encontrado nenhuma evidência de diversificação moldada por grandes mudanças na vegetação (codificadas aqui como o surgimento das primeiras pastagens de habitat aberto), é provável que nosso proxy não tivesse níveis de detalhes espaciais e temporais suficientes para capturar completamente a extensão das mudanças no bioma na diversificação proboscídeo. No entanto, nosso modelo revelou que a ligação entre características associadas à maior flexibilidade alimentar e taxas de especiação se fortaleceu em direção ao fim do Neogeno, à medida que linhagens com dentição altamente durável (caracterizada por molares de coroa alta com um alto número de cristas de esmalte) e modos mastigatórios derivados evoluíram. Interpretamos essa tendência como um efeito indireto da mudança de habitat neste clado e levantamos a hipótese de que a força ambiental poderia ter desempenhado um papel fundamental na construção de tendências evolutivas. Isso poderia ter acontecido inicialmente ao alimentar a inovação fenotípica por meio de seleção eficiente em nível populacional ( 75 ) e, então, ao exacerbar a proliferação de linhagens apresentando as morfologias derivadas mais recentes, como também inferido para outros grupos de ungulados ( 74 ).
O aspecto biogeográfico, particularmente a associação com ilhas, foi considerado importante tanto para a especiação quanto para a extinção, em linha com descobertas em vários outros táxons, como mamíferos, aves, escamados ( 76 – 79 ) e muitos outros clados [por exemplo, 80 )]. Prevê-se que esses efeitos possam surgir devido à dispersão seguida de isolamento, resultando em especiação, enquanto populações menores e recursos limitados em ilhas levam a maiores taxas de extinção ( 81 ). A maior vulnerabilidade de espécies em ilhas ampliou muito os efeitos humanos na biodiversidade das ilhas ( 79 ).
O principal fator determinante da extinção dos proboscídeos foi inferido como sendo a sobreposição com a linhagem humana, alinhando-se com o crescente corpo de evidências indicando o impacto severo dos humanos nas extinções recentes e na megafauna em particular ( 79 , 82 – 86 ). Nosso modelo e o uso de gráficos de PD nos permitem destacar o efeito dos humanos após contabilizar todos os outros fatores, sugerindo que o aumento estimado de 5 a 17 vezes na taxa atribuída aos humanos primitivos e modernos não é influenciado por outros fatores considerados aqui (por exemplo, mudanças climáticas, distribuição geográfica ou características ecomorfológicas). Descobrimos que, embora os humanos exibam o maior impacto nos últimos cerca de 120.000 anos, nossos resultados também apontam para uma influência mais fraca, porém significativa, da linhagem humana em épocas anteriores, apoiando assim outros estudos que sugerem um efeito antropogênico prejudicial duradouro na biodiversidade ( 79 , 87 – 89 ).
Em nossas análises, a mudança climática ficou em quarto lugar entre nossos preditores, sugerindo um impacto potencial do resfriamento levando a uma extinção maior, embora com um efeito inferido comparativamente pequeno. Existem, no entanto, muitos aspectos adicionais, como sazonalidade, precipitação e variação climática em pequena escala, que não podemos incorporar em nossas análises devido às limitações nos dados disponíveis, especialmente no passado mais profundo. Da mesma forma, avaliar o momento exato e a distribuição espacial das mudanças na vegetação [como a expansão da grama C4 ( 90 )] e a extensão da competição da diversificação de outras linhagens de herbívoros, como os bovinos, continua difícil ( 58 ). Assim, o efeito do clima, das mudanças ambientais e das interações bióticas na diversidade dos proboscídeos provavelmente será mais matizado do que podemos detectar aqui com os dados disponíveis. Com o crescimento contínuo na disponibilidade de ocorrências fósseis digitalizadas e características morfológicas ( 91 – 93 ), juntamente com o progresso na modelagem e medição do paleoclima e da evolução dos paleoambientes ( 94 – 96 ), esperamos obter uma imagem cada vez mais precisa da dinâmica evolutiva passada e recente. Neste contexto, nosso modelo BDNN abre caminho para uma avaliação poderosa dos processos de especiação e extinção, bem como das forças que impulsionam a ascensão e queda dos clados.
MATERIAIS E MÉTODOS
A estrutura PyRate
Nosso modelo é baseado na estrutura Bayesiana implementada no PyRate para analisar dados de ocorrência de fósseis ( 52 ). A estrutura modela os dados de ocorrência como o resultado de um processo de Poisson não homogêneo representando preservação e amostragem e um processo de nascimento-morte descrevendo a distribuição dos tempos de origem e extinção. Usando um algoritmo MCMC, o PyRate estima a partir da distribuição posterior conjunta: (i) as taxas de preservação e sua variação ao longo do tempo e entre linhagens, (ii) os tempos de origem e extinção de todas as linhagens e (iii) as taxas de especiação e extinção.
Diferentes processos de nascimento-morte podem ser usados para testar a variação da taxa ao longo do tempo com base em modelos constantes por partes ( 30 ) ou correlação com variáveis dependentes do tempo ( 84 ). Alternativamente, a variação da taxa entre linhagens pode ser estimada por meio de correlações com uma característica específica da linhagem ou dependente dos estados de uma característica categórica ( 49 ). No entanto, esses modelos atualmente não podem ser misturados para, por exemplo, incluir efeitos dependentes do tempo e da característica e assumir correlações lineares ou exponenciais simples entre as taxas e os preditores (séries temporais ou características) sem interações.
O modelo BDNN
Desenvolvemos um modelo que pode acomodar múltiplos efeitos dependentes de tempo e característica e suas interações, sem limitar os efeitos a funções simples. As características podem ser de natureza diferente, como características quantitativas, ordinais e categóricas. Espera-se que variáveis categóricas com mais de dois estados sejam codificadas one-hot em características binárias. A taxa no tempo t para a linhagem i é modelada por meio de uma rede neural feed-forward tomando como entrada o tempo e os caracteres e retornando uma taxa de especiação específica do tempo de linhagem λ ( i , t ) como saída. Uma rede independente com a mesma arquitetura é usada para obter taxas de extinção específicas do tempo de linhagem μ ( i , t ) . A entrada da rede neural é uma matriz de tamanho J , onde J é o número de preditores
(1)
(2)
Da mesma forma, a segunda camada oculta é uma matriz de tamanho L (2)onde k ∈ {1, …, L (2) } indica os nós ocultos da segunda camada e W (2) é uma matriz de pesos L (1) × L (2) . Uma terceira camada recebe h (2) como entrada e retorna uma taxa de linha de base específica da espécie-tempo, por exemplo, para taxa de especiaçãoonde W (2) é uma matriz de pesos L (2) × 1 e σ é a função de ativação softPlus ( 98 )que garante que a saída da rede, que é a taxa, seja um número positivo. Por último, a taxa de base é transformada na taxa de especiação ou extinção específica da espécie-tempo, por meio de uma camada de saídaonde φ é uma função regularizadoracom t reg ∈ [0,1] como um parâmetro regularizador e E( · ) sendo a função média aritmética que faz a média das taxas de base em todos os táxons e intervalos de tempo. Esta função não altera a taxa média ao longo do tempo e táxons em comparação com as taxas de base [ou seja, E(λ b ) = E (λ)], mas tem a propriedade de encolher as taxas em torno da média para valores de t reg mais próximos de 0, enquanto as deixa inalteradas se t reg ≈ 1. As taxas se tornam constantes ao longo do tempo e entre espécies quando t reg = 0. A rede neural pode ser configurada com diferentes números de camadas e nós ocultos e pode incluir um nó de polarização (ou interceptação) na última camada oculta. Seus parâmetros (os pesos) são compartilhados entre todas as espécies no conjunto de dados e todos os intervalos de tempo, e duas redes independentes são usadas para obter taxas de especiação e extinção com parâmetros
(3)
(4)
(5)
(6)
(7)
Ao contrário de outros modelos de nascimento-morte implementados no PyRate, os parâmetros amostrados por meio do MCMC sob este modelo não são diretamente as taxas de especiação e extinção, mas sim os pesos das duas redes neurais, W λ e W μ e o parâmetro regularizador t reg , do qual as taxas são derivadas para cada espécie e cada intervalo de tempo ( Eqs. 1 a 4 ). O número de parâmetros no modelo depende do número de preditores J (que no mínimo incluem tempo, mas podem incorporar características categóricas e contínuas e variáveis dependentes do tempo) e da arquitetura da rede (ou seja, número de camadas e nós). As redes neurais usadas no modelo BDNN não são supervisionadas, pois não são treinadas contra taxas de especiação e extinção de verdade básica, mas usadas dentro da estrutura bayesiana do PyRate, na qual as taxas de especiação e extinção são usadas para calcular a probabilidade dos dados de ocorrência fóssil observados (junto com os parâmetros do processo de preservação). A probabilidade de uma linhagem extinta i com tempo de origem estimado s i e tempo de extinção e i e preditores x i é baseada no processo de nascimento-morte amostrado ( 52 )que para linhagens existentes (ou seja, com e i = 0) se reduz a
(8)
(9)
O modelo BDNN usa o mesmo algoritmo implementado em outros modelos PyRate ( 52 ) para amostrar os tempos de origem e extinção de todas as espécies, taxas de preservação e os pesos das redes neurais de sua distribuição posterior conjunta. Atribuímos distribuições normais padrão, 𝒩(0,1), como priores regularizadores nos pesos ( 71 ), e amostramos esses pesos por meio do MCMC com kernels de proposta distribuídos normalmente. Usamos um prior exponencial truncado no parâmetro regularizador t reg ∼ Exp T = 1 (1), de modo que a maior probabilidade prior foi atribuída a t reg = 0, favorecendo assim uma hipótese nula de taxas constantes, enquanto a menor probabilidade prior diferente de zero foi atribuída a t reg = 1, o que corresponde à ausência de regularização das taxas de linha de base.
Após análises preliminares usando conjuntos de dados simulados e dados proboscídeos, encontramos diferenças insignificantes entre arquiteturas de redes neurais e escolhemos usar duas camadas ocultas com 16 e 8 nós. Essas configurações funcionaram bem para nossos conjuntos de dados com 200 a 300 espécies e 6 a 11 preditores. No entanto, encorajamos os usuários a experimentar diferentes arquiteturas de rede para avaliar a robustez dos resultados a diferentes configurações de modelo.
Pós-processamento
Implementamos uma série de etapas de pós-processamento usando a saída de uma análise BDNN para identificar se as características e os preditores de variáveis temporais afetam as taxas de especiação e extinção, a direção de seus efeitos e sua classificação de importância. Isso é importante porque a comparação de modelos normais para selecionar os preditores (por exemplo, usando fatores de Bayes) é impraticável devido ao alto número de combinações possíveis.
Usamos gráficos de PD ( 53 ) para visualizar a magnitude e a forma do efeito inferido de características e preditores de variáveis temporais nas taxas de especiação e extinção. Os gráficos de PD são amplamente usados em aprendizado de máquina interpretável ( 54 ) e mostram o efeito marginal de um ou dois preditores no resultado de um modelo de aprendizado de máquina por meio da média dos preditores restantes.
Como as redes neurais não fornecem diretamente uma medida da importância do preditor, implementamos três abordagens para classificar os preditores, quantificando diferentes aspectos da importância: influência no ajuste do modelo por meio da permutação de características, diferenças de taxa confiáveis em gráficos de PD e a mudança induzida na taxa de especiação ou extinção por um preditor por meio de SHAPs. Análises preliminares usando registros fósseis simulados com preditores conhecidos de variação de taxa mostraram que eles nem sempre podiam ser identificados corretamente usando uma única abordagem. Consequentemente, classificamos a importância dos preditores com base nas três abordagens e usamos o algoritmo QUICK ( 57 ) para obter a classificação de consenso entre eles como uma medida geral de importância. Classificamos separadamente os efeitos individuais e a força da interação entre dois preditores porque fortes efeitos individuais podem imprimir interações e, portanto, ofuscar as medidas de importância de outros preditores individuais.
Para a primeira medida de importância do preditor, quantificamos a diminuição na probabilidade de nascimento-morte após permutar os valores do preditor. Isso difere ligeiramente da implementação original da permutação de características ( 99 ) na segmentação do ajuste do modelo em vez do erro de previsão (que é tipicamente expresso como log-verossimilhança negativa). Uma interação entre dois preditores é considerada importante se a diminuição no ajuste do modelo após permutar ambos for maior do que a soma das reduções individuais quando apenas um dos preditores é permutado.
Segundo, derivamos diferenças posteriores baseadas em PD em taxas como uma medida de importância do preditor. Primeiro, identificamos o máximo e o mínimo da taxa média (ou seja, a média sobre gerações MCMC) ao longo do intervalo de valores do preditor. Então, quantificamos com que frequência as taxas no ponto máximo excederam aquelas no ponto mínimo nas amostras MCMC.
Por último, usamos SHAPs ( 55 ) para quantificar a contribuição de cada preditor para as taxas de especiação e extinção. Especificamente, obtivemos valores de SHAP usando o método baseado na integral de Choquet κ-aditiva ( 100 ) porque ele é computacionalmente mais eficiente do que o SHAP do Kernel original e quantifica a importância da interação entre dois preditores. Os valores de SHAP também indicam até que ponto cada preditor tem alavancagem na taxa de especiação e extinção de cada espécie. Isso auxilia na interpretação de taxas específicas de espécies porque, em uma rede neural, como usada para o modelo BDNN, a direção e a força da influência de um preditor em taxas individuais não podem ser facilmente obtidas a partir dos pesos da rede inferida.
Cenários de simulação
Simulamos conjuntos de dados fósseis e características de espécies em nove cenários de nascimento-morte onde as taxas de especiação e extinção são influenciadas por características, tempo, mudança ambiental ao longo do tempo ou sua combinação. Todos os cenários de diversificação foram simulados para frente no tempo em passos de tempo discretos de 0,01 Ma, começando na mesma idade raiz de 35 Ma. Para evitar conjuntos de dados extremamente pobres ou ricos em espécies, visamos uma saída de 200 a 300 espécies existentes ou extintas durante esse período. Em todos os cenários de diversificação, as espécies foram caracterizadas por uma característica contínua e uma característica discreta com dois estados, mesmo em cenários onde a especiação e a extinção não foram afetadas.
Simulamos a evolução do traço contínuo, começando com um valor de zero, por uma caminhada aleatória imparcial com uma taxa de σ 2 = 0,02 (ou seja, o equivalente ao movimento browniano em tempo contínuo). Neste modelo padrão de evolução de traço contínuo, uma espécie descendente herda o fenótipo completo de seu ancestral; portanto, não há mudança de traço em eventos de especiação cladogenética. No entanto, as espécies no registro fóssil são tipicamente definidas com base em sua morfologia [por exemplo, 101 )], o que requer algumas mudanças fenotípicas durante a especiação e relativamente poucas mudanças anagenéticas. Isso corresponde à visão microevolucionária clássica da divergência de traço na especiação ( 25 , 102 ). Embora a interconexão entre o conceito de espécie morfológica e o fenótipo torne difícil testar este modo de evolução de traço a partir do registro fóssil, mudanças de traço em eventos de especiação foram demonstradas ao longo da filogenia de espécies existentes ( 103 ). Adicionamos uma mudança fenotípica instantânea desenhando um novo valor de característica para as espécies descendentes a partir de uma distribuição normal com uma média igual ao valor da característica do ancestral e um desvio padrão (DP) de
Os cenários de diversificação diferiram na dinâmica das taxas de especiação e extinção ao longo do tempo e se estas dependem das características da espécie ou da paleotemperatura (fig. S1). O primeiro cenário apresentou taxas de especiação e extinção constantes definidas como λ = 0,2 e μ = 0,1, respectivamente.
O cenário 2 incluiu uma mudança sigmoidal nas taxas onde a especiação diminuiu de 0,4 no início do processo de diversificação para 0,1 no presente, enquanto a extinção aumentou ao longo do tempo de 0,05 para 0,4. As taxas seguiram uma função logística com ponto médio definido em 20 e 15 Ma (para especiação e extinção, respectivamente) e inclinação definida em 0,5.
O cenário 3 impôs duas mudanças instantâneas de taxa. A taxa de especiação foi definida para diminuir de 0,4 para 0,1 em 20 Ma e de 0,1 para 0,01 em 10 Ma. A taxa de extinção atingiu o pico de 0,3 de 15 a 10 Ma, enquanto antes a taxa era de 0,05 e depois de 0,01.
No cenário 4, a temperatura global ao longo do tempo exerce influência sobre a especiação e a extinção. Para explorar a capacidade da rede neural de capturar uma relação não monotônica entre um preditor e uma taxa e uma interação entre dois preditores de taxa, simulamos um link em forma de sino entre temperatura e taxa para o primeiro estado da característica categórica e uma curva de sino invertida para o segundo estado (ou seja, funções gaussianas; veja a fig. S1D). Derivamos a trajetória da temperatura para o período de tempo da simulação a partir de dados de isótopos ( 104 ) usando as equações fornecidas por Hansen et al. ( 105 ) e a escalamos para um DP de 1. O pico do sino foi definido para a temperatura média ao longo do período de tempo da simulação (ou seja, 17,7 °C) e correspondeu a uma taxa de especiação de base de 0,5 e uma taxa de extinção de 0,4. Parametrizamos a mudança de taxa com temperaturas mais quentes ou mais frias por meio de um DP dado para a função gaussiana entre a temperatura e a taxa. Efetivamente, isso transforma o sino para ser mais pontiagudo ou mais plano e menos divergente da taxa de base com um DP alto ou baixo, respectivamente. Definimos o DP da função gaussiana para 1,2. Assim, um aumento ou diminuição de 2,1 °C da média da trajetória de temperatura (ou seja, 1 DP) resultou em uma taxa 50% menor do que a linha de base, que então diminui para 3% quando a diferença de temperatura é de 2 DPs, e assintoticamente atinge uma taxa de zero com desvios de temperatura ainda maiores. Para o segundo estado, a relação entre temperatura e taxa é invertida, com uma taxa de zero na temperatura média e um aumento assintoticamente em direção à taxa de base por temperaturas mais altas e mais baixas.
No cenário de diversificação dependente de característica 5, espécies apresentando o primeiro estado da característica categórica foram caracterizadas por uma taxa de especiação de 0,1 e uma taxa de extinção que aumentou ao longo dos 35 Ma simulados de 0,01 para 0,1, resultando em taxas iguais no presente. Espécies do segundo estado tiveram taxas 5 vezes maiores, o que significa que a força do aumento na extinção foi dependente do estado e, portanto, constitui uma interação entre a característica categórica e o tempo.
No cenário 6, o traço contínuo da espécie afetou sua taxa de especiação e extinção por meio de uma função gaussiana, resultando em que ambas as taxas são mais altas em um valor de traço (arbitrário) de 0 e decrescente com o aumento de valores positivos ou negativos. Definimos a função gaussiana de modo que uma mudança de 1 ou 2 DPs do valor inicial do traço de 0 reduziu a taxa de linha de base (aqui definida como 0,5 para especiação e 0,4 para extinção), em 70 e 20%, respectivamente.
No cenário 7, as taxas de especiação e extinção foram determinadas pela interação entre características contínuas e categóricas. Enquanto as taxas para espécies com o primeiro estado da característica categórica foram obtidas da mesma função em forma de sino do cenário 6, espécies do segundo estado tiveram uma relação invertida entre característica contínua e taxa: valores de característica de 0 definiram uma taxa de 0, e taxas de linha de base são assintoticamente abordadas com valores de característica cada vez menores ou maiores.
O cenário 8 apresentou uma interação entre o tempo e a característica contínua. Estabelecemos em 15 Ma o tempo da troca entre uma função em forma de sino e uma função de sino invertida, correlacionando a característica contínua com as taxas de especiação e extinção.
Por último, o cenário 9 foi baseado nos conjuntos de dados simulados sob o cenário 6, ou seja, com taxas dependentes de uma característica contínua, que foi, no entanto, omitida da análise BDNN subsequente. Uma característica adicional foi incluída com base em um sorteio aleatório de uma distribuição normal padrão, mantendo o mesmo número de preditores.
Simulamos 100 conjuntos de dados de ocorrência fóssil em cada cenário, assumindo um processo de Poisson variável no tempo de fossilização e amostragem, com taxas independentes desenhadas para cada estágio geológico aleatoriamente de q i ∼ 𝒰[log(0.5), log(5)]. Implementamos adicionalmente a variação da taxa de preservação entre espécies desenhando multiplicadores específicos de espécies de uma distribuição gama com parâmetros de forma e taxa definidos iguais e desenhados de α ∼ 𝒰[log(0.5), log(5)]. A taxa de preservação para cada espécie é então obtida como o produto entre as taxas de base q e os multiplicadores específicos de espécies. Este modelo reflete o processo de Poisson variável no tempo com variação de taxa entre linhagens implementado em PyRate ( 52 ), com uma parametrização semelhante à que estimamos do conjunto de dados proboscídeo. Enquanto a evolução estocástica da característica contínua foi rastreada para cada espécie ao longo da simulação, o valor usado como entrada para a inferência BDNN foi definido como a média da característica apenas nos tempos de ocorrência amostrados da espécie, refletindo assim o tipo de informação morfológica observável em conjuntos de dados empíricos. Aproximamos a relação filogenética das espécies decompondo as árvores simuladas em autovetores com base na matriz de distância filogenética em pares entre as pontas da árvore. Essa abordagem, conforme implementada no pacote PVR 0.3 ( 106 ) para o ambiente de programação estatística R 4.3 ( 107 ), resulta em uma representação numérica da árvore capturando tanto a topologia quanto os tempos de ramificação ( 56 ).
Os dados de entrada usados em nossas análises incluíram (i) as ocorrências fósseis, (ii) os traços simulados e (iii) o próprio tempo e os dois primeiros autovetores filogenéticos. No cenário 4, incluímos adicionalmente a paleotemperatura como uma variável dependente do tempo. Por outro lado, os tempos de mudança de taxa implementados nos cenários 2, 3 e 8 e o traço preditivo no cenário 9 não foram incluídos entre os preditores. Esperava-se, portanto, que o tempo e os autovetores filogenéticos capturassem a variação de taxa impulsionada por traços não observados ou mudanças de taxa.
Análise de dados simulados
Analisamos a capacidade do modelo BDNN de capturar os efeitos simulados sobre a especiação e a extinção comparando os resultados com aqueles obtidos sob o modelo birth-death-shift (BDS), que infere mudanças instantâneas de taxa ao longo do tempo e é amplamente usado em análises exploratórias de dinâmica de diversificação a partir de dados fósseis [por exemplo, ( 52 )]. Para ambos os modelos, executamos 1.000.000 de iterações MCMC e amostramos cada 1.000ª iteração após um burn-in de 25% para obter a distribuição de parâmetros posteriores. Também analisamos os dados sob o método boundary-crosser ( 26 , 108 ), fornecendo uma comparação do desempenho do modelo BDNN contra uma abordagem baseada em diferentes modelos e implementações subjacentes. Obtivemos taxas de especiação e extinção boundary-crosser usando o pacote R divDyn 0.8.2 ( 109 ) usando intervalos de tempo de 2 milhões de anos. Além disso, realizamos análises em 100 conjuntos de dados bootstrap para estimar intervalos de confiança em torno das estimativas de taxas ( 26 ), que usamos para aproximar a cobertura como a fração de simulações nas quais as taxas verdadeiras foram incluídas no intervalo de confiança.
Quantificamos a precisão geral das taxas de especiação e extinção específicas de espécies inferidas como o erro relativo absoluto mediano, que foi calculado comoonde r i é a taxa real de origem ou extinção de uma espécie simulada,
(10)
Registro fóssil, características e paleoambiente dos proboscídeos
Utilizamos um conjunto de dados de ocorrência fóssil de proboscídeos recentemente compilado ( 58 ), do qual excluímos 10 espécies. Cinco espécies pertenciam a um grupo monofilético e eram representadas apenas por ocorrências únicas com mais de 10 milhões de anos de idade do que as outras espécies, resultando em taxas altamente incertas antes de 40 milhões de anos em uma análise de diversificação exploratória ( 58 ). As cinco espécies restantes não tinham informações sobre características. O conjunto de dados final incluiu 2.118 ocorrências para 175 espécies de proboscídeos. Para levar em conta as incertezas de idade em ocorrências fósseis, reamostramos as idades de ocorrência aleatoriamente a partir de distribuições uniformes abrangendo seus intervalos temporais e geramos 10 conjuntos de dados replicados ( 52 ).
Cantalapiedra et al. ( 58 ) resumiram 17 características ecomorfológicas por uma ordenação em dois eixos NMDS, que usamos em nossa análise BDNN. A durabilidade mastigatória dentária é capturada pelo primeiro eixo, que se correlaciona amplamente com a mudança de uma dieta especializada em navegação para uma dieta mais generalizada, capaz de processar uma proporção maior de grama e até mesmo madeira. O segundo eixo de características resumiu modificações craniodentais, como de um crânio plano para um crânio curto e alto. Embora a adequação de variáveis que resumem múltiplas características por meio de técnicas de redução de dimensionalidade (por exemplo, NMDS e análise de componentes principais) em análises evolutivas permaneça discutível, problemas potenciais se relacionam com nossa capacidade de modelar como essas próprias características resumidas evoluem ( 110 , 111 ). Em contraste, aqui, nós as usamos como preditores de especiação e variação da taxa de extinção sem tentar inferir os modos de sua evolução.
Devido à história evolutiva compartilhada entre espécies, características e distribuição geográfica exibem inércia filogenética ( 112 ). Não contabilizar a não independência estatística resultante poderia inflar o sinal de taxas de especiação e extinção dependentes de características ou geografia. Assim, incorporamos os dois primeiros autovetores filogenéticos, calculados a partir de uma árvore de proboscídeos vivos e extintos ( 58 ) como um proxy para parentesco filogenético, e os incluímos como características específicas da espécie, além de características e distribuição geográfica na inferência BDNN. Para contabilizar incertezas filogenéticas, recalculamos os autovetores para uma árvore diferente (amostrada aleatoriamente de sua distribuição posterior) para cada uma das 10 réplicas com idades fósseis reamostradas.
Também compilamos a distribuição geográfica para todas as espécies de proboscídeos com base em ( 58 ), com a distribuição categorizada em África, Américas e Eurásia ou se a espécie é endêmica de uma ilha. Com base na distribuição geográfica de cada espécie, obtivemos uma série temporal individual de condições paleoambientais. Como nossa primeira medida paleoambiental, registramos a origem de habitats abertos e dominados por gramíneas como um preditor binário com base em cronogramas específicos de (sub)continentes inferidos de registros paleobotânicos ( 90 ). Definimos o surgimento desses habitats na África em 15,97 Ma e na Eurásia em 20,44 Ma, exceto para espécies indianas para as quais usamos um início em 11,63 Ma. Enquanto na América do Sul os habitats abertos e dominados por gramíneas datam do Eoceno, usamos a idade norte-americana de 23,03 Ma para as Américas porque todos os proboscídeos sul-americanos são muito mais jovens. Em segundo lugar, reconstruções espacialmente explícitas do clima ao longo do tempo foram usadas para obter séries temporais específicas de espécies de paleotemperatura com base nas faixas de distribuição das espécies. Como os modelos de circulação climática (CCM) cobrem apenas uma fração da história de diversificação proboscídica, usamos uma reconstrução climática derivada da vegetação ( 94 ) com uma resolução temporal de ∼0,17 Ma e uma resolução espacial de 2° × 2°. Para o Último Máximo Glacial e o Último Interglacial, as temperaturas foram extraídas dos CCMs ( 113 , 114 ) porque uma resolução temporal mais alta é necessária para comparar os efeitos potenciais da temperatura e dos humanos nas taxas de extinção. Como há apenas um número limitado de ocorrências fósseis por intervalo de tempo, usamos todas as ocorrências fósseis de espécies, independentemente de sua idade, para aproximar as faixas de distribuição das espécies. Construímos um casco convexo ligando todos os países de um continente onde fósseis das espécies foram encontrados e não tentamos rastrear mudanças na distribuição ao longo do tempo. A paleoposição dos cascos convexos foi reconstruída para cada intervalo de tempo dos mapas paleo usando a interface rgplates 0.2.1 ( 115 ) para o aplicativo de desktop GPlates ( 116 ) e usada para extrair a paleotemperatura média através da distribuição aproximada com o pacote R terra 1.7.29 ( 117). Todos os dados espaciais foram transformados na projeção do mapa de área igual de Mollweide para evitar viés nas temperaturas devido à distorção de coordenadas geográficas em alta latitude. Usamos os seguintes limites para agrupar as paleotemperaturas: Oligoceno (33,9 Ma), Burdigaliano (20,44 Ma), Langhiano (15,97 Ma), Serravaliano (11,63 Ma), Tortoniano (7,246 Ma), Messiniano (5,333 Ma), Plioceno (2,58 Ma), Calábria (1,8 Ma), a transição do Pleistoceno Médio (∼0,8 Ma), Pleistoceno Superior (0,129 Ma) e o Holoceno (0,0117 Ma). Esses limites definem quando as taxas de especiação e extinção podem mudar ao longo do tempo. Embora tenhamos usado compartimentos com uma largura de 1 Myr em nossas simulações, optamos aqui por uma resolução temporal crescente em direção ao presente porque ela reflete a maior confiança na reconstrução da paleotemperatura, a menor incerteza nas idades fósseis e a resolução necessária para informar o modelo de aparência humana. Em um estudo anterior sobre diversificação proboscídica ao longo do tempo, não houve evidências de mudanças de taxa antes do Messiniano ( 58 ) e, portanto, a resolução relativamente grosseira usada aqui provavelmente não ocultará mudanças temporais na dinâmica de diversificação. Codificamos o tempo de origem ou chegada dos humanos primitivos e modernos em cada continente e ilhas associadas nos respectivos compartimentos de tempo usando dois níveis. Especificamente, consideramos a sobreposição proboscídica com os humanos modernos como onipresente no compartimento de tempo mais recente (ou seja, o Holoceno). Para espécies africanas e eurasianas, também incluímos a sobreposição com os humanos modernos no penúltimo compartimento (ou seja, o Pleistoceno Superior) e com os humanos primitivos em todos os outros compartimentos de tempo com menos de 1,8 Ma.
Análise BDNN dos dados proboscídeos
Transformamos em z os traços contínuos, autovetores filogenéticos e paleotemperatura para uma média de zero e um DP de 1. Para todas as 10 réplicas com idades fósseis reamostradas, executamos 50.000.000 iterações MCMC e amostramos cada 10.000ª iteração após um burn-in de 25% no PyRate. Usamos um modelo de preservação que permite heterogeneidade na amostragem de fósseis entre espécies e ao longo do tempo ( 52 ). As análises foram concluídas em aproximadamente 48 horas em unidades de processamento central AMD padrão. Após a inferência, aplicamos todas as etapas de pós-processamento descritas acima usando 1000 amostras da distribuição posterior para cada réplica. Além disso, conduzimos uma análise combinada concatenando 100 amostras de cada uma dessas réplicas.
Para determinar se a heterogeneidade da taxa inferida pelo modelo BDNN excedeu a expectativa sob um processo de nascimento-morte de taxa constante, simulamos 100 conjuntos de dados que correspondem aos dados empíricos em termos de número de características, idade do clado e número aproximado de espécies (intervalo 175 ± 52). Nós determinamos o quantil de 95% do coeficiente de variação entre as taxas específicas de espécie-tempo e o comparamos com o coeficiente de variação obtido do conjunto de dados proboscídeos.
Também reanalisamos o conjunto de dados sob um conjunto mais simples de variáveis, incluindo apenas tempo, autovetores filogenéticos, paleoclima, insularidade, massa corporal (discretizada em oito compartimentos) e sobreposição humana. Isso nos permitiu avaliar a consistência entre as duas análises e se a massa corporal sozinha poderia capturar a maior parte do sinal encapsulado pelos dois eixos NMDS.
Nota adicionada na prova : Após a aceitação, os autores tomaram conhecimento de um artigo publicado recentemente ( 118 ). Esta referência adicionada acrescenta um exemplo empírico aos fundamentos teóricos.
Agradecimentos
Os recursos computacionais foram fornecidos pelo de.NBI Cloud financiado pelo BMBF dentro da Rede Alemã para Infraestrutura de Bioinformática (de.NBI) e UBELIX ( www.id.unibe.ch/hpc ), o cluster HPC na Universidade de Berna. Agradecemos aos revisores anônimos pelo feedback construtivo. As silhuetas proboscídicas foram traçadas pelo JLC a partir de reconstruções de A. Larramendi, M. Antón, M. Romano, cisiopurple (DeviantArt) e Roman Uchytel e são disponibilizadas sob uma licença CC-BY 4.0 ( https://creativecommons.org/licenses/by/4.0/ ) no PhyloPic ( http://phylopic.org/ ).
Financiamento: TH e DS receberam financiamento da Swiss National Science Foundation (PCEFP3_187012). DS também reconhece o financiamento do Swedish Research Council (VR: 2019-04739) e da Swedish Foundation for Strategic Environmental Research MISTRA dentro da estrutura do programa de pesquisa BIOPATH (F 2022/1448). JLC foi apoiado pelo Talent Attraction Program do subsídio do Governo de Madri 2017 T1/AMB5298.
Contribuições dos autores: Conceitualização: Todos os autores. Desenvolvimento do método: TH e DS Compilação de dados: JLC e TH Análises: TH Redação: Todos os autores.
Conflitos de interesse: Os autores declaram não ter conflitos de interesse.
Disponibilidade de dados e materiais: O modelo BDNN e o tutorial correspondente estão disponíveis em https://github.com/dsilvestro/PyRate . Também desenvolvemos a biblioteca Python BDNNsim (disponível em https://github.com/thauffe/BDNNsim ) para realizar simulações estocásticas nas quais as taxas de especiação e extinção são determinadas por características em evolução, tempo e trajetórias paleoambientais. A biblioteca gera arquivos de entrada para PyRate após simular amostragem fóssil, mas também fornece a árvore filogenética correspondente, que pode ser usada para outras análises comparativas. Todos os dados necessários para avaliar as conclusões do artigo são depositados no repositório GitHub ( https://github.com/thauffe/BDNN ) e Zenodo ( https://doi.org/10.5281/zenodo.11652106 ).
