Stocking density affects transcriptome changes in the hypothalamic-pituitary-gonadal axis and reproductive performance in ducks

Abstract The hypothalamic-pituitary-adrenal (HPA) axis and hypothalamic-pituitary-gonadal (HPG) axis plays a central role in mediating physiological responses related to the reproductive system under any stressful condition. However, the molecular mechanism underlying the effects of stress on physiology still needs to be elucidated. This study demonstrated that increasing the stocking density from 4 to 8 birds/m2 during the laying period decreased the egg production rate of laying ducks by 13.04 − 63.55% and feed intake by 7.40 − 23.44%. Transcriptome analysis between high- and low-feeding-density laying ducks revealed 469, 509, 428 and 210 differentially expressed genes (DEGs) in the hypothalamus, pituitary, ovary and follicular membrane, respectively. Gene ontology (GO) and KEGG enrichment analysis showed that the DEGs in the hypothalamus and pituitary were primarily enriched in the biostimulation and dopamine secretion pathways. The major enrichment pathways in the ovarian and follicular membranes involved lipid metabolism, negative regulation of inflammatory response, and steroid hormone biosynthesis. Among the DEGs in the HPG system, POMC and GnRH1 were identified, which may be manifesting their crucial roles in regulating the stress response and reproduction. Our data showed that a high stocking density as environmental stress negatively affects the reproductive performance in ducks through transcriptional changes in the HPG axis. Highlights Raising the stocking density from 4 to 8 birds/m2 decreased the egg production rate and feed intake in laying ducks. The transcriptome indicated that stocking density affects the stress response of laying ducks through the hypothalamus-pituitary. The stress signals are subsequently transmitted to affect gene expression related to the reproduction process in laying ducks’ ovarian and follicular tissues. The hypothalamic expression of POMC and GnRH1 may play a central role in integrating stress signals and the reproductive processes of laying ducks.


Introduction
Long-term exposure to stress could cause physiological changes that negatively affect the animal health (Altan et al. 2003;Sohail et al. 2010;Simitzis et al. 2012). Therefore, elucidating the mechanism governing these stress-induced changes is important for improving the health and performance of farm animals (Daramola et al. 2013).
The effect of stress on animals' reproductive performance has been widely reported. Chrastil found that increased stocking density can significantly reduce the laying performance of hens and change the serum levels of superoxide dismutase (SOD), total antioxidant capacity (T-AOC) and malondialdehyde (MDA), which are related to antioxidant capacity and stress (Chrastil 1990). Onbasilar and Aksoy observed that in high-density cage conditions, the egg production rates of laying hens were significantly reduced, and this decrease was accompanied by an increased concentration of plasma corticosterone (Onbaşılar and Aksoy 2005). Rozenboim et al. reported that heat stress decreases egg production in White Leghorns; concurrently, they it could reduce reproductive performance by lowering plasma luteinising hormone (LH) and follicle-stimulating hormone (FSH) levels (Rozenboim et al. 2007). These studies indicated that stress could significantly affect animals' production performance, antioxidant capacity and reproductive health in poultry. Meanwhile, the underlying regulatory mechanism that integrates the stress signals and reproductive processes in animals needs further research.
It was reported that the hypothalamic-pituitaryadrenal (HPA) axis plays an essential role in response to stress. As under stressful situations, the HPA axis is activated, manifested by the secretion of corticotropin-releasing hormone (CRH) from the hypothalamus. Then, CRH can promote the secretion of adrenocorticotropic hormone (ACTH) from the pituitary gland and finally enhance the secretion of adrenocortical hormones from the adrenal glands (Star et al. 2008;Sohail et al. 2010). These stressed hormones can influence the release of gonadotropin-releasing hormone (GnRH) from the hypothalamus and inhibit the release of GnRH-induced LH from the pituitary gland. Consequently, the concentration of luteinising hormone receptor (LHR) and the sensitivity of gonads to LH were reduced, which eventually mediated changes in animal reproductive functions, including reproductive organ development, sexual maturation and egg development (Sapolsky et al. 2000;MacManes et al. 2017). As GnRH and LH (reproductive hormones) are related to the hypothalamic-pituitary-gonadal (HPG) axis. Therefore, these studies suggested that the hypothalamus and pituitary gland play essential roles in integrating stress signals and the reproductive process.
Poultry production in China is developing gradually towards intensive management and large-scale farming. Higher stocking density is an effective means of reducing the cost of livestock production and enhancing economic efficiency. In contrast, a high stocking density could cause stressful conditions in laying hens, and eventually changes in behaviour, egg production, and the fertility rate (Asghar Saki et al. 2012;Oliveira-Ferrer et al. 2014;Janczak and Riber 2015). In this study, high stocking density was used as a model for stress on laying ducks to investigate the potential roles of HPG in the relationship between stress signals and reproductive processes. At the start of the experiment, the effects of high stocking density on egg production in laying ducks was investigated, and later on, the transcriptome changes were examined in the hypothalamus, pituitary gland, ovary and follicular membranes. Most importantly, the findings of this study provide a theoretical basis for further work attempting to improve laying duck management, production performance, welfare and health conditions.

