A seleção de equilíbrio a longo prazo contribui para a adaptação na Arabidopsis e seus parentes

onde T se refere ao tempo de divergência dos dois gêneros e N A/N C se refere ao tamanho da população de A. thaliana/C. rubella, respectivamente. De acordo com nossas estimativas sob o Modelo M1, levando em consideração as mudanças no tamanho da população, esta probabilidade de identidade por descendência é da ordem de 10-9. Dado que temos < 39.275.210 locais alinhados entre as duas espécies na região gênica, esperamos que o número total de TSPs neutras seja < 1 apenas por deriva genética.

Assumimos o acasalamento aleatório em nosso modelo; entretanto, ambas as espécies são auto-suficientes e a estrutura populacional provavelmente existe dentro das espécies. No entanto, eventos demográficos recentes devem ter relativamente pouco efeito, uma vez que requeremos eventos de coalescência profunda por acaso em ambas as espécies na mesma região do genoma. Como ilustrado no estudo anterior, mesmo uma estrutura populacional profunda dentro dos humanos modernos deve ter um efeito mínimo sobre a probabilidade. Neste estudo, ambas as espécies têm um histórico de predominância de outcrossing. A. thaliana passou de outcrossing para selfing apenas um milhão de anos atrás e C. rubella transitou muito mais recentemente . Mesmo como espécie autóctone, a taxa de outcrossing das populações locais chega a 14,5% . Portanto, estruturas populacionais, se existentes, são improváveis de persistir durante uma longa escala de tempo e seu impacto na probabilidade pode ser ignorado.

Identificação de polimorfismos transespécies sob seleção de equilíbrio

TSPs podem ser distinguidas de mutações neutras porque regiões sob seleção de equilíbrio a longo prazo se agrupam por alelo, e não por espécie . Portanto, nós focamos em seguida nos 433 genes candidatos com SNPs confiáveis compartilhados na região de codificação e examinamos os haplotipos cobrindo cada SNP bi-alélico compartilhado com MAF > 0,05 nas regiões genéticas.

Para estimar o comprimento de cada segmento portador de um sinal de TSPs, nós usamos uma fórmula derivada anteriormente que se baseia em grande parte na taxa de recombinação. Do ponto de vista da coalescência, tal segmento não é decomposto por recombinação até que todas as linhagens da mesma classe alélica coalesçam ao seu mais recente ancestral comum na população ancestral . Adotando uma taxa de recombinação de 3,6 cM/Mb para ambas as espécies, o comprimento do segmento foi extremamente curto, ou seja, apenas vários pares de bases, teoricamente. Dado que ambas as espécies surgiram recentemente dos seus respectivos progenitores outcrossing e a taxa de recombinação efectiva poderia ser muito maior no passado, o comprimento esperado pode ser ainda mais curto. Esta estimativa sugere, sob as circunstâncias neutras do nosso sistema, que é altamente difícil descobrir qualquer segmento sem uma ruptura de recombinação. No entanto, quando a seleção de equilíbrio existe, a seleção pode suprimir a recombinação na região ao redor. Portanto, o comprimento do segmento deve ser maior do que a estimativa teórica sob um modelo neutro. Assim, nós escaneamos a região gênica usando uma janela de 100 bp e um passo de 1-bp.

Nos 433 genes candidatos, detectamos 975 SNPs bi-alélicos compartilhados (incluindo SNPs exônicos e intrônicos com MAF > 0,05). Similar aos estudos anteriores, procuramos em seguida janelas cobrindo pelo menos dois dos 975 SNPs que estão em forte desequilíbrio de ligação (r 2 > 0,5) nas duas espécies entre as janelas qualificadas (alinhadas a um mínimo de 95% do comprimento; ver “Métodos” para detalhes) para identificar árvores alélicas. Essas restrições podem reduzir muito os falsos positivos e produzir árvores alélicas, se elas existirem, com alta resolução. Finalmente, identificamos janelas de cinco genes, AT1G35220, AT2G16570, AT4G29360, AT5G38460 e AT5G44000, envolvendo dez locais, como candidatas TSPs sob seleção de equilíbrio a longo prazo (arquivo adicional 3: Figura S1). Nenhum dos cinco genes ortológicos que encontramos aqui está correlacionado com a variação do número de cópias (CNV) e todos eles têm apenas um acerto quando os comparamos com as referências das duas espécies, respectivamente (ver “Métodos” para detalhes).

