Séquences génomiques de Tropheus moorii et Petrochromis trewavasae , deux poissons cichlidés éco-morphologiquement divergents endémiques du lac Tanganyika

Assemblées

Sur la base des tailles de génome estimées à ~ 900 Mbp (tableau supplémentaire S11 ), nos efforts de séquençage ont donné des données de séquence avec une couverture de base moyenne de ~ 1,5 ×, ~ 88 ×, ~ 34 × et ~ 10,5 × (PT) et ~ 1,2 ×. , ~ 38×, ~ 29× et ~ 9,1× (TM) pour Roche 454, Illumina PE, Illumina MP et PacBio, respectivement (voir tableau supplémentaire S23 ). Les données de séquence filtrées ont été utilisées pour générer des assemblages primaires dérivés de différents algorithmes de reconstruction (assembleurs) et combinaisons de données (voir Méthodes). Les reconstructions finales du génome des deux espèces sont basées sur des méta-assemblages de ces ensembles d'assemblages primaires. Les méta-assemblages présentant les meilleurs scores en termes de mauvais assemblages, de contiguïté et de prédictions génétiques ont été utilisés dans des analyses ultérieures.

Petrochromis trewavasae

Les assemblages principaux présentent des tailles d'assemblage allant de ~ 779 Mbp à ~ 966 Mbp (907 Mbp PacBio uniquement ; voir le tableau supplémentaire S11 ). L'assemblage final se compose de 7261 échafaudages avec un N50 de 1,84 Mbp, 1,44 % des nucléotides sont indéterminés (N) et 90 % du génome assemblé est contenu dans 885 fragments de plus de 70 kpb. La taille totale de l'assemblage est de 917,57 Mbp (Tableau 1 ).

Tropheus moorii

Les assemblages principaux présentent des tailles d'assemblage allant de ~ 754 Mbp à ~ 952 Mbp (879 Mbp PacBio uniquement ; voir le tableau supplémentaire S11 ). L'assemblage final se compose de 7662 échafaudages avec un N50 de 1,64 Mbp, 1,29 % des nucléotides sont indéterminés (N) et 90 % du génome assemblé est contenu dans 657 fragments de plus de 192 kpb. La taille totale de l'assemblage est de 911,13 Mbp (Tableau 1 ). Les deux tailles d’assemblage se situent dans la plage attendue ; Les prédictions basées sur les spectres k-mer suggèrent des tailles de génome proches de 900 Mbp (voir tableau supplémentaire S11 ) et de 900 à 1 000 Mbp ont également été signalées pour d'autres génomes de cichlidés ^21 , 30 .

Dans ce qui suit, nous comparons nos résultats aux génomes et annotations publiés de plusieurs poissons cichlidés en mettant l'accent sur O. niloticus et M. zebra en raison de leur état bien développé. Les dernières versions (v4) de O. niloticus (44 × PacBio, nouvellement ancré) et de M. zebra (maintenant 65 × PacBio et ancré) ont été publiées par Conte et al . ²⁷ ; la tendance par rapport aux versions antérieures est claire, les qualités des séquences et des annotations sont améliorées et le nombre de structures annotées a encore augmenté. En ce qui concerne les distributions de longueur des gènes (tableau supplémentaire S1 ), les mesures de contiguïté obtenues pour PT et TM sont satisfaisantes et se situent dans la plage typique, compte tenu des technologies de séquençage appliquées et de la couverture (tableau 1 ; pour une comparaison avec les versions d'O. niloticus , voir Tableau supplémentaire S2 , et pour une comparaison générale avec les génomes de poissons publiés, voir le tableau supplémentaire S23 de Vij et al. ³¹ ).

Annotations

Structural annotation yielded ~ 40,300 (PT) and 39,600 (TM) genes and ~ 54,200 (PT) and 56,800 (TM) transcripts, respectively (Table 2); this is in line with the results of different annotation versions of ON (~ 30,200 to 42,600 genes). As to annotated features, PT and TM show similar numbers which often lie between those of version 2 and 3 of the respective ON annotations. For comparison, statistics for ON v2–v4 (the latest) are added, as ON received the most community effort and data for genome assembly and annotation of all cichlids (Supplementary Table S2). Prediction of long non-coding RNAs yielded 2782 and 2112 lncRNAs for PT and TM, respectively. With 57.7% and 63.2% a slight preference for the sense strand could be observed (Supplementary Table S3). Homology based functional annotation could be made for 41,970 (PT) and 43,918 (TM) of the coding sequences (CDSs); putative secretory signals were predicted for 5899 (PT) and 6016 (TM) of them, respectively (Table 3). Pfam domain mapping yielded 78,900 (PT) and 84,158 (TM) hits, respectively. RepeatMasker²⁷ identified 31.1% (PT) and 30.0% (TM) of the genomes as repetitive, respectively; the largest proportions of classified repeat types were held by DNA transposons, LINEs and LTR transposons with ~ 13%, ~ 7% and ~ 2% (Table 4).

Data availability and visualization

The genome and transcriptome assemblies (FASTA), the structural and functional annotations (GFF3), read mappings (BAM) and additional Integrative Genomics Viewer (IGV)³³ track files (short and long non-coding RNAs, repeats, ORFs, CpG islands, microsatellites, IPR and eggNOG domains, variant calls, read mappings, alternative splicing, and REAPR error calls; Fig. 2) are available at https://cichlidgenomes.tugraz.at.

Quality evaluation

