# Gene expression is the main driver of purifying selection in large penguin populations

Emiliano Trucchi, Piergiorgio Massa, Francesco Giannelli, Flávia Akemi Nitta Fernandes, Lorena Ancona, Nils Stenseth, Joan Ferrer-Obiol, Josephine Paris, Giorgio Bertorelle, Céline Le Bohec

► To cite this version:

Emiliano Trucchi, Piergiorgio Massa, Francesco Giannelli, Flávia Akemi Nitta Fernandes, Lorena Ancona, et al.. Gene expression is the main driver of purifying selection in large penguin populations.
2023. hal-04307140

HAL Id: hal-04307140

https://hal.science/hal-04307140v1

Preprint submitted on 25 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Gene expression is the main driver of purifying selection in large penguin populations

Emiliano Trucchi<fnref n="1" />®, Piergiorgio Massa<fnref n="2" />, Francesco Giannelli<fnref n="1" />, Flavia A. N. Fernandes<fnref n="1" />,<fnref n="3" />, Lorena Ancona<fnref n="1" />, Nils Chr Stenseth<fnref n="4" />, Joan Ferrer Obiol<fnref n="5" />, Josephine Paris<fnref n="1" />, Giorgio Bertorelle<fnref n="6" />, Céline Le Bohec<fnref n="3" />,<fnref n="7" />,<fnref n="8" />

## Affiliations

## Abstract

Purifying selection is the most pervasive type of selection, as it constantly removes deleterious mutations arising in populations, directly scaling with population size. Highly expressed genes appear to accumulate fewer deleterious mutations between divergent species' lineages, pointing towards gene expression as an additional driver of purifying selection. However, estimates of the effect of gene expression on segregating deleterious variants in natural populations are lacking, as well as an understanding of the relative contribution of population size and gene expression to overall purifying selection pressure. Here, we analyse genomic and transcriptomic data from two natural populations of closely related sister species with different demographic histories, the Emperor (*Aptenodytes forsteri*) and the King penguins (*A. patagonicus*), and demonstrate that purifying selection at the population-level depends on the level of gene expression, with larger effects than population size. Deleterious segregating variants spread less in the population when they are in genes with higher expression rate. Leveraging realistic forward simulations, we estimate that the top 10% of the most highly expressed genes in a genome experience a selection pressure corresponding to an average selection coefficient of -0.1, which decreases to a selection coefficient of -0.01 for the top 50%. Gene expression appears to be a fundamental driver of purifying selection in natural populations, also effective at small population size. We suggest gene expression could be used as a proxy for gene selection coefficients (i.e., distribution of fitness effects), which are notoriously difficult to derive in non-model species under real-world conditions.

## Keywords

Purifying selection, E-R anticorrelation, forward simulations, non-model species, *Aptenodytes*

<footnote n="1">¹Department of Life and Environmental Sciences, Marche Polytechnic University, Ancona, Italy</footnote>

<footnote n="2">²Department of Biological, Geological and Environmental Sciences, Alma Mater Studiorum University of Bologna, Ravenna, Italy</footnote>

<footnote n="3">³Université de Strasbourg, CNRS, IPHC UMR 7178, F-67000, Strasbourg, France</footnote>

<footnote n="4">⁴Centre for Ecological and Evolutionary Synthesis (CEES), Department of Biosciences, University of Oslo, Oslo, Norway</footnote>

<footnote n="5">⁵Department of Environmental Science and Policy, University of Milan, Milan, Italy</footnote>

<footnote n="6">⁶Department of Life Sciences and Biotechnology, University of Ferrara, Ferrara, Italy</footnote>

<footnote n="7">⁷Département de Biologie Polaire, Centre Scientifique de Monaco, Monaco, Monaco</footnote>

<footnote n="8">⁸Centre d'Ecologie Fonctionnelle et Evolutive (CEFE), University of Montpellier, CNRS, EPHE, IRD, Montpellier, France</footnote>

<footnote>©Corresponding author: e.trucchi@univpm.it</footnote>
Introduction

Protein evolution is constrained by purifying selection, which prevents changes in the
underlying gene sequence with a deleterious effect on organismal fitness from spreading in
natural populations. The intensity of purifying selection on weakly deleterious (e.g., most of
nonsynonymous) mutations is directly correlated with the effective size of a population
(Charlesworth 2009, Akashi et al 2012), which, in turn, is determined by species-specific life
history traits and population-specific demographic trajectories (Figuet et al 2016, Chen et al
2017). However, genes with a globally high expression rate across tissues show a slower rate
of accumulation of weakly deleterious substitutions (Duret and Mouchiroud 2000, Pal et al
2001) and lower genetic divergence between closely related species (Zhang and Yang 2015).
This inverse correlation between the rate of evolution and gene expression (E-R
anticorrelation) could be due to the strong selection acting against the toxic accumulation of
misfolded or mis-interacting proteins in cells (Yang et al 2012, Park et al 2013, Wu et al 2022).

Evidence for E-R anticorrelation was found in interspecific comparisons by estimating
nonsynonymous over synonymous (i.e., nearly neutral) fixation rates (dN/dS) in genes with
different expression levels (Zhang and Yang 2015). However, only weak marks of
anticorrelation between gene expression and genetic polymorphism (i.e., segregating
nonsynonymous over synonymous variants, pN/pS) have recently been found in
population-level genetic variation, but only in model organisms in controlled laboratory
populations (Wu et al 2022, Shibai et al 2022). On the other hand, such effects have never
been estimated in natural populations experiencing real-world selection regimes. More
importantly, the interaction between gene expression and effective population size and their
relative contribution to global purifying selection is unknown. In fact, in small populations,
purifying selection is less effective as the contribution of genetic drift to random allele
frequency change is higher: Is such a lower efficiency of purifying selection the same for any
gene regardless of its levels of expression? Or should we expect a larger accumulation of
weakly deleterious variants in highly expressed genes when population size shrinks?

Here, we use two natural populations of closely related sister species, the Emperor and the
King penguins (*Aptenodytes forsteri* and *A. patagonicus*), with different demographic histories
(Trucchi et al 2014, Cristofari et al 2016, 2018), to show that *i*) the intensity of purifying
selection on population-level nonsynonymous variation clearly scales with gene expression, *ii*)
the effect of gene expression largely overwhelms that of population size, and *iii*) for the 50%
most highly expressed genes, the average selection coefficient on any nonsynonymous
variant is stronger than -0.01, while for the top 10% of highly expressed genes, it is stronger
than -0.1. We use high-coverage whole-genome data of 24 individuals per species to estimate
patterns of genetic diversity, and whole transcriptome data of five tissues from three young
individuals per species to estimate global mRNA expression levels. Young age class was
chosen for this analysis as genes broadly expressed in early life stages have been shown to
be the most affected by purifying selection (Cheng and Kirkpatrick 2021).
## Results and Discussion

