Reduction in population size and not a shift in parasite community affect evolution of immune genes in island birds

1 Shared ecological conditions encountered by species that colonize islands often lead to the 2 evolution of convergent phenotypes, commonly referred to as “island syndrome”. Reduced 3 immune functions have been previously proposed to be part of the island syndrome, as a 4 consequence of the reduced diversity of pathogens on island ecosystems. According to this 5 hypothesis, immune genes are expected to exhibit genomic signatures of relaxed selection 6 pressure in island species. In this study, we used comparative genomic methods to study 7 immune genes in island species (N = 20) and their mainland relatives (N = 14). We gathered 8 public data as well as generated new data on innate (Toll-Like Receptors, Beta Defensins) 9 and acquired immune genes (Major Histocompatibility Complex classes I and II), but also on 10 hundreds of genes annotated as involved in various immune functions. As a control, we used 11 a set of 97 genes not involved in immune functions, to account for the lower effective 12 population sizes in island species. We used synonymous and non-synonymous variations to 13 estimate the selection pressure acting on immune genes. For the genes evolving under 14 balancing selection, we used simulation to estimate the impact of population size variation. 15 We found a significant effect of drift on immune genes of island species leading to a reduction 16 in genetic diversity and efficacy of selection. However, the intensity of relaxed selection was 17 not significantly different from control genes, except for MHC class II genes. These genes 18 exhibit a significantly higher level of non-synonymous loss of polymorphism than expected 19 assuming only drift and an evolution under frequency dependent selection, possibly due to a 20 reduction of extracellular parasite communities on islands. Overall, our results showed that 21 demographic effects lead to a decrease in the immune functions of island species, but the 22 relaxed selection caused by a reduced parasite pressure may only occur in some immune 23 genes categories. 24