Assembly quality was assessed with BUSCO³⁴ and CEGMA³⁵. BUSCO identified 98.3% and 98% of the 4584 proteins in the Actinopterygii database in complete form for PT and TM, respectively; 1.7% and 2% of the benchmarking universal single-copy orthologs (BUSCOs) were either fragmented or missing. These results compare well with those of published genomes and are generally on a par with those of the later versions of the O. niloticus genome drafts (Table 5). CEGMA identified all of the 248 core eukaryotic genes (CEGs) for both PT and TM (Table 6); CEGMA results for PT and TM transcriptome assemblies can be found in Supplementary Table S6. However, REAPR reports 17,166/11,992 (PT/TM) likely assembly errors (Supplementary Table S10); there are IGV tracks highlighting questionable regions to guide caution when analyzing in the vicinity (see Fig. 2). Completeness of conserved protein domains was assessed with DOGMA³⁶. DOGMA found 91.8% and 90.5% of the 1051 expected conserved domains at a conserved domain arrangement size of 1 for PT and TM, respectively (Table 7).

Analyse comparative

Nous avons comparé les génomes de PT et de TM en cartographiant les lectures brutes d'une espèce avec le génome de l'autre espèce. Cela a donné respectivement 4 105 604 et 4 178 777 petites variantes (SMV ; SNP et InDels) pour PT et TM. En outre, 356 428 et 577 124 SMV ont été identifiés pour PT et TM, lors de la cartographie des lectures de la même espèce sur les génomes respectifs. En moyenne, 1 variante pour ~ 220 pb (interspécifique) et 1 variante pour 2 540 pb (PT vs PT)/1 561 pb (TM vs TM) (intraespèce) a été appelée (Tableau 8 ) . Pour les deux espèces, 93 842 et 89 489 grandes variantes structurelles (SV ; insertions, délétions, duplications, inversions et translocations) entre espèces ont été détectées, la majorité étant des délétions avec respectivement 60 % et 65,6 % (Tableau 8 ) .

La distribution du SMV et du SV observée par l'analyse comparative suit largement la couverture génomique de régions structurelles/fonctionnelles particulières. Il existe des dévations petites mais notables : (1) les SMV sont (légèrement) sous-représentés dans les régions du promoteur, de la 5′ UTR, du codage, du site d'épissage, de la 3′ UTR et des régions intergéniques ; ils sont surreprésentés dans les introns ; (2) Les SV sont (légèrement) sous-représentées dans les régions du promoteur, de la 5′ UTR, du codage, de la 3′ UTR et des régions intergéniques ; ils sont surreprésentés dans les introns et les sites d'épissage (via chevauchement) (Fig.  3 A).

SNPeff ³⁷ classe les effets des variantes sur le gène en quatre groupes en fonction de l'emplacement et de la nature de la variante : « ÉLEVÉ », « MODÉRÉ », « FAIBLE » ou « MODIFICATEUR », ce dernier désignant des variantes non codantes ou des variantes affectant des variantes non codantes. -les gènes codants, pour lesquels les prédictions sont difficiles ou où il n'y a aucune preuve d'impact. Dans notre analyse, plus de 97 % des variantes identifiées sont classées comme « MODIFIER » (Tableau 9 ).

