Development of a relevant strategy using de novo transcriptome assembly method for transcriptome comparisons between Muscovy and common duck species and their reciprocal inter-specic Mule and Hinny hybrids fed ad libitum and overfed

Background: Common Pekin and Muscovy ducks and their intergeneric hinny and mule hybrids have different abilities for fatty liver production. RNA-Seq analyses from the liver of these different genetic types fed ad libitum or overfed would help to identify genes with different response to overfeeding between them. However RNA-seq analyses from different species and comparison is challenging. The goal of this study was develop a relevant strategy for transcriptome analysis and comparison between different species. Results: Transcriptomes were rst assembled with a reference-based approach. Important mapping biases were observed when heterologous mapping were conducted on common duck reference genome, suggesting that this reference-based strategy was not suited to compare the four different genetic types. De novo transcriptome assemblies were then performed using Trinity and Oases. Assemblies of transcriptomes were not relevant when more than a single genetic type was considered. Finally, single genetic type transcriptomes were assembled with DRAP in a mega-transcriptome. No bias was observed when reads from the different genetic types were mapped on this mega-transcriptome and differences in gene expression between the four genetic types could be identied. Conclusions: Analyses using both reference-based and de novo transcriptome assemblies point out a good performance of the de novo approach for the analysis of gene expression in different species. It also allowed the identication of differences in responses to overfeeding between Pekin and Muscovy ducks and hinny and mule hybrids.


Background
In waterfowls, liver steatosis can occur spontaneously as a result of energy storage before migration.
This natural ability has been exploited for thousand years to produce fatty liver or "foie gras" by overfeeding ducks or geese in a short period of time. However, fatty liver production varies according to species and breeds [1][2][3][4][5]. In France, the main producer worldwide, most of the production (90%) involves Mule ducks, i.e. intergeneric hybrids between a female common Pekin duck (Anas plathyrynchos) and a male Muscovy duck (Cairina moschata), and in a far lesser extent Muscovy ducks or geese from the Landes grey breed (Anser anser). Conversely, Pekin ducks are not involved due to their lower ability to produce fatty liver. Some studies have been conducted to better characterize hepatic steatosis development in waterfowls and analyze differences between genotypes. However, they mainly focused on biochemical levels [1][2][3][4][5]. Some studies also focused on gene expression levels but were conducted on few candidate genes only [6][7][8][9]. This was mainly due to the lack of duck speci c microarrays for transcriptome analyses at the genome level. Thus, genome-wide analyses of gene expression are still needed to better characterize hepatic steatosis in ducks and their hybrids.
The recent development of high throughput RNA sequencing now allows characterization of expressed transcripts and quanti cation in samples of any species without the need of preexisting tools such as microarrays [10,11]. Different strategies were applied for transcriptome assembly and analyses of RNA sequences referred as reference-based or de novo assemblies [12][13][14][15]. Reference-based strategy is generally conducted when a reference genome for the target transcriptome is available. In such a situation, RNA sequences are aligned on the reference genome using mapping tools like Burrows-Wheeler transform (BWA) [16] or TopHat-Cu inks [17]. This approach has also been applied in heterologous situations where reference genome and RNA sequences were from different species, for examples mapping horse, donkey and their hybrids transcriptomes on horse reference genome [18], sparrows transcriptome on zebra nch reference genome [19] or red deer transcriptome on cow reference genome [20]. In a previous study we have also conducted such a strategy on common and Muscovy ducks and their interspeci c hybrids using the common duck genome as reference [21]. This strategy was relevant to analyze differences in gene expression within a genetic type as reads from the same type would have the same mapping rate. However, a bias was expected in read mapping and counting according to different homologous, "half-homologous" and heterologous situations making this strategy not relevant to analyze differences in gene expression between genetic types. In such heterologous situations, i.e. in the absence of reference genome from the same species or when different species and hybrids are compared, de novo approaches are more generally conducted using de novo transcriptome assembling tools such as Trinity [22] or Oases [23]. Some examples of transcriptome analyses in interspeci c hybrids are the studies conducted in cyprinidae [24] or in brassica [25]. As different tools are available for reconstructing transcripts de novo, the choice of an optimal method is challenging. Surprisingly, different de novo assembly tools and methods were compared using common duck RNA sequences [26] while common duck reference genome is available [27]. Although different strategies can be applied to analyze transcriptomes from different species, no one appears as the best one. Interestingly, a de novo RNA-Seq Assembly Pipeline (DRAP) has been developed improving de novo transcriptome assemblies performed by Trinity and Oases in terms of compaction (number of contigs needed to represent the transcriptome) and quality (chimera and nucleotide error rates) [28].
In the present study, reference-based approach using Anas platyrhynchos genome as reference and de novo assembly approach were used and compared to analyze hepatic gene expression of overfed and ad libitum fed Muscovy and Pekin ducks and their reciprocal inter-speci c Mule and Hinny hybrids and to identify differentially expressed genes between feeding and duck genetic types.

