Amélioration du génome de référence de Brassica rapa par séquençage à molécule unique et technologies de capture de conformation chromosomique

Assemblage du génome

Pour guider l’assemblage du génome, nous avons estimé la taille du génome de B. rapa par cytométrie en flux en utilisant le riz comme référence. Nous avons initialement estimé que B. rapa a une taille de génome de 455 Mb (tableau supplémentaire S1). Un examen plus approfondi impliquant des calculs pour la longueur totale de la carte consensuelle générée sur la base des données BioNano a indiqué une taille du génome de 442,9 Mb (tableau supplémentaire S2). Les deux estimations étaient plus petites que la taille précédemment rapportée de 52915 ou 485 Mb2.

Nous avons assemblé le génome de B. rapa en utilisant une couverture d’environ 57 fois des sous-fils de séquençage PacBio (~25.Nous avons assemblé le génome de B. rapa en utilisant une couverture d’environ 57 fois des sous-lignes de séquençage PacBio (~25. 88 Gb), une couverture d’environ 456 fois des données BioNano (~207.70 Gb) et une couverture d’environ 164 fois des lectures Hi-C (~74.64 Gb).L’assemblage résultant était composé de 1476 contigs, avec un contig N50 de 1.45 Mb et une longueur totale de 351.06 Mb (Tableau 1). Par la suite, nous avons détecté des divergences dans 22 contigs en utilisant les lectures Hi-C (Tableau supplémentaire S4). Au lieu de supprimer ces contigs, nous les avons divisés au niveau des régions de conflit ; les données du contig01464 sont présentées à titre d’exemple (figure supplémentaire S1).

Après échafaudage et estimation de la taille des lacunes à l’aide des cartes BioNano et des lectures de paires de matrices (de BRAD, http://brassicadb.org), nous avons obtenu 1301 échafaudages avec un échafaudage N50 de 4,44 Mb (tableau 1). Pour assigner les échafaudages résultants à leurs positions chromosomiques, nous avons ancré ces échafaudages en utilisant les données Hi-C et la carte génétique améliorée (voir Méthodes). Nous avons ancré 298,19 Mb de séquence sur dix chromosomes, dont 200 échafaudages regroupés par les données Hi-C et 8 échafaudages assignés par la carte génétique. Notre assemblage final, appelé génome B. rapa v3.0, totalisait 353,14 Mb de séquence avec 396 lacunes (2,08 Mb) (Tableau 1). Le génome B. rapa v3.0 est plus long que v1.5 mais plus court que v2.5.

Pour évaluer la qualité du génome B. rapa v3.0, nous avons utilisé diverses sources de données. Tout d’abord, nous avons validé l’exhaustivité de notre assemblage en recherchant les gènes eucaryotes de base (CEGs) en utilisant CEGMA16. Au total, 247 des 248 GEC étaient complets, et 1 GEC était partiel, ce qui indique que tous les GEC ont pu être détectés dans notre assemblage (Tableau supplémentaire S6). Ensuite, la qualité du génome a été testée en faisant correspondre les séquences des étiquettes de séquences exprimées (EST) de B. rapa (téléchargées à partir de dbEST au NCBI), ce qui a montré que 99,34 % des EST pouvaient être trouvées dans le génome B. rapa v3.0 nouvellement assemblé.

Amélioration de la contiguïté

Le génome B. rapa v3.0 a amélioré la contiguïté en termes de lacunes et de tailles de contig. Le génome B. rapa v1.5 a été généré à partir de séquences Illumina, alors que davantage de lectures Illumina et une quantité relativement faible de données de séquences PacBio ont été utilisées pour l’assemblage v2.5. Ces deux assemblages présentent des limites en raison de leur fragmentation et de leur faible contiguïté (tableau 1). En combinant le séquençage unimoléculaire, la cartographie optique et la technologie Hi-C, le génome de B. rapa v3.0 représente une amélioration de la contiguïté d’environ 27 fois (contig N50 : 1446 Kb contre 53 Kb, v2.5) et d’environ 31 fois (contig N50 : 1446 Kb contre 46 Kb, v1.5) par rapport aux deux assemblages précédents (tableau 1). Il n’y avait que 396 lacunes dans la v3.0, y compris des lacunes de taille connue (122 provenant de BioNano et 74 provenant de données de paires de matrices) et inconnue (190 provenant de l’assemblage d’échafaudages Hi-C et 10 provenant de l’assemblage de cartes génétiques). Par rapport aux assemblages précédents, la v3.0 présente une amélioration d’environ 10 fois (5,89 Ko contre 60,59 Ko, v2.5) et d’environ 7 fois (5,89 Ko contre 40,09 Ko, v1.5) de la taille des lacunes par Mo par rapport aux deux assemblages précédents (tableau 1). En termes de nombre de lacunes par Mb, la v3.0 est supérieure aux v2.5 et v1.5, respectivement, avec ~23 fois (1,15 vs 25,98, v2.5) et ~35 fois (1,15 vs 40,09, v1.5) moins de lacunes par Mb (tableau 1).

Pour évaluer la contiguïté et la précision de l’ordonnancement des échafaudages des trois versions du génome de référence de B. rapa génome de référence, nous avons d’abord reconstruit les cartes génétiques basées sur les trois assemblées en utilisant le même ensemble de données de reséquençage d’une population haploïde double (DH) dérivée d’un croisement de deux lignes de chou chinois de tête17. Nous avons ensuite évalué l’emplacement des binmarkers sur les cartes génétiques en les intégrant aux cartes physiques correspondantes. Sur les 892 binmarkers de notre assemblage, 877 binmarkers (98,3 %) ont été cartographiés sur la carte génétique. Notre assemblage était en accord avec la carte génétique pour 801 binmarkers (91,3%), ce qui indique la haute qualité de la v3.0 (Fig.1 ; Tableau supplémentaire S7). Cependant, nous avons remarqué que 76 (8,7%) binmarkers sur les chromosomes A05, A08 et A09 correspondaient à des emplacements ambigus sur la carte génétique. Ces régions contenaient des séquences répétées, en particulier dans les régions centromériques, comme décrit dans l’analyse suivante. Cependant, ces régions conflictuelles ont été couvertes par les lectures PacBio et/ou les cartes BioNano ; les données du chromosome A08 dans la v3.0 sont présentées à titre d’exemple (figure supplémentaire S2).

Il y avait 1092 binmarkers sur la carte génétique de la v2.5 et 866 binmarkers sur la carte génétique de la v1.5. Cependant, nous n’avons pu cartographier que 88,7% des binmarkers (969 sur 1092) et 92,3% des binmarkers (799 sur 866) sur la carte génétique de la v2.5 et de la v1.5, respectivement (Tableau supplémentaire S7). Nous avons constaté que 15,1 % des binmarkers (166 sur 969) de la v2.5 étaient discordants, dont 146 binmarkers présentant des distances génétiques et physiques désordonnées au sein du même chromosome (intrachromosome) et 20 binmarkers présentant des distances génétiques et physiques incohérentes sur différents chromosomes (interchromosome) (figure supplémentaire S3 ; tableau supplémentaire S7). Pour la v1.5, 10,0 % des binmarkers (80 sur 799) étaient discordants, dont 71 binmarkers d’intrachromosome et 9 binmarkers d’interchromosome (Supplementary Figure S4 ; Supplementary Table S7). Cependant, la v3.0 contenait le moins de binmarkers intrachromosomiques conflictuels (8,7 %, 76 binmarkers sur 877) et aucun binmarker interchromosomique discordant (Supplementary Table S7), ce qui indique que le génome de B. rapa v3.0 présente une contiguïté plus élevée que les deux assemblages précédents. Prises ensemble, ces validations indépendantes suggèrent que B. rapa genome v3.0 a la plus grande contiguïté et le meilleur ordonnancement des échafaudages parmi les trois assemblages B. rapa.

Comparaison de l’annotation du génome

Nous avons prédit et annoté les modèles de gènes comme décrit précédemment6. Nous avons identifié un total de 45 985 modèles de gènes codant pour des protéines dans la v3.0, ce qui représentait 14,74 % de l’assemblage du génome (tableau 1). Dans notre assemblage, 98,75 % (45 411 sur 45 985) des gènes étaient annotés sur des chromosomes, et seulement 1,25 % (574 sur 45 985) étaient situés sur des échafaudages. Les gènes annotés de novo dans la v3.0 ont été nommés selon la norme de nomenclature des modèles de gènes pour les génomes de référence de Brassica (http://www.brassica.info/info/genome_annotation.php). Le nombre de modèles de gènes dans le nouvel assemblage est supérieur à celui de la v1.5 (41 020 gènes) mais inférieur à celui de la v2.5 (48 826 gènes) (Tableau 1). Pour évaluer davantage la qualité de l’annotation, une comparaison avec l’annotation des assemblages précédents a été effectuée à l’aide de BUSCO18, qui est basé sur un repère de 1440 gènes végétaux conservés. Environ 97,7 % de ces gènes végétaux conservés ont été identifiés, et 1,7 % ont été détectés comme des fragments présentés dans la v3.0 (tableau supplémentaireS11).

Une analyse de la synténie du génome a été effectuée parmi les trois assemblages en utilisant SynOrths19 pour identifier les paires de gènes synténiques et les réseaux de gènes en tandem. Un total de 2077 tableaux tandem (correspondant à 4963 gènes tandem) a été identifié dans la v3.0. Le même nombre de matrices tandem (2077 matrices correspondant à 5004 gènes) a également été détecté dans la v1.5. Une évaluation de la synténie à l’échelle du génome a indiqué que 1539 matrices en tandem (correspondant à 3757 gènes) dans la v3.0 étaient syntétiques à 1494 matrices en tandem (correspondant à 3670 gènes) dans la v1.5. Cependant, davantage de réseaux en tandem (3535 réseaux, 8002 gènes) ont été identifiés dans la v2.5 (Tableau 1). Nous avons détecté des lacunes dans les régions de gènes tandem superflus dans la v2.5, alors qu’aucune lacune n’a été trouvée dans les versions 3.0 et v1.5 (Fig. 2a). Ces lacunes peuvent être le résultat d’erreurs d’assemblage produites par la fermeture des lacunes en utilisant les lectures PacBio dans la v2.5, ce qui a conduit à l’annotation invalide des gènes tandem. Pour d’autres gènes tandem sans lacunes, nous avons observé que des gènes uniques dans la v3.0 et la v1.5 ont été annotés comme deux gènes ou plus dans la v2.5 (Fig. 2b).

En considérant chaque réseau en tandem comme un locus de gène unique, il restait 43 099 gènes dans la v3.0, 44 359 gènes dans la v2.5 et 38 093 gènes dans la v1.5 (tableau 1). Nous avons ensuite effectué une analyse de la synténie des gènes, qui a révélé que 39 858 gènes (92,48 %) dans la v3.0 servaient d’homologues à 40 442 (91,17 %) et 35 464 gènes (93,10 %) dans les v2.5 et v1.5, respectivement. Après comparaison des gènes annotés avec ceux des premières versions, nous avons identifié 3 241 gènes spécifiques de la version 3.0 par rapport aux versions 2.5 et 1.5. Parmi ceux-ci, 2380 gènes étaient soutenus par des preuves provenant de lectures d’ARNm correspondantes de B. rapa (de BRAD, http://brassicadb.org/), et 2295 gènes étaient soutenus par des séquences de protéines d’autres espèces de Brassicaceae (Tableau supplémentaire S12). Au total, 89,10% (2888 sur 3214) des gènes spécifiques de la version dans la v3.0 ont été soutenus par les données d’ARNm de B. rapa ou les séquences de protéines d’autres espèces de Brassicaceae, tandis que seulement 10,90% (326 sur 3214) des gènes n’étaient pas soutenus.

Un nouvel événement d’expansion LTR-RT identifié dans l’assemblage mis à jour

Nous avons annoté les TE dans la v3.0 en utilisant les mêmes méthodes que celles rapportées précédemment20. Un total de 235 683 TE a été identifié à partir de 1244 familles dans la v3.0, et 562 familles de TE uniques ont été trouvées par rapport aux v2.5 et v1.5. Dans la v3.0, les TEs représentent 37,51% (134 Mb) du génome assemblé, ce qui est plus élevé que dans les assemblages précédents (32,30%, 126 Mb, v2.5 ; 25,44%, 72 Mb, v1.5)2,6. Dans notre nouvel assemblage, les TE les plus abondants sont les LTR-RT, qui couvrent une longueur totale de 57,64 Mb et représentent 16,32 % du génome assemblé. Les répétitions non-LTR-RT (LINEs et SINEs) représentent 3,10% de notre assemblage (Figure supplémentaire S5). Nous avons détecté des transposons d’ADN correspondant à 26,35 Mb, qui représentent 7,46% de l’assemblage du génome (Supplementary Figure S5). Une liste complète des TE et des répétitions identifiés dans la version 3.0 se trouve dans le tableau supplémentaire S13. De plus, nous avons identifié un total de 1231 miRNA, 1281 tRNA, 2865 rRNA et 3737 snRNA dans le génome B. rapa v3.0 (tableau supplémentaire S19).

Dans notre assemblage actuel, nous avons annoté plus de LTR-RT (57 Mb) par rapport à v2.5 (44 Mb) et v1.5 (18 Mb). Nous avons identifié 51 062 LTR-RT non intacts dans la v3.0. Une analyse plus approfondie a révélé que 65,27 % (33 672 sur 51 602) des LTR non intacts étaient situés sur les dix chromosomes, tandis que 34,73 % (17 922 sur 51 602) des LTR-RT non intacts se trouvaient sur les échafaudages non ancrés. En utilisant la même méthode6, un total de 13 318 LTR-RT intacts ont été annotés dans la version 3.0. Cependant, il n’y avait que 4129 et 801 LTR-RT intacts dans la v2.5 et la v1.5, respectivement6. Une analyse plus approfondie a révélé que seulement 18,19 % des LTR-RT intacts (2423 sur 13 318) étaient situés sur les dix chromosomes, alors que la plupart (81,81 %, 10 895 sur 13 318) des LTR-RT intacts se trouvaient sur les échafaudages non ancrés dans la v3.Le temps d’insertion des LTR-RT intacts a été calculé comme décrit précédemment4, ce qui indique que le génome de B. rapa a subi trois vagues d’expansion des LTR-RT depuis sa divergence avec B. oleracea (Fig. 3). Ces LTR-RTs intacts avaient un âge d’insertion moyen de 1,88 millions d’années (MYA), avec un âge d’insertion médian de 1,59 MYA. En outre, nous avons trouvé plus de LTR-RT intacts avec des longueurs différentes dans la v3.0 par rapport à la v2.5 et à la v1.5 (figure supplémentaire S6).

Avec ces LTR-RT intacts, un nouvel événement d’expansion LTR-RT a été identifié dans le génome de B. rapa. Nous avons désigné 3155 événements d’insertion LTR-RT intacts de 0 MYA à 0,4 MYA comme une  » jeune expansion  » avec une longueur moyenne de 8135 bp et une date d’insertion moyenne de 0,2 MYA ; 2283 événements d’insertion LTR-RT intacts de 1,0 MYA à 1,4 MYA comme une  » expansion moyenne  » avec une longueur moyenne de 11 902 bp et une date d’insertion moyenne de 1,2 MYA ; et 1444 événements d’insertion LTR-RT intacts de 3.Les expansions jeunes et anciennes correspondent étroitement aux événements d’expansion identifiés précédemment ; l’expansion moyenne a été identifiée pour la première fois dans le génome de B. rapa et a un temps d’insertion similaire à celui de l’événement d’expansion LTR-RT intact chez B. oleracea. De plus, 1778 LTR-RT de type Ty1/Copia et 4179 LTR-RT de type Ty3/Gypsy ont été identifiés dans la v3.0, ce qui est beaucoup plus que ceux identifiés dans les assemblages précédents (353 Ty1/Copia et 632 Ty3/Gypsy dans la v2.5, 260 Ty1/Copia et 162 Ty3/Gypsy dans la v1.5) (Tableau supplémentaire S20 ; Figures supplémentaires S7, S8). En général, il y avait plus de LTR-RT de type Ty3/Gypsy que de LTR-RT de type Ty1/Copia (tableau supplémentaire S20). Par rapport à la v2.5 et à la v1.5, les LTR-RT de type Ty3/Gypsy de la v3.0 ont manifestement augmenté depuis 5 ans (figure supplémentaire S7), tandis que les LTR-RT de type Ty1/Copia ont augmenté depuis 2,2 ans (figure supplémentaire S8). À partir des arbres phylogénétiques, nous avons constaté que chaque groupe de LTR-RT avait plus de copies dans la v3.0 que dans la v2.5 et la v1.5 (tableau supplémentaire S21, S22 ; figure supplémentaire S9, S10).

Blocs de génome et centromères dans le génome de B. rapa

Nous avons étudié les relations des blocs de génome en utilisant l’assemblage mis à jour v3.0. Pour définir les blocs de génome et les centromères dans le génome de B. rapa v3.0, nous avons d’abord construit les trois sous-génomes (LF, MF1 et MF2) sur la base de la relation syntétique entre v3.0 et A. thaliana (figure supplémentaire S11 ; tableau supplémentaire S14). Nous avons détecté 71 des 72 (3 × 24) blocs génomiques attendus dans la v3.0, et la plupart d’entre eux étaient disposés en ligne avec ceux précédemment rapportés dans la réf.10 (Fig. 4 ; Tableau supplémentaire S15). Dans la v3.0, les deux nouveaux blocs génomiques fragmentés F (LF) et F (MF1) ont été identifiés sur les chromosomes A01 et A05 et n’ont pas été observés dans la réf.10. Nous n’avons pas pu détecter deux blocs génomiques très petits décrits précédemment dans la v3.0, le bloc C (MF2) sur le chromosome A07 et le bloc B (MF1) sur le chromosome A08 (réf.10). Les trois petits blocs génomiques adjacents (S (MF2), T (MF2) et B (MF1)) sur le chromosome A08 de la v3.0 étaient ordonnés S/T/B, alors qu’ils étaient ordonnés T/B/S dans la réf.10..

Nous avons également comparé l’orientation des blocs génomiques dans la v3.0 avec celle de la réf.10. Les blocs génomiques W (MF2) sur le chromosome A02, ainsi que G (LF) et E (LF) sur le chromosome A07, se sont révélés inversés par rapport aux autres blocs issus d’un seul chromosome ACK. Cependant, l’orientation du bloc génomique P (LF) sur le chromosome A09 et de trois blocs de V dans la v3.0 était dans le sens direct, alors qu’ils étaient inversés dans la réf.10. Ces résultats ont été étayés par les cartes génétiques de la v3.0 et de la v1.5, respectivement.

Nous avons déterminé avec précision l’emplacement des centromères de tous les chromosomes de la v3.0. En passant au crible les séquences répétées centromériques déterminées précédemment, y compris les répétitions satellites centromériques CentBr, CRB, TR238 et PCRBr21,22,23, nous avons identifié les signaux pour les 21 paléocentromères dans la v3.0, alors que trois paléocentromères n’ont pas été détectés dans la réf.10 (Fig. 4, Tableau supplémentaire S16). L’analyse des paléocentromères a indiqué que les dix centromères existants de B. rapa ont tous été hérités des 21 paléocentromères. Dans la version 3.0, les centromères des chromosomes A01, A03, A04, A05, A06, A07 et A10 avaient les mêmes blocs de génome associés flanquant les centromères correspondants, comme indiqué dans la réf.10 (Fig. 4). Cependant, le centromère du chromosome A02 était situé entre les blocs génomiques P (MF2) et V (MF1), et le centromère du chromosome A09 était situé entre les blocs génomiques P (LF) et B (LF), alors qu’ils étaient considérés comme des paléocentromères dans la réf. 10 (Fig. 4). Le centromère du chromosome A08 était situé entre les blocs génomiques T (MF2) et B (MF1), et non entre les blocs génomiques C (MF1) et T (MF2), comme indiqué dans la réf. 10. De plus, 1188 gènes ont été détectés dans les régions centromériques dans la v3.0, alors que seulement 740 gènes ont été détectés dans la réf.10 (Tableau supplémentaire S17).

Pour évaluer notre assemblage par rapport aux centromères dans la v3.0, nous avons analysé les caractéristiques de séquence des régions centromériques. Nous avons constaté qu’un nombre significativement plus élevé d’ET et de répétitions spécifiques aux centromères était cartographié dans les régions centromériques que dans les autres parties des chromosomes, et que la densité de gènes et le taux de recombinaison étaient nettement plus faibles dans les régions centromériques annotées dans la v3.0 (Fig. 5). En outre, davantage de répétitions spécifiques aux centromères ont été détectées au niveau des régions centromériques dans la v3.0 par rapport à celles rapportées dans la réf.10 (tableau supplémentaire S17).

.