Genes of interest, which are possibly affected by mutations, are highlighted in Table 10 (see Supplementary Information for details on the gene selection and GO enrichment analysis, respectively). Here, we started to analyze genes which are related to the development of the viscerocranium (untargeted gene selection) and the pharyngeal system (targeted, i.e., biased gene selection based on literature—see Table S17); however, there are other GO terms of interest which are consistently enriched over different analysis approaches such as BMP signaling, for instance. A condensed GO analysis result for an untargeted approach (A2, see Table S14b) is shown in Table 11; here, gene categories are based on variant comparison groups (within and between species groups) combined with quantile ranking and thresholding (p = 0.5, i.e., median), and variant counts (‘mutation loads’) were used as criterion. The term ‘embryonic viscerocranium morphogenesis’ is enriched in the within and the between species gene sets over all approaches (see Supplementary Table S16; genes belonging to this term were combined with genes from the targeted approach and used for further downstream analyses (see Supplementary Table S14a, Table S14b, Table S15, Table S16, Table S18 and Table S21). In the comparative analysis, biological species are coded as A (PT) and B (TM) (Table 11). The categories (AA, AB, BA and BB) refer to the within and between group comparisons. That is, there are mutations at the same genomic locations (nucleotides) which are either identical within and between species (referred to as, e.g., identical (AA) and redundantly identical (AB)) or nonidentical (referred to as, e.g., nonidentical (BB) and nonidentical (BA)); moreover, there are mutations which are unique to a group (referred to as, e.g., unique (AA) and unique (AB), i.e., at the genomic location there is only a variant in species A (unique (AA)) or there is only a variant between species A and B (unique (AB)), respectively). In the shown example, for SMV the calls for viscerocranium morphogenesis are symmetric except for the AB category (which fell below the threshold), i.e., the GO term is consistently enriched within and between species. Further analyses on the genes belonging to the term clearly verify the presence of shared and species-specific mutations in these genes (see example in Supplementary section Identification of genes putatively related to facial and jaw morphology). Hence, there is substantial variation in these genes which may drive changes in the manifestation of morphology. However, we cannot yet delineate possible effects from shared and non-shared variants.

Besides variants in the DNA structure, alternative splicing (AS) was analyzed. There are ~ 6200 AS events in ~ 2600 genes between sexes of each species and ~ 39,000 AS events in ~ 9400 genes between the two species (see Supplementary Table S13).

Assembly and annotation

Le méta-assemblage d'un ensemble d'assemblages primaires a donné des ébauches de génome de haute qualité avec 918 et 911 Mbp pour Petrochromis trewavasae et Tropheus moori , respectivement. Ceci est conforme aux tailles comprises entre 900 et 1 000 Mbp signalées pour d'autres génomes de cichlidés ^21 , 30 et aux ~ 940 Mbp estimés par notre validation d'assemblage avec REAPR ³⁸ . Le dernier assemblage d'Oreochromis niloticus s'étend sur environ 1 Gbp ; cette variation de la taille du génome peut être due à des différences biologiques, mais pourrait également indiquer que certaines parties de l'ADN répétitif dans la reconstruction génomique respective (PT/TM) ne sont pas incluses dans l'assemblage, car une conséquence de l'effondrement de la séquence (par exemple, répétitions effondrées) ³⁹ .

En règle générale, les assemblages sont basés sur les données génomiques d'un seul individu, qui provient idéalement d'une lignée consanguine. Dans ce projet, nous avons assemblé respectivement 6 (PT) et 5 (TM) individus non consanguins ; cela nécessitait une approche d’assemblage plus complexe. De plus, nous avons effectué un séquençage de novo sans aucune donnée de liaison ou de carte optique (comme le montrent les dernières versions du génome de O. niloticus et A. calliptera ) et la couverture PacBio était (avec ~ 9–10 ×) considérablement inférieure à celle utilisée. pour les assemblages de O. niloticus (44 × pour v3 ³⁰ et v4 ²⁷ ) et M. zebra (16,5 × pour v3 ⁴⁰ , ajoutés en plus des couvertures Illumina PE et MP déjà élevées, et 65 × pour v4 ²⁷ ). Néanmoins, les deux assemblées se comparent bien aux projets publiés sur le génome d’espèces comparables en ce qui concerne les mesures typiques concernant le contenu génétique. Les résultats de BUSCO (Tableau 5 ) montrent un faible taux d'environ 4 à 8 % (selon la base de données) de BUSCO dupliqués en PT et TM. Ceci est légèrement supérieur aux ~ 2 à 4 % rapportés pour d'autres génomes de cichlidés, ce qui peut être une conséquence d'un assemblage incorrect d'haplotypes provenant d'individus non consanguins. En ce qui concerne le nombre total de BUSCO identifiés, fragmentés et manquants, la reconstruction des génomes PT et TM fonctionne très bien.

Pour les deux espèces, les annotations se sont avérées valables pour les premières analyses sensées en aval, mais il existe certainement certains modèles génétiques qui pourraient nécessiter des améliorations supplémentaires, par exemple par un entraînement répété des prédicteurs génétiques ; cependant, pour AUGUSTUS, le prédicteur central, les évaluations du modèle montrent déjà de bons états d'entraînement (voir le tableau supplémentaire S12 ). Les sources les plus pertinentes de modèles génétiques insuffisants pourraient inclure les fusions, les scissions et surtout les troncatures de gènes, qui sont évidentes après un examen plus approfondi – ceci est typique des premières annotations, en particulier lorsque le pipeline d’annotations est encore en cours de développement. Nous observons une longueur moyenne et médiane relativement faible des séquences protéiques (voir Tableau 2 ) dans les deux assemblages/annotations. Cela peut refléter une erreur systématique dans le processus de génération, par exemple des InDels conduisant à des décalages de trame et, par conséquent, à des traductions erronées et à des codons d'arrêt prématurés. L'enquête sur ce phénomène a montré des InDels non triples ; cependant, ceux-ci se retrouvent également entre, par exemple, les modèles de transcription d'O. niloticus et de M.  zebra . De plus, le taux de mutations non-sens identifiées dans PT et TM est faible (Tableau 9 ). Les pipelines d'annotation NCBI et Ensembl sont à la pointe de la technologie ; de plus, la quantité et la diversité des données de séquençage d'ARN utilisées pour l'annotation, par exemple, d'O. niloticus Ã©taient beaucoup plus importantes que ce n'était le cas pour l'une ou l'autre des espèces de ce projet. Par conséquent, le plus grand nombre d’isoformes de transcription identifiées (ainsi que le nombre moyen plus élevé d’exons par transcription) peuvent être considérés comme des conséquences simples. Cependant, le nombre total d’exons chez les deux espèces est comparable aux annotations ON. Il est intéressant de noter que le nombre de modèles génétiques dans PT et TM est également comparable à celui de ON v3. Comme il n'existe pas de méthode bien établie pour évaluer l' exactitude des modèles génétiques (peut-être par une vérification de la structure générale et une notation majoritaire de similarité basée sur une base de données), il s'agit simplement d'une comparaison du nombre d'éléments. De plus, comme mentionné, il existe certaines fusions et scissions de gènes dans les ensembles de modèles de gènes PT et TM, ce qui faussera dans une certaine mesure le décompte des gènes. Comme autre mesure de qualité pour les gènes codant pour les protéines annotés, DOGMA ³⁶ et PfamScan ⁴¹ ont été utilisés ; les résultats soutiennent la notion de mauvais modèles de gènes dans l'ensemble, qui ne contiennent pas certains domaines protéiques ou seulement des fragments de ceux-ci (Tableau 7 ).

Analyse comparative

Nous avons choisi les deux espèces étudiées pour les raisons suivantes. Tropheus moorii est un brouteur d'algues très efficace que l'on trouve en grand nombre sur tous les types de rivages rocheux, tandis que Petrochromis trewavasae est un brouteur d'algues répandu sur les rivages rocheux du côté ouest du lac, vivant en sympatrie avec Tropheus . Les Tropheini comprennent 3 espèces prédatrices, une omnivore, 10 brouteurs d'algues et 15 brouteurs d'algues. Les brouteurs d'algues ont des dents en forme de ciseau pour mordre les algues filamenteuses du substrat rocheux, tandis que les brouteurs d'algues ont des dents en forme de peigne sur plusieurs rangées pour éliminer les algues unicellulaires et les détritus des roches. En raison de l'âge avancé de la tribu Tropheini, qui s'élève à environ 2 à 6,5 millions d'années au début de leur rayonnement ⁶ , le degré de divergence éco-morphologique est plus grand que chez les équivalents éco-morphologiques beaucoup plus jeunes du lac Victoria, mais comparable à l'éco-morphoespace couvert par l'ensemble du troupeau du lac Malawi. Il est intéressant de noter que le genre Tropheus comprend environ 120 populations, pour la plupart allopatriques et distinctes en termes de couleur, ainsi que des espèces sÅ“urs morphologiquement similaires. Ils sont tous restés dans la même niche trophique sur toutes les rives rocheuses du lac. Petrochtomis trewavasae ne présente pas beaucoup de variation de couleur, a une répartition restreinte sur la rive sud-ouest du lac et fait partie d'une lignée de brouteurs complexe et morphologiquement distincte comprenant le complexe d'espèces beaucoup plus diversifié de P. polyodon . Si l’on considère l’ensemble de la lignée, il a subi une trajectoire évolutive similaire à celle de Tropheus . Il convient de noter ici que le nombre d'espèces généralement beaucoup plus faible dans le lac Tanganyika par rapport aux lacs Malawi et Victoria résulte également des différents concepts d'espèces utilisés, dans la mesure où plusieurs entités allopatriques sont traitées comme des espèces dans les lacs Victoria et Malawi, alors que comme des variétés géographiques. dans le rayonnement plus ancien du lac Tanganyika.

L'analyse comparative présentée ici a donné, comme prévu, un grand nombre de régions variantes entre les deux espèces et même un nombre considérable au sein de chaque espèce. La grande quantité de variation au niveau intraspécifique peut en fait être due à notre approche consistant à utiliser plusieurs individus F1 non consanguins d'une même population échantillonnée dans l'environnement naturel, mais reflète mieux la diversité intra-population et, finalement, l'âge évolutif ancien de l'espèce. lignée. Nous avons utilisé GATK ⁴² et DELLY ⁴³ , deux outils bien établis, pour l'appel de variantes ; cependant, l'appel de variantes n'est toujours pas un problème bien résolu avec souvent peu de chevauchement entre les résultats des différentes routes algorithmiques (par exemple, voir ^44 , 45 ). En ce qui concerne les statistiques rapportées sur les effets variants, il est connu que l'état de l'annotation structurelle et l'annotateur d'effet variant utilisé influencent fortement les résultats ⁴⁶ . Les résultats d'analyse présentés ici reflètent l'état des reconstructions du génome (v1).

Le nombre relativement important de SV et de SMV triés réciproquement parmi les deux espèces étudiées est remarquable et pourrait refléter le temps de divergence relativement ancien entre les deux espèces étudiées, s'élevant à environ 2,5 à 6 Mya pour les deux clades ⁶ . En fait, on s’attend à ce que les mutations structurelles affectant les informations codantes mettent plus de temps à évoluer que les mutations régulatrices. Ainsi, lorsque l’on compare des espèces provenant des lacs Victoria et du Malawi, beaucoup plus jeunes, on ne s’attendrait pas à un degré aussi marqué de variation de codage réciproquement distincte. Le SV et le SMV peuvent également être interprétés à la lumière de l'hypothèse de la tige flexible ^4 , 18 . La tige flexible des radiations cichlidés est formée d’espèces écologiquement et phénotypiquement flexibles adaptées aux habitats fluviaux saisonnièrement instables. Une fois qu'ils ont semé des radiations lacustres, ils peuvent rapidement s'adapter à des niches vides dans cet environnement plus stable en raison de leur grande plasticité phénotypique ¹⁸ . Par la suite, la population phénotypiquement plastique est subdivisée en phénotypes adaptatifs alternatifs, puis les facteurs génétiques adaptatifs sont triés au cours de la spéciation pour continuer via l'accommodation génétique et l'assimilation génétique ⁴⁷ . La plasticité phénotypique ou développementale fait référence à la capacité d'un seul génotype à produire plusieurs phénotypes dans différentes conditions environnementales. L'hypothèse de la tige flexible postule que la plasticité d'une population peut influencer la direction de l'évolution en exposant une variation génétique cryptique à la sélection dans un nouvel environnement. Selon ce modèle, des sous-ensembles d'une population ancestrale exploitent des niches écologiques distinctes dans un nouvel habitat, comme différents types d'aliments. Au cours d'une seule génération, la plasticité de l'anatomie peut conduire à une amélioration de la condition physique, par exemple une capture ou une transformation plus efficace des aliments, dans chaque niche. Les variations phénotypiques nouvellement exposées seront ciblées par la sélection, et si le nouvel environnement est stable, les phénotypes plastiques pourront être canalisés par assimilation génétique. L'hypothèse est que les mécanismes moléculaires de la réponse plastique sous-tendent également l'évolution des phénotypes clés, c'est-à-dire que la variation génétique des mêmes molécules/voies de signalisation, qui permettent la plasticité, est ciblée par sélection et fixée afin de canaliser le phénotype. Dans une étude récente, le rôle de la signalisation hérisson (Hh) dans la plasticité cranio-faciale chez les téléostéens a été mis en évidence, démontrant que les niveaux de Hh ajustent la sensibilité aux signaux mécaniques liés aux conditions d'alimentation - où les changements morphologiques adaptatifs dans les structures immédiatement affectées, par exemple, le les os pharyngés, peuvent propager des changements morphologiques à d'autres structures cranio-faciales ⁴⁸ .

Des variantes ont été appelées dans pratiquement toutes les régions génétiques. Environ 99 % ont au moins une variante possible, selon les paramètres appliqués, dans le corps du gène ou 5 kb en amont/en aval (Fig.  3 B). Les gènes avec au moins une mutation ont été soumis à une analyse d'ontologie génétique (GO) pour obtenir des indications sur d'éventuels groupes fonctionnels intéressants affectés par davantage de variantes - c'est-à-dire que le nombre de variantes (ou « charge de mutation ») a été utilisé comme indicateur de la probabilité d'une mutation efficace. changements. La justification de cette approche était l'hypothèse de l'exactitude du modèle infinitésimal ou du modèle omnigénique ⁴⁹ , respectivement. On peut s'attendre à ce que les changements phénotypiques observés ne soient pas dus à quelques variantes à fort impact (généralement une région codante), mais plutôt à plusieurs variantes à « impact moindre » (dans les catégories utilisées, probablement les « variantes modificatrices » qui représentent généralement > 90 % des charge de mutation). Même si à ce stade, la pertinence de la variation dans les gènes sélectionnés n'est pas claire, tous les gènes répertoriés ont de multiples appels concernant le SMV et le SV (Fig. 3  B ), ce qui peut augmenter les chances d'influencer efficacement les phénotypes. Compte tenu des fonctions qui leur sont attribuées et rapportées dans d'autres organismes (tableau 10 ), ces gènes méritent cependant d'être étudiés. Par exemple, cinq gènes liés à la définition de la forme du nez et du menton ( DCHS2 , RUNX2 , GLI3 , PAX1 et EDAR ) ont récemment été identifiés dans une étude GWAS humaine ⁵⁰ ; plusieurs variantes de tous ces gènes ont également été trouvées entre les deux espèces. De plus, les membres des familles PAX3 , KCTD15 et TBX ( TBX1 et TBX10 , mais pas TBX15 comme indiqué précédemment) sont dans le jeu de résultats ; ces gènes ont été liés à la morphologie du visage chez l'homme dans deux autres études GWAS récentes ^51 , 52 (Tableau 10 ). Les futures analyses en aval devraient se concentrer particulièrement sur les gènes présentant des différences stables dans l’expression des gènes entre les espèces étudiées. Comme indiqué précédemment, nous nous concentrons sur les différences de formes faciales et pharyngées (voir Fig. S1 supplémentaire).). Il est intéressant de noter que cette méthode simple de comptage de variantes impartiales («charge de mutation») produit de manière reproductible les termes GO liés à la morphogenèse du viscérocrâne (voir Informations supplémentaires), sans donner une longue liste plutôt non spécifique de termes GO. Du résultat GO découle la mise en évidence de plusieurs voies de signalisation importantes : signalisation BMP (par exemple, bmp2, bmp4), signalisation Hedgehog (Hh) (par exemple, Shh, famille Gli, famille Sec, smo, med12, plcb3), signalisation de l'endothéline (par exemple , edn1, furine, famille dlx), signalisation de l'acide rétinoïque (RA) (par exemple, rere, rerea) et signalisation du facteur de croissance des fibroblastes (FGF) (par exemple, fgf8, fgf20b) (voir le tableau 10 et le tableau supplémentaire S21 ) . Tous ces réseaux de signalisation sont connus pour jouer un rôle dans la régulation de la morphogenèse faciale des vertébrés et interagissent. Il existe, par exemple, de fortes interactions coopératives et fonctionnelles entre Shh et l'acide rétinoïque ^{53 , 54 , 55 , 56 , 57 , 58} . Une analyse comparative plus approfondie de la distribution des variantes génétiques observées entre les deux espèces et de leurs phénotypes respectifs n'a pas été réalisée à ce stade ; ce sera une tâche importante pour les études de suivi.