Animals and experimental design
Animals, experimental design, RNA preparation and sequencing have been described in previous publications [5,6,8,21]. Brie y, male ducks corresponding to common Pekin (Anas platyrhynchos) and Muscovy (Cairina Moschata) duck species and to their two reciprocal interspeci c Mule and Hinny hybrids were fed ad libitum or overfed 14 days with corn and corn meal at Experimental Station for Waterfowl Breeding (INRA, Artiguères, France). After electronarcosis, ducks were euthanized by neck sectioning and liver samples were collected and frozen in liquid nitrogen [5]. RNA were extracted from liver samples [6], cloned in libraries and paired-ends sequenced [21]. The whole transcriptomic project included 80 Illumina libraries for a total of 2 252 041 541 paired-end   reads with a read length of 100 bp. The four genetic types (Anas platyrhynchos, Cairina moschata, Hinny  and Mule hybrids) correspond to 20, 19, 20 and 21 libraries for 554 893 036, 513 608 258, 544 599 927 and 638 940 320 paired end reads in total for each of them, respectively. RNA sequences are accessible in the NCBI sequence read archive (SRA) under the accession number SRP144764.
Transcriptome assembly using reference genome After adapter trimming, mapping of RNA sequences and transcripts reconstruction on the common duck (Anas platyrhynchos) reference genome BGI_duck_1.0 (INSDC Assembly GCA_000355885.1) [27] were performed using TopHat2, Cu inks and Cuffmerge tools of the Tuxedo suite [17] as previously described [21].
De novo transcriptome assemblies RNA-seq data were alternatively assembled using de novo approaches. The raw datasets for each of the 4 species were assembled independently using Trinity v2.3.2 [22,29], including a reads cleaning process with Trimmomatic (using defaults parameters) [30] and an in silico reads normalization step (provided within the Trinity suite). Parental transcriptome (mixing Pekin and Muscovy normalized reads) and hybrid transcriptome (mixing Hinny and Mule duck normalized reads) were assembled with Trinity and Oases (m 25 -M 65 -s 10) [23]. Composite meta-transcriptome assembly of single genetic type transcriptomes was also conducted using DRAP v1.91 [28] for compacting and correcting Trinity assemblies and constructing a reference transcriptome. By default, DRAP produced 4 assemblies which were ltered for low coverage with Fragments Per Kilobase per Million (FPKM) thresholds set at 1, 3, 5 and 10 respectively.

Assembly assessment
Quanti cation for the abundances of transcripts was performed by library using pseudo-aligner Kallisto [31] (trinityrnaseq-2.4.0/util/align_and_estimate_abundance.pl) in no strand-speci c paired-end mode. Transrate [32] assembly evaluation which was regrouped in DRAP package was set to version 1.0.1.

Results
Transcriptome assembly using reference genome Transcriptome assemblies of RNA sequences from Pekin (Anas platyrhynchos), Muscovy (Cairina moschata), Hinny and Mule ducks were performed as in a previous study [21] by a reference-based approach using the common duck BGI_duck_1.0 assembly as reference genome. Heterologous mapping rate of Muscovy and to a lesser extend of Hinny and Mule transcriptomes (40%, 58%, 59%, respectively) were lower than homologous mapping rate of Pekin transcriptome (71%) on this common duck reference genome, clearly indicating an important bias in heterologous mapping when compared to homologous mapping. Direct mapping of reads on reference genome produced similar results (data not shown). This mapping bias prompted us to test the relevance of a de novo approach.

