Skip to main content
Download PDF

Immune profiles of male giant panda (Ailuropoda melanoleuca) during the breeding season

BMC Genomics volume 22, Article number: 143 (2021) Cite this article

Abstract

Background

The giant panda (Ailuropoda melanoleuca) is a threatened endemic Chinese species and a flagship species of national and global conservation concern. Life history theory proposes that reproduction and immunity can be mutually constraining and interrelated. Knowledge of immunity changes of male giant pandas during the breeding season is limited.

Results

Here, we researched peripheral blood gene expression profiles associated with immunity. Thirteen captive giant pandas, ranging from 9 to 11 years old, were divided into two groups based on their reproductive status. We identified 318 up-regulated DEGs and 43 down-regulated DEGs, which were enriched in 87 GO terms and 6 KEGG pathways. Additionally, we obtained 45 immune-related genes with altered expression, mostly up-regulated, and identified four hub genes HSPA4, SUGT1, SOD1, and IL1B in PPI analysis. These 45 genes were related to pattern recognition receptors, autophagy, peroxisome, proteasome, natural killer cell, antigen processing and presentation. SUGT1 and IL1B were related to pattern recognition receptors. HSP90AA1 was the most up-regulated gene and is a member of heat shock protein 90 family. HSP90 contributes to the translocation of extracellular antigen. KLRD1 encodes CD94, whose complex is an inhibitor of the cytotoxic activity of NK cells, was down-regulated. IGIP, which has the capability of inducing IgA production by B cells, was down-regulated, suggesting low concentration of IgA in male giant pandas. Our results suggest that most immune-related genes were up-regulated and more related to innate immune than adaptive immune.

Conclusions

Our results indicated that breeding male giant pandas presented an immunoenhancement in innate immunity, enhanced antigen presentation and processing in cellular immunity compared to non-breeding males. The humoral immunity of male giant pandas may show a tendency to decrease during the breeding season. This study will provide a foundation for further studies of immunity and reproduction in male giant pandas.

Background

Reproductive activity, with a high metabolic cost, is associated with body and immunological conditions [1]. Life history theory proposes that reproduction and immunity can be mutually constraining and interrelated due to the optimal allocation of limited nutrients and energy [2]. An immune response is responsible for a substantial energetic costs [3]. Energy investment in reproduction leads to a corresponding decrease in immune investment, resulting in a trade-off between the two [4]. It is challenging to elucidate the underlying reproduction and immunity trade-off mechanisms, while it is easy to observe and record immune traits during reproduction [5]. Many male-focused studies on a variety of species have documented that males have reduced innate immunity and lower cellular immunity during reproductive periods [3, 6, 7]. However, several studies have found that during reproduction there is immunoenhancement of cellular immunity and higher resistance against bacteria [8, 9].

The giant panda (Ailuropoda melanoleuca), known as China’s national treasure, is a Chinese flagship species, a threatened species of global concern and a worldwide symbol of conservation [10]. Our research group previously showed that maternal giant panda immunity undergoes dramatic changes during estrus, early pregnancy, and late pregnancy [11]. Males also undergo reproductive changes, although not as dramatic as females. Male giant pandas reach sexual maturity at approximately 8 years old and undergo increases in testes volume, androgen concentrations and sperm production each breeding season thereafter [12]. According to previous studies, male giant pandas’ reproductive cycle was divided into three periods: breeding season (February–May), prebreeding season (October–January) and non-breeding season (June–September) [13]. However, most studies focused on the reproductive behaviors [13,14,15] and little is known about the immune change of male giant panda during the breeding season. The greater understanding of giant panda maternal immunity has allowed better health care and disease prevention and therefore there is a real need for male reproductive immune system research.

Peripheral blood is the important vehicle of the immune system. Many methods can quantify and measure immune traits from blood samples, such as total leukocyte counts [16], cytokines [17], complement and lysozyme activity [18] and phagocytosis activity [19]. Alternatively, immunocompetence can be quantified by phytohaemagglutinin challenge test [8] and sheep red blood cells injection test [20]. Measuring only one or few immune parameters, which treats the immune system like a “black box” [4], may be insufficient to reflect the overall immune system status and change. Transcriptome analysis is a robust tool that investigates the immune system by assessing transcript abundance changes in blood on a genome-wide scale, earning its place in immune function study [21]. Therefore, we aimed to compare the immune profiles of breeding and non-breeding male giant pandas. Transcriptome analysis was used to quantitatively evaluate transcript levels and identify the immune-related differentially expressed genes (DEGs) and pathways. This study will provide a foundation for further reproductive immunity and disease prevention studies on breeding male giant pandas.

Results

Reads sequencing and processing

Raw Illumina RNA-seq data were converted into clean reads data. All raw data has been deposited at NCBI Sequence Read Archive under the project accession no. PRJNA631846. A total of 95.49 Gb of paired-end clean data were generated. FastQC showed that the percent of Q30 was above 85%.

HISAT2 mapping results revealed all sample overall alignment rates were between 86 and 90%. Read summarization counted by program featureCounts was converted into a numerical matrix. PCA results based on normalized matrix demonstrated that the 13 samples were divided into two groups from different dimensions (Fig. 1). Giant pandas in the breeding season were clustered into one group, while the non-breeding individuals were clustered into another group.

Fig. 1
figure1

PCA analysis of 13 samples

Full size image

Identification of DEGs

We detected 1128 genes with changed expression level, given the FDR threshold (Additional file 1: Table S1). By setting a cutoff of log2FC, 318 up-regulated DEGs and 43 down-regulated DEGs were identified in the breeding season compared to the non-breeding season. Two hundred seventy-five in 318 upregulated genes had annotations, while 33 in 43 down-regulated genes had annotations.

In the top 10 up-regulated DEGs, eight genes were involved in genetic information processing, mainly in transcription. HSP90AA1 is a member of heat shock protein 90 family. HSP90AA1 participates in numerous immune processes, such as antigen processing and presentation, Th17 cell differentiation, and NOD-like receptor signaling pathway. PSMD7 encodes proteasome 26S subunit. Proteasome plays a great role in innate and adaptive immune responses.

In top 10 down-regulated DEGs, five genes were related to genetic information processing, such as transcription, translation, and protein export. IGIP (Immunoglobulin A inducing protein) belongs to the Immunoglobulin A regulatory factors family. KLRD1 (killer cell lectin-like receptor subfamily D member 1) is associated with natural killer cell immunity.

Gene ontology enrichment of DEGs

Up-regulated DEGs were enriched in 69 GO terms, being 22 terms in biological process, 39 terms in cellular component and 8 terms in molecular function (Fig. 2). Down-regulated DEGs were enriched in 18 GO terms, which were 8 terms in cellular component and 10 terms in molecular function (Fig. 3). All GO term enrichments are shown in Additional file 2: Table S2. There were some overlap top-level cellular component terms between up-regulated DEGs and down-regulated DEGs, such as protein-containing complex (GO:0032991), cell (GO:0005623), cell part (GO:0044464) and organelle (GO:0043226). For down-regulated DEGs, the most significantly enriched molecular function term was cytochrome-c oxidase activity (GO:0004129). For up-regulated DEGs, the enriched GO terms in molecular function included gene expression (GO:0010467) and RAGE receptor binding (GO:0050786) which was associate with immune and inflammatory responses.

Fig. 2
figure2

Partial GO enrichment of up-regulated DEGs

Full size image
Fig. 3
figure3

GO enrichment of down-regulated DEGs

Full size image

KEGG pathway enrichment of DEGs