Pour résumer, les deux nouveaux projets de génomes ajoutent deux espèces clés monophylétiques et éco-morphologiquement divergentes qui comblent une lacune phylogénétique importante. De plus, ils représentent la première ramification des cichlidés haplochromines dits modernes, la lignée la plus riche en espèces de cichlidés d’Afrique de l’Est. Tandis que les Tropheini rayonnaient dans les confins du lac Tanganyika, leurs alliés se répandaient sur plusieurs rivières pour semer des radiations supplémentaires comme celles des lacs Malawi et Victoria, où celles-ci atteignaient une diversité éco-morphologique comparable.

Espèces étudiées

Les spécimens échantillonnés de T. moorii sont des descendants F2 d'individus capturés dans la nature dans la partie zambienne de la rive sud-ouest du lac Tanganyika (08°38′ S 30°52′ E) près du village de Nakaku, qui ont été amenés à l'Université de Graz. en 2005. Les spécimens de P. trewavasae utilisés dans cette étude sont des descendants F1 de poissons sauvages également de la côte sud-ouest, mais plus au nord-est près du village de Katete (08°20′S 30°30′E) et ont été obtenus à partir d'un poisson d'ornement. importateur. La collecte de la génération parentale de poissons a été réalisée dans le cadre d'un protocole d'accord entre le Département des pêches, le Ministère de l'agriculture et des coopératives de Zambie, le Département des sciences biologiques de l'Université de Zambie à Lusaka, le Département de zoologie de l'Université de Graz, en Autriche, le Département d'écologie comportementale de l'Université de Berne, en Suisse, et le Département de zoologie de l'Université de Bâle, en Suisse, en vertu du permis de recherche délivré au CSt par le ministère zambien de l'Intérieur (numéro de permis :SP006515). Les données de séquence présentées ici sont basées sur des extractions d'ADN de 6 P. trewawasae et 5 T . individus moorii ; les spécimens comprenaient les deux sexes et étaient âgés d'environ un an.

