Nutrigenomics in Animal Feeding: Digital Gene Expression Analysis in Poultry Fed Tenebrio molitor Larvae Meal

The purpose of this study was to investigate the effects of high levels of Tenebrio molitor dietary inclusion (15%) on molecular mechanisms that influence poultry health in a broiler chicken diet. The global gene expression of four tissues (breast, liver, jejunum, and caecum) was evaluated using the RNA-Seq approach. The analysis of differentially expressed genes suggested that the use of Tenebrio molitor leads to the overexpression of genes related to protein elongation required for tissue growth and development in the gut and liver. It would also appear to contain nutrients that reduce the expression of genes related to the immune system and inflammation of the mucosa. The dietary inclusion of Tenebrio molitor in poultry could also lead to a possible inactivation of the growth factor and a reduction of tissue free-radicals. No genes alterations have been detected in liver RNA expression that would discourage the use of larvae in feeding broilers.


Introduction
In the modern world, human and animal populations experience exponential growth, and the consequent scarcity of protein sources leads to an increasing interest in insect meal for animal and human nutrition. Insect production requires less energy, less land utilization, and consequently, exerts a lower environmental impact [1][2][3][4]. Several studies highlighted the opportunity to use insects, containing a large amount of high-quality protein, as a standard ingredient in animal feeds [5][6][7][8][9][10]. Furthermore, insect meal as a protein source represents a unique opportunity for animal nutrition due to its low competitiveness with human nutrition [11]. In particular, in poultry production (both egg and meat), which is rising due to the absence of cultural or religious obstacles, dietary protein sources represent the primary production cost [12,13]. Moreover, wild birds naturally consume insects as freerange poultry. Therefore, insect meal currently appears as a promising protein source in poultry nutrition [3,4,6,11,12]. In the European Union, insect-derived proteins are currently authorized for feeding fish, pets, poultry, and pigs, but not for ruminants [1]. Tenebrio molitor (TM), a species of darkling beetle, is the first insect approved by the European Food Safety Authority as a novel food [1,12]. Recent studies investigated the use of larva meal supplementation in broiler chicken feed, and TM resulted as an acceptable protein source

Experimental Animals
The experimental protocol and diet formulation are described in detail in Biasato et al., 2018 [12]. Briefly, 80 one-day-old male broiler chicks (Ross 708) were raised in 10 pens (8 birds/pen) and fed with two diets (5 pens/diet): a control diet (CT) and an isonitrogenous and isoenergetic TM diet, obtained by partially replacing soybean meal, corn gluten meal, and soybean oil with 15% levels of full-fat TM larva meal (Gaobeidian Shannong Biology Co., Ltd., Gaobeidian, China). Feed ingredients and the proximate composition of the experimental diets are reported in Table 1, and detailed in Biasato et. 2018 [12]. The growth performance of the broiler chickens was also evaluated as reported in detail by Biasato et al. [12]: the live weight and the average daily feed intake, as well as the feed conversion ratio of the birds, increased along with the increasing levels of dietary TM meal inclusion. The experimental trial lasted 53 days: during the first 3 weeks, the animals were heated by infrared lamps to maintain the suitable temperature, according to standard breeding practices. The lighting schedule was 23 h light:1 h darkness until day 3, and then 18 h light:6 h darkness until slaughter age. Tissues were sampled from 20 broilers (10 for each dietary treatment, 2 birds/pen). Intestinal tissue collection included approximately 5 cm length of the jejunum proximal tract adjacent to the Meckel's diverticulum (J), and the apex of the caecum (C). The liver parenchyma (L; central right lobe) and breast (B; central Pectoralis major) were collected on the external surface.

Sample Preparation, and RNA Extraction
In total, 80 tissue samples were collected, 10 for each diet (CT and TM)/tissue (breast, liver, jejunum, and caecum). They were immediately stored in RNAlater solution (Ambion) and then kept at −80 • C until the RNA extraction. The total RNA of each homogenized sample was extracted using the TRIzol ® method (Invitrogen, Carlsbad, CA, USA). The RNA concentrations were quantified using Qubit 2.0 Fluorometer (Invitrogen, Carlsbad, CA, USA). The integrity of the total RNA samples was evaluated using an Agilent 2100 (Agilent Technologies Inc., Santa Clara, CA, USA), and only the samples with an RNA integrity number (RIN) near 7.0 were used for RNA sequencing. Each group of 10 RNAs were pooled before library construction and the 8 pools, one for each tissue and diet, were used for RNA-seq analysis.

