Characterization of avian β-defensin genes in Galliformes reveals widespread evolutionary diversification and distinct evolutionary relationships with infection risk

Abstract

Background

Avian β-defensins (AvBDs) represent a key family of antimicrobial host defense peptides in birds. Accumulating evidence suggests that the evolutionary trajectory of β-defensin genes is specific to the gene, timescale, and species involved, implying that species-specific ecological and life-history differences drive divergent selective pressures on these genes. However, their evolutionary dynamics, particularly the interactions with ecological factors and life-history traits, remain insufficiently explored.

Results

Through a comprehensive survey of 25 species spanning all major clades of Galliformes, 354 AvBD genes were identified. Comparative sequence analysis, genomic organization, and phylogenetic studies collectively reveal significant evolutionary diversification characterized by gene duplication, pseudogenization, and gene loss across these species. Notably, chicken AvBD3 exhibits significant differences in its coding regions, while AvBD6 and AvBD7 appear to have copy number variations, with species-specific paralogs of AvBD6 being especially prominent. Moreover, positive selection was more frequently observed in recently diverged gene lineages compared to ancestral ones. Using 70 samples from eight galliform species, the study further identified the prevalence of species-specific amino acid alleles. Phylogenetic comparative analysis demonstrated that the evolution of nine AvBD genes (AvBD2, -4, -5, -8, -9, -10, -11, -12, and -14) is significantly associated with specific ecological factors and life-history characteristics. Additionally, the evolutionary rates of these genes showed distinct relationship with inferred infection risk, likely reflecting the multifunctionality of β-defensins and potential trade-offs between immune defense and other biological functions.

Conclusions

This cross-species identification and systematic evolutionary analysis of AvBDs in Galliformes deepen our understanding of the co-evolution of host defense peptides, offering valuable insights into their natural biology and evolution, and paving the way for future applications as alternatives to traditional antibiotics.

Background

Host defense peptides, also referred to as antimicrobial peptides, comprise a group of cationic amphipathic short peptides ubiquitously expressed across multicellular eukaryotes. These molecules play a pivotal role in immune defense, exhibiting potent antimicrobial activity against a wide array of pathogens, including viruses, bacteria, fungi, protozoa, and parasitic microbes [1, 2]. Beyond their broad-spectrum antimicrobial functions, recent studies have highlighted additional immunomodulatory roles for host defense peptides [3, 4]. Their strong antimicrobial efficacy, reduced environmental persistence, multifunctionality within the immune system, and lower likelihood of inducing bacterial resistance make them promising candidates for alternatives to conventional antibiotics [5,6,7]. However, emerging evidence underscores the necessity for a deeper foundational understanding of these peptides to avoid replicating the resistance crisis currently plaguing traditional antibiotics [8, 9]. Among the critical areas of focus, systematically characterizing their evolutionary history is deemed essential for enhancing their efficacy and sustainability in translational applications, as well as for predicting their future utility [10].

Defensins, a prominent family of cysteine-rich host defense peptides, are present in fungi, plants, and animals. They are categorized into two major subfamilies—cis- and trans-defensins—based on secondary structure arrangement, tertiary conformation, and disulfide bond connectivity [11]. In vertebrates, most defensins are further subdivided into four families: α-, β-, θ- and ovo-defensins [12,13,14]. Of these, only β-defensins and ovo-defensins have been identified in birds. Avian β-defensins (AvBDs), initially identified in chicken leukocytes [15], are characterized by a six-cysteine motif, C-X_4–6-C-X_3–5-C-X_9–10-C-X_5–6-CC, where Xn indicates the number of amino acids between cysteines. Genomic analyses reveal that avian genomes encode multiple β-defensin genes across 14 gene lineages (AvBD1-14), flanked by CTSB and TRAM2 [16]. These genes are densely clustered within a complex genomic region ranging from 66 kb to approximately 200 kb, and are prone to structural variation [17, 18].

Despite their ancient phylogenetic origins, β-defensins remain evolutionarily dynamic. Instances of rapid evolution, including positive selection, extensive copy number variation, and duplication with subsequent divergence, have been observed across a broad spectrum of vertebrates [19,20,21,22]. These events have led to amino acid substitutions, alterations in gene expression, and sequence variations, all of which can significantly impact their biological functions, particularly regarding microbial specificity [18, 23]. Moreover, accumulating evidence suggests that the evolutionary trajectory of β-defensin genes is specific to the gene, timescale, and species involved [24,25,26], implying that species-specific ecological and life-history differences drive divergent selective pressures on these genes. Ecological traits and life-history factors associated with infection risk have recently been recognized as major drivers of immune-related gene evolution [27,28,29]. Additionally, the life history trade-off hypothesis posits that trade-offs between immunity and fitness costs shape species-specific balances between immune function and other biological processes [30, 31]. These frameworks provide insights into how pathogen-mediated selection shapes the evolution of defensins.