Both Emperor and King penguins feature single, large and quasi-panmictic populations (Cristofari et al 2016, 2018), but they show different levels of genetic diversity (Fig. 1A), corresponding to their different ecological adaptations and past demographic dynamics (Cristofari et al 2016, 2018, Cole et al 2022). As a consequence of the historically larger effective population size in the Emperor penguin, this species has a higher proportion of segregating variants and a lower proportion of fixed derived variants (Fig. 1B), a lower proportion of segregating nonsynonymous over synonymous variants (pN/pS, Fig. 1B), but only a minor difference in the proportion of fixed nonsynonymous over synonymous differences (dN/dS, Fig. 1B) given the relatively short time since its divergence from the King penguin. Both gene-by-gene estimates of diversity (nucleotide diversity: $\pi$) and expression rate (normalised as transcripts per million, TPM) are highly correlated between the two species (Fig. 1C,D), thus minimising any confounding effect of sequence and expression divergence in our downstream analyses.

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Patterns of genetic diversity and gene expression in Emperor (E, teal) and King (K, gold) penguins. A. Distribution of nucleotide diversity ($\pi$) and Tajima's D in 50 kb genomic windows; B. Proportion of derived alleles as segregating variants (p) or fixed differences (d) at synonymous and nonsynonymous sites (left panel) and estimates of pN/pS and dN/dS (right panel); C. Per gene comparison of nucleotide diversity between King and Emperor penguins; D. Per gene comparison of expression rate between King and Emperor penguins, quantified as transcripts per million (up to TPM = 1100; see Supp. Fig. 5 for the whole expression range).</figcaption></figure>
Purifying selection more efficiently removes nonsynonymous segregating variants in genes while expression rate increases. Average pN/pS across genes binned in 5% percentiles of expression rate shows a declining trend of pN/pS with expression, with the bin average dropping by ca. 80% across the whole range of gene expression in both species (average pN/pS: from > 1.1 in the bottom 5%, to < 0.4 in the top 5% of expression rate; Fig. 2A). Conversely, the difference in pN/pS between the two species (i.e., likely due to the effect of the population size) spans from 2% to 37% across the whole range of gene expression, suggesting that effective population size as a less important determinant of pN/pS in these species. Similar results are obtained if we analyse the expression rate per segregating site (mRNA sequencing coverage per site normalised as count per million reads - CPM) in order to take into account heterogeneous expression rate among exons. Again, counts of nonsynonymous over synonymous variants in bins of 0.05 CPM, from 0 to 5 CPM, are inversely correlated with expression rate (Supp. Fig. 6). The decline of pN/pS with increasing gene expression rate is due to the decreasing count of nonsynonymous variants in highly expressed genes, whereas the count of synonymous variants is stable across the whole gene expression range in both species (Fig. 2B,C). More importantly, the difference in the counts of synonymous variants between the two species is also stable and always significant (Kolmogorov-Smirnov test *p*-value << 0.005), whereas the difference in the counts of nonsynonymous variants decrease with increasing gene expression, with this difference disappearing in the upper 50-60% of gene expression rate (Fig. 2C). This result supports the hypothesis that gene expression overrides effective population size in driving purifying selection. Even if gene expression has been suggested to be one of the causes of non-neutrality in synonymous variants in yeast (Shen et al 2022), or that codon usage bias is more intense in highly expressed genes (Frumkin et al 2018), gene expression rate does not appear to perturb synonymous variation in our datasets from two vertebrate species.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Increasing purifying selection with gene expression in Emperor (teal) and King (gold) penguins.
Distribution of pN/pS (A), and average number of synonymous (B) and nonsynonymous (C) segregating variants (normalised by coding sequence - CDS - length) in genes binned by 5% percentiles of expression rate (normalised as TPM). Median (solid white line) and mean (white triangle) is shown in each boxplot. Statistical significance for the difference in the distribution of synonymous and nonsynonymous variants (Kolmogorov-Smirnov test *p*-value) per percentile between the two species is shown.</figcaption></figure>
Purifying selection more efficiently prevents nonsynonymous segregating variants from increasing in frequency in genes with higher expression rate. The derived allele frequency spectrum of nonsynonymous variants with expression rate higher than 0.3 CPM is depleted in medium-high frequency categories, while there is no difference in the derived allele frequency spectrum of synonymous variants across the whole expression range (Fig. 3). Changing the arbitrary threshold to discriminate between low and high expression rate, or using more than two categories of expression rate (low: < 0.3 CPM, medium: 0.3-2 CPM, high: > 2 CPM) does not change the observed pattern (Supp. Fig. 7). The pattern holds when all nonsynonymous and synonymous variants are used in the allele frequency spectrum estimate (Supp. Fig. 7) as well as when one nonsynonymous and one synonymous variant are randomly sampled from each gene (Fig. 3), thus excluding the possibility that few genes with many variants (i.e., pseudoreplication) drive our observation.

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Nonsynonymous variants in highly expressed genes segregate at lower frequency in Emperor (teal) and King (gold) penguins. Site frequency spectra of ten random resampling (95% distribution) of one nonsynonymous (upper panels) and one synonymous (lower panels) variant per gene with lower (dark shade) or higher (light shade) than 0.3 CPM mRNA expression. The relative frequency of each count class is log10 transformed. Note that nonsynonymous variants in highly expressed genes show a higher frequency at low derived allele count classes than nonsynonymous variants in lowly expressed genes.</figcaption></figure>
Purifying selection in the top 10% of highly expressed genes largely exceeds the effect of 100,000-individuals effective population size. In simulated populations, under either Wright-Fisher or more realistic non Wright-Fisher models, median pN/pS across genes declines from 1.8 to 0.9, while population size increases from 1,000 to 100,000 individuals (Fig. 4). Such values of pN/pS are much higher than the mean or median values observed for genes in the top 10% of expression rate (Fig. 2A). In these models, the effect of population size on purifying selection was explored by simulating a set of realistic values for mutation and recombination rate, synonymous to nonsynonymous ratio, selection and dominance coefficient distributions, coding sequence length and gene numbers. In particular, new mutations were given a selection and dominance coefficient (*h*-mix) based on a nearly neutral prior distribution (Kim et al 2017, Kyriazis et al 2021), meaning that most of the mutations are weakly deleterious. To reproduce the pN/pS values that we observed in the highly expressed genes (i.e., mean pN/pS < 0.4 in the top 10%) in both penguin species, we designed a more extreme selection scenario: all nonsynonymous mutations appearing in a gene were given a fixed selection coefficient of -0.1, -0.01 or -0.001 (100 replicated genes per selection coefficient) and a dominance coefficient derived from the *hs* relationship (Henn et al 2016). In models with a selection coefficient of -0.01, pN/pS decreases below 0.6 (Fig. 4) while a selection coefficient of -0.1 results in pN/pS below 0.4, as observed in genes in the top 50% and 10% of expression rate, respectively. Such strong selection coefficients are expected to be effective even when the population size is small (i.e., s >> 1/N, per N = 1,000), however we observe more variance in simulations with smaller population sizes (Supp. Fig. 8).

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Population size and gene-specific extreme selection coefficient explain low observed pN/pS values in simulations. Distribution of pN/pS across 1000 genes simulated under nonWrightFisher (nonWF, left, solid border) and WrightFisher (WF, right, dashed border) models with effective population size from 1,000 to 100,000 (darker grey background) and across 100 genes with selection coefficient from -0.001 to -0.1 (lighter grey background). Note that the dominance coefficient is set according to the *h*-mix or *hs* models in simulations testing different population sizes or selection coefficients, respectively. More efficient purifying selection in nonWF models, where effective population size tends to be lower than in WF models, can be explained by the fact that, in such models, individuals with high fitness can survive and reproduce for multiple generations (Haller and Messer 2019).</figcaption></figure>
Gene expression can be used as a proxy of the distribution of gene selection coefficients in natural populations of non-model species. Variants with highly deleterious effects on individual fitness are expected to be immediately lost in natural populations. Consistent with this expectation, the highly deleterious variants (< 1000 HIGH effect SNPs per species) predicted by SNPeff (Cingolani et al 2012) in each penguin population show a much lower average expression level than weakly deleterious (MODERATE effect) and nearly-neutral (LOW effect and synonymous) variants (Tab. 1). This observation means that HIGH effect variants are mainly present in lowly expressed genes with limited impact on fitness. Expression level is even lower in the very few fixed differences with HIGH effect (Tab. 1), thus supporting our hypothesis. As site-specific expression of highly deleterious variants (mainly start/stop codons loss/gain and splice acceptor/donor variants) could be biassed in mature mRNA sequencing, we also estimated the expression of highly deleterious variants as the expression of the gene they belong to. Even applying a rather conservative test, the expression of genes with predicted highly deleterious variants is on average three times lower than all genes (Kolmogorov-Smirnov test *p*-value = 0.00095) in the King penguin and slightly lower, even if not significant, in the Emperor penguin. As previously suggested for the distribution of dominance coefficients in a model plant species (Huber et al 2018), gene expression should be taken into account when using predictions of fitness effects and, more generally, when using such predictions to calculate the genetic load in populations of conservation concern (Bertorelle et al 2022). In fact, predicted highly deleterious variants could be on lowly expressed genes, thus with little contribution to individual or population fitness.

