Molecular evolution in immune genes across the avian tree of life

Abstract

All organisms encounter pathogens, and birds are especially susceptible to infection by malaria parasites and other haemosporidians. It is important to understand how immune genes, primarily innate immune genes which are the first line of host defense, have evolved across birds, a highly diverse group of tetrapods. Here, we find that innate immune genes are highly conserved across the avian tree of life and that although most show evidence of positive or diversifying selection within specific lineages or clades, the number of sites is often proportionally low in this broader context of putative constraint. Rather, evidence shows a much higher level of negative or purifying selection in these innate immune genes – rather than adaptive immune genes – which is consistent with birds' long coevolutionary history with pathogens and the need to maintain a rapid response to infection. We further explored avian responses to haemosporidians by comparing differential gene expression in wild birds (1) uninfected with haemosporidians, (2) infected with Plasmodium and (3) infected with Haemoproteus (Parahaemoproteus). We found patterns of significant differential expression with some genes unique to infection with each genus and a few shared between ‘treatment’ groups, but none that overlapped with the genes included in the phylogenetic study.

Introduction

Our understanding of the evolution of the avian genome has recently been updated by the landmark avian phylogenomic paper by Jarvis et al. (2014). This phylogeny includes all extant orders from Neoaves (e.g. songbirds, parrots, pigeons) and species from several clades in the more basal orders including Anseriformes (waterfowl), Galliformes (landfowl), Tinamiformes (tinamous) and Struthioniformes (ostrich). The phylogeny supports the rapid radiation of the most derived order Passeriformes (i.e. perching birds) and their sister taxon, the parrots (Psittaciformes). This diversification occurred around 60–55 million years ago during the Late Paleocene in Gondwana (Gill, 1995; Jarvis et al., 2014). Passeriformes account for more than half of all avian species, are found on every continent except for Antarctica, and occupy a diverse range of habitats (Gill, 1995). Although birds as a whole appear to have evolved at a conservative molecular evolutionary rate when contrasted with other large groups of organisms (notably mammals), passeriform songbirds have an average mutation rate that is close to twice that of other avian species (Jarvis et al., 2014; Zhang et al., 2014).

Immune genes are shaped by selection pressure from an onslaught of coevolving pathogens (Hendrick, 1998; Acevedo-Whitehouse and Cunningham, 2006; Piertney and Oliver, 2006). The genes of the major histocompatibility complex (MHC) have been the most useful for understanding the adaptive immune response in wildlife, particularly in birds (Acevedo-Whitehouse and Cunningham, 2006; Turner et al., 2012), and many other systems (natural and model or experimental) have been used to explore gene-for-gene coevolutionary processes (reviewed in Brockhurst and Koskella, 2013). However, the role of selection in shaping diversity in other aspects of the immune system, such as innate immunity, remains poorly understood (Acevedo-Whitehouse and Cunningham, 2006; Vinkler and Albrecht, 2009; Grueber et al., 2012, 2014).

To date, the only non-MHC immune genes to be extensively studied in birds are Toll-like receptors (TLR). TLRs are membrane-bound sensors of the innate immune system that recognize invariant and distinctive molecular features of invading microbes and are essential for initiating adaptive and innate immune responses in vertebrates. These types of receptors are part of a large family of proteins called pattern-recognition receptors (PRRs). The genetic variation at TLR genes has been directly related to differential pathogen responses in birds, humans and livestock (Alcaide and Edwards, 2011; Grueber et al., 2012; Grueber et al., 2014; Bateson, et al., 2016; Raven et al., 2017; Kannaki et al., 2018; Nelson-Flower et al., 2018). Prior research has shown that vertebrate TLR evolution is mostly characterized by purifying (negative) selection (Barreiro et al., 2009; Alcaide and Edwards, 2011; Nelson-Flower et al., 2018) and, to a lesser extent, balancing selection (Ferrer-Admetlla et al 2008). However, Grueber et al. (2014) also found evidence of episodic positive selection in avian TLRs. This discovery was aided by the doubling of avian transcriptomic and genomic data made available at the time of their research (Grueber et al., 2014). These new data allowed Grueber et al. (2014) to evaluate up to 23 bird species per gene and gave added resolution in discerning patterns of selection.