Using an overrepresented analysis, we performed KEGG enrichment analysis for further understanding of DEGs. Up-regulated DEGs and down-regulated DEGs were enriched in four and two KEGG pathways respectively (Fig. 4). Up-regulated genes were enriched in ribosome (aml03010), spliceosome (aml03040), oxidative phosphorylation (aml00190) and thermogenesis (aml04714) pathways. Ribosome and spliceosome pathways were associated with genetic information processing. Thermogenesis was the child term of environmental adaptation pathway. Oxidative phosphorylation was the downstream term of thermogenesis. When focusing on down-regulated genes, we found the protein export (aml03060) and ribosome (aml03010) pathway were significantly enriched. Protein export was the child term of genetic information processing pathway.

Fig. 4
figure4

KEGG enrichment of up-regulated and down-regulated DEGs

Full size image

We identified the biological impact of the breeding stage and the direction of the impact using a Dynamic Impact Approach (DIA). The summary of KEGG main categories and sub-categories is shown in Fig. 5. Among the main categories of KEGG, the category “Genetic Information Processing” was the most impacted, followed by “Organismal Systems” and “Cellular Processes”. Except for inhibition of “Membrane Transport” and “Digestive System”, the flux values of sub-categories were activated. The sub-category “Transcription” was the most impacted, followed by “Sensory System”. The top 20 most-impacted pathways are shown in Fig. 6. The most impacted pathway was “Fatty acid elongation in mitochondria” followed by “Progesterone-mediated oocyte maturation”. “Notch signaling pathway” was the only inhibited pathway. Among the top 20 pathways, “NOD-like receptor signaling pathway” and “Antigen processing and presentation” were associated with the immune system.

Fig. 5
figure5

Summary of the main categories and subcategories of KEGG pathways analyzed by DIA. On the right are the bar denoting the overall impact (in blue) and the shade denoting the effect on the pathway (from green (inhibited)—to red (activated)). Darker the color larger the activation (if red) or inhibition (if green) of the pathway. “B” represents the “breeding season”. “non-B” represents the “non-breeding season”

Full size image
Fig. 6
figure6

The 20 most impacted pathways analyzed by DIA. On the right are the bar denoting the overall impact (in blue) and the shade denoting the effect on the pathway (from green (inhibited)—to red (activated)). Darker the color larger the activation (if red) or inhibition (if green) of the pathway. “B” represents the “breeding season”. “non-B” represents the “non-breeding season”

Full size image

Expression of immune-associated genes

We obtained 45 immune-related genes and clustered them into 12 key categories according to KEGG annotation (Fig. 7). We also plotted the heatmap of immune-related genes to visualize their expression in all samples (Fig. 8). These categories were roughly divided into innate immune entries and adaptive immune entries. Innate immune system entries consisted of C-type lectin receptor, NOD-like receptor, autophagy, peroxisome, proteasome, natural killer cell, cytokine and chemokine, and TNF signaling pathway. Adaptive immune entries consisted of antigen processing and presentation, T cell receptor signaling pathway, Th17 cell differentiation, and IL-17 signaling pathway.

Fig. 7
figure7

Chord diagram of categories of immune-associated genes. Genes were divided into twelve branches according to KEGG pathways: A –Antigen processing and presentation; B –Autophagy; C –C-type lectin receptor signaling pathway; D –Cytokine and Chemokine; E –IL-17 signaling pathway; F –Natural killer cell mediated cytotoxicity; G –NOD-like receptor signaling pathway; H –Peroxisome; I –Proteasome; J –T cell receptor signaling pathway; K –Th17 cell differentiation; L –TNF signaling pathway

Full size image
Fig. 8
figure8

Heat map plot of 45 immune-related genes. The expression values of 13 pandas are presented after being centered and scaled in the row direction. Each column represents a specimen and each row represents a gene. Red color indicates genes which were up-regulated and green color indicates genes which were down-regulated

Full size image

The expression trends of 45 genes were consistent, mostly up-regulated, while KLRD1 (killer cell lectin-like receptor subfamily D member 1), IL15 (interleukin 15) and TRAF1 (TNF receptor-associated factor 1) were down-regulated. CLEC4E (C-type lectin domain family 4 member E) is a member of the C-type lectin receptor signaling pathway. NAMPT (nicotinamide phosphoribosyltransferase) and GABARAPL1 (GABA type A receptor associated protein like 1) participate in the NOD-like receptor signaling pathway. BECN1 (Beclin1), PRDX5 (peroxiredoxin 5) and PSME1 (PA28 alpha) shows the great function in autophagy, peroxisome and proteasome respectively. VAV1 (guanine nucleotide exchange factor) and PLCG2 (phosphatidylinositol phospholipase C gamma-2) are linked to natural killer cell. IL15 together with IL1R2 (interleukin 1 receptor type 2) are two important cytokines. Lastly, CD3D (T-cell surface glycoprotein CD3 delta chain) and CD3G (T-cell surface glycoprotein CD3 gamma chain) are associated with T cell receptor signaling pathway and Th17 cell differentiation.

Protein-protein interaction network of immune-associated genes

All immune-associated genes were converted into proteins by STRING. A total of 64 interaction edges between 36 nodes were extracted from the database after removing 9 isolated nodes. What’s more, we calculated the hub genes by using cytoHubba. The score of 36 genes were calculated by topological analysis methods. (Additional file 3: Table S3). We plotted the network diagram to illustrate interaction among proteins (Fig. 9). HSPA4 (heat shock 70 kDa protein 4), SUGT1 (SGT1 homolog), SOD1 (superoxide dismutase 1), and IL1B (interleukin 1 beta) were at important position within the interaction network.

Fig. 9
figure9

Protein-protein interaction network of DEGs. The nodes represent proteins and edges represent pair-wise interactions. The size of the nodes is score of the protein calculated by cytoHubba. The red nodes represent up-regulated proteins. The blue nodes represent down-regulated proteins

Full size image

Real-time quantitative PCR (qRT-PCR) validation

Twelve DEGs (MFAP1, HSP90AA1, PSMD7, S100A9, SOD1, CD3D, RPL9, KLRD1, IGIP, SEC61B, PHF5A, VMA21) were selected for verification. As shown in Fig. 10, the results of qRT-PCR indicated similar expression tendencies with transcriptome sequencing. The qRT-PCR validation further improves the reliability of the present study.

Fig. 10
figure10

The expression level of genes verified by qRT-PCR. “*” represents the P-value< 0.05, “**” represents the P-value< 0.01, “***” represents the P-value< 0.001. Data were shown as mean ± SD. The left axis represents gene expression levels verified by qRT-PCR. The right axis represents the expression levels in TPM units of RNA-seq

Full size image

Discussion

Animals in nature need to balance resource allocation between reproduction and self-maintenance, and immunity is a major component of self-maintenance [22]. The reproduction and conservation of giant pandas has been, and continues to be, of global concerns [12, 23]. Our research group’s previous work investigated immune changes at four key phases of female giant panda reproduction [11]. However, the immune performance of male giant pandas during reproduction has been little studied. Here we investigated the immune changes in 8 male giant pandas over the breeding season compared with 5 males in the non-breeding season. We monitored the expression of immune-related genes based on peripheral blood transcriptome. We identified 45 immune-related genes with altered expression, mostly up-regulated, in breeding males compared to non-breeding males.