Experimental animals
The Jinding ducks used, in this study, were obtained from the Huaning Agricultural Company, Xichang City, Sichuan Province. The Jinding duck breed was selected due to the higher egg production rate and is already famous as a breed of the South China region (Yang 2011). The selected 360 ducks were divided into five groups with respect to their densities of 4-8 birds/m 2 , respectively, and each group was repeated twice (36 ducks per replicate). Hence, a total of ten pens were used with available space that is equal to the number of ducks divided by density per grid. All ducks had free access to water and feed. Drip-nipple water supply lines provided fresh water. The pelleted meals were provided in feed troughs on one side of each pen. All ducks were kept in the same environmental control farm. The ducks were fed on wire-floor pens and exposed to natural lighting and temperature throughout this study. The laying nest was set as one nest for 12 ducks. All laying ducks were female, and the sex was identified on the day of hatching. According to the experimental design, all selected ducks (a total of 360 ducks) were 21 weeks old at the start of the experiment because the Jinding ducks laying period started at the age of 140 d approximately. All experimental procedures were carried out according to the regulations of the animal ethics and welfare committee (AEWC) of Sichuan Agricultural University China (approval code: AEWC2016, 6 January 2016).

Measured traits
Ducks were individually weighed every week from 21 to 35 weeks of age. Eggs were collected and weighed daily. Feed consumption of each density treatment group was recorded on weekly basis from 21 to 35 weeks. The feed intake and egg mass ratio were measured during 21-35 weeks.

Sampling
At 30 weeks of age, laying ducks were slaughtered to collect samples. Three ducks were randomly selected from the low-density group (n ¼ 72, 4 birds/m 2 ) and the high-density group (n ¼ 72, 8 birds/m 2 ). These six ducks were slaughtered to enable tissue samples to be collected. The specific sample collection details, including the hypothalamus was collected at the bifurcation point of the septopallio-mesencephalic tract and ended after the cerebellum became apparent. This reference was used to determine the specific location of the hypothalamus and pituitary (Cheng Jiayue 1994). Similarly, we determined the location of the ovaries based on the book 'Colour Anatomy of Ducks' and subsequently obtained the ovarian matrix while collecting the largest follicles (F1) on the ovaries. The largest follicle is designated as the F1 follicle, which is the most mature and is destined for the next ovulation. This follicle was punctured with forceps, rinsed with PBS, and placed in a cryopreservation tube. In brief, the hypothalamus, pituitary, ovary and follicular membrane (diameter range: 8 ± 10 mm) were collected separately, rapidly frozen in liquid nitrogen and then stored at À80 C for further analysis.

