Improvado genoma de referência da Brassica rapa por seqüenciamento monomolécula e tecnologias de captura de conformação cromossômica

Montagem do genoma

Para orientar a montagem do genoma B. rapa, estimamos o tamanho do genoma B. rapa por citometria de fluxo usando o arroz como referência. Estimamos inicialmente que a B. rapa tem um tamanho de genoma de 455 Mb (Tabela Suplementar S1). Outras investigações envolvendo cálculos para o comprimento total do mapa de consenso gerado com base em dados BioNano indicaram um tamanho de genoma de 442,9 Mb (Tabela Suplementar S2). Ambas as estimativas foram menores do que o tamanho previamente reportado de 52915 ou 485 Mb2.

Nós montamos o genoma B. rapa usando uma cobertura de ~57 vezes as subrelações do seqüenciamento PacBio (~25.88 Gb), ~456 vezes a cobertura dos dados BioNano (~207,70 Gb), e ~164 vezes a cobertura das leituras Hi-C (~74,64 Gb). A montagem resultante consistiu de 1476 contigs, com um contig N50 de 1,45 Mb e um comprimento total de 351,06 Mb (Tabela 1). Posteriormente, detectamos discrepâncias dentro de 22 contigs usando as leituras Hi-C (Tabela Suplementar S4). Ao invés de remover essas contigs, nós as dividimos nas regiões de conflito; os dados para Contig01464 são mostrados como exemplo (Figura Suplementar S1).