RNA Sequencing and Data Analysis
RNA-seq library construction and RNA high-throughput sequencing were performed at the IGA Institute (Udine, Italy). The TruSeq Stranded mRNA Sample Prep Kit (Illumina, San Diego, CA, USA) was used for library preparation following the manufacturer's instructions, starting with 1µg of RNA. Final libraries were quantified by using the Qubit 2.0 Fluorometer (Invitrogen, Carlsbad, CA, USA) and quality tested using the Agilent 2100 Bioanalyzer RNA Nano assay (Agilent Technologies, Santa Clara, CA, USA). Libraries were analyzed for RNA-seq using the Illumina HiSeq 2000 System (Illumina, San Diego, CA, USA), according to the manufacturer's specifications. Lower quality bases and adapters were removed (trimming) using Cutadapt software [19].
The expression levels of the transcripts in each sample were estimated as the number of fragments per kilobase of exon model per million reads mapped (FPKM) [17]. The following equation was used: FPKM = (C)/(N × L) ×109, where C is the uniquely mapped counts determined from the high-quality category in million, L is the length of the cDNA for the longest splice variant for a particular gene model, and N is the total mappable reads, which were determined as the sum of the high-quality reads and the highly repetitive reads [17]. The Cuffdiff option was used to perform comparisons of genes/transcripts expression between CT and TM groups in each tissue and the log2 fold change (FC) was defined [18]. Geometric normalization was applied. Transcripts with a FPKM < 5 were considered as not expressed. Transcripts with a delta (∆FPKM = absolute value of the FPKM difference between transcripts CT and TM,) below 5 were discarded.
Differentially expressed genes (DEGs) were therefore identified using a Fisher's exact test with the Bonferroni procedure [21,22] and significant differences were considered for p < 0.05; in this study, DEGs are the only differentially expressed gene with ∆FPKM > 5 and p< 0.05. The complete list of DEGs was used in functional analysis.