Procédures de séquençage et de laboratoire

Nous avons séquencé l'ADN génomique extrait des échantillons ci-dessus à l'aide de plusieurs technologies de séquençage : Illumina HiSeq paire d'extrémités 2 × 101 pb (taille de fragment de 300 pb et 600 pb), Illumina Nextera mate-pair 2 × 100 pb (taille de fragment de 1 à 6 kpb). ), 454 Life Sciences (~ 350 pb de longueur de lecture moyenne ; taille de fragment de 8 et 20 kbps) et technologie de séquençage de molécule unique en temps réel (SMRT) de Pacific Biosciences (PacBio) (~ 8 000 à 9 000 pb de longueur de lecture moyenne après correction) .

Les méthodes liées au laboratoire (extraction d'ADN, préparation de bibliothèques et séquençage) ont, en partie, été décrites précédemment dans l'article d'accompagnement sur les génomes mitochondriaux ⁵⁹ . De plus, nous avons effectué deux cycles de séquençage en utilisant la technologie de séquençage de deuxième génération de Pacific Biosciences, basée sur un individu par espèce. L'extraction de l'ADN a été réalisée à Graz, la préparation de la bibliothèque et le séquençage au Centre de technologies génomiques de Lausanne : l'ADN a été cisaillé dans un g-TUBE Covaris (Covaris, Woburn, MA, USA) pour obtenir des fragments de 20 kpb. Après cisaillement, la distribution de la taille de l'ADN a été vérifiée sur un analyseur de fragments (Advanced Analytical Technologies, Ames, IA, USA). 5 µg d'ADN cisaillé ont été utilisés pour préparer une bibliothèque SMRTbell avec le kit de préparation de modèles PacBio SMRTbell 1 (Pacific Biosciences, Menlo Park, Californie, États-Unis) conformément aux recommandations du fabricant. La bibliothèque résultante a été sélectionnée en taille sur un système BluePippin (Sage Science, Inc. ; Beverly, MA, USA) pour les molécules de plus de 11 kpb. La bibliothèque récupérée a été séquencée sur treize/seize cellules SMRT (TM/PT) avec la chimie P6/C4 et MagBeads sur un système PacBio RSII (Pacific Biosciences, Menlo Park, CA, USA) pendant une durée de film de 240 minutes.

Pour l'ARN-seq, l'ARN total d'un individu mâle et femelle par espèce (regroupé à partir des tissus suivants : foie, rate, cerveau, cÅ“ur et muscle squelettique) a été extrait avec Trizol comme suit : les tissus ont été homogénéisés avec MagnaLyser et incubés avec Trizol- tube 5 min à température ambiante (RT) ; 200 µl de chloroforme (/ml de Trizol) ont été ajoutés et secoués vigoureusement pendant 15 s, incubés pendant 2 à 3 min/RT et centrifugés à 12 200 tr/min/4 °C/15 min ; le surnageant a été transféré dans un nouveau tube de 1,5 ml et 500 µl d'isopropanol (/ml de Trizol) ont été ajoutés ; après vortex, incubation pendant 10 min/RT, centrifugation à 12 200 tr/min/4 °C/10 min, le surnageant a été jeté et le culot placé immédiatement sur la glace. Les culots ont été lavés 2 fois : ajouter 1 ml d'EtOH à 80 % (-20 °C), centrifuger : pleine vitesse/4 °C/5 min, éliminer le surnageant et enfin séchés à 37 °C. Les culots séchés ont été remis en suspension dans 20 µl d'eau distillée. Les bibliothèques d'ARN-seq ont été dérivées d'ARN total appauvri en ARNr, normalisé et séquencé sur une seule voie Illumina HiSeq 2500 par espèce.

(Pré)traitement général des données

Tout le pipeline et le traitement de niveau supérieur ont été effectués avec R/Bioconductor, quelques pipelines mineurs dans Bash et certaines fonctionnalités de pointe ont été écrites en C ++ (appelées depuis R). Pour plus de détails sur les réglages des paramètres pour les étapes/outils importants, voir le tableau supplémentaire S22 .

FastQC v0.10.1 ⁶⁰ a été utilisé pour l'évaluation de base de la qualité de lecture. Une approche personnalisée basée sur le spectre k-mer utilisant JELLYFISH v2.0 ⁶¹ (en conjonction avec une base de données de séquences techniques connues) et une approche basée sur De Bruijn (implémentée dans Minion à partir du package Kraken v13-274 ⁶² ) ont été utilisées pour le identification automatique des contaminants techniques et des séquences suspectes (en fonction des fréquences attendues). De plus, FastQScreen v0.4.4 ⁶³ a été utilisé pour l'identification spécifique à l'espèce de la contamination biologique et DeconSeq v0.4.3 ⁶⁴ pour sa suppression. Cutadapt v1.5 ⁶⁵ a été utilisé pour l'élimination des contaminants techniques, Scythe v0.994 ⁶⁶ pour un découpage supplémentaire de l'adaptateur 3', CLC Quality Trim v4.2 ⁶⁷ pour un découpage de lecture basé sur le score de qualité et Reaper v13-274 ⁶² pour une meilleure qualité. et un filtrage basé sur la complexité. BBmerge v33.40 ⁶⁸ a été utilisé pour la fusion de lectures appariées superposées et FastUniq v1.1 ⁶⁹ pour la suppression des doublons. Nextclip v1.2 ⁷⁰ a été utilisé pour le filtrage et la classification de lecture des paires de partenaires Nextera. 454 ensembles de données ont également été filtrés avec sffToCA (utilitaire Celera Assembler). BAMtools v2.4.0 ⁷¹ , SAMtools/BCFtools/HTSlib v1.4 ⁷² et les outils Picard v1.119 ⁷³ ont été utilisés pour le mappage et les manipulations de fichiers de séquence telles que l'indexation, la fusion, le tri et la génération de sous-ensembles, la suppression des lectures en double et la suppression. de contamination PE à partir de bibliothèques MP dans des fichiers de séquence. Proovread v2.13.10 ⁷⁴ a été utilisé pour la correction de lecture PacBio en utilisant toutes les données Illumina PE disponibles et les unitigs créés par MaSuRCA v2.3.2 ⁷⁵ . SEECER v0.1.3 ⁷⁶ et Rcorrector v1.0.2 ⁷⁷ ont été utilisés pour RNA-seq et Musket v1.1 ⁷⁸ pour la correction des appels de base ADN-seq. Les ensembles de données DNA-seq et RNA-seq ont été prétraités en utilisant le même pipeline (avec des paramètres différents) ; en général, deux régimes de filtrage ont été appliqués à chaque ensemble de données (« strict »/« standard » et « détendu ») en préparation de différents cas d'utilisation en aval (voir le tableau supplémentaire S22 ). La taille du génome a été estimée par une approche basée sur le spectre k-mer mise en Å“uvre dans GCE v1.0.2 ⁷⁹ .

Assemblage du génome

Du point de vue du méta-assemblage réalisé, l'algorithme implémenté dans MaSuRCA v2.3.2 ⁷⁵ (utilise Celera Assembler v6.5 ⁸⁰ ) a servi de procédure d'assemblage de base ; Tous les ensembles de données disponibles à l'heure actuelle (c'est-à-dire Illumina PE et MP, Illumina Nextera MP et 454 MP et SE) ont été utilisés. Celera Assembler v8.3rc2 (CA) ⁸⁰ a été utilisé pour les assemblages « PacBio uniquement ». Comme plusieurs individus par espèce (tous des diploïdes non consanguins) ont été séquencés dans ce projet, l'hétérozygotie était une préoccupation. Par conséquent, des algorithmes d'assemblage spécialement conçus pour mieux gérer la divergence ont été incorporés dans l'approche de reconstruction : Platanus v1.2.1 ⁸¹ est un assembleur récent conçu pour traiter de manière plus judicieuse les problèmes d'hétérozygotie dans les données génomiques (5 itérations ; tous les ensembles de données Illumina ont été utilisés) ; Redundans v0.12c ⁸² (utilise SSPACE3 ⁸³ , GapCloser ⁸⁴ , bwa ⁸⁵ et last ⁸⁶ ) vise également à fournir des assemblages plus précis et contigus de génomes hautement hétérozygotes (5 itérations ; tous les ensembles de données Illumina ont été utilisés). La suite PBJelly v15.8.24 ⁸⁷ (utilise BLASR ⁸⁸ ) a été utilisée pour incorporer les lectures à séquence longue (PacBio) dans un processus d'assemblage guidé par référence dans les ébauches établies (5 itérations). L'ensemble diversifié de projets de génome générés a été soumis à Metassembler ⁸⁹ dans le but de générer des séquences consensus de haute qualité. Un algorithme personnalisé, qui prend en compte plusieurs mesures sur les erreurs d'assemblage probables, la contiguïté et les prédictions génétiques (tirant des informations de QUAST ⁹⁰ et REAPR ³⁸ ), a été appliqué pour déterminer le meilleur ordre des méta-assemblages successifs.