Although the galliform order is nearly globally distributed, individual species tend to inhabit restricted regions, resulting in a wide range of ecological characteristics and life-history traits [32, 33]. Furthermore, the availability of complete genomes and related sequences for these species enables comprehensive characterization of the AvBD gene family. This presents an excellent opportunity to investigate the evolutionary relationships of AvBDs and to explore whether their evolution is shaped by potential pathogen exposure and associated fitness costs. By integrating systematic evolutionary analyses, this study characterized the full repertoire of AvBD genes from 25 species spanning major clades of Galliformes. It also investigated the relationship between the evolution of avian β-defensins and key ecological factors, along with life-history traits related to infection risk and host investment. Our findings reveal widespread gene duplication, pseudogenization, and gene loss across galliform genomes, along with amino acid substitutions and other variations arising during species divergence. Notably, most positively selected sites were located in the mature peptide region, with these sites appearing more frequently in recently evolved gene lineages than in older ones. Additionally, our results indicate that the evolutionary rates of AvBD genes display divergent correlation patterns with ecological factors and life-history traits, suggesting distinct evolutionary relationships with infection risk.

Materials and methods

Genomic identification and full-length coding sequences prediction

To identify homologous sequences of the AvBD gene family, the amino acid sequences of the 14 chicken AvBD mature peptides (AvBD1-14) were used as reference queries (Additional file 1). TBLASTN searches were performed against the genomes and available Sequence Read Archives (SRAs) of the corresponding species, as well as SRAs from representative tissues of 24 Galliformes species (last accessed: January 19, 2025). Notably, due to the lack of support for the 3-coding-exon structure of AvBD13 in chickens [34], and in many other birds [12] and reptiles [26], only the amino acids encoded by the first two exons of chicken AvBD13 (NM_001001780.1) were used for homology identification in this study. The 24 species included the red-legged partridge (Alectoris rufa), Chinese bamboo-partridge (Bambusicola thoracicus), Indian peafowl (Pavo cristatus), green peafowl (Pavo muticus), Gunnison grouse (Centrocercus minimus), greater sage-grouse (Centrocercus urophasianus), greater prairie chicken (Tympanuchus cupido), lesser prairie-chicken (Tympanuchus pallidicinctus), white-tailed ptarmigan (Lagopus leucura), rock ptarmigan (Lagopus muta), turkey (Meleagris gallopavo), brown eared-pheasant (Crossoptilon mantchuricum), silver pheasant (Lophura nycthemera), ring-necked pheasant (Phasianus colchicus), Mikado pheasant (Syrmaticus mikado), California quail (Callipepla californica), scaled quail (Callipepla squamata), northern bobwhite (Colinus virginianus), marbled wood quail (Odontophorus gujanensis), helmeted guineafowl (Numida meleagris), white-crested guan (Penelope pileata), Australian brush-turkey (Alectura latham), as well as those of previously studied species such as the golden pheasant (Chrysolophus pictus), and Japanese quail (Coturnix japonica) (Additional file 2). For each defensin sequence identified in the genomes, 5,000-bp sequences upstream and downstream were retrieved to predict the full-length coding sequences using SignalP 6.0 [35] and GeneWise [36] with default settings, and subsequently adjusted based on known transcript sequences. Sequences containing mutations of the initiation codon, premature stop codons or frameshift mutations (insertions or deletions) before the last conserved cysteine were classified as pseudogenes, while sequences missing any part of the typical AvBD gene structure were considered partial genes. AvBD genes from different species were named according to a nomenclature system based on the first two letters of the genus and species (Additional file 2), as previously proposed [37]. For instance, Gaga AvBD1 refers to the AvBD1 gene from chicken (Gallus gallus). Duplicates within the same subfamily were numbered according to their deduced genomic positions.

Sequence alignment and gene phylogenetic construction

All deduced AvBD nucleotide sequences were aligned codon-by-codon based on exon structure using the MUSCLE v5 program (Edgar, 2004), with default settings and manual adjustments where necessary. Since the transversional (TVM) model + gamma distribution (G) + proportion of invariable sites (I), recommended by jModelTest 2 [38], is not supported by the Bayesian Evolutionary Analysis by Sampling Trees (BEAST) software package version 1.8.4 [39], this study employed the closest over-parameterized model, General Time Reversible (GTR) + I + G, for Markov chain Monte Carlo (MCMC) analysis [40]. The MCMC iteration ran for 100 million generations with a sampling frequency every 1,000 steps, using a random starting tree, lognormal Yule birth rate, and lognormal uncorrelated relaxed clock. After discarding the first 30% of generations as burn-in, all effective sample sizes exceeded 200. The maximum clade credibility tree was then generated and visualized using FigTree v.1.4.3.

Evolutionary selective pressure analysis