<table><caption>**Table 1. Expression rate by predicted fitness effect.** LOWsyn: low effect and synonymous; MDR: moderate effect; HIGH: high effect.</caption><thead><tr><th>PREDICTED FITNESS EFFECT</th><th>LOWsyn</th><th>LOWsyn</th><th>MDR</th><th>MDR</th><th>HIGH</th><th>HIGH</th></tr><tr><th>Species</th><th>Emperor</th><th>King</th><th>Emperor</th><th>King</th><th>Emperor</th><th>King</th></tr></thead><tbody><tr><td>Total variants</td><td>73501</td><td>57088</td><td>47183</td><td>41352</td><td>934</td><td>840</td></tr><tr><td>→Average (stdev) Z-normalised CPM</td><td>1.5 (5.02)</td><td>1.3 (4.56)</td><td>0.78 (2.85)</td><td>0.78 (2.86)</td><td>0.24 (1.97)</td><td>0.14 (0.88)</td></tr><tr><td>Fixed differences</td><td>1846</td><td>3500</td><td>1166</td><td>2229</td><td>16</td><td>44</td></tr><tr><td>→Average (stdev) Z-normalised CPM</td><td>1.41 (4.12)</td><td>1.22 (4.14)</td><td>1.4 (6.31)</td><td>0.93 (3.39)</td><td>0.08 (0.54)</td><td>0.02 (0.57)</td></tr><tr><td>Average count of derived alleles in segregating variants</td><td>6.17</td><td>6.25</td><td>4.75</td><td>4.73</td><td>5.14</td><td>5.48</td></tr></tbody></table>

## Concluding remarks

Overall, our study provides evidence that gene expression is a critical driver of purifying selection in natural populations, and to a higher extent than population size. About half of the genes in a genome, which are likely responsible for basic cellular and molecular functions
(Boyle et al 2017), are under a strong selective constraint preventing deleterious sequence changes even when population size declines to about 1,000 individuals. Below this order of magnitude, random effects could prevail in the population evolutionary trajectory. Importantly, gene expression can be used as a proxy of the gene selection coefficient, which is notoriously difficult to study in natural populations of non-model species (Huber et al 2017). Gene expression data are easier to collect than selection coefficients and are usually highly conserved across closely related species (Fig. 1D), so that they can be used to refine estimates of genetic load (Bertorelle et al 2022) in natural populations of conservation concern.

**Data availability**

Genomic and transcriptomic raw reads are publicly available at ENA database with Project accession number PRJEB64484 and XXX (genomic raw reads to be submitted).

**Acknowledgements**

We thank Fabrizio Mafessoni, Paolo Gratton and Robin Cristofari for helpful comments and suggestions on the analyses, and Federica Pirri for RNA data production. ET was supported by PNRA_16 00164 ("Programma Nazionale di Ricerca in Antartide"). Bando PNRA 5 aprile 2016, n. 651. - Linea B "Genomica degli adattamenti estremi alla vita in Antartide"). The study was supported by the Institut Polaire Français Paul-Emile Victor (IPEV) within the framework of the Program 137-ANTAVIA (PI: CLB), by the Centre Scientifique de Monaco with additional support from the LIA-647 and RTPI-NUTRESS (CSM/CNRS-UNISTRA), by the Centre National de la Recherche Scientifique (CNRS) through the Programme Zone Atelier de Recherches sur l'Environnement Antarctique et Subantarctique (ZATA). This study was approved by the French ethics committee (last: APAFIS#29338-2020070210516365) and the French Polar Environmental Committee, and permits to handle animals and access breeding sites were delivered by the "Terres Australes et Antarctiques Françaises" (TAAF). The sequencing service was provided by the Norwegian Sequencing Centre (www.sequencing.uio.no), a national technology platform hosted by the University of Oslo and supported by the "Functional Genomics" and "Infrastructure" programs of the Research Council of Norway and the Southeastern Norway Regional Health Authorities. Bioinformatic analyses were performed on the HPC clusters of the Department of Life and Environmental Sciences ("HappyComputing@DiSVA"), Marche Polytechnic University, and of the Department of Life Sciences and Biotechnology, University of Ferrara.

**Author contributions**

ET designed the study and secured the fundings together with GB, extracted the DNA samples, analysed the genomic and transcriptomic data, and wrote the manuscript; PM and FG built, performed and analysed the genomic simulations, and wrote the relevant section in the Methods; NCS contributed to the pilot study leading to this work; MB, CC collected the samples for RNA analyses and contributed to the manuscript; FANF, LA, JFO, JP and GB discussed the results and contributed to the manuscript; CLB coordinated the project logistics and the samples collection and the associated fundings, discussed the results and contributed to the manuscript. All authors contributed to the development of the paper.
## References

Akashi H, Osada N, Ohta T (2012). Weak selection and protein evolution. *Genetics*, **192**, 15-31

Bertorelle G, Raffini F, Bosse M et al (2022). Genetic load: genomic estimates and applications in non-model animals. *Nature Reviews Genetics*, **23**(8), 492-503.

Boyle EA, Li YI, Pritchard JK (2017). An expanded view of complex traits: from polygenic to omnigenic. *Cell*, **169**(7), 1177-1186.

Charlesworth B (2009). Fundamental concepts in genetics: Effective population size and patterns of molecular evolution and variation. *Nature Reviews Genetics*, **10**(3), 195-205.

Chen J, Glémin S, Lascoux M (2017). Genetic diversity and the efficacy of purifying selection across plant and animal species. *Molecular biology and evolution*, **34**(6), 1417-1428.

Cheng C, Kirkpatrick M (2021). Molecular evolution and the decline of purifying selection with age. *Nature communications*, **12**(1), 1-10.

Cingolani P, Platts A, Wang LL et al (2012). A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3. fly, 6(2), 80-92.

Cole TL, Zhou C, Fang M et al (2022). Genomic insights into the secondary aquatic transition of penguins. *Nature Communications*, **13**(1), 3912.

Cristofari R, Bertorelle G, Ancel A et al (2016). Full circumpolar migration ensures evolutionary unity in the Emperor penguin. *Nature communications*, **7**(1), 11842.

Cristofari R, Liu X, Bonadonna F, Cherel Y et al (2018). Climate-driven range shifts of the king penguin in a fragmented ecosystem. *Nature Climate Change*, **8**(3), 245-251.

Duret L, Mouchiroud D (2000). Determinants of substitution rates in mammalian genes: expression pattern affects selection intensity but not mutation rate. *Molecular biology and evolution*, **17**(1), 68-70.

Figuet E, Nabholz B, Bonneau M et al (2016). Life history traits, protein evolution, and the nearly neutral theory in amniotes. *Molecular biology and evolution*, **33**(6), 1517-1527.

Frumkin I, Lajoie MJ, Gregg CJ et al (2018). Codon usage of highly expressed genes affects proteome-wide translation efficiency. *Proceedings of the National Academy of Sciences*, **115**(21), E4940-E4949.
Henn BM, Botigué LR, Peischl S et al (2016). Distance from sub-Saharan Africa predicts mutational load in diverse human genomes. *Proceedings of the National Academy of Sciences*, 113(4), E440-E449.

Huber CD, Kim BY, Marsden CD, Lohmueller KE (2017). Determining the factors driving selective effects of new nonsynonymous mutations. *Proceedings of the National Academy of Sciences*, 114(17), 4465-4470.

Huber CD, Durvasula A, Hancock AM, Lohmueller KE (2018). Gene expression drives the evolution of dominance. *Nature communications*, 9(1), 2750.

Kim BY, Huber CD, Lohmueller KE (2017). Inference of the distribution of selection coefficients for new nonsynonymous mutations using large samples. *Genetics*, 206(1), 345-361.

Kyriazis CC, Wayne RK, Lohmueller KE (2021). Strongly deleterious mutations are a primary determinant of extinction risk due to inbreeding depression. *Evolution letters*, 5(1), 33-47.

Pál C, Papp B, Hurst LD (2001). Highly expressed genes in yeast evolve slowly. *Genetics*, 158(2), 927-931.

Park C, Chen X, Yang JR, Zhang J (2013). Differential requirements for mRNA folding partially explain why highly expressed proteins evolve slowly. *Proceedings of the National Academy of Sciences*, 110(8), E678-E686.

Shen X, Song S, Li C, Zhang J (2022). Synonymous mutations in representative yeast genes are mostly strongly non-neutral. *Nature*, 606(7915), 725-731.

Shibai A, Kotani H, Sakata N, Furusawa C, Tsuru S (2022). Purifying selection enduringly acts on the sequence evolution of highly expressed proteins in Escherichia coli. *G3*, 12(11), jkac235.

Trucchi E, Gratton P, Whittington JD et al (2014). King penguin demography since the last glaciation inferred from genome-wide data. *Proceedings of the Royal Society B: Biological Sciences*, 281(1787), 20140528.

Yang JR, Liao BY, Zhuang SM, Zhang J (2012). Protein misinteraction avoidance causes highly expressed proteins to evolve slowly. *Proceedings of the National Academy of Sciences*, 109(14), E831-E840.