De novotranscriptome assemblies
After adapters removing, quality trimming and ltering of reads, 545 691 171, 504 591 419, 534 967 150 and 627 575 995 pairs were conserved for Pekin, Muscovy, Hinny and Mule ducks, respectively. Subsequently, 37 073 089, 36 822 075, 43 650 600 and 46 866 094 pairs were selected during normalization. They were assembled de novo using Trinity in four independent "single genetic type" transcriptomes, in a "mixed parental" transcriptome (with Pekin and Muscovy reads) and in a "mixed hybrids" transcriptome (with Hinny and Mule reads). A mixed hybrids transcriptome was also assembled using Oases for comparison. As shown in Table 1, single parental species Pekin or Muscovy transcriptome assemblies with Trinity were of higher quality when compared to single hybrid Hinny or Mule transcriptome assemblies with greater average lengths, higher N50 values and lower numbers of transcripts. The quality of the mixed parental (Pekin+Muscovy) transcriptome assembly was also lower, similar to those of single hybrids transcriptome assemblies. The quality of the mixed hybrids (Hinny+Mule) transcriptome assembly was even lower. When mixed hybrids transcriptome assembly was conducted with Oases, very long transcripts with a great N50 value were produced suspected to be chimeras.
Unique meta-transcriptomes were nally assembled with DRAP using single genetic type Pekin, Muscovy, Mule and Hinny assemblies and after ltering at 4 different FPKM thresholds. This ltering on FPKM values resulted in a reduction of transcript number by removing lowly expressed transcripts. For example, 175 572 sequences were conserved after ltering of FPKM > 1 and 64 862 with FPKM > 3, representing 28 6879 845 and 127 395 003 residues, respectively. We chose to keep the transcriptome assembly produced by DRAP after FPKM > 1 ltering for further analyses. This assembly was of high quality, with a N50 value of 2807 bp. Pseudo-alignments of reads on this meta-transcriptome and quanti cation were performed using Kallisto. As shown on Fig. 2, pseudo-alignment rates of sample reads on DRAP metaassembly were very similar to pseudo-alignment rates on Pekin Trinity assembly, with low "heterologous biases" and variability across samples (84.56 ± 0.98 % and 85.63 ± 1.10, respectively).
Annotation completeness in terms of gene content of de novo transcriptome assembled with DRAP was assessed with BUSCO for the presence/absence of the conserved eukaryotic single copy orthologous genes. To refer to our previous work [21], completeness of transcriptome assembly using mapping approach with Cu inks and BGI_duck_1.0 as reference genome was also assessed. As shown in Table 2, DRAP de novo approach produced more orthologues (97.1 % completeness), provided a more complete (0 % missing) and less fragmented (only 3%) catalog of orthologues when compared to the mapping approach using Cu inks (81.5 %, 7.3 % and 11.2 %, respectively). Same BUSCO analyses were also conducted on BGI_duck_1.0 genome assembly and on two more recent genomes assembled at the chromosome level and expected to be more complete, i.e. CAU_duck1.0 (GCA_002743455.1) and ASM874695v1 (GCA_008746955.1) ( Table 2). These three genomes contain more complete and single copy orthologous genes (81.2%, 82.8% and 90.1%, respectively) than de novo transcriptome assembly with DRAP (59.1%) but less complete orthologues in total (82.9 %, 83.5% and 90.8% versus 97.1%) with more genes missing (7.2 %, 7.3% and 6.9% versus 0.0%) suggesting that our RNA sequencing and de novo transcriptome assembly allowed identifying new genes that are missing in genome assemblies.