A codon-based maximum likelihood approach implemented in PAML 4.7 [41] was employed to assess positive selection at individual sites for each AvBD gene using site models. Likelihood ratio tests (LRTs) compared three model pairs: M1a (nearly neutral) vs. M2a (positive selection), M7 (neutral, β-distribution of ω < 1) vs. M8 (positive selection, β-distribution of ω > 1), and M8 vs. M8a (β-distribution of ω = 1). When LRT results indicated significance, Bayes empirical Bayes (BEB) was used to identify codons under selection in M2a and M8, with a posterior probability > 0.95. Concurrently, the mixed effects model of evolution (MEME), fixed effects likelihood (FEL), fast unconstrained Bayesian approximation (FUBAR), and single-likelihood ancestor counting (SLAC) methods, all implemented through the HyPhy via the Datamonkey platform [42] web interface (last accessed 26 December 2024), were used to detect indications of positively selected sites. Results from the FEL, FUBAR, and SLAC models were also employed to identify sites under negative selection. A significance threshold of p < 0.1 for MEME, FEL, and SLAC, and a posterior probability > 0.9 for FUBAR, was applied. Sites were considered positively selected when all tests in PAML (M1a vs. M2a, M7 vs. M8,and M8 vs. M8a) or HyPhy (MEME, FEL, FUBAR, and SLAC) indicated significance, while residues were classified as under negative selection when FEL, FUBAR, and SLAC collectively indicated significance.

Net charge and hydrophobicity estimation

For each AvBD gene, the mature peptide was deduced from the amino acid alignment of corresponding chicken counterpart. Net charge and hydrophobicity were estimated using online tools PepCalc.com and bioinformatics.org, respectively.

Ecological and life-history variables collection

In this study, nine continuous variables and three categorical variables (Additional file 3), potentially associated with pathogen exposure and fitness costs, were collected or calculated from previous records. The continuous variables included the Shannon–Wiener index of dietary diversity [43], sexual maturity age (mean value for males and females, days) [44, 45], longevity (years) [44, 46], weight (average body mass of males and females, grams) [44], clutch size [44, 47], incubation period (days) [44], evolutionary distinct and globally endangered score [48], breeding range cell count [48], and habitat elevation (average, meters) [45, 47, 49]. The categorical variables included breeding system (1: biparental care, 2: uniparental care) [47, 49], family system (1: no family living, 2: family living) [46, 47, 50], and movement type (1: resident, 2: migrant) [46, 47].

Association analyses between AvBD gene evolution as well as ecological factors and life-history traits

Phylogenetic generalized least squares (PGLS) [51], implemented in the R package ‘caper’, were used to assess the potential relationship between the evolutionary rate of each deduced single-copy AvBD subfamily (AvBD1-5 and AvBD8-14) and nine continuous factors, with log₁₀-transformed root-to-tip ω as the predictor variable. Briefly, 10,000 species phylogenetic trees for each AvBD gene dataset were downloaded from BirdTree [52], and the maximum clade credibility trees were constructed using TreeAnnotator in BEAST 1.8.4 [39] with a 25% burn-in. Evolutionary rates (root-to-tip ω) were estimated using the free-ratio model in CODEML 4.7.4 [41] by summing the dN values (dN_s) and dS values (dS_s) values from the root to the terminal species (tip), then calculating the dN_s/dS_s ratio as previously described [48, 53]. To ensure that root-to-tip ω values were appropriate for gene-trait associations, any cases where dN_s or dN_s were less than 0.0002 were excluded from further analysis. Both root-to-tip ω and the nine continuous variables were log₁₀-transformed for normalization. The maximum-likelihood method was used to estimate the lambda (λ) value, ranging from 0 to 1, indicating the presence of phylogenetic signal in the traits. Simultaneously, Markov Chain Monte Carlo generalized linear mixed models (MCMCglmms), in the MCMCglmm package [54], were employed to investigate the potential association between the root-to-tip ω of the 12 AvBD genes (AvBD1-5 and AvBD8-14) and the three categorical variables. Species phylogenetic relationships were incorporated as the covariance structure, with covariance V = 1 and degree of belief parameter (nu) = 0.002 as inverse-Wishart priors. MCMC chains were run for 140,000 iterations, with a 14,000 burn-in and sampling every 500 iterations. Autocorrelation was checked to ensure values were < 0.1. Variables with 95% credible intervals (CI) not overlapping 0 and pMCMC < 0.05 were considered statistically significant. Additionally, to assess whether AvBD repertoire evolution and species-specific duplication of AvBD6 and AvBD7 were related to the 13 variables, log₁₀-transformed total numbers of AvBD genes, as well as AvBD6 and AvBD7 gene counts in each species, were used in regression analyses following the same methods.

Three-dimensional structure visualization

The three-dimensional structures of mature AvBD2, AvBD7, and AvBD11 peptides were retrieved from the Protein Data Bank [55], with corresponding entry IDs 2LG5, 5LCS, and 6QEU, respectively. Tertiary structures of other AvBD sequences were predicted using SWISS-MODEL [56], employing templates with the highest identity scores. The six conserved cysteines and positions under selections were highlight by PyMOL v3.1.

Sample collection PCR amplification and sequencing

