Effects of Chronic Heat Stress on Lung and Growth Performance in Cherry-Valley Ducks under Different Ambient Temperatures

1. Introduction

The growth of the global population requires large amounts of protein production and puts increasing pressure on animal production, including poultry in tropical and subtropical regions, which makes a significant contribution to global food security [1,2]. The poultry industry is an important subsector of livestock production and contributes immensely to economic growth [3,4]. However, due to global warming, the sustained hyperthermia in summer has become one of the most significant challenges which are the industry facing, resulting in reduced productivity and economic losses [5]. This disastrous effect, or called heat stress, has been worldwide studied [6].

Heat stress is defined as an imbalance between the acquisition and dissipation of heat, which induces an increase in body temperature of poultry [7-9]. It is normally classified as intense or chronic heat stress (CHS), associated with short or prolonged exposure to high ambient temperatures, respectively [10]. Unfortunately, heat stress, whether intense or chronic, can harm avian health [11,12]. Similar consequence can also be observed in all domestic animals. Particularly, poultry is more susceptible to it as partly because they do not possess sweat glands and a large part of their body surface is covered with feathers. Therefore, lungs in poultry take on more responsibility heat dissipation [13,14]. When the ambient temperature is within a certain range known as the thermal comfort zone, birds can maintain body temperature with minimal effort. Typically, the optimum temperature for growing Pekin ducks is 18 °C to 20 °C [15]. Once the temperature exceeds this range, birds make behavioral, physiological, and metabolic changes to maintain body temperature and mitigate the damage caused by high temperatures [16-18]. These changes are characterized by panting, reduced food intake, and weight loss, which can lead to increased feed conversion ratios (FCRs), stunted growth, and poor meat quality. If they are unable to dissipate body heat efficiently, the ensuing multi-organ dysfunction can be fatal [19-21].

With advances in high-throughput screening technology in recent years, the transcriptome sequencing is routinely utilized to identify and quantify RNAs of different tissues and cells [22]. RNA molecules carry transcribed information encoded in genes, which can be translated into proteins or affect the expression levels of genes directly or indirectly [23]. The number of different transcripts under the temperature variation can provide insights into the cell state and stress mechanisms [22,24]. Heat shock protein (HSPs) and phosphoinositide 3-kinase (PIK3) are key players in heat stress acclimation, and the expression levels of genes encoding these two proteins were found to be significantly upregulated in heat-treated Shan Ma, Fujian shelducks, Pekin, and Muscovy ducks [25-29]. However, the differentially expressed genes (DEGs) associated with duck growth under heat stress conditions are unstudied.

The Cherry-Valley duck (CVds) is one kind of dual-purpose duck, occupying a big market share of the global poultry industry on account of its comprehensive superiority [30]. Here, we mainly investigated the effects of the CHS on the growth performance and carcass traits of the CVds and the transcriptomics response in lung. In other words, the present study aimed to explore the effects of different ambient temperatures on the growth performance of CVds and characterize the possible mechanisms responding to heat stress using Next-Generation Sequencing technology (NGS). This work should be of great importance for expanding our knowledge of the optimal temperature for CVds breeding and how high temperatures affect their growth performance. Also, it should provide some guidelines and referential value for modern broiler producers.

2. Materials and Methods

2.1. Animals, Animal Ethics, and Experimental Design