Wu Z, Cai X, Zhang X et al (2022). Expression level is a major modifier of the fitness landscape of a protein coding gene. *Nature Ecology & Evolution*, 6(1), 103-115.

Zhang J, Yang JR (2015). Determinants of the rate of protein sequence evolution. *Nature Reviews Genetics*, 16(7), 409-420.
# Extended methods and supplementary information

## 1. Genomic data

### 1.1 Samples collection and storage

Genomic DNA extractions of 24 King penguins (*Aptenodytes patagonicus* from South Georgia, Crozet archipelago, and Heard island), and 24 Emperor penguins (*Aptenodytes forsteri* from Terre Adélie and Dronning Maud Land, Antarctica), were selected from samples used in Cristofari et al (2016) and Cristofari et al (2018). In addition, four Gentoo penguin (*Pygoscelis papua* from Crozet archipelago) and two Adelie penguins (*Pygoscelis adeliae* from Terre Adélie, Antarctica, 2007) were collected during the field campaigns of the French IPEV programme 137 in 2017 and 2007, respectively (Supp. File 1). Samples were stored in ETOH (muscle biopsy) or in Queen Lysis buffer (blood samples), frozen at -80°C from the field to the lab.

### 1.2 DNA extraction, pooled library preparation and sequencing

DNA extraction was performed using DNeasy Blood and Tissue kit (Qiagen) following manufacturer's instructions. Whole genome sequencing libraries were prepared and sequenced at the Norwegian Sequencing Centre, Oslo, using Illumina pcr-free single or dual-indexing kits. To minimise batch effects, genomic samples from different species and different localities were randomised in six libraries and sequenced on 1-3 lanes of Illumina HiSeq2500 and HiSeq4000 aiming at 20X coverage depth. Raw reads are publicly available at ENA database (XXX, to be made available).

### 1.3 Variant calling, filtering and annotation

After Illumina adapters trimming and quality filtering with Trimmomatic (Bolger et al 2014) and URQT (Modolo and Lerat 2015), respectively, fastq reads were mapped to the *Aptenodytes forsteri* reference genome (ASM69914v1; RefSeq assembly accession: GCF_000699145.1) using *bwa mem* (v0.7.15; Li 2013), converted to bam files and sorted with *samtools* (v0.1.19; Li et al 2009), keeping only reads with phred-scaled mapping quality higher than 10. Duplicated reads were removed with *picard-tools* (v1.98) and bam files were assigned to individual samples by adding ID read groups with picard-tools.

Small variants (SNPs, indel and MNPs) were called with *freebayes* (Garrison and Marth 2012) using reference genome scaffolds longer than 100Kb (476 in total), all samples grouped per species (*--populations* flag), and a minimum phred-scaled mapping quality of 20. Resulting vcf files (one per scaffold) were then filtered: i) MNPs were first broken down into SNPs using *vcfallelicprimitives* script in *vcflib* (with -k -g flags; Garrison et al 2022); ii) variants were then filtered for quality (QUAL > 30), strand bias (SAF > 0 & SAR > 0), read placement bias (RPL > 0 & RPR), and type of variants (TYPE = snp) with *vcffilter* script in *vcflib*; iii) SNPs were finally filtered using *vcftools* (v0.1.15; Danecek et al 2011) for minimum coverage depth of 3 reads per individual (*--minDP* 3; individual genotypes discarded if below threshold), mean maximum coverage depth of 50 (*--max-meanDP* 50 which is ca. three times the average individual coverage depth; whole locus discarded if above threshold), and retained only if biallelic across all samples. A total of 44 sex-linked scaffolds were identified and removed from the dataset by running *samtools idxstats* on all bam files and the SATC (Nursyifa et al 2022)
Rscript on the resulting data. SNPs were annotated using SNPeff (Cingolani et al 2012) and the GCF_000699145.1 genome annotation. Annotated vcf files are available upon request. Mean coverage depth of filtered SNPs per allele was 6.93X, standard deviation 1.07X.

## 1.4 SNPs polarisation in ancestral and derived alleles

