The Sequence Read Archive (SRA) repository of the National Center for Biotechnology Information (https://www.ncbi.nlm.nih.gov/sra) was screened precisely to find the appropriate RNA-seq datasets that address the research question of the current work using the keywords “Gallus gallus,” “chicken,” “liver” and “acute heat stress”. We could find only four RNA-seq datasets in which the AHS was the main treatment. Detailed information of the four selected datasets can be found in Table 1. In addition, the accession numbers of the used samples are reported in In Supplementary Table S1 .

After retrieving the RNA-seq datasets, fastq-dump tool of SRA-Toolkit version 2.9.6 (Staff, 2011) was employed to convert the SRA files into the FASTQ format. The data quality was assessed using the fastQC tool version 0.12.1 (Andrews, 2010) (available at https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). The low quality reads were eliminated using the Trimmomatic (version 0.39) software (Bolger et al., 2014) (available at http://www.usadellab.org/cms/?page=trimmomatic) with ILLUMINACLIP, SLIDING WINDOW (3-5: 20-28), CROP (3-10), AVGQUAL (20-25) and MINLEN (40-45) options. The implemented values varied among the datasets according to their quality metrics. The trimmed reads were mapped onto the reference genome Gallus_gallus.GRCg6a (https://asia.ensembl.org/Gallus_gallus/info/index) using HISAT2 (version 2.2.1) software (Kim et al., 2019) (available at https://daehwankimlab.github.io/hisat2). The expression count matrix was generated using HTSeq-count (version 0.9.1) (Anders et al., 2015). Then, the DESeq2 package (version 3.16) (Love et al., 2014) was used to identify the DEGs using the default parameters. To characterize the gene symbols, biotypes, and positions of each transcript, we uploaded the DESeq2 results and the gtf file (http://ftp.ensembl.org/pub/release-104/gtf/gallus_gallus/Gallus_gallus.GRCg6a.104.gtf.gz) to galaxy (available at https://usegalaxy.eu/) and we used the Annotate DESeq2/DEXSeq output tables (version 1.1.0) for annotation.

The metaRNASeq package version 1.0.2 (Marot and Rau, 2013) in R was used to identify the meta-genes. First, raw p-values and log2-fold change values of all expressed genes of all four datasets were gathered in a new file. Then, p-values were combined via the fishcomb function. Comparing the sign (positive or negative) of the log2-fold change of the genes across the four datasets revealed the consistency or inconsistency of the expression of genes in all datasets. Finally, genes with consistent expression in all datasets, with adjusted p-value ≤0.05 in at least one dataset, and meta-analysis p-value ≤0.05 were considered as meta-genes.

Meta-genes were interpreted based on Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathways by DAVID web base software (https://david.ncifcrf.gov). Terms with p-value ≤0.05, FDR ≤0.2, and a fold enrichment >2 were considered as significant. Furthermore, protein-protein interaction (PPI) network was created using the STRING database (https://string-db.org/) to construct the network and identify the hub-genes. Cytoscape software (version 3.7.2) (Shannon et al., 2003) was used to visualize the retrieved networks.

The R package Weighted Gene Co-expression Network Analysis (WGCNA) (version 1.71) (Langfelder and Horvath, 2008) was used to identify the co-expressed networks, modules, and hub-genes. To alleviate the influence of noise when calculating the correlations based on the read counts, we filtered out the genes with read counts less than 10 in more than 90% of the samples, in order to reduce the sampling differences. The variance-stabilizing transformations (VSTs) output from Deseq2 were also used as input. However, in the case of using VST data, the variations caused by batch effects or other covariates could not be accounted for. Therefore, the function “removeBatchEffect” was used to eliminate the batch variations by package limma (version 3.50.3). In the first step, excessive missing values and outlier samples were examined with the options “goodSamplesGenes” and “hclust”. Soft threshold power of 0.9 was selected based on the scale-free topology index (R ²) (Zhang and Horvath, 2005) and, thus, the “pickSoftThreshold” option was used to calculate the adjacency matrix. The Pearson correlation coefficients between each pair of genes were calculated. The adjacency matrix was converted to a Topological Overlap Matrix (TOM) to minimize the effects of noise. Then, the corresponding dissimilarity matrix (1-TOM) was generated. A dynamic tree cut (DTC) algorithm was used to detect and construct the gene co-expression modules with the following parameters; cut height of 0.975, minimum module size of 30 genes, DeepSplit of 2, and hybrid method.

To identify the hub-genes, the module eigengenes were calculated using the “moduleEigengenes” function, and the “intramodularConnectivity” process was used to calculate both the intramodular connectivity (k_within) and the total connectivity (k_total). Then, the “corPvalueStudent” process was used for the identification of hub-genes based on the p-values (Degli Esposti et al., 2019). In addition, the “chooseTopHubInEachModule” and “chooseOneHubInEachModule” options were used to identify the hub-genes within each module.

After identifying the meta-genes and significant modules, the common genes between them were identified using venn diagram (http://bioinformatics.psb.ugent.be/webtools/Venn/). The common significant meta-genes were interpreted by enrichment analysis methods based on the KEGG pathway and GO using the ClueGO plugin of Cytoscape software (version 3.7.2) (Bindea et al., 2009), and subsequently, only pathways with Bonferroni Step Down corrected p-values <0.01 were considered as significantly enriched.

The identified meta-genes with their corresponding expression values were submitted to Machine-Learning for the selection of the most important HS-related genes. Additionally, the chicken breed, age at sample collection, and duration of the HS were used as three attributes. This was accomplished based on seven weighting algorithms, including Uncertainty, Chi-Squared, Relief, Gini index, Rule, Information Gain, and Gain Ratio. Meta-genes with average weighting values above 0.7 across all seven algorithms were analyzed via Rapid Miner software (version 9.9) for 10-fold cross-validation using stratified sampling with Decision Tree (Accuracy, Gain Ratio, Gini Index, and Information Gain criterion), Random Forest (Accuracy and Gain Ratio criterion), Deep Learning (Tanh and Rectifier criterion), and Naive Bayes.

For SNP calling, the Galaxy platform (available at https://usegalaxy.eu/) was used with Genotype-variants outline (available at https://github.com/cfarkas/Genotype-variants). After sorting the BAM files with samtools sort (galaxy version 2.0.4) (https://samtools.github.io/hts-specs/), we used FreeBayes (galaxy version 1.3.6+ galaxy0) (https://github.com/ekg/freebayes) for variant calling. We then used VCFfilter (galaxy version 1.0.0_rc3+galaxy3) (https://github.com/ekg/vcflib) to keep only the SNPs with a depth of more than 25 reads and a quality of more than 30. The transcript IDs harboring each variant were retrieved using VCFannotate (galaxy version 1.0.0_rc1+galaxy0) (GitHub. https://github.com/ekg/vcflib). For expression based genome wide association study (eGWAS), we used 22,197 eSNPs located on all analyzed transcripts. The expression values (read counts) of the transcripts were required to be above 10 in at least 10% of the samples to be included in the eGWAS. The quality assessment of the SNP genotypes was carried out by PLINK software (version 1.9) (Purcell et al., 2007). The nine significant modules that were output by WGCNA software were considered as nine traits, and the eigenvalues of the 32 samples for each of the significant modules were considered as dependent variables (phenotypes). A linear model was employed in the “assoc” analysis of PLINK. The significant eSNPs were visualized using an R package called CMplot (version 4.2.0) (Yin et al., 2021). Figure 1 illustrates the flowchart of the analysis steps of the four RNA-seq datasets in the present study. It should be mentioned that, all gene feature were reported in the current manuscript and supplementary files were according to the Gallus_gallus.GRCg6a genome assembly of domestic chicken.

Two additional RNA-seq datasets related to chicken livers under chronic heat stress (CHS) were used to validate the identified 77 meta-genes. The mentioned datasets were downloaded from the NCBI website with accession numbers ERP014602 and SRP100368 with 8 and 16 samples, respectively. Half of the samples in each of the mentioned datasets were from the chickens treated with 35 degrees of Celsius for 7-8 h and for 7 days, whereas the next halves were considered as control groups. The accession numbers of the used runs for the validation analysis can be found in Supplementary Table S2. The analysis pipeline was as same as mentioned above for the analysis of the AHS datasets.

We examined four selected RNA-seq datasets and identified 41, 115, 39, and 12 DEGs for SRP268422-A, SRP268422-B, ERP014602-A, and ERP014602-B datasets, respectively. Of which, 24, 83, 10, and 5 were downregulated while 17, 32, 29, and 7 were upregulated DEGs, respectively. In Supplementary File S1, detailed information about the results of the four RNA-seq datasets analyses are provided. The primary aim of analyzing four similar, related datasets was to discover the DEGs that were commonly present in all of the studied datasets. Therefore, we were most interested in the common DEGs across the results of four individual datasets. We, found, however, a little number of common DEGs. In Supplementary Figure S1, the Venn diagram of the DEGs identified in the four individual experiments has been shown.

As mentioned above, the results obtained from the four RNA-seq datasets did not support each other. Therefore, the implementation of meta-analysis seemed necessary. As such, we combined the results of the four abovementioned datasets and carried out a meta-analysis. Finally, a total of 77 significant meta-genes were identified. Ten out of the 77 meta-genes did not have gene symbols, while the remaining 67 meta-genes possessed known gene symbols. Supplementary File S2 includes the detailed results of the significant meta-genes from the meta-analysis of four RNA-seq datasets.

We used the DAVID database to perform KEGG pathway and GO analyses on meta-genes to characterize their molecular function (MF), biological process (BP), and cellular component (CC). The results revealed 3, 1, and 3 significant GO terms for BPs, CCs, and MFs, respectively. The terms “posttranslational protein targeting to membrane translocation,” “endoplasmic reticulum unfolded protein response” and “ubiquitin-dependent ERAD pathway” were three significant BP terms. The “endoplasmic reticulum membrane” CC and there MF terms “identical protein binding,” “ribosome binding,” and “misfolded protein binding” were significantly enriched by meta-genes. Through this process, four significant KEGG pathways including “Protein processing in endoplasmic reticulum,” “Protein export,” “Biosynthesis of nucleotide sugars,” and “Amino sugar and nucleotide sugar metabolism” were revealed to associate with AHS. A thorough information on the KEGG pathways and GO terms have been provided in Table 2. The PPI network of the meta-genes was also studied and the resulted network was visualized using the Cytoscape. In Figure 2 the visualized PPI network of the meta-genes is shown.

The expression values of 7829 genes were used as input data for WGCNA (Supplementary File S3). In the first step, excessive missing values and outlier samples were examined, and outlier samples were removed ( Supplementary Figure S2A ), and a value of 6 was determined as a soft threshold power (Supplementary Figure S2B). We used the dynamic tree cut algorithm to construct the modules. Input genes were grouped into 22 modules, ranging in size from 49 to 492 genes, along with 3940 genes that could not be grouped in any of the modules and therefore were called as unassigned (Figures 3–A). The hierarchical clustering of genes using the topological overlap matrix (TOM) is shown in Figures 3–B. Nine out of the 22 modules including green (r = −0.75, p-value = 8e-7), pink (r = −0.43, p-value = 0.01), royal blue (r = −0.43, p-value = 0.01), light cyan (r = −0.48, p-value = 0.006), purple (r = +0.36, p-value = 0.04), brown (r = +0.5, p-value = 0.003), turquoise (r = +0.37, p-value = 0.04), tan (r = +0.42, p-value = 0.02), and yellow (r = +0.44, p-value = 0.01) modules were identified as particularly significant. In Figure 4 and (Supplementary File S4) the genes counts and names within each of the 22 modules, and in Figure 5 the Pearson correlation coefficient and p-values of the identified modules are reported. Thereafter, hub-genes were identified per module based on the p-values using the “corPvalueStudent” (Supplementary File S5). In Table 3 the identified five hub-genes per module are reported.

Meta-analysis relieved the sample size limitation, while the WGCNA identified genes with high expression correlation. For functional enrichment analysis, we assigned significant common genes between the modules and meta-genes (Supplementary File S6 ) using the ClueGO plugin of Cytoscape software. As a result of pathway enrichment analysis, “Protein processing in endoplasmic reticulum,” and “Protein export” were found to be significant KEEG pathways, and “rough endoplasmic reticulum” and “rough endoplasmic reticulum membrane” were two significant CC terms. “Protein transmembrane transport,” “intracellular protein transmembrane transport,” “response to endoplasmic reticulum stress,” “response to unfolded protein,” “endoplasmic reticulum to cytosol transport,” “endoplasmic reticulum unfolded protein response,” “negative regulation of response to endoplasmic reticulum stress,” “regulation of response to endoplasmic reticulum stress,” “cellular response to unfolded protein,” “negative regulation of cellular protein catabolic process,” “protein exit from endoplasmic reticulum,” “ERAD pathway,” “retrograde protein transport, ER to cytosol,” “negative regulation of proteolysis involved in cellular protein catabolic process,” “regulation of ERAD pathway,” “negative regulation of ERAD pathway,” “negative regulation of proteasomal protein catabolic process” and “ubiquitin-dependent ERAD pathway” were significant terms for BP category. The connections among the terms are illustrated in Figure 6. In addition, the detailed information about the significantly enriched terms are reported in Supplementary File S7.

Seventeen meta-genes showed an average weighting values greater than 0.7. For breed, age at sampling, and duration of HS, however, the average weighting values were less than 0.7. In Supplementary Table S3, the weights of the 17 meta-genes obtained from the 7 different machine learning algorithms are reported. Ten models were applied to the dataset. The performances of the models are reported in Table 4. According to the cross validation results, the Decision Tree with the Gain Ratio criterion, the Random Forest with the Accuracy criterion, and Deep Learning with the Maxout parameter showed acceptably high accuracy. The Decision Tree identified key meta-genes that classify the chickens under AHS or control conditions based on the expression values. Here, three genes including EIF5B, HSPA5, and SEC23B are introduced as biomarkers for AHS as they were included in all Decision Trees with the Gain Ratio, Accuracy, Gini Index, and Information Gain criteria. The accuracy range of the mentioned models were 67.50%+/-37.26%, 57.50%+/-40.64%, 65.00%+/-43.23%, and 55.00%+/-42.61%, respectively. Thus, based on the Gain Ratio, if the expression value of EIF5B was greater than 6.164, the samples would fell into the AHS, otherwise, the expression value of HSPA5 must be taken into account (Figure 7; Supplementary Figure S3). The important role of these three genes on AHS was also identified earlier in the WGCNA analysis, as all of them included in the significant green module as hub-genes. The rediscovery of the mentioned three genes by the Decision Tree models confirms the important association of them with AHS.

To identify significant SNPs associated with nine significant modules, eGWAS was conducted. As a result, we identified 406 significant eSNPs (p-values ≤0.05). It is noteworthy to mention that most of the significant eSNPs were associated with two modules, including yellow and green. In Supplementary File S8 significant levels of all 406 eSNPs on all nine significant modules are provided. A total of 134 significant eSNPs located on the genes that included within the significant modules, and 52 of them were significantly associated with all nine significant modules. In Figure 8, a circular Manhattan plot illustrating the association of eSNPs with nine modules (layers of circles) is shown. We identified one eSNP on gene SERP1 (position chr9:23805255) and four eSNPs on gene SESN1 (positions chr3:67138413, 67135835, 67137358, 67137810). The SERP1 was a meta-genes which consistently downregulated in AHS, while SESN1 was a meta-genes which consistently upregulated in AHS as compared with the control group.

To summarize our results from the five employed approaches and make relevant conclusions, we decided to focus on genes that were commonly observed as meta-genes and included within the significant terms, pathways, and modules, especially those that were identified as hub-genes. Furthermore, the meta-genes with significant eSNPs that highlighted by Machine Learning were considered as more important. In Supplementary Figure S4 and Supplementary File S9, the resulting Venn diagram has been reported. Based on the summary of our results, we introduce important genes associated with AHS including; HSPA5, SSR1, SDF2L1, SEC23B, SERP1, and SESN1.

Fifteen and two meta-genes out of the 77 meta-genes were identified as differentially expressed (p-value <0.05) in the two additional datasets (i.e., ERP014602 and SRP100368, respectively). Surprisingly, only one of the 17 genes showed no match with the results of the current work. In other words, the log2-fold change of the 16 meta-genes in all four AHS datasets were as in the same direction as in the validation datasets (CHS). In Supplementary File S10 the log2-fold change of the 77 meta-genes in two validation datasets are provided. In Figure 9, the expression concordance of the 17 identified meta-genes between the AHS datasets and CHS datasets were shown. We also compared the expression pattern of the six most associated meta-genes we introduced above (i.e., HSPA5, SSR1, SDF2L1, SEC23B, SERP1, and SESN1). The expression pattern of five of them were matched between the AHS and CHS datasets.