The GO term enrichment of “translation”, “peptide biosynthetic process” and “structural constituent of ribosome” and KEGG pathway enrichment of “ribosome” were observed in up-regulated genes. DIA revealed that “Genetic Information Processing” was the most impacted pathway and was overall, strongly activated. These results suggest an increased requirement for protein synthesis in breeding male giant pandas. The amplification of protein synthesis was also reported in male freshwater spotted snakehead (Channa punctatus) during reproductive phases [24]. The enrichment of the ribosome pathway is consistent with findings in sheep testes, and indicates that the normal function of the ribosome plays an essential role in spermatogenesis [25]. The dramatically up-regulated genes were enriched in spliceosome, which removes noncoding introns from transcribed mRNA precursors, suggesting spliceosome is very important in producing necessary gene products related to male sexual development [26]. Oxidative phosphorylation was another enriched pathway in our study. This pathway is an important ATP-related metabolic pathway and provides energy for male reproduction [26]. For the ‘Sensory System’ subcategory, we uncovered a clear upregulation of “Olfactory transduction” and “Phototransduction”. The male giant pandas may use olfactory and visual cues to assess their sexual partner during the breeding season [27]. The most impacted pathway from DIA analysis was “Fatty acid elongation in mitochondria”. Fatty acid elongation is associated with fatty acid synthesis [28]. The sperm polyunsaturated fatty acid content increases during the breeding season and sperm characteristics are affected by fatty acid composition [29, 30]. The need of fatty acid for spermatogenesis might in part explain the upregulation of fatty acid elongation in breeding male giant pandas. During spermatogenesis, Sertoli cells and germ cells can respond to Notch signaling that is crucial for germ cell differentiation and germ stem cell pool maintenance and migration [31, 32]. During mouse folliculogenesis, Notch signaling was reported to plays an important role in follicular development and angiogenic growth [33]. Moreover, two hub genes HSPA4 and SOD1 were 3.36 and 3.25-folder higher in breeding males than non-breeding males, respectively. The expression of HSPA4 is higher in germ cells of prenatal gonads [34] and SOD1 activity is higher in stallion during the breeding season [35]. This suggests HSPA4 and SOD1 are involved in spermatogenesis and antioxidant protection of sperm in male giant pandas. The up-regulated genes and enriched pathways may indicate that male giant panda reproductive systems prepare for breeding by triggering protein synthesis, energy generating and spermatogenesis.

Innate immune changes

The innate immune subsystem typically includes pattern recognition receptors, autophagy, antimicrobial peptides, and many cell types (e.g. dendritic cells, macrophages and natural killer cells), establishing the first line of defense against a wide range of invading pathogens [36, 37]. Moreover, the innate immune subsystem is responsible for the activation of the adaptive immune subsystem [36]. During the breeding season, captive male tree lizards reduced innate immunity [3], while the innate immunity showed no change in male Eurasian tree sparrows and temperate bats [6, 38]. Yet, Arabian and Thoroughbred horses presented an increased innate immunity [17]. Here, we explored the alteration of innate immunity from several aspects and found an enhanced innate immunity in male giant pandas.

We found several key genes referred to as pattern recognition receptors (PRRs) were upregulated, including CLEC4E (also known as Mincle), SUGT1 (SGT1 homolog), HSP90AA1, IL1B, and GABARAPL1 (LC3 paralog). Pattern recognition receptors mainly include Toll-like receptors (TLRs), C-type lectin receptors (CLRs), and NOD-like receptors (NLRs) [39]. CLRs were found to recognize microorganisms such as viruses, bacteria, and fungi, and then regulate the production of proinflammatory cytokines [39]. CLEC4E encodes macrophage-inducible C-type lectin (Mincle) who is a member of the CLRs family [40]. Mincle has been known to recognize dead cells and bacteria [40]. Evidence showed that Mincle was strongly up-regulated after skin injury and irritation, and mediated a severe inflammatory response and the production of IL1B [41]. The up-regulation of Mincle and IL1B in our study may suggest activation of the host immune responses and protection of giant pandas from skin infectious diseases. DIA revealed that NOD-like receptor signaling pathway was activated in male giant pandas during reproduction. NLRs exert function in inflammatory responses and tissue homeostasis [42, 43]. NOD1 as a member of the NLRs family recognizes invasive microbial pathogens which threaten homeostasis by specific peptidoglycans [39]. The SGT1 was reported to positively regulate NOD1 activation and depletion of SGT1 blocks multiple cellular responses caused by NOD1 activation [44]. Besides, HSP90, which is an evolutionarily conserved molecular chaperone, protects NOD1 from degradation and functions as a stabilizer [44, 45]. NLRs can interact indirectly with LC3 through a signaling cascade to regulate autophagy [45, 46]. Autophagy-associated genes (BECN1, GABARAPL1) were also up-regulated in our analysis. What’s more, we also found the up-regulation of NAMPT (also known as visfatin), which is a cytokine hormone [47]. The NAMPT was up-regulated fourfold in adult chicken testis compared with prepubertal chickens, suggesting the critical functions in spermatogenesis and steroidogenesis [47]. In the mouse testis, the expression of visfatin was found to be significantly associated with antioxidant enzymes activities when using dexamethasone treatment [48]. NAMPT inhibition resulted in the decreased production of many proinflammatory mediators in mouse macrophages [49]. Therefore, the up-regulation of SGT1 and NAMPT in breeding male giant pandas may regulate tissue homeostasis, antioxidant enzymes activities, and proinflammatory mediators.

Autophagy is a fundamental intracellular bulk degradation process with multiple roles in innate immune responses and cellular stress [50, 51]. Beclin1 and LC3 encoded by BECN1, GABARAPL1, respectively, were both up-regulated. Mammalian core autophagy-related proteins mainly involve several functional units, including the PI3K complex that is composed of Beclin1, the LC3 conjugation system and so on [52]. LC3 conjugation system regulates the elongation of the phagophore and promotes the completion of autophagosome formation [46]. During the breeding stage of testicular recovery, the expression of BECN1 and LC3 began to increase in South American plains vizcacha [53]. Tabecka-Lonczynska et al. observed an increase in beclin1 and LC3 synthesis and confirmed the function of autophagy in adult reproductive male European bison [51]. Up-regulated expressions of BECN1 and LC3 suggests the increased demand for maintaining homeostasis in male giant pandas during the period of reproductive activity.

Peroxisomes are crucial metabolic organelles which play central roles in lipid metabolism and ROS turnover [54]. Accumulating evidence suggests a new function for peroxisomes in microbial infection resolution and antiviral response [54, 55]. What’s more, peroxisomes have been pointed to play an important role in cell type-specific metabolic function in the testis and spermiogenesis [56]. Two antioxidant genes SOD1 and PRDX5 were up-regulated in the present study. ROS which has recently emerged as a signal factor in innate immune responses, is influenced by the disruption of redox balance in enzymes and subcellular compartments [55, 57]. Kwang et al. demonstrated that SOD1 tightly regulated the generation of ROS during virus infection [57]. Excess ROS production can induce lipid peroxidation, disrupting membrane characteristics of sperm [58]. Elevated expression of the PRDX5 and SOD1 improves the quality of porcine oocytes by modulating the ROS level [59]. Bernard et al. reported that human PRDX5 interacted with, or bound to, PRRs to activate a proinflammatory response [60]. PRDX5 can trigger the expression and release of IL1B [60]. The expression of PRDX5 and IL1B were both up-regulated in giant pandas, confirming an association between PRDX5 and IL1B. Collectively, peroxisomes are essential for the activation of the innate immune system and the normal function of testis in giant pandas.

The proteasome is responsible for poly-ubiquitinated substrates recognition and intracellular proteins degradation [61]. The proteasome system and autophagy are closely interconnected [62]. The proteasome is a multi-subunit protein complex, consisting of a 20S core particle and 19S regulatory particles [62]. Standard 20S proteasomes can be replaced by immunoproteasomes which are activated by PA28 complex in conditions of infection, inflammation, and an intensified immune response [62]. PSME1 encoding PA28 alpha, one of PA28 complex, was up-regulated about 2.8 fold in giant pandas during the breeding season. Furthermore, the proteasomes generate spliced peptides from major histocompatibility complex type I (MHC class I) molecules and PA28 enhances the presentation of several viral epitopes [63]. Proteasome subunits were reported to increase the immune tolerance of the rhesus monkey during early pregnancy [64]. The upregulation of PA28 may indicate the enhancement of immunity in giant pandas.