More generally, PRRs are highly conserved and found in a diverse range of animal species. They provide the foundation upon which the hosts' immune response to pathogens has evolved (e.g. Bagheri and Zahmatkesh, 2018). PRRs account for many other families of sensors beside TLRs (Acevedo-Whitehouse and Cunningham, 2006). These include soluble components (e.g. lipopolysaccharide-binding protein, collectins and pentraxins) and cell-associated components (e.g. TLRs, C-type lectins, NOD-like receptors, RIG-I like receptors, scavenger receptors, formyl-peptide receptors and various intracellular receptors; Chen et al., 2013). Pattern recognition is considered to occur at three different levels. These are interactions involving soluble extracellular PRRs, membrane-bound PRRs and PRRs found in the cytoplasm (see also Zhang et al., 2018). Most major pathogen-response groups are included in one or more of these three levels. The redundancy in recognition of particular pathogen-associated molecular patterns and coordination of different PRRs to sense any single invading class of microorganism is an important feature of the PRR network. The outcome of PRR-mediated responses is to initiate downstream immune cascades that represent the first stage in establishing the immune response. These cascades are in turn mediated by a large group of regulatory proteins called cytokines. Cytokines include interleukins, interferons, tumour necrosis factors and chemokines (Chen et al., 2013). Variation in cytokine genes has been associated with certain pathogens in humans, but relatively little research is available for birds, and what is available is primarily limited to fowl (e.g. Chhabra et al., 2015; Kannaki et al., 2018).

PRR immune genes have shown signatures of selection that can be associated with specific bird species or evolutionary families throughout the avian phylogenetic tree (Grueber et al., 2014). For example, chickens, unlike ducks, geese and finches, lack RIG-I, the PRR sensor used to detect highly pathogenic avian influenza. Instead, they express a gene called the melanoma differentiation-associated gene 5 (MDA5), which compensates for the missing RIG-I (Chen et al., 2013; Magor et al., 2013). Some PRR-induced interferon genes are missing from birds completely. These include interferon genes ISG15, ISG54 and ISG56 (see Velova et al., 2018). Birds as a whole have a reduced repertoire of immune genes when compared to mammals (Magor et al., 2013), and because of this, we could predict that avian PRRs are evolutionarily constrained (Chen et al., 2013; Magor et al., 2013; Grueber et al., 2014).

To explore PRR evolution in birds, we conducted an in silico analysis of bird gene sequences in the GenBank. Additionally, we examined gene expression differences in Northern cardinals uninfected by haemosporidians, infected by Plasmodium haemosporidians, or infected by Haemoproteus haemosporidians. The computational analyses suggest that some avian PRRs are under positive, diversifying selection, while the differential gene expression analysis highlighted the differences in response among naturally-occurring ‘treatment’ groups.

Methods

Sequence mining and alignment