The experiment was conducted in compliance with the Regulations for the Administration of Affairs Concerning Experimental Animals (Ministry of Science and Technology, Beijing, China, revised in June 2004) and was approved by the Institutional Animal Care and Use Committee of China West Normal University, Sichuan, China (Approval No. CWNU2022D030). A total of 70 Galliformes samples were used to amplify the mature peptides of 12 deduced single-copy genes (AvBD1-5, and AvBD8-14). These included 26 Indian peafowls, 20 Chinese monals (Lophophorus lhuysii), 9 silver pheasants, 7 Reeves's pheasants (Syrmaticus reevesii), 4 ring-necked pheasants, 2 local dwarf chickens (Gallus gallus domestic), 1 golden pheasant (Chrysolophus pictus), and 1 Temminck's tragopan (Tragopan temminckii). Of the 20 Chinese monal samples, 17 were derived from eggshell membranes of captive populations, while 3 were muscle samples from wild individuals that died of natural causes, as previously described [57]. The remaining 50 samples were collected from eggshell membranes and specimens at the Changle Wild Animal Breeding Farm (Sichuan, China). Genomic DNA was extracted using DNAiso Reagent (TaKaRa, Beijing, China) according to the manufacturer’s protocol. DNA quality and concentration were verified before use in PCR reactions. Species identification for all samples was first confirmed by sequencing the Mitochondrially Encoded NADH Dehydrogenase 2 (MT-ND2) gene. Primers used in this study are listed in Additional file 4. PCR reactions were performed in 25 μl volumes, consisting of 12.5 μl of TSINGKE 2xT5 Super PCR Mix, 1 μl each of forward and reverse primers, 8.5 μl of sterilized distilled water, and 2 μl of DNA template. The PCR program was set as follows: 94 °C for 3 min, followed by 35 cycles of 94 °C for 30 s, annealing at the optimal temperature for 30 s, and extension at 72 °C for 30 s, with a final extension at 72 °C for 5 min. All amplified products were sequenced by Sangon Biotechnology Co., Ltd. (Shanghai, China). For target sequences longer than 800 bp, both forward and reverse sequencing was performed, and the sequences were assembled. The deduced mature peptides were predicted with reference to chicken sequences. All PCR products and functional AvBDs identified from the 25 genomes were aligned using the MUSCLE program [58].

Results

Identification of AvBD gene repertoire in Galliformes

Based on a homology search of the 14 AvBD mature peptides, a total of 364 sequence fragments containing the β-defensin motif were identified in the aforementioned database across 24 species of Galliformes. Among these, 360 sequences were identified in the genomes, while the remaining four fragments were exclusively found in the SRA database. Using a combination of predictive strategies, together with the known 14 chicken AvBD genes, 354 putative AvBD genes were identified, comprising 333 intact genes, 14 partial genes, and 7 pseudogenes (Fig. 1, Additional file 5 and Additional file 6). Of the 7 pseudogenes, the deduced sequences indicate that the pseudogenization of AlruAvBD3-pseudo and CojaAvBD6-pseudo is caused by mutations in the initiation codon, while PamuAvBD10-pseudo is pseudogenized due to a premature stop codon. The pseudogenization of CemiAvBD14-pseudo, CemiAvBD10-pseudo, CrmaAvBD7-pseudo, and LonyAvBD7-pseudo results from frame-shift mutations. The number of AvBDs varied among species, ranging from 12 in the Japanese quail and the Australian brush-turkey to 24 in the lesser prairie-chicken.

Fig. 1

Distribution of putative AvBD genes across 25 galliform genomes. Species are arranged based on phylogenetic relationships. Circle size represents the number of genes in the 14 AvBD lineages. Purple, blue, and gray circles indicate intact genes, partial genes, and pseudogenes, respectively

Sequences comparison of AvBD genes

In line with the typical defensin gene structure, the first two coding exons of each AvBD gene identified in this study are separated by a phase-1 intron. Notably, all three intron phases were observed between the second and third coding exons, following the strict GT/AG rule and known transcripts across the 25 species. Among the 11 genes (AvBD1-11) containing three coding exons, the intron phase between the last two coding exons of the AvBD8 sub-clade is phase-0, whereas all AvBD11 genes exhibit a phase-1 intron in corresponding regions. For the remaining nine AvBD genes (AvBD1-7, AvBD9, and AvBD10), a consistent phase-2 intron structure is present before the final coding exon (Additional file 6).

Multiple sequence alignment of AvBDs reveals that most coding regions are highly conserved, with the signal peptide region composed of approximately 20 amino acid residues and the mature peptide containing a β-defensin motif. The six cysteine residues in this motif are largely invariantly spaced, following the pattern C-X_2–4-G-X_1–2-C-X_3–5-C-X_2–8-G-X_1–6-C-X_5–6-CC (Additional file 5). Additionally, two glycine residues, positioned between the first and second cysteines and between the third and fourth cysteines, are conserved across the β-defensin motif, suggesting their importance for proper folding, conformation, and function. Interestingly, length variation was observed in the terminal sequences of certain AvBDs. For instance, chicken AvBD3 contains an additional 20 amino acids beyond the last conserved cysteine compared to other species.