Finition du génome

Pour une autre série de fermeture de l'espace entre les échafaudages, GMcloser ⁹¹ (utilise Nucmer ⁹² / BLAST ⁹³ et​​Bowtie2 ⁹⁴ ) a été appliqué sur les méta-assemblages avec les données PacBio et Illumina PE. Enfin, Sealer ⁹⁵ (utilise Konnector , une partie du pipeline assembleur ABYSS ^{96 ) a été utilisé avec les bibliothèques Illumina PE (libérales) pour le remplissage final des lacunes et une finition personnalisée du génome 42} basée sur GATK (via la rétro-cartographie Illumina PE et le rappel de consensus) a été appliquée.

Validation du génome

REAPR v1.0.18 ³⁸ (utilisant SMALT v0.7.0.1 ⁹⁷ ) a été utilisé avec les bibliothèques Illumina Nextera mate-pair (6 kpb) et Illumina PE (600 pb) pour évaluer l'exactitude des assemblages et QUAST v4.1 90 a ^été appliqué pour les statistiques de contiguïté et de prédiction génétique. L'exhaustivité des assemblages a été évaluée à l'aide de CEGMA v2.5 ³⁵ (en utilisant GeneWise v2.4.1 ⁹⁸ , HMMER v3.0 ⁹⁹ et NCBI BLAST + v2.2.29 +  ⁹³ ) avec optimisation des paramètres pour les génomes de vertébrés (–vrt) et BUSCO v3.0.2 ³⁴ (en utilisant NCBI BLAST + v2.2.29 + , HMMER v3.1 ⁹⁹ et AUGUSTUS v3.2.1 ¹⁰⁰ ).