Para verificar as regiões identificadas, primeiro determinamos todos os haplótipos nas regiões identificadas de cada população e ressequenciamos acessos representativos para cada haplótipo (ver arquivo adicional 1: Tabela S4 para os primers). Como esperado, todos os locais de TSP candidatos nos cinco genes foram validados e as sequências das duas espécies nas regiões candidatas agrupadas por alelo, em vez de espécies (Fig. 4). No gene AT1G35220, os dois sites candidatos à TSP estavam em desequilíbrio de ligação completa em uma região intrônica; esta região pode ser o alvo de seleção de equilíbrio ou ligada a um site TSP codificado não detectado.

Todas as regiões candidatas nos cinco genes produzem uma árvore alélica, em vez de uma árvore espécie

Embora os haplótipos de cada região agrupados por alelo, em vez de espécies, a partilha de haplótipos entre as duas espécies foi raramente detectada, excepto em AT2G16570 (Col-0 partilhou o seu haplótipo com vários C. rubéola; Fig. 4). Isto não é surpreendente dado o longo tempo de divergência; a partilha extensiva de haplótipos geralmente aparece a uma escala de tempo muito menor e é induzida por eventos como a recente introgressão entre espécies estreitamente relacionadas.

Estudos de simulação neutra validam os cinco genes candidatos

Para ver se as janelas observadas poderiam ser geradas aleatoriamente sob evolução neutra, resultando em falsos positivos, fizemos simulações adicionais baseadas nos parâmetros demográficos estimados usando fastsimcoal2 (arquivo adicional 4: Texto S1). Além das mutações neutras recorrentes, o fluxo gênico também pode resultar em SNPs compartilhados. Assim, fizemos simulações tanto sob o Modelo M1 (sem fluxo gênico) quanto sob o M2 (com fluxo gênico antigo), embora nossa análise demográfica tenha indicado que M1 se encaixa ligeiramente melhor nos dados. Em ambas as simulações, consideramos heterogeneidade nas taxas de mutação para diferentes classes de mutações, notadamente a maior taxa de mutação em locais CpG, o que pode resultar em falsos positivos (arquivo adicional 1: Tabela S5, arquivo adicional 4: Texto S1). Usando fastsimcoal2 , geramos 1.000.000 segmentos neutros de 100 bp sob cada modelo e procuramos por aqueles com dois ou mais SNPs compartilhados e cluster por alelo enquanto pesquisávamos por TSPs.

Para ambos os modelos, nenhuma das 1.000.000 execuções deu origem a uma janela que atendesse aos nossos critérios (Arquivo Adicional 1: Tabela S6). Apesar da existência de SNPs neutros compartilhados, nenhuma janela simulada deu origem a uma árvore alélica, uma vez que todas as janelas com SNPs compartilhados foram acompanhadas de diferenças muito mais fixas entre as duas espécies, implicando em níveis de divergência mais elevados do que a diversidade. Este resultado sugere que estes SNPs neutros simulados compartilhados são mutações recorrentes, ao invés de TSPs, e mais importante, os cinco genes que encontramos não são consistentes com a evolução neutra e, portanto, provaram ser TSPs reais sob seleção de equilíbrio. Os locais e genes finais de TSP estão listados na Tabela 1. Além disso, junto com o estudo demográfico acima mencionado, nossos resultados implicam que mesmo que ocorresse um fluxo de genes antigos, sob evolução neutra, as TSPs se perderiam por deriva neste sistema.

Propriedades dos genes sob seleção de equilíbrio

A seguir calculamos a diversidade de nucleotídeos (π) para todas as regiões de TSP nos cinco genes de cada espécie e usamos as sequências neutras simuladas sob M1 para determinar os níveis de diversidade de fundo. Todas as regiões nos cinco genes exibiram valores significativamente maiores π do que os níveis de fundo tanto em C. rubella quanto em A. thaliana (Wilcoxon-Mann-Whitney test, corrigido por FDR P < 0.05, Tabela 2, arquivo adicional 3: Figura S2A), exceto AT5G38460 em A. thaliana. Além disso, os alelos destes genes mostraram uma tendência para frequências intermediárias (Wilcoxon-Mann-Whitney test, P = 0,0752/0,03474 para A. thaliana/C. rubella; Arquivo adicional 3: Figura S2B). Entretanto, uma freqüência intermediária é uma indicação de seleção de equilíbrio, mas não uma evidência definitiva, uma vez que a distribuição da freqüência alelo de locais ligados a um polimorfismo balanceado é esperado que exiba uma mudança em direção ao equilíbrio de freqüência, que pode estar em qualquer freqüência alelo .