Genomic organization of the AvBD gene cluster

In addition to the reported gene locations in chicken, golden pheasant, and Japanese quail, the genomic structures of 22 other Galliformes species were strategically analyzed. As illustrated in Fig. 2, all galliform genomes appear to encode a single cluster of AvBD genes, though some gaps remain in a few genomes. The order and transcriptional orientation of orthologous AvBDs are largely conserved across species, with the exception of turkey. Nine genes (AvBD1, −2, −5, and 8–13) exhibit one-to-one orthologs among the 25 species, although some subfamilies include pseudogenes or gene fragments from unassembled sequences. AvBD3, −6, −7, and −14 were absent in a few genomes, while AvBD6 and AvBD7 experienced significant expansion in multiple species, suggesting that these AvBD gene lineages have independently undergone birth-and-death processes. Despite two copies of AvBD4 being identified in the Chinese bamboo partridge, AvBD4 was treated as a single-copy gene in this study, as BathAvBD4B is located in a short fragment, and no gaps were found near BathAvBD4A in the reference genome.

Fig. 2

Genomic organization of AvBD gene clusters in 25 galliform species. The position of each gene is marked by the location of the β-defensin motif. Red genes are transcribed opposite to green ones, while pseudogenes are marked in blue. A broken line indicates a gap in the genomic sequence

Phylogenetic analysis of AvBDs

To investigate the evolutionary relationships among AvBDs, a phylogenetic tree was constructed using 333 intact AvBD genes from 25 avian species via the BEAST software package (Fig. 3 and Additional file 7). All AvBDs were grouped into 14 gene lineages with strong posterior probabilities of 1.0, except for the AvBD7 lineage, which had a posterior probability of 0.85. AvBD6, characterized by numerous species-specific paralogs, was grouped with AvBD7 with a posterior probability of 1.0, supporting the notion that AvBD6 duplicated from AvBD7 prior to the divergence of these species. Further classification of AvBD9/10/11/13/14 and AvBD1/3 into distinct clusters, with posterior probabilities of 0.90 and 0.91, respectively, suggests that each cluster originated from a common ancestral duplication. However, the phylogeny among other subgroups remains poorly supported, indicating that many of these genes have undergone diversifying selection during evolution.

Fig. 3

Phylogenetic relationships of AvBD genes across 25 species. Bayesian posterior probabilities for the 14 AvBD lineages are shown beside the nodes

Net charge and hydrophobicity estimation of the deduced mature peptides

The deduced mature peptide subgroups exhibit positive net charges, ranging from 1.28 in AvBD12 to 7.56 in AvBD1 (Fig. 4A), consistent with the theory that positively charged amino acid residues are essential for their biological functions. In terms of hydrophobicity, the deduced peptides displayed indexes ranging from −1.006 to 0.626 (Fig. 4B). While most lineages have negative or neutral hydrophobicity indexes, the AvBD3 subfamily stands out with an average index of 0.294, suggesting that this increased hydrophobicity may influence its antimicrobial mechanism, particularly in microbial membrane disruption.

Fig. 4

Comparison of net charge (A) and hydrophobicity index (B) of putative AvBD mature peptides. Dots represent average values, while the lines depict the ranges

Positive selection in different AvBD clusters

In the AvBD7, −8, −12, and −13 gene lineages, no significant positive selection was fully supported by the methods employed in either PAML or HyPhy (Additional file 8). The Bayes Empirical Bayes approach also indicated that AvBD5 contains only two positively selected sites, both located in or near the signal peptide region. In contrast, at least one positively selected site was identified in the mature peptide regions of the remaining nine gene lineages (Fig. 5 and Additional file 8). The number of positively selected sites varied significantly among AvBD subfamilies, with AvBD1 exhibiting the highest numbers (9 sites). Notably, site models identified positively selected sites in the C-terminal sequences encoded by the last exon in the AvBD1 and AvBD6 subfamilies. The results inferred from HyPhy also indicated that there are at least two sites under purifying selection in each of the 14 genes. In addition to the highly conserved six cysteines, a large number of positions were found to be invariant (identical at the nucleotide level), particularly in the AvBD13 and AvBD14 lineages (Fig. 5 and Additional file 9).

Fig. 5

Three-dimensional structure and location of positively selected sites in mature peptides from different AvBD lineages. Structures were modeled based on chicken amino acid sequences, except for AvBD3-Pamu. Positively selected sites are shown in red, purifying selected sites in green, and the six conserved cysteines are highlighted in yellow

Amino acid alleles of mature peptides