This study was carried out strictly following the regulations of the Administration of Affairs Concerning Experimental Animals (Decree No. 63 of the Jiangsu Academy of Agricultural Science on 8 July 2014). All animal experiments were approved by the Research Committee of the Jiangsu Academy of Agricultural Sciences (Nanjing, China). Refer to [15,31], a total of 36 newborn CVds with the same genetic background were purchased from a commercial farm, then averagely and randomly divided into three groups (i,e., every group has three biological samples (three individual CVds per sample)). For short, W20, W29 and W320 were referred hereinafter. In the first stage of feeding, all of the ducks were housed in one environmental control chamber and watered and fed ad libitum (Figure 1b). During the first 4 days, the room temperature was set at 35 °C, dropping to 34 °C for the next three days. During the period from 8 d to 12 d, the temperature in W29 room dropped by 1 °C per day, finally reaching 29 °C and lasting for thirty days (blue line in Figure 1a). Meanwhile, the temperature in W20 and W320 room was adjusted downward by 1 °C per day from 8 d to 21 d. After that, a constant temperature of 20 °C was set for W20 (Figure 1a, red line), while W320 (green line in Figure 1a and Figure 1c) was moved outside for bearing the big temperature change caused by the alternation of day and night. During the whole breeding process, the ducks were free to feed from the commercial standard diets in pellet form (11.27 MJ metabolizable energy/kg and 175 g crude protein/kg) and water from drip-nipple water supply line. The feed consumption was recorded every day for the calculation of feed conversion ratio as follows: FCR = feed given / animal weight gain. In simple terms, our interests were to investigate the physiological changes at CVds caused by three rearing modes only with different ambient temperature settings (Figure 1c and 1d). Finally, at day 43 of the market-age, nine samples from each group were stunned with an electrified pool (36 V) and euthanized by jugular vein dissection. The second-largest organ for heat dissipation, lungs, were carefully separated, frozen in liquid nitrogen, and stored at -80 °C before sent to the sequencing company.

2.2. Measurements of Carcass Traits

Body length, weight gain, and feed conversion ratio (FCR) were measured for each group of Cherry Valley ducks on 7 d, 16 d, 25 d, 34 d, and 43 d. Following the example of the reference [32], body length, weight gain, and feed conversion ratio were multiply measured on 7 d, 16 d, 25 d, 34 d, and 43 d. After weighing on 49 d, keel length (from the anterior to the posterior edge of the keel) and body oblique length (between the tuberosity of the shoulder joint and posterior tuberosity of the ischium) were also measured using a tape measure with an accuracy of 1 mm. In addition, in order to examine the other body changes, the depth and width of breast, the back width, the weights of liver and spleen, abdominal fat, breast muscle, leg muscle were also measured.

2.3. Microscopic Structure of Lung

The microexamination with x20 magnification were implemented for every individual duck in the three groups. The lung samples excised from the CVds in W20, W29, and W32 were fixed in 4% paraformaldehyde (PFA) for 24 h and then dehydrated with a graded alcohol series. The specimens were then infiltrated and embedded in paraffin, cut into thin seriate slices of approximately 5 μm, mounted onto glass slides, stained with hematoxylin and eosin, and observed with an Olympus microscope (BX53).

2.4. Next-Generation Sequencing and Bioinformatics Pipeline

The RNAs extracted from the lung samples from W20, W29, and W320 were isolated using TRIzol® reagent (Invitrogen, Carlsbad, CA, USA) and treated with DNase I RNase-free (TaKara) to remove genomic DNA. The concentration and integrity of RNA were measured by a NanoDrop 2100 spectrophotometer and an Agilent Bioanalyzer 2100, respectively. For each group, a NEBNext® Ultra™ Directional RNA Library Prep Kit was used to prepare an mRNA library using 3 µg of RNA per sample. The obtained mRNA libraries were sequenced on an Illumina Novaseq 6000 platform.