Um dos cinco genes sob seleção de equilíbrio a longo prazo neste estudo, AT1G35220, tem uma função desconhecida, mas exibe fosforilação protéica sob tratamento com etileno . Entre outros, AT2G16570 é uma enzima chave na via de biossíntese de nucleotídeos puros e é importante para a divisão celular, biogênese cloroplástica e germinação de sementes; AT4G29360 é uma proteína da família O-glycosyl hydrolase 17, envolvida em respostas de defesa ; AT5G38460 é uma glicosiltransferase e catalisa a transferência de um grupo glicosil de um composto (doador) para outro (aceitante) e está envolvido em diversas funções, incluindo estresse biótico ; AT5G44000 é uma glutationa S-transferase, que geralmente está envolvida na resposta ao estresse abiótico e biótico . Aparentemente, esses genes estão potencialmente envolvidos na resposta ao estresse biótico ou abiótico (AT4G29360, AT5G38460 e AT5G44000) ou funções bioquímicas fundamentais (AT2G16570).

Como esperado, os genes sob seleção de equilíbrio foram funcionalmente importantes e todos os homólogos dos cinco genes já existiam no ancestral comum mais recente das plantas verdes. Como indicado na Tabela S7 (arquivo adicional 1: Tabela S7), homólogos (ortologues ou paralogues) podem ser encontrados até mesmo nas espécies mais basais de plantas verdes, Chlamydomonas reinhardtii, para todos os cinco genes, exceto AT4G29360, que podem ser rastreados até Physcomitrella patens.

No entanto, loci que são amplamente aceitos como estando sob seleção de equilíbrio, como os genes S-locus ou R , não se destacaram neste estudo. Isto é esperado, já que estes loci são muito variáveis para serem identificados com base em leituras curtas. Por exemplo, os genes R são demasiado dinâmicos para chamar SNPs; o S-locus não existe na última anotação do genoma Arabidopsis e apenas um haplótipo S-locus é mantido em C. rubella desde a transição do outcrossing para o selfing e a quebra da auto-incompatibilidade . Além disso, o S-locus não está mais sob seleção de equilíbrio, uma vez que ambas as espécies são agora autocompatíveis. Em contraste, os genes que identificamos aqui, embora antigos, não foram estudados de forma abrangente e podem fornecer uma visão dos tipos de genes sob seleção de equilíbrio.

A seleção de equilíbrio contribuiu para a adaptação a habitats divergentes

Para ver se as variantes alélicas sob seleção de equilíbrio a longo prazo estão associadas à diversificação ecológica, investigamos a divergência com respeito a 48 fatores ecológicos (arquivo adicional 5: Tabela S8A). Devido à falta de informação GPS e ao pequeno tamanho da amostra de C. rubella, esta análise só foi possível para as amostras de A. thaliana. A estrutura populacional é geralmente altamente correlacionada com a diversificação ecológica e pode, portanto, confundir os nossos resultados. Nós primeiro verificamos se algum site de TSP estava correlacionado com a estrutura da população nas amostras de A. thaliana, embora tal estrutura não afeta a probabilidade de observar a árvore da espécie A. thaliana e C. rubella. Usando ADMIXTURE , verificamos que as 80 amostras de A. thaliana podem ser classificadas em dois grupos (arquivo adicional 3: Figura S3; arquivo adicional 6: Tabela S9) e apenas as classificações alélicas dos dois locais do gene AT5G38460 estão significativamente correlacionadas com a estrutura populacional (teste qui-quadrado, corrigido por FDR P < 0.05,; arquivo adicional 1: Tabela S10). Assim, excluímos AT5G38460 de análises ecológicas subsequentes.