Introduction majority of studies that focused on island avifauna have found ambiguous results, with either 66 no support for a reduced immune response on island species (Matson 2006;Beadell et al. 67 2007), or contrasted results, such as a lower humoral component (total immunoglobulins) on 68 islands, but a similar innate component (haptoglobin levels) between island and mainland 69 species ). The use of immune parameters as proxies of immune function is 70 fraught with difficulties ). The study of molecular evolution of immune genes 71 therefore represents an alternative strategy to tackle this question. However, it is necessary 72 to distinguish neutral effects, the demographic effects resulting from island colonization, from 73 selective ones, the potential relaxation of selection pressures due to the changes in the 74 pathogen community. 75 First, the bottleneck experienced by species during island colonization leads to a decrease in 76 genetic variability (Frankham 1997). A reduced genetic diversity at loci involved in immunity 77 should have a direct implication on immune functions (Hawley et al., 2005; Hale & Briskie, 78 2007 but see Spurgin et al. 2011). Second, small population sizes increase genetic drift, which 79 may counteract the effect of natural selection of low-effect mutations, especially weakly 80 deleterious mutations (Ohta 1992). Several recent studies found a greater load of deleterious 81 mutations in island species (Loire et  Toll-Like Receptors (TLR; transmembrane proteins) trigger a chain reaction leading to the 87 production of various substances, including antimicrobial peptides such as beta-defensins 88 (BD) that have active properties in pathogen cell lysis (Velová et al. 2018). On the other hand, 89 the acquired immune system allows a specific response, characterized by immune memory. 90 Major Histocompatibility Complex (MHC) genes code for surface glycoproteins that bind to 91 antigenic peptides, and present them to the cells of the immune system; class I and II genes 92 ensure the presentation of a broad spectrum of intra-and extracellular-derived peptides, 93 respectively (Klein 1986  island birds compared to continental species exhibit a genome-wide reduction in genetic 136 diversity and efficacy of selection (Kutschera et al. 2020; Leroy et al. 2021). Therefore, we 137 expect a similar reduction in immune genes diversity even without any change in the parasite 138 pressure. 139 To disentangle the effect of population size from a change in parasite pressure and estimate 140 the impact of demography on the efficacy of selection, we randomly selected protein-coding 141 genes (i.e., control genes) implied in various biological functions (Fijarczyk et al. 2016;Leroy 142 et al. 2021). If a reduced parasite pressure on islands directly impacts the evolution of 143 immunity genes, the genetic diversity of immunity genes is expected to show a larger variation 144 between island and continental species than the control genes. More specifically, for genes 145 under purifying selection, non-synonymous weekly deleterious mutations, normally eliminated 146 under strong selection, would be maintained, leading to an increase of genetic diversity. By   Alignments of Coding DNA Sequences (CDS) of individuals from 24 species were obtained 183 from Leroy et al., (2021). In addition, data for ten other species (six and four from islands and 184 mainland, respectively) were newly generated for this study by targeted-capture sequencing. 185 Blood samples and subsequent DNA extractions were performed by different research teams. 186 The complete dataset consisted of 34 bird species (20 and 14 insular and mainland species 187 respectively; Table 1; Figure 1). We filtered alignments in order to retain only files containing 188 a minimum of five diploid individuals per site (Table 1)  198 . CC-BY-NC 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in  This allowed us to identify Phylloscopus trochilus as an outlier. Unlike for all other species 283 (e.g. Fringilla coelebs, Figure S3), synonymous polymorphism level was very dependent on 284 the number missing data tolerated in P. trochilus alignments ( Figure S3). As such, we 285 excluded P. trochilus from further analysis. 286 The mean Pn/Ps, calculated from the concatenated sequences of genes from the same gene 287 class (control genes; BD; TLR; MHC I; MHC II; Database-group; Sma3s-group), was 288 estimated for each bird species. Alternative transcripts were identified based on the genomic 289 position in the GFF file. If several transcripts were available, one transcript was randomly 290 selected. Pn/Ps estimates based on less than four polymorphic sites were excluded from the 291 analysis, as were those with no polymorphic non-synonymous sites. 292

Statistical analyses 293
To estimate the impact of demographic history on genome-wide polymorphism of island 294 species and the potentially reduced constraints on their immunity genes, we computed the 295 ratio of non-synonymous genetic diversity over synonymous genetic diversity (Pn/Ps). We 296 distinguished the part due to the island or mainland origin of species from the one due to the 297 gene category (i.e., immunity versus control genes). A linear mixed model was performed, 298 using the Pn/Ps ratio as dependant variable and, as explanatory variables, the mainland or 299 . CC-BY-NC 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in  Kuznetsova et al. 2017)). In order to take into account the phylogenetic effect, the 301 taxonomic rank "family" was included as a random effect in the model. Five linear mixed 302 models were defined i) null model, ii) model with only the origin parameter, iii) model with only 303 the gene category parameter, iv) model using both origin and gene category parameters, and 304 finally v) model including those two parameters and the interaction effect. In some cases, the 305 phylogenetic effect was difficult to estimate because the number of species per family was 306 reduced to one. In that case, we chose to reduce the number of families by grouping Turdidae 307 with Muscicapidae, Nectariniidae, and Estrildidae with Ploceidae and Fringillidae within 308 Thraupidae. The results obtained with these family groupings were similar to the original model 309 (Table S1), except when stated. The categories Database-group and Sma3s-group were 310 tested separately from the Core group because they contained hundreds of genes annotated 311 using the automatic pipeline that were only available for species with genome wide data. 312 Database-group and Sma3s-group were not analysed simultaneously because they contained 313 a partially overlapping set of genes. Finally, genes evolving under purifying selection and 314 genes evolving under balancing selection were also analysed separately. Model selection was 315 based on two methods. First, we use the difference in corrected Akaike Information Criterion 316 (ΔAICc) calculated using the qpcR package (Spiess and Spiess 2018). Second, a model 317 simplification using an ANOVA between models was also performed. 318 We also tested an alternative model using the difference between Pn/Ps of immunity genes 319 and control genes ( Pn/Ps) as dependent variable, and species origin as explanatory variable. 320 Under the hypothesis of a relaxation in selection pressure on islands due to a change in the 321 parasite community, we expect the Pn/Ps to be higher on island species compared to the 322 mainland ones and, therefore, the species origin (i.e., mainland or island) to be significant. In . CC-BY-NC 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in

336
For the ten species (N = 150) for which we generated new data by targeted capture 337 sequencing, an average of 3.3 millions paired-ends reads per individual was generated (Table  338   S1). After mapping, genotyping and cleaning, we analysed 112.5 control and 16.4 immunity 339 genes on average per species, out of the 141 targeted genes (120 control and 21 immunity 340 related genes; Table S3). For the species with whole-genome sequences, we analysed 133 341 control and 20 immunity genes on average per species, out of the 141 targeted genes, and 342 904 and 785 genes on average in the Database-group and Sma3s-group respectively (Table  343 S4). 344