Assemblage du transcriptome et cartographie des lectures d'ARN-seq

Les assemblages du transcriptome ont été réalisés avec Trinity v2.3.2 ^101 , 102 et le pipeline PASA2 v2.0.2 ¹⁰³ (en utilisant GMAP v2014-12-06 ¹⁰⁴ , BLAT v36.1 ¹⁰⁵ et MySQL v5.7.12 ¹⁰⁶ ) ; ^{Transdecoder v3.0.1 102} a également Ã©té appliqué pour identifier les régions codantes candidates (utilisées avec MAKER3 ¹⁰⁷ ). Les alignements de lecture d'ARN-seq pour d'autres analyses ont généralement été effectués avec STAR v2.4.2a ¹⁰⁸ en utilisant les paramètres par défaut.

Annotation du génome

Des annotations structurelles ont été réalisées sur la base de données expérimentales provenant d'ensembles de données d'ARNm-Seq. De plus, des informations ont été tirées de modèles de transcription et de protéines provenant d'ensembles de données sélectionnés accessibles au public ( Danio  rerio , H. burtoni, M. zebra, N. brichardi, O. niloticus et P. nyererei ) et d'autres modèles dans UniProt | , nr/nt et UniRef90|teleost. L'annotation fonctionnelle a été principalement réalisée via des comparaisons basées sur BLAST avec les bases de données mentionnées et via une multitude de bases de données coordonnées par InterProScan 5 (voir Tableau 2 ).