Para obter uma compreensão completa da divergência ecológica, usamos 1135 genomas recentemente publicados A. thaliana . Primeiro, aplicamos um processo de “desbaste” para garantir que cada amostra fosse altamente representativa do seu habitat natural, o que deixou 584 amostras (ver “Métodos”). Em segundo lugar, para cada gene, classificamos as 584 adesões de A. thaliana em dois grupos com base nos haplótipos faseados para os dois sítios TSP (Ficheiro adicional 5: Tabela S8B, C, algumas amostras foram removidas porque não puderam ser faseadas). Em seguida, avaliamos a divergência entre os dois grupos de acedências em relação aos 48 fatores ecológicos para cada um dos quatro genes. Curiosamente, todos estes quatro genes foram associados com a divergência de alguns parâmetros ecológicos específicos. AT1G35220 e AT4G29360, em particular, mostraram divergência significativa com respeito à maioria dos fatores ecológicos relacionados à temperatura (arquivo adicional 5: Tabela S8 A, teste Wilcoxon-Mann-Whitney, corrigido por FDR P < 0,05).

A seguir modelamos os nichos ecológicos para todos os quatro genes. Aparentemente, os dois grupos de amostras para cada gene, como indicado pela estatística I de Warren que mede a similaridade de nichos, exibiram uma identidade de nicho observada significativamente menor do que 100 permutações aleatórias (teste t de uma amostra, corrigido por FDR P < 0,01; Fig. 5a, arquivo adicional 5: Tabela S8 D). Em outras palavras, os dois grupos alélicos de amostras exibem divergência significativa de nicho. Além disso, as amostras de cada tipo alélico para cada gene foram dispersas, ao invés de serem isoladas em uma pequena área local (arquivo adicional 3: Figura S4). Estes resultados sugerem que todos estes loci estão correlacionados com a adaptação.

Divergência Ecológica e de Expressão. a Divergência ecológica significativa entre os dois tipos de amostras para cada um dos quatro genes, conforme indicado pelo escore I observado (I O) e escore I simulado (I S). b Divergência de expressão do gene AT5G44000. c Esquerda: Modelagem do nicho com alta probabilidade (≥ 0,5) dos dois tipos de amostras para AT5G44000. Direita: Resultados significativos sob diferentes estratégias de permutação (para nichos com probabilidade ≥ 0,5; I O = 0,673, 100 permutações)

Examinamos também a diferenciação da expressão dos quatro genes entre os dois grupos correspondentes com base nos haplótipos faseados nos dois sítios da TSP, escolhendo 84 transcriptomas de folhas extraídas de A. thaliana (uma amostra foi sequenciada para cada adesão e o nível de expressão foi medido como fragmentos por kilobase de exon por milhão de fragmentos mapeados ) como nosso estudo anterior . Um gene, AT5G44000, apresentou diferença significativa de expressão (teste Wilcoxon-Mann-Whitney, corrigido por FDR P < 0,05, Fig. 5b) entre os dois grupos haplótipos.

Por isso, realizamos uma modelagem profunda de nicho de AT5G44000 (Fig. 5c) e examinamos a diversificação dos dois grupos de amostras (503 vs 75). Primeiro comparamos a identidade do nicho entre os dois grupos de haplótipos AT5G44000, restringindo a nossa análise a nichos com elevada probabilidade (≥ 0,5) e obtivemos resultados semelhantes (Fig. 5c, Ficheiro adicional 5: Tabela S8 D). Para ver se o tamanho desequilibrado da amostra poderia afetar os resultados, usamos outra estratégia de permutação, restringindo a análise ao mesmo tamanho de amostra (75) para ambos os conjuntos em cada repetição (com probabilidade > 0,5). Conforme apresentado na Fig. 5c, quando a permutação foi realizada para os grupos de amostra reais (simulação 1), o valor observado I (0,673) não apresentou diferença significativa (teste t de uma amostra, P = 0,166), indicando que o valor observado foi confiável, independentemente da diferença de tamanho da amostra. Quando os dois grupos reais foram misturados e dois grupos aleatórios de tamanhos reais foram selecionados (simulação 2) ou dois grupos aleatórios de tamanho igual (75) foram selecionados (simulação 3), a diferença entre o valor observado e as permutações foi significativa novamente (teste t de uma amostra, P = 1,9 × 10-75 para a simulação 2 e P = 2,6 × 10-75 para a simulação 3). Estes resultados implicam que os dois grupos haplótipos funcionalmente diferenciados de AT5G44000 adaptaram-se a habitats ecológicos divergentes.