Selenium concentration in blood samples
Chicks were supplemented with selenium-deficient diet caused decrease selenium concentration in blood samples compared with chicks supplemented with basal diet (p < 0.05) (Fig. 2)
Transcriptome Data Generation And Processing
High-quality data of 3.35 and 2.15 Gb were generated for the control and selenium deficient groups on Illumina Next Seq500 using 2 x 150 bp chemistry. The high-quality reads of each of the samples were mapped to the reference genome of Gallus gallus GRCg6a using Tophat at default parameters. The read mapping percentage were found to be 75.6% and 77.6% for the control and selenium deficient groups respectively (Supplementary Table S1 & Table S2).
Identification Of Differential Expressed Genes (Degs)
Individual transcriptomes were assembled and then differential expressed genes (DEGs) were identified using cufflinks and cuffdiff software’s. DEGs were analyzed in oviduct tissues between selenium-deficient and control groups. Significant DEGs were identified based on the statistical significance (p ≤ 0.05) and ≥ 2 log2 fold changes (FC) for Up-regulated, and ≤ -2 log2 FC values for Down-regulated genes. In oviduct samples, a total of 10,266 genes were differentially expressed, of which significant DEGs were 213 (Up-regulated) and 230 (Down-regulated). A complete list of Up and Down regulated DEGs is provided in Supplementary Table S3 & Table S4. The DEGs results are summarized in the Supplementary Table S5. In addition, the raw data of DEGs was provided in the Supplementary Table S6-S10.
Figure 3a showed that Heat maps of the top 50 DEGs, where red and green colours represent up-regulated and down-regulated DEGs, respectively in selenium-deficient group compared with control group. In addition, scatter plot (Fig. 3b) and volcano plot show significant up-regulated and down-regulated DEGs in oviduct samples (Fig. 3c).
Go And Keegs Pathway
g: Profiler is a public web server (http://biit.cs.ut.ee/gprofiler/) was used to gain insight into various GO (Gene Ontology) terms of the up and down regulated genes in chickens. The official gene symbol of the upregulated genes was uploaded to the functional annotation tool in the g: Profiler web server and Gallus gallus was selected as the reference genome16. Of the 210 genes uploaded, 194 genes were annotated into three GO terms: biological process (BP), cellular component (CC), and molecular function (MF) (Fig. 4a). Like-wise, 230 down-regulated genes uploaded, 146 genes were annotated into three GO terms: BP, CC, and MF (Fig. 4b).
In the GO analysis, a total of up-regulated 213 DEGs were annotated to 90 GO terms, covering 30 biological processes, 2 molecular functions, and 30 cellular components. The GO terms for biological processes were “cellular component movement”, “regulation of cell motility”, regulation of “signaling”, “regulation of cell migration”, “positive regulation of chemotaxis”, “regulation of signal transduction”, “regulation of locomotion”, “cell communication”. Cellular components GO terms included “endoplasmic reticulum” (GO:0005783), “endomembrane system”, “membrane”, and “vesicle”. Molecular functions enriched with “protein binding” and “protein-containing complex binding” was shown in Table 1.
Table 1
Gene Ontology (GO) terms for up-regulated differentially expressed genes (DEGs). Note: ‘BP’, ‘CC’, and ‘MF’ represent biological process, cellular component and molecular function, respectively. Adjusted P-Value.
Category
|
GO ID
|
GO Term
|
Adjusted P-value
|
BP
|
GO:2000145
|
Regulation of cell motility
|
5.50E-05
|
BP
|
GO:0030334
|
Regulation of cell migration
|
8.49E-05
|
BP
|
GO:0050921
|
Positive regulation of chemotaxis
|
0.000132202
|
BP
|
GO:0009966
|
Regulation of signal transduction
|
0.032977079
|
BP
|
GO:0010646
|
Regulation of cell communication
|
0.01941158
|
BP
|
GO:0023051
|
Regulation of signalling
|
0.021871511
|
BP
|
GO:0040012
|
Regulation of locomotion
|
0.001326849
|
BP
|
GO:0051270
|
Regulation of cellular component movement
|
0.029940223
|
CC
|
GO:0005783
|
Endoplasmic reticulum
|
0.000678557
|
CC
|
GO:0012505
|
Endomembrane system
|
0.003640007
|
CC
|
GO:0016020
|
Membrane
|
0.005992789
|
CC
|
GO:0031982
|
Vesicle
|
0.01733706
|
MF
|
GO:0005515
|
Protein binding
|
0.011364633
|
MF
|
GO:0044877
|
Protein-containing complex binding
|
0.042467714
|
On other hand, a total of 237 down-regulated DEGs between the control and selenium deficient groups were assigned into 90 enriched GO terms consisting of 30 biological processes (BP), 30 molecular functions (MF) and 30 cellular components (CC). The significant GO terms for the biological process (BP) were, GO:0007166 (cell surface receptor signalling pathway), GO:0007165 (signal transduction), GO:0007155 (cell adhesion), GO:0007154 (cell communication), GO:0007049 (cell cycle). The most relevant terms for molecular function were GO:0005488 (binding), GO:0030695 (GTPase regulator activity), GO:0022890 (inorganic cation transmembrane transporter activity), and GO:0030246 (carbohydrate binding). Under cellular component GO terms was, GO:0032432 (actin bundle filaments), GO:0015629 (actin cytoskeleton), and GO:0030670 (phagocytic vesicle membrane) (Table 2).
Table 2
Gene Ontology (GO) terms for down-regulated differentially expressed genes (DEGs). Note: ‘BP’, ‘CC’, and ‘MF’ represent biological process, cellular component and molecular function, respectively. Adjusted P-Value.
Category
|
GO ID
|
GO Term
|
Adjusted P-value
|
BP
|
GO:0007166
|
Cell surface receptor signalling pathway
|
1
|
BP
|
GO:0007165
|
Signal transduction
|
1
|
BP
|
GO:0007155
|
Cell adhesion
|
1
|
BP
|
GO:0007154
|
Cell communication
|
1
|
BP
|
GO:0007049
|
Cell cycle
|
1
|
BP
|
GO:0007259
|
Receptor signalling pathway via JAK-STAT
|
1
|
BP
|
GO:0040012
|
cell-cell signalling
|
1
|
BP
|
GO:0007267
|
Ras protein signal transduction
|
1
|
CC
|
GO:0032432
|
Actin bundle filaments
|
0.140050699
|
CC
|
GO:0015629
|
Actin cytoskeleton
|
0.475423994
|
CC
|
GO:0030670
|
Phagocytic vesicle membrane
|
0.581274167
|
MF
|
GO:0005488
|
Binding
|
0.201762502
|
MF
|
GO:0031267
|
GTPase regulator activity
|
1
|
MF
|
GO:0022890
|
Inorganic cation transmembrane transporter activity
|
1
|
MF
|
GO:0030246
|
Carbohydrate binding
|
1
|
We also analyzed the pathways enrichment of up and down regulated genes of DEGs in chickens using Kyoto Encyclopedia of Genes and Genomes (KEGGs) database. The 10,266 DEGs identified in the control versus selenium-deficient groups were classified into 23 functional pathway categories in KEGG database. These pathways were divided into five specific pathways, namely, metabolism (A), genetic information processing (B), environmental information processing (C), cellular processes (D), and organismal system (E) as shown in Table 3. In the environmental information processing, the “signal transduction” is the most abundant subcategory (gene count: 1,280). While in cellular processes, the predominant group was “transport and catabolism” (gene count: 482). In genetic information processing category, “translation” was largest subgroup and its gene count was 312. In the metabolism category, there were two categories much more common than others, called as “lipid metabolism” (gene count: 260) and “carbohydrate metabolism” (gene count: 208). Finally, in the organismal systems, there the most abundant group called “environmental adaption” (gene count: 222). Moreover, the KEGGs annotation statistics was provided in the Supplementary Table S11.
Table 3
KEGG Pathway Classification
Pathways
|
Control Vs Selenium-Deficient groups
|
|
Gene Counts
|
A. Metabolism
|
|
Carbohydrate metabolism
|
208
|
Energy metabolism
|
147
|
Lipid metabolism
|
260
|
Nucleotide metabolism
|
113
|
Amino acid metabolism
|
185
|
Metabolism of other amino acids
|
67
|
Glycan biosynthesis and metabolism
|
184
|
Metabolism of cofactors and vitamins
|
135
|
Metabolism of terpenoids and polyketides
|
27
|
Biosynthesis of other secondary metabolites
|
28
|
Xenobiotics biodegradation and metabolism
|
69
|
B. Genetic Information Processing
|
|
Transcription
|
119
|
Translation
|
312
|
Folding, sorting and degradation
|
309
|
Replication and repair
|
113
|
C. Environmental Information Processing
|
|
Membrane transport
|
30
|
Signal transduction
|
1,280
|
Signalling molecules and interaction
|
599
|
D. Cellular Processes
|
|
Transport and catabolism
|
482
|
Cell growth and death
|
411
|
Cellular community - eukaryotes
|
366
|
Cell motility
|
129
|
E. Organismal Systems
|
|
Environmental adaptation
|
220
|
Table 4. Composition of Basal Diet
Ingredients
|
Percentage
|
|
|
Corn
|
63
|
Soybean Meal
|
21
|
Limestone
|
7
|
Bran
|
5
|
Vitamin-mineral premixa
|
4
|
|
|
Nutrient level
|
|
|
|
Crude Protein %
|
16.5
|
Crude Fibre %
|
3.35
|
Lysine %
|
0.7
|
Methionine %
|
0.4
|
Cystine %
|
0.3
|
Digestible energy (ME)/kJ.kg-1
|
11.4
|
Phosphorus %
|
0.7
|
Calcium %
|
3.5
|
Selenium per mg. Kg-1
|
0.012
|
a Vitamin-mineral premix (per Kg Diet): 26,300 IU vitamin A, 8000 IU vitamin D3, 370 IU vitamin E, 40 mg vitamin K3, 35 mg vitamin B1, 100 mg vitamin B2, 30 mg vitamin B6, 0.6 mg vitamin B12, 50 mg niacin, 12 mg pantothenic acid, 13 mg folic acid, 0.8 mg biotin, 110 mg choline chloride, 50 mg Fe, 1.5 mg Cu, 50 mg Mn, 12 Zn, 0.012 mg Se, 0.095 mg I.
Table 5. Composition of Selenium-Deficient Diet
Ingredients
|
Percentage
|
|
|
Corn
|
64
|
Soybean Meal
|
21
|
Limestone
|
7
|
Bran
|
5
|
Vitamin-mineral premixa
|
3
|
|
|
Nutrient level
|
|
|
|
CP %
|
16
|
Crude Fibre %
|
4
|
Lysine %
|
0.8
|
Methionine %
|
0.4
|
Cystine %
|
0.5
|
Digestible energy (ME)/kJ.kg-1
|
12
|
Phosphorus %
|
0.6
|
Calcium %
|
3.5
|
Selenium per mg. Kg-1
|
0.002
|
a Vitamin-mineral premix (per Kg Diet): 24,500 IU vitamin A, 6000 IU vitamin D3, 570 IU vitamin E, 30 mg vitamin K3, 45 mg vitamin B1, 111 mg vitamin B2, 30 mg vitamin B6, 0.5 mg vitamin B12, 40 mg niacin, 15 mg pantothenic acid, 12 mg folic acid, 0.6 mg biotin, 100 mg choline chloride, 30 mg Fe, 1.6 mg Cu, 30 mg Mn, 15 Zn, 0.002 mg Se, 0.090 mg I.
Ppi Network And Module Analysis
The STRING database analysis depicted that 324 of DEGs PPI network comprised of 324 nodes connected with 439 different edges after applying relevant filters. Expected number of edges was observed to be 305 while the average node degree score was found to be 2.71 i.e., one node had at least 2.71 interacting nodes. Average local clustering coefficient was predicted to be 0.373 and PPI enrichment value was observed as 3.56e-13 (Fig. 5a).
Next, we imported all 324 DEGs into the STRING database loaded the STRING network in Cytoscape App software. We found that yellow color node represents up-regulated DEGs, whereas the green color node denotes down-regulated DEGs (Fig. 5b). Furthermore, we conducted functional module analysis using the MCODE plug-in in Cytoscape software (Supplementary Table S12-S20). The results showed that a total of 8 functional modules are shown in Fig. 6a-h. Among these, module 1 included 7 nodes and 21 edges with highest score of 7 (Fig. 6a). Module 2 included 6 nodes and 12 edges with a score of 4.8 (Fig. 6b). Module 3 comprised 19 nodes and 35 edges with a score of 3.889 was shown in Fig. 6c). The genes in Module 4 consists of ANKRD2, CSRP2, and FHL2 with 3 nodes and 3 edges with a score 3 (Fig. 6d). Module 5 (SEC16B, ETV5, and TMEM18) consists of score 3 (3 nodes and 3 edges) was shown in Fig. 6e. Module 6 has ERO1L, CLGN, and TXNDC5 (3 nodes, and 3 edges; score: 3) (Fig. 6f), Module 7 (CLDN10, CLDN1, and OCLN) has 3 nodes and 3 edges with score of 3 and 8 contained 3 nodes and 3 edges with a score 3 (Fig. 6g and Fig. 6h).
Hub Gene Identification And Enrichment Analysis
According to the Degree algorithm ranking in the CytoHubba, the top ten highest scored genes were selected as the hub genes (Fig. 7a), including FRK, JUN, PTPRC, ACTA2, MST1R, SDC4, SDC1, CXCL12, MX1, and EZR was shown in Fig. 7b. Furthermore, ten hub genes were submitted to CancerGeneNet network analysis resulted in three types of graphs displayed. Figure 8a showed that direct connection between query hub genes (Level 1), Fig. 8b revealed that connecting hub cancer genes and also interacting with two cancer genes, all to form first neighbours (Level 2). The protein-to-protein (PPI) that pointed to external stimulus such as angiogenesis, metastasis, apoptosis, monocyte differentiation, cell migration, brown adipogenesis, and proliferation reported in Fig. 8c (Level 3). Additionally, the enrichment of hub genes involved in cancer pathways were provided in Table 6.
Table 6
Gene enrichment analysis of hub genes
Cancer Module
|
Hits
|
P-values
|
AML1-ETO in AML
|
4/18
|
< 1.0E-4
|
RTKs in cancer
|
1/12
|
< 1.0E-4
|
ASXL1 in AML
|
1/13
|
< 1.0E-4
|
FLT3 in AML
|
1/40
|
< 1.0E-4
|
Triple mutant AML
|
1/18
|
< 1.0E-4
|
KIT in AML
|
1/30
|
< 1.0E-4
|
AML_TRIPLETS
|
1/22
|
< 1.0E-4
|
Sara
|
1/67
|
< 1.0E-4
|
BCR-ABL in AML
|
1/20
|
< 1.0E-4
|
ErbB receptors in Cancer
|
1/12
|
< 1.0E-4
|
WNT/FLT3
|
3/27
|
< 1.0E-4
|
mTOR in Cancer
|
1/15
|
< 1.0E-4
|
Onco-fusion protein in AML
|
1/28
|
< 1.0E-4
|
MLL fusion protein in AML
|
1/21
|
< 1.0E-4
|
miRNA in AML
|
3/35
|
< 1.0E-4
|
Wnt in Cancer
|
1/10
|
< 1.0E-4
|
FLT3-ITD in AML
|
3/36
|
< 1.0E-4
|
Haematopoiesis Transcription Control
|
2/16
|
< 1.0E-4
|
NPM1_ new
|
1/25
|
< 1.0E-4
|