Construction libraries and RNA-seq sequencing
Total RNA was extracted from tissues using TRIzol reagent (Beyotime, Shanghai, China) according to the manufacturer's instructions. The degree of RNA degradation was analysed by agarose gel electrophoresis (BIO-RAD, Hercules, CA), and the samples were separated by electrophoresis in 1.5% agarose gel. A NanoDrop 2000 spectrophotometer (NanoVue TM Plus Spectrophotometer, Holliston, MA) was used to detect the purity of RNA (OD 260/280 ! 1.8, OD 260/ 230 ! 1.0), the concentration of RNA was detected with a Qubit version 2.0 fluorometer (QuantiFluorV R ssDNA System, USA) and the integrity of RNA was measured with an Agilent 2100 Bioanalyzer (RIN value !8.0, 28S/18S !1.5) (Agilent, Waldbronn, Germany). Following the manufacturer's recommendations, sequencing libraries were generated using the NEBNext V R UltraTM RNA Library Prep Kit for Illumina V R (NEB, Ipswich, MA) and index codes were added to attribute sequences to each sample. Briefly, mRNA was purified from total RNA using oligo (dT) magnetic beads. Using divalent cations under elevated temperature, fragmentation was performed in NEBNext First Strand Synthesis Reaction Buffer (5X). The first-strand cDNA was synthesised using a random hexamer primer and M-MuLV Reverse Transcriptase (RNase H). Second strand cDNA synthesis was performed using DNA Polymerase I and RNase H. The remaining overhangs were blunt and ended by exonuclease/ polymerase activities. After adenylation of the 3 0 ends of DNA fragments, NEBNext Adaptor with a hairpin loop structure was ligated to prepare the DNA for hybridisation. cDNA fragments measuring 250-300 bp in length were selected by purifying the library fragments with the AMPure XP system (Beckman Coulter, Beverly, MA). Then, 3 lL of the USER enzyme (NEB, Ipswich, MA) was used with size-selected, adaptorligated cDNA at 37 C for 15 min, followed by 5 min at 95 C before PCR analysis. PCR was then performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers and the Index (X) Primer. Finally, PCR products were purified (AMPure XP system), and the library quality was assessed on the Agilent Bioanalyzer 2100 system. The library preparations were sequenced on an Illumina Novaseq platform.

RNA-seq data analysis
High-throughput sequencing data were originally presented as raw image files (in the format BCL). After base recognition by bcl2fastq Conversion Software version 1.8.4 (Illumina, USA), the data were converted into raw data. In the quality control step, clean reads were obtained by removing reads containing adapter, poly-N and low-quality read from raw data. All downstream analyses were based on high-quality clean data. The clean data were matched to the reference genome [duckbase.RefSeq.v4] using Hisat2 version 2.1.0 (Amazon Web Services, ). The transcripts were assembled using Stringtie version 1.3.3b, and the genes were quantified. The differentially expressed genes (DEGs) were calculated by the edgeR package in R version 3.6.1 software (Ross Ihaka, New Zealand). The screening criteria were jLogFCj !1, p value < .05. Gene expression levels were quantified using CPM.