Using annotations of immune genes from the GenBank, we used BLAST against the known genomes and transcriptomes of birds and aligned 47 individual gene datasets (i.e. fasta files). We included most PRRs that we could align with at least 250 base pairs (bp), some non-PRR genes and some from the adaptive immune system with at least 250 bp; we did not include all the TLRs because they had been analysed in another study (Grueber et al., 2014). We then trimmed each down to contiguous alignments for each taxon that were aligned at a first codon position and ended at a third codon position using webPrank (https://www.ebi.ac.uk/goldman-srv/webprank/). There were no ambiguities, stop codons or missing data in any alignment, and all genes were checked for recombination using RDP4 (Martin et al., 2015). Numbers of taxa per dataset ranged from 18 to 57, and gene regions ranged in length from 279 to 2346 bp. Only coding regions were used in these analyses. All gene alignments and MEME results are freely available in the Dryad Depository (doi:10.5061/dryad.r5594f0).

Characterization of gene evolution

We analysed each dataset to determine patterns of selection using several methods within the HyPhy platform (via the Adaptive Evolution server at DataMonkey.org; Kosakovsky Pond and Frost, 2005; Delport et al., 2010). The single-likelihood ancestor counting (SLAC; Kosakovsky Pond and Frost, 2005) method was used to determine the overall d _N/d _S per gene dataset. The fixed-effects likelihood (FEL; Kosakovsky Pond and Frost, 2005) method was used to determine the number of sites (codons) under negative selection; FEL fits a site-specific likelihood model on a per-codon basis and then determines whether d _N is significantly greater than d _S. FEL reconstructs the phylogeny using maximum likelihood while optimizing parameters. The mixed effects model of evolution (MEME; Murrell et al., 2012) was then used to determine the number of positively-selected sites; MEME essentially fits a likelihood model assuming a background d _N/d _S and determines, on a site-by-site basis, whether d _N is greater than the background rate. Like FEL, MEME reconstructs the phylogeny using maximum likelihood while optimizing parameters. Finally, even with the limitations recently pointed out by Venkat et al. (2018), we used adaptive branch-site REL (aBSREL; Kosakovsky Pond et al., 2011; Smith et al., 2011), a branch-based rather than site-based model, to determine if and where episodic diversifying selection has occurred in branches within trees reconstructed with distance methods (per the programme).

Comparative transcriptome analyses

In a previous study (Walstrom and Outlaw, 2016), we collected blood samples from Northern cardinals (Cardinalis cardinalis). All samples were adult (after hatch-year) females. We determined whether each animal was (1) uninfected with haemosporidians, (2) infected with Plasmodium or (3) infected with Haemoproteus (Parahaemoproteus) parasites. As only a small amount of blood was required to determine the haemosporidian infection status, excess blood for each cardinal was stored at −80°C. In the current study, we extracted RNA from blood samples of three birds uninfected with haemosporidians, three birds infected with Plasmodium (two with OZ03OZ01 and one with SIAMEX01) and three birds infected with Haemoproteus (two with OZ45_MEX19 and one with TUTI233) for a total of nine samples. RNA was extracted using a chloroform/ethanol protocol followed by a DNase digestion, and then cleaned using a Qiagen RNeasy Mini kit (Qiagen.com). cDNA libraries were prepared using a NEBNext^® Ultra™ II RNA Library Prep Kit for Illumina and run on a HiSeq (Illumina.com). Salmon (v0.9.1, Patro et al., 2017) generated transcript expression estimates for each library by mapping the raw Illumina reads to the zebra finch (Taeniopygia guttata) transcriptome (v3.2.4, Warren et al., 2010). We recovered ~11 200 genes from each ‘treatment’. Gene expression estimates were produced by combining the transcript expression estimates from the same gene for each sample (Sonenson et al., 2015). Genes expressed in fewer than three of the nine samples were removed from further analysis (Bourgon et al., 2010). The gene expression estimates were normalized based on sequencing depth of each library. EdgeR (v3.12.1; Robinson et al., 2009) was used to find differentially expressed genes. Generalized linear models were used to test for differential expression (DE) based on contrasting each treatment against each other, i.e. Plasmodium-infected vs uninfected, Haemoproteus-infected vs uninfected and Plasmodium-infected vs Haemoproteus-infected. Genes with an FDR adjusted P value ≤ 0.05 were considered differentially expressed. DE results are freely available in the Dryad Depository (doi:10.5061/dryad.r5594f0).

Results

Gene evolution

The number of taxa in each aligned immune-associated gene dataset ranged from 18 to 57 with an average number of 39 species. The number of bps included in each dataset ranged from 279 to 2346 with an average of 940 bps. The number of positively selected sites ranged from 0 to 29 (detected using MEME) and the number of negatively selected sites (detected using FEL) ranged from 0 to 398. Across all datasets d _N/d _S ranged from 0.040 to 0.789. Eight phylogenetic trees generated for specific gene datasets showed statistically significant evidence of positive, diversifying selection and 39 showed no evidence of positive, diversifying selection (detected using aBSREL). These results are listed in Table 1. Two trees with evidence of positive selection – one with high taxon and high bp values and one with low taxon and low bp values, and two trees with no evidence of positive selection – one with high taxon and high bp values and one with low taxon and low bp values are shown in Fig. 1 (from aBSREL analyses). Note that MEME and aBSREL analyses returned similar topologies.

Fig. 1. Gene trees reconstructed with aBSREL. Branch lengths are proportional to the number of codon substitutions. Songbirds (Passeriformes) are the purple clade and game birds (Anseriformes and Galliformes) are the green clade.

Table 1. Summary of results for each gene from MEME (# + sites), FEL (# − sites), SLAC (d _N/d _S) and aBSREL (diversifying selection) analyses in HyPhy

^a DEXH-box helicase 58 is involved in nucleic acid binding and helicase activity (https://www.genecards.org/cgi-bin/carddisp.pl?gene=DHX58&keywords=DEXH-box).

^b NFIL3 is involved in DNA-binding transcription factor activity and transcription corepressor activity (https://www.genecards.org/cgi-bin/carddisp.pl?gene=NFIL3&keywords=NFIL3).

^c NOD-like receptor P3 is involved in peptidoglycan binding (https://www.genecards.org/cgi-bin/carddisp.pl?gene=NLRP3&keywords=NOD-like,P3).

^d TLR7 is involved in transmembrane signalling receptor activity and double-stranded RNA binding (https://www.genecards.org/cgi-bin/carddisp.pl?gene=TLR7&keywords=TLR7).

Trees from genes in bold are shown in Fig. 1.

Differential gene expression

All samples in this study were collected from free-living birds, and although our polymerase chain reaction-based tests were repeated several times and have been shown in our lab and in others to be highly reliable, we did not control or specifically account for infection with other pathogens. Relevant values for significantly differentially expressed genes are listed in Table 2. Six genes were uniquely downregulated and four were uniquely upregulated in samples classified as being infected solely with Haemoproteus parasites. Nine genes were uniquely downregulated and seven genes were uniquely upregulated in samples putatively infected solely with Plasmodium parasites. Three genes were upregulated in both Haemoproteus-infected samples and Plasmodium-infected samples. The first of these, SIX Homeobox 4 is found to be involved in eye development in Drosophila and is a transcription factor in mouse development (NCBI gene). The second, immunoglobulin superfamily2C member 11, is involved in cell adhesion (NCBI gene), and the third, tetratricopeptide repeat domain 37, is involved in protein–protein interactions (NCBI gene). The small sample sizes almost certainly bias our results and affect the substantiation of our interpretation, but regardless of these limitations, DE is evident between groups (see also Supplementary Table 1). Surprisingly, none of the genes included in the phylogenetic analyses was found in the DE analyses. Neither did our DE analyses include overlap with the results from Videvall et al. (2015), in which the authors examined patterns of DE in response to experimental infection of birds with Plasmodium.

Table 2. Comparative patterns of gene expression between uninfected birds and those infected with Plasmodium or Parahaemoproteus parasites

Upregulated genes are in light grey and downregulated genes are in darker grey. See also Supplementary Table 1.

Discussion

Recently, the genomic sequences of 45 species have been made available from across the avian tree of life (Zhang, 2014). These data should provide a more complete picture of immune gene evolution in birds, including both TLRs and other non-MHC-related immune genes. To date, the only examination of molecular evolution in avian immune genes has been researched by making pairwise comparisons of orthologous genes between bird species (see Ekblom and Galindo, 2010; Ekblom et al., 2011, but see Minias et al., 2018). Here we evaluated signatures of natural selection in avian PRR gene families using a phylogenetic framework.

Like humans, birds are highly susceptible to malaria (haemosporidian) parasite infection and are confronted with no fewer than four genera of malaria parasites (Valkiunas, 2004; see also Bichet et al., 2015). Because malaria has had such a profound effect on the human genome (see Kwiatkowski, 2005) it would seem necessary that we look for the effects of molecular evolution in avian genomes. However, given that it is relatively unknown which genes are turned on or off in response to infection with malaria parasites, studies – like this one – must first compare birds with different types of parasite infections. Many genetic mechanisms of malaria resistance in humans were shaped by evolution in such a way that they have even become harmful to the constituent populations. The most well-known example is heterozygote advantage in the HBB (haemoglobin beta) gene, which confers such resistance to human Plasmodium falciparum malaria that the deleterious allele remains at 10% frequency or higher even though the homozygous recessive is lethal (Kwiatkowski, 2005). Such costly mechanisms have arisen independently many times and in response to different species of human Plasmodium (summarized in Kwiatkowski, 2005).

General evolutionary trends

Innate immune genes of birds appear to evolve in a highly conserved manner under strong purifying selection (see Table 1) which is entirely consistent with their evolutionary role. In the eight phylogenetic trees with at least one branch under positive selection, selection occurred mostly on terminal branches and there was no identifiable pattern with which taxa were under selection.

Haemosporidian malaria parasites elicited different transcriptomic responses from our wild-caught Northern cardinal samples. None of the significantly differentially expressed genes overlap with the immune genes in the first part of the study. Of note, a few of the differentially expressed genes in the birds are known to be involved in immune responses and signal transduction. With further sampling of other bird species, differentially expressed genes may be candidates for further studies of avian–haemosporidian interactions. However, considerable caution is warranted as the birds were caught in the wild and not tested for infection with other parasites, and the sample size is very small. Thus, the observed DE may reflect the effects of many other environmental factors.

Adaptive strategy

The adaptive strategy of avian innate immune genes is their conservation: a d _N/d _S less than 1. Overwhelmingly, genes of the innate immune system have been under intense selection to stay the same, although it is likewise clear than many sites are under positive selection even against the background of constraint. Therefore, why then, have other studies found immune genes to be rapidly evolving? The most comprehensive study was conducted between a galliform (chicken) and a passerine (zebra finch), two branches of the avian tree of life separated by ~40% genomic sequence divergence (Ekblom et al., 2011). In a follow-up study, Ekblom et al. (2011) summarized across gene ontology groups within the zebra finch genome and found d _N/d _S values about half of what we see here. On the other hand, analyses of selection in the MHC genes show rampant evidence of positive selection, a d _N/d _S greater than 1 within passerines (Minias et al., 2018), which is much higher than we found. More recently, big data from across the avian tree of life have shown that songbirds radiated in a big bang – even with low overall values of d _N/d _S (Zhang et al., 2014) – as we have seen in our gene trees (Fig. 1a–d). This big bang would at least partially explain the elevated rates that Ekblom et al. (2011) found between these very divergent lineages, except that their overall estimate of selection was closer to neutral, with positive selection showing up at specific residues rather than across whole genes.

Our analyses paint a more complex picture given that we could leverage more data than in previous studies. Immune genes are a blend of both positively and negatively selected sites in a constrained background, and their complex patterns of selection reflect their complex roles. The complexity of evolutionary interactions between hosts and pathogens is becoming increasingly apparent, particularly within the plant literature (e.g. Thrall et al., 2016 and references therein), and mathematical models to understand these interactions increasingly refer to them as ‘tangled’. Detailed analyses of sites under selection and the interpretation thereof are beyond our expertise, but the myriad parasites that differentially affect avian clades suggest that this pursuit would unlock many secrets about the evolution of avian immune genes.

Although studies of parasites in wild birds are extremely patchy in geographic location and host taxonomy, there is ample evidence to suggest that some parasites infect all birds but that different birds are affected by different parasites in many cases. Plasmodium (Atkinson, 2008), Toxoplasma (Dubey, 2008) and capillarid nematodes (Yabsley, 2009) are cosmopolitan and seemingly ubiquitous across birds (and other vertebrate taxa). On the other hand, Trichomonas primarily infects doves, raptors and parrots (Forrester and Foster, 2008), and the Dispharynx/Echinuria/Streptocara complex and tracheal worms (Fernando and Barta, 2008) are most commonly found in songbirds and galliformes. There are many other examples of host-specificity (at the order and family level) in bird parasites (see Wobeser, 2008; see also Kumar et al., 2017) and coevolutionary relationships between these hosts and parasites presumably over millions of years have shaped the adaptive responses within and between avian clades. Differences between the innate and adaptive immune systems are profound, and therefore, our expectations of patterns of selection between these groups of genes are likewise very different. Here, we focused on (primarily) genes of the innate immune system, and found an overall signature of purifying selection. However, where there are clear distinctions between groups of birds that are infected by different parasite groups, there is evidence of molecular evolution that may correspond with these differences.