L'annotation structurelle des gènes codants et des ARNt a été générée à l'aide des pipelines MAKER v3.0 ¹⁰⁷ (en utilisant les chercheurs de gènes GeneMark-ES v4.32 ¹⁰⁹ , AUGUSTUS v3.2.1 ¹⁰⁰ , SNAP v2013-11–29 ¹¹⁰ et tRNAscan v1.3.1 ¹¹¹ ) , Funannotate v0.5.5-v0.7.0 ¹¹² (FA) et BRAKER1 v1.9 ¹¹³ (utilisant GeneMark-ET v4.32 ¹¹⁴ et AUGUSTUS v3.2.1 ¹⁰⁰ ) ; BRAKER1 a également été utilisé pour la formation AUGUSTUS. De plus, des modèles génétiques ont été créés avec StringTie v1.3.2d ¹¹⁵ et Cufflinks v2.2.1 ¹¹⁶ . Tous les modèles ont été combinés par EVidenceModeler v1.1.1 ¹¹⁷ (EVM) sous le contrôle de MAKER3. Pour les ARN non codants, Infernal v1.1.2 ¹¹⁸ , Rfam v12.1 ¹¹⁹ et FEELnc v0.1.0 ¹²⁰ ont été utilisés. L'ensemble de formation d'ARNm pour FEELnc a été dérivé des données d'annotation FA/MAKER, où des modèles de gènes présumés « bons » avec une structure similaire aux modèles précédemment publiés ont été sélectionnés ; l'ensemble de formation lncRNA a été généré par mélange des séquences d'ARNm. Les microsatellites ont été appelés avec MISA v1.0 ¹²¹ , les îles CpG avec EMBOSS v6.6.0 ¹²² cpgplot et les ORF avec EMBOSS v6.6.0 getorf (et post-traitement R). Les répétitions ont été déterminées à l'aide de RepeatMasker v4.0.6 ³² (avec RepBase v20160321 ¹²³ et de bibliothèques spécifiques aux espèces générées avec RepeatModeler v1.0.8 ¹²⁴ ), RepeatScout v1.0.5 ¹²⁵ et TRF v406 ¹²⁶ .

L'annotation fonctionnelle a été réalisée à l'aide d'InterProScan v5.24–63.0 ¹²⁷ (en utilisant les bases de données CDD-3.14, Coils-2.2.1, Gene3D-3.5.0, Hamap-201605.11, MobiDBLite-1.0, PANTHER-11.1, Pfam-30.0, PIRSF- 3.01, PRINTS-42.0, ProDom-2006.1, ProSitePatterns-20.119, ProSiteProfiles-20.119, SFLD-2, SMART-7.1, SUPERFAMILY-1.75, TIGRFAM-15.0 et TMHMM-2.0c). De plus, sous le contrôle de FA, les bases de données eggNOG v4.5.1 ¹²⁸ (fiNOG), MEROPS v12.0 ¹²⁹ , dbCAN v5.0 ¹³⁰ et BUSCO vertebrata v3 ³⁴ ont été utilisées pour les recherches de similarité et SIGNALP v4.1 ¹³¹ pour l'identification de l'emplacement cible. séquences de signaux.

L'intégration finale de toutes les annotations a été réalisée avec R 3.4.3/Bioconductor 3.6 en utilisant les packages data.table 1.12.2, GenomicFeatures 1.30.3, VariantAnnotation 1.24.5 et leurs dépendances.

Cartographie de lecture de séquences d'ADN

Les lectures prétraitées ont été alignées en mode paire avec BWA mem ⁸⁵ en utilisant les paramètres par défaut avec les indicateurs ‐M et ‐R. Les lectures alignées ont été triées par coordonnées avec Picard SortSam v1.119 ⁷³ et indexées avec l'index SAMtools v1.4 ⁷² . Les doublons ont été supprimés avec Picard MarkDuplicates v1.119. La qualité des cartographies a été évaluée avec QualiMap v2.0 ¹³² .

Analyse comparative â€“ appel de petites (SMV) et de variantes structurelles (SV) â€“ prédiction de l'effet des variantes

La boîte à outils d'analyse du génome (GATK) v3.7 a été utilisée pour le réalignement local des lectures ainsi que pour la détection et le filtrage des variantes SNP/InDel (appelées petites variantes, SMV) ⁴² , comme recommandé par la documentation GATK ; le HaplotypeCaller a été appliqué - avec un score minimum pour l'émission de variantes de 10, pour l'appel de 30 et un élagage minimum de 10. Les SMV avec un score de qualité ≥ 30 ont été inclus dans des analyses plus approfondies. DELLY v0.7.7 ⁴³ a été appliqué pour appeler des variantes structurelles (SV, insertions, suppressions, duplications, inversions et translocations) avec une taille limite d'insertion de 3 (pour les suppressions) et une qualité minimale de mappage des extrémités appariées de 20. Toutes les variantes avec un minimum de 5 paires de lectures brisées prenant en charge la variante ainsi qu'une longueur minimale de 300 pb (pour les suppressions, les inversions et les duplications) ont été incluses dans des analyses plus approfondies, comme recommandé par la documentation DELLY. Des effets de variantes présumés ont été appelés avec SNPeff v4.3r ³⁷ . Whippet v0.11.1 ¹³³ a été utilisé pour l'appel d'événements d'épissage alternatifs. Les analyses comparatives ont été effectuées dans R 3.4.3/Bioconductor 3.6 en utilisant les packages data.table 1.12.2, GenomicFeatures 1.30.3, VariantAnnotation 1.24.5 et leurs dépendances.

Analyse GO

Pour affiner la liste des gènes candidats, une analyse d'enrichissement GO a été réalisée sur les régions génétiques portant des variants à l'aide du package R topGO v2.30.1 ¹³⁴ ; les annotations GO personnalisées ont été générées sur la base des mappages InterProScan. La topologie GO a été prise en compte ( poids de la méthode ) et l'enrichissement a été évalué via un test exact de Fisher avec un seuil de p  ≤ 0,001. Voir les détails de l'analyse GO dans les informations supplémentaires.

Approbation éthique et consentement à participer

Le traitement des animaux rapporté dans cet article est conforme aux normes de la loi autrichienne sur la protection des animaux et de la directive communautaire 86/609 de la Communauté européenne. Les poissons ont été conservés dans notre aquarium certifié à l'Institut de biologie de l'Université de Graz. Les individus ont été échantillonnés par CSt et SK, euthanasiés à l'aide d'une surdose d'huile de clou de girofle et décapités conformément à la législation autrichienne sur le bien-être des animaux. Conformément à la loi autrichienne sur l'expérimentation animale (TVG, BGBI. Nr. 501/1989, modifiée en dernier lieu par BGBI. I Nr. 162/2005), l'approbation n'était pas requise car aucun traitement expérimental n'était effectué.

Consentement à la publication

N'est pas applicable.