Functional Analysis
A differential expression pattern analysis of DEGs was performed using Reactome Knowledgebase (https://reactome.org (accessed on 31 January 2022)) considering Gallus gallus, but the information stored in the chicken database is quite small, so the analysis was carried out using the most complete human database in order to obtain more information. The complete list of differentially expressed coding genes was used as a reference set in the analysis. The functional gene ontology (GO) annotation was subsequently extracted from the database. Pathways with p-value and false discovery rate (FDR) > 0.05 were removed and the expression level of genes involved in significant biological pathways were evaluated.
In total, around 40 × 10 6 reads (excluding the adaptor sequences) were sequenced from each library and on average, 70% were aligned. The number of distinct aligned reads was similar among each library (liver 26 × 10 6 , 74%; jejunum 22 × 10 6 , 57%; ceacum 21 × 10 6 , 58%; breast 21 × 10 6 , 69%) ( Figure S1). The coverage over the gene body, defined as the ratio of the base number in a gene covered by unique mapping reads to the total base numbers of that gene, is high (more of 50.000×) and dropped dramatically only after the 95 • percentile.

DEGs Analysis
Based on the number of reads, the expression levels of the transcripts in each sample were estimated as FPKM.
Around 10,000 transcripts were shown to have a FPKM below 5 ( Figure 1); these were considered as not expressed; consequently, only 40% of the transcripts tested were expressed. Only 1% of transcripts had a FPKM > 1000, and they are listed in supplementary material.
the ratio of the base number in a gene covered by unique mapping reads to the total base numbers of that gene, is high (more of 50.000×) and dropped dramatically only after the 95° percentile.

DEGs Analysis
Based on the number of reads, the expression levels of the transcripts in each sample were estimated as FPKM.
Around 10,000 transcripts were shown to have a FPKM below 5 ( Figure 1); these were considered as not expressed; consequently, only 40% of the transcripts tested were expressed. Only 1% of transcripts had a FPKM > 1000, and they are listed in supplementary material. In most transcripts, little variation of expression was observed, and only the 25% had a ΔFPKM above 5 FPKM. The number of DEGs in which the expression was significantly (p < 0.05) modified by TM feeding ranged between 234, in liver tissue, to 308, in caecum tissue ( Table 2).  In most transcripts, little variation of expression was observed, and only the 25% had a ∆FPKM above 5 FPKM. The number of DEGs in which the expression was significantly (p < 0.05) modified by TM feeding ranged between 234, in liver tissue, to 308, in caecum tissue ( Table 2). The results in Table 2 showed that 117, 118, and 182 DEGs were upregulated and 120, 116, 168, and 126 were downregulated in breast muscle, liver, jejunum and caecum tissues, respectively, and less than half have a FC greater than 2.5 in absolute value. The fold-change and chromosome distribution of the DEGs are reported in Figure S2. Table 3 reports the type of transcripts identified in the DEGs.
In all tissues, the highest DEGs were related to small nuclear RNA (snoRNA or SNOR) and micro-RNA (miRNA) transcripts, represented almost completely by a FC > 2.5 in absolute value. The SNORs belong to the family of non-coding RNA (ncRNA). The ncRNAs, although not codifying any protein, are involved in chromatin structure regulation, replication mechanisms, DNA transcription, and epigenetic mechanisms of regulation [23,24]. For each tissue, specific SNOR were modulated by TM feeding (Figure 2).

Caecum DEGs
The DEGs with a fold change below −2.5 were 52 in the caecum. Downregulated coding genes in the caecum included avian beta-defensin 9 (AvBD9, FC = −2.6), also called gallinacin 9 (GAL9), which is a host defense peptide related to immune response against Gram-negative bacteria [25,26]. Host defense peptides are important effector molecules of the innate immune system in vertebrates [27]. Another coding gene was apolipoprotein A4 (APOAIV, FC = −2.6) associated with innate immune response in mucosa, but also with cholesterol biosynthesis and lipid transport [28]. Twenty-nine transcripts with a fold change above +2.5 were observed, including exclusively non-coding genes.

Jejunum DEGs
The DEGs with a FC lower than −2.5 were 40 in the jejunum. The downregulated coding genes were chymotrypsin-C precursor (CTRC, FC = −3.7), a member of the peptidase S1 family of genes that encode serum calcium-decreasing factor proteins with chymotrypsinlike protease activity [29]; chymotrypsin-like elastase family member 1 (CELA1, FC = −3.9), and serine-type peptidase activity (CELA2A, FC = −4.4), a subfamily of serine proteases that hydrolyze many proteins in addition to elastin [29]. The other DEGs were ncRNA or unidentified genes. Twenty-five transcripts with a fold change above +2.5 were observed but, as in the liver, we included exclusively SNOR, miRNA non-coding genes, or unidentified transcripts. The results in Table 2 showed that 117, 118, and 182 DEGs were upregulated and 120, 116, 168, and 126 were downregulated in breast muscle, liver, jejunum and caecum tissues, respectively, and less than half have a FC greater than 2.5 in absolute value. The foldchange and chromosome distribution of the DEGs are reported in Figure S2. Table 3 reports the type of transcripts identified in the DEGs. In all tissues, the highest DEGs were related to small nuclear RNA (snoRNA or SNOR) and micro-RNA (miRNA) transcripts, represented almost completely by a FC > 2.5 in absolute value. The SNORs belong to the family of non-coding RNA (ncRNA). The ncRNAs, although not codifying any protein, are involved in chromatin structure regulation, replication mechanisms, DNA transcription, and epigenetic mechanisms of regulation [23,24]. For each tissue, specific SNOR were modulated by TM feeding (Figure 2).

Breast DEGs
Transcripts with a fold change below −2.5 were 40 and only 1 coding gene was downregulated: ankyrin repeat domain 2 (ANKRD2, FC = −3.1), belongs to the conserved muscle ankyrin repeat protein (MARP) family. The expression of MARPs is induced in response to physiologic stress, injury, and hypertrophy [30,31]. The number of DEGs with an FC higher than 2.5 in breast tissue was 55, including upregulated coding genes such as HPS5 biogenesis of lysosomal organelles complex 2 subunit 2 (HPS5, FC = +4.6), which may regulate the intracellular vesicular trafficking in fibroblasts and may be involved in the regulation of the general functions of integrins, synthesis, and the function of lysosomes and of highly specialized organelles [29,32,33]. The fibrillin 1 (FBN1, FC = +2.6), proteoglycan 4 (PRG4, FC = +4.0), and avidin gene (AVD, FC = +3.45) were also upregulated. The other DEGs are ncRNA or unidentified genes.

Functional Analysis of DEGs
A functional analysis was performed to explore which functional components were regulated by TM meal inclusion. The analysis showed that only half of the DEGs identified were implicated in known biological processes. The DEGs distribution in biological network is illustrated in Figure 3. Results showed that TM inclusion mainly modulates the expression of genes related to gene expression, the immune system, signal transduction, metabolism, and proteins metabolism. Results showed that TM inclusion mainly modulates the expression of genes related to gene expression, the immune system, signal transduction, metabolism, and proteins metabolism.

Expression Pathway Analysis in Breast
Ninety-one genes, among the observed 181 DEGs, were annotated in the GO category and used for functional analysis (Table 2). Among the enriched significant biological pathways, different interactions were observed for signal transduction, particularly with homeostasis, the immune system, and gene expression. The last also includes an interaction with protein metabolism.
Striated muscle contraction (Table S1) was the most significant pathway that was influenced by TM dietary inclusion in breast tissues. Fifty-six known molecules were involved in the muscle contraction pathway, which occurs in the cytoplasm and in the cellular membrane. Five genes, coding 11 molecules in this pathway, were downregulated in the TM group (Figure 4).   Myosin light chain 1 (MYL1) expressed in fast contraction muscular fibers codifies for the alkali light chain that constitutes the myosin, a hexameric ATPase protein [34]. The troponin genes (TNN) codify for the troponin subunits; this is an important regulatory protein of striated muscle contraction, and together with tropomyosin, is located on the actin filament [34].
Groups of genes implicated in translational processes and ribosome biogenesis were downregulated in breast tissues ( Figure 5). The decrease in ribosome biogenesis and translation, essential to cell growth, might lead to decreased tissue development. Moreover, ribosomal protein large P1 (RPLP1), which was downregulated in breast tissues (FC = −0.47), but unaltered in other tissues, is related to human tumor development [35].

Expression Pathway Analysis in Gut and Liver
On the list of DEGs, 85/185 (liver), 116/223 (jejunum) and 117/135 (caecum) were annotated in the GO category. According to the functional analysis, as in breast tissue, TM inclusion led to significant dysregulation of the "signal transduction" pathway, together with different interactions concerning homeostasis, the immune system, and gene expression. The last influenced protein metabolism in particular. In the liver and caecum, the cellular cycle and DNA replication pathways were also dysregulated by TM dietary inclusion.  In all three tissues, the most significant pathway dysregulation was peptide chain elongation (Table S1). This pathway was mainly studied in humans [36], while it has not yet been described in the Gallus species. Pathways related to ribosome biogenesis and translation processes, which are essential for cell growth and tissue development, were also dysregulated by TM dietary inclusion. All these pathways were associated with the upregulation of the ribosomal protein large subunit (RPL) family gene. While RPL genes codify the structural proteins of the 80S ribosomes, the elongation factor (EEF1A and EEF2) genes codify proteins involved in protein elongation ( Figure 6).  These genes are involved in many pathways related to protein turnover, and some were significantly over-expressed. Both the RPL family factors and elongation factors were upregulated in the liver and gut of poultry TM, except for RPL39, which was downregulated (FC = −0.35) in the caecum.

Protein Metabolism Expression Pathway Analysis
Considering the protein content of TM meal, the protein metabolism pathways were taken into consideration. The analysis identified 15 molecules in breast, 22 in liver, 38 in the jejunum, and 46 in the caecum tissues involved on 1244 known pathways; the main pathways were related to ubiquitin protein, coding by the ubiquitin B gene (UBB). The expression analysis of this gene was downregulated in breast (FC = −0.8), upregulated in liver (FC = +0.6), slightly upregulated in the caecum (FC = +0.05), and was not dysregulated in the jejunum tissues (Figure 7). taken into consideration. The analysis identified 15 molecules in breast, 22 in liver, 38 in the jejunum, and 46 in the caecum tissues involved on 1244 known pathways; the main pathways were related to ubiquitin protein, coding by the ubiquitin B gene (UBB). The expression analysis of this gene was downregulated in breast (FC = −0.8), upregulated in liver (FC = +0.6), slightly upregulated in the caecum (FC = +0.05), and was not dysregulated in the jejunum tissues (Figure 7).

Discussion
Over the last decades, several efforts were made by the scientific community to propose and develop new protein sources for animal nutrition. Protein sources represent one of the major costs in animal production and their scarcity is related both to limited natural resources and high commercial competition with human nutrition [1][2][3]. The aim of this study was to understand, through the use of the RNA-seq tool, the effects of dietary TM meal inclusion on the broiler chicken's health in four different tissues (breast, liver, jejunum, and caecum).
The results showed that only 13% of the transcripts were differentially expressed. Among these, one-quarter resulted in ncRNA, one-quarter showed no annotated sequences, and only one-half resulted in coding genes. In particular, the greatest number of DEGs was observed in caecum and jejunum tissues.
In this study, the attention was directed towards the most dysregulated genes (a threshold FC = 2.5 was selected) to highlight the main effects of a diet containing high levels of TM meal.
The transcripts with the highest FC > 2.5 were ncRNA, or microRNA, involved in the processes of transcript regulation, RNA maturation, and translation; in particular, SNOR

Discussion
Over the last decades, several efforts were made by the scientific community to propose and develop new protein sources for animal nutrition. Protein sources represent one of the major costs in animal production and their scarcity is related both to limited natural resources and high commercial competition with human nutrition [1][2][3]. The aim of this study was to understand, through the use of the RNA-seq tool, the effects of dietary TM meal inclusion on the broiler chicken's health in four different tissues (breast, liver, jejunum, and caecum).
The results showed that only 13% of the transcripts were differentially expressed. Among these, one-quarter resulted in ncRNA, one-quarter showed no annotated sequences, and only one-half resulted in coding genes. In particular, the greatest number of DEGs was observed in caecum and jejunum tissues.
In this study, the attention was directed towards the most dysregulated genes (a threshold FC = 2.5 was selected) to highlight the main effects of a diet containing high levels of TM meal.
The transcripts with the highest FC > 2.5 were ncRNA, or microRNA, involved in the processes of transcript regulation, RNA maturation, and translation; in particular, SNOR carry out several roles in ribosome synthesis, translation, and the regulation of alternative splicing and oxidative stress [23,24]. It would be interesting to focus on the role of these transcripts because they could be related to the mechanisms of body adaptation to TM dietary inclusion; ncRNA may lead to the alternative splicing of genes, which are not modulated at the RNA expression level. Additionally, it would be interesting to use the proteome analysis to test this aspect.
The results showed that few known genes had a high fold-change; these were not involved in significant biological processes assessed by the pathway enrichment analysis. It is possible that these genes could be involved in biological processes that are currently unknown. These genes are related to the defense response, but also to the cholesterol biosynthetic process and lipid transport (AvBD9, alias GAL9, APOAIV), and peptidase activity (CTRC, CELA1, and CELA2A) [25][26][27][28][29]. Host defense peptides such as AvBD9, important effector molecules of the innate immune system of vertebrates [26,27], were downregulated, with a FC below −2.5 in gut tissues.
Groups of genes implicated in the translational processes and ribosome biogenesis (RPL family), which are essential for tissue growth, were upregulated in gut tissues due to TM feeding. Biasato et al., 2017 [15] reported a greater morphological development of jejunum tissues when compared to the ileum, thus representing the physiological gut devel-opment. In fact, the significant modulation of genes associated with intestinal epithelium health [37][38][39] were not highlighted in our study, as mucine or claudine gene expression was not dysregulated by TM dietary inclusion. In addition, the modulation of proteins involved in the transport of the amino acids was not influenced by dietary treatment. This is in agreement with Biasato et al. [4,6,15], Dabbou et al. [10], and Gariglio et al. [9], who did not find any significant morphological gut alterations in broiler chickens or ducks fed diets including TM or Hermetia illucens, respectively.
Acute-phase response protein, fibrillin1, proteoglycan4, and avidin genes (HPS5, FBN1, PRG4, AVD) in breast tissue had a FC > 2.5. Fibrillin-1 is an extracellular matrix glycoprotein that serves as a structural component of calcium-binding microfibrils [40]. These microfibrils provide force-bearing structural support in elastic and nonelastic connective tissues [40]. In addition, microfibrils store transforming growth factor beta (TGF-β), which is a critical growth factor [40]. Microfibrils help to regulate the availability of TGF-β [40]. The TGF-β is inactivated when stored in and activated when released from microfibrils [40]. The upregulation of fibrillin1 might therefore inactivate TGF-β, causing slow breast muscle growth. This is also corroborated by the downregulation of the RPL genes implicated in the translational processes and ribosome biogenesis. These results are in agreement with Biasato et al., 2018 [12], who reported a lower body weight and breast weight in birds fed with feed containing 15% TM. The protein metabolism analysis highlighted that the ubiquitin B gene (UBB) was downregulated in breast tissue, thus suggesting a decrease in non-lysosomal intracellular protein degradation. This could be related to a slow protein turnover.
The biological function of avidin is not known, whereas the ankyrin repeat domain 2 (ANKRD2) gene codes for a protein with chromatin binding and transcription cofactor activity, and it is also involved in muscle remodeling accompanied by the change in content of slow and fast muscle fibers [31]. The UBB, fibrillin1, RLP genes, and ANKRD2 patterns may also partially explain the above-mentioned worsening in growth and slaughtering performance of the TM birds [4,6,15]. The downregulation of genes involved in muscle contraction (troponin and myosin) in breast tissue can be related to a reduction in the free radicals [30,31,41]. In liver tissues, the non-coding genes were markedly dysregulated by TM dietary inclusion. Protein metabolism analysis highlighted that UBB was upregulated in liver tissue, whereas the expression of this gene did not change in the intestinal mucosa. This increased expression of UBB might suggest an increase in non-lysosomal intracellular protein degradation in liver tissue, and the over-expression of the RPL genes implicated in the translational processes and ribosome biogenesis suggested an increasing of protein turnover. The UBB is also involved in the maintenance of the chromatin structure, the regulation of gene expression, and the stress response [29]. Nevertheless, the dysregulation of the stress (HSP70), apoptosis (BcL2), and necrosis (NF-kB) genes was not observed in liver tissues. Therefore, this result, and the ncRNA activation, could potentially reflect mechanisms of hepatic adaptation to TM utilization.
Interestingly, the selenocysteine synthesis and selenoamino acid metabolism pathways highlighted in the Reactome analysis was linked to expression of the RPL gene. The concentration of selenium-enriched soybean, for instance, and TM is 0.18-0.4 µg/g and 54.21 ± 1.25 µg/g, respectively [42]. In particular, selenocysteine is encoded from the codon UGA, which is normally a stop codon. Nevertheless, in presence of a particular segment of mRNA, the UGA sequence is interpreted as a constitutive element. Selenocysteine was found in 25 human selenoproteins and selenoenzymes that are fundamental for the cellular processes of homeostasis and metabolic rate [43]. Both the synthesis of selenocysteine and its incorporation into a selenoprotein requires an elaborate synthetic and translational apparatus, which does not resemble the canonical enzymatic system employed for the 20 standard amino acids [42]. Thus, the expression of the RPL family genes and ncRNA transcription could be linked to the availability of selenium. Emerging evidence suggests that enzymes responsible for selenocysteine formation and the decoding of the selenocysteine UGA codon are essential for the development and health of the human organism [44]. In the absence of selenium, the mRNAs of the selenium proteins run into NMD (nonsense-mediated decay).

Conclusions
This study showed that dietary TM meal inclusion can influence metabolic pathways in poultry. More specifically, this work is the first step towards understanding how TM utilization could shape poultry health. While TM activates the protein elongation pathway required for tissue growth and development in the gut and liver, it also contains nutrients that improve gut health and reduce the mucosa inflammation. However, high levels of TM meal inclusion in poultry diets could also lead to a possible inactivation of the growth factor, along with a reduction in tissue free-radicals. No alterations in RNA expression were detected in liver that would discourage the use of the yellow mealworm in broiler chicken nutrition. In conclusion, these results confirm previous data concerning the safety of TM meal use in poultry diets.
Further investigations can be carried out to evaluate the differential expression of genes identified with a low level of expression and the role of ncRNA. The splicing variations of transcribed genes and the use of different levels of TM inclusion (5% or 10%) might also be investigated to study the modulation of the DEGs under different dietary conditions. Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/poultry1010003/s1, Figure S1: Distribution of categories reads; Figure S2: Distribution of DEGS on Chromosomes; Table S1: DEGs Pathways analysis.