Gene ontology and KEGG pathway analysis
Gene ontology (GO) enrichment analysis of DEGs was performed using the David website (https://david. ncifcrf.gov/home.jsp). This analysis includes three parts: molecular function, biological process and cellular component. A specific p value was employed to determine the GO enrichment of DEGs. Pathway enrichment analysis was assessed using KOBAS version 3.0 (Bejing, China) (http://kobas.cbi.pku.edu.cn/). The Cluster Profiler R package was used to test the statistical enrichment of DEGs in KEGG pathways.

Protein-protein interaction analysis
The protein-protein interaction networks were performed using the interaction relationship in the STRING database (http://string-db.org/). Visual editing of differential protein-protein interaction network data files was performed by Cytoscape software version 3.7.2 (http://www.cytoscape.org/).

Statistical analysis
Data were analysed by a repeated-measures ANOVA procedure of SPSS version 20.0 software (SPSS Inc., Chicago, IL). When the density treatment was significant (p < .05), the means were compared using the LSD multiple comparison procedure with SPSS software.

Production performance
As shown in Figure 1, the duck feed intake decreased with increasing stocking density. During the whole experimental period (21-35 weeks of age), compared with the low-density group (4 birds/m 2 ), the feed intake of ducks in the high-density group (8/m 2 ) decreased by 7.40 À 23.44%, with an average decrease of 15.15%. The egg-laying rate reached a peak at 27 weeks of age in the low-density group (4 birds/m 2 ) and gradually decreased with increased age. In each experimental group, increased stockingrearing density was accompanied by a decrease in the egg production rate of duck flocks. Meanwhile, during 21-35 weeks of age, the egg production of the high-density group (8 birds/m 2 ) decreased by 13.04 À 63.55% compared with the low-density group (4 birds/m 2 ), with an average decrease of 35.51%. By calculating the feed/egg ratio, we found that the ducks at 27 weeks had the lowest feed/egg ratio among all groups. After 32 weeks of age, the feed/egg ratio increase further supported the downward trend in egg production. The feed/egg ratio increased by 60.55% in the high-density group (8 birds/m 2 ) compared with the low-density group (4 birds/m 2 ). Through repeated measurement ANOVA, the results are shown in Table 1, a p value of time <.001 indicates that the difference in data at each time point is significant, and a p value of time Ã group >.05 suggests no interaction between time and group. These results explain that the effect of the time factor (i.e. 21-35 weeks of age) does not vary with the grouping (i.e. 5 density groups). The above-mentioned results indicate that high-density treatment may decrease ducks' feed intake and egg production rate, as well as, poor feed to egg ratio, was found in highdensity groups. Consequently, this approach provided a suitable model for investigating the role of the hypothalamic-pituitary in integrating stress signals and the reproductive process.

Overview of the RNA-seq data
Approximately 380 million clean reads were obtained from the hypothalamus, pituitary, ovary and 8-10 mm follicular tissue of all sequenced individuals with three biological replicates. More than 80% of clean reads were aligned to the reference genome for each sample. The detailed data for each sample are presented in Table S1. The Pearson correlation coefficient indicated that the choice of experimental samples was consistent and reliable ( Figure S1).

Analysis of differentially expressed genes
The DEGs were identified between the low-and high-density treatment groups (low-vs. high-density group) in four tissues under jLog 2 FCj ! 1 and p value < .05. A total of 469, 509, 4284 and 210 DEGs were obtained from the hypothalamus, pituitary, ovarian and follicular membranes, respectively (Table 2 and Figure 2).
Through David (https://david.ncifcrf.gov/home.jsp) (Huang et al. 2009), 14 significant GO terms were enriched in the hypothalamus (Figure 3(A)), and the top 3 were an integral component of membrane (GO:0016021), G-protein coupled receptor activity (GO:0004930) and acute-phase response (GO:0006953). A total of 84 significant GO terms were enriched in the ovary (Figure 3(B)), and the top 3 were integral components of the membrane (GO:0016021), calcium ion binding (GO:0005509) and cell communication (GO:0007154). Furthermore, 22 significant GO terms were enriched in the pituitary (Figure 3(C)), and the top 3 were extracellular space (GO:0005615), an integral component of membrane (GO:0016021) and immune response (GO:0006955). A total of 6 terms were enriched with significant p values (p <.05; Figure  3(D)). The top three terms were integral components of the membrane (GO:0016021), positive regulation of cell differentiation (GO:0045597) and cholesterol homeostasis (GO:0042632).
By using KOBAS (http://kobas.cbi.pku.edu.cn/) (Xie et al. 2011), KEGG enrichment analysis was conducted based on all screened DEGs. In the hypothalamus, the most significant pathway (the screening standard is p values < .05) was neuroactive ligand-receptor interaction; thus, in the pituitary, the DEGs were primarily enriched in cytokine cells, factor receptor interactions, RIG-I-like receptor signalling pathways and cell adhesion molecules (CAMs). In the ovary, the most significant pathway was steroid hormone biosynthesis, which may be related to the synthesis of reproductive hormones. The most prominent pathway in the follicular membrane was the peroxisome proliferatoractivated receptors (PPARs) signalling pathway involved in lipid metabolism (Figure 4).

Expression of genes related to the stress process is altered in the hypothalamus and pituitary tissues
We screened 8 stress-related pathways from the hypothalamus and pituitary, including 31 DEGs. In the adipocytokine signalling pathway, we observed three downregulated genes: POMC (logFC À2.23, p value 1.12E-03), AGRP (logFC À2.88, p value 4.39E-04) and ACSL5 (logFC À1.29, p value 3.29E-02). A total of 11 pathways related to the stress process were screened in the pituitary, including 50 DEGs. In protein processing in the endoplasmic reticulum, we detected 1 downregulated gene, EIF2AK2 (logFC À2.23, p value 1.12E-03). The stress-related pathways involved in the hypothalamus and pituitary are listed in Table 3, and their DEGs are listed in Table S2.

Expression of genes related to reproductive processes is altered in ovarian tissues
We screened 9 pathways related to the reproduction process from the ovarian tissues, including 121 DEGs. In steroid hormone biosynthesis, we detected 14 DEGs, 8 of which were upregulated and 6 downregulated. Among these genes, HSD17B2 (logFC À3.53, p value 3.91E-06) was related to oestrogen activity (Cornel et al. 2019). A total of 8 pathways related to the reproductive process were screened from the follicular membrane tissues, and these pathways covered 13 DEGs. In focal adhesion, we observed 3 downregulated genes, most important was PGF (logFC À1.72, p value 4.83E-02), which has been reported to be involved in ovulation (Takahashi et al. 2018). The reproduction-related pathways involved in the ovary and follicular membrane are listed in Table 4, and the enriched DEGs are listed in Table S2.

Genes related to stress processes were clustered with genes involved in the reproduction process
The hypothalamus and pituitary regulation of other tissues and organs is primarily achieved through neuroendocrine mechanisms. Therefore, we investigated the secreted proteins in the hypothalamus and pituitary tissue of laying ducks. Then, we compiled all the DEGs of the KEGG pathway (Table 3) and compared these DEGs with the Human Protein Atlas database (http://www.proteinatlas.org/). A total of nine genes encoding secreted proteins from the hypothalamus and pituitary gland were screened out (Table 5, marked with red). Then, we analysed the relations between these genes and the genes expressed in ovary tissues and follicle membranes (Table 5). We performed cluster analysis on the expression levels of all these genes ( Figure 5(A)). The results showed that POMC, DRD3 and AVPR1B were grouped into one cluster, and NPVF and OPRM1 were grouped into another cluster ( Figure 5(A)). In addition, the protein-protein interaction analysis was performed by STRING (http:// string-db.org/), and the results were optimised using Cytoscape ( Figure 5(B,C)). In Figure 5(B), the results calculated by the cytoHubba plugin showed that the degree of POMC was the most, indicating that POMC might play an essential role in the stress response caused by stocking density treatment. Similarly, in Figure 5(C), GNRH1 and POMC are highlighted, suggesting that they play key roles among these DEGs.   Figure 2. Volcano plot of differentially expressed genes. Horizontal coordinates represent the fold change in genes between the low-density group (4 birds/m 2 ) and the high-density group (8 birds/m 2 ). The ordinate represents the statistical significance of changes in gene expression. The smaller the p value is, the larger À log 10 p is. Each dot in the image represents a gene; dots marked in grey indicate no significant differences in genes, and dots marked in blue and red indicate significant differences.

Discussion
It is generally believed that a high stocking density directly affects animals' production performance and health status (Feddes et al. 2002;Xie et al. 2014;Tsiouris et al. 2015). In this study, the effects of five different stocking densities (4, 5, 6, 7 and 8 birds/m 2 ) on the laying performances of ducks were compared under the net rearing mode. Compared with the low-density group (4 birds/m 2 ), the laying rate and feed intake were decreased by 35.51% and 15.15% on average, respectively, in the high-density group (8 birds/m 2 ), indicating that the increase in density could negatively affect the production status of laying ducks. Meanwhile, the results of this study demonstrated that the ducks in the 8 birds/m 2 group had the worst feed intake, egg production rate and feed/egg ratio. These results were consistent with the findings of previous studies, which also indicated that the production performance of laying ducks was significantly reduced under high-density rearing ( Skrbi c et al. 2011).   GCG# FOS# EGF# GNG13# GRP# DRD3" GRM1# GRM5" GRM7" OPRM1" MRAP2" Ovary ADRA2A" CRHR2" ESR2# CCKAR# GABRR1" AVPR1B# STAR# THBS1" VWF" FOS" LUM" CCL24" EOMES# B2M" FBLN1" ELFN1" CALR3# LHFPL2" IL13RA1" EGR1" Follicular membrane MMP27"MGLL# PGF# SOST" FGF19" Arrows indicate low-density group vs. high-density group differential gene up-and down-direction.   Table 5.
We performed a subsequent transcriptome analysis to determine whether the decline in production performance is related to stress. In the hypothalamus and pituitary, 78 DEGs, including NPVF, TBX19, POMC, MCHR1 and EIF2AK2, were screened from stress-related pathways. Among these genes, the mRNA expression level of the POMC gene was upregulated in the high-density group (8 birds/m 2 ), and it has been reported to play key functions in increasing the satiety of the stomach and suppressing appetite, which was consistent with our findings (Chu et al. 2014;Liu et al. 2018). TBX19, another DEG in the stress-related pathway, was proven to activate POMC expression . According to the UniProt database, ACTH is one of the proteins encoded by the POMC gene. Moreover, the MCHR1 gene, which affects feeding and anxiety, was also upregulated in the high-density group (Roy et al. 2006;Luthin 2007). Meanwhile, EIF2AK2 exhibited a similar expression mode as the above-mentioned DEGs, and it has been reported to play an important role in the integrated stress response (ISR) (Pakos-Zebrucka et al. 2016). With the differential expression of the genes related to the stress pathway, it could be confirmed that the laying ducks had a stress response under high-density feeding conditions.
Reproduction, a complex and sophisticated physiological process, requires coordination of the HPG axis, which promotes the growth and maturation of germ cells (Ye et al. 2019) and the synthesis of gonadal steroids in both males and females gender (Oyola and Handa 2017). In the HPG axis, GnRH is widely recognised as an important hormone controlling reproduction (Stohn et al. 2016). We determined that the expression of the GnRH1 gene in the hypothalamus was downregulated in the high-density group. Studies have also shown that dopamine in the hypothalamus plays an important role in poultry reproduction (Kagya-Agye et al. 2012). However, we did not observe any differential expression levels of genes related to dopamine in the hypothalamus. Only the pituitary, the dopamine type 3 receptor (DRD3), was observed. Hernandez et al. showed that a lack of active DIO3 in mice could lead to excessive T3 in the brain, whether developmental or adulthood (Hernandez et al. 2006(Hernandez et al. , 2010. This phenomenon is associated with decreased anxiety, depression-like behaviour, and increased hyperactivity (Stohn et al. 2016). In this study, the downregulation of DIO3 gene expression in the lowdensity group was consistent with the findings of previous studies, thereby reducing the stress behaviour of laying ducks, and the feed intake and egg production rate were higher than those in the high-density group. Therefore, the downregulation of such genes as GnRH1 in the 8 birds/m 2 group indicated that highdensity feeding could decrease the egg production rate by affecting the function of HPG. Meanwhile, in the ovaries and follicular membranes, we screened a total of 132 DEGs enriched in reproductive pathways, including the GnRH signalling pathway (14 DEGs), the steroid hormone biosynthesis pathway (7 DEGs), the progesterone-mediated oocyte maturation pathway (22 DEGs) and the oocyte meiosis pathway (19 DEGs). The HSD17B2 gene is related to oestrogen activity (Cornel et al. 2019). Furthermore, the EGR1 and MMP27 genes, reportedly associated with follicular development, were downregulated in the high-density group (Goldman and Shalev 2004;Hatzirodos et al. 2014). The FOS gene has been widely studied concerning steroid hormone synthesis in the ovaries (Oliveira-Ferrer et al. 2014). The expression changes of these genes indicated that the follicular development and ovarian functions of laying ducks were affected under high-density feeding conditions. These results helped explain the decrease in egg production observed in the high-density population.

Conclusions
We demonstrated that increasing the stocking density from 4 to 8 birds/m 2 during the laying period decreased the egg production rate of laying ducks by 13.04 À 63.55% and decreased the feed intake by 7.40 À 23.44%. The transcriptome analysis indicated that stocking density affects the stress response of laying ducks through the hypothalamus-pituitary, and the stress signals are subsequently transmitted to affect gene expression related to the reproduction process in the ovarian and follicular tissues of laying ducks. In addition, the hypothalamic expression of POMC and GnRH1 may play a central role in integrating stress signals and the reproductive processes of laying ducks.

Author contributions
All experiments were designed by Fajun Pu, Hehe Liu, and Xia Xiong. Yang Xi, Shengchao Ma, and Yanying Li performed the genome data analysis. Chaowu Yang, Rongping Zhang, Lili Bai, Liang Li, and Fajun Pu participated in the sample collection and preparation. Hehe Liu, Xia Xiong, Fajun Pu, and Jianmei Wang completed the manuscript writing. All authors have read and approved the manuscript.

Disclosure statement
The authors declare that they have no competing interests.