Immunity genes evolving under purifying selection 345
We first focused on a restricted set of genes unambiguously involved in immunity function, 346 namely the BD and TLR genes. At control genes, insular species had, on average, higher 347 Pn/Ps ratios than the mainland ones (0.12 and 0.2 respectively). Model selection based on AICc identified two models as similarly performant at explaining 353 variation in Pn/Ps across species (ΔAICc < 2; Table 2). The first one only includes the gene 354 category. The second one includes the origin (i.e., mainland or island) and gene category 355 without interaction (Table 2). A model selection approach based on simplification with ANOVA 356 identified the latter as the best (Table 2, p < 0.05). In this model, island origin of species is 357 . CC-BY-NC 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for this this version posted November 22, 2021. ; https://doi.org/10.1101/2021.11.21.469450 doi: bioRxiv preprint associated with a greater Pn/Ps (0.12 vs. 0.09; Table 3; p < 0.01). Gene categories 358 corresponding to TLRs and BDs showed a significantly higher Pn/Ps than control genes (Table  359 3; p < 0.001). In all cases, the best models have no interaction between origin and gene 360 categories invalidating the hypothesis of a reduced parasite communities on island (Figure 2). 361   Figure 4). For genes of the Sma3s-group, the category 374 . CC-BY-NC 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for this this version posted November 22, 2021. ; https://doi.org/10.1101/2021.11.21.469450 doi: bioRxiv preprint parameter was also identified by simplification with ANOVA, associated with a reduction of the 375 Pn/Ps of about 0.02 compared to control genes (p < 0.05; Table S9). 376 The alternative statistical approach using the difference between Pn/Ps of immunity genes 377 and control genes ( Pn/Ps) as dependent variable, and species origin as explanatory variable 378 under a PGLS framework lead to similar results. Island was never associated to a statistically 379 higher Pn/Ps (table S2)  Using simulations under frequency dependence selection, model selection identifies the two 400 models as equivalent, first the model with origin and category parameters and the full model 401 (Table 4). However, the full model is not significantly different from the model with origin and 402 category using the method by simplification with ANOVA (Table 4). 403 Using simulations under the overdominance, model selection identifies the model with origin 404 as the best, contrary to the method by simplification with ANOVA which identified the full model 405 therefore including significant interaction between origin and genes category (Table 4). This 406 interaction effect is significant for the MHC II (p < 0.01, Table S12) but not for MHC I. As 407 expected, island species have a significantly lower Pn/Ps in MHC genes compared to 408 mainland species (p < 0.5; except for the full model based on control genes evolving under 409 overdominance Table S12). 410

417
On oceanic islands, the depauperate parasite community is expected to lead to a relaxation 418 of selection on the immune system. In this study, we found support for such an effect, but only 419 on MHC class II genes and under a specific simulation model (i.e, overdominance), which 420 evolves under balancing selection. No effect was detected for MHC class I genes nor for innate 421 immunity genes (TLRs and BDs), evolving under purifying selection. On these gene sets, 422 increased drift effects on island populations limit the efficacy of selection in accordance with 423 the nearly-neutral theory (Ohta 1992). The ability to distinguish between the selective and 424 nearly-neutral processes (relaxed selection due to environmental change vs. drift) could only 425 be achieved by our approach of using random genes (i.e., "control genes'') to estimate the 426 genome-wide effect of potential variation in effective population size between populations. 427

Effects of effective population size variation 428
Our results support the nearly-neutral theory of evolution for those genes under purifying 429 selection, whereby strong genetic drift acting on small island populations reduces the efficacy 430 of natural selection, leading to an increase in non-synonymous nucleotide diversity compared 431 to the mostly neutral, synonymous nucleotide diversity (i.e., Pn/Ps) (Ohta 1992). This is 432 Our results rely on simulations that may be affected by the choice of the parameter values. 471 First, we performed simulations using a fixed effective population size (Ne) estimated from the 472 polymorphism data. Using others values of Ne had a weak impact on the relative difference 473 between island and mainland species for the overdominance type of selection, but had a more 474 noticeable impact for the frequency dependent type of selection ( Figure S4, S5). Secondly, 475 we simulated two types of selection, namely overdominance (Doherty and Zinkernagel 1975) 476 and frequency dependence (Slade and McCallum 1992), but it has been argued that the 477 maintenance of MHC polymorphism could be the result of fluctuating selection (Hill 1991). 478 Additionally, recombination and gene conversion has also been put forward as a mechanism 479 responsible for generating diversity (Spurgin et al. 2011). Therefore, our results for the MHC 480 II, which is based on the relative difference between Pn/Ps of island and mainland species 481 comparing empirical and simulated data, should be taken cautiously as their significance can 482 be dependent on the specific parameters that we used, although we did our best to select a 483 realistic range of parameters. Our comparative population genomics study has investigated the combined effects of drift and 527 selection on immunity genes from island and mainland passerines. The study of synonymous 528 and non-synonymous polymorphism of these genes confirmed that island species, with 529 smaller population sizes than their mainland counterparts, were more impacted by drift, which 530 induces a load of weakly deleterious mutations in their genome. Indeed most of the genes 531 studied here involved in the immune response do not show a statistically different pattern from 532 control genes. Only MHC II genes, involved in the recognition of extracellular pathogens, 533 showed a reduction in their non-synonymous polymorphism in island species. This response, 534 which may be attributed to reduced selection pressures on these genes, could be associated 535 with the suspected reduced parasitic communities on islands. The increased load of 536 deleterious mutations as well as the potential relaxed selection pressures on MHC II support 537 the reduced immune functions of island species, which could be added to the list of other 538 convergent responses of the island syndrome. 539 CC-BY-NC 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for this this version posted November 22, 2021. ; https://doi.org/10.1101/2021.11.21.469450 doi: bioRxiv preprint