Gene expression analyses
Gene expressions after DRAP de novo approach were analyzed with edgeR. Signi cant (p-value < 0.01) differentially expressed genes (fold change ≥ 2) were determined. In total, 13898 different genes were found up-or down-regulated by overfeeding in the four genetic types (Table 3). Much less differentially expressed genes (DEG) were identi ed in Pekin than in other genetic types (3749, versus 6167, 8920 and 7696, respectively). Among all 13898 DEG, 903 were identi ed in the four genetic types. Similar results were found when analyses were conducted with DESeq2 (data not shown).
Expression levels of DEG were analyzed by principal component analysis (PCA) and hierarchical clustering (HC) to cluster samples according to similarities in gene expression. As shown in Fig. 3A, the rst principal component (PC 1) of PCA summarized 50% of the whole variability and discriminated samples according to genetic type, pure species being extreme and hybrids intermediate. The second principal component (PC2) summarized 16% of the whole variability and discriminated samples according to feeding. The cluster corresponding to overfed Pekin ducks appeared more dispersed than the other clusters. Two clusters were de ned for Mule and Hinny samples according only to feeding (ad libitum and overfed) without any distinction according to genetic type. When expressions of DEG were analyzed by HC (Fig. 3B), two clusters were rst de ned: one including all Pekin samples, whatever the feeding was, and the other including all other samples. This second cluster included two groups. One corresponded to ad libitum and overfed Muscovy samples. The second included two groups: on one hand ad libitum and on the other hand overfed mule and hinny hybrid samples as was observed in PCA. Finally, 4 different clusters were de ned again according rst to genetic type: a Pekin cluster more distant from the 3 others and including ad libitum and overfed samples; a Muscovy cluster also including ad libitum and overfed samples; a fed ad libitum Hinny and Mule cluster; and an overfed Hinny and Mule cluster. As observed in PCA, these two latter clusters suggest that Hinnies and Mules in the same feeding status cannot be distinguished according to differential gene expression.
The presence of down-and up-regulated DEG in the four genetic types was shown with Venn diagrams (Fig. 4A). Few down-and up-regulated DEG (364 and 539, respectively) were identi ed in the four genetic types. However, fold change in expressions of these common DEG were very different according to genetic type (Fig. 4B). DEG expression levels in Pekin were more different from the three other genetic types with less down-and up-regulated DEG visualized, suggesting that DEG responses to overfeeding are less important in Pekin ducks. Conversely, many DEG were not found in each of the 4 genetic types or speci cally found in one genetic type only (821, 2171, 1281, and 2125 in Pekin, Muscovy, Hinny and Mule ducks, respectively) indicating genetic type effects (Fig. 4A). Some examples taken at random of such feeding and genetic type effects interactions are shown in additional le 1. They illustrate different types of response to overfeeding (fold change and/or up and down regulation) between the 4 genetic types. Hybrids shared the highest number of DEG with 2677 down-and 3061 up-regulated DEG found in both mule and hinny ducks. Only two genes were found with interaction, i.e. down-regulated in one hybrid and up-regulated in the other hybrid. These results indicate that few differences were observed in response to overfeeding between the two hybrids or in other words few genetic type effects and therefore few interactions with feeding effect. For comparison, less down-and up-regulated DEG were found in both mule and Muscovy ducks (1640 and 1936, respectively), even less in both mule and Pekin ducks (956 and 1611, respectively) and only 415 down-and 598 up-regulated DEG in both Pekin and Muscovy ducks indicating more differences in feeding effect between the two duck "pure" species.

Functional enrichments
Functional enrichments in GOBP associated to DEG identi ed after DRAP de novo method were analyzed and compared to those identi ed after reference based method with BGI_duck_1.0 genome as conducted previously [21]. Many GOBP enrichments (535) were found whatever the method used (additional les 2 and 3A). However, some differences were observed between methods as evidenced with enrichment annotation pro les (Fig. 5). Enrichments were less important with de novo method. Difference between Pekin and the 3 other genetic types was more important with de novo method and many GOBP enrichments found with reference based method (744) were not found with DRAP method (additional le 3C). Interestingly, many GOBP were found enriched with DRAP method (1257), but most of them corresponded to few DEG (additional le 3B) and/or were found in 1 or 2 genotypes only (additional le 2).
As expected, most of the enriched biological processes found with the 2 methods were related to metabolism of lipids but also to many other processes related to mitosis and cell cycle, transmembrane transport, ion homeostasis and in ammatory response to cite the more enriched (additional le 2). To better characterize lipid metabolism, we focused on GOBP terms including lipid or fatty acid words. As shown in Fig. 6, metabolism and regulation of fatty acids and neutral lipids, lipid signaling and response to lipid biological processes were enriched in Muscovy ducks and the two hybrids and much less or not in Pekin ducks. Again, difference between Pekin and the 3 other genetic types was more important with de novo method.