After the above deep sequencing, the raw paired-end reads were preprocessed using SeqPrep (https://github.com/jstjohn/SeqPrep) and Sickle (https://github.com/najoshi/sickle) by removing adaptor sequences and filtering low-quality reads. Clean reads were aligned to the high-resolution reference genome of Peking duck using HISAT [33,34]. The mapped reads were assembled by StringTie in a reference-based approach [35]. Also, the number of reads mapped to each transcript was calculated using RSEM, and the transcripts per million reads (TPM) was estimated to measure the transcript level of each gene [36]. DESeq2 was used to identify differentially expressed genes (DEGs) between W20 and W29/W320 [37]. Genes with adjusted P values <0.05 and a >1.5-fold change in expression were considered statistically differentially expressed. To understand their functions, Gene Ontology (GO) functional enrichment and KEGG pathway analysis were carried out through the Goatools and KOBAS, respectively [38,39]. The DEGs were also imported into the STRING database (version 11.5) [40]. The results obtained were imported into Cytoscape software to establish protein–protein interactions (PPI) network, and the cytoHubba plug-in was used to screen important PPI modules [41,42]. The GO and KEGG enrichment analysis of co-action targets was conducted by the BINGO plug-in and the R package clusterProfiler [43,44].

2.5. RNA Extraction, Quantitative Real-Time PCR Validation

The total RNA was extracted from lung samples using the RNA isolator Total RNA Extraction Reagent (Vazyme, Nanjing, China) according to the manufacturer’s instructions. The concentration and purity of RNA were measured with an ultra-microspectrophotometer at 260/280 nm. The integrity of RNA was estimated by agarose gel electrophoresis. RNA was further reverse transcribed into cDNA using HiScript III RT SuperMix with gDNA wiper (Vazyme, Nanjing, China) and a thermal cycler according to manufacturer’s instructions. The RT reaction was conducted in 20 μl of the reaction mixture at 37 °C for 15 min and 85 °C for 5 s, followed by cooling at 4 °C.

The two-step quantitative real-time PCR (RT-qPCR) was performed using SYBR Green Master Mix (Yeasen, Shanghai, China) and ABI 7,500 Sequence Detector (Applied Biosystem, Carlsbad, CA, USA) following the manufacturer’s recommendation. The reactions were performed with 10 μl of first-strand cDNA, 0.4 μM of forward and reverse primers and 10 μl of 2 × SYBR Green Master Mix in a final volume of 20 μl. Primers were synthesized by Sangon Biotechnology (Shanghai, China) according to the sequences described in GenBank and the reference genome of Anas platyrhynchos (GCF_015476345.1)[34,45]. The specific sequences of primers are listed in Table 2. The procedure of real-time PCR analyses was as follows: 95 °C for 2 min, followed by forty cycles of 95 °C for 10 s and 60 °C for 30 s. All of the PCR reactions were performed in triplicate corresponding to the three groups in each group. To confirm specific product amplification, a melting curve analysis was carried out. The relative gene expression was calculated using the 2−^ΔΔCT method.

2.6. Statistical Analysis

All experiments designed for different groups were independently carried out in triplicate as every group has three independent samples. The differences in the comparison of W20 with W29 and W20 with W320 were analyzed by the student’s t-test. The data are presented as the mean ± standard deviation. The significance levels were set as ^*P < 0.05 and ^**P < 0.01.

3. Results

3.1. Growth Performance and Carcass Traits

The effects of ambient temperature on the growth performance of CVds is presented in Figure 2. At the beginning, the whole-body weight and the daily intake in all of the three groups increased rapidly. The same trend was also observed in the daily intake of the ducks. As a result, at the harvest time-node of day 43, W20 had significantly higher body weight than both of W29 (P < 0.05) and W320 (P < 0.01). The average daily intake of ducks in W20 was much more than that of ducks in the other two groups since day 16 (P < 0.01). The measurement of the body size was carried out on day 43. The calculated results shows that W20 had longer keel than W29 (P < 0.05), but no statistical difference comparing with W320. Meanwhile, the body oblique of W20 was significantly longer than that of the other two groups (P < 0.01). These outcomes indicated that the growth rate of CVds would slow down when the temperature increases from 20 °C to 29 °C, which is coincide with the conclusion of [4]. On the other hand, although there was a slow upward trend in the FCR of the three groups of ducks from day 7 to day 43, there was no significant difference between them except on day 16 (data not shown). It seems that high temperature may not have much impact on the FCR of CVds.

Regarding to the other carcass traits, as the Figure 3 shows, there was no significant change between any two groups at breast depth, back width and abdominal Fat. But W20 had bigger indexes than W29 and W320 with statistical significance at breast width, liver weight, spleen weight and leg muscle weight. In addition, W20 had more the breast muscle only than W29. On the other hand, although there was a slow upward trend in the FCR of the three groups of ducks from day 7 d to day 43 d, there was no significant difference between them except on day 16 d. It seems that high temperature may not have much impact on the FCR of CVds.

3.2. Lung Appearance

The lung tissues dissected from W20 were generally reddish without oozing. On the whole, the alveoli were observed to be intact and smooth under the microscopy, without thickening of the alveolar walls but with little visible inflammatory exudate in the alveolar space. In addition, there was no congestion or edema in the pulmonary interstitial area (Figure 4). In contrast, large oozing points could be observed clearly in in the lung tissue of CVds from W29. The shape of the alveoli became elongated and irregular, which differed greatly from their normal morphology. Moreover, as the Figure 4 shows, the thickening alveolar wall, the occlusion of alveolar capillaries, the inflammatory cellular infiltration in the alveolar lumen and the edema in the pulmonary interstitial area are at the sight of on seeing. In the lung tissues in W320, it was observed that some of the alveoli were intact, while the shape of some other alveoli was close to that of W29. Furthermore, the alveolar infiltration and interstitial congestion were alleviated compared to W29 (Figure 4). Overall, W29 had the most damaged lungs, W320 had the second most damaged lungs, and W20 had the least, perhaps indicating that high temperatures cause lung lesions in CVds.

3.3. Bioinformatics Outcome

To investigate the genes associated with heat tolerance in CVds, the transcriptome analysis was performed using the Illumina Novaseq 6000 platform on three individual lung samples from W20, W29, and W320, independently. A total of 140,085,820, 152,129,548, 141,677,462, 132,793,402, 140,933,246, 118,826,918, 128,259,234, 118,574,410, and 122,680,820 clean reads with 17,575,757,205, 19,161,362,732, 18,577,850,816, 17,260,232,547, 17,798,316,784, 15,501,757,897, 16,516,406,070, 15,726,318,356, and 16,005,030,545 nucleotides were obtained from the nine cDNA libraries, respectively. More than 91 % of the clean reads from each library can be mapped to the reference genome of Anas platyrhynchos. A total of 16,074 expressed genes were detected, including 16,058 known genes and 16 predicted novel genes. In addition, a total of 1,968 novel transcripts were identified (Table 1).

The expression level of each transcript was calculated with the transcript per million (TPM) method. When studying the differences in the expression levels of genes, the expression profiles of their longest transcript isoforms were used. Those with a fold change ≥ $\frac{3}{2}$ or ≤ $\frac{2}{3}$and a q-value < 0.05 were determined as DEGs. The number of DEGs from the comparison of W20 and W29 (hereafter refer to W20/W29) was 307, while the number was 364 when comparing W20 and W320 (hereafter refer to W20/W320). To be specific, the number of up-regulated genes in W29 and W320 was respectively 163 and 178, while the number of down-regulated genes was 144 and 186, respectively (Figure 5a and Supplemental Figure S1). Of these DEGs, 76 were expressed in both groups, 231 were only expressed in W29 while the number in W320 is 288 (Figure 4b).

The Goatools software was used to identify GO terms that were remarkably enriched in DEGs [38]. The remarkably enriched GO terms in DEGs were screened out with Fisher’s exact test (Benjamini-Hochberg false discovery rate < 0.05). Unfortunately, according to the GO functional enrichment analysis, there were no significantly enriched terms comparing either in W20/W29 or W20/W320. Nevertheless, the GO functional annotations indicated that DEGs in both W20/W29 and W20/W320 were involved in three biological processes: ‘cellular response to external stimulus’, ‘response to extracellular stimulus’ and ‘cellular response to extracellular stimulus’. The genes may play important roles in responding to HS [46-48]. In the category of molecular function, only the term ‘epidermal growth factor binding’ was enriched in both W20/W29 and W20/W320. As for the category of cellular component, ‘respirasome’ as well as ‘mitochondrial respirasome’ were enriched for the W20/W29 DEGs and the W20/W320 DEGs. This result suggested that the genes involved in respiration and the ATP synthesis process, might be affected by CHS.

The significantly enriched pathways to which DEGs belong to were determined by the KEGG Orthology-Based Annotation System (KOBAS) [39]. According to the KEGG pathway analysis, no pathways were observed as significantly enriched either W20/W29 or W20/W320 (Supplemental Figure S3). Nevertheless, we found that these DEGs were associated with multiple metabolic pathways, such as ‘fatty acid biosynthesis’ [map00061], ‘fatty acid degradation’ [map00071], ‘oxidative phosphorylation’ [map00190], ‘pyruvate metabolism’ [map00620], ‘purine metabolism’ [map00230], and ‘lysine degradation’ [map00310]. This indicated that these metabolic processes may play important roles in response to heat. Additionally, ‘thermogenesis’ [map04714], ‘growth hormone synthesis, secretion and action’ [map04935], and ‘insulin secretion’ [map04911] were also enriched for DEGs in both W20/W29 and W20/W320, suggesting that the secretion of certain hormones in CVds may be altered at high temperatures, which would affect thermogenesis. CVds may be altered at high temperatures, which would affect thermogenesis [16].

The PPI network was constructed with the W20/W29 DEGs using the STRING database, as shown in Figure 7a [40]. Then, the obtained results were imported into Cystoscope software and analyzed by the cytoHubba plugin to screen out hub genes [41,42]. Based on the MCC algorithm, the top 15 genes were identified as potential hub genes, which were ENSAPLP00000000396, ENSAPLP00000006645, ENSAPLP00000006901, MYD88, NDUFA2, NDUFA5, NFKB1, NFKBIA, TRAF3, TRAF6, UQCR10, UQCRB, UQCRC1, UQCRC2, and UQCRFS1 (Supplemental Figure S4a). The GO analysis and pathway enrichment analysis showed that these DEGs were involved in oxidative phosphorylation (Supplemental Figure S5a and Supplemental Figure S6a). The same analysis was applied to the W20/W320. The PPI network of these genes is shown in Figure 7b, and the hub genes were RPL12, RPL14, RPL24, RPL27, RPL27A, RPL32, RPL34, RPL35, RPL7A, RPS16, RPS23, RPS27A, RPS27L, RPS3A, and RPS7 (Supplemental Figure S4b). The results of GO analysis and the pathway enrichment analysis suggested that these genes may affect the translation process and even gene expression due to their involvement in the biosynthesis of ribosomal large subunits (Supplemental Figure S5ad Supplemental Figure S6b).

3.4. Real-Time Quantitative Transcription-Polymerase Chain Reaction (RT-qPCR) Validation

To experimentally confirm that the DEGs obtained in this study were credible, six DEGs were selected and their expression pattern was examined via RT-qPCR. First of all, taking W20 as a standard of comparison, IGF1 expression decreased in both of W29 and W320 with statistical significance, while NR1H3 expression inversely increased. ACSL3 had increased expression only in W29, and ACSL6 and HSPA8 had decreased expression only in W320. It is notable that FGFR1 had opposite expression change in the two comparisons. Two points are summed up: (1) CHS leads to up-regulated genes and down-regulated genes; (2) diverse rearing modes result in different gene expression panel, even disparate results.

3.5. Potential Gene Markers for CHS

We found several DEGs that were potentially associated with the growth and development in CVds, including GPAM, AGPAT3, SREBF1, ACACA, ACSL3, ACSL6, NR1H3, PYGL, FGFR1, and CACNA1H. Of the ten genes, the first six were related to lipid metabolism. GPAM and AGPAT3 act on the first and second steps, respectively, in the glycerophosphate pathway, which is the main pathway for triglyceride synthesis [49]. In W320, the expression of GPAM was significantly increased. On the contrary, in W29, the expression of AGPAT3 was significantly decreased. The triglyceride synthesis could be accelerated by SREBF1, the expression of which was significantly decreased in W29 [50]. Within cells, triglycerides are the main storage form of fatty acids [51]. ACACA regulates the synthesis of long-chain fatty acids, while the ACSL gene family plays an essential role in the activation of long-chain fatty acids [52,53]. High mRNA expression of ACACA and ACSL3 was observed in W29. However, in W320 low mRNA expression of ACSL6 was observed. Similar to SREBF1, NR1H3 is a transcription factor. It had a critical role in cholesterol homeostasis and was highly expressed in W320 [54]. In addition, the expression of PYGL encoding glycogen phosphorylase was found to be significantly up-regulated in W320 [55]. This may indicate an increased energy expenditure in ducks at high temperatures since the glycogen phosphorylase catalyzes glycogenolysis. Finally, the expression of FGFR1 and CACNA1H were also remarkably up-regulated in W320. These two genes in chicken have been associated with myogenesis and body weight, respectively [56-58]. Although the mechanisms of the aforementioned genes in heat-stressed ducks need to be studied in more detail, the results of their differential expression analysis and qRT-PCR together suggested that these genes may be used as markers of HS in CVds.

4. Discussion

The negative effects in physiology caused by CHS have been widely reported for poultry and livestock industry, including but not limited to the decreases in feed intake, feed efficiency, growth performance, the production and quality of meat and egg, and survival ability [15,24,31,32,59-66]. In addition, inflammatory response, dysbiosis, reactive oxidative stress (ROS), signal reactions and energy metabolism have been further examined in a number of studies [67-75]. Specifically, the effects of increasing ambient temperatures on the health status and well-being of ducks were investigated, focusing on granulosa cells and the jejunum [25,76]. Moreover, CHS-related genes and their mediate pathways were revealed, such as alpha-Enolase (ENO1, also HSP40), 70-kDa heat shock protein (HSP70, also HSPA9, U3IP36), HSP47, HSP60, HSP90, HSP100, IRF4, IRF8, ATP-independent pathway, rapamycin and PI3K/Akt pathway [46,47,77-81]. Plus, various strategies were proposed to alleviate these adverse effects of CHS: (1) better housing, ventilation and cooling systems (e.g., little rearing system and cage rearing system) [82,83]; (2) dietary supplementation (e.g., vitamin A, vitamin C, vitamin E, Glutamine and Herbs) [84-88]; (3) feed additive (e.g., probiotics, prebiotics, polyphenols and palm oils) [89,90]; (4) others (e.g., feed restrictions and genetic selection for heat tolerance) [91-95]. In the present study, we actually exhibited three patterns of duck rearing with the focus on the effects of ambient temperature in the industrial feedlots. From this, the big advantages of using thermostatic chambers in the process of duck breeding were illustrated. It is worth noting that we verified the suitable temperature (~20-25℃) in doors for the following similar researches. On the contrary, the outcome of W320 demonstrated the natural rearing method after early-period artificial feeding does not have superiority than the two other methods. In addition, we replicated the observation of the growth performance of ducks and investigated the transcriptomics response in the lungs. To our knowledge, it is the first comprehensive study to explore the CHS-induced effects with transcriptome profiles. Based on the bioinformatics analysis, we found some suspicious CHS-related genes, which potentially enrich the heat marker list. In sum, we have paved the way to enhance the harvest of poultry and other animals such as pigs and rabbits.

Admittedly, there are some shortcomings in this study. Firstly, only morphological changes in lung tissue were studied. Future studies should include biochemical tests of lung tissue and morphological tests of skin tissue, as the latter is the largest heat dissipation system. Second, the short-read transcriptome sequencing with data analysis and the followed validation cannot uncover the full view. The DEG list obtained here contain a few unreported genes, which probably also play potential roles in the CHS-induced response of the ducks and is worthy to be characterized in the next studies. Further, the study of the mechanisms regulating the expression of these genes is also necessary, as well as the use of other omics approaches (multi-omics studies). In addition, more examinations could be added in the present study, such as the Amino Acid composition, nutritional value, protein digestibility [96].

In conclusion, like other poultry, the CVds are also sensitive to high ambient temperature and the environmental control chamber has the advantage to improve the life quality [20,97,98]. If climate change continues as predicted in most models [19,99], further studies are urgent to reveal the response mechanism and regulatory network of ducks to CHS.

Number of Pages:31

Number of Figures: 7

Number of Tables: 2