Annotated vcf files were parsed using a custom python script (vcf2missenseFreq.2d.py; https://github.com/emitruc/ExpressionLoad). For each SNP position we summed up the joint derived allele counts for King and Emperor penguin samples using Adelie and Gentoo penguins samples as outgroups. The resulting summary table (daf.joint; Supp. Tab. 1) included the following information:

- counts of derived alleles and total alleles in Emperor penguin (der1; tot1), King penguin (der2; tot2), and Adélie and Gentoo penguins (der_out; tot_out) samples;

- average allele coverage (avgCov);

- ancestral (ref) and derived (alt) allele;

- genomic site type (vartype) based on the first annotation by SNPeff as missense if containing the word 'missense', synonymous if containing the word 'synonymous', intergenic if containing the word 'intergenic', intronic if containing the word 'intron', else otherwise;

- genomic site predicted effect (effect) based on the first annotation description as HIGH, MODERATE, LOW, MODIFIER;

- a flag whether the polymorphic site was originally called as MNP or SNP by freebayes (fyellow bellied toadlagQual: haplo, snp);

- a flag to track how the derived allele was called (flagPol) on the basis of the options shown in Supplementary Figure 1.

<table><caption>**Supplementary Table 1.** Example of SNPs recorded in the summary table daf.joint</caption><thead><tr><th>scaffold</th><th>position</th><th>flagPol</th><th>flagQual</th><th>der1</th><th>tot1</th><th>der2</th><th>tot2</th><th>der_out</th><th>tot_out</th><th>avgCov</th><th>ref</th><th>alt</th><th>vartype</th><th>effect</th></tr></thead><tbody><tr><td>NW_008793941.1</td><td>85516</td><td>InFixAnc</td><td>snp</td><td>0</td><td>48</td><td>0</td><td>44</td><td>6</td><td>8</td><td>7.2</td><td>G</td><td>A</td><td>intron</td><td>MODIFIER</td></tr><tr><td>NW_008793941.1</td><td>85528</td><td>unfolded</td><td>snp</td><td>3</td><td>48</td><td>46</td><td>46</td><td>0</td><td>8</td><td>6.6</td><td>G</td><td>A</td><td>intron</td><td>LOW</td></tr><tr><td>NW_008793941.1</td><td>85558</td><td>allFix</td><td>snp</td><td>0</td><td>48</td><td>0</td><td>46</td><td>8</td><td>8</td><td>6.4</td><td>T</td><td>C</td><td>synonymous</td><td>LOW</td></tr><tr><td>NW_008793941.1</td><td>85657</td><td>unfolded</td><td>snp</td><td>0</td><td>48</td><td>1</td><td>48</td><td>0</td><td>12</td><td>6</td><td>C</td><td>T</td><td>synonymous</td><td>LOW</td></tr><tr><td>NW_008793941.1</td><td>85734</td><td>unfolded</td><td>snp</td><td>2</td><td>48</td><td>0</td><td>48</td><td>0</td><td>12</td><td>6.9</td><td>C</td><td>T</td><td>intron</td><td>MODIFIER</td></tr><tr><td>NW_008793941.1</td><td>85746</td><td>allFix</td><td>snp</td><td>0</td><td>48</td><td>0</td><td>48</td><td>12</td><td>12</td><td>6.8</td><td>G</td><td>A</td><td>intron</td><td>MODIFIER</td></tr><tr><td>NW_008793941.1</td><td>85762</td><td>unfolded</td><td>snp</td><td>5</td><td>48</td><td>0</td><td>48</td><td>0</td><td>10</td><td>6.6</td><td>C</td><td>T</td><td>intron</td><td>MODIFIER</td></tr><tr><td>NW_008793941.1</td><td>85788</td><td>InFixAnc</td><td>snp</td><td>0</td><td>48</td><td>0</td><td>46</td><td>4</td><td>8</td><td>7.1</td><td>A</td><td>C</td><td>intron</td><td>MODIFIER</td></tr><tr><td>NW_008793941.1</td><td>85790</td><td>unfolded</td><td>snp</td><td>0</td><td>48</td><td>15</td><td>46</td><td>0</td><td>8</td><td>7.1</td><td>T</td><td>C</td><td>intron</td><td>MODIFIER</td></tr><tr><td>NW_008793941.1</td><td>85800</td><td>unfolded</td><td>snp</td><td>2</td><td>48</td><td>0</td><td>48</td><td>0</td><td>8</td><td>6.9</td><td>C</td><td>T</td><td>intron</td><td>MODIFIER</td></tr><tr><td>NW_008793941.1</td><td>85804</td><td>unfolded</td><td>snp</td><td>1</td><td>48</td><td>0</td><td>48</td><td>0</td><td>8</td><td>6.9</td><td>C</td><td>T</td><td>intron</td><td>MODIFIER</td></tr><tr><td>NW_008793941.1</td><td>85829</td><td>unfolded</td><td>snp</td><td>0</td><td>48</td><td>11</td><td>48</td><td>0</td><td>8</td><td>6.7</td><td>T</td><td>C</td><td>intron</td><td>MODIFIER</td></tr><tr><td>NW_008793941.1</td><td>85836</td><td>allFix</td><td>snp</td><td>0</td><td>48</td><td>0</td><td>48</td><td>8</td><td>8</td><td>6.6</td><td>A</td><td>G</td><td>intron</td><td>MODIFIER</td></tr><tr><td>NW_008793941.1</td><td>85845</td><td>allFix</td><td>snp</td><td>0</td><td>48</td><td>0</td><td>46</td><td>8</td><td>8</td><td>6.8</td><td>A</td><td>G</td><td>intron</td><td>MODIFIER</td></tr><tr><td>NW_008793941.1</td><td>85864</td><td>InFixAnc</td><td>snp</td><td>0</td><td>48</td><td>0</td><td>44</td><td>4</td><td>10</td><td>6.1</td><td>T</td><td>C</td><td>intron</td><td>MODIFIER</td></tr><tr><td>NW_008793941.1</td><td>85870</td><td>unfolded</td><td>snp</td><td>0</td><td>46</td><td>5</td><td>46</td><td>0</td><td>10</td><td>6.6</td><td>G</td><td>A</td><td>intron</td><td>MODIFIER</td></tr></tbody></table>
<figure><img src="image_6.png" /><figcaption>**Supplementary Figure 1.** Schematic representations of the potential distribution of ancestral and derived alleles at a biallelic site across two ingroup and one outgroup samples. Ingroup1: Emperor penguin; ingroup2: King penguin; outgroup: *Pygoscelis* (Adélie + Gentoo) penguins. Blue and yellow circles represent copies of the ancestral and derived alleles. The number of allele copies per sample is only given as an example. Missing data are represented as no circles. The yellow star represents the branch and timing the mutation to the derived allele likely occurred at. For ALLFIX and INFIXOUTMISS configurations, the allele of ingroup1 is arbitrarily considered as ancestral. Ingroup1 or ingroup2 missing data configurations are not shown. This algorithm is implemented in the python script *vcf2missenseFreq.2d.py* (available at https://github.com/emitruc/ExpressionLoad).</figcaption></figure>

**1.5 Distribution of derived alleles frequencies**

After removing monomorphic sites across Emperor and King penguin samples (e.g., INFIXANC and ALLFIX configurations in Supp. Fig. 1) from the dataset, we checked whether the derived alleles distribution was in line with basic expectations from population genetics. To avoid downstream normalisation, we selected only sites without missing data in the target species. In addition, we choose sites which were properly polarised (flagPol=UNFOLDED) based on at least four alleles present in the outgroup and with coverage depth range between 6X and 8X, i.e., one standard deviation (ca. 1X) from the mean (ca. 7X). Such a narrow coverage range was applied to mitigate as much as possible the inclusion of sex-chromosome related regions, which could have been incorrectly assembled within autosomal scaffolds, and multiple-copy regions, which were not properly assembled in the reference genome sequence (see Supp. Fig. 2 for a comparison at different coverage ranges). Both types of regions are also
characterised by marked deviation from Hardy-Weinberg equilibrium (HWE). Applying a coverage threshold as a filter, we could retain sites which are not in HWE due to other processes (e.g., population structure, selection). Using all intergenic sites with no missing data for both species, we summarised the joint derived allele frequency (DAF) distribution (Supp. Fig. 2). As expected, a very small proportion (note the log scale of the heatmap in Supp. Fig. 2) of derived alleles appear as segregating in both species, with a minimal covariance in the low and in the very high frequency classes. The latter suggests some mis-polarisation due to sites hit by multiple mutations (see below). Nevertheless, when setting the coverage range between 6X and 8X, co-segregating derived allele frequencies largely appear as uncorrelated between the two species as expected in case of incompletely sorted ancient variation.

<figure><img src="image_7.png" /><figcaption>**Supplementary Figure 2.** Joint derived allele frequency distributions in King and Emperor penguins samples at different ranges of coverage.</figcaption></figure>

DAF were separately summarised according to the predicted annotation as intergenic, intronic, synonymous and nonsynonymous in each species (Supp. Fig. 3). In general, the shape of observed DAF is consistent with the expectations from population genetics theory. The slight (note the log scale on the y-axis in Supp. Fig. 3) increase in very high frequency variants has been commonly observed (also in human population data; Marchi and Excoffier 2020) and it is likely due to mis-polarisation of ancestral/derived alleles at sites hit by multiple mutations (Hernandez et al 2007), rather than to migration from a “ghost” population (Marchi and Excoffier 2020). DAF from intergenic sites represent the best approximation to neutrality, with a shape depending on past population demography only. On the contrary, DAF from missense sites are expected to be enriched in low frequency and depleted in medium-high frequency classes as deleterious alleles are less likely to increase in frequency in the population due to negative selection. Such a pattern is clearly appearing in both Emperor and King penguins. Synonymous sites DAF show small deviations from neutrality, likely due to linked selection within missense sites in exons. On the other hand, DAF from intronic and intergenic sites are fully overlapping with each other, suggesting very limited linked selection spanning from exons to introns.
<figure><img src="image_8.png" /><figcaption>**Supplementary Figure 3.** Derived allele frequency distributions for different types of variants (Intergenic, intronic, synonymous and missense as annotated by SNPeff)</figcaption></figure>

**1.6 Population statistics**

Nucleotide diversity ($\pi$) and Tajima's D were estimated in non-overlapping windows of 10 thousand base pairs (kbp) in Emperor and King penguin samples using vcftools from 432 scaffolds. As we are interested in the species-level long-term adaptation processes, we included in our samples individuals from multiple locations to get an accurate representation of the whole species diversity. Previous studies, based on reduced-representation genome sequencing of hundreds of individuals per species, described both Emperor and King penguin as quasi-panmictic species, with only one very large population each and very little differentiation among colonies (Cristofari et al 2016, 2018). By estimating Weir and Cockerham (1984) $F_{ST}$ in 10 kbp non-overlapping windows using vcftools (Supp. Fig. 4 inset), we confirmed previous results finding negligible genetic differentiation between King penguin samples from Crozet and Heard, and South Georgia ($F_{ST}$ mean 0.003, std 0.024) and between Emperor penguin samples from Terre Adélie and Dronning Maud Land ($F_{ST}$ mean 0.004, std 0.020). Principal component analysis (SNPRelate - Zheng et al 2012, 500 bp pruning for linkage disequilibrium) on a subset of ca. 20,000 SNPs from scaffold NW_008796188 showed no genetic structure in any of the two species (Supp. Fig. 4).
<figure><img src="image_9.png" /><figcaption>**Supplementary Figure 4.** No signature of genetic differentiation among distant colonies of Emperor (left) and King penguin (right). Main plots show the results of Principal Component Analysis run in SNPRelate using 20,256 and 18,256 SNPs not closer than 500 bp from scaffold NW_008796188, respectively. Histograms of *F*<sub>ST</sub> estimated in non-overlapping 10 kbp windows across the whole genome are shown as insets. DML: Dronning Maud Land; ADE: Terre Adélie; CRO: Crozet; HEA: Heard; GEO: South Georgia.</figcaption></figure>
## 2. Gene expression data

### 2.1 Sample collection and storage

During June-September 2016 field campaigns at Dumont d'Urville Station in Antarctica and at the Alfred Faure station in Crozet archipelago, samples of five different tissues (brain, liver, kidney, skin and muscle) were collected from three freshly-predated 3-7 months-old chicks from natural populations of Emperor and King penguins, respectively (Pirri et al in preparation). All tissue samples were collected immediately after death, directly fixed in RNA later (Applied Biosystems, Warrington, UK) and frozen at -80°C until RNA extraction.

### 2.2 RNA extraction, pooled library preparation and sequencing

Total RNA was isolated from 40 mg of each tissue sample by a standard laboratory-based chloroform extraction after homogenization in 500 ul of TRIzol® reagent (Invitrogen, ThermoFisher Scientific). Samples were added to 100 µl of chloroform, vortexed, and centrifuged at 12,000×g for 15 min at 8°C; the upper aqueous phase was collected and transferred to a new tube for precipitation with isopropanol by centrifugation at 12,000 g for 10 min at 8 °C; the RNA pellet was washed with 75% ethanol and centrifuged at 7,500×g for 5 min at 8 °C, the ethanol removed and the RNA pellet resuspended in RNase free water to be stored at -80°C. As TRIzol-based extraction from skin and muscle yielded poor RNA quality and quantity, likely due to large amount of proteins, connective tissue, and collagen in these tissues, RNA from these two tissues was extracted using the RNeasy Fibrous Tissue Mini Kit (Qiagen) according to the manufacturer's instructions. Also in this case, isolated RNA was dissolved in RNase free water and stored at -80°C. Concentration and purity (i.e., the A260/A280 ratio) of each RNA sample was assessed by Nanodrop 2000 (Thermo Fisher Scientific) and Qubit 4.0 fluorometer (ThermoFisher Scientific) while RNA integrity was evaluated by capillary electrophoresis on Agilent 2100 Bioanalyzer (Agilent technologies, Santa Clara, CA). As the target of our study was to estimate the global level of gene expression (across tissues), a total of six RNA pools (three pools of five tissues per three individuals for each species) were assembled starting from 15 RNA samples per species, after concentration was normalised. RNA-seq library preparation and sequencing was carried out by BMR Genomics Service (Padova, Italy). Libraries were synthesised using the TruSeq Stranded mRNA Sample Prep kit (Illumina, San Diego, CA), according to the manufacturer's instructions. Poly-A mRNA was fragmented for 3 minutes at 94°C, and each purification step was carried out with 1 × Agencourt AMPure XP beads. Paired-end sequencing (100 bp from each end) was then performed on the Novaseq 6000 (Illumina, San Diego, CA) at a sequencing depth of 100 million reads per library. Raw reads are publicly available at ENA database as one pool of reads per species (Project ID: PRJEB64484, sample accession ID King penguin: ERS16093259; sample accession ID Emperor penguin, ERS16093260).

### 2.3 RNA mapping, base-pair and gene expression rate estimates

RNAseq reads from the two penguin species were mapped to the same reference genome (i.e., *A. forsteri* reference genome - RefSeq assembly accession: GCF_000699145.1) as for the genomic data. In particular, after standard filtering and trimming with Trimmomatic, RNA reads were mapped using STAR v.2.7.9a (Dobin et al 2013) and resulting bam files indexed with SAMtools. From bam files, counts of reads overlapping each gene were estimated with HTseq (Anders et al 2015) using the available genome annotation for GCF_000699145.1 reference
genome. Multi-mapped and overlapping multiple expression features reads were discarded.
For each gene, genomic coordinates, exon and CDS length were extracted from the
annotation of GCF_000699145.1 reference genome with the following bash commands and
merged with the RNAseq expression counts from both Emperor and King penguin samples:

while read i; do echo -ne $i '\t' ; grep -e 'gene='$i';' GCF_000699145.1_ASM69914v1_genomic.gff | grep -e ' gene ' | cut -f 1,4,5 | sort -k1,1 -k2n | bedtools merge ; done < gene.list > geneCoord

while read i e ; do echo -ne $i '\t' ; grep -e 'gene='$i';' GCF_000699145.1_ASM69914v1_genomic.gff | grep exon | cut -f 1,4,5 | sort -k1,1 -k2n | bedtools merge | awk '{sum+=$3-$2} END {print sum}'; done < gene.list > geneLength

while read i e ; do echo -ne $i '\t' ; grep -e 'gene='$i';' GCF_000699145.1_ASM69914v1_genomic.gff | grep -e " CDS" | cut -f 1,4,5 | sort -k1,1 -k2n | bedtools merge | awk '{sum+=$3-$2} END {print sum}'; done < gene.list > cdsLength
# 3. Testing the effect of gene expression on purifying selection

## 3.1 Gene-based analyses

Counts of synonymous and missense polymorphic sites were summarised per gene using a custom python script (FinalPipeline.ipynb; https://github.com/emitruc/ExpressionLoad). After applying the stringent filters for coverage (from 6X to 8X) and missing data (no missing genotype), we used the genomic coordinates of all genes to subset the list of SNPs in the *daf.joint* dataset and count the number of synonymous and missense polymorphic sites per gene. RNAseq counts per gene from HTseq were then normalised by CDS length and total number of reads as transcript per million (TPM), separately per Emperor and King penguin samples (Supp. Fig. 5).

<figure><img src="image_10.png" /><figcaption>**Supplementary Figure 5.** Expression rate as transcripts per million (TPM) per gene in King (yellow) and Emperor (blue) penguins. Whole range of expression (left panel) and up to 6000 TPM (right panel).</figcaption></figure>

Counts of synonymous and missense polymorphic sites per gene were also normalised by CDS length. To estimate the intensity of purifying selection in each gene we used the ratio of missense (N) over synonymous (S) polymorphic sites (pN/pS).

We used these data to test two hypotheses:

First, if purifying selection is driven by gene expression rate, we should observe a decline in pN/pS with increasing expression and such decline should be determined by a corresponding decline in missense polymorphism. To test this hypothesis, we studied the distribution of normalised synonymous and missense counts and pN/pS values in genes grouped by 5% percentiles of TPM.

Our second question is more complex and concerns the relative weight of population size and gene expression, respectively, in driving purifying selection. When comparing two populations with different historical effective population size, we expect the smaller
population to show lower diversity at both neutral and deleterious sites, but higher pN/pS because of larger drift which, in turn, reduces the efficacy of purifying selection. If population size is the main driver of purifying selection, we expect that the difference in diversity, and of pN/pS, between the two populations of different size will be the same across the whole range of gene expression. If gene expression is instead the main driver of purifying selection, we expect that the difference in diversity between the two populations of different size will decline with increasing gene expression rate at deleterious sites but not at neutral sites. To answer this question we tested whether the difference in the distributions of normalised counts of synonymous and missense sites between the two species was statistically significant by a 2-samples Kolmogornov-Smirnov test for groups of genes characterised by increasing expression rate (5%-percentiles of TPM).

Python scripts used for data handling, parsing, plotting and statistical testing are available in the *FinalPipeline.ipynb* (https://github.com/emitruc/ExpressionLoad).

## 3.2 Polymorphic sites-based analyses

For each polymorphic SNP in our genomic dataset (*daf.joint.no00*), we separately estimated the RNA reads coverage in King and Emperor penguin samples using SAMtools module *depth* and a bed file listing all polymorphic sites in the genomic dataset. Per site RNA read coverage was then added to the genomic dataset, as one value (total number of reads) per species (*daf.joint.no00.rnaCov*; https://doi.org/10.6084/m9.fishshare.23863503.v1). Raw mRNA coverage per site was normalised as counts per million reads (CPM) dividing this value by the sum of mapped RNA reads x 1 million, separately per species:

emp_counts = [62589631 + 85819833 + 72587896] #mapped reads in each of the three pools
king_counts = [66803162 + 64732729 + 58348853] #mapped reads in each of the three pools

emp_CPM = emp_RnaCov / emp_counts * 1000000

king_CPM = king_RnaCov / king_counts * 1000000

After applying the same stringent filters for coverage (from 6X to 8X), outgroup missing data (at least 4 alleles present), and King and Emperor penguin missing data (no missing genotype), we generated different datasets of synonymous and nonsynonymous variants separately per species. Next, after capping the maximum CPM value to 5, we grouped the variants in each dataset applying the same 100 bins of CPM. We then estimated the ratio of nonsynonymous over synonymous variants in each bin separately per species (Supp. Fig. 6). To investigate the effect of gene expression on the allele frequency of synonymous and nonsynonymous variants, we estimated the derived allele frequency spectra grouping the variants by discrete values of CPM (Supp. Fig. 7): i) CPM < 0.3, CPM > 0.3; ii) CPM < 0.5, CPM > 0.5; iii) CPM < 0.3, 0.3 < CPM < 2, CPM > 2. In order to exclude the possibility that a few genes were driving the observed pattern (i.e., pseudoreplication) in the comparison of the site frequency spectra at different expression rates, we replicated ten times the analysis by grouping variants with CPM < 0.3 or CPM > 0.3 after randomly subsampling one synonymous and one nonsynonymous variant per gene. Out of ten replicas, we estimated the 95% intervals for each derived allele count in the site frequency spectra (Figure 3 in main text).
<figure><img src="image_11.png" /><figcaption>**Supplementary Figure 6.** Ratio of missense to synonymous segregating sites per 0.05 intervals of counts per million (CPM) mRNA coverage of each site. Total mRNA read coverage (normalised as counts per million reads) across five tissues from three specimens has been scored for all nonsynonymous and synonymous sites. Note these are density histograms where the total sums up to 1 in each species.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Supplementary Figure 7.** Site frequency spectra of all nonsynonymous (left column) and synonymous variants (right column) across all genes of Emperor (blue) and King (yellow) penguins with mRNA expression partitioned by different values of CPM. The relative frequency of each count class is log10 transformed.</figcaption></figure>
## 4. Testing the relative weight of effective population size and gene expression in realistic forward simulations

The genomic and transcriptomic data analysed in the previous section show that i) purifying selection (estimated as pN/pS) on any gene appears as correlated with its level of expression, and that ii) the effect of gene expression on purifying selection overrides the effect of population size (i.e., King and Emperor penguin populations show similar trends concerning pN/pS while their demographic histories were different; Trucchi et al 2014, Cristofari et al 2016, 2018). To understand the relative contribution of gene expression and population size on purifying selection, we devised a forward in time genomic simulation framework in SLiM v4.0.1 (Haller and Messer 2023). In particular, we used both general Wright-Fisher (WF) and penguin-specific non Wright-Fisher (nonWF) models to study the effect of different population sizes (i.e., genetic drift) and the effect of different selection coefficients on pN/pS. Our hypothesis was that the effect of demography alone isn't strong enough to generate the pN/pS values as observed in our King and Emperor penguin data, but a much stronger selection coefficient on highly expressed genes is needed in the model.