In order to briefly explore the distribution of mature peptide alleles across Galliformes, we sequenced 70 individuals across 8 galliform species, including the Indian peafowl, Chinese monal, silver pheasant, Reeves's pheasant, ring-necked pheasant, chicken, golden pheasant, and Temminck's tragopan (Additional file 10 and Additional file 11). Alongside the deduced peptides from the genomes of 25 species, a total of 28 Galliformes species were used to investigate variations in AvBD mature peptides across the order. Among the mature peptide sequences obtained, the number of amino acid alleles ranged from 7 in AvBD13 to 26 in AvBD1 (Fig. 6 and Additional file 12). Although certain alleles were shared across species at amino acid levels, widespread sharing was only observed in AvBD13, with the most common allele identified in 18 members of Phasianidae. On the other hand, amino acid substitutions were noted in a few cases, despite the majority of samples being obtained from captive populations (Fig. 6 and Additional file 12). The prevalence of species-specific amino acid alleles strongly suggests adaptive evolution driven by diverse ecological factors and life-history traits.

Fig. 6

Distribution of amino acid alleles for AvBD1 (A) and AvBD13 (B). Species are vertically aligned based on phylogenetic relationships from BirdTree. Shaded boxes indicate shared alleles among species. The number of shaded boxes from each branch tip reflects the count of alleles found in that species

Associations between AvBD gene evolution and related factors

To explore whether the prevalence of species-specific AvBD sequences is influenced by ecological and life-history variables, PGLS (Fig. 7 and Additional file 13) and MCMCglmm (Fig. 8 and Additional file 14) regressions were performed for continuous and categorical variables, respectively. Given that variations in the signal peptides or pro-segments can affect gene expression and peptide secretion and potentially lead to pseudogenization, the 12 deduced single-copy genes with intact sequences from the 25 genomes were included in the analysis. A significant negative correlation between log₁₀(root-to-tip ω) and incubation time was detected for AvBD2 (Fig. 7A), while a negative association was found between the evolutionary rates of AvBD10 and body weight (Fig. 7D). Additionally, longevity exhibited a significant relationship with the evolutionary rates of AvBD5 (Fig. 7B) and AvBD11 (Fig. 7E), although these two genes demonstrated opposite trends. Notably, the molecular evolution of AvBD9 (Fig. 7C), AvBD12 (Fig. 7F), and AvBD14 (Fig. 7G) was significantly correlated with habitat elevation, with AvBD9 showing a positive association and AvBD12 and AvBD14 displaying negative associations. The Japanese quail, the only migratory species in this study, exhibited significantly higher log₁₀(root-to-tip ω) values for AvBD2 and AvBD12, and lower values for AvBD4 compared to the other species, as revealed by MCMCglmm analysis (Fig. 8). These results suggest that the evolution of most AvBD genes has been shaped by habitat and/or life-history traits. However, none of the 13 variables were significantly related to the total number of AvBD genes, or to the number of AvBD6 or AvBD7 genes in our dataset (Additional file 13 and Additional file 14).

Fig. 7

Scatterplots showing significant relationships between AvBD evolutionary rates (root-to-tip dN/dS) and traits. Dots represent log₁₀-transformed root-to-tip dN/dS values, with lines representing phylogenetically controlled regression. R² and p-values from PGLS analysis are shown in the top right corner

Fig. 8

Boxplots illustrating significant relationships between AvBD evolutionary rates (root-to-tip dN/dS) and traits. The x-axis shows migrant or resident species, with dots representing log₁₀-transformed values. p-values from MCMCglmm analysis are displayed

Discussion

By searched through the sequences database of 25 species spanning major clades of Galliformes, 354 AvBD genes were identified in the present study. Sequences comparison indicate that C-terminal sequences are not fully conserved, with amino acid substitutions and length variations occurring since the divergence of these 25 species (Supplementary File 5). Previously, the extended C-terminal of chicken AvBD3 was attributed to a two-nucleotide insertion based on cDNA from chickens and turkeys [59]. However, cross-species genomic analysis revealed that this unusual structural feature of chicken AvBD3 (Fig. 5) most likely arose from a C-to-G mutation, creating a novel splicing site (AG) and causing a frameshift, leading to the long C-terminal extension (Supplementary File 15). These findings suggest that terminal untranslated sequences can serve as critical resources for defensin evolution, similar to observations in pig cathelicidins [60]. While it remains unclear which of the two mutations occurred first, further studies are needed to clarify the evolutionary timeline and functional divergence of these sequences. Additionally, a potential G-to-A substitution was found near the splice site in the coding region, indicating that if chicken AvBD3 was spliced at the common site observed in the other species (Additional file 15), the hypothetical peptide (Gaga-AvBD3h) would be two amino acids shorter than the typical form. Interestingly, the AvBD3 amino acid sequences in the California quail, the scaled quail, and the northern bobwhite are also two amino acids shorter than those in other species. However, these premature stop codons likely emerged independently, given the significant nucleotide-level differences.