NK cells comprise 5–10% of lymphocytes in peripheral blood and vary with age [65]. Natural killer (NK) cells play an immensely significant role in innate immunity by defending against virus infections [36]. CD94, encoded by KLRD1, was down-regulated. CD94-NKG2A receptor complex which recognizes MHC class I, is an inhibitor of the cytotoxic activity of NK cells [66]. CD94/NKG2A inhibitory receptor was reported to down-regulate invariant natural killer T cells responses in mice [67]. PLCG2 (PLC-gamma2) which encodes phospholipase C-gamma2 belongs to PLC-gamma proteins family and was up-regulated. PLC-gamma proteins, serving as cytoplasmic enzymes, are involved in NK cell activation [68]. NK cell cytotoxicity was completely abrogated in PLC-gamma2-deficient cells by 4-h 51Cr-release assay [69]. Integration of down-regulated inhibitor and up-regulated activators may imply the partial activation states of NK cells [65]. Shigeru et al. observed the increase of NK cells in pregnant mouse uterus and identified the roles of NK cells in the maintenance of pregnancy [70]. In the testis of macaque and rat, NK cells have an ability for the maintenance of immune privilege and the surveillance for pathological antigens [71]. The increased expression level of NK cells was also reported in breeding horses, suggesting an increased innate immunity of giant pandas during the breeding season [17].

Considered together, these results suggest an enhanced innate immunity in male giant pandas during the breeding season, which is consistent with previous finding [17]. The energy investment in reproduction does not lead to a corresponding decrease innate immune investment. One possible explanation is that the captive pandas are not in a resource-limited environment [3]. Another possible explanation is that the enhancement of innate immunity and cellular immunity can compensate for the decline in the humoral immunity [72, 73].

Adaptive immune changes

The two typical cellular subsets T and B cells comprise the adaptive immune system [74]. In terms of cellular immunity, male ruffs showed a decreased immunity while tree frogs showed an increased immunity during the breeding season via phytohaemagglutinin challenge test [8, 75]. When it comes to humoral immunity, many studies have documented the changes in the immune system. Male bank voles and Eurasian tree sparrows had lower humoral immunocompetence [7, 38], While the immunoglobulin concentration of the Great Tit increased during breeding in accordance with previous studies on birds [16]. In this study, we also observed some changes in cellular immunity and humoral immunity of male giant pandas.

Several key genes involved in antigen presentation and processing were up-regulated. The DIA results also confirmed this finding. The antigen processing pathway is required for proteasomes which produce peptide fragments of MHC class I ligands [63]. The activation of the proteasome relies on PA28 who enhances the liberation of immunopeptidome [63]. Not only was PA28 up-regulated, but also HSP70 and HSP90 were also up-regulated in our study. HSP70 stimulates antigen cross-presentation of dendritic cells and the immune response of activated NK cells [76]. HSP90 contributes to the translocation of extracellular antigen and associates with peptides implicated as precursors of MHC class I ligands [77]. Our data indicates that male giant pandas may have better antigen presentation and processing compared to non-breeding males.

T cell receptors consist of an antigen-binding subunit (TCRαβ) and three dimers of protein CD3 signaling subunit assemble in a coordinated way [78]. CD3D and CD3G coding genes showed elevated transcript levels in the current study, which are involved in TCR activation [79]. Moreover, Aykut et al. also documented that male horses had higher CD3 expression level during the breeding season [17]. The upregulation of IL1R2 was found in T-cell activation [79], and we found that IL1R2 was up-regulated in the male giant panda. However, we found the down-regulation of IL15. IL15 is an important cytokine in lymphocyte survival [79] as well as T cell proliferation and differentiation [73]. Moreover, some co-receptors are indispensable for the activation of T cells [79]. Therefore, T cells did not show increased proliferative and differentiation in male giant pandas during the breeding season.

B cells can differentiate into plasma cells and secrete immunoglobulins against the pathogen [73]. For male temperate bats, reproduction did not influence the concentration of immunoglobulin G (IgG) [6]. However, for male Eurasian tree sparrows, birds during the breeding stage had lower IgA levels [38]. In our study, we found the down-regulation of IGIP. IGIP has the capability of inducing IgA production by B cells [80]. IGIP is primarily produced by dendritic cells and acts as a switch or differentiation factor to regulate IgA [80, 81]. In our study, we found the down-regulation of IGIP. The down-regulation of IGIP may indicate the low concentration of IgA and a reduced humoral immunity in male giant pandas during reproduction.

Conclusions

The present study is the first RNA-seq report on immune profiles of male giant pandas during the breeding season. We identified 45 immune-related genes with altered expression, mostly up-regulated. These genes were related to pattern recognition receptors, autophagy, peroxisome, proteasome, natural killer cell, antigen processing and presentation. Our results suggest an enhanced innate immunity and cellular immunity in male giant pandas during the breeding season, while low humoral immunity was also found. This study provides a foundation for further studies on reproductive immunity in male giant pandas.

Methods

Samples

Peripheral blood samples were collected from 13 captive adult male giant pandas housed at the China Conservation and Research Center for Giant Panda, Sichuan Province, China. The eight male giant pandas peripheral blood samples were collected in April when they were in the breeding season. The five male samples were collected in August when they were in the non-breeding season. A routine physical examination of giant pandas was conducted together with blood sampling. Thirteen pandas with ages ranging from 9 to 11 years old were divided into two groups depending on whether or not they were in the breeding season (Table 1).

Table 1 Information of samples, including identification name, age and group
Full size table

No giant pandas were harmed as a result of this research and continued their normal captive existence after completion of the sampling. Before initiation of this project, the China Conservation and Research Center for Giant Panda had consented to provide blood samples and assisted us in this research.

Library preparation and sequencing

Total RNA from fresh blood was prepared by TRIzol reagent (Invitrogen) and RNeasy kit (Qiagen). We used Nanodrop 8000 Spectrophotometer (Thermo scientific) to evaluate the purity and concentration of RNA. RNA integrity was then checked by using an RNA PicoChip with Agilent 2100 Bioanalyzer (Agilent Technologies). The extracted RNA samples were used for the cDNA synthesis. Double-stranded cDNA was amplified. Sequencing libraries construction, quality control and quantification were performed as per the manufacturer’s instructions using kits from the Illumina Company. The cDNA library was sequenced on the Illumina sequencing platform (HiSeq 2000) using standard procedures. Image analysis and base-calling were performed by the Genome Analyzer Pipeline version 2.0 with default parameters. The 150-bp paired-end reads were generated.

Quantification and mapping