Firstly, to investigate the correlation between population size (Ne) and pN/pS, we designed a WF (WF Pop Size Effect Model) and a nonWF model (nonWF Pop Size Effect Model), which we ran with different carrying capacities (Ne = 1,000, 10,000 and 100,000). The main difference in the nonWF model is that we modelled realistic King penguin life history traits (Céline Le Bohec, personal communication) with population size as a non-fixed parameter which can fluctuate around the set carrying capacity. The genomic model is the same in WF and nonWF simulations and it is implemented as a reduced version of the whole CDS of the King penguin with the coding component of 1,000 genes of length 2,400 bp (i.e., the mean value for the coding sequence per gene in King penguin). The recombination rate is set at 1e-8 within genes and at 4.8e-4 between genes (1e-8 rate scaled to the average intergene length). Mutation rate is set at 1e-8 and the ratio between the occurrence of deleterious and neutral mutations is 2.31:1 (Kim et al 2017). The selection coefficient is assigned to each deleterious mutation using a random value from a gamma distribution with mean -0.01314833 and shape 0.186 (Kim et al 2017). This distribution has been used before in humans and other mammals but we believe it could approximate the distribution of fitness effect also in our target species. For the dominance coefficient we used a *h*-mix model (Kyriazis et al 2021), where weakly deleterious mutations (*s* ≥ -0.01) are partially recessive (*h* = 0.25), while strongly deleterious mutations (*s* < -0.01) are totally recessive (*h* = 0).