Discussion
The aim of this study was to analyze and compare gene expressions in four different duck genetic types, common Pekin and Muscovy duck species and their reciprocal mule and hinny hybrids using a relevant approach. As mentioned by Moreton et al. [14], choice for using reference-based or de novo transcriptome assembly approach for gene expression analyses is generally based on the question whether or not a reference genome is available. However, when gene expressions from different species are to be analyzed and compared, the question of strategy is more complex and cannot be taken up with the question of availability of a reference-genome and neither strategy seems a priori better than the other.
In a previous study, we have conducted a reference-based approach using common duck (Anas platyrhynchos) BGI_duck_1.0 genome assembly as reference for the four genetic types [21]. As indicated, mapping rates were very different when homologous and heterologous mapping were considered (71% with Pekin duck reads and 41% with Muscovy duck reads). In the study of Wang et al. [18], homologous mapping rate of horse reads on horse reference genome (61.70%) was very similar to heterologous mapping rate of donkey reads on horse reference genome (62.07%). The reason of this difference of heterologous mapping rate between equids and ducks was not directly investigated but we can speculate that it is partly due to the difference of evolution distance between equids and duck species (estimated time according to http://www.timetree.org/, [45], 7.7 and 22.8 MYA, respectively), but also to different mapping parameters in the two studies. In our study, these parameters were set by defaults and were not optimized for heterologous mapping. As described in our previous study [21], this bias of mapping was not a problem to analyze differences in gene response to feeding in ducks from the same genetic type. It was however too important for gene expression comparisons between the four genetic types and to analyze genetic effect in gene expression. Therefore, de novo transcriptome assembly approaches were then conducted. As shown in Table 1 and Fig.1A, Trinity "single" transcriptome assemblies from Pekin or Muscovy duck sequences were very similar and better in terms of length and alignments of reads than assemblies from hinny and mule hybrid sequences. This result is probably due to the fact that in hybrids transcripts are expressed from two different genomes making assembly less e cient according to the presence of some polymorphisms in transcripts from the same orthologs between the two species. This was also observed when assemblies were performed with sequences from two different species (Fig. 1B). We have also shown that more reads were aligned on these de novo assembled transcriptomes than on common duck reference genome BGI_duck_1.0. Furthermore, homologous mapping of Pekin reads on de novo Pekin transcriptome assembly (around 85%) was more e cient than that on reference genome (71%). We can speculate that some reads correspond to genes that are missing in this reference genome, corresponding to new or misassembled genes, or to genes that are in unknown regions, i.e. unassembled and unassigned regions. Thus, these reads are not mapped on reference genome but are assembled together in transcripts by de novo approach and map on theses transcripts.
As Trinity and Oases were not well suited for "complex" transcriptome assembly, we have tested de novo RNA-Seq Assembly Pipeline (DRAP) [28] to perform a meta-transcriptome assembly of transcriptomes from the four duck genetic types. As shown in Fig. 2, this meta-transcriptome assembly with DRAP was very e cient, with a mapping rate of sample reads similar to the best rates previously observed, i.e. homologous mapping rate of Pekin or Muscovy duck reads on de novo assembled Pekin or Muscovy transcriptomes. Further analysis with BUSCO [33] con rmed the quality of this meta-transcriptome assembled with DRAP (Table 2). In our study, more than 97% of expected single-copy orthologues were assembled by DRAP as complete (i.e. full expected transcript length), 3% were fragmented (assembled with a shorter length than expected) and nally no (0%) were missing. With reference-based approach using Cu inks and BGI_duck_1.0 assembly, only 81.5% single-copy orthologues were assembled to full length, 11.2% were fragmented and 7.3% were missing. This clearly indicates that de novo assembly of meta-transcriptome using DRAP is of higher quality when compared to assembly using reference-based approach as performed earlier [21]. Again, this is probably due to the missing of genes in thisgenome. Moreover, when BUSCO was applied directly to this genome, proportions of complete and missing orthologues were similar (82.9% and 7.2%, respectively) to those found after reference-based approach for transcriptome assembly using Cu inks. BUSCO analyses were also conducted on CAU_duck1.0 and ASM874695v1 common duck genomes assembled more recently. As expected, completeness of these two genomes was more important (83.5% and 90.8%, respectively). However, completeness were lower than that of DRAP transcriptome assembly and no improvement on missing orthologues was observed (7.3% and 6.9%), again suggesting that de novo transcriptome assembly with DRAP was more complete.
It is interesting to note that lower proportions of complete and single copy orthologs were found in DRAP assembled transcriptome and Cu inks transcriptome assembled with BGI_duck_1.0 genome (59.1% and 42.9 %) than in all three different assemblies of the genome (81.2%, 82.8% and 90.1%). The fact that most of orthologues were complete and as expected in a single copy in the genome and not in transcriptomes whatever the approach used is probably the result of the presence of only one isoform for one orthologue in the genome (by de nition because BUSCO applies only on single-copy orthologues in related species) but of multiple alternative transcripts from one single-copy orthologue making assembly in one transcript more tricky. This clearly indicates that development of methods to improve transcriptome assembly and reduce duplication and more generally redundancy remains a concern today. Nevertheless, our results show that the de novo approach with DRAP is more relevant than the referencebased approach to assemble a meta-transcriptome from reads of the four duck genetic types.
Gene expression analyses were then conducted using this meta-transcriptome as reference, assuming that mapping bias were substantially reduced and allowed us to analyze DEG within and between genetic types. Many genes were found up-and down-regulated by overfeeding and some of them (903) were found in the four genetic types. These results indicate some similarities but also some genetic typespeci c responses to overfeeding. This was further con rmed by principal component analyses and hierarchical clustering, clearly showing that Pekin ducks are different from the three other genetic types.
Functional enrichments in GO biological processes were nally conducted on DEG identi ed by DRAP de novo method and compared to those of DEG found with reference based method using BGI_duck_1.0 genome. Most of the most enriched biological processes identi ed after reference based method, notably lipid metabolism, were also observed with DRAP method, indicating that most of our previous conclusions remain valid when using de novo method. However, it is important to note that enrichments appeared overestimated with reference based method when compared to de novo method. According to mapping rates on BGI_duck_1.0 genome and DRAP meta transcriptome, we assume that de novo method is more relevant to compare DEG and associated functions between the four genetic types. It would have be interesting to compare the two methods using a more complete duck genome assembly, ASM874695v1 for example, to apply the reference based method. However, as this genome is less complete than meta-transcriptome assembled de novo and more importantly is always fully speci c to Pekin duck only, we think that mapping biases would be similar to those observed with BGI_duck_1.0 genome and the results as well.
Many enriched biological processes speci cally found with DRAP method were represented by very few DEG (additional le 3B). It is di cult to draw de nitive conclusion about this result without further studies on related genes. Conversely, no unexpected biological process was highly enriched and speci cally found using DRAP method. Different conclusions can be drawn from these results. First, responses to overfeeding in hinny and mule hybrids were very similar and indicate that expectation that hinnies are not used for fatty liver production due to a different and lower response to overfeeding is wrong. In fact, the reason essentially lies in the di culty to produce viable hinnies due to a lower e ciency of arti cial insemination of Muscovy females and their lower fertility when compared to Pekin females. Second important conclusion is that many genes are differentially expressed in the four genetic types, including those related to lipid metabolism, but also related to many other biological processes. As previously observed [21],this suggests that response to overfeeding and consequently development of hepatic steatosis is a complex trait involving many genes and functions and is not well described by the expression of few candidate genes [6][7][8][9], even if expressions of these genes are well correlated to fatty liver weight [6]. Third, although some differentially expressed genes were found in the four genetic types, some other were not found in Pekin ducks or expressed at a lower level, including those involved in fatty acid and lipid metabolism making this species less similar and more distant to the three other genetic types. Furthermore, individual variability of responses in Pekin ducks was very important making some overfed Pekin samples near to Pekin ducks fed ad libitum samples. This is directly linked to the greater variability and lower ability of Pekin ducks to produce fatty liver in response to overfeeding when compared to hinny and mule hybrids and Muscovy ducks. To conclude on this point, response to overfeeding in ducks is very complex, different in some extent from a genetic type to the other and involves many genes suggesting that genetic improvement of duck selection for fatty liver production should be undertaken at the whole genome level.

Conclusions
We have shown that differential gene expression analyses in different duck genetic types and comparisons are more accurate when a meta-transcriptome is assembled with DRAP from single transcriptomes of the four different genetic types assembled by Trinity. Using this strategy, we have shown that differences in responses to overfeeding can be identi ed between the four genetic types and that Pekin ducks greatly differ from ducks of the 3 other genetic types. This work provides new information to describe hepatic steatosis in different duck genetic types in a much more holistic way.
This information could also be exploited in the context of genomic selection programs that will certainly be developed in ducks in the coming years.

Declarations
The dataset supporting the conclusions of this article is available in the NCBI Sequence Read Archive (SRA), under the accession number SRP144764 (https://www.ncbi.nlm.nih.gov/sra/SRP144764).
Authors' contributions CD and EC conceived and coordinated the study. XL and FH performed the analyses. XL, FH, EC and CD interpreted the results and drafted the manuscript. All authors read and approved the nal manuscript.
Ethics approval and consent to participate All experimental procedures were performed in accordance with the French National Guidelines for the care and use of animals for research purposes (http://ethique.ipbs.fr/charteexpeanimale.pdf). Ethics Committee authorization to experiment on living animals was not sought for this study as samples were previously collected by another study.

Consent for publication
Not applicable.

Competing interests
The authors declare that they have no competing interests.