Após a montagem de andaimes e a estimativa do tamanho das lacunas usando mapas BioNano e leituras de pares mate (do BRAD, http://brassicadb.org), obtivemos 1301 andaimes com um andaime N50 de 4,44 Mb (Tabela 1). Para atribuir os andaimes resultantes às suas posições cromossómicas, ancorámos estes andaimes usando os dados Hi-C e o mapa genético melhorado (ver Métodos). Ancorámos 298,19 Mb de sequência em dez cromossomas que incluíam 200 andaimes agrupados por dados Hi-C e 8 andaimes atribuídos pelo mapa genético. Nossa montagem final, denominada genoma B. rapa v3.0, totalizou 353,14 Mb de sequência com 396 lacunas (2,08 Mb) (Tabela 1). O genoma B. rapa v3.0 é mais longo que o v1.5 mas mais curto que o v2.5.

Para avaliar a qualidade do genoma B. rapa v3.0, usamos várias fontes de dados. Primeiro, validamos a completude da nossa montagem procurando por genes eucarióticos (CEGs) usando o CEGMA16. Um total de 247 dos 248 CEGs foi completo, e 1 CEG foi parcial, indicando que todos os CEGs puderam ser detectados em nossa montagem (Tabela Suplementar S6). Em seguida, a qualidade do genoma foi testada combinando as sequências de ESTs (express sequence tags) do B. rapa (baixado do dbEST no NCBI), que mostrou que 99,34% dos ESTs podiam ser encontrados no genoma B. rapa v3.0.

A melhoria da contigüidade

O genoma B. rapa v3.0 melhorou a contigüidade em termos de lacunas e tamanhos de contigüidade. O genoma B. rapa v1.5 foi gerado a partir das sequências Illumina, enquanto mais leituras do Illumina e uma quantidade relativamente pequena de dados da sequência PacBio foram utilizados para a montagem v2.5. Estas duas montagens têm limitações devido à sua fragmentação e baixa contiguidade (Tabela 1). Ao combinar sequenciamento monomolécula, mapeamento óptico e tecnologia Hi-C, o genoma B. rapa v3.0 representa uma melhoria de ~27 vezes (contig N50: 1446 Kb vs. 53 Kb, v2.5) e ~31 vezes (contig N50: 1446 Kb vs. 46 Kb, v1.5) na contiguidade em relação às duas montagens anteriores (Tabela 1). Houve apenas 396 lacunas na v3.0, incluindo lacunas de tamanhos conhecidos (122 de BioNano e 74 de dados de pares de pares de mate) e desconhecidos (190 de junção de andaimes Hi-C e 10 de junção de mapas genéticos). Em comparação com as montagens anteriores, a v3.0 tem ~10 vezes (5.89 Kb vs. 60.59 Kb, v2.5) e ~7 vezes (5.89 Kb vs. 40.09 Kb, v1.5) melhoria no tamanho das lacunas por Mb em relação às duas montagens anteriores (Tabela 1). Em termos do número de lacunas por Mb, v3.0 é superior a v2.5 e v1.5, respectivamente, com ~23 vezes (1.15 vs. 25.98, v2.5) e ~35 vezes (1.15 vs. 40.09, v1.5) menos lacunas por Mb (Tabela 1).

Para avaliar a contiguidade e precisão da ordenação de andaimes das três versões do B. rapa, primeiro reconstruímos os mapas genéticos baseados nas três montagens usando o mesmo conjunto de dados de ressequenciamento de uma população haplóide (DH) duplicada derivada de um cruzamento de duas linhas de couve chinesa17. Em seguida, avaliamos a localização dos binmarkers nos mapas genéticos, integrando-os com os mapas físicos correspondentes. Dos 892 binmarkers em nossa montagem, 877 binmarkers (98,3%) foram mapeados no mapa genético. Nossa montagem concordou com o mapa genético de 801 binmarkers (91,3%), indicando a alta qualidade da v3.0 (Fig.1; Tabela Complementar S7). Entretanto, notamos que 76 (8,7%) marcadores de binômeros nos cromossomos A05, A08 e A09 foram mapeados para locais ambíguos no mapa genético. Essas regiões continham sequências repetidas, especialmente em regiões centrômicas, conforme descrito na análise a seguir. Entretanto, estas regiões conflitantes foram cobertas por mapas PacBio e/ou BioNano; os dados para o cromossomo A08 na v3.0 são mostrados como exemplo (Figura Complementar S2).

Existiam 1092 marcadores no mapa genético da v2.5 e 866 marcadores no mapa genético da v1.5. No entanto, só conseguimos mapear 88,7% dos marcadores de lixo (969 de 1092) e 92,3% dos marcadores de lixo (799 de 866) no mapa genético da v2.5 e v1.5, respectivamente (Tabela Suplementar S7). Verificamos que 15,1% dos marcadores (166 de 969) na v2.5 eram discrepantes, incluindo 146 marcadores com distâncias genéticas e físicas desordenadas dentro do mesmo cromossomo (intracromossomo) e 20 marcadores com distâncias genéticas e físicas inconsistentes em diferentes cromossomos (intercromossomo) (Figura Suplementar S3; Tabela Suplementar S7). Para a v1.5, 10,0% dos marcadores (80 de 799) foram discrepantes, incluindo 71 marcadores intracromossômicos e 9 marcadores intercromossômicos (Figura S4 Suplementar; Tabela S7 Suplementar). Entretanto, a v3.0 continha os binmarkers intracromossômicos menos conflitantes (8,7%, 76 de 877 binmarkers) e nenhum binmarkers intercromossômicos discrepantes (Tabela Suplementar S7), indicando que o genoma B. rapa v3.0 tem uma contiguidade maior do que as duas montagens anteriores. Tomadas em conjunto, estas validações independentes sugerem que o genoma B. rapa v3.0 tem a maior contiguidade e a melhor ordenação de andaimes entre as três montagens de B. rapa.

Comparação da anotação do genoma

Prevemos e anotamos os modelos genéticos como descrito anteriormente6. Identificamos um total de 45.985 modelos de genes codificadores de proteínas na v3.0, que representaram 14,74% do conjunto genômico (Tabela 1). Em nossa montagem, 98,75% (45.411 de 45.985) dos genes foram anotados em cromossomos, e apenas 1,25% (574 de 45.985) foi anotado em andaimes. Os genes de novo anotados na v3.0 foram nomeados seguindo o padrão de nomenclatura do modelo de genes para os genomas de referência da Brassica (http://www.brassica.info/info/genome_annotation.php). O número de modelos de genes na montagem do novo é superior ao da v1,5 (41.020 genes) mas inferior ao da v2,5 (48.826 genes) (Tabela 1). Para avaliar melhor a qualidade da anotação, foi feita uma comparação com a anotação de montagens anteriores usando BUSCO18, que se baseia em um benchmark de 1440 genes vegetais conservados. Aproximadamente 97,7% desses genes vegetais conservados foram identificados, e 1,7% foram detectados como fragmentos apresentados na v3.0 (Tabela ComplementarS11).

Uma análise do genoma sintético foi realizada entre as três assembléias usando SynOrths19 para identificar pares de genes sintéticos e matrizes de genes tandem. Um total de 2077 matrizes de genes tandem (correspondentes a 4963 genes tandem) foram identificados na v3.0. O mesmo número de arrays tandem (2077 arrays correspondentes a 5004 genes) também foi detectado na v1.5. Uma avaliação do genoma synteny indicou que 1539 arrays tandem (correspondentes a 3757 genes) em v3.0 eram sintéticos para 1494 arrays tandem (correspondentes a 3670 genes) em v1.5. Entretanto, mais arrays tandem (3535 arrays, 8002 genes) foram identificados na v2.5 (Tabela 1). Detectamos lacunas nas regiões de genes tandem supérfluos na v2.5, enquanto não foram encontradas lacunas na 3.0 ou v1.5 (Fig. 2a). Estas lacunas podem ser o resultado de erros de montagem produzidos pelo fechamento de lacunas usando o PacBio lido na v2.5, o que por sua vez levou à anotação inválida dos genes tandem. Para outros genes tandem sem lacunas, observamos que genes individuais na v3.0 e v1.5 foram anotados como dois ou mais genes na v2.5 (Fig. 2b).

Ao tomar cada matriz tandem como um único gene locus, havia 43.099 genes restantes na v3.0, 44.359 genes na v2.5, e 38.093 genes na v1.5 (Tabela 1). Realizamos então uma análise de síntese de genes, que revelou que 39.858 genes (92,48%) no v3,0 serviram como contrapartidas para 40.442 (91,17%) e 35.464 genes (93,10%) no v2,5 e no v1,5, respectivamente. Após a comparação dos genes anotados com os das versões iniciais, identificamos 3241 genes específicos da versão em v3.0, em comparação tanto com v2.5 quanto com v1.5. Destes, 2380 genes foram suportados por evidências de correspondência de mRNA de B. rapa (de BRAD, http://brassicadb.org/), e 2295 genes foram suportados por sequências proteicas de outras espécies de Brassicaceae (Tabela Suplementar S12). No total, 89,10% (2888 de 3214) dos genes específicos da versão em v3.0 foram suportados pelos dados do mRNA de B. rapa ou sequências proteicas de outras espécies de Brassicaceae, enquanto apenas 10,90% (326 de 3214) dos genes não foram suportados.

Um novo evento de expansão LTR-RT identificado na montagem atualizada

Anotamos TEs na v3.0 usando os mesmos métodos reportados anteriormente20. Um total de 235.683 TEs foram identificadas de 1244 famílias na v3.0, e 562 famílias TE únicas foram encontradas em comparação com a v2.5 e v1.5. Na v3.0, as EC representando 37,51% (134 Mb) do genoma montado, que foi maior do que nas montagens anteriores (32,30%, 126 Mb, v2.5; 25,44%, 72 Mb, v1.5)2,6. Em nossa nova assembléia, as TEs mais abundantes são as LTR-RT, que cobrem um comprimento total de 57,64 Mb e representam 16,32% do genoma montado. As repetidas não LTR-RT (LINEs e SINEs) representam 3,10% da nossa montagem (Figura Complementar S5). Detectamos transposições de DNA correspondentes a 26,35 Mb, que representam 7,46% da montagem do genoma montado (Figura Suplementar S5). Uma lista completa de ET identificadas e repetidas em v3.0 pode ser encontrada na Tabela Suplementar S13. Além disso, nós identificamos um total de 1231 miRNAs, 1281 tRNAs, 2865 rRNAs, e 3737 snRNAs no genoma B. rapa v3.0 (Tabela Suplementar S19).

Na nossa montagem atual, nós anotamos mais LTR-RTs (57 Mb) comparado com a v2.5 (44 Mb) e v1.5 (18 Mb). Identificamos 51.062 LTR-RTs não antigos na v3.0. Análise adicional revelou que 65,27% (33.672 de 51.602) dos LTR-RTs não ancestrais estavam localizados nos dez cromossomos, enquanto 34,73% (17.922 de 51.602) dos LTR-RTs não ancestrais estavam localizados nos andaimes não ancorados. Usando o mesmo método6, um total de 13.318 LTR-RTs intactos foram anotados na v3.0. No entanto, foram encontrados apenas 4129 e 801 LTR-RT intactos na v2.5 e v1.5, respectivamente6. Análises posteriores revelaram que apenas 18,19% dos LTR-RTs intactos (2423 de 13.318) estavam localizados nos dez cromossomos, enquanto a maioria (81,81%, 10.895 de 13.318) dos LTR-RTs intactos foram encontrados nos andaimes não ancorados na v3.0.O tempo de inserção de LTR-RTs intactos foi calculado como descrito anteriormente4, o que indicou que o genoma de B. rapa sofreu três ondas de expansão de LTR-RT por divergir de B. oleracea (Fig. 3). Estas LTR-RT intactas tinham uma idade média de inserção de 1,88 milhões de anos atrás (MYA), com uma idade média de inserção de 1,59 MYA. Além disso, encontramos mais LTR-RTs intactos com diferentes comprimentos na v3.0 comparados com na v2.5 e v1.5 (Figura Complementar S6).

Com estes LTR-RTs intactos, um novo evento de expansão LTR-RT foi identificado no genoma B. rapa. Designamos 3155 eventos de inserção de LTR-RT intactos de 0 MYA a 0.4 MYA como uma “expansão jovem” com um comprimento médio de 8135 bp e uma data de inserção média de 0.2 MYA; 2283 eventos de inserção de LTR-RT intactos de 1.0 MYA a 1.4 MYA como uma “expansão média” com um comprimento médio de 11.902 bp e uma data de inserção média de 1.2 MYA; e 1444 eventos de inserção de LTR-RT intactos de 3.0 MYA a 3,4 MYA como uma “expansão antiga” com comprimento médio de 9823 bp e data de inserção (Fig. 3). As expansões jovem e antiga correspondem estreitamente aos eventos de expansão previamente identificados; a expansão média foi identificada pela primeira vez no genoma B. rapa e tem um tempo de inserção semelhante ao do evento de expansão LTR-RT intacto em B. oleracea. Além disso, 1778 Ty1/Copia-like LTR-RTs e 4179 Ty3/Gypsy-like LTR-RTs foram identificados na v3.0, que é muito mais do que aqueles identificados nas montagens anteriores (353 Ty1/Copia e 632 Ty3/Gypsy na v2.5, 260 Ty1/Copia e 162 Ty3/Gypsy na v1.5) (Tabela suplementar S20; Figura suplementar S7, S8). Em geral, havia mais Ty3/Gypsy-like LTR-RTs do que Ty1/Copia-like LTR-RTs (Tabela Suplementar S20). Em comparação com a v2.5 e v1.5, os LTR-RTs semelhantes aos Ty3/Gypsy na v3.0 foram obviamente aumentados desde 5 MYA (Figura Suplementar S7), enquanto os LTR-RTs semelhantes aos Ty1/Copia foram aumentados desde 2.2 MYA (Figura Suplementar S8). A partir das árvores filogenéticas, descobrimos que cada grupo de LTR-RT tinha mais cópias em v3.0 do que em v2.5 e v1.5 (Tabela Suplementar S21, S22; Figura Suplementar S9, S10).

Blocos e centrômeros do genoma B. rapa

Nós investigamos as relações dos blocos do genoma usando a montagem atualizada v3.0. Para definir os blocos e centrômeros do genoma B. rapa no genoma v3.0, primeiro construímos os três subgêneros (LF, MF1 e MF2) baseados na relação sintética entre v3.0 e A. thaliana (Figura Suplementar S11; Tabela Suplementar S14). Detectamos 71 dos 72 (3 × 24) blocos genômicos esperados na v3.0, e a maioria deles foram dispostos de acordo com os anteriormente relatados na ref.10 (Fig. 4; Tabela Suplementar S15). Na v3.0, os dois novos blocos genómicos fragmentados F (LF) e F (MF1) foram identificados nos cromossomas A01 e A05 e não foram observados na ref.10. Não foi possível detectar dois blocos de genoma previamente descritos, muito pequenos na v3.0, bloco C (MF2) no cromossomo A07 e bloco B (MF1) no cromossomo A08 (ref.10). Entretanto, em nossa montagem, os blocos de genoma N/M (MF1), O/P (LF) e A/C(LF) foram dispostos nos cromossomos A01, A09 e A10, respectivamente, enquanto que eles foram dispostos em lados opostos na ref.10. Os três pequenos blocos de genoma adjacentes(S (MF2), T(MF2) e B(MF1) nos cromossomos A08 da v3.0 foram ordenados S/T/B, enquanto que estes foram dispostos como T/B/S na ref.10.

Comparamos também a orientação dos blocos do genoma na v3.0 com a da ref.10. Os blocos genómicos W (MF2) no cromossoma A02, assim como G (LF) e E (LF) no cromossoma A07, foram encontrados invertidos em relação aos outros blocos que se originaram de um único cromossoma ACK. Entretanto, a orientação do bloco do genoma P (LF) no cromossomo A09 e três blocos de V na v3.0 estavam na direção direta, enquanto estes estavam invertidos na ref.10. Estes resultados foram ainda suportados pelos mapas genéticos da v3.0 e v1.5, respectivamente.

Determinamos com precisão a localização dos centrômeros de todos os cromossomos na v3.0. Através da triagem de sequências de repetição centroméricas previamente determinadas, incluindo as repetições centroméricas por satélite CentBr, CRB, TR238 e PCRBr21,22,23, identificamos os sinais para todas as 21 regiões paleocentroméricas na v3.0, enquanto três regiões paleocentroméricas não foram detectadas na ref.10 (Fig. 4, Tabela Suplementar S16). A análise dos paleocentromeros indicou que os dez centromeros de B. rapa existentes foram todos herdados dos 21 paleocentromeros. Na v3.0, os centrômeros dos cromossomos A01, A03, A04, A05, A06, A07 e A10 tinham os mesmos blocos genômicos associados, flanqueando os centrômeros correspondentes, conforme relatado na ref.10 (Fig. 4). Entretanto, o centrômero no cromossomo A02 estava localizado entre os blocos genômicos P (MF2) e V (MF1), e o centrômero no cromossomo A09 estava situado entre os blocos genômicos P (LF) e B (LF), enquanto que estes foram considerados paleocentrômeros na ref.10 (Fig. 4). O centrômero no cromossomo A08 estava localizado entre os blocos genômicos T (MF2) e B (MF1), ao invés de entre os blocos genômicos C (MF1) e T (MF2), como relatado na ref.10. Além disso, foram detectados 1188 genes dentro das regiões centrômicas na v3.0, enquanto apenas 740 genes foram detectados na ref.10 (Tabela suplementar S17).

Para avaliar nossa montagem em relação aos centrômeros na v3.0, analisamos as características da seqüência das regiões centrômicas. Verificamos que um número significativamente maior de TEs e repetições específicas dos centrômeros foi mapeado para as regiões centrômicas do que para outras partes dos cromossomos, e a densidade gênica e a taxa de recombinação foram marcadamente menores nas regiões centrômicas anotadas na v3.0 (Fig. 5). Além disso, mais repetições centrômicas específicas foram detectadas nas regiões centrômicas na v3.0 em comparação com aquelas relatadas na ref.10 (Tabela Suplementar S17).