We designed similar WF and nonWF models (WF Gene Expression Effect Model and nonWF Gene Expression Effect Model, respectively) to investigate the effect of gene expression on pN/pS. As the pN/pS in the largest simulated population size (Ne = 100,000) of the Pop Size Effect models did not reach small values as observed in our King and Emperor penguin data, we implemented a more extreme selection scenario. In the previous Pop Size Effect models the gamma distribution used to randomly assign the selection coefficient of a novel deleterious mutation resulted in most of the mutations being weakly deleterious (*s* >= -0.01). Here, we assigned a fixed selection coefficient (*s* = -0.001, -0.01, and -0.1) to every gene so that all nonsynonymous mutations appearing in a gene have the same selection coefficient, then we simulated 100 genes for each selection coefficient resulting in a set of 300 genes. We
followed the same strategy for the dominance coefficient assignment, thus it is assigned to every gene deriving it from its fixed selection coefficient following the *hs* relationship (1) (Kyriazis et al 2021), so that dominance and selection coefficients result to be inversely proportional.

$$ (1) \quad h(s) = \frac{1/2}{(1 + 7071.07 \times s)} $$

The recombination coefficient is set at 1e-8 within genes and 0.5 between genes in order to make genes independent to one another. Since most of the genes resulted in a low number of mutations at the end of the simulated generations and to keep the computational running time tractable, we instead increased the gene length from 2,400 to 34,000 bp (the average total gene length in King penguin genome) in the final models, we then estimated the pN/pS per each selection coefficient category of genes under different carrying capacities (Ne = 1,000, 10,000, 100,000).

In total, we designed four models, each of them testing three carrying capacities. For each of these 12 scenarios, we ran three replicas of 10*N* generations each (except for the Ne = 100,000 model where we ran for Ne generations due to computational time), as suggested by SLiM authors in order to make the population reach an equilibrium state (i.e., burn-in). At the end of the simulations, we estimated the pN/pS based on 24 individuals, as in our genomic and transcriptomic real data, repeating the estimate 100 times by randomly resampling 24 individuals. Slim scripts are available at XXX github repository (to be published).

<figure><img src="image_13.png" /><figcaption>**Supplementary Figure 8.** Selection coefficients in highly expressed genes. Distributions of pN/pS estimated in 100 genes where each deleterious mutation is assigned a fixed selection coefficient of -0.001, -0.01, or -0.1 in simulated populations of 1,000, 10,000, or 100,000 individuals, plotted by population size (left panels) or selection coefficient (right panels). As expected, stronger selection coefficients (i.e., -0.1) are effective also when population size is small.</figcaption></figure>
## 5. Expression rate per fitness effect class predicted by SNPeff

Genetic load is the cost paid by any population to its potential further evolution. It is an inherent feature of populations evolving by random mutations which more often have deleterious than advantageous effects on fitness. In conservation genomics, genetic load is getting growing attention as a more appropriate measure of a population's genomic health (Bertorelle et al 2022). However, genetic load is difficult to estimate in non-model species, especially when relying on genomic data only, without information on mutations' fitness effect. In such a case, genetic load can be estimated using *i*) either the predicted effect of a mutation on the amino acid sequence (Cingolani et al 2012), *ii*) or the evolutionary conservation of a certain allele at orthologous sites across multiple species (i.e., GERP; Davydov et al 2010). The latter has the advantage that it can be applied even outside coding regions, but on the other hand, it requires large multi-species genomic alignments, which are extremely computation intensive and error-prone. Even if it can only be applied to coding sequences, a genomic load proxy based on gene expression could be more accurate and easier to compute, given that mRNA expression data from multiple tissues of the target species (or one closely related) are available.

After excluding intergenic variants from the *daf.joint.no00.rnaCov* file with a simple bash command (awk '$14 != "intergenic"' daf.joint.no00.rnaCov > daf.joint.no00.rnaCov.noIntergenic) and applying the same filters for coverage (between 6X and 8X) and missing data in ingroup (no missing) and outgroup samples (at least 4 alleles present for polarisation) as before, we applied an additional Z-normalisation to the CPM for clarity of interpretation only.

**#CPM normalization**

eCPMscal = sum of rnaCov_emp / 1000000

kCPMscal = sum of rnaCov_king / 1000000

eCPM = rnaCov_emp / eCPMscal

kCPM = rnaCov_king / kCPMscal

**#CPM Z-normalisation**

eCPMstdz = eCPM - eCPM_mean / eCPM_std

kCPMstdz = kCPM - kCPM_mean / kCPM_std

Next we selected the variants based on their fitness effect predicted by SNPeff (Cingolani et al 2012): LOW-synonymous (we excluded all of the LOW which are present in introns, mostly as splice-region-variants, as these are not covered by mature mRNA seq data), MODERATE, and HIGH effect.