We failed to identify AvBD3 gene in the Japanese quail, the Chinese bamboo partridge, and the golden pheasant. Defects at the AvBD3 locus have been observed in various lines of Japanese quail [20], while the phylogenetic relationship between the Japanese quail and the Chinese bamboo partridge further supports this observation. Additionally, our primers failed to amplify the defensin domain of AvBD3 in a golden pheasant sample, with only the third coding exon successfully amplified (data not shown). These findings suggest that AvBD3 may have been lost in certain Galliformes species. It has been proposed that the annotation of multiple genes situated within GC-rich regions poses significant challenges in avian genomes, primarily due to the inherent difficulties in sequencing and assembling these GC-rich segments [61]. Supporting this notion, our CG content analysis revealed that the coding sequences of AvBD genes exhibit a relatively high GC bias, with an average GC content of 53.4% (Additional file 16). Therefore, we cannot exclude the possibility that the absence of certain AvBD genes in our current analysis may be attributed to incompleteness in the available genome sequences. The evolutionary trajectories of the AvBD6 and AvBD7 gene lineages are notably complex. Although phylogenetic analysis suggests that AvBD6 and AvBD7 emerged before the divergence of the 25 species studied (Fig. 3), AvBD6 is exclusively found within the Phasianidae family (Fig. 1). Its absence in the Odontophoridae, Cracidae, and Numididae families suggests that AvBD6 was lost following the divergence of these families from Phasianidae. Given that duplications of AvBD6 [20] and AvBD7 [19] have been reported in galliform species, the possibility that some duplications result from copy number variations cannot be ruled out. Notably, AvBD6 duplications were observed in only eight genomes, all of which share a close phylogenetic relationship, with some AvBD6 paralogs suggested to have arisen much earlier than others (Additional file 7). This suggests that certain AvBD6 genes were fixed in the genome millions of years ago, while some species-specific paralogs may have originated from gene conversion followed by homologous recombination, similar to the evolutionary scenario observed in chicken AvBD6 and AvBD7 [62]. In Japanese quail, AvBD6 duplication is accompanied by pseudogenization of the locus [20], indicating that certain pseudogenes in reference genomes may result from pseudo-loci or null alleles.

Consistent with previous studies [16], the results indicate that AvBD genes are predominantly under negative selection, with episodic diversifying selection acts on a few specific sites. Considering the phylogenetic relationships of defensin genes from reptiles to birds, it is speculated that selection patterns in mature peptides may be linked to their homology with reptilian defensins. Among the seven AvBD genes that emerged before the divergence of birds and crocodiles (AvBD5, −8, −10, −12, −13, and −14) [26] or originated from ancient reptilian defensins (AvBD11) [63], a higher proportion of conserved sites is observed. In contrast, only a single positively selected site was identified in the mature peptides of AvBD11 and AvBD14, whereas multiple positively selected sites were identified in the bird-specific β-defensin cluster. Cross-genome analyses of 53 bird species, spanning Struthioniformes to Passeriformes, suggest that many AvBD3 paralogous genes originated at various points in avian evolution, with a significant number of sites under positive selection following gene duplication [16]. Aside from the conservative SLAC method, the other six model consistently indicated that AvBD3 contain at least two positively selected sites in the deduced mature peptide region in this study. Additionally, no duplication events were detected in the 25 bird genomes analyzed. These observations suggest that extensive diversification, rather than rapid gene duplication, has driven AvBD3 evolution in galliform species. Conversely, while duplications were observed in both AvBD6 and AvBD7, the two subgroups exhibited inconsistent evolutionary trends. Seven positively selected sites were identified in AvBD6, whereas no site in AvBD7 was fully supported by either PAML or HyPHy as being under positive selection. Surprisingly, two positions (site 45 and site 64) inferred to be under positive selection in AvBD6 showed evidence of purifying selection in AvBD7. These results support a neofunctionalization model, where the ancestral gene copy undergoes purifying selection to retain its original function, while the duplicate evolves under positive selection to acquire a new function [64]. Nevertheless, the species-specific evolution of AvBD6 and AvBD7 remains difficult to fully explain. Sequence changes are a hallmark of adaptive evolution, driving functional diversification. However, the precise reasons why selection has favored the evolution of certain distinctive and diversified defensin repertoires remain elusive. Further experimental studies are needed to explore functional differences among the various positively selected sites, to better understand AvBD diversity.

Unlike reports in Anseriformes [24], no widely shared amino acid alleles were identified in AvBD4, −5, and −10 across galliform species, which may explain the significant differences in selective pressure between Galliformes and Anseriformes for these genes. In contrast, identical amino acid alleles of AvBD13 were widely observed across both Galliformes and Anseriformes, likely due to its evolutionary conservation and lack of positively selected sites in both Galliformes and Anseriformes. These results further underscore the gene-, species-, and timescale-specific evolution of β-defensin genes. Widespread alleles across multiple species likely reflect evolutionary advantages tailored to the biological functions of closely related species. On the other hand, only a few nonsynonymous mutations were detected within species, most likely due to the limited number of individuals and species available for genotyping, with the majority of samples originating from captive populations.

