Mejora del genoma de referencia de Brassica rapa mediante tecnologías de secuenciación de molécula única y de captura de la conformación del cromosoma

Ensamblaje del genoma

Para guiar el ensamblaje del genoma, estimamos el tamaño del genoma de B. rapa mediante citometría de flujo utilizando el arroz como referencia. Inicialmente estimamos que B. rapa tiene un tamaño de genoma de 455 Mb (Tabla Suplementaria S1). Una investigación posterior que incluía cálculos para la longitud total del mapa de consenso generado en base a los datos de BioNano indicó un tamaño del genoma de 442,9 Mb (Tabla Suplementaria S2). Ambas estimaciones fueron más pequeñas que el tamaño previamente reportado de 52915 o 485 Mb2.

Ensamblamos el genoma de B. rapa usando una cobertura de ~57 veces de subredes de secuenciación PacBio (~25.88 Gb), ~456 veces la cobertura de los datos de BioNano (~207,70 Gb), y ~164 veces la cobertura de las lecturas Hi-C (~74,64 Gb).El ensamblaje resultante consistió en 1476 contigs, con un contig N50 de 1,45 Mb y una longitud total de 351,06 Mb (Tabla 1). Posteriormente, detectamos discrepancias en 22 contigs utilizando las lecturas Hi-C (Tabla Suplementaria S4). En lugar de eliminar estos contigs, los dividimos en las regiones conflictivas; los datos del Contig01464 se muestran como ejemplo (Figura Suplementaria S1).