To perform quality control on the data, all raw reads were processed with an adapter trimming and reads filtering by NGS QC Toolkit version 2.3.3 [82]. Quality reports were generated to ensure that clean data were subjected to further analysis by using FastQC version 0.11.5 (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). A genome index was built and reads data were mapped on the giant panda reference genome by using HISAT2 version 2.1.0 [83]. The reference genome (v90) and reference annotation were downloaded from the Ensembl website (ftp://ftp.ensembl.org/pub/release-90/).

Calculating differentially expressed genes (DEGs)

SAM files were generated from alignment program HISAT2 and sorted by SAMtools version 1.7 [84]. SAM files stored mapping information. Sample reads counts file was obtained from BAM files by featureCounts version v1.6.2 [85] which counted the features faster. The expression value of transcripts per million (TPM) was calculated from reads counts. We performed principal components analysis (PCA) on the TPM data matrix of the number of reads by R function prcomp. The PCA plot that showed the clustering information was drawn by R package ggbiplot (https://github.com/vqv/ggbiplot). We then used the counts file as the input to R package edgeR [86]. edgeR performed a method based on the poisson model to infer genes with significant expression differences. The expression fold change (FC) and false discovery rate (FDR) were computed for looking for differentially expressed genes between the two conditions. We set a cut-off of 0.05 for FDR to filter genes for immune-related genes analysis. A cut-off of 1.5 and 0.05 were respectively set for absolute value of log2FC and FDR. The filtered genes were DEGs for further analysis.

Analysis of gene enrichment

To deeply understand the function of these gene sets, genes were annotated and enriched by using the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (http://www.genome.jp/kegg) databases. The GO system includes three categories, molecular function, biological process, and cellular component. A web server g:Profiler [87] was used to cluster genes to GO terms and computed the adjusted P-value (p.adj) by g:SCS algorithm. The KEGG pathway map analysis, presenting biological interpretation of higher-level systemic functions, was performed by KOBAS 2.0 [88]. A threshold of 0.05 was set for p.adj.

To understand the impact and direction of the pathways, KEGG enrichment analysis was also performed by the Dynamic Impact Approach [89]. The estimate of the dynamic change was represented by the “Impact”. The overall direction of the dynamic change was represented by the “Flux”. Detailed methodology has been reported previously [89]. The data set included Entrez Gene ID, FDR, expression ratio, and P-value. We uploaded the data set and chose Human as annotation species. An FDR < 0.05 and a P-value < 0.05 between the breeding and non-breeding were used as cut-off.

Analysis of protein-protein interaction network

STRING [90] is a database of known and predicted protein-protein interactions. We input immune-related genes into STRING to obtain protein-protein interaction networks. The output of protein-protein interactions was simple texts in tabular form. Cytoscape [91] was used to visualize molecular interaction. We used cytoHubba [92] which is the plugin of Cytoscape to calculate hub genes. The genes were ranked by radiality algorithm. Three additional topological algorithms of betweeness, closeness, and degree were used for verification [93].

Real-time quantitative PCR (qRT-PCR)

We performed the qRT-PCR to confirm the expression changes of genes. A dozen of top-ranked up-regulated and down-regulated genes were selected. The reference gene was GAPDH, which was same as previous studies [73]. We had six samples (B2, B5, B6, N1, N2, N3) because several samples were lost. The sequences of primers were predicted using Primer-BLAST [94] and listed in Table 2. Primers were synthesized at TSINGKE. Each qPCR reaction system was 10 μL containing 5 μL of 2x M5 Hiper SYBR Premix EsTaq (mei5, Beijing, China), 0.2 μL cDNA, 0.2 μL of forward and reverse primer and nuclease-free water. All reactions were performed on Bio-Rad CFX96 Touch in triplicates. The following program was used: 95 °C for 30 s, followed by 40 cycles of 95 °C/10 s, 60 °C/15 s; then 72 °C for 10 s. The optimized comparative Ct (2-ΔΔCt) value method was applied to calculate gene expression levels. We used GraphPad Prism 8.0.2 software to analyze data.

Table 2 The primer information of qRT-PCR
Full size table

Availability of data and materials

Raw sequence data are accessible at NCBI under the BioProject accession number PRJNA631846 (https://www.ncbi.nlm.nih.gov/sra/PRJNA631846). The giant panda reference genome (v90) (http://ftp.ensembl.org/pub/release-90/fasta/ailuropoda_melanoleuca/) and reference annotation (http://ftp.ensembl.org/pub/release-90/gtf/ailuropoda_melanoleuca/) were downloaded from the Ensembl website. The accession numbers in Additional file 1 Table S1 and Additional file 2 Table S2 correspond to Ensembl website (https://asia.ensembl.org/index.html).

Abbreviations

DEGs:

Differentially expressed genes

TPM:

Transcripts per million

PCA:

Principal components analysis

FC:

Fold change

FDR:

False discovery rate

GO:

Gene Ontology

KEGG:

Kyoto Encyclopedia of Genes and Genomes

PRRs:

Pattern recognition receptors

TLRs:

Toll-like receptors

CLRs:

C-type lectin receptors

NLRs:

NOD-like receptors

MHC class I:

Major histocompatibility complex type I

NK:

Natural killer

IgG:

Immunoglobulin G

References

  1. 1.

    Graña Grilli M, Pari M, Ibañez A. Poor body conditions during the breeding period in a seabird population with low breeding success. Mar Biol. 2018;165(9):142.

    Article  Google Scholar 

  2. 2.

    Stearns SC. The evolution of life histories: OUP Oxford; 1992.

    Google Scholar 

  3. 3.

    French SS, Moore MC. Immune function varies with reproductive stage and context in female and male tree lizards, Urosaurus ornatus. Gen Comp Endocr. 2008;155(1):148–56.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  4. 4.

    Segner H, Verburg-van Kemenade BML, Chadzinska M. The immunomodulatory role of the hypothalamus-pituitary-gonad axis: proximate mechanism for reproduction-immune trade offs? Developmental & Comparative Immunology. 2017;66:43–60.

    CAS  Article  Google Scholar 

  5. 5.

    Schwenke RA, Lazzaro BP, Wolfner MF. Reproduction-immunity trade-offs in insects. Annu Rev Entomol. 2016;61(1):239–56.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  6. 6.

    Ruoss S, Becker NI, Otto MS, Czirják GÁ, Encarnação JA. Effect of sex and reproductive status on the immunity of the temperate bat Myotis daubentonii. Mamm Biol. 2019;94(1):120–6.

    PubMed  Article  PubMed Central  Google Scholar 

  7. 7.

    Saino N, Canova L, Fasola M, Martinelli R. Reproduction and population density affect humoral immunity in bank voles under field experimental conditions. OECOLOGIA. 2000;124(3):358–66.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  8. 8.

    Desprat JL, Lengagne T, Dumet A, Desouhant E, Mondy N. Immunocompetence handicap hypothesis in tree frog: trade-off between sexual signals and immunity? Behav Ecol. 2015;26(4):1138–46.

    Article  Google Scholar 

  9. 9.

    Gupta V, Ali ZS, Prasad NG. Sexual activity increases resistance against Pseudomonas entomophila in male Drosophila melanogaster. BMC Evol Biol. 2013;13(1):185.

    PubMed  PubMed Central  Article  Google Scholar 

  10. 10.

    Wei F, Hu Y, Yan L, Nie Y, Wu Q, Zhang Z. Giant pandas are not an evolutionary cul-de-sac: evidence from multidisciplinary research. Mol Biol Evol. 2015;32(1):4–12.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  11. 11.

    Wu W, Wu H, He M, Zhang L, Huang Y, Geng Y, Liu J, Wang Q, Fan Z, Hou R, et al. Transcriptome analyses provide insights into maternal immune changes at several critical phases of giant panda reproduction. Dev Comp Immunol. 2020;110:103699.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  12. 12.

    Martin-Wintle MS, Kersey DC, Wintle N, Aitken-Palmer C, Owen MA, Swaisgood RR. Comprehensive breeding techniques for the Giant panda. Adv Exp Med Biol. 2019;1200:275–308.

    PubMed  Article  PubMed Central  Google Scholar 

  13. 13.

    Aitken-Palmer C, Hou R, Burrell C, Zhang Z, Wang C, Spindler R, Wildt DE, Ottinger MA, Howard J. Protracted reproductive seasonality in the male giant panda (Ailuropoda melanoleuca) reflected by patterns in androgen profiles, ejaculate characteristics, and selected behaviors. Biol Reprod. 2012;86(6):195.

    PubMed  Article  CAS  PubMed Central  Google Scholar 

  14. 14.

    Zhou W, Nie Y, Hu Y, Swaisgood RR, Zhang Y, Liu D, Wei F. Seasonal and reproductive variation in chemical constituents of scent signals in wild giant pandas. Sci China Life Sci. 2019;62(5):648–60.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  15. 15.

    Charlton BD, Owen MA, Zhou X, Zhang H, Swaisgood RR. Influence of season and social context on male giant panda (Ailuropoda melanoleuca) vocal behaviour. PLoS One. 2019;14(11):e225772.

    Article  CAS  Google Scholar 

  16. 16.

    Pap PL, Vágási CI, Tökölyi J, Czirják GÁ, Barta Z. Variation in Haematological indices and immune function during the annual cycle in the great tit Parus major. ARDEA. 2010;98(1):105–12.

    Article  Google Scholar 

  17. 17.

    Uner AG, Sulu N, Altinsaat C, Ergun A. Blood levels of selected metabolic factors, cytokines, and lymphocyte subpopulations in Arabian and thoroughbred horses during the longest and shortest days of the year. J EQUINE VET SCI. 2013;33(11):969–76.

    Article  Google Scholar 

  18. 18.

    O'Brien KA, Waterman JM, Anderson WG, Bennett NC. Trade-offs between immunity and testosterone in male African ground squirrels. J Exp Biol. 2018;221(16):b177683.

    Article  Google Scholar 

  19. 19.

    Millet S, Bennett J, Lee KA, Hau M, Klasing KC. Quantifying and comparing constitutive immunity across avian species. Developmental & Comparative Immunology. 2007;31(2):188–201.

    CAS  Article  Google Scholar 

  20. 20.

    Cutrera AP, Zenuto RR, Luna F, Antenucci CD. Mounting a specific immune response increases energy expenditure of the subterranean rodent Ctenomys talarum (tuco-tuco): implications for intraspecific and interspecific variation in immunological traits. J Exp Biol. 2010;213(5):715–24.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  21. 21.

    Chaussabel D, Pascual V, Banchereau J. Assessing the human immune system through blood transcriptomics. BMC Biol. 2010;8(1):84.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  22. 22.

    Kiss J, Rádai Z, Rosa ME, Kosztolányi A, Barta Z. Seasonal changes in immune response and reproductive investment in a biparental beetle. J Insect Physiol. 2020;121:104000.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  23. 23.

    Gocinski BL, Knott KK, Roberts BM, Brown JL, Vance CK, Kouba AJ. Changes in urinary androgen concentration indicate that male giant pandas (Ailuropoda melanoleuca) respond to impending female oestrus during and outside the typical spring breeding season. Reprod Fertil Dev. 2018;30(2):399–408.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  24. 24.

    Roy A, Basak R, Rai U. De novo sequencing and comparative analysis of testicular transcriptome from different reproductive phases in freshwater spotted snakehead Channa punctatus. PLoS One. 2017;12(3):e173178.

    Article  CAS  Google Scholar 

  25. 25.

    Bai M, Sun L, Zhao J, Xiang L, Cheng X, Li J, Jia C, Jiang H. Histological analysis and identification of spermatogenesis-related genes in 2-, 6-, and 12-month-old sheep testes. SCI NAT-HEIDELBERG. 2017;104(9–10):84.

    Article  CAS  Google Scholar 

  26. 26.

    Jin S, Hu Y, Fu H, Sun S, Jiang S, Xiong Y, Qiao H, Zhang W, Gong Y, Wu Y. Analysis of testis metabolome and transcriptome from the oriental river prawn (Macrobrachium nipponense) in response to different temperatures and illumination times. Comp Biochem Physiol Part D Genomics Proteomics. 2020;34:100662.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  27. 27.

    Guillaume D, Moussu C, de Geoffroy F, Chesneau D, Keller M. Olfactory stimulation or inhibition of sexual behavior of stallions in non-breeding season. Physiol Behav. 2018;186:1–9.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  28. 28.

    Torella JP, Ford TJ, Kim SN, Chen AM, Way JC, Silver PA. Tailored fatty acid synthesis via dynamic control of fatty acid elongation. Proc Natl Acad Sci. 2013;110(28):11290.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  29. 29.

    Samadian F, Towhidi A, Rezayazdi K, Bahreini M. Effects of dietary n-3 fatty acids on characteristics and lipid composition of ovine sperm. ANIMAL. 2010;4(12):2017–22.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  30. 30.

    Aurich C, Ortega FC, Peña VF, Schrammel N, Morcuende D, Aurich J. Seasonal changes in the sperm fatty acid composition of Shetland pony stallions. THERIOGENOLOGY. 2018;107:149–53.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  31. 31.

    Huang Z, Rivas B, Agoulnik AI. NOTCH1 gain of function in germ cells causes failure of spermatogenesis in male mice. PLoS One. 2013;8(7):e71213.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  32. 32.

    Garcia TX, Hofmann M. NOTCH signaling in Sertoli cells regulates gonocyte fate. Cell cycle (Georgetown, Tex.). 2013;12(16):2538–45.

    CAS  Article  Google Scholar 

  33. 33.

    Di R, He J, Song S, Tian D, Liu Q, Liang X, Ma Q, Sun M, Wang J, Zhao W, et al. Characterization and comparative profiling of ovarian microRNAs during ovine anestrus and the breeding season. BMC Genomics. 2014;15(1):899.

    PubMed  PubMed Central  Article  Google Scholar 

  34. 34.

    Held T, Barakat AZ, Mohamed BA, Paprotta I, Meinhardt A, Engel W, Adham IM. Heat-shock protein HSPA4 is required for progression of spermatogenesis. REPRODUCTION. 2011;142(1):133–44.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  35. 35.

    Koziorowska-Gilun M, Gilun P, Mietelska K, Kordan W. Determination of the activity and relative abundance of mRNA for antioxidant enzymes in stallion testicular and epididymal tissues: a comparison between two breeding seasons. Anim Reprod Sci. 2018;196:230–8.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  36. 36.

    Nicholson LB. The immune system. Essays Biochem. 2016;60(3):275–301.

    PubMed  PubMed Central  Article  Google Scholar 

  37. 37.

    Mair KH, Sedlak C, Käser T, Pasternak A, Levast B, Gerner W, Saalmüller A, Summerfield A, Gerdts V, Wilson HL, et al. The porcine innate immune system: an update. Dev Comp Immunol. 2014;45(2):321–43.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  38. 38.

    Zhao Y, Li M, Sun Y, Wu W, Kou G, Guo L, Xing D, Wu Y, Li D, Zhao B. Life-history dependent relationships between body condition and immunity, between immunity indices in male Eurasian tree sparrows. Comp Biochem Physiol A Mol Integr Physiol. 2017;210:7–13.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  39. 39.

    Takeuchi O, Akira S. Pattern recognition receptors and inflammation. CELL. 2010;140(6):805–20.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  40. 40.

    Lu X, Nagata M, Yamasaki S. Mincle: 20 years of a versatile sensor of insults. Int Immunol. 2018;30(6):233–9.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  41. 41.

    Kostarnoy AV, Gancheva PG, Lepenies B, Tukhvatulin AI, Dzharullaeva AS, Polyakov NB, Grumov DA, Egorova DA, Kulibin AY, Bobrov MA, et al. Receptor Mincle promotes skin allergies and is capable of recognizing cholesterol sulfate. Proc Natl Acad Sci. 2017;114(13):E2758–65.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  42. 42.

    Van Gorp H, Kuchmiy A, Van Hauwermeiren F, Lamkanfi M. NOD-like receptors interfacing the immune and reproductive systems. FEBS J. 2014;281(20):4568–82.

    PubMed  Article  CAS  PubMed Central  Google Scholar 

  43. 43.

    Kufer TA, Sansonetti PJ. NLR functions beyond pathogen recognition. Nat Immunol. 2011;12(2):121–8.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  44. 44.

    Da SCJ, Miranda Y, Leonard N, Ulevitch R. SGT1 is essential for Nod1 activation. Proc Natl Acad Sci U S A. 2007;104(16):6764–9.

    Article  Google Scholar 

  45. 45.

    Mayor A, Martinon F, De Smedt T, Pétrilli V, Tschopp J. A crucial function of SGT1 and HSP90 in inflammasome activity links mammalian and plant innate immune responses. Nat Immunol. 2007;8(5):497–503.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  46. 46.

    Matsuzawa-Ishimoto Y, Hwang S, Cadwell K. Autophagy and inflammation. Annu Rev Immunol. 2018;36:73–101.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  47. 47.

    Ocón-Grove OM, Krzysik-Walker SM, Maddineni SR, Hendricks GR, Ramachandran R. NAMPT (visfatin) in the chicken testis: influence of sexual maturation on cellular localization, plasma levels and gene and protein expression. REPRODUCTION. 2010;139(1):217–26.

    PubMed  Article  CAS  PubMed Central  Google Scholar 

  48. 48.

    Annie L, Gurusubramanian G, Kumar RV. Dexamethasone mediated downregulation of PGC-1α and visfatin regulates testosterone synthesis and antioxidant system in mouse testis. Acta Histochem. 2019;121(2):182–8.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  49. 49.

    Cameron AM, Castoldi A, Sanin DE, Flachsmann LJ, Field CS, Puleston DJ, Kyle RL, Patterson AE, Hässler F, Buescher JM, et al. Inflammatory macrophage dependence on NAD + salvage is a consequence of reactive oxygen species-mediated DNA damage. Nat Immunol. 2019;20(4):420–32.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  50. 50.

    Kuo CJ, Hansen M, Troemel E. Autophagy and innate immunity: insights from invertebrate model organisms. AUTOPHAGY. 2018;14(2):233–42.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  51. 51.

    Tabecka-Lonczynska A, Mytych J, Solek P, Koziorowski M. Autophagy as a consequence of seasonal functions of testis and epididymis in adult male European bison (Bison bonasus, Linnaeus 1758). Cell Tissue Res. 2020;379(3):613–24.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  52. 52.

    Shibutani ST, Saitoh T, Nowag H, Münz C, Yoshimori T. Autophagy and autophagy-related proteins in the immune system. Nat Immunol. 2015;16(10):1014–24.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  53. 53.

    González CR, Muscarsel IM, Vitullo AD. The balance between apoptosis and autophagy regulates testis regression and recrudescence in the seasonal-breeding south American plains vizcacha, Lagostomus maximus. PLOS ONE. 2018;13(1):e191126.

    Google Scholar 

  54. 54.

    Ferreira AR, Marques M, Ribeiro D. Peroxisomes and Innate Immunity: Antiviral Response and Beyond. Int J Mol Sci. 2019:20(15).

  55. 55.

    Di Cara F, Sheshachalam A, Braverman NE, Rachubinski RA, Simmonds AJ. Peroxisome-mediated metabolism is required for immune response to microbial infection. IMMUNITY. 2017;47(1):93–106.

    PubMed  Article  CAS  PubMed Central  Google Scholar 

  56. 56.

    Nenicu A, Lüers GH, Kovacs W, Bergmann M, Baumgart-Vogt E. Peroxisomes in human and mouse testis: differential expression of Peroxisomal proteins in germ cells and distinct somatic cell types of the Testis1. Biol Reprod. 2007;77(6):1060–72.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  57. 57.

    Jung KI, Pyo CW, Choi SY. Influenza a virus-induced autophagy contributes to enhancement of virus infectivity by SOD1 downregulation in alveolar epithelial cells. Biochem Biophys Res Commun. 2018;498(4):960–6.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  58. 58.

    Aitken RJ. Reactive oxygen species as mediators of sperm capacitation and pathological damage. Mol Reprod Dev. 2017;84(10):1039–52.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  59. 59.

    Kim WJ, Lee SE, Park YG, Jeong SG, Kim EY, Park SP. Antioxidant hesperetin improves the quality of porcine oocytes during aging in vitro. Mol Reprod Dev. 2019;86(1):32–41.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  60. 60.

    Knoops B, Becker S, Poncin MA, Glibert J, Derclaye S, Clippe A, Alsteens D. Specific interactions measured by AFM on living cells between Peroxiredoxin-5 and TLR4: relevance for mechanisms of innate immunity. CELL CHEM BIOL. 2018;25(5):550–9.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  61. 61.

    Kammerl IE, Meiners S. Proteasome function shapes innate and adaptive immune responses. Am J Physiol Lung Cell Mol Physiol. 2016;311(2):L328–36.

    PubMed  Article  PubMed Central  Google Scholar 

  62. 62.

    Feist E, Burmester GR, Krüger E. The proteasome - victim or culprit in autoimmunity. Clin Immunol. 2016;172:83–9.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  63. 63.

    Sijts A, Sun Y, Janek K, Kral S, Paschen A, Schadendorf D, Kloetzel PM. The role of the proteasome activator PA28 in MHC class I antigen processing. Mol Immunol. 2002;39(3–4):165–9.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  64. 64.

    Zhao H, Sui L, Miao K, An L, Wang D, Hou Z, Wang R, Guo M, Wang Z, Xu J, et al. Comparative analysis between endometrial proteomes of pregnant and non-pregnant ewes during the peri-implantation period. J ANIM SCI BIOTECHNO. 2015;6(1):18.

    Article  CAS  Google Scholar 

  65. 65.

    Pegram HJ, Andrews DM, Smyth MJ, Darcy PK, Kershaw MH. Activating and inhibitory receptors of natural killer cells. Immunol Cell Biol. 2011;89(2):216–24.

    PubMed  Article  PubMed Central  Google Scholar 

  66. 66.

    Mesecke S, Urlaub D, Busch H, Eils R, Watzl C. Integration of activating and inhibitory receptor signaling by regulated phosphorylation of Vav1 in immune cells. Sci Signal. 2011;4(175):a36.

    Article  CAS  Google Scholar 

  67. 67.

    Ota T, Takeda K, Akiba H, Hayakawa Y, Ogasawara K, Ikarashi Y, Miyake S, Wakasugi H, Yamamura T, Kronenberg M, et al. IFN-gamma-mediated negative feedback regulation of NKT-cell function by CD94/NKG2. BLOOD. 2005;106(1):184–92.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  68. 68.

    Matalon O, Fried S, Ben-Shmuel A, Pauker MH, Joseph N, Keizer D, Piterburg M, Barda-Saad M. Dephosphorylation of the adaptor LAT and phospholipase C-γ by SHP-1 inhibits natural killer cell cytotoxicity. Sci Signal. 2016;9(429):a54.

    Article  CAS  Google Scholar 

  69. 69.

    Caraux A, Kim N, Bell SE, Zompi S, Ranson T, Lesjean-Pottier S, Garcia-Ojeda ME, Turner M, Colucci F. Phospholipase C-gamma2 is essential for NK cell cytotoxicity and innate immunity to malignant and virally infected cells. BLOOD. 2006;107(3):994–1002.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  70. 70.

    Saito S, Shima T. Nakashima a: [regulatory T cells and regulatory NK cells play essential roles for maintenance of pregnancy]. Nihon Rinsho Meneki Gakkai Kaishi. 2012;35(5):424–8.

    CAS  PubMed  Article  Google Scholar 

  71. 71.

    Duan YG, Gong J, Yeung W, Haidl G, Allam JP. Natural killer and NKT cells in the male reproductive tract. J Reprod Immunol. 2020;142:103178.

    CAS  PubMed  Article  Google Scholar 

  72. 72.

    Du L, Liu Q, Shen F, Fan Z, Hou R, Yue B, Zhang X. Transcriptome analysis reveals immune-related gene expression changes with age in giant panda (Ailuropoda melanoleuca) blood. Aging. 2019;11(1):249–62.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  73. 73.

    Geng Y, Shen F, Wu W, Zhang L, Luo L, Fan Z, Hou R, Yue B, Zhang X. First demonstration of giant panda's immune response to canine distemper vaccine. Developmental & Comparative Immunology. 2020;102:103489.

    CAS  Article  Google Scholar 

  74. 74.

    den Haan JMM, Arens R, van Zelm MC. The activation of the adaptive immune system: Cross-talk between antigen-presenting cells, T cells and B cells. Immunol Lett. 2014;162(2, Part B):103–12.

    Article  CAS  Google Scholar 

  75. 75.

    Lozano GA, Lank DB. Seasonal trade-offs in cell-mediated immunosenescence in ruffs (Philomachus pugnax). Proceedings of the Royal Society of London. Series B: Biological Sciences. 2003;270(1520):1203–8.

    Google Scholar 

  76. 76.

    Multhoff G, Pockley AG, Schmid TE, Schilling D. The role of heat shock protein 70 (Hsp70) in radiation-induced immunomodulation. Cancer Lett. 2015;368(2):179–84.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  77. 77.

    Imai T, Kato Y, Kajiwara C, Mizukami S, Ishige I, Ichiyanagi T, Hikida M, Wang J, Udono H. Heat shock protein 90 (HSP90) contributes to cytosolic translocation of extracellular antigen for cross-presentation by dendritic cells. P NATL ACAD SCI USA. 2011;108(39):16363–8.

    CAS  Article  Google Scholar 

  78. 78.

    Guy CS, Vignali KM, Temirov J, Bettini ML, Overacre AE, Smeltzer M, Zhang H, Huppa JB, Tsai Y, Lobry C, et al. Distinct TCR signaling pathways drive proliferation and cytokine production in T cells. Nat Immunol. 2013;14(3):262–70.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  79. 79.

    Wang M, Windgassen D, Papoutsakis ET. Comparative analysis of transcriptional profiling of CD3+, CD4+ and CD8+ T cells identifies novel immune response players in T-cell activation. BMC Genomics. 2008;9:225.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  80. 80.

    Austin AS, Haas KM, Naugler SM, Bajer AA, Garcia-Tapia D, Estes DM. Identification and characterization of a novel regulatory factor: IgA-inducing protein. J Immunol. 2003;171(3):1336–42.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  81. 81.

    Endsley MA, Njongmeta LM, Shell E, Ryan MW, Indrikovs AJ, Ulualp S, Goldblum RM, Mwangi W, Estes DM. Human IgA-inducing protein from dendritic cells induces IgA production by naive IgD+ B cells. Journal of immunology (Baltimore, Md.: 1950). 2009;182(4):1854–9.

    CAS  Article  Google Scholar 

  82. 82.

    Patel RK, Jain M. NGS QC toolkit: a toolkit for quality control of next generation sequencing data. PLoS One. 2012;7(2):e30619.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  83. 83.

    Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12(4):357–60.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  84. 84.

    Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R. The sequence alignment/map format and SAMtools. BIOINFORMATICS. 2009;25(16):2078–9.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  85. 85.

    Liao Y, Smyth GK. Shi W: featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014;30(7):923–30.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  86. 86.

    Robinson MD, McCarthy DJ. Smyth GK: edgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–40.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  87. 87.

    Reimand J, Arak T, Adler P, Kolberg L, Reisberg S, Peterson H. Vilo J: g:profiler-a web server for functional interpretation of gene lists (2016 update). Nucleic Acids Res. 2016;44(W1):W83–9.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  88. 88.

    Xie C, Mao X, Huang J, Ding Y, Wu J, Dong S, Kong L, Gao G, Li CY, Wei L. KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res. 2011;39(Web Server issue):W316–22.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  89. 89.

    Bionaz M, Periasamy K, Rodriguez-Zas SL, Hurley WL, Loor JJ. A novel dynamic impact approach (DIA) for functional analysis of time-course omics studies: validation using the bovine mammary transcriptome. PLoS One. 2012;7(3):e32455.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  90. 90.

    Szklarczyk D, Gable AL, Lyon D, Junge A, Wyder S, Huerta-Cepas J, Simonovic M, Doncheva NT, Morris JH, Bork P, et al. STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. 2019;47(D1):D607–13.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  91. 91.

    Smoot ME, Ono K, Ruscheinski J, Wang PL, Ideker T. Cytoscape 2.8: new features for data integration and network visualization. BIOINFORMATICS. 2011;27(3):431–2.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  92. 92.

    Chin CH, Chen SH, Wu HH, Ho CW, Ko MT, Lin CY. cytoHubba: identifying hub objects and sub-networks from complex interactome. Bmc Syst Biol. 2014;8(Suppl 4):S11.

    PubMed  PubMed Central  Article  Google Scholar 

  93. 93.

    Anupama R, Sajitha LS, Mukherjee A, Babu S. Cross-regulatory network in Pseudomonas aeruginosa biofilm genes and TiO2 anatase induced molecular perturbations in key proteins unraveled by a systems biology approach. Gene. 2018;647:289–96.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  94. 94.

    Ye J, Coulouris G, Zaretskaya I, Cutcutache I, Rozen S, Madden TL. Primer-BLAST: a tool to design target-specific primers for polymerase chain reaction. BMC Bioinformatics. 2012;13:134.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

Download references

Acknowledgments

We acknowledge Dr. Megan Price for providing language help.

Funding

This work was supported by the Open Project of Key Laboratory of State Forestry and grassland administration (KLSFGAGP2020.010) and Project of Key Laboratory of State Forestry and grassland administration (CCRCGP181913). The funder had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Author information

Author notes

  1. Haibo Shen and Caiwu Li contributed equally to this work.

Affiliations

  1. Key Laboratory of Bio-resources and Eco-environment, Ministry of Education, College of Life Science, Sichuan University, No. 24 South Section 1, Yihuan Road, Chengdu, 610065, Sichuan, China

    Haibo Shen & Xiuyue Zhang

  2. Key Laboratory of State Forestry and Grassland Administration on Conservation Biology of Rare Animals in The Giant Panda National Park, China Conservation and Research Center for the Giant Panda, Dujiangyan, 611830, Sichuan, PR China

    Caiwu Li, Ming He, Yan Huang, Jing Wang & Minglei Wang

  3. Sichuan Key Laboratory of Conservation Biology on Endangered Wildlife, College of Life Sciences, Sichuan University, Chengdu, 610064, PR China

    Bisong Yue

Authors
  1. Haibo Shen
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Caiwu Li
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Ming He
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Yan Huang
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Jing Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Minglei Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Bisong Yue
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Xiuyue Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

HBS and CWL coordinated and performed the research. HBS analyzed the data, prepared all figures and wrote the manuscript. MH, YH, JW, and MLW provided the blood samples, and performed qRT-PCR validation. BSY and XYZ designed the research. All authors have read and agreed to the published version of the manuscript.

Corresponding author

Correspondence to Xiuyue Zhang.

Ethics declarations

Ethics approval and consent to participate

This study was carried out according to the Regulation on the Administration of Laboratory Animals (2017 Revision) published by the Ministry of Science and Technology of the People’s Republic of China. All study procedures and animal care activities were approved by the Academic and Ethics Committee of Sichuan University.

Consent for publication

Not applicable.

Competing interests

The authors declare no conflict of interest.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Table S1.

Differentially expressed genes between the breeding season and the non-breeding season. 1128 genes changed in expression level were detected with a given FDR threshold. The table contains three columns, including GENEID, log2FC and FDR.

Additional file 2: Table S2.

GO enrichment of differentially expressed genes. Up-regulated DEGs were enriched in 69 GO terms and down-regulated DEGs were enriched in 18 GO terms.

Additional file 3: Table S3.

The score of 36 genes calculated by topological analysis methods. The table contains five columns, including Gene Symbol, Radiality, Betweenness, Closeness, and Degree.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Shen, H., Li, C., He, M. et al. Immune profiles of male giant panda (Ailuropoda melanoleuca) during the breeding season. BMC Genomics 22, 143 (2021). https://doi.org/10.1186/s12864-021-07456-x

Download citation

Keywords

  • Male giant panda
  • Immune change
  • Breeding season
  • RNA-seq

BMC Genomics

ISSN: 1471-2164

Contact us