We then summarised the Z-normalised CPM mean and standard deviation across all private variants (segregating and fixed) or only in private fixed variants in each species per SNP effect (Table 1 in the main text).

As site-specific expression of HIGH deleterious variants (mainly start/stop codons loss/gain and splice acceptor/donor variants) could be biased in mature mRNA sequencing due to
their sequence position, we also estimated their expression using the expression rate of the gene they occurred in. Using the gene-based expression data generated before, we tested whether the expression rate of the group of genes with HIGH deleterious derived alleles was significantly lower than all the rest of the genes. Even if we applied a rather conservative test, the expression of genes with predicted highly deleterious variants is on average three times lower than all genes (Kolmogorov-Smirnov test *p*-value = 0.00095) in the King penguin and slightly lower, even if not significant, in the Emperor penguin.

MODERATE variants show the same average expression rate in both Emperor and King penguin samples (0.78) while a higher average expression in private fixed variants, with a larger difference in Emperor (1.41) than King (0.93) penguin sample. The variance in expression is also larger in fixed than in segregating variants (Table 1 in the main text). Instead of a major role of purging, which could anyway still be part of the process, derived nonsynonymous (which mainly contribute to the MODERATE effect variants) alleles which are fixed at highly expressed genes could actually be truly advantageous variants fixed by positive selection. The lower average expression in King penguin is in line with the expected larger effect of random drift in this population leading to fixation of a larger number of derived nonsynonymous variants (2229 in the King penguin as compared to 1166 in the Emperor penguin), but in genes with lower expression rate and, hence, lower effect on individual fitness. More purging could be again suggested in this case, limiting the fixation of deleterious missense in highly expressed genes. We also suggest that the more intense purging of deleterious variants in the King penguin could be the cause of the lower average expression of synonymous variants which were reduced in more expressed genes by background selection (Table 1 in the main text).

To test our hypothesis that private fixed derived MODERATE variants could actually be truly advantageous, we screened the Emperor and King penguin genome data for selection signatures using Sweepfinder2 (De Giorgio et al 2016) and OmegaPlus (Alachiotis et al 2012) in windows of 10 kb, following the authors' instructions with default settings. Fixed derived alleles are in regions showing higher signatures of selection and lower diversity in both species, and higher differentiation between the two species, hallmarks of highly conserved and lowly recombining genomic regions (Supp. Tab. 3). This has to be considered as a preliminary indication of the potentially positive effect of some of the private fixed MODERATE variants in each species. More targeted analyses are, however, necessary to conclude on the fitness effect of fixed differences with a MODERATE effect on fitness.
<table><caption>**Supplementary Table 3.** Signatures of selection (SF2: Sweepfinder2; OP: Omega+), nucleotide diversity ($\pi$) and differentiation between King (k) and Emperor (e) penguins ($F_{ST}$) in 10 kb genomic regions with fixed differences (Fix) or segregating variants (Seg). All comparisons between the distribution of the statistics across regions with fixed differences and segregating variants for both species are significant (Kolmogorov-Smirnov test $p$-value < 0.05).</caption><thead><tr><th></th><th>eFix</th><th>eSeg</th><th>kFix</th><th>kSeg</th></tr></thead><tbody><tr><td>eSF2</td><td>29.953265</td><td>3.886132</td><td>21.077176</td><td>5.764172</td></tr><tr><td>kSF2</td><td>11.175642</td><td>2.338862</td><td>10.622353</td><td>1.844183</td></tr><tr><td>eOP</td><td>3.166764</td><td>2.367568</td><td>2.759333</td><td>2.46341</td></tr><tr><td>kOP</td><td>3.613646</td><td>2.965366</td><td>3.808706</td><td>2.965772</td></tr><tr><td>eπ</td><td>0.002055</td><td>0.003088</td><td>0.002386</td><td>0.002906</td></tr><tr><td>kπ</td><td>0.001484</td><td>0.00196</td><td>0.001539</td><td>0.001929</td></tr><tr><td>F<sub>ST</sub></td><td>0.644212</td><td>0.509277</td><td>0.612177</td><td>0.521541</td></tr></tbody></table>
Supplementary References

Alachiotis N, Stamatakis A, Pavlidis P (2012). OmegaPlus: a scalable tool for rapid detection of selective sweeps in whole-genome datasets. *Bioinformatics*, 28(17), 2274-2275.

Anders S, Pyl PT, Huber W (2015). HTSeq—a Python framework to work with high-throughput sequencing data. *bioinformatics*, 31(2), 166-169.

Bertorelle G, Raffini F, Bosse M et al (2022). Genetic load: genomic estimates and applications in non-model animals. *Nature Reviews Genetics*, 23(8), 492-503.

Bolger AM, Lohse M, Usadel B (2014). Trimmomatic: a flexible trimmer for Illumina sequence data. *Bioinformatics*, 30(15), 2114-2120.

Cristofari R, Bertorelle G, Ancel A et al (2016). Full circumpolar migration ensures evolutionary unity in the Emperor penguin. *Nature communications*, 7(1), 11842.

Cristofari R, Liu X, Bonadonna F et al (2018). Climate-driven range shifts of the king penguin in a fragmented ecosystem. *Nature Climate Change*, 8(3), 245-251.

Danecek P, Auton A, Abecasis G et al (2011). The variant call format and VCFtools. *Bioinformatics*, 27(15), 2156-2158.

Davydov EV, Goode DL, Sirota M et al (2010). Identifying a high fraction of the human genome to be under selective constraint using GERP++. *PLoS computational biology*, 6(12), e1001025.

DeGiorgio M, Huber CD, Hubisz MJ, Hellmann I, Nielsen R (2016). SweepFinder2: increased sensitivity, robustness and flexibility. *Bioinformatics*, 32(12), 1895-1897.

Dobin A, Davis CA, Schlesinger F et al (2013). STAR: ultrafast universal RNA-seq aligner. *Bioinformatics*, 29(1), 15-21.

Garrison E, Marth G (2012). Haplotype-based variant detection from short-read sequencing. *arXiv preprint arXiv:1207.3907*.

Garrison E, Kronenberg ZN, Dawson ET, Pedersen BS, Prins P (2022). A spectrum of free software tools for processing the VCF variant call format: vcflib, bio-vcf, cyvcf2, hts-nim and slivar. *PLOS Computational Biology*, 18(5), e1009123.

Haller BC, Messer PW (2023). SLiM 4: multispecies eco-evolutionary modeling. *The American Naturalist*, 201(5), E000-E000.

Hernandez RD, Williamson SH, Bustamante CD (2007). Context dependence, ancestral misidentification, and spurious signatures of natural selection. *Molecular biology and evolution*, 24(8), 1792-1800.
Kim BY, Huber CD, Lohmueller KE (2017). Inference of the distribution of selection coefficients for new nonsynonymous mutations using large samples. *Genetics*, *206*(1), 345-361.

Kyriazis CC, Wayne RK, Lohmueller KE (2021). Strongly deleterious mutations are a primary determinant of extinction risk due to inbreeding depression. *Evolution letters*, *5*(1), 33-47.

Li H, Handsaker B, Wysoker A et al (2009). The sequence alignment/map format and SAMtools. *bioinformatics*, *25*(16), 2078-2079.

Li H (2013). Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. *arXiv preprint arXiv:1303.3997*.

Marchi N, Excoffier L (2020). Gene flow as a simple cause for an excess of high-frequency-derived alleles. *Evolutionary applications*, *13*(9), 2254-2263.

Modolo L, Lerat E (2015). UrQt: an efficient software for the Unsupervised Quality trimming of NGS data. *BMC bioinformatics*, *16*, 1-8.

Nursyifa C, Brüniche-Olsen A, Garcia-Erill G, Heller R, Albrechtsen A (2022). Joint identification of sex and sex-linked scaffolds in non-model organisms using low depth sequencing data. *Molecular Ecology Resources*, *22*(2), 458-467.

Trucchi E, Gratton P, Whittington JD et al (2014). King penguin demography since the last glaciation inferred from genome-wide data. *Proceedings of the Royal Society B: Biological Sciences*, *281*(1787), 20140528.

Weir BS, Cockerham CC (1984). Estimating F-statistics for the analysis of population structure. *evolution*, *1358-1370*.

Zheng X, Levine D, Shen J et al (2012). A high-performance computing toolset for relatedness and principal component analysis of SNP data. *Bioinformatics*, *28*(24), 3326-3328.