Tras el andamiaje y la estimación de los tamaños de los huecos utilizando los mapas de BioNano y las lecturas de parejas (de BRAD, http://brassicadb.org), obtuvimos 1301 andamios con un N50 del andamio de 4,44 Mb (Tabla 1). Para asignar los andamios resultantes a sus posiciones cromosómicas, anclamos estos andamios utilizando los datos Hi-C y el mapa genético mejorado (ver Métodos). Anclamos 298,19 Mb de secuencia en diez cromosomas que incluían 200 andamios agrupados por los datos Hi-C y 8 andamios asignados por el mapa genético. Nuestro ensamblaje final, denominado genoma de B. rapa v3.0, tenía un total de 353,14 Mb de secuencia con 396 huecos (2,08 Mb) (Tabla 1). El genoma de B. rapa v3.0 es más largo que v1.5 pero más corto que v2.5.

Para evaluar la calidad del genoma de B. rapa v3.0, utilizamos varias fuentes de datos. En primer lugar, validamos la integridad de nuestro ensamblaje mediante la búsqueda de genes eucariotas centrales (CEG) utilizando CEGMA16. Un total de 247 de 248 CEGs estaban completos, y 1 CEG era parcial, indicando que todos los CEGs podían ser detectados en nuestro ensamblaje (Tabla Suplementaria S6). A continuación, se comprobó la calidad del genoma cotejando las secuencias de las etiquetas de secuencias expresadas (EST) de B. rapa (descargadas de dbEST en el NCBI), lo que demostró que el 99,34% de las EST podían encontrarse en el genoma de B. rapa v3.0 recién ensamblado.

Mejora de la contigüidad

El genoma de B. rapa v3.0 ha mejorado la contigüidad en términos de huecos y tamaños de contig. El genoma de B. rapa v1.5 se generó a partir de secuencias de Illumina, mientras que para el ensamblaje v2.5 se utilizaron más lecturas de Illumina y una cantidad relativamente pequeña de datos de secuencias de PacBio. Estos dos ensamblajes tienen limitaciones debido a su fragmentación y baja contigüidad (Tabla 1). Al combinar la secuenciación de una sola molécula, el mapeo óptico y la tecnología Hi-C, el genoma de B. rapa v3.0 representa una mejora de ~27 veces (contig N50: 1446 Kb frente a 53 Kb, v2.5) y ~31 veces (contig N50: 1446 Kb frente a 46 Kb, v1.5) en la contigüidad respecto a los dos ensamblajes anteriores (Tabla 1). En la v3.0 sólo había 396 huecos, incluyendo huecos conocidos (122 de BioNano y 74 de datos de parejas) y desconocidos (190 de la unión de andamios Hi-C y 10 de la unión de mapas genéticos). En comparación con los montajes anteriores, la v3.0 presenta una mejora de ~10 veces (5,89 Kb frente a 60,59 Kb, v2.5) y ~7 veces (5,89 Kb frente a 40,09 Kb, v1.5) en el tamaño de los huecos por Mb con respecto a los dos montajes anteriores (Tabla 1). En cuanto al número de huecos por Mb, la v3.0 es superior a la v2.5 y a la v1.5, respectivamente, con ~23 veces (1,15 frente a 25,98, v2.5) y ~35 veces (1,15 frente a 40,09, v1.5) menos huecos por Mb (Tabla 1).

Para evaluar la contigüidad y la precisión del ordenamiento de los andamios de las tres versiones del genoma de referencia de B. rapa, primero reconstruimos los mapas genéticos basados en los tres ensamblajes utilizando el mismo conjunto de datos de resecuenciación de una población doblemente haploide (DH) derivada de un cruce de dos líneas de col china17. A continuación, evaluamos la ubicación de los marcadores binarios en los mapas genéticos integrándolos con los mapas físicos correspondientes. De los 892 marcadores binarios de nuestro ensamblaje, 877 marcadores binarios (98,3%) estaban mapeados en el mapa genético. Nuestro ensamblaje coincidió con el mapa genético para 801 binmarkers (91,3%), lo que indica la alta calidad de la v3.0 (Fig.1; Tabla Suplementaria S7). Sin embargo, observamos que 76 (8,7%) marcadores binarios de los cromosomas A05, A08 y A09 correspondían a lugares ambiguos en el mapa genético. Estas regiones contenían secuencias repetidas, especialmente en las regiones centroméricas, como se describe en el siguiente análisis. Sin embargo, estas regiones conflictivas estaban cubiertas por las lecturas de PacBio y/o los mapas de BioNano; los datos del cromosoma A08 en la v3.0 se muestran como ejemplo (Figura Suplementaria S2).

Hubo 1092 marcadores en el mapa genético de v2.5 y 866 marcadores en el mapa genético de v1.5. Sin embargo, sólo pudimos mapear el 88,7% de los marcadores binarios (969 de 1092) y el 92,3% de los marcadores binarios (799 de 866) en el mapa genético de v2.5 y v1.5, respectivamente (Tabla Suplementaria S7). Encontramos que el 15,1% de los binmarkers (166 de 969) en la v2.5 eran discrepantes, incluyendo 146 binmarkers con distancias genéticas y físicas desordenadas dentro del mismo cromosoma (intracromosoma) y 20 binmarkers con distancias genéticas y físicas inconsistentes en diferentes cromosomas (intercromosoma) (Figura Suplementaria S3; Tabla Suplementaria S7). Para la v1.5, el 10,0% de los binmarkers (80 de 799) eran discrepantes, incluyendo 71 binmarkers de intracromosoma y 9 binmarkers de intercromosoma (Figura Suplementaria S4; Tabla Suplementaria S7). Sin embargo, la v3.0 contenía el menor número de marcadores intracromosómicos conflictivos (8,7%, 76 de 877 marcadores) y ningún marcador intercromosómico discrepante (Tabla Suplementaria S7), lo que indica que el genoma de B. rapa v3.0 tiene una mayor contigüidad que los dos ensamblajes anteriores. En conjunto, estas validaciones independientes sugieren que el genoma de B. rapa v3.0 tiene la mayor contigüidad y el mejor ordenamiento de los andamios entre los tres ensamblajes de B. rapa.

Comparación de la anotación del genoma

Predijimos y anotamos los modelos de genes como se describió previamente6. Identificamos un total de 45.985 modelos de genes codificadores de proteínas en la v3.0, que representaban el 14,74% del ensamblaje del genoma (Tabla 1). En nuestro ensamblaje, el 98,75% (45.411 de 45.985) de los genes estaban anotados en cromosomas, y sólo el 1,25% (574 de 45.985) estaba localizado en andamios. Los genes anotados de novo en la v3.0 se nombraron siguiendo el estándar de nomenclatura de modelos de genes para los genomas de referencia de Brassica (http://www.brassica.info/info/genome_annotation.php). El número de modelos de genes en el nuevo ensamblaje es mayor que el de la v1.5 (41.020 genes) pero menor que el de la v2.5 (48.826 genes) (Tabla 1). Para evaluar aún más la calidad de la anotación, se realizó una comparación con la anotación de ensamblajes anteriores utilizando BUSCO18, que se basa en una referencia de 1440 genes vegetales conservados. Aproximadamente el 97,7% de estos genes vegetales conservados fueron identificados, y el 1,7% fueron detectados como fragmentos presentados en la v3.0 (Tabla SuplementariaS11).

Se realizó un análisis de sintenia del genoma entre los tres ensamblajes utilizando SynOrths19 para identificar pares de genes sintéticos y conjuntos de genes en tándem. En la v3.0 se identificaron un total de 2077 matrices en tándem (correspondientes a 4963 genes en tándem). El mismo número de matrices en tándem (2077 matrices correspondientes a 5004 genes) se detectó también en la v1.5. Una evaluación de la sintenia en todo el genoma indicó que 1539 matrices en tándem (correspondientes a 3757 genes) en la v3.0 eran sintéticas a 1494 matrices en tándem (correspondientes a 3670 genes) en la v1.5. Sin embargo, en la v2.5 se identificaron más matrices en tándem (3535 matrices, 8002 genes) (Tabla 1). Detectamos lagunas en las regiones de genes superfluos en tándem en la v2.5, mientras que no se encontraron lagunas ni en la 3.0 ni en la v1.5 (Fig. 2a). Estos huecos pueden ser el resultado de errores de ensamblaje producidos por el cierre de huecos utilizando lecturas de PacBio en la v2.5, lo que a su vez condujo a la anotación inválida de genes en tándem. Para otros genes en tándem sin huecos, observamos que genes individuales en v3.0 y v1.5 fueron anotados como dos o más genes en v2.5 (Fig. 2b).

Cuando se tomó cada matriz en tándem como un único locus génico, quedaron 43.099 genes en la v3.0, 44.359 genes en la v2.5 y 38.093 genes en la v1.5 (Tabla 1). A continuación, realizamos un análisis de sintenia de genes, que reveló que 39.858 genes (92,48%) de la v3.0 servían de contrapartida a 40.442 (91,17%) y 35.464 genes (93,10%) de la v2.5 y la v1.5, respectivamente. Tras comparar los genes anotados con los de las primeras versiones, identificamos 3241 genes específicos de la versión v3.0 en comparación con la v2.5 y la v1.5. De estos, 2380 genes fueron apoyados por evidencia de lecturas de ARNm coincidentes de B. rapa (de BRAD, http://brassicadb.org/), y 2295 genes fueron apoyados por secuencias de proteínas de otras especies de Brassicaceae (Tabla Suplementaria S12). En total, el 89,10% (2888 de 3214) de los genes específicos de la versión en v3.0 fueron apoyados por los datos de ARNm de B. rapa o las secuencias de proteínas de otras especies de Brassicaceae, mientras que sólo el 10,90% (326 de 3214) de los genes no fueron apoyados.

Un nuevo evento de expansión LTR-RT identificado en el ensamblaje actualizado

Anotamos TEs en v3.0 utilizando los mismos métodos que se informaron previamente20. Se identificaron 235.683 TEs de 1244 familias en la v3.0, y se encontraron 562 familias de TEs únicas en comparación con la v2.5 y la v1.5. En la v3.0, los TEs representaron el 37,51% (134 Mb) del genoma ensamblado, lo que fue mayor que en los ensamblajes anteriores (32,30%, 126 Mb, v2.5; 25,44%, 72 Mb, v1.5)2,6. En nuestro nuevo ensamblaje, los TEs más abundantes son los LTR-RT, que cubren una longitud total de 57,64 Mb y representan el 16,32% del genoma ensamblado. Las repeticiones no LTR-RT (LINEs y SINEs) representan el 3,10% de nuestro ensamblaje (Figura Suplementaria S5). Detectamos transposones de ADN correspondientes a 26,35 Mb, que constituyen el 7,46% del ensamblaje del genoma (Figura Suplementaria S5). La lista completa de TEs y repeticiones identificadas en la v3.0 se puede encontrar en la Tabla Suplementaria S13. Además, identificamos un total de 1231 miRNAs, 1281 tRNAs, 2865 rRNAs, y 3737 snRNAs en el genoma de B. rapa v3.0 (Tabla Suplementaria S19).

En nuestro ensamblaje actual, anotamos más LTR-RTs (57 Mb) en comparación con v2.5 (44 Mb) y v1.5 (18 Mb). Identificamos 51.062 LTR-RTs no intactos en la v3.0. Un análisis más detallado reveló que el 65,27% (33.672 de 51.602) de los LTRs no intactos estaban localizados en los diez cromosomas, mientras que el 34,73% (17.922 de 51.602) de los LTR-RTs no intactos se encontraban en los andamios no anclados. Utilizando el mismo método6, se anotaron un total de 13.318 LTR-RT intactas en la v3.0. Sin embargo, sólo había 4129 y 801 LTR-RT intactas en la v2.5 y la v1.5, respectivamente6. Un análisis más detallado reveló que sólo el 18,19% de las LTR-RT intactas (2.423 de 13.318) estaban localizadas en los diez cromosomas, mientras que la mayoría (81,81%, 10.895 de 13.318) de las LTR-RT intactas se encontraban en los andamios no anclados en la v3.El tiempo de inserción de las LTR-RTs intactas se calculó como se describió previamente4 , lo que indicó que el genoma de B. rapa sufrió tres oleadas de expansión de LTR-RTs desde que divergió de B. oleracea (Fig. 3). Estos LTR-RTs intactos tenían una edad media de inserción de 1,88 millones de años (MYA), con una edad media de inserción de 1,59 MYA. Además, encontramos más LTR-RTs intactos con diferentes longitudes en la v3.0 en comparación con la v2.5 y la v1.5 (Figura Suplementaria S6).

Con estos LTR-RTs intactos, se identificó un nuevo evento de expansión de LTR-RT en el genoma de B. rapa. Designamos 3155 eventos de inserción LTR-RT intactos de 0 MYA a 0,4 MYA como una «expansión joven» con una longitud media de 8135 pb y una fecha de inserción media de 0,2 MYA; 2283 eventos de inserción LTR-RT intactos de 1,0 MYA a 1,4 MYA como una «expansión media» con una longitud media de 11.902 pb y una fecha de inserción media de 1,2 MYA; y 1444 eventos de inserción LTR-RT intactos de 3.Las expansiones jóvenes y antiguas se corresponden estrechamente con los eventos de expansión previamente identificados; la expansión media fue identificada por primera vez en el genoma de B. rapa y tiene un tiempo de inserción similar al del evento de expansión LTR-RT intacto en B. oleracea. Además, en la v3.0 se identificaron 1778 LTR-RT similares a Ty1/Copia y 4179 LTR-RT similares a Ty3/Gypsy, lo cual es mucho más que los identificados en los ensamblajes anteriores (353 Ty1/Copia y 632 Ty3/Gypsy en la v2.5, 260 Ty1/Copia y 162 Ty3/Gypsy en la v1.5) (Tabla suplementaria S20; Figura suplementaria S7, S8). En general, había más LTR-RT de tipo Ty3/Gypsy que LTR-RT de tipo Ty1/Copia (Tabla Suplementaria S20). En comparación con la v2.5 y la v1.5, los LTR-RT de tipo Ty3/Gypsy en la v3.0 aumentaron obviamente desde 5 MYA (Figura Suplementaria S7), mientras que los LTR-RT de tipo Ty1/Copia aumentaron desde 2,2 MYA (Figura Suplementaria S8). A partir de los árboles filogenéticos, descubrimos que cada grupo de LTR-RTs tenía más copias en v3.0 que en v2.5 y v1.5 (Tabla Suplementaria S21, S22; Figura Suplementaria S9, S10).

Bloques genómicos y centrómeros en el genoma de B. rapa

Investigamos las relaciones de los bloques genómicos utilizando el ensamblaje actualizado v3.0. Para definir los bloques genómicos y los centrómeros en el genoma de B. rapa v3.0, primero construimos los tres subgenomas (LF, MF1 y MF2) basándonos en la relación sintética entre v3.0 y A. thaliana (Figura suplementaria S11; Tabla suplementaria S14). Detectamos 71 de los 72 (3 × 24) bloques genómicos esperados en la v3.0, y la mayoría de ellos estaban dispuestos en línea con los reportados previamente en la ref.10 (Fig. 4; Tabla Suplementaria S15). En la v3.0, los dos nuevos bloques genómicos fragmentados F (LF) y F (MF1) se identificaron en los cromosomas A01 y A05 y no se observaron en la ref.10. No pudimos detectar dos bloques genómicos muy pequeños descritos anteriormente en la v3.0, el bloque C (MF2) en el cromosoma A07 y el bloque B (MF1) en el cromosoma A08 (ref.10). Sin embargo, en nuestro ensamblaje, los bloques genómicos N/M (MF1), O/P (LF) y A/C(LF) estaban ordenados en los cromosomas A01, A09 y A10, respectivamente, mientras que en la ref.10 estaban ordenados en lados opuestos.Los tres pequeños bloques genómicos adyacentes (S (MF2), T(MF2) y B(MF1)) en el cromosoma A08 de la v3.0 estaban ordenados S/T/B, mientras que en la ref.10 estaban ordenados T/B/S.

También comparamos la orientación de los bloques del genoma en la v3.0 con la de la ref.10. Los bloques del genoma W (MF2) en el cromosoma A02, así como G (LF) y E (LF) en el cromosoma A07, resultaron estar invertidos en relación con los otros bloques que se originaron en un solo cromosoma ACK. Sin embargo, la orientación del bloque genómico P (LF) en el cromosoma A09 y tres bloques de V en la v3.0 estaban en la dirección hacia adelante, mientras que éstos estaban invertidos en la ref.10. Estos resultados fueron respaldados por los mapas genéticos de la v3.0 y la v1.5, respectivamente.

Determinamos con precisión la ubicación de los centrómeros de todos los cromosomas en la v3.0. Mediante el cribado de las secuencias de repetición centromérica previamente determinadas, incluyendo las repeticiones satélite centroméricas CentBr, CRB, TR238 y PCRBr21,22,23, identificamos las señales de las 21 regiones paleocentroméricas en la v3.0, mientras que tres regiones paleocentroméricas no se detectaron en la ref.10 (Fig. 4, Tabla Suplementaria S16). El análisis de los paleocentrómeros indicó que los diez centrómeros existentes de B. rapa fueron todos heredados de los 21 paleocentrómeros. En la v3.0, los centrómeros de los cromosomas A01, A03, A04, A05, A06, A07 y A10 tenían los mismos bloques genómicos asociados flanqueando los correspondientes centrómeros como se informó en la ref.10 (Fig. 4). Sin embargo, el centrómero del cromosoma A02 estaba situado entre los bloques genómicos P (MF2) y V (MF1), y el centrómero del cromosoma A09 estaba situado entre los bloques genómicos P (LF) y B (LF), mientras que estos fueron considerados paleocentrómeros en la ref.10 (Fig. 4). El centrómero del cromosoma A08 estaba situado entre los bloques genómicos T (MF2) y B (MF1), en lugar de entre los bloques genómicos C (MF1) y T (MF2), como se indica en la ref.10. Además, se detectaron 1188 genes dentro de las regiones centroméricas en la v3.0, mientras que en la ref.10 sólo se detectaron 740 genes (Tabla Suplementaria S17).

Para evaluar nuestro ensamblaje con respecto a los centrómeros en la v3.0, analizamos las características de la secuencia de las regiones centroméricas. Encontramos que un número significativamente mayor de TEs y de repeticiones específicas del centrómero fueron mapeadas en las regiones centroméricas que en otras partes de los cromosomas, y la densidad génica y la tasa de recombinación fueron notablemente menores en las regiones centroméricas anotadas en la v3.0 (Fig. 5). Además, se detectaron más repeticiones específicas del centrómero en las regiones centroméricas en la v3.0 en comparación con las reportadas en la ref.10 (Tabla Suplementaria S17).