Pathogen-mediated selection can have contrasting effects on the immune system. One view, posits that both the innate and adaptive immune systems are highly plastic and continuously evolve to counter rapidly changing pathogen fauna [28, 65]. In contrast, another view suggests that highly efficient gene types have already been selected through natural selection, and effective alleles exhibit low tolerance for variants [66], though constraints may relax under low infection pressure [29]. Furthermore, the trade-offs hypothesis proposes that immune responses are energetically costly and can contribute to autoimmune processes, imposing fitness costs or becoming toxic when pathogen pressure diminishes [29, 31]. Comparative analyses revealed that the evolution of nine AvBD genes (AvBD2, −4, −5, −8, −9, −10, −11, −12, and −14) is significantly correlated with five of the 13 ecological factors and life-history traits examined. However, no uniform pattern of response emerged, as certain factors exert gene-specific effects on AvBDs. AvBD2, one of the earliest expressed host defense peptides, is found in the chicken yolk sac as early as embryonic day 7 [67]. Since the embryo lacks adaptive immunity, host defense peptides serve as a critical first line of defense by enhancing innate immunity [68]. The negative correlation between AvBD2's evolutionary rate and incubation time may be explained by prolonged pathogen exposure during extended incubation, which applies stronger selection pressure on the innate immune system. Natural selection has optimized AvBD2's antimicrobial activity to combat pathogens prevalent during incubation. However, highly efficient variants may also incur fitness costs when infection pressure is low. As a result, species experiencing lower pathogen pressure during incubation may tolerate higher mutation rates, allowing variants with compromised antimicrobial activity but beneficial physiological functions to become fixed, reflecting a trade-off between innate immunity and other physiological processes. Birds in high-altitude environments face increased basal metabolic rates and limited food availability, creating ongoing energy challenges [69]. Additionally, higher altitudes typically reduce pathogen richness and disease transmission due to lower temperatures and increased ultraviolet radiation [70, 71]. Therefore, the positive correlation between AvBD9 evolution and habitat elevation similarly suggests that highly efficient AvBD9 alleles may impose fitness costs in sterile environments. Given that long lifespan and migration increase pathogen exposure, AvBD11 and AvBD4 are thought to primarily function in infection clearance, aligning with the evolutionary patterns observed in efferent molecules of the innate immune system. Migratory pressures may also impose strong purifying selection on AvBD4 in migratory waterfowl [24], contributing to the distinct differences in positively selected sites between Galliformes and Anseriformes. Moreover, larger animals are increasingly recognized as facing greater pathogen exposure and may adaptively evolve diverse strategies to boost immunity [72, 73]. This suggests that AvBD10 may have evolved to maintain robust pathogen tolerance, particularly in larger galliform species.

The other significant associations identified between AvBD evolution and ecological and life-history characteristics—such as the positive correlation between AvBD5 and longevity, the negative correlation with habitat elevation for AvBD12 and AvBD14, and the higher evolutionary rates of AvBD2 and AvBD12 in migratory species—are likely driven by factors beyond direct pathogen elimination [74, 75]. Many AvBDs, such as AvBD2 and AvBD5, are expressed in various immune cell types and have the capacity to modulate immune signaling and antigen presentation [76]. Additionally, host defense peptides have been shown to induce tumor cell death, potentially extending lifespan by exerting anti-cancer effects during aging. In such cases, AvBDs may function as afferent molecules, with genetic diversity conferring advantages by allowing hosts to recognize and target a wide range of receptors and aberrant host cells. Furthermore, host defense peptides are integral to maintaining a balanced microbiome, and many can be induced by metabolites produced by symbiotic bacteria [77, 78]. Certain amino acid alleles may help prevent microbiome dysbiosis, particularly in migratory species where high energy consumption could otherwise lead to imbalance. Alternatively, birds living at high altitudes, where basal metabolic rates are significantly elevated [69], may reduce the evolutionary rates of specific genes to tradeoff between immune investment and other life-history traits. Although the antimicrobial activity of a single allele can be broad, it is important to consider that amino acid substitutions often lead to shifts in microbial specificity [79]. Specific genes may need to retain greater variability to keep pace with the co-evolutionary arms race between host and pathogen, as well as to explore a wider functional space to combat emerging diseases.

Conclusion

In conclusion, the identification and analysis of AvBD genes within the primary clade of living galliform species reveal widespread duplication, pseudogenization, gene loss, and sequence variation, contributing to substantial diversity in both gene sequences and numbers across species. The evolution of AvBDs exhibits varied responses to distinct evolutionary pressures and life-history traits, underscoring unique evolutionary relationships with infection risk and emphasizing the multifunctional roles of β-defensin genes in birds. This study provides valuable insights into the adaptive immune responses of hosts to diverse ecological and life-history factors, laying the foundation for further exploration of their evolutionary trajectories before considering them as potential novel antimicrobial agents against drug-resistant pathogens. Additionally, further mechanistic and evolutionary research is needed to fully understand why ecological factors and life-history traits differentially influence the evolution of diversified defensin groups.