Maria Del Mar

Dr Alba Soler

Pr Guillaume Michael Diss

Dr Stéphanie Knop

Blandin

An astonishing phenotypic diversity can be observed in natural populations. One of the major goals of modern biology is to unravel the genetic origins of this phenotypic landscape. Gene expression is known to be a main determinant of the relationship between genotypes and phenotypes. In recent decades, several analytical and technical advances have made it possible to study gene expression at every step of the expression process (e.g., transcriptome and proteome) and at very large scales. However, a complete exploration of gene expression across the entire process and at the population scale is still lacking. The goal of this dissertation is to get a more comprehensive view of how each layer of gene expression varies, influences each other, and is related to the natural genetic diversity observed within species. To this end, we analysed the transcriptomes and proteomes of a large natural population of S. cerevisiae (bringing together more than 1,000 individuals) and found unsuspected differences between mRNA and protein abundance regulation. Simultaneously, we studied the gene expression process at three different molecular levels (transcriptome, translatome and proteome) and found that important buffering mechanisms underlie the expression variation between individuals.

Introduction

Comprendre l'origine de l'importante diversité phénotypique observée au sein des populations naturelles est au coeur de la biologie moderne. Plus spécifiquement, la détermination des origines génétiques sous-jacentes aux variations de traits observées chez des individus génétiquement distincts est un prérequis indispensable dans de nombreux domaines tels que la médecine, l'agro-alimentaire ou les sciences environnementales. Les liens entre la diversité génétique et phénotypique, aussi appelés relation génotype-phénotype, sont le résultat d'une multitude de facteurs (environnementaux ou internes à l'organisme) influençant un organisme sur plusieurs échelles (moléculaire, cellulaire, tissulaire, par exemple). Le processus d'expression de gènes est quant à lui l'un des principaux moteurs de la relation génotypephénotype [START_REF] Aguet | Molecular quantitative trait loci[END_REF]Albert and Kruglyak, 2015;Tak and Farnham, 2015). En effet, de nombreux travaux ont montré l'impact important des modifications de l'expression génique sur de nombreuses pathologies humaines (Corbett, 2018;Lee and Young, 2013). Ainsi, comprendre comment varie l'expression des gènes à travers les individus est essentiel dans l'exploration de la relation génotype-phénotype. En dépit de sa nature linéaire (l'ADN est transcrit en ARN qui est traduit en protéine), le processus d'expression des gènes est incroyablement complexe. Chacune de ses étapes (transcription, traduction, dégradation des ARNm et des protéines…) est finement régulée via plusieurs centaines de facteurs cellulaires et l'aspect hiérarchique du processus cache en fait de nombreuses interactions entre ARN et protéines (Buccitelli and Selbach, 2020;Liu et al., 2016). Aussi, les mécanismes de variations de l'expression des gènes à l'échelle des populations sont à ce jour grandement incompris.

Au cours des deux dernières décennies, de nombreux progrès techniques et analytiques ont débouché sur une large gamme d'outils permettant à la fois des mesures fines et à grande échelle de l'expression des gènes, mais aussi de relier les variations de cette expression au fond génétique des individus. Sur le plan technique il est possible de citer la mise en place du RNAsequencing (Stark et al., 2019) pour quantifier le transcriptome des individus, le développement à plus large échelle de la spectrométrie de masse en tandem pour leur protéome (Demichev et al., 2020;Messner et al., 2023Messner et al., , 2022)), et enfin la création de techniques plus spécifiques comme le ribosome-profiling permettant de mesurer avec précision le processus de traduction (Ingolia et al., 2019). Au niveau analytique, l'extension à grande échelle des techniques précédemment citées a été accompagnée par le développement des analyses d'associations pangénomiques (Dehghan, 2018) (Genome-Wide Association Studies -GWAS). En utilisant de large cohortes d'individus, ces dernières permettent d'associer des variants génétiques, généralement des variants nucléotidiques, aux variations d'un phénotype précis quantifiées dans une population.

Les variants génétiques ainsi détectés sont appelés QTL (pour Quantitative Trait Loci) et dans le cas de GWAS visant à étudier les variations des niveaux d'ARNm ou de protéines, on parle respectivement d'eQTL (expression QTL) et de pQTL (protein QTL) [START_REF] Aguet | Molecular quantitative trait loci[END_REF].

Cependant et malgré ces différentes avancées, de nombreux aspects concernant les variations d'expression de gènes restent méconnus, notamment à l'échelle des populations. Par exemple, il n'y a toujours pas de consensus sur le degré de similarité entre les variations des transcriptomes et des protéomes à travers de larges groupes d'individus (Buccitelli and Selbach, 2020;Fortelny et al., 2017;Liu et al., 2016). De plus, la similitude entre les eQTL et les pQTL est encore grandement débattue car les résultats des études à ce sujet sont contradictoires.

C'est dans ce contexte que prend place mon projet de thèse. Il s'agit d'explorer simultanément plusieurs étapes de l'expression des gènes en utilisant une population naturelle de la levure S. cerevisiae pour laquelle le génome a été entièrement séquencé via la technique Illumina (Peter et al., 2018) et où un large jeu de données lié au transcriptome a déjà été généré (Caudal et al., 2023). Cependant, l'exploration à grande échelle de l'expression génique est encore difficile à cause de certaines limitations techniques, notamment au niveau protéomique où il est difficile de quantifier avec précision un grand nombre de gène lorsqu'on travaille sur de nombreux individus. De ce fait, ma thèse s'est articulée autour de deux principaux projets : l'un portant sur la quantification précise de l'expression des gènes à travers chacunes des étapes de cette dernière et donc quantifiant un grand nombre de gènes à travers un nombre limité d'individus (8 isolats naturels de S. cerevisiae), l'autre portant sur l'exploration à travers un très grand nombre d'individus du transcriptome et du protéome de S. cerevisiae mais en ayant une couverture génique plus faible.

Resultats

La variation de traduction à travers différents fonds génétiques révèle une signature de l'atténuation post-transcriptionnelle chez la levure Dans ce projet, nous avons quantifié avec précision l'expression de plus de 4 344 gènes dans 8 isolats naturels de S. cerevisiae provenant d'environnements très différents. La quantification de l'expression de gènes s'est faite via du RNA-sequencing, du ribosome-profiling et de la spectrométrie de masse en tandem, permettant d'avoir une vision globale des dynamiques de régulation tout au long du processus d'expression.

Plusieurs points ont été explorés au cours de ce projet. Tout d'abord, en comparant les données de transcription et de traduction (donc, de RNA-sequencing et de ribosome-profiling fait en collaboration avec le Riken Institute -Tokyo), nous avons observé que les variations d'expression de gènes semblaient principalement liées aux préférences trophiques de chacune des souches de levures. Celles-ci ayant des origines écologiques très variées, elles ont probablement adapté la régulation des gènes liés au métabolisme en conséquence. Si ces résultats sont partagés entre la régulation transcriptionnelle et traductionnelle, ces dernières sont pourtant assez différentes en termes d'intensité. En effet, nous avons observé que les variations d'abondance d'ARNm sont plus importantes que les variations observées sur les données de ribosome-profiling. Ceci est dû à un phénomène appelé atténuation posttranscriptionnelle (post-transcriptional buffering) déjà décrit dans des études antérieures (Artieri and Fraser, 2014;Blevins et al., 2019;McManus et al., 2014;Wang et al., 2020). Ce phénomène suggère que les étapes avancées de l'expression de gènes sont plus conservées et donc plus contraintes évolutivement parlant. Nous avons aussi détecté que l'atténuation posttranscriptionnelle affecte préférentiellement certains types de gènes, comme les gènes essentiels, les gènes liés à des complexes protéiques ou même les gènes ayant un faible niveau d'expression. Ceci n'avait pas été démontré jusqu'à présent et permet d'éclaircir les mécanismes sous-jacents à l'atténuation post-transcriptionnelle qui sont encore très méconnus.

Enfin nous avons utilisé les données d'abondances transcriptionnelles et traductionnelles pour explorer comment certains gènes présents chez S. cerevisiae mais issus d'espèces différentes de levures sont exprimés. Ces gènes sont issus de mécanismes d'échange de matériel génétiques entre espèces (comme les introgressions ou les transferts horizontaux de gènes) et la régulation de leur expression a rarement été explorée [START_REF] D'angiolo | A yeast living ancestor reveals the origin of genomic introgressions[END_REF]Marsit et al., 2015;Novo et al., 2009;Peter et al., 2018). Nous avons pu observer différents profils d'expression en fonction de l'origine de ces gènes. Les introgressions avait par exemple des niveaux d'expression similaires en comparaison avec leur orthologues alors que les gènes issus de transferts horizontaux sont pour leur part moins traduit que les autres gènes de S. cerevisiae.

Nous avons ensuite étendu l'exploration de l'expression de gènes de ces 8 souches au niveau protéique (en collaboration avec le Weizmann Institute of Science -Israël). Plusieurs résultats précédemment obtenus sont retrouvés au niveau protéique. Tout d'abord les variations d'abondance protéique sont principalement liées aux gènes du métabolisme, ce qui confirme les résultats obtenus en utilisant des données de RNA-sequencing et de Ribosome-profiling.

De plus, en comparant le niveau de variation au sein de chacune des étapes d'expression de gènes, nous avons aussi pu observer le phénomène d'atténuation post-transcriptionnelle, les variations d'abondance protéique étant les plus faibles, et donc cette étape du processus d'expression est de fait la plus conservée à travers les 8 isolats. Ces résultats confirment la présence du phénomène d'atténuation post-transcriptionnelle. Nous nous sommes par la suite penchés plus précisément sur les différences de contraintes évolutives entre chacune des étapes du processus d'expression. Dans cette optique, nous avons utilisé une base de données décrivant plusieurs centaines de caractéristiques des gènes de S. cerevisiae. À l'aide d'une quantification précise de la vitesse évolutive associée à chacune de ces caractéristiques (Wang et al., 2020), nous avons pu montrer que l'évolution de la régulation de l'abondance génique, bien qu'ayant de nombreuses spécificités pour chacune des étapes de l'expression, suit plusieurs principes communs, notamment que les gènes ayant une place centrale dans les interactions entre protéines ou ayant un rôle fondamental dans le fonctionnement cellulaire verront leur expression grandement conservée à travers les individus. À l'inverse et en adéquation avec les résultats précédemment cités, les gènes liés au métabolisme et aux capacités respiratoires des cellules ont une expression qui évolue beaucoup plus rapidement que les reste des gènes.

Les transcriptomes et protéomes quantitatifs à l'échelle de l'espèce révèlent un contrôle génétique distinct de la variation de l'expression génique chez la levure

Conclusion

L'ensemble des travaux réalisés pendant ma thèse ont permis d'éclaircir les mécanismes de régulations qui façonnent les variations d'expression de gènes à l'échelle d'une population. Le message principal étant que chacune des étapes du processus d'expression suit des dynamiques de régulations propres, avec une tendance globale à la diminution de la variation d'expression observée entre chaque individu au fur et à mesure de l'expression génique. Dans le cadre de la compréhension de la relation génotype-phénotype, il est donc primordial de considérer l'impact d'un variant tout au cours du processus d'expression génique pour comprendre comment il peut impacter la diversité des traits observés au sein d'une population naturelle. exploration ................................................................ regulation ........................................................................... Submitted .................................................................................................................................................. In preparation .......................................................................................................................................... Oral ............................................................................................................................................................ Poster .......................................................................................................................................................

Résumé graphique du projet

A more exhaustive transcriptome and proteome

Population-wide exploration of translation

STATE OF THE ART

The genotype-phenotype relationship: it's complicated… An exceptional phenotypic diversity can be observed within all species. On both macroscopic and molecular scales, individuals differ in a myriad of observable and quantifiable traits. The origin of this diversity has long been questioned and understanding the mechanisms underlying the phenotypic landscape observed in a population has been a central and long-standing challenge in biology. In fact, as humanity moved from nomadism to sedentism, relying on both agriculture and livestock, it became crucial to improve the traits of interest in domesticated plants and animals. Thus, the control of phenotypic diversity quickly became a central concern.

A complex relationship for complex traits

A long-standing interest Throughout the history of science, the relationship between heredity and traits has been questioned. Traces of primitive explorations of pedigree have been observed on a 6000 year old Babylonian tablet, which may describe horse breeding, suggesting early hypotheses on the hereditary nature of livestock characteristics (Coonen, 1952). Later, during the Ancient Greek period (400-300 BCE), physicians and philosophers also established early theories of trait heredity and reproduction [START_REF] Bazopoulou-Kyrkanidou | Genetic concepts in Greek literature from the eighth to the fourth century B[END_REF]. However, in the late 19 th century, Mendel established the first law of heredity [START_REF] Mendel | Versuche über Pflanzen-Hybriden[END_REF] and, with the rediscovery of his work in the early 20 th century by De Vries, Correns and Tschermak, genetics emerged as one of the major disciplines in the life sciences. At the same time, several fundamental concepts in genetics were proposed and established by Bateson and Johannsen, such as "gene", "genotype" and "phenotype" [START_REF] Bateson | Mendel's principles of heredity[END_REF][START_REF] Johannsen | The Genotype Conception of Heredity[END_REF]. These terms respectively described the "unit factor" of heredity, the "sum of all the genes" and the "direct inspection […] or direct measures of assessment" [START_REF] Johannsen | The Genotype Conception of Heredity[END_REF]. Johannsen supported the theory of a direct relationship between genotype and phenotype. The highly valuable implications of a precise dissection of the genotype-phenotype relationship in medicine, agriculture, food production, or industry quickly led to intensive research efforts. At the same time, several fundamental discoveries in the 20 th century led to a clearer view of the genetic mechanisms underlying heredity. One example is the work of Morgan, Sturtevant, Muller and Bridges on the fruit fly, which combined Mendelian and chromosomal theories [START_REF] Morgan | Poor transcript-protein correlation in the brain: negatively correlating gene products reveal neuronal polarity as a potential cause[END_REF]. This work was a real keystone at the time, on which several other major discoveries were based, such as the concept of mutation [START_REF] Muller | The Production of Mutations by X-Rays[END_REF]. Later, in the mid-20 th century, the rise of molecular biology and biochemistry also revolutionized genetics. For example, Beadle and Tatum demonstrated that genes are involved in biochemical reactions [START_REF] Beadle | Genetic Control of Biochemical Reactions in Neurospora[END_REF] through the action of an enzyme. In addition, Rosalind Franklin's exploration of the structure of the deoxyribonucleic acid (DNA) molecule and the results published by Crick and Watson in the 1950s [START_REF] Watson | The Structure of Dna[END_REF] clearly established the role of DNA as being the vector of heredity. At the end of this period, a major discovery introduced the concept of gene expression by messenger ribonucleic acid (mRNA) molecules and the regulation of gene expression [START_REF] Jacob | Genetic regulatory mechanisms in the synthesis of proteins[END_REF]. From this point on, the relationship between the amazing diversity observed in natural populations and its biological origin became clearer. Genetic variations in each individual, through the process of gene expression (transcription and translation), affect the function or quantity of a given protein, which ultimately affects one or more traits of the individual. This link between genotype and phenotype is commonly referred to as the genotypephenotype relationship and is still the subject of considerable research in modern biological science.

The decomposition of a phenotype

Despite its seemingly simple nature, the genotype-phenotype relationship is a truly complex and subtle process. First, the number of genes influencing a trait can range from 1 to several thousand. The distribution of the phenotypes within a population according is usually a good approach to determine if the trait is control by one or several genes (Figure 1A,B). If the population exhibits a phenotypic bimodal distribution (Figure 1A), this usually indicates that the phenotype is controlled by a single gene. In humans, about 6,000 thousand diseases [START_REF] Condò | Rare Monogenic Diseases: Molecular Pathophysiology and Novel Therapies[END_REF], such as cystic fibrosis, neurofibromatosis or Duchenne muscular dystrophy, are caused by a single defective gene. Several non-pathological characteristics are also determined by a single gene, such as the ABO blood group [START_REF] Yamamoto | Molecular genetic basis of the histoblood group ABO system[END_REF]. Conversely, and for most quantifiable traits, the phenotypic distribution within a population follows a normal distribution (Figure 1B) These traits are considered as complex. Autism, Alzheimer's disease, or human height [START_REF] Akiyama | Characterizing rare and low-frequency height-associated variants in the Japanese population[END_REF][START_REF] Nikolac Perkovic | Genetic Markers of Alzheimer's Disease[END_REF][START_REF] Ramaswami | Genetics of autism spectrum disorder[END_REF] are notable examples of complex traits for which hundreds or even thousands of genetic variants have been found to influence the appearance or intensity of the trait. Deciphering and capturing all the genetic factors involved in such traits is a tedious, but crucial, task for predicting and treating complex diseases. It is worth noting that both monogenic and complex traits are, of course, not solely due to genetic factors. An individual's environmental background can influence a trait as much, if not more, than their genetic background. In the case of autism, for example, several external factors such as zinc deficiency, prenatal and perinatal stress, or parental age are known to be associated with the onset of autism spectrum disorders [START_REF] Grabrucker | Environmental Factors in Autism[END_REF]. Another and perhaps more famous example of the impact of environmental factors on a phenotype is the association between lung cancer and smoking [START_REF] Hecht | Cigarette smoking: cancer risks, carcinogens, and mechanisms[END_REF].

σ 2 P = σ 2 G + σ 2 E + σ 2 GxE A C B σ 2 Additivity + σ 2 Dominance + σ 2 Epistasis X X X H 2 h 2
The phenotypic variance of a trait (s 2 P) is therefore considered as equal to the sum of the genetic variance (s 2 G), the environmental variance (s 2 E) as well as the genetic and environmental interaction variance (s 2 GxE) (Figure 1C). The genetic variance itself is a highly complex composition of additive and non-additive effects (Figure 1C). The fraction of phenotypic variance controlled by the genetic variance is called broad-sense heritability (H 2 ) while the narrow-sense heritability (h 2 ) includes only the fraction of the phenotypic variance associated with additive genetic effects (Figure 1C). The non-additive effects are related to phenomena such as epistasis and dominance [START_REF] Su | Estimating Additive and Non-Additive Genetic Variances and Predicting Genetic Merits Using Genome-Wide Dense Single Nucleotide Polymorphism Markers[END_REF]. Briefly, epistasis describes the interaction between two or more loci that causes a phenotype different from that expected from an additive effect, while dominance describes the situation where one allele at a heterozygous locus masks the phenotype attributed to the alternative allele. Recent estimates in the budding yeast Saccharomyces cerevisiae showed that trait variation is mostly due to additive effects, about 55% of phenotypic variation, while non-additive effects account for 29% of phenotypic variation (Bloom et al., 2015(Bloom et al., , 2013;;Fournier et al., 2019).

From molecular changes to macroscopic traits

The relationship between a specific DNA sequence modification and a phenotypic variation depends on multiple mechanisms that affect individuals at several scales. For example, a mutation will influence molecular reactions that affect one or more cellular processes, resulting in tissue, organ, or overall macroscopic trait variation. Studying the molecular mechanisms underlying phenotypic changes is therefore a critical step to fully understand the genotypephenotype relationship.

Genetic diversity

Identifying and characterizing the genetic variants that cause a change in phenotype is the first step in understanding the molecular processes involved. The simplest and most common type of DNA variant is the modification of a single base (Figure 2A), often called a Single Nucleotide Polymorphism (SNP). Extensive exploration of the genomes of more than 2,500 human individuals has for example resulted in the discovery of more than 88 million SNPs (The 1000 Genomes Project Consortium, 2015). SNPs are often used to describe the genetic differences within species. In fact, the number of SNPs in each individual (compared to a reference genome), and its integration at a population-scale level is a good, although incomplete, indicator of the genetic diversity within this population. In humans, each individual carries on average 4.5 million SNPs, which corresponds to one SNP every 1,000 bp. In microorganisms, intraspecific genetic diversity can be much higher. For example, in the budding yeast S. cerevisiae, an average of 1 SNP per 200 bp has been reported between individuals from a natural population (Peter et al., 2018). The second most common type of genetic variant is insertion or deletion of a few bases, commonly referred to as "indels" (Figure 2A). More than 3.6 million indels have been identified in the 2,500 human genomes (The 1000 Genomes Project Consortium, 2015). Both SNPs and indels can be detected efficiently using short-read sequencing techniques such as Illumina sequencing.

Structural variants (SVs) comprise diverse type of genetic changes that are much larger (at least >50 base pairs) than SNPs and indels. SVs include large chromosomal alterations such as deletions or insertions, translocations, aneuploidies, copy number variants (CNV), inversions and duplications (Figure 2B). The SVs are usually much more difficult to detect than SNPs or indels. Although the detection of these variants is theoretically possible using classical shortread sequencing methods (such as Illumina paired-end sequencing), their detection is more reliable using specific sequencing techniques, such as long-read sequencing (Logsdon et al., 2020;[START_REF] Shi | Long-read sequencing and de novo assembly of a Chinese genome[END_REF][START_REF] Wenger | Accurate circular consensus long-read sequencing improves variant detection and assembly of a human genome[END_REF]. Although they are less common, about 60,000 SVs have been detected in the 2,500 human genomes, they affect a larger number of bases compared to SNPs and indels. On average, an individual will differ from the human reference genome at about 4.5 million genetic positions due to SNPs and indels, while the individual's SVs will affect an average of 20 million bases (The 1000 Genomes Project Consortium, 2015). It is important to note that until recently, the vast majority of large-scale genomic explorations were based on short-read sequencing methods, resulting in poor exploration and characterization of SVs across individuals. Today, the development and democratization of long-read sequencing technologies, such as Oxford Nanopore Sequencing or PacBio technologies, are filling this gap [START_REF] Audano | Characterizing the Major Structural Variant Alleles of the Human Genome[END_REF]He et al., 2023;Zhou et al., 2022). Overall, the nature of a variant will, of course, be a major determinant of its molecular effects, as well as several other characteristics such as its genetic location, its homo-or heterozygous state, or its size. The frequency at which a genetic variant occurs in a population, often measured as minor allele frequency (i.e., the percentage of individuals carrying the 2 nd most common allele), is often low (Figure 3). Overall, the vast majority of variants present in an individual are rare (Peter et al., 2018;The 1000Genomes Project Consortium, 2015). Interestingly, rare variants are thought to play an important role in complex traits, even if detecting their actual contribution to the phenotype can be laborious [START_REF] Bodmer | Common and rare variants in multifactorial susceptibility to common diseases[END_REF][START_REF] Cirulli | Uncovering the roles of rare variants in common disease through wholegenome sequencing[END_REF][START_REF] Gibson | Rare and common variants: twenty arguments[END_REF]The UK10K Consortium, 2015). In yeast, for example, rare variants have been shown to play a major role in several growth phenotypes [START_REF] Bloom | Rare variants contribute disproportionately to quantitative trait variation in yeast[END_REF]Fournier et al., 2019) and deleterious variants tend to be enriched among the rare variants (Figure 3) (Peter et al., 2018). The rare variants (MAF<0.05) account for 92% of the SNPs in a natural population of S. cerevisiae.

Among them, deleterious variants tend to be enriched (adapted from Peter et al., 2018).

DNA and molecular changes for macroscopic consequences

The molecular changes induced by genetic variants are diverse and include a wide range of phenotypes. For example, a change in the amino acid composition of a protein or the apparition of a premature stop codon is a possible consequence of non-synonymous substitutions that can be induced by SNPs (Figure 2A). Even though non-synonymous SNPs are less common than synonymous SNPs, they are the most frequently associated with human diseases [START_REF] Yip | Annotating single amino acid polymorphisms in the UniProt/Swiss-Prot knowledgebase[END_REF]. Such a change is obviously a strong modifier of several protein properties. The enzymatic function itself can be affected as it is the case for the human gene ALDH1L1, a tumor suppressor gene involved in folate metabolism [START_REF] Krupenko | The Role of Single-Nucleotide Polymorphisms in the Function of Candidate Tumor Suppressor ALDH1L1[END_REF] whose catalytic activity is reduced in individuals with a specific SNP causing an amino acid change [START_REF] Frosst | A candidate genetic risk factor for vascular disease: a common mutation in methylenetetrahydrofolate reductase[END_REF] leading to an increased risk of cardiovascular disease, birth defects, and cancer.

Protein stability is also a known molecular property that can be altered by SNPs or indels [START_REF] Bromberg | Correlating protein function and stability through the analysis of single amino acid substitutions[END_REF][START_REF] Casadio | Correlating disease-related mutations to their effect on protein stability: A large-scale analysis of the human proteome[END_REF]. Severe diseases such as amyotrophic lateral sclerosis [START_REF] Lindberg | Systematically perturbed folding patterns of amyotrophic lateral sclerosis (ALS)-associated SOD1 mutants[END_REF][START_REF] Ling | ALS-associated mutations in TDP-43 increase its stability and promote TDP-43 complexes with FUS/TLS[END_REF] are associated with an increased or decreased protein stability [START_REF] Randles | Using Model Proteins to Quantify the Effects of Pathogenic Mutations in Ig-like Proteins[END_REF][START_REF] Wang | SNPs, protein structure, and disease[END_REF]. Finally, modification of protein amino acid sequence can also affect the protein-protein interaction ability of a particular protein [START_REF] Cheng | Comprehensive characterization of protein-protein interactions perturbed by disease mutations[END_REF][START_REF] Porta-Pardo | A Pan-Cancer Catalogue of Cancer Driver Protein Interaction Interfaces[END_REF]. Sickle-cell disease, an autosomal recessive pathology, is a famous case of protein interaction modification where a unique SNP that is located in the coding sequence of the β-globin gene [START_REF] Ingram | Gene Mutations in Human Haemoglobin: the Chemical Difference Between Normal and Sickle Cell Haemoglobin[END_REF][START_REF] Rees | Sickle-cell disease[END_REF][START_REF] Sundd | Pathophysiology of Sickle Cell Disease[END_REF] leads to an abnormal protein aggregation into large polymers which ultimately causes an abnormal red cell shape [START_REF] Bunn | Pathogenesis and Treatment of Sickle Cell Disease[END_REF]. In addition, it is worth noting that synonymous SNP (i.e., a nucleotide substitution that does not change the final amino acid sequence of a protein, figure 2A), despite their apparently neutral effect, are also known to be associated with human disease [START_REF] Sauna | Synonymous Mutations as a Cause of Human Genetic Disease[END_REF]. An elegant example of this is in cystic fibrosis, where a synonymous SNP associated with the disease alters the mRNA structure of the CFTR gene, resulting in misfolding of the cognate protein and its dysfunction [START_REF] Bartoszewski | A Synonymous Single Nucleotide Polymorphism in ΔF508 CFTR Alters the Secondary Structure of the mRNA and the Expression of the Mutant Protein[END_REF]. Although less studied than SNPs, SVs can also induce specific molecular changes. A very recent example of this is the case of rust resistance in wheat, which has been linked to a translocation linked to an introgression event (the transfer of genetic material from one species to another through hybridization), and results in a fused kinase that drives resistance to the pathology (Wang et al., 2023).

The case of missing heritability

In the recent decades, many studies have attempted to unambiguously identify the genetic variants that influence the onset or intensity of human diseases or phenotypes by testing the association of each variant (usually SNPs) with the phenotypic values quantified in a population. This is usually done using analytical tools such as genome-wide association studies or linkage mapping studies, which will be described later. This has become particularly relevant with the development of short-read sequencing techniques (e.g., Illumina sequencing) since 2005, allowing for more and more large-scale variant characterizations. For example, the influence of genetic factors on human height has been extensively studied (Wood et al., 2014;[START_REF] Yang | Genetic variance estimation with imputed variants finds negligible missing heritability for human height and body mass index[END_REF][START_REF] Yang | Common SNPs explain a large proportion of the heritability for human height[END_REF][START_REF] Yengo | Meta-analysis of genome-wide association studies for height and body mass index in ∼700000 individuals of European ancestry[END_REF] and the most recent large-scale study using data from 5.4 million individuals identified more than 12,000 SNPs significantly associated with height (Yengo et al., 2022). However, the fraction of phenotypic variation explained by this large number of genetic variants is at most 40% (depending on the ancestry of the individuals), suggesting that the dissection of the origin of height variation is far from complete.

Interestingly, for many complex human traits, examination of the genetic fraction of phenotypic variation fails to explain most of the trait variation. For example, in the case of familial colorectal cancer, less than half of the heritability of this disease is associated with clearly identified genetic variants [START_REF] Schubert | The missing heritability of familial colorectal cancer[END_REF]. Another example is the case of Alzheimer's disease, where the various identified genetic variants account for only 7.78% of the phenotypic variance [START_REF] Ridge | Alzheimer's Disease: Analyzing the Missing Heritability[END_REF], while the estimated heritability of this disease reaches 58% [START_REF] Gatz | Role of Genes and Environments for Explaining Alzheimer Disease[END_REF]. This phenomenon is better known as the "missing heritability" and has been widely discussed and investigated (Génin, 2020;[START_REF] Maher | Personal genomes: The case of the missing heritability[END_REF]Manolio et al., 2009;[START_REF] Owen | Explaining the missing heritability of psychiatric disorders[END_REF][START_REF] Young | Solving the missing heritability problem[END_REF]. Several hypotheses have been formulated to explain the gap between the total estimated heritability and the proportion experimentally associated with genetic variants. Rare genetic variants, which are poorly considered when studying the relationship between a trait and its genetic origin, epigenetic factors as well as genetic interactions could explain part of this missing heritability. Moreover, since SVs have so far been poorly characterized across individuals, their effects on phenotypes have been less studied, which may also explain some of the missing heritability. This large gray area is still far from being resolved, and it highlights how complex and still poorly understood the genotypephenotype relationship is.

Gene expression as a driver of the genotype-phenotype relationship

Interestingly, despite the strong influence of SNPs in coding regions at the molecular level, the majority of SNPs detected as influencing a specific trait in humans are located in non-coding or intronic regions [START_REF] Aguet | Molecular quantitative trait loci[END_REF]Tak and Farnham, 2015). This suggests that the effect of these non-coding SNPs is more likely to be on regulatory processes, causing changes in gene expression, which ultimately induce cellular and macroscopic phenotypic changes.

All molecular steps of gene expression can be affected by a genetic variant. The initial accessibility of a gene to the transcriptional machinery is strongly influenced by nearby or distal DNA modifications. Chromatin organization, for example, is a tightly regulated process that is highly sensitive to SNPs or indels [START_REF] Degner | DNase I sensitivity QTLs are a major determinant of human expression variation[END_REF][START_REF] Delaneau | Chromatin three-dimensional interactions mediate genetic effects on gene expression[END_REF]. Similarly, DNA methylation alteration, which has long been associated with cancer [START_REF] Das | DNA Methylation and Cancer[END_REF][START_REF] Koch | Analysis of DNA methylation in cancer: location revisited[END_REF], is also a known regulatory step that is tightly regulated and influenced by numerous genetic variants [START_REF] Hawe | Genetic variation influencing DNA methylation provides insights into molecular mechanisms regulating genomic function[END_REF]Min et al., 2021). Overall, several diseases are linked to transcriptional modification, as several types of cancer are known to be associated with aberrant transcriptional regulation, such as prostate cancer [START_REF] Demichelis | Identification of functionally active, low frequency copy number variants at 15q21.3 and 12q21.31 associated with prostate cancer risk[END_REF], melanoma [START_REF] Huang | Highly recurrent TERT promoter mutations in human melanoma[END_REF] or leukemia [START_REF] Sanda | Core transcriptional regulatory circuit controlled by the TAL1 complex in human T cell acute lymphoblastic leukemia[END_REF]. Dysregulation of mRNA degradation is also a known factor influencing human disease [START_REF] Saramago | The Implication of mRNA Degradation Disorders on Human DISease: Focus on DIS3 and DIS3-Like Enzymes[END_REF]. Furthermore, targeting specific mRNA degradation pathways is a promising avenue for the development of anticancer therapies [START_REF] Bokhari | Targeting nonsense-mediated mRNA decay in colorectal cancers with microsatellite instability[END_REF][START_REF] Huang | Antisense suppression of the nonsense mediated decay factor Upf3b as a potential treatment for diseases caused by nonsense mutations[END_REF][START_REF] Lindeboom | The impact of nonsense-mediated mRNA decay on genetic disease, gene editing and cancer immunotherapy[END_REF]. The molecular effects of genetic variants may also involve post-transcriptional mechanisms. For example, RNA splicing has been shown to be an important link between DNA variation and disease (Y. I. Li et al., 2016). Finally, alterations in translational regulation itself can be a major source of phenotypic variation. A good example in humans is the Fragile X Mental Retardation Protein (FMRP), encoded by the FRM1 gene [START_REF] Verkerk | Identification of a gene (FMR-1) containing a CGG repeat coincident with a breakpoint cluster region exhibiting length variation in fragile X syndrome[END_REF]. This protein normally regulates and represses translation through various mechanisms [START_REF] Li | The fragile X mental retardation protein inhibits translation via interacting with mRNA[END_REF][START_REF] Richter | Dysregulation and restoration of translational homeostasis in fragile X syndrome[END_REF].

When FMRP is not expressed, this leads to a global and abnormal translation of many different mRNA [START_REF] Udagawa | Genetic and acute CPEB1 depletion ameliorate fragile X pathophysiology[END_REF], resulting in the Fragile X syndrome, a common inherited form of intellectual disability (Corbett, 2018).

Therefore, because of the central role of gene expression in the genotype-phenotype relationship, this link between the genetic origins of a trait and its establishment is inherently tedious to dissect and understand. The complete expression process of each gene (i.e., transcription, translation, transcript or protein maturation, and degradation) is the combination of tenths of tightly regulated mechanisms, with a plethora of interactions and retro-controls between each step (Figure 4). The multifaceted nature of the genotype-phenotype relationship implies that the consequences of a genetic variant, and thus the biological mechanism underlying a particular phenotype, may affect any of the gene expression steps. Gene expression is a tightly regulated process in which each step interacts with the others, resulting in a truly complex regulatory network (adapted from Buccitelli and Selbach, 2020)

Determining the role of gene expression in the genotype-phenotype relationship

As described above, alteration of gene expression is one of the main mechanisms that translates genotypes into phenotypes. The accurate study of gene expression at any level (transcriptome, proteome...) is a constant challenge for modern biology, as is the development of bioinformatic tools to link gene expression and genetic variants. Combined, these two aspects of gene expression studies are crucial to dissect the impact of changes in transcript and protein abundance on the phenotypic landscape of natural populations.

Quantification of gene expression

Together with sequencing techniques, the precise quantification of gene expression at each step of the expression process has been a keystone of modern biology, leading to considerable advances in medicine, agriculture, and biotechnology. Over the last two decades in particular, several tools have allowed a steady increase in the precision and scale of quantification of either mRNA or protein abundance.

Quantification of mRNA

Methods for quantifying mRNA molecules in an individual have been developed since the early days of molecular biology. Northern blotting was developed in the late 1970s and at the time was a robust technique for relative quantification of mRNAs of interest [START_REF] Reue | mRNA Quantitation Techniques: Considerations for Experimental Design and Application[END_REF][START_REF] Sambrook | Molecular cloning: a laboratory manual[END_REF]. The principle was based on the transfer of fractionated and separated (by denaturing gel electrophoresis) mRNAs onto a membrane. The mRNA of interest was then revealed on the membrane by a hybridization step using a cDNA probe, either radiolabeled or linked to an enzyme, and relatively quantified by comparing the label intensity to a control label. Also based on hybridization, the ribonuclease protection assay (RPA) involved a liquid mixture between the total mRNAs and a specific probe [START_REF] Azrolan | A solution hybridization/RNase protection assay with riboprobes to determine absolute levels of apoB, A-I, and E mRNA in human hepatoma cell lines[END_REF]. Once hybridization is achieved, single-stranded RNAs are degraded, leaving only the hybrid, which is then electrophoresed on a denaturing gel, allowing relative or even absolute quantification of mRNA using titration reactions [START_REF] Reue | mRNA Quantitation Techniques: Considerations for Experimental Design and Application[END_REF]. However, the precise quantification of mRNA in an absolute manner has been improved and achieved mainly by polymerase chain reaction (PCR) techniques, in particular the combination of reverse transcriptase (RT) and PCR, known as quantitative RT-PCR [START_REF] Bustin | Absolute quantification of mRNA using real-time reverse transcription polymerase chain reaction assays[END_REF][START_REF] Foley | Quantitation of RNA using the polymerase chain reaction[END_REF]. Theoretically, RT-PCR can quantify a single RNA molecule in a sample. However, the practical limit is closer to ten molecules due to RT inefficiency [START_REF] Reue | mRNA Quantitation Techniques: Considerations for Experimental Design and Application[END_REF]. Absolute quantification using quantitative RT-PCR was developed in the late 1980s by adding an exogenous transcript standard [START_REF] Becker-André | Absolute mRNA quantification using the polymerase chain reaction (PCR). A novel approach by a PCR aided transcipt titration assay (PATTY)[END_REF][START_REF] Gilliland | Analysis of cytokine mRNA and DNA: detection and quantitation by competitive polymerase chain reaction[END_REF][START_REF] Wang | Quantitation of mRNA by the polymerase chain reaction[END_REF]. Despite the significant advances in mRNA quantification achieved with these techniques, one of their major drawbacks was their limited suitability for probe multiplexing. In other words, the simultaneous quantification of different mRNA was severely limited. This limitation was overcome with the development of DNA microarrays in the early 1990s [START_REF] Bumgarner | DNA microarrays: Types, Applications and their future[END_REF]. Briefly, the principle is based on the detection of hybridization of DNA fragments on a surface containing probes corresponding to genic regions in the case of mRNA quantification. The hybridization is detected and quantified by the prior labeling of DNA fragments extracted from the sample or obtained after reverse transcriptase of mRNA. This method allows relative mRNA quantification in a much higher throughput than previous techniques and has therefore been widely used to analyze gene expression [START_REF] Schena | Quantitative Monitoring of Gene Expression Patterns with a Complementary DNA Microarray[END_REF][START_REF] Tarca | Analysis of microarray experiments of gene expression profiling[END_REF].

However, the development and the increased accessibility of RNA-sequencing (or RNA-seq, figure 5) from the late 2000s [START_REF] Emrich | Gene discovery and annotation using LCM-454 transcriptome sequencing[END_REF][START_REF] Lister | Highly integrated single-base resolution maps of the epigenome in Arabidopsis[END_REF][START_REF] Nagalakshmi | The transcriptional landscape of the yeast genome defined by RNA sequencing[END_REF] has led to a decline of the use of DNA microarrays. RNA-seq, that is originally based on shortread sequencing (mainly Illumina sequencing), allows a theoretical absolute quantification of the cell's transcripts, giving a global view of all mRNAs from a sample, a tissue or an individual (such quantifications are called transcriptome) and is therefore less biased than microarray methods (Stark et al., 2019). In addition, where microarrays required prior knowledge of the genome of the species being studied to construct probes, RNA-seq has no such requirement and can be performed on all species. The classical workflow (Figure 5) is based on an RNA extraction step, followed by a ribosomal RNA depletion or mRNA enrichment, followed by cDNA synthesis and short-read sequencing after a final sequencing adapter ligation (Hrdlickova et al., 2017). The proper quantification of the expression of each gene then relies on bioinformatics pipelines [START_REF] Corchete | Systematic comparison and assessment of RNA-seq procedures for gene expression quantitative analysis[END_REF]. Classically, the sequencing reads are filtered and aligned, with low quality or multi-mapped reads typically discarded after these steps. The expression quantification of each gene is then performed. The resulting data set is typically a raw count of the number of reads aligned to each gene. To formally compare mRNA abundance across datasets, the raw counts are transformed and normalized using, for example, the transcripts per million (TPM) unit [START_REF] Conesa | A survey of best practices for RNA-seq data analysis[END_REF]. The RNA-seq technique has been widely used to study various gene expression processes, such as RNA splicing [START_REF] Wang | Alternative isoform regulation in human tissue transcriptomes[END_REF] or RNA-mediated gene regulation (W. Li et al., 2016;[START_REF] Morris | The rise of regulatory RNA[END_REF]. Due to its versatility and its accessibility, RNA-seq has been an important tool to understand or diagnose gene expression alterations in various pathologies [START_REF] Byron | Translating RNA sequencing into clinical diagnostics: opportunities and challenges[END_REF][START_REF] Doebele | An Oncogenic NTRK Fusion in a Patient with Soft-Tissue Sarcoma with Response to the Tropomyosin-Related Kinase Inhibitor LOXO-101[END_REF][START_REF] Hong | RNA sequencing: new technologies and applications in cancer research[END_REF][START_REF] Wirka | Advances in Transcriptomics[END_REF] and to explore the transcriptome at the population level (Caudal et al., 2023;The GTEx Consortium, 2015). However, even though RNA-seq has considerably advanced the study of RNA abundance, the fact that this strategy was until recently solely based on short-read sequencing has made it to suffer from several limitations.

For example, in the case of mRNA isoforms, short reads sequencing prevents from accurate individual quantification [START_REF] Djebali | Landscape of transcription in human cells[END_REF]Stark et al., 2019). In organisms where transcripts can be very long and variable, such as humans, where more than half of the transcripts are longer than 2,500 bp [START_REF] Frankish | GENCODE reference annotation for the human and mouse genomes[END_REF], this issue is highly relevant and specific methods have emerged as powerful alternatives to short-read-based RNA-seq to account for mRNA isoforms.

Indeed, the development of new techniques to sequence long fragments of DNA (also known as long-read sequencing) is leading to new information on mRNA abundance (Figure 5). These methods, namely Pacific Biosciences (PacBio) and Oxford Nanopore Technologies (ONT) [START_REF] Rhoads | PacBio Sequencing and Its Applications[END_REF][START_REF] Wang | Nanopore sequencing technology, bioinformatics and applications[END_REF] allow the capture of the entire mRNA molecule, mainly through cDNA sequencing. They are particularly efficient for mRNA isoform detection and de novo transcriptome analysis (Stark et al., 2019). In addition, ONT can also be used to sequence native mRNA molecules, which allows for more precise exploration of mRNA base modifications. However, the throughput of these long-read sequencing methods is still limited.

Quantification of protein

Similar to mRNA, protein abundance quantification has constantly evolved since its inception in the second half of the 20 th century. Although some quantification methods are similar to what can be done for mRNA abundance, such as Western blotting [START_REF] Burnette | Western Blotting": Electrophoretic transfer of proteins from sodium dodecyl sulfatepolyacrylamide gels to unmodified nitrocellulose and radiographic detection with antibody and radioiodinated protein A[END_REF][START_REF] Towbin | Electrophoretic transfer of proteins from polyacrylamide gels to nitrocellulose sheets: procedure and some applications[END_REF], most of the quantification techniques are specific to protein abundance due to the very different chemical nature of peptide chains compared to RNA. Early protein quantification methods were based on spectrophotometric measurements, such as the Bradford or Lowry protein assay [START_REF] Bradford | A rapid and sensitive method for the quantitation of microgram quantities of protein utilizing the principle of protein-dye binding[END_REF][START_REF] Lowry | PROTEIN MEASUREMENT WITH THE FOLIN PHENOL REAGENT[END_REF]. However, these methods are not compatible for global quantification of all the different proteins in a sample. One of the first methods to Long-read mRNA sequencing Short-read mRNA sequencing allow for the quantification of the total set of proteins expressed in a cell, a tissue, or whole organism (such sets of protein are called proteomes) was based on 2-dimensional electrophoresis gels [START_REF] Magdeldin | Basics and recent advances of two dimensional-polyacrylamide gel electrophoresis[END_REF][START_REF] O'farrell | High Resolution Two-Dimensional Electrophoresis of Proteins[END_REF]. The principle is simple: all proteins in a sample are extracted and successively separated by two properties on the two dimensions of a gel (usually a polyacrylamide gel). The first dimension of the gel resolves protein molecules according to their isoelectric point (using a pH gradient), while the second dimension resolves them according to their molecular weight. The resulting gel consists of several separated dots that can be excised and quantified using, for example, a coupled mass spectrometry device. This technique has been repeatedly used to quantify the global proteome of many organisms, such as bacteria [START_REF] Wasinger | Progress with gene-product mapping of the Mollicutes: Mycoplasma genitalium[END_REF], yeast (Gygi et al., 1999) or humans [START_REF] Friedman | Proteome analysis of human colon cancer by two-dimensional difference gel electrophoresis and mass spectrometry[END_REF].

Since the 2000s, however, the evolution of both molecular biology tools and mass spectrometry methods has led to more precise and comprehensive proteome acquisitions. For example, tagbased protein quantifications have been widely used in the last decade. These techniques rely on the construction of libraries in which each protein is individually fused with a tag, such as green fluorescent protein (GFP) or a high-affinity epitope. The tagged proteins are quantified by immuno-or photo-detection. In yeast, several studies have quantified the proteome using either GFP microscopy [START_REF] Breker | A novel single-cell screening platform reveals proteome plasticity during yeast stress responses[END_REF][START_REF] Chong | Yeast Proteome Dynamics from Single Cell Imaging and Automated Analysis[END_REF][START_REF] Dénervaud | A chemostat array enables the spatio-temporal analysis of the yeast proteome[END_REF][START_REF] Mazumder | Genome-wide single-cell-level screen for protein abundance and localization changes in response to DNA damage in S. cerevisiae[END_REF][START_REF] Tkach | Dissecting DNA damage response pathways by analysing protein localization and abundance changes during DNA replication stress[END_REF][START_REF] Yofe | One library to make them all: streamlining the creation of yeast libraries via a SWAp-Tag strategy[END_REF], GFP flow cytometry [START_REF] Davidson | The proteomics of quiescent and nonquiescent cell differentiation in yeast stationary-phase cultures[END_REF][START_REF] Lee | Global protein expression profiling of budding yeast in response to DNA damage[END_REF][START_REF] Newman | Single-cell proteomic analysis of S. cerevisiae reveals the architecture of biological noise[END_REF] or tag immunodetection [START_REF] Ghaemmaghami | Global analysis of protein expression in yeast[END_REF]. However, these methods suffer from major shortcomings when it comes to largescale or unbiased quantification. First, library construction can be tedious, and the quantification itself requires each protein to be measured independently. Second, the fusion of a tag to a protein is obviously not biologically neutral, and this can lead to misfunction or dysregulation of abundance. Therefore, other methods related to mass spectrometry techniques have become prominent tool to explore and precisely resolve proteomes in an almost exhaustive manner. For example, liquid chromatography coupled with tandem mass spectrometry (LC-MS/MS) is particularly suited to multiplexed sample analysis, making it a perfect tool for proteome exploration. The method is generally based on double analysis of peptides (obtained after trypsin digestion and separated by LC) using two coupled mass spectrometers.

Quantification and identification of each protein is performed by analyzing the mass-to-charge (m/z) ratio of the peptide fragments on the two spectrometers. Over the past decade, several technological advances have improved the reliability of mass spectrometry proteomics. One notable example is the rise of data-independent acquisition (DIA) [START_REF] Chapman | Multiplexed and data-independent tandem mass spectrometry for global proteome profiling[END_REF][START_REF] Gillet | Mass Spectrometry Applied to Bottom-Up Proteomics: Entering the High-Throughput Era for Hypothesis Testing[END_REF], a method that allows for a broader range of protein quantification [START_REF] Li | Data-independent acquisition (DIA): An emerging proteomics technology for analysis of drug-metabolizing enzymes and transporters[END_REF].

Several software packages, such as DIA-NN, are specifically designed for large-scale proteomics experiments (Demichev et al., 2020). More importantly, these technological advances have enabled proteome exploration at a much higher throughput and scale in the recent years (Figure 6) (Messner et al., 2023(Messner et al., , 2022;;Muenzner et al., 2022), although proteomes at population level remain largely uncharacterized, especially compared to transcriptomes at population level.

The limitations of LC-MS/MS have been studied extensively, and one of the major limitations is missing data. Missing peptides is a common and long-standing issue in LC-MS/MS studies, especially for large-scale exploration [START_REF] Karpievitch | Liquid Chromatography Mass Spectrometry-Based Proteomics: Biological and Technological Aspects[END_REF]Muenzner et al., 2022). There are several reasons for this. For example, the abundance level of a protein is a major determinant of its detectability: low abundance peptides are often missed by LC-MS/MS, resulting in a biased quantification towards highly abundant proteins. Also, the chemical and physical characteristics of some proteins make them prone to be trapped in the LC column.

Sample n

Protein abundance bioinformatic tools such as DIA-NN (Demichev et al., 2020)

Ribosome profiling

Despite the interconnected nature of mRNA and protein abundance, transcriptomic and proteomic regulation can be very different. Therefore, understanding the quantitative relationship between the two expression levels requires a precise analysis of the translation process. While there are several techniques that allow the exploration of translation [START_REF] Dermit | Methods for monitoring and measurement of protein translation in time and space[END_REF], such as polysome profiling [START_REF] Arava | Genome-wide analysis of mRNA translation profiles in Saccharomyces cerevisiae[END_REF], a more precise and robust method was developed more than a decade ago, called ribosome profiling or ribo-seq (Figure 7A) (Ingolia et al., 2009). The principle is to extract intact polysomes from a cell and subject the cell extract to RNase digestion. The mRNA fragments that are actively being translated are protected from the RNase treatment by the translating ribosomes, resulting in a pool of ≈28 nucleotides of mRNA fragments. These fragments are then sequenced to determine which parts of the transcriptome are being translated. Usually, ribo-seq and RNA-seq experiments are performed on the same sample.

Thanks to the mechanisms underlying the translation process, the ribosome profiling results

have some specific features. First, the average length of the reads corresponds to the size of the mRNA covered by the ribosome, i.e., 28 nucleotides (Ingolia, 2010). Second, due to the codonwise movement of the ribosome along the mRNA, ribosome profiling libraries exhibit a characteristic 3-nucleotide periodicity (Figure 7B), supporting that single codon resolution is achievable with this technique (Ingolia et al., 2009). In addition, the distribution of reads along the gene sequence can be indicative of several translational features, such as changes in elongation rate and ribosomal frameshifting or stalling (Figure 7C) (Ingolia et al., 2019;[START_REF] Michel | Observation of dually decoded regions of the human genome using ribosome profiling data[END_REF][START_REF] Napthine | Protein-directed ribosomal frameshifting temporally regulates gene expression[END_REF]. Because of this versatility, ribosome profiling has been a powerful tool for the precise study of translation (Brar and Weissman, 2015). In fact, studies have been conducted on quantitative, mechanistic, and spatial aspects of the translation process [START_REF] Guydosh | Dom34 rescues ribosomes in 3' untranslated regions[END_REF]G.-W. Li et al., 2014;[START_REF] Williams | Targeting and plasticity of mitochondrial proteins revealed by proximity-specific ribosome profiling[END_REF]. Finally, ribosome profiling has also been a powerful tool for exploring functional genome evolution, and notably how de novo genes are expressed and fixed (Blevins et al., 2021;[START_REF] Wacholder | A vast evolutionarily transient translatome contributes to phenotype and fitness[END_REF].

Currently, ribosome-profiling still faces several limitations (Brar and Weissman, 2015). At the experimental level, rapidly stopping translating ribosomes to obtain an accurate snapshot of translation is a critical step. Translation elongation inhibitors, such as cycloheximide, have been widely used, but they are known to induce ribosome distribution biases, particularly around the translation start site [START_REF] Guydosh | Dom34 rescues ribosomes in 3' untranslated regions[END_REF][START_REF] Hussmann | Understanding Biases in Ribosome Profiling Experiments Reveals Signatures of Translation Dynamics in Yeast[END_REF]Ingolia et al., 2009). In this regard, flash freezing with liquid nitrogen is a robust alternative to effectively stop the ribosome movement [START_REF] Ingolia | The ribosome profiling strategy for monitoring translation in vivo by deep sequencing of ribosome-protected mRNA fragments[END_REF]. Other experimental biases are known, such as contamination of mRNA fragments that are not ribosome-protected fragments, but rather structured RNA. In silico data processing is usually required to address such issues [START_REF] Ingolia | Ribosome profiling: new views of translation, from single codons to genome scale[END_REF]. Finally, one of the main limitations of ribosome profiling is that due to the higher level of sample processed (compared to mRNA), large amounts of cellular material are required to accurately quantify translation [START_REF] Ingolia | The ribosome profiling strategy for monitoring translation in vivo by deep sequencing of ribosome-protected mRNA fragments[END_REF]. As a result, ribosome profiling remains challenging to scale for single-cell approaches or high-throughput studies. However, recent advances in the sensitivity of the method and in the analytical workflow, especially with the incorporation of machine learning steps, have enabled the development of single-cell ribosome profiling [START_REF] Vaninsberghe | Singlecell Ribo-seq reveals cell cycle-dependent translational pausing[END_REF]. Regarding the adaptation of ribosome profiling for largescale studies, the linear amplification of ribosome-protected fragments after their extraction is a promising solution and may allow the exploration of translation at the population scale (Mito et al., 2023).

Analytical link between genotype and phenotype variations

Accurately determining the genetic loci involved in phenotypic variation has been a constant and arduous challenge over the past 50 years. The development of molecular genetics and biology, together with advances in DNA sequencing through the Sanger technique [START_REF] Sanger | DNA sequencing with chain-terminating inhibitors[END_REF], led to classical approaches to link genotypes and phenotypes, such as random mutagenesis. In this technique, chemical or physical mutagenesis agents were applied to individuals such as mice or yeast to induce point mutations. These latter were then associated to wide ranges of phenotypic changes, such as cell division in the fission yeast Schizosaccharomyces pombe [START_REF] Nurse | Genetic control of the cell division cycle in the fission yeast Schizosaccharomyces pombe[END_REF][START_REF] Nurse | Regulatory genes controlling mitosis in the fission yeast Schizosaccharomyces pombe[END_REF] or motility in the nematode Caenorhabditis elegans [START_REF] Brenner | The genetics of Caenorhabditis elegans[END_REF]. In parallel, the advent of the PCR amplification techniques allowed more targeted genetic disruption, especially in model organisms with efficient recombination capacities, such as S. cerevisiae [START_REF] Shortle | Lethal Disruption of the Yeast Actin Gene by Integrative DNA Transformation[END_REF][START_REF] Wach | New heterologous modules for classical or PCRbased gene disruptions in Saccharomyces cerevisiae[END_REF]. In this case, systematic deletion (or knock-out, KO) of each gene in the genome allows precise characterization of the effect of each gene [START_REF] Giaever | The Yeast Deletion Collection: A Decade of Functional Genomics[END_REF].

More recently, targeted mutagenesis was deeply developed with the introduction of CRISPR-CAS9 technologies [START_REF] Cong | Multiplex genome engineering using CRISPR/Cas systems[END_REF][START_REF] Doudna | The new frontier of genome engineering with CRISPR-Cas9[END_REF][START_REF] Jinek | A programmable dual-RNA-guided DNA endonuclease in adaptive bacterial immunity[END_REF].

CRISPR-CAS9 allows for highly efficient and nucleotide resolution modifications that are difficult to achieve in human cells, for example.

However, while these methods are highly efficient for characterizing the cellular or molecular effects of particular genetic variants, they are not well suited for exploring the influence of natural genetic variation on phenotype, especially for complex traits. Such investigations require either a large cohort of individuals or a model organism for which large-scale segregant generation is possible. Genetic association studies are the main strategy for such investigations.

Two main types of association studies are commonly used: linkage analysis (or linkage mapping) and genome-wide association studies (GWAS). The goal of these strategies is to detect and associate genetic regions with quantitative phenotypic variation (e.g., human height, mRNA or protein abundance, growth on a particular medium). The results of these studies are called quantitative trait loci (or QTL). Depending on the method used to link genetic regions to quantitative phenotypes, the resolution of QTLs can vary from large chromosomal regions to single nucleotide variants. Gene expression, specifically mRNA and protein abundance, can be used as a phenotype to investigate the genetic origins of gene expression variation between individuals. In the case of mRNA abundance, the detected loci are usually referred to as eQTL (for expression QTL) [START_REF] Nica | Expression quantitative trait loci: present and future[END_REF]. In the case of protein abundance, they are referred as pQTL (for protein QTL) (Ferkingstad et al., 2021).

Linkage mapping studies

Linkage mapping studies are based on the generation of a large number of segregants from an original cross between two genetically and phenotypically distinct individuals (Figure 8A). Due to meiotic recombination, the parental genetic variants are shuffled among the offspring, resulting in unique genotype in each segregant. By combining the phenotypic measurement of both parents and segregants with their genotypes (Figure 8B), it is possible to recover the causal regions of the phenotype under study. In fact, individuals with a similar phenotype will most likely share the genomic regions that carry the causal variants, while the vast majority of the other parental variants will be randomly distributed along the genome (Figure 8C). While early linkage mapping studies relied on genetic markers such as restriction fragment length polymorphisms [START_REF] Botstein | Construction of a genetic linkage map in man using restriction fragment length polymorphisms[END_REF], most linkage mapping analyses are based on SNPs (Albert et al., 2018;[START_REF] Brem | Genetic dissection of transcriptional regulation in budding yeast[END_REF]. More recently, SVs were also integrated in such analyses (Weller et al., 2023). The budding yeast S. cerevisiae is an important tool for linkage mapping studies because it combines several characteristics that are crucial for an efficient QTL detection: a small genome, a short sexual generation time and, more importantly, a high meiotic recombination rate [START_REF] Fay | The molecular basis of phenotypic variation in yeast[END_REF][START_REF] Liti | Advances in quantitative trait analysis in yeast[END_REF]. (A) In linkage mapping studies, 2 individuals with different genotypes and phenotypes are crossed and the phenotype of their offspring is associated with genetic markers (M1, M2... M5) distributed throughout the genome (e.g., SNPs). (B) A genetic marker associated with the phenotype difference between the parents will also show this association in the segregants. (C) All markers ultimately have an association score (e.g., LOD, logarithm of odds) that must exceed a certain threshold to be considered significant.

Disentangling the genetic origins of inter-individual variation in gene expression has often been investigated using linkage mapping. Studies of mRNA abundance have been performed in several organisms, ranging from plants [START_REF] Jansen | Genetical genomics: the added value from segregation[END_REF] to animals [START_REF] Schadt | Genetics of gene expression surveyed in maize, mouse and man[END_REF] and of course, yeast (Albert et al., 2018;[START_REF] Brem | Genetic dissection of transcriptional regulation in budding yeast[END_REF], resulting in the discovery of thousands of eQTL. Protein abundance has also been explored and genetically mapped using linkage mapping studies. For example, in yeast, linkage mapping and proteomics have been combined twice (Albert et al., 2014;Foss et al., 2007), but due to technical limitations, the number of proteins included in these studies was limited. In the case of the earlier study (Foss et al., 2007), the LC-MS/MS used at the time had a high signal-to-noise ratio, which made it difficult to repeatedly cover a large number of proteins in the samples. In the end, 221 proteins were used for linkage mapping. In the latter study (Albert et al., 2014), the proteins were quantified using a single-cell measurement, where the proteins are fused to GFP and the signal is measured using fluorescence-activated cell sorting. To ensure good statistical power and the possibility to control the results, only 160 proteins were ultimately used.

While linkage mapping approach is a powerful method to determine the genomic region associated with a phenotype, it has some inherent limitations. First, the genetic diversity captured by this method is limited to that found in the two parental individuals which does not recapitulate the complete genetic diversity of the species. Also, depending on the recombination rate of the organism under study, the resolution of linkage mapping may be limited, especially if the study focuses on the first generation of offspring [START_REF] Flint | Strategies for mapping and cloning quantitative trait genes in rodents[END_REF]. As a result, large regions are associated with the phenotype under study, and precise identification of the causal variant can be tedious. Several tools have been developed specifically to address these issues.

For example, in mice, an entire lineage has been generated from 8 inbred individuals (Collaborative Cross Consortium, 2012) leading to an increase in both genetic diversity and resolution [START_REF] Gatti | Quantitative Trait Locus Mapping Methods for Diversity Outbred Mice[END_REF]. In yeast, the generation of large population by crossing 16

genetically distinct individuals has also been achieved to tackle the genetic diversity limitation [START_REF] Bloom | Rare variants contribute disproportionately to quantitative trait variation in yeast[END_REF]. Given the limited resolution of linkage mapping studies, one possible solution is to map the genetic origin of a phenotype using the F6 generation of the parental cross [START_REF] Jakobson | Molecular Origins of Complex Heritability in Natural Genotype-to-Phenotype Relationships[END_REF].

Genome-wide association studies

The development of large-scale sequencing projects over the last two decades has led to a clearer view of intraspecific genetic diversity. In the mid-2000s, using large cohorts of individuals whose allelic status is clearly defined across their genome, a new method for linking genotypes to phenotypes was developed: genome-wide association studies (GWAS) [START_REF] Dewan | HTRA1 promoter polymorphism in wet age-related macular degeneration[END_REF][START_REF] Klein | Complement factor H polymorphism in age-related macular degeneration[END_REF]Wellcome Trust Case Control Consortium, 2007). This method is based on testing the association of hundreds of thousands or even millions of variants with phenotypes of interest. Until recently, the tested genetic variants were mostly SNPs [START_REF] Uffelmann | Genome-wide association studies[END_REF] sometimes completed with CNVs data (Wellcome Trust Case Control Consortium, 2010). Overall, the experimental workflow of a GWAS follows several steps (Figure 9A). First, a large population is genotyped to resolve the allelic status of each individual. Historically, microarray data were used to resolve the genotype of each individual, but nowadays whole-genome sequencing is preferred and has the advantage to capture nearly all genotypic variation across the genome. Each of the variants (generally the bi-allelic SNPs) is then tested for association with the phenotype of interest. There are several ways to perform the association, among which the used of either linear mixed regression models (Figure 9A) or logistic mixed regression models, depending on the nature of the phenotype (either continuous or discrete) [START_REF] Uffelmann | Genome-wide association studies[END_REF]. As population structure or familial relatedness is a major confounding factor in GWAS [START_REF] Balding | A tutorial on statistical methods for population association studies[END_REF][START_REF] Kang | Variance component model to account for sample structure in genome-wide association studies[END_REF][START_REF] Zhang | Mixed linear model approach adapted for genome-wide association studies[END_REF][START_REF] Zhao | An Arabidopsis Example of Association Mapping in Structured Samples[END_REF], linear mixed models (LMM) are powerful statistical methods to correct for such confounding. As LMM can be computationally intensive, several tools have been developed to increase the accessibility of these methods, such as TASSEL or FaST-LMM (Lippert et al., 2011;[START_REF] Zhang | Mixed linear model approach adapted for genome-wide association studies[END_REF].

Because associations are made across thousands to millions of SNPs, false discovery rate correction is a critical part of GWAS. Depending on the population and the species, different types of significance correction are used. In humans, a common method is to use a Bonferroni corrected p-value to detect significant associations, resulting in a p-value threshold of 5×10 -8 [START_REF] Uffelmann | Genome-wide association studies[END_REF]. For more complex GWAS focusing on multiple phenotypes simultaneously, a trait-specific p-value can be defined by performing permutation tests (Caudal et al., 2023;Peter et al., 2018). The results of a GWAS can be easily visualized using a Manhattan plot (Figure 9A), where for each SNP across the chromosomes (x-axis), an association score (usually the -log10 transformation of the association p-value) is plotted (yaxis). Since its development, GWAS have been used to study a wide variety of phenotypes. As of April 2023, more than 6,000 human GWAS have been published, for which more than 500,000 associations has been highlighted (Figure 9B) [START_REF] Macarthur | The new NHGRI-EBI Catalog of published genome-wide association studies (GWAS Catalog)[END_REF]. Human diseases have naturally caught lot of attention and examples include cancer [START_REF] Sud | Genome-wide association studies of cancer: current insights and future perspectives[END_REF], type 2 diabetes (Zhao et al., 2017), and psychological or mental disorders [START_REF] Duncan | Significant Locus and Metabolic Genetic Correlations Revealed in Genome-Wide Association Study of Anorexia Nervosa[END_REF][START_REF] Hyde | Identification of 15 genetic loci associated with risk of major depression in individuals of European descent[END_REF][START_REF] Jansen | Genome-wide analysis of insomnia in 1,331,010 individuals identifies new risk loci and functional pathways[END_REF][START_REF] Li | Genome-wide association analysis identifies 30 new susceptibility loci for schizophrenia[END_REF]. Non-pathological complex phenotypes can also be investigated with GWAS such as body mass index [START_REF] Yengo | Meta-analysis of genome-wide association studies for height and body mass index in ∼700000 individuals of European ancestry[END_REF], educational attainment [START_REF] Lee | Gene discovery and polygenic prediction from a genome-wide association study[END_REF], or even musical beat synchronization [START_REF] Niarchou | Genome-wide association study of musical beat synchronization demonstrates high polygenicity[END_REF].

human phenotype. The colors represent the type of phenotype, see www.ebi.ac.uk/gwas/ for more details. Plot and data were generated and obtained from the NHGRI-EBI catalog [START_REF] Macarthur | The new NHGRI-EBI Catalog of published genome-wide association studies (GWAS Catalog)[END_REF].

Of course, gene expression itself has been studied with GWAS. In this case, a hundred to a few thousand phenotypes are analyzed simultaneously, i.e., mRNA or protein abundance.

Technically, GWAS focusing on gene expression can be tedious to perform, as gene expression exploration (e.g., RNA sequencing or LC-MS/MS) has to be performed in each individual. In addition, as mentioned above, specific significance thresholds need to be set in such studies, as the risk of false-positive discovery is high [START_REF] Liu | Abundant associations with gene expression complicate GWAS follow-up[END_REF]. A common way to visualize GWAS performed on mRNA or protein abundance is to plot the genomic location of the eQTL or pQTL against the location of its associated trait (i.e., the affected gene) (Figure 10A). Several studies have focused on the genetic origins of mRNA abundance, the most famous in human being the Genotype-Tissue Expression (GTEx) (The GTEx Consortium, 2020[START_REF] Auesukaree | Molecular mechanisms of the yeast adaptive response and tolerance to stresses encountered during ethanol fermentation[END_REF], 2015), in which transcript abundance was monitored in 49 human tissues from 838 postmortem donors.

Recently, a catalog of human eQTL has been established [START_REF] Kerimov | A compendium of uniformly processed human gene expression and splicing quantitative trait loci[END_REF]. Plants [START_REF] Lan | AtMAD: Arabidopsis thaliana multi-omics association database[END_REF] and of course, yeast (Caudal et al., 2023) have also been used for eQTL exploration by GWAS. Fewer large-scale studies have been performed for pQTL exploration with GWAS [START_REF] Suhre | Genetics meets proteomics: perspectives for large population-based studies[END_REF], mainly because large-scale quantification of protein abundance has been a major limiting factor. Notable studies focused on the plasma and serum proteome (Ferkingstad et al., 2021;[START_REF] Gudjonsson | A genome-wide association study of serum proteins reveals shared loci with common diseases[END_REF]. Until now, finding a good trade-off between the number of samples in a proteomic GWAS and the number of proteins included has been tedious. Despite the significant progress that has been made with GWAS, the technique still suffers from limitations [START_REF] Tam | Benefits and limitations of genomewide association studies[END_REF]. First, due to the stringency applied to avoid false-positive associations, genetic variants are likely to be missed. Therefore, only a fraction of the

Gene position

Local Distant A B heritability of a complex trait is captured by GWAS [START_REF] Dudbridge | Estimation of significance thresholds for genomewide association scans[END_REF]Manolio et al., 2009) and this fraction will be biased towards high effect variants. Conversely, small effect variants are more difficult to detect. In addition, a large fraction of the heritability of traits is also missed because several sources of genetic variants are poorly considered (Manolio et al., 2009;[START_REF] Zuk | Searching for missing heritability: Designing rare variant association studies[END_REF]. For example, SVs are rarely accounted for in GWAS, mainly because of limitations for their accurate detection that were only recently overcome, but also because of the difficulty to integrate these data in GWAS. The democratization of long-read sequencing and the development of pangenome graph-based GWAS are promising solutions to account for the SV effect (He et al., 2023;Li et al., 2022;Logsdon et al., 2020;Zhou et al., 2022). Rare variants are also poorly considered in GWAS because their low representation in the cohorts studied makes it difficult to detect their association with sufficient statistical power. In this regard, the construction of diallel crosses allows to artificially increase the allele frequency of rare variants (Fournier et al., 2019). Finally, due to linkage disequilibrium (the non-random association between two or more alleles on different loci), finding the correct causal variant can be difficult [START_REF] Altshuler | Genetic mapping in human disease[END_REF], especially when the QTL is in a non-coding region, and sometimes requires additional studies to confirm the effect of a genetic variant. For this reason, the predictive power of GWAS in the clinical context remains limited [START_REF] Janssens | Genome-based prediction of common diseases: advances and prospects[END_REF][START_REF] Loos | Predicting Polygenic Obesity Using Genetic Information[END_REF]. Increasing the size of cohorts is a good solution to overcome such difficulties [START_REF] Tam | Benefits and limitations of genomewide association studies[END_REF]. In this context, the development of high-throughput molecular phenotyping (Caudal et al., 2023;Messner et al., 2022) is a promising approach to overcome some of the current limitations of GWAS.

Gene expression and SNP

As described above, both linkage mapping approaches and GWAS have been used to dissect the genetic origins of variations in transcript and protein abundance across individuals. At the transcriptomic level, extensive efforts have been made to understand which and how genetic variants, more specifically SNPs, affect mRNA abundance (Caudal et al., 2023;Foss et al., 2007;[START_REF] Gan | Multiple reference genomes and transcriptomes for Arabidopsis thaliana[END_REF]Ghazalpour et al., 2011;The GTEx Consortium, 2015). These studies have led to the discovery of thousands of eQTL that affect gene expression through a variety of mechanisms. Some eQTL affect gene expression of nearby genes, while others affect their target genes in a distant manner (Figure 11). The former, also known as local-or cis-eQTL (or local-and cis-pQTL in the case of protein abundance), usually affects the ability of the transcription machinery to bind to the promoter of its target gene (Figure 11A). Any element involved in the direct regulation of a gene's transcription can be affected, such as the core promoter [START_REF] Lubliner | Core promoter sequence in yeast is a major determinant of expression level[END_REF][START_REF] Tirosh | Promoter architecture and the evolvability of gene expression[END_REF], enhancer regions [START_REF] Garieri | The effect of genetic variation on promoter usage and enhancer activity[END_REF][START_REF] Kikuchi | Enhancer variants associated with Alzheimer's disease affect gene expression via chromatin looping[END_REF], nearby chromatin accessibility [START_REF] Keele | Integrative QTL analysis of gene expression and chromatin accessibility identifies multi-tissue patterns of genetic regulation[END_REF], or terminal regions [START_REF] Hill | Molecular and evolutionary processes generating variation in gene expression[END_REF]. In addition, the local-eQTL tend to have a greater effect on their target genes compared to the distant-eQTL (Albert et al., 2018;Caudal et al., 2023). In humans, late lactose tolerance is a famous case of local DNA variants inducing a change of gene expression, resulting in a persistent lactase expression. Briefly, while populations with the ancestral allele have a decreased expression of the lactase (LCT) gene after childhood, single mutations (either alone or in combination) in an enhancer region of LCT create a new binding site for a transcription factor leading to a non-downregulated pathway for LCT expression [START_REF] Enattah | Identification of a variant associated with adult-type hypolactasia[END_REF][START_REF] Fang | The human lactase persistence-associated SNP -13910*T enables in vivo functional persistence of lactase promoter-reporter transgene expression[END_REF][START_REF] Lewinsky | T-13910 DNA variant associated with lactase persistence interacts with Oct-1 and stimulates lactase promoter activity in vitro[END_REF][START_REF] Olds | Lactase persistence DNA variant enhances lactase promoter activity in vitro: functional role as a cis regulatory element[END_REF]. Graphically, cis regulation is easily observed as a diagonal line when plotting the eQTL position vs the target gene position (reflecting a similar location of the QTL and the trait on the genome, Figure 10A).

Conversely, eQTL or pQTL described as distant (or trans-) can be located anywhere in the genome, either on the same or different chromosomes (Figure 11B). The effects of distant-QTL are usually achieved through the proteins or RNA involved in transcriptional or translational regulation such as mRNA binding protein, transcription factors or non-coding RNA [START_REF] He | Genome-wide pQTL analysis of protein expression regulatory networks in the human liver[END_REF][START_REF] Lutz | DNA variants affecting the expression of numerous genes in trans have diverse mechanisms of action and evolutionary histories[END_REF]. Interestingly, distant-QTL tend to be located in hotspots (Albert et al., 2018;[START_REF] Qu | Lead Modulates trans-and cis-Expression Quantitative Trait Loci (eQTLs) in Drosophila melanogaster Heads[END_REF][START_REF] Yao | Dynamic Role of trans Regulation of Gene Expression in Relation to Complex Traits[END_REF] and have more pleiotropic effects than local regulatory variants [START_REF] Lemos | Dominance and the evolutionary accumulation of cis-and trans-effects on gene expression[END_REF][START_REF] Prud'homme | Emerging principles of regulatory evolution[END_REF]. On a plot of QTL position versus affected gene position, distant QTL hotspots are graphically indicated by the vertical accumulation of points (Figure 10B). Due to their multiple effects, distant regulatory variants tend to be more deleterious and less beneficial than local regulatory variants [START_REF] Coolon | Molecular Mechanisms and Evolutionary Processes Contributing to Accelerated Divergence of Gene Expression on the Drosophila X Chromosome[END_REF][START_REF] Emerson | Natural selection on cis and trans regulation in yeasts[END_REF][START_REF] Schaefke | Inheritance of gene expression level and selective constraints on trans-and cis-regulatory changes in yeast[END_REF]. Therefore, they are an important driver of the relationship between diseases and their genetic origins [START_REF] Westra | Systematic identification of trans eQTLs as putative drivers of known disease associations[END_REF], as it is the case for several autoimmune pathologies such as type 1 diabetes or systemic lupus erythematosus [START_REF] Han | Genome-wide association study in a Chinese Han population identifies nine new susceptibility loci for systemic lupus erythematosus[END_REF][START_REF] Heinig | A trans-acting locus regulates an anti-viral expression network and type 1 diabetes risk[END_REF]. regulatory variants affect the expression of genes located on another chromosome or far away on the same chromosome. This is usually achieved through the action of an intermediary protein.

Gene expression and SV

Structural variants are also major contributors to variability in mRNA or protein abundance.

Down syndrome (or trisomy 21) is a notable example, where individuals carry a full or partial extra copy of the 21st chromosome. This results in a huge macroscopic phenotypic impact and increased risk for a wide range of diseases [START_REF] Patterson | Molecular genetic analysis of Down syndrome[END_REF]. At the molecular level, the consequences of this aneuploidy are extensive, affecting gene expression of genes located on all chromosomes [START_REF] Letourneau | Domains of genome-wide gene expression dysregulation in Down's syndrome[END_REF][START_REF] Prandini | Natural Gene-Expression Variation in Down Syndrome Modulates the Outcome of Gene-Dosage Imbalance[END_REF]. In plants, a large-scale study of 100 tomato lines revealed that SVs (which in this organism are mainly associated with transposons) are a major source of variation in gene expression across the species (Alonge et al., 2020). Another elegant example of this is the case of the sulfite tolerance in [START_REF] Ribéreau-Gayon | The Microbiology of Wine and Vinifications[END_REF]. This enhanced tolerance was associated with a reciprocal translocation between chromosomes VIII and XVI, resulting in an overexpression of the SSU1 gene, a sulfite pump whose promoter is altered by the translocation (Pérez-Ortín et al., 2002). Later studies showed that three different chromosomal rearrangements (two translocations and one inversion) can induce SSU1 overexpression in different yeast isolates (García-Ríos et al., 2019;Yuasa et al., 2004;Zimmer et al., 2014). However, due to the tedious nature of detecting and characterizing SVs in large populations, their impact on gene expression, especially at the population level, remains largely unexplored. Recently however, CNV has been taken into consideration to explore variation in transcript abundance in yeast (Caudal et al., 2023).

The transcriptome and proteome relationship

In the last two decades, both technological (e.g., RNA sequencing, large-scale LC-MS/MS) and analytical (e.g., GWAS) developments have led to a better understanding of the molecular mechanisms underlying the genotype-phenotype relationship. However, when both transcript and protein abundance studies are numerous, their conclusions can be very contradictory, especially when it comes to the relationship between the transcriptome and the proteome. As mentioned previously, the apparent linear and hierarchical nature of gene expression hides a complex and tightly regulated phenomenon (Buccitelli and Selbach, 2020;Liu et al., 2016;[START_REF] Vogel | Insights into the regulation of protein abundance from proteomic and transcriptomic analyses[END_REF]. Thus, the final protein abundance at any-given time in a cell is the result of a subtle balance between transcription rate, mRNA half-life, translation rate, protein half-life, and the cell cycle progression [START_REF] Baum | Of Gene Expression and Cell Division Time: A Mathematical Framework for Advanced Differential Gene Expression and Data Analysis[END_REF]Buccitelli and Selbach, 2020).

Therefore, many aspects of the relationship between the transcriptome and the proteome are poorly understood. For example, how well protein abundance can be predicted from transcript abundance is still an ongoing debate. Similarly, it is still unclear if and how proteome variation reflects transcriptome variation. Large-scale exploration of mRNA and protein abundance is a promising tool to address these grey areas, as it allows precise quantification of expression variation and deep exploration of the genetic origin of these variations. However, at the population level, studies focusing on the relationship between the transcriptome and the proteome are sparse.

Transcript-protein correlation

A good way to explore the interaction between the transcriptome and the proteome is to focus on the correlation between mRNA and protein abundance. An important distinction must be made, as two types of mRNA-protein correlations can be calculated (Buccitelli and Selbach, 2020;Liu et al., 2016): the first focuses on the correlation of all genes in a sample (e.g., a tissue, a cell, a strain...), while the other highlights the variation in mRNA and protein abundance across different conditions. The first type is referred to as the "across-gene" correlation while the second is referred to as the "within-gene" correlation. The confusion between these two types of correlations is common, even in the scientific literature (Fortelny et al., 2017).

Correlations are often calculated using either Spearman or Pearson coefficients, and the resulting value is highly dependent on the type of correlation (i.e., across-or within-gene).

Across-gene correlation

As explained above, the across-gene correlation focuses of the comparison between the mRNA and protein levels within a single sample. Graphically, this is often analyzed by plotting the transcript and peptide abundance together (Figure 12A), where each point corresponds to a gene. The across-gene correlation has been extensively studied in several species such as humans (Battle et al., 2015;Edfors et al., 2016;Gautier et al., 2016;Salovska et al., 2020;Wang et al., 2019;Wilhelm et al., 2014;Zhang et al., 2014), mice and rats (Aydin et al., 2023;J. J. Li et al., 2014;[START_REF] Moritz | Poor transcript-protein correlation in the brain: negatively correlating gene products reveal neuronal polarity as a potential cause[END_REF]Schwanhäusser et al., 2011), fruit flies (Becker et al., 2018), maize (Ponnala et al., 2014) and, of course, yeast (Gygi et al., 1999;Ingolia et al., 2009;Marguerat et al., 2012). In the vast majority of these studies, across-gene correlation typically shows medium-high to high correlation indexes (0.4-0.8) (Figure 12B). This observation is consistent through all the species in which this has been examined. Overall, this means that highly abundant transcripts encode for highly abundant proteins. It is important to note that across-gene correlation may be sensitive to the time point at which the mRNA and proteome abundances were surveyed. Indeed, the hierarchical nature of gene expression will introduce a delay between transcription and translation [START_REF] Fournier | Delayed Correlation of mRNA and Protein Expression in Rapamycin-treated Cells and a Role for Ggc1 in Cellular Sensitivity to Rapamycin*[END_REF][START_REF] Gedeon | Delayed Protein Synthesis Reduces the Correlation between mRNA and Protein Fluctuations[END_REF].

Therefore, when biological samples are exposed to changing conditions, a change in gene expression will follow, and for most genes, this will first affect transcription and then translation. Ultimately, two measures of mRNA and protein abundance during the steady state and transition phase will most likely have different correlation levels, with the steady state correlation index being higher than the transition state correlation index. Also, biological features such as the group of genes being monitored or the context of the study (e.g., tissue, cell type...) are known to influence the mRNA-protein correlation level (Buccitelli and Selbach, 2020). For example, in a study of several cell lines and tissues targeting specific proteins, the Pearson coefficient between mRNA and protein abundance varied from 0.39 to 0.79 (Edfors et al., 2016). Overall, this type of correlation is a good indicator to empirically quantify the impact of mRNA abundance on the final protein abundance (Liu et al., 2016). Investigations on mouse data showed that between 56% and 84% of the variance in protein levels is explained by mRNA levels (J. J. Li et al., 2014). This highlights that at steady state and within a sample, protein abundance is primarily explained by mRNA abundance.

Within-gene correlation

While across-gene correlation provides valuable information about how mRNA and protein levels are coupled in a specific sample, it does not provide information about how transcript and protein abundance variations are reflected at larger scales. Within-gene correlation is a powerful tool to interrogate and investigate gene expression changes across multiple samples (tissues, cell types, strains, growth conditions...). In this approach, the mRNA-protein correlation is calculated for each gene using the sample transcript and protein levels (Figure 13). Interestingly, there is no clear consensus across studies for the within gene correlation (Archer et al., 2018;Aydin et al., 2023;Battle et al., 2015;Chick et al., 2016;Ghazalpour et al., 2011;Huang et al., 2017;[START_REF] He | Genome-wide pQTL analysis of protein expression regulatory networks in the human liver[END_REF]Mertins et al., 2016;Mirauta et al., 2020;Mun et al., 2019;Vasaikar et al., 2019;Wang et al., 2019;Zhang et al., 2014Zhang et al., , 2016)). Within gene correlation indexes range from 0.14 to 0.59 (Battle et al., 2015;Upadhya and Ryan, 2022), so there is no clear consensus on whether mRNA and protein changes are correlated or not. In the within-gene correlation, the transcript and protein abundances of all the individuals are correlated for each gene. Graphically, a plot can be drawn for each gene where each dot corresponds to a sample. The correlation indexes are calculated for each gene and used to get an overall correlation median.

There are several reasons for this uncertainty. Technical biases can be an influencing factor, as proteins that are preferentially captured by previous proteomic methods will show up at higher correlation levels (Alam et al., 2016;Upadhya and Ryan, 2022). The type of data used in the correlation calculation (i.e., absolute, or relative mRNA and protein quantification) is also a determinant of the overall correlation levels. Indeed, absolute transcript and protein levels span several orders of magnitude, while the relative expression change of protein across samples remains in a much narrower range (Marguerat et al., 2012;Messner et al., 2023). Similarly, the magnitude of variation across the samples also strongly influences the mRNA-protein correlation: genes with important expression variation are more likely to have a high withingene correlation because the changes are more likely to affect both expression levels (Buccitelli and Selbach, 2020;Wang et al., 2019). In part because of this, the cellular function of the gene is also an important determinant of within-gene correlation. For example, metabolism-related genes tend to be associated with high levels of correlation (Buccitelli and Selbach, 2020;Wang et al., 2019). It is worth noting that these genes are known to have highly variable expression across individuals (Caudal et al., 2023). Conversely, ribosomal genes tend to show no correlation or slight anticorrelation (Buccitelli and Selbach, 2020;Wang et al., 2019).

Expression noise is also a confounding factor for mRNA-protein correlation, especially for low expressed genes where biological signal variations are of the same magnitude as the noise. More generally, the precision of proteome and transcriptome exploration will strongly influence the mRNA-protein correlation (Buccitelli and Selbach, 2020). Finally, the aforementioned investigation included at most a few hundred samples (192 mouse samples in Chick et al., 2016). To truly explore both proteome and transcriptome variation at the population level, larger cohorts of individuals are required. Until recent advances (Messner et al., 2022), this was technically difficult to achieve.

Post-transcriptional buffering

Although tightly regulated, gene expression is a noisy process, especially at the transcriptional level. External factors can induce gene expression noise through surface receptors, as cells live in a highly fluctuating environment where important changes in condition must be distinguished from rapid and noisy signals (Liu et al., 2016). Similarly, internal factors such as random transcription initiation can also lead to expression noise [START_REF] Chalancon | Interplay between gene expression noise and regulatory network architecture[END_REF][START_REF] Gandhi | Transcription of functionally related constitutive genes is not coordinated[END_REF]. Erroneous or inappropriate gene expression can ultimately lead to proteome imbalance and is obviously detrimental to proper cellular homeostasis. If the external factors can be compensated by annealing incorrect cellular signals before they affect gene expression [START_REF] Chalancon | Interplay between gene expression noise and regulatory network architecture[END_REF][START_REF] Hornung | Noise Propagation and Signaling Sensitivity in Biological Networks: A Role for Positive Feedback[END_REF], the internal factors are more likely to activate transcription. Cells must therefore cope with expression noise. More generally, the overall effect of gene expression variations can alter several key cellular functions.

Interestingly, some of these central functions tend to show robustness to expression variation [START_REF] Félix | Pervasive robustness in biological systems[END_REF]. Post-transcriptional buffering is a good example of the cellular mechanisms that deal with deleterious transcriptional noise or variation.

Different contexts, one phenomenon

The phenomenon of post-transcriptional buffering describes the fact that transcriptional variation tends to be buffered as the gene expression process progresses (Figure 14A). Over the past decades, this phenomenon has been repeatedly observed in several contexts and has thus emerged as a crucial determinant of the relationship between transcripts and proteins. Several early investigations of post-transcriptional buffering were made by comparing proteome and transcriptome changes associated with CNV in cancer cells [START_REF] Geiger | Proteomic Changes Resulting from Gene Copy Number Variations in Cancer Cells[END_REF][START_REF] Stingele | Global analysis of genome, transcriptome and proteome reveals the response to aneuploidy in human cells[END_REF].These studies, along with later investigations on larger sample in both human cancer cells and yeast (Dephoure et al., 2014;Gonçalves et al., 2017;[START_REF] Liu | Systematic proteome and proteostasis profiling in human Trisomy 21 fibroblast cells[END_REF]Zhang et al., 2014) highlighted that the proteome composition is not as sensitive to gene dosage variation as the transcriptome, which is known to typically reflects CNVs [START_REF] Fehrmann | Gene expression analysis identifies global gene dosage sensitivity in cancer[END_REF][START_REF] Schlattl | Relating CNVs to transcriptome data at fine resolution: assessment of the effect of variant size, type, and overlap with functional regions[END_REF]. Accordingly, early interspecies comparisons between proteome and transcriptome also highlighted that proteome variation were more constrained than transcriptome variation [START_REF] Khan | Primate Transcript and Protein Expression Levels Evolve Under Compensatory Selection Pressures[END_REF][START_REF] Laurent | Protein abundances are more conserved than mRNA abundances across diverse taxa[END_REF][START_REF] Schrimpf | Comparative Functional Analysis of the Caenorhabditis elegans and Drosophila melanogaster Proteomes[END_REF].

Taken together, these results indicate that protein abundance is more constrained and conserved than transcript abundance. Conceptually, this fits with the scheme that proteome variation will more directly affect the final phenotypic landscape compared to the transcriptome. Thus, changes in protein abundance are more likely to be deleterious and will be under stronger selective pressure. In addition to mRNA and protein comparisons, several ribosome profiling experiments showed that transcript abundance fluctuations between individuals, species or even different conditions are also buffered at the translational level (Artieri and Fraser, 2014;Blevins et al., 2019;McManus et al., 2014;Wang et al., 2015Wang et al., , 2020)). This highlights that posttranscriptional buffering is a multilayered phenomenon, suggesting that multiple mechanisms underlie the phenomenon.

Mechanisms underlying the post-transcriptional buffering

Although post-transcriptional buffering is frequently observed when comparing transcription, translation, and protein abundance together, the cellular mechanisms underlying this phenomenon remain elusive and poorly understood. Since post-transcriptional buffering is observed in both ribosome profiling and proteomic data, it is likely that multiple and distinct processes are involved.

Conceptually, post-transcriptional buffering requires some kind of feedback in which the cell detects and copes with expression fluctuations. Autoregulation may be a simple mechanism that allows this type of feedback (Figure 14B), as many transcription factors and RNA-binding proteins can inhibit their own expression [START_REF] Grönlund | Transcription factor binding kinetics constrain noise suppression via negative feedback[END_REF][START_REF] Müller-Mcnicoll | Auto-regulatory feedback by RNAbinding proteins[END_REF].

E3 ubiquitin ligases are also under strong autoregulation as they are able to recognize and target their excess proteins, thereby inducing proteasome-mediated degradation (de Bie and Ciechanover, 2011). However, autoregulation is limited to proteins that can influence their own expression or abundance (transcription factors, mRNA-binding proteins, ubiquitin ligases…)

and is obviously not sufficient to explain the extensive and frequently observed buffering.

Moreover, it is mechanistically difficult to imagine a global mechanism capable of detecting fluctuations in protein abundance and driving translation in response to these variations.

Accordingly, a previous study in aneuploid yeast has not found feedback mechanisms on protein synthesis in the case of expression variation for molecular complex-related proteins (Taggart and Li, 2018). Yet, the proteins with the most robust abundance to variation seem to be precisely molecular complex proteins (Dephoure et al., 2014;Gonçalves et al., 2017;Liu et al., 2017;[START_REF] Stingele | Global analysis of genome, transcriptome and proteome reveals the response to aneuploidy in human cells[END_REF]. Indeed, protein degradation of the unassociated complex components may play an important role in post-transcriptional buffering at the proteome level (Gonçalves et al., 2017;Juszkiewicz and Hegde, 2018;Taggart et al., 2020). The strong robustness of the complex-related protein is consistent with the numerous investigations highlighting that gene dosage imbalance can be highly deleterious (Deutschbauer et al., 2005;Morrill and Amon, 2019;Ohnuki and Ohya, 2018;Veitia and Potier, 2015). Again, these findings represent only a subset of genes. What other processes enable the cell to cope with expression variation at the proteome level is still unknown.

At the level of translation, knowledge about post-transcriptional buffering is even more scarce.

Although the phenomenon has been detected in several studies, no clear description of posttranscriptional buffering at this level has been made. The few insights into the mechanisms underlying translational buffering suggest that variations in mRNA abundance at the translation level tend to be attenuated by modulation of translation efficiency (i.e., the ratio between the quantification of translation for a gene and its mRNA abundance) (McManus et al., 2014).

Overall, although post-transcriptional buffering plays a major role in expression variation across multiple scales (intra-and interspecific, across conditions...) and levels (translation and protein abundance), it is still poorly characterized.

Overlap in the genetic origins of transcript and protein abundance

The final abundance of both mRNA and protein is the consequence of tight genetic regulation.

Both GWAS and linkage mapping have been used to explore the genetic origin of gene expression resulting in the discovery of thousands of eQTL and pQTL. The question of whether pQTL mirror eQTL has of course been explored on several occasions, but to date, no clear consensus has been reached. Due to the different correlation levels of the across-and withingene correlation (Buccitelli and Selbach, 2020), it is puzzling to answer to this question intuitively. On the one hand, the high degree of similarity between mRNA abundances within the sample (i.e., the across-gene correlation) suggests that protein abundance is largely the result of mRNA abundance and therefore the genetic regulation of the proteome could only result from that of the transcriptome. Conversely, the lower and still undetermined mRNA-protein correlation across samples (i.e., the within-gene correlation) could emphasize that the abundance variations between individuals are highly expression layer specific, and thus the genetic origins of the transcriptome and proteome should be as well.

A debated similarity

The simultaneous exploration of mRNA and protein abundance has several prerequisites in addition to those necessary for GWAS or linkage mapping studies. First, the conditions between transcriptomic and proteomic exploration need to be as similar as possible to ensure that the observed gene expression variation is not related to condition biases and is mostly due to genetic variation between individuals. Moreover, GWAS or linkage mapping exploration must be performed on the same cohort of individuals to obtain comparable genetic origins.

A significant number of studies have questioned the overlap between the genetic origin of a protein and its transcript abundance. In yeast, for example, two major studies focused on this issue (Albert et al., 2014;Foss et al., 2007). In both cases, these studies were performed on segregants from a cross between two isolates (a laboratory strain (BY) and a wine strain (RM11)). The early study found that the comparison between eQTL and pQTL was very modest: about 10% of the eQTL also affected the abundance of the related protein. This was consistent with later findings focusing on protein networks, which showed that protein coregulation across natural yeast isolates was mostly different from mRNA co-regulation [START_REF] Foss | Genetic Variation Shapes Protein Networks Mainly through Non-transcriptional Mechanisms[END_REF]. However, the more recent eQTL/pQTL exploration in yeast (Albert et al., 2014) showed very contrasting results: in this work, more than 60% of the eQTL had a corresponding pQTL. It is important to note that the two studies were technically different, as the latter was performed and based on a single-cell approach with quantification based on a GFP tag, whereas the previous study used LC-MS/MS quantification.

The similarity between eQTL and pQTL has also been investigated in the mouse model on several occasions (Chick et al., 2016;Ghazalpour et al., 2011). Here, the two studies are technically similar for the protein quantification (based on tandem mass-spectrometry) but slightly different for mRNA quantification (the latest used RNA-seq while the first one used RNA-microarray). It is also worth noting that in one case the eQTL/pQTL was performed on an inbred population (Ghazalpour et al., 2011) while in the other case, an outbred population was used (Chick et al., 2016). Again, there is a large discrepancy between the two studies: in the case of the inbred mouse set, approximately 5 to 6% of the eQTL have a corresponding pQTL, while in the case of the outbred mouse set, this value reaches 33%. Despite the drastically different conclusions of these studies, they agree on the fact that the local eQTL and pQTL tend to be more shared than the distant ones. For example, the overlap between the two types of QTL in Chick et al. (2016) is almost exclusively related to local QTL (1,392 local eQTL out of a total of 1,401 overlapping eQTL). Another contrasting result from these two studies needs to be emphasized: while the two studies had a similar within-gene correlation (between 0.25 and 0.30), suggesting that the mRNA-protein correlation across a population may not be an accurate predictor of the similarities between the genetic origins of mRNA and peptide abundance, one of them highlights that the genes with overlapping eQTL and pQTL tend to have higher within-gene correlations (Chick et al., 2016). Interestingly, a survey on human lymphoblastoid cell lines revealed a similarly high overlap: 33% of the eQTL replicated in the set of pQTL (Battle et al., 2015). Consistent with the eQTL/pQTL overlap, the correspondence between the regulatory hotspots affecting mRNA or protein abundance is often inconsistent across the studies (Albert et al., 2014;Foss et al., 2007), especially since pQTL or eQTL hotspots are not always detected (Ghazalpour et al., 2011).

Limitation of the explorations

The important difference in the overlap between eQTL and pQTL may be due to several reasons. First, the precision of quantification: in the early studies (Foss et al., 2007;Ghazalpour et al., 2011), the quantification of mRNA abundance was based on microarray technology, which is known to be less sensitive to subtle changes in mRNA abundance compared to RNA sequencing [START_REF] Mantione | Comparing Bioinformatic Gene Expression Profiling Methods: Microarray and RNA-Seq[END_REF]. However, the study with the higher overlap (≈60%) was also microarray-based (Albert et al., 2014;[START_REF] Smith | Gene-Environment Interaction in Yeast Gene Expression[END_REF]. Thus, quantification precision is unlikely to be the main reason for the contrasting results.

Another reason that could affect the reliability of the previous studies and therefore cause the observed discrepancies is the size of the populations in which mRNA and protein abundance was monitored. Overall, the more recent the study, the larger the cohort: the first studies on yeast and mice included 94 and 97 individuals, respectively (Foss et al., 2007;Ghazalpour et al., 2011) whereas the later studies included 114 and 192 individuals, respectively (Albert et al., 2014;Chick et al., 2016). The study on human lymphoblastoid cell lines focused on a total of 62 lines (Battle et al., 2015). Despite the tendency to study larger cohorts, these numbers are relatively small compared to the dimensionality of the problem, as only a very small fraction of the genetic diversity of the species is studied. Moreover, the genetic backgrounds of both the yeast and mouse strains in the aforementioned studies suffer from major limitations in terms of natural diversity: the yeast strains were for the most part generated from a simple cross between only two isolates (Albert et al., 2014;Foss et al., 2007), and the 97 and 192 mice are either inbred lines (Ghazalpour et al., 2011) or outbred lines derived from 8 inbred individuals (Chick et al., 2016). Again, this is a major limitation in terms of natural genetic diversity. It is therefore difficult to extend these results to larger scales. Ideally, a reliable approach to explore the similarities between the genetic origins of mRNA and protein abundance at the species level would be based on a larger scale exploration of gene expression, which is now possible as both RNA-seq and LC-MS/MS have been developed to reach population scales (Caudal et al., 2023;Messner et al., 2023Messner et al., , 2022;;The GTEx Consortium, 2015). Equally important, to ensure truly reliable genetic diversity, the population studied should consist of natural rather than constructed strains or isolates.

Saccharomyces cerevisiae, a powerful model to explore gene expression variation

As mentioned in the previous chapter, large-scale exploration of gene expression is a fundamental step in dissecting the genotype-phenotype relationship. Although technological advances have made population-scale studies more accessible, this type of exploration remains laborious. In this context, S. cerevisiae is a convenient and safe organism on which most molecular techniques have either been developed or adapted. This yeast is an ascomycetous fungus with a 12 Mb nuclear genome, distributed across 16 chromosomes resulting from an ancestral whole genome duplication likely caused by an hybridization event and forming an allopolyploid (Marcet-Houben and Gabaldón, 2015; [START_REF] Wolfe | Molecular evidence for an ancient duplication of the entire yeast genome[END_REF]. The genome is very compact compared to other eukaryotes as 70% of the genome corresponds to coding sequences and only 2% of the protein coding genes have an intron [START_REF] Hooks | Intron Evolution in Saccharomycetaceae[END_REF]. This yeast can be found is a very large diversity of natural and anthropized natural environment. The natural history of this species has obviously been strongly influenced by its extensive use in the context of anthropized fermented substrates [START_REF] De Guidi | Yeast domestication in fermented food and beverages: past research and new avenues[END_REF].

S. cerevisiae as a model to explore genome, transcriptome and proteome variations

A deeply characterized genome S. cerevisiae, which was the first eukaryote to be fully sequenced in 1996 [START_REF] Goffeau | Life with 6000 Genes[END_REF], has been a central model organism for biological science. There are more than 6,000 genes in the S. cerevisiae genome, although this number varies slightly between isolates (Peter et al., 2018). As of today (June 2023), complete genome sequencing of the S. cerevisiae genome with

Illumina sequencing has been conducted on more than 3,000 isolates from a vast diversity of geographical and ecological origins [START_REF] Basile | Large-scale sequencing and comparative analysis of oenological Saccharomyces cerevisiae strains supported by nanopore refinement of key genomes[END_REF][START_REF] Duan | The origin and adaptive evolution of domesticated populations of yeast from Far East Asia[END_REF][START_REF] Gallone | Domestication and Divergence of Saccharomyces cerevisiae Beer Yeasts[END_REF][START_REF] Lee | Extensive sampling of Saccharomyces cerevisiae in Taiwan reveals ecology and evolution of predomesticated lineages[END_REF]Peter et al., 2018;[START_REF] Strope | The 100-genomes strains, an S. cerevisiae resource that illuminates its natural phenotypic and genotypic variation and emergence as an opportunistic pathogen[END_REF]. Among these studies, and to date, the more complete genomic exploration of S. cerevisiae has been performed by fully sequencing the genomes of 1,011 natural isolates with Illumina technology (Peter et al., 2018). The studied population includes both wild and domesticated strains, with diverse ecological and geographical origins. Indeed, the strains were sampled from all 5 continents and come from different isolation sources such as clinical, wild (e.g., flower, soil, tree, water), wine, bread, bioethanol production. In this population, more than 1.6 million SNPs and 125,000 indels were detected, highlighting a high nucleotide diversity within this species, reaching up to 1.8% between the most distantly related isolates. The vast majority of SNPs are present at low frequency within the population, as 92% of these polymorphic positions have a MAF less than 5%. The genetic diversity observed within this population results from specific evolutionary events and determinants that have shaped the genomes of the strains during the evolution of the species. Though the construction of a neighbor-joining tree based on the complete dataset of bi-allelic SNPs, 26 subpopulations were identified, recapitulating for most of them the ecological origin of the isolates (Figure 15) and allowing a clear distinction between wild and domesticated subpopulations. This neighbor-joining tree constructed from the biallelic SNPs, highlights the subpopulations and their ecological origins (domesticated, wild, human, and unknown). Interestingly, this population also displays an important diversity in terms of CNVs. Indeed, each gene has a duplicated or deleted copy at least once across the 1,011 isolates, and CNVs were highlighted in each isolate. This study allowed to define the pangenome of the species, i.e., the complete set of genes found in the species based on this population, that is composed of 7,796 genes, as well as the core genome, i.e., the set of genes shared by all the isolates, that gathered 4,940 genes. In-between, the 2,856 remaining genes are considered as accessory, as only a fraction, that can be large or small, of the population carries them [START_REF] Li | The pan-genome of Saccharomyces cerevisiae[END_REF]Peter et al., 2018). Accessory genes are particularly interesting because they reveal much about the evolutionary history of S. cerevisiae. For example, many accessory genes originate from introgressions with the closely related species Saccharomyces paradoxus [START_REF] D'angiolo | A yeast living ancestor reveals the origin of genomic introgressions[END_REF]Peter et al., 2018). The uneven distribution of introgressed accessory genes across subpopulations allowed to trace different hybridization events between the two species. The Mexican agave, Alpechin, and French Guyana subpopulations have, on average, more introgressed genes than the other subpopulations. Horizontal gene transfer is also a mechanism that has led to the accumulation of several accessory genes, especially in wine isolates, where several genes coming from Zygosaccharomyces bailii and Torulaspora microellipsoides confer evolutionary advantage to the recipient isolates, especially in the winemaking environment (Marsit et al., 2015;Novo et al., 2009).

Another fundamental aspect of the S. cerevisiae genome is that the isolates show very different levels of ploidy. While the majority of natural isolates are in a diploid state, haploid or higher ploidy levels (3n, 4n and 5n) are common place (Peter et al., 2018). More specifically, some anthropized subpopulations, such as the beer ones, are enriched in polyploids. Among the 1,011 studied isolates, more than 200 isolates have at least one chromosome in an aneuploid state.

The aneuploidies are both related to chromosome gain and loss and are unevenly distributed across the genome: the smallest chromosomes, i.e., 1 st and 9 th , are preferentially affected.

Interestingly, despite the known preference for asexual reproduction in S. cerevisiae, about 63% of isolates (mainly domesticated isolates) show heterozygosity punctuated by loss of heterozygosity (LOH) regions of varying size, depending on the subpopulation. For example, the Sake subpopulation, that is mostly diploid, has an average of 80% of its genome affected by LOH.

Even if a large diversity of genetic variants was detected in this population, it is important to emphasize that the genotyping relied on a short-read sequencing approach (Peter et al., 2018), which missed a large fraction of SVs (translocations, inversions, insertions, deletions, etc.).

Since this type of variant is known to have a strong influence on yeast phenotype and gene expression [START_REF] Gorkovskiy | The Role of Structural Variation in Adaptation and Evolution of Yeast and Other Fungi[END_REF][START_REF] Hou | Chromosomal rearrangements as a major mechanism in the onset of reproductive isolation in Saccharomyces cerevisiae[END_REF]Zimmer et al., 2014), further genomic characterizations using long-read sequencing techniques at large scale will be an essential step to have a comprehensive view of genetic diversity among the S. cerevisiae species.

Population-scale gene expression exploration in S. cerevisiae

Because of its safe and easy manipulation, combined with rapid growth capabilities and a wide range of available technics for yeast, S. cerevisiae is indeed a very good model to study gene expression. As already mentioned, several studies of gene expression have been carried out using S. cerevisiae, for each step of gene expression (Albert et al., 2014;Artieri and Fraser, 2014;[START_REF] Brem | The landscape of genetic complexity across 5,700 gene expression traits in yeast[END_REF][START_REF] Foss | Genetic Variation Shapes Protein Networks Mainly through Non-transcriptional Mechanisms[END_REF]Gygi et al., 1999;[START_REF] Khan | Protein quantification across hundreds of experimental conditions[END_REF]McManus et al., 2014;[START_REF] Smith | Gene-Environment Interaction in Yeast Gene Expression[END_REF]. However, the large-scale exploration of gene expression is a very recent advance for this species (Caudal et al., 2023;Messner et al., 2023). The 1,011 population described above (Peter et al., 2018) was recently used for an extensive transcriptomic survey (Caudal et al., 2023). In this study, high-quality transcriptomes of 969 isolates were generated in a synthetic complete medium. The detected transcripts comprised 6,445 open reading frames (ORF), of which 4,977 belonged to the core genome and 1,468 to the accessory genome. This is one of the largest population transcriptome explorations to date.

Due to its very large scale, this study allowed for an in-depth characterization of the transcriptome variation between individuals. Surprisingly, the results showed that, overall, accessory genes have a very specific transcriptional behavior, being less expressed than the other genes, but more variable in expression across the individuals. However, this behavior was also related to the type of accessory genes. For example, there was a large variability in the mRNA abundance of accessory ORF originating from HGT evens across the population, while this was not the case for the ORF that were acquired through introgression events. In this particular case, it was even possible to compare the expression of the genes introgressed from S. paradoxus with the expression of their orthologs in the context of allelic heterozygosity (i.e., for an ORF, an isolate has one copy of the S. cerevisiae allele and one copy of the S. paradoxus allele). Using allele-specific expression, no difference was found between the expression of the S. paradoxus allele and the S. cerevisiae allele. By focusing on the subpopulation level, it was possible to detect specific differentially expressed genes (DEG) for most subpopulations (Caudal et al., 2023). The DEG were strongly associated with the environment from which the strains were isolated. For example, the GAL pathway was overexpressed in dairy fermentationrelated isolates, even in the presence of glucose. This type of metabolic switch represents a key adaptation to lactose-rich media [START_REF] Boocock | Ancient balancing selection maintains incompatible versions of the galactose pathway in yeast[END_REF][START_REF] Duan | Reverse Evolution of a Classic Gene Network in Yeast Offers a Competitive Advantage[END_REF], and highlights the central role of gene expression modulation in adaptation to anthropized processes. Finally, this investigation allowed for a precise association between mRNA abundance and genetic variants (SNPs and CNVs) present in this population. Using GWAS, it was possible to detect 7,273 SNP-eQTL and 2,197 CNV-eQTL affecting a total of 3,471 genes. Several fundamental aspects underlying the variation in mRNA abundance were observed with this GWAS. First, both cis and trans-eQTL were detected, with cis having a larger effect on their target genes. Regarding the CNV-eQTL, many of them were in fact related to the aneuploidies of the 1 st , 3 rd , 8 th , 9 th , and 11 th chromosomes. Interestingly, a difference was found between the effect of CNV-eQTL and SNP-eQTL, with the latter having a more important effect on mRNA abundance.

On the proteomic side, a recent study was carried out on the same set of isolates (Muenzner et al., 2022): the proteomes of 613 isolates grown on a synthetic minimal medium were accurately monitored using a high-throughput LC-MS/MS, resulting in the quantification of 1,563 proteins. The main focus of this study was to investigate the effect of aneuploidies on the proteome and to compare this with the transcriptome. Consistently with the post-transcriptional buffering phenomenon, a general dosage compensation was observed at the proteome level. In fact, the mRNA-abundance reflected more accurately the chromosomal imbalance than the protein abundance. In addition, the ubiquitin-proteasome system was likely a major process in variation buffering. Yet, due to the differences in the culture conditions, it is difficult to extend the comparison between the population's proteome and the transcriptome to other aspects of the gene expression population, as the medium difference could be a major confounding factor, especially when comparing GWAS results. Therefore, a large-scale and exact comparison between the population transcriptome and proteome is still lacking.

The S. cerevisiae domestication and its consequences on its evolutionnary history

The aforementioned investigations of both the genomes and transcriptomes of the S. cerevisiae population have emphasized the drastic impact of domestication on the S. cerevisiae evolutionnary history. As this species is likely to be the most widely used in food making because of its efficient fermentation capabilities, the domestication of S. cerevisiae has been extensively studied. However, due to the large number of strains that are used in different industrial contexts, the impact of domestication, particularly at the molecular level, have yet to be fully characterized.

Domestication and industrial use of S. cerevisiae

Evidence of deliberate fermentation has been found in early human history, dating back to prehistorical time [START_REF] Gallone | Domestication and Divergence of Saccharomyces cerevisiae Beer Yeasts[END_REF][START_REF] Mcgovern | Fermented beverages of preand proto-historic China[END_REF][START_REF] Michel | Chemical evidence for ancient beer[END_REF][START_REF] Samuel | Investigation of Ancient Egyptian Baking and Brewing Methods by Correlative Microscopy[END_REF]. The clear detection and association between yeast and fermentation came as modern science was in its early stage, when Louis Pasteur described the role of S. cerevisiae in alcoholic fermentation [START_REF] Pasteur | Nouveaux faits concernant l'histoire de la fermentation alcoolique[END_REF]. The case of the domestication of microorganisms is somewhat peculiar, since human selection, until the beginning of the last century, was mainly based on the indirect assessment of metabolic capacities, whereas the domestication and selection of animals or plants is based on visual and more quantifiable phenotypes (e.g., biomass production, size) [START_REF] De Guidi | Yeast domestication in fermented food and beverages: past research and new avenues[END_REF]. Because of this, it is expected that the domesticated strains will have very specific metabolic capacities depending on their isolation origin.

In the 1,011-population described above, the ecological backgrounds of the domesticated isolates are diverse (Figure 15, 16) (Peter et al., 2018). Alcoholic beverages are the main source of domesticated isolates. As observed through the structure of the population: the largest subpopulation is related to European wine isolates and 3 subpopulations mainly include beer isolates (African beer, mosaic beer, ale beer). In addition, 3 subpopulations include strains from other alcoholic beverages (sake, African palm wine and Mexican agave-related beverage).

Isolates related to food fermentation also grouped in specific subpopulations such as the ones related to cheese making (French dairy) while bakery related strain are mostly dispersed within the population. Finally, 35 domesticated isolates come from bioethanol production sites and are grouped in a specific subpopulation (Brazilian bioethanol). It is worth noting that, although being selected for fermentative purposes, the ecological niches of the domesticated isolates are very different. The microbial interactions, the temperature, the consistency of the medium and many other fundamental factors that are known to affect cell biology are specific to each ecological origin. Therefore, characterization of the molecular impact of domestication may require subpopulation-specific investigations. It is still unclear whether a single or multiple events led to the domestication of S. cerevisiae.

The general consensus states that the species originated from eastern China [START_REF] Duan | The origin and adaptive evolution of domesticated populations of yeast from Far East Asia[END_REF]Peter et al., 2018). However, while the 1,011-isolates collection supports the fact that multiple independent events shaped the domestication of S. cerevisiae (Peter et al., 2018), a study of Asian strains showing that domesticated isolates form two main groups (depending on the liquid or solid state on which they were isolated) suggests that domestication resulted from a single bottleneck event [START_REF] Duan | The origin and adaptive evolution of domesticated populations of yeast from Far East Asia[END_REF]. Regarding the Chinese domesticated isolates, their population structure indeed seems to support the bottleneck hypothesis [START_REF] Duan | The origin and adaptive evolution of domesticated populations of yeast from Far East Asia[END_REF].

However, the larger study of the 1,011 isolates includes a greater genetic diversity and may therefore be more reliable in accurately describing the domestication history of the species as a whole.

Human, clinical (138) Lab strain ( 2) Probiotic ( 2) Unknown ( 28) Flower ( 14)

Fruit ( 47)

Insect ( 20)

Nature ( 52) Soil ( 38)

Tree ( 64)

Water ( 19) Bakery ( 37) Beer ( 59) Bioethanol ( 27) Cider ( 17) Dairy ( 27) Distillery ( 29) Fermentation ( 36) Industrial (30) Palm wine (30) Sake ( 47)

Wine (248)
Phenotypic and molecular impact of the domestication Domestication usually results in a profound modification of the phenotypic landscape of a species. In the case of yeast, a wide range of trait variation can be observed when comparing wild and domesticated isolates. On the metabolic side, the fermentative capacity of the domesticated isolates are increased [START_REF] Bell | Comparison of fermentative capacities of industrial baking and wild-type yeasts of the species Saccharomyces cerevisiae in different sugar media[END_REF] and, more globally, there is a shift towards fermentation rather than respiration (Lahue et al., 2020). Interestingly, each domesticated subpopulation has undergone specific metabolic evolution due to the diversity of the ecological niches. As mentioned above, the diary subpopulation has a particular shift between glucose and galactose metabolism, which conferred the cell an increased ability to ferment lactose [START_REF] Boocock | Ancient balancing selection maintains incompatible versions of the galactose pathway in yeast[END_REF]Caudal et al., 2023;[START_REF] Duan | Reverse Evolution of a Classic Gene Network in Yeast Offers a Competitive Advantage[END_REF]. Beer isolates have on their side a better fitness when grown with maltose as a carbon source [START_REF] Gallone | Domestication and Divergence of Saccharomyces cerevisiae Beer Yeasts[END_REF].

However, the adaptation of carbon or energy metabolism is only a small part of the drastic changes that result from domestication. Several other cellular or molecular phenotypes have been associated with the domesticated isolates of S. cerevisiae. The general signatures (other than metabolic adaptation) of the domesticated isolate include improved osmotic stress tolerance and reduced sporulation (despite the higher proportion of heterozygotes) [START_REF] De Guidi | Yeast domestication in fermented food and beverages: past research and new avenues[END_REF]. The Sake isolates for instance went through morphological changes [START_REF] Ohnuki | Phenotypic Diagnosis of Lineage and Differentiation During Sake Yeast Breeding[END_REF] and are now highly resistant to high ethanol concentration [START_REF] Shiroma | Enhancement of ethanol fermentation in Saccharomyces cerevisiae sake yeast by disrupting mitophagy function[END_REF][START_REF] Watanabe | Ethanol fermentation driven by elevated expression of the G1 cyclin gene CLN3 in sake yeast[END_REF]. Stress resistance is also an important domestication trait in the wine subpopulation, as these isolates tend to be more resistant to the presence of copper or sulfite [START_REF] Brandolini | Saccharomyces cerevisiae wine strains differing in copper resistance exhibit different capability to reduce copper content in wine[END_REF]Yuasa et al., 2004). Interestingly, several SVs underlie this increase in stress resistance in wine isolates. For copper tolerance, this is associated to an increased copy number of the CUP1 gene [START_REF] Fogel | Tandem gene amplification mediates copper resistance in yeast[END_REF]Peter et al., 2018;[START_REF] Steenwyk | Copy Number Variation in Fungi and Its Implications for Wine Yeast Genetic Diversity and Adaptation[END_REF]. As mentioned above, the sulfite resistance is associated to large chromosomal rearrangements such as translocations and inversions (García-Ríos et al., 2019;Pérez-Ortín et al., 2002;Yuasa et al., 2004;Zimmer et al., 2014). Chromosomal alterations is also frequently observed in beer isolates, as many of them are polyploids (3n, 4n and 5n) [START_REF] Gallone | Domestication and Divergence of Saccharomyces cerevisiae Beer Yeasts[END_REF]Peter et al., 2018;[START_REF] Saada | Phased polyploid genomes provide deeper insight into the multiple origins of domesticated Saccharomyces cerevisiae beer yeasts[END_REF], which is interesting as high ploidy levels are also a hallmark of plant domestication [START_REF] Purugganan | The nature of selection during plant domestication[END_REF].

(Japan), we performed ribosome profiling and RNA sequencing on 8 natural isolates of S. cerevisiae. Our results covered a large number of genes, 3,755 in total. We found that the transcriptional and translational variations were associated with metabolism-related genes. We detected post-transcriptional buffering in our dataset and found that modulation of translational efficiency is an important mechanism underlying this phenomenon. Interestingly, essential genes, protein complex-related genes as well as less transcribed genes were preferentially affected by post-transcriptional buffering. In addition, we investigated the translation of a subset of the S. cerevisiae pangenome, the accessory genes, focusing on the introgressed and Horizontal Gene Transfers (HGT) ORFs. We found that the introgressed genes were translated similarly to their orthologous ORF in the other isolates, whereas the HGT ORFs showed a lower translation efficiency. Overall, our results provide insights into the mechanisms underlying post-transcriptional buffering at the translational level and its specificity. For example, the cellular systems that cope with complex imbalance toxicity could be one of the drivers of the phenomenon already at the translational level, as it is at the proteomic level.

We sought to extend these findings at the proteome level by measuring protein abundances in these 8 isolates, which is described in the second chapter (Figure 1A). We analyzed the aforementioned RNA sequencing and ribosome profiling data together with proteomic data obtained by an LC-MS/MS approach. This dataset was generated in collaboration with the Weizmann Institute of Science (Israel) and includes 3,635 proteins. The 3-layer quantification (transcriptome, translatome and proteome) was possible for a total of 2,840 genes. We found that protein abundance variations were also mainly associated with metabolism-or respirationrelated genes. Again, we found that post-transcriptional buffering was a major determinant of protein abundance variation across our isolate. We even observed that the more advanced the gene expression process, the stronger this phenomenon was. Therefore, protein abundance variations are more constrained, and different from what can be found at the transcriptome or translatome level. The difference in variation can also be observed by looking at the correlation between each of the gene expression steps. Despite being thought as a proxy for protein abundance, the ribosome profiling data was only slightly more correlated to protein data in comparison with RNA-seq data for the across-gene correlation, meaning that the proteome is slightly better reflected by the translatome than the transcriptome. When looking at the withingene correlation, both RNA-seq and ribosome profiling data had a mediocre correlation with the protein abundance, meaning that inter-individual proteome variation are barely captured by transcriptome and proteome. Taken together, this suggests that each gene expression layer may be subject to different evolutionary constraints. We sought to verify this by quantifying gene expression evolution at each step for more than 700 features characterizing each gene. We found that although each gene expression layer has some specificities, and that there are some general rules underlying gene expression evolution. For example, genes that are central to cellular networks (i.e., that interact a lot with other proteins or play a fundamental process for the cell) tend to have more constrained gene expression regulation, while genes related to metabolism tend to have faster gene expression evolution. This is strongly consistent with the previous findings on post-transcriptional buffering.

Finally, we explored the relationship between the transcriptome and proteome, this time at the population level, in order to get an accurate view of how the transcriptome and proteome vary at the species level (Figure 1B). In the third and final chapter, I describe the comparison between two large-scale surveys related to the transcriptomes and the proteomes of the 1,011 S. cerevisiae strains from our collection. While the transcriptome dataset was generated in our laboratory for a previous project, the proteome dataset was generated in collaboration with the Charité University of Medicine in Berlin and the Francis Crick Institute in London. The combined proteome and transcriptome dataset covers 629 genes from 889 isolates. As our data was one of the largest mRNA and protein abundance comparisons to date, we explored several gene expression phenomena that have never been studied at this scale in yeast, such as gene coexpression networks, post-transcriptional buffering, or domestication-related proteomic signatures. More importantly, we detailed the within-gene correlation in our population and showed that it was rather weak, around 0.16. Although the level of this correlation has been debated, this is significantly lower than what has been found in previous studies. Interestingly, the correlation level tended to be gene-dependent, with metabolism-and respiration-related genes showing high correlation levels, while ribosome-related genes tended to be uncorrelated or anticorrelated. The overall weak correlation between mRNA and protein abundance variation suggested that the genetic origins of the transcriptome and proteome are different. We performed GWAS to unravel the association between genome variation (taking both SNPs and CNVs into account) and proteome or transcriptome variation. We found that the overlap between eQTL and pQTL was very modest, especially when looking at the SNP-based GWAS.

Only 3.6% of SNP-QTLs were shared between transcriptome and proteome. The CNV-based GWAS, on the other hand, showed a higher similarity between the genetic origins of mRNA and protein abundance, but this was mainly related to the presence of large aneuploid segments on specific chromosomes, which have a strong impact on gene expression. No overlap between CNV-eQTL or pQTL nonrelated to aneuploidy was observed. Taken together, our results show that at the population level, the proteome and the transcriptome are two very distinct layers of gene expression, with very specific mechanisms underlying inter-individual variation. (A) A first project focuses on gene expression variation across 8 isolates for which RNA-seq, ribo-seq and LC-MS/MS data are available. In a first chapter I will focus on the comparison between transcriptional and translational variation and on post-transcriptional buffering at the translational level. In a second chapter, I will combine these datasets with new proteomic data from one of the 8 isolates and explore the evolutionary constraints on gene expression evolution across the gene expression process. (B) In a third and final chapter, I will present the population-level exploration of both transcriptome and proteome across 889 isolates. The main focus and findings will be related to within-gene mRNAprotein correlation and the genetic origins of transcript and protein abundances. steps. Some of these technologies include DNA high-throughput sequencing methods such as RNA sequencing (Hrdlickova et al., 2017), and mass spectrometry [START_REF] Lu | Absolute protein expression profiling estimates the relative contributions of transcriptional and translational regulation[END_REF], which enable a global description of transcriptomic and proteomic dynamics. By associating these data with genomic data, we can greatly improve our understanding of the mRNA and protein abundance regulation at the population level [START_REF] He | Genome-wide pQTL analysis of protein expression regulatory networks in the human liver[END_REF]Kita et al., 2017;Messner et al., 2022;[START_REF] Suhre | Genetics meets proteomics: perspectives for large population-based studies[END_REF]The GTEx Consortium, 2015). Furthermore, the core of the translational process can be precisely dissected with the development of ribosome profiling (or Ribo-Seq) (Ingolia et al., 2019(Ingolia et al., , 2009)). This strategy relies on the sequencing of mRNA fragments covered by the ribosomes during the translation process, revealing which parts of the transcriptome are actively being translated. This method can quantify translation in mRNAwise (number of fragments of the corresponding mRNA) and also the behavior of the ribosomes at the given codon (density of the ribosomes along the mRNA) [START_REF] Ingolia | Ribosome profiling: new views of translation, from single codons to genome scale[END_REF].

The budding yeast Saccharomyces cerevisiae has been a powerful model for ribosome profiling experiments, as this technique was developed on this organism (Ingolia et al., 2009).

Translational variation in yeast has been explored with ribosome profiling on several occasions (Albert et al., 2014;Artieri and Fraser, 2014;Blevins et al., 2021Blevins et al., , 2019;;McManus et al., 2014;Wang et al., 2015). Interestingly, several of these studies highlighted that the transcriptional variations tended to be buffered when looking at the translational variations (Artieri and Fraser, 2014;Blevins et al., 2019;McManus et al., 2014;Wang et al., 2015). This phenomenon is known as post-transcriptional buffering and has also been observed when comparing transcriptomic and proteomic datasets (Gonçalves et al., 2017;Kustatscher et al., 2017).

However, despite recurring observations, the mechanisms underlying this post-transcriptional buffering are still poorly understood. Moreover, while transcription and protein abundance have been extensively monitored, translation itself has been considerably less studied, and no clear description of the phenomenon has been made at this level. More globally, translational variation remains largely unexplored, and several known sources of expression variation, such as accessory ORF (open reading frames), have yet to be investigated at the translational layer.

Here, we conducted ribosome profiling and RNA sequencing in the same conditions on eight S. cerevisiae natural isolates coming from very diverse ecological environments and being genetically different (Peter et al., 2018). We first compared the transcriptional and translational variations, and found that they had similar functional patterns. Metabolism-related genes tended to be more variable across the eight isolates while essential genes and genes involved in molecular complexes had more conserved transcription and translation regulation.

Interestingly, we found that the transcriptional profiles were less correlated to each other compared to the translational profiles. Accordingly, Euclidean distances and expression variations (quantified using the absolute log2 transformed foldchanges for each gene in each isolate pairwise comparisons) were approximately 10% higher in the transcriptomic data, indicating that post-transcriptional buffering is a strong determinant of the translational variations. More importantly, we found that this phenomenon has a specific signature in terms of affected genes. We observed that essential genes and protein complex-related genes as well as lowly transcribed genes tended to be preferentially buffered. Furthermore, we investigated the transcription and translation of accessory open reading frames (ORFs) present in the eight isolates, particularly those acquired through introgression or horizontal gene transfer (HGT) events. We observed that introgression-related ORFs were similarly transcribed and translated compared to their orthologs, while HGT-related ORFs displayed a significantly lower translation efficiency than the rest of genes. Together, our results provide an overview of translational variations as well as an accurate description of post-transcriptional buffering.

Results

Ribosome profiling and RNA sequencing across eight natural isolates

We performed both RNA sequencing (RNA-seq) and ribosome profiling (Ribo-seq) on eight genetically diverse S. cerevisiae isolates (Table S1), which were cultivated and harvested in the exact same condition. These isolates were selected to represent the genetic diversity of the species (Figure S1) and were grown on a synthetic complete medium. All the genomes of the isolates were all previously sequenced (Peter et al., 2018), and in addition to their very different genetic backgrounds, they also came from very diverse environments. After TPM (transcript per million) normalization of RNA-seq and Ribo-seq raw counts (see Methods), we computed a translation efficiency (TE) value of each gene in each isolate by dividing its Ribo-seq TPM value by its RNA-seq TPM value. In total, we analyzed 3,755 genes and our results showed a strong correlation between RNA-seq and Ribo-seq data (Spearman correlation test between 0.769 and 0.867), highlighting the relationship between transcription and translation. To gain a global view of intraspecific variation, we performed pairwise Spearman correlation tests on the two datasets for each strain (Figure 1A). The RNA-seq and Ribo-seq correlation matrices showed similar patterns, indicating that transcriptional variations were largely reflected at the translational level. However, the correlation coefficients were generally higher in the Ribo-seq matrix than the average value in the RNA-seq matrix (Figure 1B), suggesting that translational profiles were more similar than transcriptional profiles.

Next, we sought to identify genes that did not follow these correlation trends in RNA-seq and Ribo-seq data to detect genes with variable regulation of transcription and translation. To achieve this, we used a combination of two different methods to make a pairwise comparison of all isolates. The first method was a Mahalanobis distance calculation to detect outliers in the pairwise comparison (See Methods, Figure S2A) (Ho et al., 2018). The second method relied on the selection of genes displaying residual (obtained from a linear regression model computed on the isolate pairwise comparison), which are in 2.5% highest or 2.5% lowest residual quantiles. We selected the genes that overlapped between the two methods (Figure S2B) and identified a total of 357 genes in Ribo-seq data and 352 genes for RNA-seq, with 179 overlaps between the two expression layers (Figure S3A andB). These genes are later mentioned as "variable genes". Most of these genes were only detected once in the 28 pairwise comparisons (128 out of 357 Ribo-seq variables genes, 118 out of 352 RNA-seq variables genes, Figure 1C-D). The number of variable genes detected in the pairwise comparisons against the 7 other strains ranged from 125 to 204 for Ribo-seq (median=172) and from 128 and 204 to RNA-seq (median=177.5).

To investigate functional enrichment among these identified genes, we conducted a Gene Ontology analysis (GO) (Ashburner et al., 2000;Gene Ontology Consortium, 2021) on both Ribo-seq and RNA-seq variable genes, as well as on the 179 overlapping gene set. The RNAseq variable genes yielded a relatively low number of terms compared to the other dataset (20 terms, Table S2). However, the majority of its terms were detected in the Ribo-seq variable genes, suggesting that despites the difference between the two genes groups, the functions of the genes are mostly shared. We also found that 63 terms were shared between the Ribo-seq variable genes (out of 79 terms, Table S3) and the overlapping variable genes (out of 86 terms, Table S4). Many of the shared results were related to metabolism terms (such as "glycolytic process", "pentose-phosphate shunt" and "glucose metabolic process"). This observation may be explained by the fact that the eight isolates used in this study were obtained from distinct environments (Table S1) and might have adapted their regulation of several metabolic functions to different trophic conditions.

Interestingly, our analysis revealed that the GO term displaying the lowest p-value was the depletion of the "protein-containing complex" term (Ribo-seq variable genes: 1.08e-17, RNAseq variable genes: 7.55e-09, overlapping variable genes: 6.35e-18), suggesting that genes encoding protein involved in protein complexes are underrepresented in the variable gene sets.

To confirm this observation, we checked if genes previously annotated as related to protein complexes (Pu et al., 2009) were indeed significantly depleted from the variable genes set.

Similarly, we explored if genes characterized as essential (Dowell et al., 2010;Giaever et al., 2002) were enriched or depleted in the variable genes since these characteristics (essentiality and being involved in protein complexes) have been cited in previous studies exploring evolutionary constraints (Morrill and Amon, 2019;[START_REF] Pál | An integrated view of protein evolution[END_REF][START_REF] Rancati | Emerging and evolving concepts in gene essentiality[END_REF]. Using

Fisher's exact test (FET), we found that essential genes were strongly depleted in the variable gene set (RNA-seq: odds ratio = 0.40, p-value = 1.55e-10; Ribo-seq: odds ratio = 0.22, p-value = 8.89e-16). The results were very similar when examining protein complex-related genes (RNA-seq: odds ratio = 0.24, p-value = 1.81e-28; Ribo-seq: odds ratio = 0.17, p-value = 1.79e-38). Interestingly, the depletion was lower with the Ribo-seq variable genes for both essential genes and protein complex-related genes, suggesting that these genes were less likely to exhibit variable regulation at the translational level compared to the transcriptional level.

Together, our results highlight that expression variation is unequal among the genes. While metabolism related genes display important transcriptional and translational variations among the eight isolates, essential genes and protein complex related genes are related to a lower expression variation.

Post-transcriptional buffering at the translation level across isolates

Several results suggest lower variability at the translational level compared to the transcriptomic level, such as the higher correlation in the Ribo-seq dataset (Figure 1B). To confirm these observations, we first computed the Euclidean distances and checked if the distances between each profile were higher in RNA-seq or Ribo-seq log10 data. We found that the distance between the strains were approximately 10% higher in RNA-seq data (t-test, p-value = 0.036) (Figure 2A), suggesting that the transcriptional variations tended to be higher than the translational ones. Consistently, the expression variance of each gene across the 8 isolates was significantly higher at the transcription level (Figure S4).

We also quantified gene expression variation in the two expression layers by computing the absolute value of the log2 transformed foldchange (|log2(FC|, see Methods) for each gene in the pairwise comparisons of the isolates (n = 105,140). The higher this value will be, the stronger the transcriptional or translational variation will be. We found that |log2(FC| was approximately 10% higher in RNA-seq data (mean = 0.975) compared to Ribo-seq data (mean = 0.872, Figure 2B), suggesting that overall, the transcriptional variations tended to be buffered.

These results are consistent with the phenomenon of post transcriptional buffering that has been observed previously (Artieri and Fraser, 2014;Blevins et al., 2019;McManus et al., 2014;Wang et al., 2015). In their study, they showed that the buffering of transcriptional variation may be linked to modification of translational efficiency. We sought to confirm this by comparing the RNA-seq, Ribo-seq and TE log2(FC) values in each pairwise comparison. We observed that RNA-seq log2(FC) values were strongly anti-correlated with TE log2(FC) values (Figure 2C, S5A, Table S5). Conversely, we found that the comparisons between RNA-seq and Ribo-seq log2(FC) always showed a positive correlation (Figure 2D, S5B, Table S5). Finally, we found no correlation between Ribo-seq and TE log2(FC) (Figure S5C, Table S5). Together, our results suggest that even though transcriptional and translation variations are similar in direction, their strength are strongly attenuated at the translation level by an opposite change in translation efficiency.

Signature of the post-transcriptional buffering at the translation level

Despite several observations of post-transcriptional buffering, this phenomenon remains largely unknown, especially at the functional level (Artieri and Fraser, 2014;Blevins et al., 2019;McManus et al., 2014;Wang et al., 2015). We sought to further characterize the general rules underlying this phenomenon by looking for genes that would be preferentially affected by the post-transcriptional buffering at the translational level.

With that in mind, we split each of the RNA-seq vs TE log2(FC) pairwise comparisons in orthogonal spaces using a TE and RNA-seq fold-change threshold of 1.5 (Figure 3A). This method made it possible to distinguish two categories of genes. First, we detected genes with transcriptional variation using the RNA-seq log2(FC) threshold. Second, using the TE log2(FC) threshold, we were able to identify genes with buffered variation from those not affected by the buffering (Figure 3A). The number of genes affected by post-transcriptional buffering ranged from 854 to 1,319 across the pairwise comparisons, with a mean value of 1,051 per comparison (Figure S6). No genes were buffered in all pairwise comparisons, and the proportion of buffered genes among the genes with transcriptional variation averaged 46.8%, ranging from 37.2% to 58.8% (Figure S7).

We then selected genes whose transcriptional variation was recurrently buffered or unbuffered (see Methods). We detected 361 and 507 genes whose transcriptional variation was recurrently buffered and unbuffered, respectively. We searched for a functional signature among the buffered genes but found no significant enrichment using GO analysis. We then focused on the content of essential (Dowell et al., 2010;Giaever et al., 2002) and protein complex-related genes (Pu et al., 2009) in recurrently buffered genes in comparison to recurrently unbuffered genes. The proportion of essential genes was significantly higher in the buffered gene set compared to the unbuffered genes (Figure 3B). Similarly, protein complex-related genes were also in higher proportion among the buffered gene (Figure 3C). Together, these results support that the genes are unequally affected by the post-transcriptional buffering at the translational level, with essential genes and protein complex-related genes preferentially buffered.

Interestingly, essentiality and protein-protein interactions are features debated for their influence on protein sequence evolution [START_REF] Fraser | Evolutionary rate in the protein interaction network[END_REF][START_REF] Pál | An integrated view of protein evolution[END_REF][START_REF] Rancati | Emerging and evolving concepts in gene essentiality[END_REF][START_REF] Zhang | Determinants of the rate of protein sequence evolution[END_REF]. Another major determinant of sequence evolution is gene expression level, as highly expressed proteins are known to be more conserved [START_REF] Drummond | A Single Determinant Dominates the Rate of Yeast Protein Evolution[END_REF][START_REF] Drummond | Why highly expressed proteins evolve slowly[END_REF][START_REF] Rocha | The quest for the universals of protein evolution[END_REF][START_REF] Zhang | Determinants of the rate of protein sequence evolution[END_REF]. We therefore sought to see if the transcription level was also involved in preferential buffering. Surprisingly, we found that recurrently buffered genes tended to be less transcribed (median TPM = 49.8) than recurrently unbuffered genes (median TPM = 230.1, Figure 3D). The same results were observed when comparing the recurrently buffered genes to the rest of the genes (Figure S8A). The results were also similar when looking at the average expression levels of buffered and not buffered genes in the pairwise comparisons (Figure S9). Together, these results highlight that transcription level is also a determinant of the phenomenon of post-transcriptional buffering, since buffered genes tend to be less transcribed.

We sought to confirm these results by exploring the codon usage bias of the buffered genes since this feature is known to be related to expression level [START_REF] Coghlan | Relationship of codon bias to mRNA concentration and protein length in Saccharomyces cerevisiae[END_REF][START_REF] Plotkin | Synonymous but not the same: the causes and consequences of codon bias[END_REF] and essentiality [START_REF] Dilucca | Codon Bias Patterns of E. coli's Interacting Proteins[END_REF]. We computed a codon usage bias index for each gene using tRNA Adaptation Index (tAI) [START_REF] Dos Reis | Solving the riddle of codon usage preferences: a test for translational selection[END_REF][START_REF] Dos Reis | Unexpected correlations between gene expression and codon usage bias from microarray data for the whole Escherichia coli K-12 genome[END_REF], and we confirmed that this index correlated with RNA-seq or Ribo-seq data (Figure S10). We then compared the tAI values of the buffered group to those of the unbuffered group. We observed a significantly lower tAI in the buffered group (median tAI = 0.34) compared to the unbuffered group (median tAI = 0.38) (Figure 3E). Similar results were observed when comparing the recurrently buffered genes with the rest of the genes (Figure S8B). This observation supported the previous results of expression level difference between the two groups.

Overall, these results highlight the fact that the phenomenon of post-transcriptional buffering preferentially affects essential, protein complex-related genes, or genes with lower transcription levels and therefore has a specific signature. This behavior toward some specific categories of genes has never been shown before and highlights evolutionary constraints affecting the translational regulation of these genes.

Transcription and translation variation of accessory genes

Recent advances in S. cerevisiae population genomics have highlighted the presence of more than 1,700 variable ORFs (accessory genes) in this species (Peter et al., 2018). Our translation exploration across multiple individuals from very different genetic and environmental origins is a unique opportunity to explore the translation of such ORFs, which was omitted in previous yeast ribosome profiling investigations. In our eight isolates, the number of these ORFs varies between 63 to 215, corresponding to a total of 446 unique accessory ORFs (median = 94 accessory ORFs per strain) (Figure S11), but depending on the isolate, between 36% and 72% of them were expressed (Figure S11). We observed that the unexpressed ORFs tended to be smaller (on average 148.81 bp) than the expressed ones (on average 365.2 bp, Figure S12). Overall, our eight isolates displayed variable profiles in terms of accessory ORFs origins (Figure S11). Two strains differed notably to the others in their compositions: the CPI isolate due to a very high number of ORFs acquired by introgression and the BPL isolate due to genes acquired through horizontal gene transfer (HGT).

Regarding the CPI isolate, this strain was originally isolated in Mexico and has been described

as part of the "Mexican agave clade" (Peter et al., 2018). This subpopulation has a high number of introgressed ORFs coming from the yeast Saccharomyces paradoxus (median = 161 ORFs per strain vs 25.75 in the overall population). The CPI isolate had 87 expressed ORFs coming from introgression events, and 45 ORFs had known orthologs in S288C. In order to strictly explore the impact on transcription and translation, we focused on 18 out of the 45 ORFs that were homozygous for the S. paradoxus allele in the CPI isolate and found no expression difference between these ORFs and their orthologs in the 7 other strains (again at the transcriptional and translational level) (Figure 4A,B). These results imply that the transcriptional and translational regulation of ORFs acquired by introgressions from S.

paradoxus is similar to their regulation of their orthologs.

We then focused on the expression of the 16 accessory ORFs coming from HGT in the case of the BPL isolate (Table S7). This strain has already been described as part of a wine subpopulation (Peter et al., 2018) and the occurrence of HGT events in this type of strain has already been observed (Marsit et al., 2015;Novo et al., 2009). Briefly, the coexistence of S.

cerevisiae with other yeast species in the wine environment led to gene transfer that can confer evolutionary advantage in the winemaking environment. We compare the expression of these ORFs to other genes in the BPL isolate (Figure S13). Surprisingly, HGT ORFs almost all showed lower Ribo-seq values than RNA-seq, suggesting a lower TE than the rest of the genes.

We then compared HGT ORFs with the rest of the genes and found a significant difference where HGT ORFs were less efficiently translated than other genes in the BPL isolate (mean HGT TE = 0.545, mean HGT other genes = 0.865) (Figure 4C). This observation clearly shows that HGT-related ORFs exhibited significantly lower translation efficiency compared to the rest of the genome.

Discussion

Translational variation is a major determinant of the transcriptome-proteome relationship, and therefore plays a central role in the phenotypic diversity observed in natural populations.

However, the translation process itself remains largely unexplored, and the central mechanisms driving translation variations are still poorly understood. In this study, we have precisely monitored translation variations across natural isolates of S. cerevisiae using ribosome profiling.

Gene expression is known to differ across individual. We observed that the translational and transcriptional variations were functionally similar, with metabolism-related genes displaying the greatest variation while the translation regulation of essential genes and protein complexrelated genes were more conserved. These results are consistent with recent large-scale exploration of mRNA abundance (Caudal et al., Submitted) and highlight that gene expression plasticity might be driven by the metabolism preferences between isolates coming from very different environments (Hodgins-Davis et al., 2012). Conversly, expression variation of genes with central and essential functions is likely to be deleterious and tend to be therefore more conserved [START_REF] Fraser | Noise Minimization in Eukaryotic Gene Expression[END_REF].

Our dataset also allowed to get better insight into the phenomenon of post-transcriptional buffering at the translational level (Artieri and Fraser, 2014;Blevins et al., 2019;McManus et al., 2014;Wang et al., 2015). Using |log2(FC)|, we quantitatively measured gene expression variation across both layers of expression (transcription and translation), and found that the median transcriptional value |log2(FC)| was 10% higher than the median translational |log2(FC)|. Together with the differences in gene expression variance and correlations coefficients between isolates, this clearly indicates that post-transcriptional buffering is detected in our dataset. As previously suggested (McManus et al., 2014), we observed that TE modulation plays a central role in compensating for variations in mRNA abundance.

Interestingly, we found that the post-transcriptional buffering has a specific signature. It preferentially affects genes that are essential or related to protein complexes, as well as genes with low transcript levels. Several reasons could underlie this preferential buffering. For complex-forming protein, it is well established that the imbalance of complex components can be deleterious for (Deutschbauer et al., 2005;Ohnuki and Ohya, 2018;Veitia and Potier, 2015), partly for stoichiometric reasons (Morrill and Amon, 2019). More generally, protein complexrelated genes are known to have stronger regulatory control at the protein level rather than mRNA level (Jüschke et al., 2013) and programmed translation of complex components precisely proportional to stoichiometry was not only found in yeast (Taggart and Li, 2018), but also in bacteria and plants (Chotewutmontri and Barkan, 2016;Li et al., 2014;Lukoszek et al., 2016;Trösch et al., 2018). Essential genes are known to carry a central and highly conserved function in the cell [START_REF] Coghlan | Relationship of codon bias to mRNA concentration and protein length in Saccharomyces cerevisiae[END_REF], which can lead to higher constraints on gene expression evolution (Wang et al., 2020). However, the conditional nature of essentiality [START_REF] Larrimore | The conditional nature of gene essentiality[END_REF][START_REF] Papp | Metabolic network analysis of the causes and evolution of enzyme dispensability in yeast[END_REF]) and the ongoing debate on the importance of essentiality on evolutionary constraint [START_REF] Pál | An integrated view of protein evolution[END_REF][START_REF] Rancati | Emerging and evolving concepts in gene essentiality[END_REF] suggest that the link between expression conservation and essentiality remains unclear. Regarding the fact that buffered genes tend to be less transcribed than other genes, this is surprising since very abundant proteins are usually well conserved [START_REF] Drummond | A Single Determinant Dominates the Rate of Yeast Protein Evolution[END_REF][START_REF] Drummond | Why highly expressed proteins evolve slowly[END_REF][START_REF] Rocha | The quest for the universals of protein evolution[END_REF][START_REF] Zhang | Determinants of the rate of protein sequence evolution[END_REF]. Thus, we would have expected that the regulation of highly expressed genes would also be highly conserved as well, and therefore preferentially affected by the phenomenon of post-transcriptional buffering.

Finally, working with genetically distinct natural isolates allowed to explore the translation of a part of the S. cerevisiae accessory genome, which is something that has been barely investigated so far. Accessory genes had very different translation dynamics depending on their origins of acquisition. While introgression-related ORFs displayed similar levels of translation compared to their orthologs, HGT-related ORF were less translated, resulting in low translation efficiency. These results on the pangenome nevertheless remain limited due to the low representation of the entire pangenome of S. cerevisiae (Peter et al., 2018). More generally, a broader view of S. cerevisiae population translation would improve our understanding of the post-transcriptional buffering phenomenon.

Overall, our results highlight the importance of the post-transcriptional buffering at the translation level, as well as its specific signature. Moreover, they give one of the first insight into the translation dynamics of a specific part of the genome such as accessory genes.

Materials and methods

Strain, culture, and flash freezing

The complete list of isolates used in this study is available and described in Table S1. The strains were grown in liquid SC medium (Yeast Nitrogen Base with ammonium sulfate 6.7 g.l -1 , MPbio, OH, USA; amino acid mixture 2 g.l -1 , MPbio; glucose 20 g.l -1 , Euromedex, France).

The culture was maintained until the strains reached their growth mid log phase using an optical plate reader (Tecan infinite F200 pro). The cells were then filtered using 0.45 µm MCE membrane (Merk Millipore, France). The filters were then plunged into a 50 mL tube containing liquid nitrogen and stocked in a -80°C freezer before being used for ribosome profiling and RNA sequencing experiment.

Ribosome profiling and RNA sequencing

The library preparation for ribosome profiling was performed as previously described with modifications [START_REF] Mcglincy | Transcriptome-wide measurement of translation by ribosome profiling[END_REF][START_REF] Mito | Protocol for Disome Profiling to Survey Ribosome Collision in Humans and Zebrafish[END_REF]. Cells on the filters were mixed with frozen droplets of 600 µL lysis buffer (20 mM Tris-HCl pH 7.5, 150 mM NaCl, 5 mM MgCl2, 1 mM dithiothreitol, 100 µg/mL cycloheximide, and 1% Triton X-100) and crushed using Multi-beads Shocker (Yasui Kikai, Japan). Lysate containing 20 µg of total RNA was digested with 10 U of RNase I (Lucigen, WI, USA) for 45 min at 25ºC. Ribosomes were precipitated by sucrose cushion and ultracentrifugation, suspended into EDTA lysis buffer, in which 5 mM MgCl2 in lysis buffer was substituted with 5 mM ethylenediaminetetraacetic acid (EDTA), and transferred to Amicon Ultra-0.5 Ultracel-100 (Merck Millipore, MA, USA) to separate footprints from ribosome subunits (the details will be described elsewhere [M.M. and S.I., unpublished data]). RNAs ranging from 17 to 34 nt were excised from a polyacrylamide TBE-Urea gel. The rRNA was depleted using RiboMinus Transcriptome Isolation Kit (yeast) (Thermo Fisher Scientific, MA, USA).

For RNA-seq, total RNA was purified using TRIzol LS reagent (Thermo Fisher Scientific) and Direct-zol RNA Miniprep Kit (Zymo research, CA, USA). Following the removal of rRNA by RiboMinus Transcriptome Isolation Kit (yeast), the sequencing library was prepared with TruSeq Stranded mRNA Library Prep Kit (Illumina, CA, USA). The ribosome profiling and RNA-Seq libraries were sequenced on a HiSeq 4000 platform (Illumina) with a single-end 50 bp.

Sequence data alignment, quantification, and normalization

Alignment and quantification of ribosome profiling and RNA-Seq data were performed as previously described with modifications [START_REF] Mcglincy | Transcriptome-wide measurement of translation by ribosome profiling[END_REF]. After the removal of the linker sequence and the splitting based on sample barcode, we removed reads that mapped to non-coding RNA (ncRNA) sequences using STAR 2.7.0a [START_REF] Dobin | STAR: ultrafast universal RNA-seq aligner[END_REF]. Despite the fact that we used an rRNA depletion method for both RNA-seq and Ribo-seq, our libraries were highly contaminated with ncRNA (Table S7), resulting in a relatively low reads number input for the alignment: between 92,676 and 350,221 reads for RNA-seq and between 248,228 and 1,049,718 reads for Ribo-seq. Remaining reads were aligned to the S288C S. cerevisiae genome using STAR 2.7.0a [START_REF] Dobin | STAR: ultrafast universal RNA-seq aligner[END_REF]. For the analysis of the accessory ORFs (open reading frames), the reads were also aligned to all the ORF detected in the pangenome of S. cerevisiae (Peter et al., 2018). The A-site offsets of ribosome footprints were determined according to the location of the 5′ end of reads mapped to start codons. For RNA-seq, offsets were set to 15 for all mRNA fragments. Reads corresponding to the first and last five codons of each coding sequence (CDS) were excluded from the analysis. For calculation of transcript per million (TPM) values for each CDS, we normalized read counts by CDS length minus 10 and adjusted sum of all normalized values to one million. The custom scripts will be available upon requests.

We finally calculated a translation efficiency (TE) value by dividing the Ribo-seq TPM value by the RNA-seq TPM value for each gene in each isolate.

Expression variation analysis

The TPM normalized datasets for Ribo-seq and RNA-seq were log10-transformed. A first general overview of the strain variation was obtained using Spearman correlation test. Next, on each strains vs. strain pairwise comparisons, we applied Mahalanobis distance (Ho et al., 2018) using the check_outlier() function (from the R package "performance"), where the genes with a distance higher than 10.59 (χ 2 distribution, with a 0.005 alpha level and 2 degrees of freedom) were selected. These genes were then filtered using a linear model on the strains vs. strain pairwise comparison: the residuals coming from the linear model were used to select the genes

Introduction

Gene expression is the main driver of the relationship between the phenotypic landscape observed in natural populations and the genetic mechanisms underlying this large phenotypic diversity. Modification of gene expression, and more specifically protein abundance, is a known source of phenotypic variation (Albert and Kruglyak, 2015;Maurano et al., 2012). In recent decades, the increased accessibility of various gene expression quantification methods, such as RNA sequencing or LC-MS/MS, has allowed several studies to be conducted to examine mRNA or protein levels across healthy or natural individuals (Battle et al., 2015;Ferkingstad et al., 2021;The GTEx Consortium, 2015). In addition, the development of techniques such as ribosome profiling has led to the precise study of the translation process (Ingolia et al., 2009).

This technique is based on the sequencing of mRNA fragments protected by ribosomes during translation and is considered a better proxy for protein abundance than transcript abundance (Brar and Weissman, 2015;Ingolia, 2010). While several pathologies are known to be associated with changes in gene expression (Corbett, 2018;Lee and Young, 2013), physiological and benign changes in either transcript or protein abundance are repeatedly observed across individuals (Battle et al., 2015;[START_REF] He | Genome-wide pQTL analysis of protein expression regulatory networks in the human liver[END_REF][START_REF] Niu | Plasma Proteome Variation and its Genetic Determinants in Children and Adolescents[END_REF]The GTEx Consortium, 2015).

However, it is well known that the natural variation in gene expression between individuals differs depending on the level at which gene expression is considered. In fact, several studies have highlighted that expression variation tends to decrease as the gene expression process progresses. This phenomenon, called post-transcriptional buffering, has been observed when both proteomic and ribosome profiling data are compared with transcriptomic data (Artieri and Fraser, 2014;Blevins et al., 2019;Dephoure et al., 2014;Gonçalves et al., 2017;McManus et al., 2014;Wang et al., 2020). The mechanisms underlying this phenomenon remain unclear.

Several cellular processes such as autoregulation or stoichiometry control are thought to be involved in such a phenomenon, but they are not sufficient to fully explain the extent of posttranscriptional buffering (Buccitelli and Selbach, 2020). Overall, this phenomenon suggests that protein abundance is more conserved than transcript abundance and is therefore subject to different evolutionary constraints. Yet, the comparison between the determinants of gene expression evolution between each layer of expression is still lacking, mainly because few studies have focused on all expression levels together.

To explore what are the determinants of gene expression evolution across individuals, we precisely quantified the protein abundance of 8 natural isolates of Saccharomyces cerevisiae previously used for RNA sequencing and ribosome profiling experiments (Teyssonniere et al., Submitted). We combined these 3 datasets and observed that the variation in protein abundance was mainly related to metabolic genes, which is consistent with the variation observed at the transcriptome and translatome levels (Teyssonniere et al., Submitted). Consistent with the previous observation of post-transcriptional buffering, we observed that protein abundance was less variable than both transcript abundance and ribosome-protected mRNA fragments (RPF)

abundance. Interestingly, we observed that the ribosome profiling data was only slightly more correlated with protein abundance than the transcriptomic data when looking at the mRNAprotein correlation of each individual. In addition, when looking at the gene-wise mRNAprotein correlation (i.e. the correlation between the abundance variation across the isolate), RPF

was not more correlated with protein abundance than mRNA abundance. Finally, we explored the determinants of gene expression evolution at each level (mRNA, RPF and protein abundance) and found that if the proteome tends to have specific constraints, several general rules shaped the evolution of each gene expression level. Consistent with the results above, metabolism-related genes tended to have the faster evolutionary rate across all expression layers. Conversely, highly interacting genes or genes involved in central cellular processes were associated with higher evolutionary constraints in each layer. Taken together, our results provide a more accurate picture of how gene expression evolves at each step of the expression process.

Strain specific proteome variations

Protein abundance is known to vary among individuals and is a major driver of phenotypic variation and adaptation to environmental fluctuations. Since the 8 strains used in our study were isolates from very different environments, we sought to detect possible adaptation in protein abundance. To this end, we performed a differential expression analysis for each strain.

Quantified protein abundance

Based on the replicated protein level data, we compared the abundance of each protein between isolates using the linear regression-based method LIMMA [START_REF] Ritchie | limma powers differential expression analyses for RNA-sequencing and microarray studies[END_REF][START_REF] Smyth | limma: Linear Models for Microarray Data[END_REF].

We used a one-vs-all strategy, which allowed us to accurately determine the strain-specific under-or over-abundance of each protein (Figure 2A, Figure S7). We considered a gene to be differentially expressed if the fold change and the FDR-adjusted p-value reported by LIMMA were respectively higher than 1.2 and lower than 10 -5 . Overall, across the 8 isolates, we detected 317 cases of differentially expressed proteins (DEP), corresponding to 244 unique proteins (Table S6, Figure S8). In addition, protein expression has been quantified in all samples for 68% of the DEP (i.e., no imputed expression), and for 92% of the DEP in half of the samples at least and at most completely missing from one or two strains. Hence, our imputation of missing values should have a low impact on the result of the differential expression analysis, since we also average the expression between all strains. Surprisingly, we found that the DEP number greatly varied across the isolate. It ranged between 9 for the AMH isolate and 193 for the BPL isolate. The high number of DEP in BPL is in line with its different proteome profile (Figure 1C). We looked for functional enrichment among the DEP using the gene ontology (GO) annotation and found that the enriched features were mainly related to either amino-acid metabolism or respiration function (Table S9, Figure S9). This suggests that intraspecies variability in terms of protein abundance seems to be the result of metabolic adaptations to different trophic conditions, as it is observed using the RNA-seq data (Caudal et al., 2023).

We then focused on isolate-specific adaptation of protein abundance. It was possible to relate some of the DEP to the environmental origin of the isolates. For example, BPL showed a very important overexpression of PDC5, a minor isoform of pyruvate decarboxylase that is essential for alcoholic fermentation. Since BPL is a wine strain, such an adjustment in protein abundance may reflect how selection for ethanol synthesis performance can affect protein abundance.

However, to get a broader view of protein abundance specificities for each strain, we performed GSEA on the Log2(FC) value obtained for each gene in each one-vs-all comparison (Table S8).

We found several signatures that recapitulated the environmental condition of the isolate. For example, BTT, a bioethanol strain, tended to overexpress the lipid metabolic pathway genes (Figure 2B), which is consistent with the central role of lipid adaptation and turn-over for ethanol tolerance in S. cerevisiae [START_REF] Eardley | Yeast Cellular Stress: Impacts on Bioethanol Production[END_REF][START_REF] Ma | Mechanisms of ethanol tolerance in Saccharomyces cerevisiae[END_REF][START_REF] Vanegas | Role of Unsaturated Lipid and Ergosterol in Ethanol Tolerance of Model Yeast Biomembranes[END_REF]. However, a large proportion of the enriched features in all the isolates were related to ATP metabolism and respiration (Table S8). The switch between respiration and fermentation being one of the signatures of S. cerevisiae domestication (Lahue et al., 2020), we checked if this signature was observable among the domesticated (namely, BPL, BTT, CQC and CPI). Interestingly, we found that domesticated isolates had different DEP pattern for the respiration related genes (Figure 2C). If BPL and BTT showed a clear underexpression of respiration-related genes, CPI showed no particular enrichment and CQC even showed an overexpression of many respiration-related genes. This could be consistent with the fact that, unlike wine or bioethanol production (for BPL and BTT), where yeasts are inoculated either by backsloping or addition of starter culture, cocoa fermentation in West Africa (the ecological origin of CQC) is a spontaneous fermentative process [START_REF] De Vuyst | Sourdough production: fermentation strategies, microbial ecology, and use of non-flour ingredients[END_REF][START_REF] De Vuyst | The cocoa bean fermentation process: from ecosystem analysis to starter culture development[END_REF][START_REF] De Vuyst | Functional role of yeasts, lactic acid bacteria and acetic acid bacteria in cocoa fermentation processes[END_REF][START_REF] Díaz-Muñoz | Phylogenomics of a Saccharomyces cerevisiae cocoa strain reveals adaptation to a West African fermented food population[END_REF][START_REF] Fernández Maura | The environmental and intrinsic yeast diversity of Cuban cocoa bean heap fermentations[END_REF][START_REF] Leroy | Lactic acid bacteria as functional starter cultures for the food fermentation industry[END_REF]. This likely explains why CQC has not shifted its metabolism contrarily to BTT or BPL and still primarily bases its energy production on respiration rather than fermentation. Taken together, our results show that protein abundance signature measured as DEP is a marker of environmental adaptation, but also of evolutionary history. (A) For each isolate (here CQC is used as an example), we detected DEP by applying a gene-wise linear model using the R package LIMMA. This allowed a one-against-all strain comparison to accurately detect isolate-specific over-or under-expressed genes. The p-value and Log2(FC) are

Transcriptional variation across the isolate is buffered across the gene expression

Interestingly, the proteome variations (mostly related to metabolism and respiration genes) seem at first highly similar to what was observed when looking at transcriptome and translatome variations (Teyssonniere et al., Submitted). We sought to compare more precisely the variations across the expression layers. We explored the inter-strain expression variations at each layer using two main approaches. First, we looked at the 28 pairwise correlations (corresponding to the correlations of each isolate against another) and observed a significant increase similarity between the expression profiles as long as the expression process goes on (median 𝜌transcriptome < median 𝜌Ribo-seq < median 𝜌proteome, figure 3A). This was also observed using the nonnormalized abundance (Figure S10A). Consistently, we quantified the variations level using an absolute log2 transformed fold change for each gene in each isolate pairwise comparison (|log2(FC)|, see methods). Briefly, the more this value increases, the more variable is the expression of a gene between two isolates. We found that the |log2(FC)| median value across the 28 pairwise comparisons were significantly decreasing as long as the expression process progresses (median |log2(FC)|transcriptome > median |log2(FC)|Ribo-seq > median |log2(FC)|proteome, figure 3B). Again, this was also observed using the non-normalized data (Figure S10B). The tendency of variation diminution at each step was also observed using Euclidean distances and gene-wise variance (Figure S10C,D). Taken together, these findings imply that gene expression is more constrained and therefore more conserved at the later steps of the process. This is in line with a phenomenon called post-transcriptional buffering that has been observed using both proteomics and ribo-seq data (Artieri and Fraser, 2014;Blevins et al., 2019;Dephoure et al., 2014;Gonçalves et al., 2017;McManus et al., 2014;Wang et al., 2020). Our data confirms that individuals. Furthermore, it suggests that the constraints shaping proteome variation might be different from those shaping transcriptome and proteome variation.

Gene expression evolution is gene specific

We sought to explore the determinants of gene expression evolution at each stage of the expression process. We used all categorical features from the Yeastomics dataset, which collects 3,685 gene characteristics from 27 studies (see full list and citations on https://github.com/benjamin-elusers/yeastomics). These characteristics include many gene and protein features related to e.g., chromosomal location, cellular function, interaction capabilities.

The overlap between this dataset and our expression dataset reached 2,308 genes. We selected the characteristics that affected at least 10 genes, resulting in a list of 793 features (Table S6).

To test whether each of the characteristics had an impact on the evolution of gene expression, we adapted a method previously published on mammalian expression data (Wang et al., 2020).

Briefly, this method relies on the construction of Euclidean distance trees based on gene expression across the 8 isolates (Figure 4A). The total length of the resulting tree serves as a measure of gene expression evolution. Using all the normalized expression of all genes, we observed that the resulting trees had different size: the proteomic based tree was the shortest while the transcriptomic based tree was the longest (Figure 4A). This is in line with our aforementioned observations on post-transcriptional buffering.

For each category, we then compared the length of the tree resulting from the included genes to a randomly generated length (see Methods), as the length of the tree is strongly correlated with the number of genes used for its computation (Figure 4B, Figure S12). For example, the category cat_genomics.sgd.chr_A (genes located on chromosome 1) encompassed 30 genes.

Using normalized abundance, the resulting tree has a total branch length of 21.72 (arbitrary unit), which is not significantly different from the lengths obtained with randomly generated trees (using 30 genes) (Figure 4C). The ratio between the median random tree lengths and the computed tree length are used as a measure of gene expression evolution. We tested the significance of the difference between the length and using a corrected (FDR) p-value threshold of 0.001, we detected 59 features influencing gene expression variation in at least one step of gene expression (mRNA abundance, translation, protein abundance) (Table S7), revealing that gene expression evolution is unequal across the genes. Several categories were based on previous expression variation exploration [START_REF] Lahtvee | Absolute Quantification of Protein and mRNA Abundances Demonstrate Variability in Gene-Specific Translation Efficiency in Yeast[END_REF] and were associated respectively with high and low expression variation in our dataset, supporting the reliability of the dataset and the tree-based exploration of expression evolution constraints.

We found that the overlap between the expression layer was overall small. Only 3 features influenced both mRNA and protein expression evolution. Similarly, 4 features influenced both translation and protein expression evolution while 22 features influenced both transcription and translation regulation evolution. At first glance, this suggests that the gene expression evolutionary constraints are layer specific, even if mRNA and RPF abundance seem to face similar constraints. However, we observed several trends that were conserved across each expression step. For instance, we found that the categories related to metabolism (Figure 4D)

were associated with faster gene expression evolution compared to the rest of the genes.

Consistently, in each dataset, the category associated with the fastest gene expression evolution across our 8 isolates was always related to metabolism (Table S7). This is in line with the DEP detected previously and other exploration in yeast highlighting that metabolism genes are usually among the most variable genes across individuals (Caudal et al., 2023). This important plasticity in gene expression is most likely a mechanism allowing for an optimal adaptation to different trophic specificities across our 8 isolates and reflects the environmental differences between the ecosystem in which each strain naturally occurs. Inversely, we found that several categories were associated with strong evolutionary constraints. For example, the feature associated with protein interactions resulted in the construction of short expression trees in each layer (Figure 4E). Accordingly, the most conserved categories were all related to interaction features. More generally, we found that central and essential cell functions were associated with evolutionary constraints. For instance, genes annotated as the core essential group were associated with constrained protein abundance evolution, while cytoplasmic translation related genes displayed a more conserved translation regulation. (A) Trees constructed from the gene expression-based Euclidean distances and their respective lengths. Consistent with post-transcriptional buffering, the proteomic tree is smaller than the others.

The error bars correspond to the bootstrapping (n=1000) performed on each tree. (B) The number of genes is a strong determinant of the length of the trees. The lines are constructed by smoothing the 1000 lengths obtained for each number of genes in each expression layer (see full plot in figure S7).

(C) For each gene feature (here cat_genomics.sgd.chr_A, i.e., the genes located on chromosome I) in each expression level (here RNA-seq), the length of the tree resulting from the selected gene is calculated (the dashed purple line). This length is then compared to the 1000 tree lengths generated using exactly the same number of genes (here n=30) if it is significantly shorter or longer.

Discussion

Gene expression is one of the major determinants of phenotypic variation observed between individuals. It is a complex phenomenon in which each step is tightly regulated. Here, we deeply quantified the proteome of 8 natural isolates of S. cerevisiae and combined these data with previously generated transcriptomes and translatomes on the same set of isolates (Teyssonniere et al., Submitted), allowing for a precise exploration of the expression variation throughout the gene expression process.

Our results highlight that protein abundance variation tended to be mainly related to metabolic pathways. Furthermore, we observed that protein abundance signatures of each isolate (detected as DEPs) often included genes related to respiration, which is a known marker of metabolic adaptation in the domestication of S. cerevisiae. These signatures adequately matched the wild or domesticated origin of the strain, and even highlighted differences within domesticated isolates related to their different uses. These observations are consistent with previous transcriptome and translatome surveys (Caudal et al., 2023;Teyssonniere et al., Submitted).

When comparing the variation between different expression levels, we observed that the expression variation between isolates tends to be buffered as the gene expression process progresses. This phenomenon, known as post-transcriptional buffering (Artieri and Fraser, 2014;Blevins et al., 2019;Dephoure et al., 2014;Gonçalves et al., 2017;McManus et al., 2014;Wang et al., 2020), highlights that the constraints on gene expression evolution are greater in the final stages of the expression process, which is consistent with the fact that proteins are the drivers of cellular machinery and functions. Abnormal changes in protein abundance can therefore be highly deleterious, and post-transcriptional buffering is often considered as a coping mechanism for faulty expression regulation (Buccitelli and Selbach, 2020;Liu et al., 2016). Our results show that this phenomenon is multilayered and its effect increases with the course of gene expression. Interestingly, we observed that although ribosome profiling is known to be a proxy for protein abundance (Brar and Weissman, 2015), the translatome was only slightly better at reflecting protein abundance within isolates than the transcriptome. When looking at the similarity between transcriptome, translatome and proteome variation across the expression layer in a gene-wise manner, we found that both translatome and transcriptome poorly reflected the proteome variation observed across isolates.

Finally, we focused on exploring the constraints that shape gene expression evolution. By adapting an analysis method previously developed (Wang et al., 2020), we could detect several gene characteristics that seemed to be involved in either fast or slow gene expression evolution.

gprofiler2 [START_REF] Kolberg | gprofiler2 --an R package for gene list functional enrichment analysis and namespace conversion toolset g:Profiler[END_REF][START_REF] Raudvere | Profiler: a web server for functional enrichment analysis and conversions of gene lists (2019 update)[END_REF] to functionally characterize the genes encompassed in our dataset. We used semantic similarity to reduce the number of detected GO annotations (biological process) with the rrvgo R package (Sayols, 2022).

Quantitative analysis of proteomes and differentially expressed proteins

The quantitative analysis of detected proteins was performed in R. First, we applied a logtransformation on expression intensities with a normalization of the data by the median of each sample to obtain relative protein expression within and across samples. To minimize information loss, we imputed missing protein expression values using the Bayesian Principal Component Analysis algorithm [START_REF] Oba | A Bayesian missing value estimation method for gene expression profile data[END_REF] as long as the protein was detected in 2 replicates of the same strain and among at least 4 strains. We computed the CV for each protein across the 8 strain and performed a GSEA using the R package FGSEA (REF) and with the GO annotation to explore which genes tended to display expression variability.

The LIMMA package [START_REF] Ritchie | limma powers differential expression analyses for RNA-sequencing and microarray studies[END_REF][START_REF] Smyth | limma: Linear Models for Microarray Data[END_REF] was used for identifying differentially regulated proteins, using the average across all four replicates, since the variance within strains was found to be much lower than between samples (cf Heatmap). We corrected the statistical significance for multiple testing of differential expression using the False-Discovery Rate procedure [START_REF] Benjamini | Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing[END_REF]. Then, a protein was considered differentially expressed if the fold change was higher than 1.2 and the adjusted p-value was lower than 10-5, as reported by LIMMA. For each isolate, we used the log2 transformed fold change of each gene to compute GSEA to detect strain-specific over or under expressed cellular pathway.

Expression variation exploration and comparison

We sought to explore gene expression variation all along the gene expression process. We combined our proteomic data with data previously generated on the same set of strains and under the same culture conditions (Teyssonniere et al., Submitted). This data, generated using RNA-seq and ribo-seq, comprises the precise measurement of 3,755 genes at the transcriptome and translatome level. The overlap between the two datasets covers 2,840 genes. We applied a LOESS normalization on the RNA-seq, ribo-seq and proteomic data together in order to strictly compare biological variations.

Using this dataset, we computed and compared the pairwise correlation levels of each isolate versus another using Spearman correlation test in each dataset. And we also computed the Euclidean distances between each isolate and computed the gene wise variance as well and checked if the resulting values were different depending on the expression level. Additionally, we computed for each gene, and in each isolate pairwise comparison, the absolute value of the log2 transformed fold change (|log2(FC)|) between two isolates as follow:

|𝑙𝑜𝑔2(𝐹𝐶)|𝑔𝑒𝑛𝑒 𝑋 𝑖𝑛 𝑐𝑜𝑚𝑝𝑎𝑟𝑖𝑠𝑜𝑛 𝑠𝑡𝑟𝑎𝑖𝑛 1 𝑣𝑠 𝑠𝑡𝑟𝑎𝑖𝑛 2 = |𝑙𝑜𝑔2( normalized expression gene X strain 1 normalized expression gene X strain 2 )|
For both pairwise correlations and absolute log2 fold change, the analyses were also performed on the non-normalized data.

Correlation between the expression levels

We sought to explore the proximity between Ribo-seq data and proteomic data. To do so, we computed two types of gene expression correlations: the across-gene correlation and the withingene correlation (Buccitelli and Selbach, 2020;Liu et al., 2016). The across-gene correlation is computed for each isolate (so in our case, 8 correlations per expression level) and it is based on the comparison between the gene expression data of all genes (in our case, the correlation between 2,840 values) in each expression level correlation (i.e., transcriptome vs translatome, transcriptome vs proteome and translatome vs proteome). We computed the correlation using a Spearman correlation test. The within-gene correlation is computed for each gene (in our case, 2,840 correlation per expression level) and is based on the gene expression values across the 8 isolates. Again, this was computed using a Spearman correlation for each expression level correlation: transcriptome vs translatome, transcriptome vs proteome and translatome vs proteome.

Gene expression evolution constraints

We looked for the determinants of gene expression evolution. We did so by using the Yeastomics database available online (https://github.com/benjamin-elusers/yeastomics for the full list, its construction, and citations). This database encompasses 3,685 gene characteristics that are either numeric (ex: specific codon composition, gene length, variance level in previous studies…) or boolean (ex: is on chromosome A, part of specific GO annotation, essentiality…).

We used the boolean characteristics for which at least 10 genes fulfilled the characteristics, which resulted in a list of 793 characteristics. The overlap between the gene characterized by this data and our 2,840 genes reached 2,308 genes. We used these characteristics to explore the constraints on gene expression evolution by adapting a tree-based approach previously published (Wang et al., 2020). This method relies on the construction of trees using the Euclidean distances generated from the 3 gene expressions. The total length of the resulting tree branches is used as a measure of gene expression evolution. Using all genes, we observed a difference in tree length depending on the expression level (statistically confirmed by a 1,000step bootstrapping) which was in line with the post-transcriptional buffering phenomenon and was a first control to support the reliability of this method. The tree length being highly dependent on the number of genes used for the tree construction, we computed for each number the lengths expected by chance for each characteristic by generating 1,000 trees constructed with randomly selected genes. For each gene characteristic, we selected the corresponding genes and constructed a tree for each gene expression level. We compared the resulting tree length with the 1,000 random lengths using a simple normal density probability test as the random lengths were normally distributed. The p-values obtained were FDR-corrected using and a threshold of 0,001 was considered to detect characteristics significantly associated with gene expression evolution. The ratio between the computed tree-length and the median random tree length was used to detect characteristics associated with expression evolutionary constraint (ratio <1) or with fast expression evolution (ratio>1). (A) The density area represents the distribution of protein expression determined from the average logarithmic intensity (log10) of detected peptides for each strain. A compact representation is also shown as a one-dimensional marginal distribution horizontally below each plot. The number of proteins quantified in each strain is written on the top left corner. A black density line indicates the distribution of protein expression after the imputation of missing values for proteins that were partially undetected among samples. The number of proteins partially detected is indicated in the top right corner. (B) The position of the circles relative to the x-axis is given by the average difference between imputed and quantified protein expression (Δ intensity) in each strain. For each strain, two arms extend in opposite directions to highlight the range of Δ intensity within one standard deviation from the mean represented by the circles.

Figure S3

: GO enrichment analysis of the encompassed proteins. Graphical representation of the GO enrichment analysis on the included genes from the proteomic data. The GO categories (in white) were obtained using semantic similarities on the terms detected using the gprofiler2 package in R. The semantic similarity was performed using the rrvgo package in R. (A) Heatmap of the Pearson correlation coefficients for pairwise comparisons of protein expression between proteomics samples. A hierarchical clustering based on the complete Euclidean distances between samples' expression profiles was applied to both rows and columns. The colors represent each isolate as shown in figure 1. (B) Principal component analysis of samples expression profile. The expression profiles of all samples are plotted as points scatterplot against the first (xaxis) and second (y-axis) principal components, which capture 52% of the variability between them. Points are colored according to their matching strains. The shaded regions delimit the range of variability among samples of the same strain. Volcano plot for each isolate. The points highlighted in purple are the DEP. The detection was performed using the LIMMA R package. For each gene number, 1,000 trees were constructed by randomly selecting the corresponding number of genes and the total length was computed.

CQC_1_1_10 CQC_1_2_12 BPL_2_2_12 BPL_2_1_10 BPL_1_2_12 BPL_1_1_10 CMP_1_1_10 CMP_1_2_12 CQC_2_2_10 CMP_2_1_10 CMP_2_2_10 BAN_1_2_10 BAN_2_1_10 BAN_1_1_12 BAN_2_2_10 AMH_2_1_12 AMH_2_2_12 BED_2_2_10 CQC_2_1_10 BED_2_1_12 BED_1_2_10 BED_1_1_12 BTT_1_1_12 BTT_1_2_10 BTT_2_2_12 BTT_2_1_12 CPI_1_2_10 CPI_2_1_12 CPI_2_2_12 CPI_1_1_12 AMH_1_1_12

Introduction

Understanding the genetic basis of phenotypic variation in natural populations is one of the main goals of modern biology. Gene expression differs among individuals and is known to be a main determinant of phenotypic variation (Albert and Kruglyak, 2015;Maurano et al., 2012).

In humans, the onset and development of numerous diseases have been linked to abnormal regulation of gene expression [START_REF] Cookson | Mapping complex disease traits with global gene expression[END_REF]. It is therefore essential to understand how genomic information is expressed through the different layers of gene regulation (i.e., transcriptomes and proteomes). Over the past decades, the development of methods for highthroughput quantification of mRNA and protein abundance has made it possible to explore both the proteome and the transcriptome on a larger scale (Messner et al., 2022b;[START_REF] Moyerbrailean | A highthroughput RNA-seq approach to profile transcriptional responses[END_REF]. These approaches facilitated the detection of numerous genetic loci (quantitative trait loci, QTL) affecting either transcript (eQTL) or protein (pQTL) levels (Chick et al., 2016;Ferkingstad et al., 2021;[START_REF] Folkersen | Genomic and drug target evaluation of 90 cardiovascular proteins in 30,931 individuals[END_REF]The GTEx Consortium, 2020[START_REF] Auesukaree | Molecular mechanisms of the yeast adaptive response and tolerance to stresses encountered during ethanol fermentation[END_REF], 2015).

However, the relationship between transcript and protein levels remains debated and poorly understood at the population level (Buccitelli and Selbach, 2020).

The transcript-protein correlation provides a first global view of the dependency of the two gene expression layers. Two types of mRNA-protein correlation can be determined, across-and within-gene, reflecting very different dynamics (Buccitelli and Selbach, 2020;Fortelny et al., 2017;Liu et al., 2016). The across-gene correlation analysis focuses on the overall correlation of a large set of genes coming from the same sample under a given condition to find out how well the absolute abundances of mRNAs and proteins are correlated. This correlation has been widely investigated in several species, such as human (Battle et al., 2015;Edfors et al., 2016;Gautier et al., 2016;Salovska et al., 2020;Wang et al., 2019;Wilhelm et al., 2014;Zhang et al., 2014), rats and mice (Aydin et al., 2023;Li et al., 2014;[START_REF] Moritz | Poor transcript-protein correlation in the brain: negatively correlating gene products reveal neuronal polarity as a potential cause[END_REF]Schwanhäusser et al., 2011), flies (Becker et al., 2018), plants (Ponnala et al., 2014) or yeast (Gygi et al., 1999;Ingolia et al., 2009;Marguerat et al., 2012). Across-gene correlations are consistently high and range from 0.4 to 0.8, suggesting that the absolute number of transcripts and proteins are globally correlated. Therefore, very abundant transcripts generally lead to very abundant proteins and vice versa.

However, the relationship between the transcript and protein abundance at the population level is explored via their variation across samples (e.g., individuals, tissues, or cell lines). Withingene correlation analysis gives a view on how the protein level of each gene tracks its mRNA level in a population. Different studies have investigated this within-gene correlation in different contexts and organisms, but they often show divergent results. Several surveys of tumors, normal human tissues, as well as pluripotent stem cells have highlighted this discrepancy in estimates with median within-gene correlation coefficients ranging from 0.14 to 0.59 (Archer et al., 2018;Aydin et al., 2023;Battle et al., 2015;Huang et al., 2017;[START_REF] He | Genome-wide pQTL analysis of protein expression regulatory networks in the human liver[END_REF]Mertins et al., 2016;Mirauta et al., 2020;Mun et al., 2019;Upadhya and Ryan, 2022;Vasaikar et al., 2019;Wang et al., 2019;Zhang et al., 2014Zhang et al., , 2016)). Similarly, the overlap of the detected loci influencing mRNA (eQTL) and protein (pQTL) abundance greatly differed across the datasets. It ranges from a very weak overlap of 5.5% in a study on 97 inbred and recombinant mice to nearly 35% in human (n = 62) and mice (n = 192) (Battle et al., 2015;Chick et al., 2016;Ghazalpour et al., 2011).

Part of the diverging results might have been driven by technical limitations. For instance, it has been shown that by selecting the most representative peptides in prior proteomic methods, the overall correlation of global transcript and mRNA abundance improves significantly (Alam et al., 2016;Upadhya and Ryan, 2022). A key difference is also whether the goal of the survey is to correlate absolute number of transcripts and proteins, or relative changes in protein or mRNA levels, which differ between samples. While the absolute number of transcripts and proteins spans several orders of magnitude, the relative expression differences of any individual protein across samples varies within a much narrower range (Marguerat et al., 2012;Messner et al., 2022a). Finally, a main limitation of these studies is that the sample size is much lower than the dimensionality of the problem.

To determine to which extent differences in relative changes in mRNA and protein levels are correlated and the genetic origins of their abundance variation are shared, a large-scale population survey exploring these two facets in a quantitative way was therefore necessary.

Here, we took advantage of the 1,011 yeast Saccharomyces cerevisiae population we genomesequenced and for which we have a species-level understanding of the natural genetic and phenotypic diversity (Peter et al., 2018). In order to be fully able to compare and analyze at unprecedented detail the relationship between these two layers of gene regulation, we therefore generated 942 quantitative proteomes in which cells were also cultured in synthetic complete medium supplemented with amino acids using high-throughput mass-spectrometry. We found that protein levels are molecular traits that exhibit considerable variation between individuals and specific signatures related to certain subpopulations. This large available population also makes it possible to generate a detailed map of loci involved in the variation of protein abundance (pQTL) at the species level, via genome-wide association studies (GWAS).

Interestingly, local pQTL are less frequent than distant ones (8% of the total set of pQTL) but they have a higher impact on their respective traits. Integration of proteomic and transcriptomic datasets acquired in parallel under similar conditions allowed comparison of accurate quantification of the mRNA and protein abundance of 629 genes across 889 natural isolates (Caudal et al., 2023). Based on these unique datasets, we clearly demonstrated that the degree of within-gene correlation between protein and mRNA abundance is very low (Rho = 0.165).

Consistently, we found that the genetic variants influencing protein and mRNA abundance are very dissimilar. Our study highlights that population-scale proteomes are essential and add a new dimension to the characterization of the genotype-phenotype relationship when integrated with genomic and transcriptomic information.

out of 629 consistently quantified proteins fall into the 20% highest transcribed genes (n = 1,304). In total, 537 out of 629 quantified proteins were found in the two highest abundance deciles as defined in a recent yeast protein abundance meta-analysis (Ho et al., 2018) (Figure S1E). Overall, proteins related to essential genes and involved in molecular complexes were both significantly enriched in the set of proteins quantified by Scanning SWATH (odd-ratio = 3.5 and 2.2 respectively, Fisher's exact test, p-values < 2.2x10 -16 ) (Dowell et al., 2010;Giaever et al., 2002;Pu et al., 2009). Function-wise, we found that metabolism-related genes were overrepresented among the 629 genes included in our study (Table S4).

We then investigated the level of variation in protein abundance by calculating the coefficient of variation (CV) for each protein using the non-normalized dataset. We found an average CV of 31%, varying between 12% and 98% and one high outlier reaching 300% (PDC5, a pyruvate decarboxylase). The precursor-level CVs across quality control samples (15.15%) were much lower than the precursor-level CVs across the natural isolate samples (34.21%), confirming that a biological signal was observed across the isolates (see Methods). Gene set enrichment analyses (GSEA) were performed using the CVs and significant enrichment of genes related to amino acid metabolism, respiration or pyruvate metabolism was found for proteins with a high CV, indicating that they vary the most (Table S5). By contrast, proteins with a low CV were significantly related to genes involved in tRNA aminoacylation or protein degradation. A. The proteomic dataset was generated on isolates grown in synthetic complete (SC) medium with amino acids using a semi-automated sample preparation workflow, and Scanning-SWATH MS (see Methods). The overlap between this dataset and the recently generated transcriptomic dataset on the same population in the same condition (Caudal et al., 2023) resulted in 629 protein/transcript abundances across 889 isolates. B. Phylogenetic trees of the isolates used in this study. Colors correspond to previously defined subpopulations (Peter et al., 2018). C. Gene-wise correlation coefficients (Spearman correlation test) between the proteome and the transcriptome. D. and E. mRNA-protein within-gene correlation across isolates for the RPL38 and GLR1 genes (ρ corresponds to the Spearman correlation coefficient with p-values of 4.8x10 -7 and 2.3x10 -153 , respectively).

Transcript and protein abundances are weakly correlated at the gene level across isolates

As proteomes and transcriptomes were obtained using the same growth media, our dataset allowed us to characterize the different types of correlation between mRNA and protein abundance across a natural population. We first determined the across-gene correlation, i.e. the concordance between protein and transcript abundance for each isolate, and found a very high correlation (median rho = 0.53, interquartile range of 0.06, Figure S2), which is consistent with what was previously described (Battle et al., 2015;Becker et al., 2018;Edfors et al., 2016;Gautier et al., 2016;[START_REF] Moritz | Poor transcript-protein correlation in the brain: negatively correlating gene products reveal neuronal polarity as a potential cause[END_REF]Ponnala et al., 2014;Salovska et al., 2020;Wang et al., 2019;Wilhelm et al., 2014;Zhang et al., 2014). We next computed the correlation between the protein and mRNA normalized abundance for each gene across the 889 natural isolates (Figure 1C-D-E, Table S6). While the across-gene correlation levels were in line with previous explorations, we found an overall very low within-gene correlation level (median ρ = 0.165, interquartile range of 0.17). This value is much lower than the one determined with smaller samples in mice (approximately 0.25) (Chick et al., 2016;Ghazalpour et al., 2011) and in human healthy tissues (0.35 and 0.46) [START_REF] He | Genome-wide pQTL analysis of protein expression regulatory networks in the human liver[END_REF]Wang et al., 2019), but it is in line with what was found in human lymphoblastoid cell lines (0.14) (Battle et al., 2015). For a total of 385 out of the 629 quantified proteins, the level is significantly correlated with RNA level (Bonferroni corrected p-value < 0.05). Out of these 385 proteins, only 3 show a negative correlation: Rps13, Asc1 and Rpl38 (Figure 1D), all ribosomal related proteins. This observation is consistent with previous surveys pointing out that some ribosome-related proteins are negatively correlated with their cognate transcripts (Buccitelli and Selbach, 2020;Wang et al., 2019). But overall, this correlated set of 385 proteins/transcripts is significantly enriched of genes related to several metabolism pathways (Table S7). Moreover, the most strongly correlated set of proteins/transcripts (n = 33) show functional enrichment of genes related to mitochondrial respiration (Table S8) (see Methods). Interestingly, it points out that this specific pathway has similar gene regulation at both levels. Finally, we observed that four genes with very high mRNA-protein correlation were located outside the main correlation index distribution (Figure 1C). These genes all have correlation coefficients greater than 0.6: SFA1 (alcohol dehydrogenase), HBN1 (unknown function), GLR1 (glutathione oxidoreductase, Figure 1E) and YLR179C (unknown function). Such a high correlation clearly points to common regulatory mechanisms and genetic bases underlying the two levels of variation, as we have seen below.

Gene expression is more constrained at the proteome level

By combining these proteomic and transcriptomic datasets, we are in a position to simultaneously explore and compare the variation of these two gene expression layers at the population level. We therefore computed the absolute Log2(fold change) value (i.e., |Log2(FC)|) for each gene in each pair of isolates and found that this value is 32% lower on average for the proteome (Figure 2A), suggesting that protein abundance is less variable and more constrained than mRNA abundance. Furthermore, a higher correlation was observed between proteomes (rho = 0.92) compared to transcriptomes (rho = 0.83) (Figure 2B). Finally, the variance observed for each gene was lower for the proteomic data (Figure S3A) and the Euclidean distances between each isolate were smaller when computed with the protein abundance dataset (Figure S3B). Overall, these observations reflect and highlight the presence of a global post-transcriptional buffering of the transcriptome variations.

Despite recurrent observations (Blevins et al., 2019;Kustatscher et al., 2017;McManus et al., 2014;Muenzner et al., 2022;Wang et al., 2020), the post-transcriptional buffering phenomenon remains largely functionally uncharacterized and poorly understood. We sought to better understand this phenomenon at the genetic level by examining the cellular functions that tended to be most affected by post-transcriptional buffering. Briefly, we constructed neighbor-joining trees using the proteome or transcriptome Euclidean distances between each isolate (Figure 2C and see Methods) (Wang et al., 2020). Total branch length was used as a measure of expression variation and evolution at the species level. We then calculated the ratio between the lengths of the proteome and transcriptome tree branches to quantify the strength of the post-transcriptional buffering phenomenon. The lengths of branches from the proteome-based tree were shorter than those from the transcriptome-based tree, resulting in a length ratio of 0.93 (Figure 2C, Figure S3C). This observation is consistent with the differences in Euclidean distances observed previously (Figure 2B). We then applied the same procedure to 101 sets of genes, representing central biological processes obtained from a reduced list of gene ontology (GO) annotations (Table S9). We found that a total of 16 sets display a ratio lower than 0.93 and a significant difference between the proteome and transcriptome branch lengths, meaning that these sets are strongly affected by the phenomenon of post-transcriptional buffering (Figure 2D, Table S10).

Interestingly, 6 out of the 16 sets include genes with functions related to protein production and maturation (Figure 2D), highlighting that the evolution of the cellular machinery involved in protein production and maturation is highly constrained. The other set of genes are related to several metabolism processes and detected as strongly buffered, despite being highly variable in the proteomic data (Table S5). This observation could be due to the fact that metabolismrelated genes are among the genes with the greatest variation in mRNA abundance at the species level (Caudal et al., 2023). This variation is largely attenuated at the proteome level but remains important, reflecting differences in metabolic preferences within the population. Moreover, we also found 3 sets with a ratio higher than 0.93 and a significant difference between the proteome and transcriptome trees, which means that the expression variation of these genes is greater at the proteome level (Figure 2D, Table S10). Interestingly, all of them are related to protein catabolism, highlighting a difference in post-transcriptional mechanism for this specific functional category. Taken together, these results provide new insights into post-transcriptional buffering as well as its functional impact.

displaying a ratio lower than 0.93 and a significant difference in terms of branch lengths (see Methods)

were considered as strongly buffered. The biological processes displaying a ratio higher than 1 and a significant difference in term of branch lengths had an enhanced abundance variation at the proteome level. D. Biological processes detected as strongly buffered or with an enhanced variation using the procedure detailed in C.

Architecture of the proteome landscape

Using these datasets, we then sought to understand the main determinants shaping the proteome architecture at the population level. The S. cerevisiae yeast species exhibit a clear population structure, which potentially can impact the proteome landscape (Peter et al., 2018) (Table S1).

We performed a principal component analysis (PCA) with the protein abundance data and found that no clear grouping emerged from the subpopulations when plotting together the 6 first principal components (Figure S4A-B-C). The same results were observed for transcriptomes (Figure S4D-E-F). To confirm this, we also computed the Euclidean distance across transcript and protein levels between every pair of isolates and used these to construct a neighbor-joining tree (Figure 3B-C). We observed that none of the subpopulations present in the genetic-based tree merged in either the proteome-or transcriptome-based tree (Figure 3A-B-C). Together, these results highlight that population structure does not impact transcriptomes and proteomes in the S. cerevisiae species.

One potential determinant of the proteome organization could be related to co-expression networks that strongly influence the coordination of gene expression or various cellular processes. Using Weighted Gene Co-Expression Network Analysis (WGCNA) [START_REF] Zhang | A general framework for weighted gene co-expression network analysis[END_REF] on the normalized protein abundance data, we detected seven co-expression modules (Figure 3D, Table S11). Each of these modules corresponds to a specific biological function (Table S12, Figure S5) and encompasses between 38 (Cellular amino acid biosynthetic process) and 114 (Ribosome biogenesis) genes. Interestingly, very similar modules were found applying the same procedure on the mRNA normalized data. Five co-expression modules were detected (Figure 3E, Figure S6, Table S13, Table S14), and all of them were detected in the seven proteomic modules, suggesting that co-expression patterns recapitulate central cell functions are conserved across the two expression layers (Figure 3D-E).

Insight into subpopulation-specific protein expression

We further wanted to explore and determine the presence of subpopulation-specific signatures.

We therefore sought to identify differential protein expression patterns by comparing each clade to the rest of the population and we detected a total number of 1,129 differentially expressed proteins (DEPs) (corresponding to 465 unique proteins, Figure S7, Table S15). An average of 59 DEPs was found per clade, ranging from 218 for the Wine clade to 0 for wild Asian clades represented by a small sample (e.g. CHN, Taiwanese and Far East Russian) (Figure S8A).

Several DEPs were adequately related to the ecological origin of the different subpopulations.

Proteasome complex

For example, several subpopulations related to alcoholic fermentation show overexpression of alcohol dehydrogenases, such as ADH4 in Wine and Brazilian bioethanol clades as well as ADH3 in the Sake subpopulation. In the French Dairy subpopulation, we also observed an underexpression of SEC23, a GTPase-activating protein involved in the COPII related vesicle formation, which could reflect an adaptation to this secretory pathway to the cheese-making environment [START_REF] Celińska | Filamentous fungi-like secretory pathway strayed in a yeast system: peculiarities of Yarrowia lipolytica secretory pathway underlying its extraordinary performance[END_REF]. Overall, these observations suggest that domestication and more generally, ecological constraints are drivers of the proteomic landscape evolution in a natural population. We then performed GSEA based on differential expressed proteins in each subpopulation and found significant enrichments for various biological processes (Figure 4A, Table S16). Many enriched functional categories were associated with respiration related genes (e.g. "respiratory electron chain transport"). Interestingly, we observed that while most wild clades (8 out of 13) tend to have overexpression of respirationrelated proteins, these are underexpressed in domesticated subpopulations (5 out of 7).

We therefore further explored the impact of domestication on the proteome at the population level. Using the same DEP detection method, we assessed the proteome differences between the domesticated and wild isolates (Peter et al., 2018) and found a total of 133 DEPs (Table S17). Among these proteins, other alcohol dehydrogenases such as SFA1 and ADH3 were highly abundant in domesticated isolates. A GSEA performed on this set of DEPs clearly shows an enrichment of underexpressed respiration-related proteins in domesticated clades (Figure 4B, Table S18). Unlike wild isolates, domesticated isolates were selected for fermentation purposes, likely leading to this specific signature. This observation is in line with the previous finding pointing out that the switch from a preference between respiration and fermentation is one of the hallmarks of domestication in yeast (Lahue et al., 2020). In addition, significant enrichment of the functional category "chaperon mediated protein folding" points to overexpression of this set of proteins in the domesticated isolates (Figure 4B), which may be an adaptative response to long-term exposure to ethanol, known to induce protein denaturation [START_REF] Auesukaree | Molecular mechanisms of the yeast adaptive response and tolerance to stresses encountered during ethanol fermentation[END_REF]. By performing the same analysis on transcriptomic data (Figure S8B, Table S17), similar results, showing overexpression of respiration-related genes in domesticated clades, were obtained (Table S19).

The genetic bases of protein abundance at the population scale

To uncover the genetic origins of the proteome variation at the population-scale, we performed genome-wide association studies (GWAS) and considered both SNPs and CNVs that were characterized previously (Peter et al., 2018). We focused on isolates for which both proteomic and transcriptomic data were available, resulting in a set of 889 isolates. In this population, a total of 84,633 SNPs and 1,019 CNVs were considered, with a minor allele frequency higher than 5%. We performed GWAS using the raw protein abundances of the genes for which we have both levels of expression (i.e., 629 genes). Overall, we detected a total of 598 SNP-pQTL after colliding SNP affected by linkage disequilibrium (R2 > 0.6), and 4,528 CNV-pQTL corresponding to 501 and 520 loci and affecting 300 and 93 genes, respectively (Figure 5A-B, Table S20, Table S21, data file 1).

Among the SNP-pQTL, 8% (n = 50) were local-pQTL, showing that regulation of protein abundance is primarily achieved through trans regulation. This fraction is consistent with previous exploration in yeast (Foss et al., 2007) and lower than what is usually found at the transcriptome level (Albert et al., 2018;Caudal et al., 2023). Nonetheless, we observed that the local SNP-pQTL have a higher effect size compared to trans SNP-pQTL (Figure 5C) and tend to be located near the transcription starting site of the gene (Figure S9). We found no strong SNP-pQTL hotspots, suggesting that most of the distant pQTL are evenly distributed throughout the genome (Figure 5D).

In contrast, CNVs impacting protein abundance had a biased location on the chromosomes 1, 3, 8, 9 and 11 (Figure 5B). Out of 4,528 CNV-pQTL, a total of 4,303 were located on these chromosomes and affected a gene on their respective chromosome. This observed bias is due to the presence of aneuploidies on these chromosomes in our population (Peter et al., 2018).

These CNV-pQTL have also a higher impact on the protein abundance variation compared to the other CNV-pQTL, suggesting that aneuploidies represent a major source of proteome variation at the population level (Figure S10). Only 24 local CNV-pQTL out of 4,528 were detected, and no significant effect size between local and distant CNV-pQTL was found (Figure 5C).

We then looked at the extent to which the genetic bases of protein abundance are common with those underlying the abundance of transcripts. We performed GWAS using the transcriptomic dataset and detected 596 SNP-eQTL and 4,877 CNV-eQTL (Figure S11, data file 2), which is of the same order of magnitude as the GWAS proteome results. Surprisingly, the overlap between the SNP-pQTL and the SNP-eQTL is very low, with only 3.6% of shared SNP-QTL (n = 22). Interestingly, 18 out of 22 were related to local regulation, meaning that 36% of the local SNP-pQTL (18 out of 50) also impact the cognate transcripts of their target protein. This observation is consistent with previous findings showing that the common regulation between mRNA and protein abundances is mainly related to local regulation (Chick et al., 2016;Ghazalpour et al., 2011). Overall, we observed that genes with a strong correlation between transcript and protein abundance, such as the top four most correlated genes previously mentioned (SFA1, HBN1, GLR1 and YLR179C), tend to have a shared pQTL and eQTL (Figure S12). Additionally, we found that the SNP-pQTL distribution across the genome did not match the SNP-eQTL distribution, where a QTL hotspot could be detected around the CTT1 gene (Caudal et al., 2023;[START_REF] Stuecker | Linkage mapping of yeast cross protection connects gene expression variation to a higher-order organismal trait[END_REF]. The reasons for the weak overlap are likely multifactorial, but protein-specific regulation, such as protein degradation, may play a central role. We sought to confirm this by looking at the average protein turnover (Muenzner et al., 2022) of the proteins with and without overlapping pQTL and eQTL (Figure S13A, see Methods). We found that proteins, for which an overlap between pQTL and eQTL was detected, show a lower turnover rate compared to the other proteins. Consistently, the half-life of proteins with an overlapping SNP-QTL was higher than the rest of the proteome (Figure S13B). This observation suggests that protein degradation is probably involved in the large differences observed between the genetic origins of mRNA and protein abundance.

In contrast, the overlap between the two sets of CNV-QTL is much higher, as 3,097 QTLs were shared between the transcriptome and proteome, i.e., approximately 68% of the CNV-pQTL.

However, these shared CNV-QTLs are all aneuploidy-related CNVs, suggesting that the effect of aneuploidies is persistent through the expression layers (Muenzner et al., 2022). None of the non-aneuploidy CNV-QTL (22 CNV-eQTL and 216 CNV-pQTL) were shared. Together, our results highlight that the genetic bases underlying population-level protein abundance are very distinct from those underlying mRNA abundance.

Discussion

Quantifying transcripts and proteins expressed in a large natural population is fundamental for having a better understanding of the genotype-phenotype relationship. In this study, we have quantitatively analyzed the proteome of 942 natural isolates of S. cerevisiae, allowing in-depth exploration of protein abundance and precise characterization of the genetic origins of its variation at the species level.

The S. cerevisiae species is characterized by a complex population structure, with domesticated and wild subpopulations (Peter et al., 2018). Structured populations are also observed in a large number of other species, such as humans, and their impact on the proteome remains unexplored.

In our dataset, the population structure had no significant impact on the proteomic landscape.

This observation is consistent with previous results obtained with the transcriptomes of S.

cerevisiae isolates (Caudal et al., 2023;Kita et al., 2017). In fact, most subpopulations are characterized by specific signatures related to a small set of genes but not to a general pattern.

This dataset allowed us to have better insight into the architecture of the species-wide proteome variation. First, we found that the co-expression network captures main biological functions and is globally conserved across the species. Second, we detected differential protein expression signatures specific to subpopulations, reflecting an adaptation to specific ecological conditions, such as domesticated environments. Similar expression signatures can be also observed using transcriptomic data (Caudal et al., 2023;Hodgins-Davis et al., 2012), highlighting that gene expression plasticity at both levels is a key mechanism of environmental adaptation.

The species-wide proteomes and transcriptomes obtained in the same condition represent a unique opportunity to compare the gene regulation at both levels. The overall agreement between protein and transcript within each isolate appears to be high and this in the whole population, showing again that very abundant transcripts generally lead to very abundant proteins and vice versa (Battle et al., 2015;Becker et al., 2018;Edfors et al., 2016;Gautier et al., 2016;[START_REF] Moritz | Poor transcript-protein correlation in the brain: negatively correlating gene products reveal neuronal polarity as a potential cause[END_REF]Ponnala et al., 2014;Salovska et al., 2020;Wang et al., 2019;Wilhelm et al., 2014;Zhang et al., 2014). However, our data allow for the first time to have an accurate estimation of the correlation per gene at the population level and we found that this gene-wise correlation is very weak with a median of 0.165, which is lower than most previous estimates based on much smaller human and mice populations (Chick et al., 2016;Ghazalpour et al., 2011;[START_REF] He | Genome-wide pQTL analysis of protein expression regulatory networks in the human liver[END_REF]Upadhya and Ryan, 2022;Wang et al., 2019). Consistent with this result, genome-wide association studies also highlighted that SNPs related to variation in protein (pQTL) and transcript (eQTL) levels poorly overlap (3.6%), with mostly common local QTL. This result is consistent with one of the first eQTL/pQTL comparisons (Foss et al., 2007) but unlike other studies, showing a higher overlap (Albert et al., 2014). However, we should emphasize that we were not able to map the genetic basis of the entire S. cerevisiae proteome and therefore the eQTL/pQTL overlap might be biased and underestimated.

Mechanistically, our results suggest that the regulation of protein degradation has an impact on the variation of the proteome, and therefore on its genetic basis. Proteins with a high turnover rate will be more affected by proteome-specific regulation and will therefore show a weaker correspondence with the transcriptome. Conversely, proteins with a low turnover rate are more likely to be impacted by variation in transcript abundance. They will therefore likely reflect variation in mRNA abundance.

Although mass spectrometers are highly sensitive, it should be noted the limitation that proteomic methods are biased towards quantification of highly abundant proteins. Indeed, the fraction of the proteome quantified constitutes the vast majority of the total proteomic mass of a cell and is enriched for essential genes as well as in genes most connected in functional networks. Our dataset captures many of the fundamental processes. Yet, results related to low abundant proteins are missed by this approach.

Overall, our study clearly highlights that the dependency between transcript and protein levels is complex, pointing to the importance of post-transcriptional regulation of protein abundance.

Proteome and transcriptome are indeed two distinct layers of gene regulation, which need to be further explored to understand the genotype-phenotype relationship. As gene function is ultimately executed by the proteome, while mRNA is the messenger, more proteomic approaches will be needed to create a better understanding of the phenotypic diversity. Our study provides a first species-wide insight into the genetics that underlies both proteome and transcriptome diversity in a natural population.

LC-MS/MS measurements

In brief, peptides were separated on a 3-min high-flow chromatographic gradient and recorded by mass spectrometry using Scanning SWATH [START_REF] Messner | Ultra-fast proteomics with Scanning SWATH[END_REF] using an online coupled 1290 Infinity II LC system (Agilent) -6600+ TripleTOF platform (Sciex). 5 µg of sample were injected onto a reverse phase HPLC column (Luna ® Omega 1.6µm C18 100A, 30 × 2.1 mm, Phenomenex) and resolved by gradient elution at a flow rate of 800 µL/min and column temperature of 30 °C. All solvents were of LC-MS grade. The gradient program used 0.1% formic acid in water (Solvent A) and 0.1% formic acid in acetonitrile (Solvent B) and was as follows: 1% to 40% B in 3 min, increase to 80% B at 1.2 mL over 0.5 min, which was maintained for 0.2 min and followed by equilibration with starting conditions for 1 min. For mass spectrometry analysis, the scanning swath precursor isolation window was 10 m/z; the bin size was set to 1/5th of the window size, the cycle time was 0.7 s, the precursor range 400 m/z to 900 m/z, the fragment range 100 m/z to 1500 m/z as previously described [START_REF] Messner | Ultra-fast proteomics with Scanning SWATH[END_REF]). An IonDrive TurboV source was used with ion source gas 1 (nebulizer gas), ion source gas 2 (heater gas), and curtain gas set to 50 psi, 40 psi and 25 psi respectively. The source temperature was set to 450 °C and the ion spray voltage to 5500 V.

Data processing

The mass spectrometry files were processed following the approach previously described (Muenzner et al., 2022). Briefly, an experimental spectral library obtained using the S288c was filtered to reduce the search space to peptides well shared across the strains. This library was then used with the software DIA-NN (Demichev et al., 2020) (Version 1.8) and the following parameters: missed cleavages: 0, Mass accuracy: 20, Mass accuracy MS1: 12, scan windows:

6. The option 'MBR' was used to process the data. As the peptides selected were not necessarily present ubiquitously in all the strains, an additional step was required to remove false positives (entries where a peptide is detected in a strain where it should be absent). This represents only ~1% of the total entries of the report. S22). Batch correction was carried out at the precursor level using median batch correction, which consists in bringing the median value of the precursors in the different batches to the same level. Proteins were then quantified from the peptide abundance using the maxLFQ (Cox et al., 2014) function implemented in the DIA-NN R package. The resulting dataset consists of 630 proteins for 942 strains. We imputed the missing value for further exploration using the KNN imputation method from the impute R package [START_REF] Hastie | impute: impute: Imputation for microarray data[END_REF].

Combination of transcriptomic and proteomic data

Unless specified, all the analysis performed below were conducted using R version 4.1.2. The transcriptomic data was generated previously (Caudal et al., 2023). We used the log2 transcript per million (TPM) data, where the overlap with proteomic data was encompassing 629 genes across 889 isolates, for the genome wide association studies (see later for the method). For the exploration of gene expression variation, subpopulation related DEG and gene expression network, we used the variance stabilized data obtained directly from the log2 TPM data. In this case one gene was removed from the analysis and the reference strain data was not considered, which resulted in an overlap of 628 genes across 888 isolates. To only focus on real expression variation difference between the expression layers, we normalized the proteomic and transcript abundance using quantile normalization. Unless specified, all the analyses described below use the quantile normalized transcriptomic and proteomic data. We recomputed the raw protein abundance coefficient of variation (CV) of each gene by dividing the standard deviation by the mean (using the non-normalized abundance) and transformed it to a percentage. Based on the CV, we performed a functional exploration by gene set enrichment analysis (GSEA) (Subramanian et al., 2005) using the fgsea R package (Korotkevich et al., 2021) for the gene ontology annotation (Ashburner et al., 2000;Gene Ontology Consortium, 2021) to detect cellular pathways with a conserved regulation across the population. The within-and acrossgene mRNA-protein correlation was performed for each gene or each isolate using a Spearman correlation test. We selected the genes with a mRNA-protein correlation index higher than 0.42 (> 95% percentile) and performed gene ontology (GO) enrichment analysis using the biological process (BP) database using the topGO R package [START_REF] Alexa | topGO: Enrichment Analysis for Gene Ontology[END_REF]. For the was significantly different after 10 bootstrapping steps (Bonferroni corrected Wilcoxon test pvalue < 0.001).

Transcriptome and proteome landscape exploration

We sought to check if the genetic structure of the population had an impact on the transcriptome and proteome structure. We obtained the genetic distances from (Peter et al., 2018) between pairs of isolates and compared them to the pairwise isolate correlation (Spearman correlation test) obtained with the normalized transcript or protein abundances. We also used both normalized protein and mRNA abundance data to perform principal component analysis (PCA) using the prcomp function from the stats R package. For the 2 PCA (transcriptomic and proteomic), we plotted the 6 first principal components (PC) together (PC1-PC2, PC3-PC4 and PC5-PC6) and looked for eventual grouping according to the subpopulation as defined previously (Peter et al., 2018). We then computed a Weighted Gene Co-Expression Network Analysis (WGCNA) using the WGCNA R package [START_REF] Langfelder | WGCNA: an R package for weighted correlation network analysis[END_REF] to detect co-expression module in both mRNA and peptide normalized abundance. To do so, we generated a Topological Overlap Matrix (TOM) using the blockwiseModules function. The TOM were calculated based on a signed adjacency matrix with the power of 9 for the mRNA abundance data and 5 for the peptide abundance data. The blockwiseModules automatically detected the co expression modules by generating a clustering from a dissimilarity matrix (1-TOM) using the following option: detectCutHeight = 0.995; minModuleSize = 30. This resulted in the detection of 5 and 7 transcriptome and proteome modules respectively. We computed an overrepresentation analysis for each co-expression module with the GO terms as annotation and using the mod_ora function from the CEMiTool R package [START_REF] Russo | CEMiTool: a Bioconductor package for performing comprehensive modular coexpression analyses[END_REF] and used the most representative GO terms as the final annotation for each detected module. The two coexpression networks were generated for plotting by computing an adjacency matrix from the TOM matrix (generated previously) and ultimately plotted using the ggnet2 function from the GGally R package.

Transcriptome and proteome differentially expressed gene detection

We used the normalized protein abundance to detect subpopulation-specific (Peter et al., 2018) differentially expressed proteins (DEPs). The goal was to detect either over-or underexpressed genes by comparing the normalized expression of all the isolates from a subpopulation against the rest of the population using a Wilcoxon test for each gene. The p-value of the test was corrected using a Bonferroni correction with the p.adjust function in R. A gene was considered as differentially expressed if the corrected p-value of the Wilcoxon test was below 0.05. We computed as well a log2 transformed fold change (log2(FC)) value for each gene in each subpopulation using the mean expression of the subpopulation divided by the mean expression of the rest of the population. To further characterize the detected DEPs, we performed a functional exploration using GSEA (with the fgsea function from the fgsea R package) using the log2(FC) value from the DEP exploration as score rankings. In order to have a global view of the pathways that were significantly differentially expressed in each subpopulation, we used the 16 co-expression modules detected and defined previously using the population transcriptome data (Caudal et al., 2023) as biological function annotations for the GSEA. We performed the same procedure but this time comparing the domesticated against the wild isolate using the clade wise annotation from (Peter et al., 2018). This time, the test was performed on both normalized protein and transcript abundances.

Proteome and transcriptome genome-wide association studies

We computed GWAS with a linear mixed model-based method as described previously (Caudal et al., 2023;Peter et al., 2018) using FaST-LMM (Lippert et al., 2011). We performed the GWAS using either the transcriptome log2 transformed TPM data or the protein abundance.

For each dataset, we performed two separated GWAS, one based on SNP as genotype, and one based on the CNV as genotype. The SNP GWAS was run with total of 84,633 SNP displaying a minor allele frequency (MAF) > 5% and that were not located in the telomeric regions (< 20kb away from the chromosome ends). The CNV GWAS was run on a total of 1,019 CNV (MAF > 5%). We used the SNP matrix for both SNP and CNV GWAS, thus evaluating the kinship between the isolate to account for the population structure. We set a phenotype-specific p-value threshold using 100 permutation tests where the phenotypes were randomly permutated between the isolates. We use the 5% lowest p-value quantile from these permutation tests to define the significance threshold. We finally scaled the significance thresholds of the CNV GWAS to account for the size difference between the SNP and CNV matrices.

Regarding the SNP GWAS, the detected QTL were filtered to avoid false positives detection due to linkage disequilibrium among the SNP as described previously (Caudal et al., 2023).

This resulted in the filtration of 81 eQTL and 131 pQTL (out of respectively 677 and 729 QTL).

The QTL were considered as "local" QTL when they were located 25 kb around their affected (A) number of isolates encompassed in the proteomic datasets, in the transcriptomic dataset (Caudal et al., submitted) and in the overall population (Peter et al., 2018). The x-axis corresponds to the clades (or subpopulations) as defined previously (Peter et al., 2018). (B, C) Expression values of 5 randomly selected isolates for protein and transcript abundance before (B) and after (C) quantile normalization.

(D) mRNA levels of the gene encompassed or not in the proteomic data (*** = p-value < 2.2x10 -16 , Wilcoxon test). (E) Protein levels (as defined in Ho et al., 2018) of the genes encompassed by our proteomic data. Functional enrichment of each co-expression module detected using WGCNA on protein abundance data. The enrichment was performed using the CEMiTool package. The dotted lines on each graph represent the significance threshold. Volcano plots for each subpopulation highlighting the DEPs. The blue points correspond to underexpressed gene in a subpopulation while the red points correspond to over-expressed genes.

CV density group

QCs Samples aspects such as the chromatograph time and liquid flow has been drastically reduced. This was an adaptation of a previously published study that already allowed a theoretical throughput of 180 samples per day (Messner et al., 2020). Covering a set of 630 proteins in more than 900 isolates at this rate is a novelty in itself. However, several improvements are still needed, especially in peptide signal acquisition or proteomic data handling, as this is still a limited number of genes covering only 10% of the theoretical proteome, that is moreover biased toward highly expressed genes. Increasing the proteome coverage would be crucial, in particular to capture some more population-specific trends. For example, the LAC and MAL genes that were related to strong subpopulation-specific transcription signals in the French dairy and beers populations, respectively (Caudal et al., 2023) were not quantified in our population-scale proteome exploration.

Exploration of other culture conditions

Throughout this work, yeast cultures have been performed on complete synthetic medium (at 30°C) to provide a nutrient-rich and controlled environment. However, such culture conditions obviously do not reflect the natural environment for the vast majority of the 1,011 isolates we worked with, especially since S. cerevisiae is ubiquitously distributed on earth across both human and wild niches, and thus face a vast diversity of trophic conditions [START_REF] Bai | The Ecology and Evolution of the Baker's Yeast Saccharomyces cerevisiae[END_REF]Peter et al., 2018;[START_REF] Wang | Surprisingly diverged populations of Saccharomyces cerevisiae in natural environments remote from human activity[END_REF]. Therefore, an interesting continuation of this work could be the study of gene expression in various growth culture conditions, closer to the environmental constraints faced by certain subpopulations, for which subpopulation-specific gene expression could be much stronger. Indeed, the strains adapted to specific conditions would most probably have an improved gene expression network to cope with stresses they commonly face in their natural environment, such as high copper or sulfite concentrations for wine isolates, or lactose-rich environment for dairy isolates. This could also allow to better characterize the impact of the S. cerevisiae pangenome on gene expression, as several accessory genes are known to be advantageous in specific contexts such as vinification for HGT-related ORF in wine isolates (Marsit et al., 2015). Finally, performing GWAS on gene expression data coming from different culture conditions would be a promising way to explore in depth the overall genotype-phenotype relationship, especially when GWAS based on growth data was already performed using the 1,011 collection (Peter et al., 2018).

Transcript and peptide degradation

Gene expression is a complex mechanism where the final abundance of transcripts and peptides in a cell results from the combination of several factors, including the rate of synthesis, but also the rate of degradation (Buccitelli and Selbach, 2020). Therefore, both transcript and protein degradation can be considered as a determinant mechanism underlying gene expression, and their study could represent an interesting follow-up to this work. Protein degradation has been shown to be highly important for buffering mechanisms and partly underlies the phenomenon of post-transcriptional buffering (de Bie and Ciechanover, 2011;Gonçalves et al., 2017;Juszkiewicz and Hegde, 2018;Taggart et al., 2020). We have also shown that the rate of degradation has a significant influence on the overlap between the genetic origin of mRNA and protein abundance: proteins with high turnover tend to be more associated with proteomespecific regulation.

However, to date there is no high-throughput method for both transcript and protein degradation measurement, and current methods are still too laborious to be applied to more than a thousand samples.

Population-wide exploration of translation regulation

Thor-Ribo-Seq

The study of translation regulation is of particular interest as it represents the central link between transcriptome and proteome. Defective translation regulation can have phenotypic consequences, and has for example been shown to be implicated in the pathogenesis of many diseases such as cancer [START_REF] Robichaud | Translational Control in Cancer[END_REF]. However, compared to the transcriptome (Caudal et al., 2023;The GTEx Consortium, 2015) or proteome (Ferkingstad et al., 2021;Messner et al., 2023), the translatome has been poorly explored, especially at the population level where no large-scale studies have been conducted so far, most probably due to technical limitations.

Recently, a team developed a high-throughput approach for ribosome profiling called Thor-Ribo-Seq (Mito et al., 2023). This method relies on the use of a small amount of substrate by linear amplification of mRNA fragments covered by the ribosomes (Figure 1). This is a major advance in the field as it allows for ribosome profiling scalable to large number of samples.

This newly developed Thor-Ribo-Seq technique allows to consider generating the translatome for the whole collection of 1,000 natural isolates.

Expected insights

A large-scale ribosome profiling approach would help revealing several fundamental aspects of translation. First, the genetic origins of translation variation across individuals could be explored in depth and would help gain a more precise view of the mechanisms underlying translation regulation (e.g., local, or distant regulation, presence of regulatory hotspots).

Equally important will be the comparison of these results with those obtained on both transcriptomes and proteomes (Caudal et al., 2023;Muenzner et al., 2022). The comparative analysis of these datasets would represent an incredible opportunity to decipher the interactions between each expression layer and would provide an exhaustive view of the gene expression process. Finally, translation-specific regulation, such as ribosome velocity diversity and population-scale frameshift catalogs, will be achievable for the first time with this type of data.

Accounting for missing heritability

We performed genome-wide associations to identify the genetic origins of variation in both mRNA and protein abundance. Our study focused on SNPs and CNVs, yet other types of variants can impact phenotypes in general, and in gene expression in particular. Among them, structural variants (SVs) are central in modifying gene expression (Alonge et al., 2020). Their frequency and effects have already been studied in yeast (Dephoure et al., 2014;Muenzner et al., 2022;[START_REF] O'donnell | Telomere-to-telomere assemblies of 142 strains characterize the genome structural landscape in Saccharomyces cerevisiae[END_REF], and the largest catalog of SVs available for this species so far was constructed from a set of 142 natural isolates mostly out of the 1,011 population [START_REF] O'donnell | Telomere-to-telomere assemblies of 142 strains characterize the genome structural landscape in Saccharomyces cerevisiae[END_REF]. In this work, SVs were observed as impacting gene expression, especially near their breakpoints. This is in line with previous findings showing that SV-like With the Thor-Ribo-Seq method, inputs with low mRNA quantities are amplified via in vitro transcription after a fusion of the mRNA fragments with a T7 promoter. The resulting ribosome profiling data is highly similar to data obtained with standard ribosome profiling. Figure adapted from (Mito et al., 2023) Cell lysate with polysome RNA treatment

Linker ligation

In vitro transcription RT-PCR &

DNA sequencing T7 promoter inversions or translocations can directly affect the promoter of a gene and thus lead to a modification of gene expression, as has been shown in the case of sulfite resistance in some wine-related S. cerevisiae isolates (García-Ríos et al., 2019;[START_REF] Marullo | SSU1 Checkup, a Rapid Tool for Detecting Chromosomal Rearrangements Related to the SSU1 Promoter in Saccharomyces cerevisiae: An Ecological and Technological Study on Wine Yeast[END_REF]Pérez-Ortín et al., 2002;Yuasa et al., 2004;Zimmer et al., 2014). Yet, 142 isolates represent only a fraction of the total 1,011 population, suggesting that a large part of the SVs is certainly missed. A project aiming to sequence the entire 1,011-population with long-read sequencing method (Nanopore sequencing technology) is currently ongoing is our laboratory. As the experimental part is now completed, an exhaustive catalog of the SVs across the 1,011 isolates should be released within the next few months. The best methodology to perform GWAS on SVs is still debated, however pangenome graph-based association studies are a promising approach (He et al., 2023;Li et al., 2022;Logsdon et al., 2020;Zhou et al., 2022) Finally, rare variants (with a frequency in the population below 5%) have also been poorly considered in association studies so far (Génin, 2020;Manolio et al., 2009). The exploration of low-frequency variants in a natural population can be performed by artificially increasing the frequency of the variants using a diallel cross strategy (Fournier et al., 2019). A recent project in our lab investigated mRNA abundance within such a diallel cross and emphasized the significant impact of rare variants on gene expression [START_REF] Tsouris | Diallel panel reveals a significant impact of low-frequency genetic variants on gene expression variation in yeast[END_REF]. The phenotypic characterization of this population could for example be expanded to protein abundance which could help clarifying the effect of rare variants on the proteome.

Figure 1 :

 1 Figure 1 : (A) À l'aide de 8 souches naturelles de S. cerevisiae, nous avons quantifié l'abondance d'ARNm et de protéines et la traduction afin d'élucider comment l'expression de gènes varie à chaque étape. Ceci nous a permis de caractériser avec précision le phénomène d'atténuation post-transcriptionnelle et les contraintes évolutives spécifiques à chaque étape du processus d'expression de gènes. (B) Nous avons combiner des jeux de donnés de transcriptomes et de protéomes obtenus à très large échelle pour comparer les variations d'abondance d'ARNm et de protéines. Nous avons observé que ces variations étaient très différentes et que cela peut être expliqué par l'importante différence entre les origines génétique des abondances des ARNm et des protéines.

A

  larger coverage of the proteome .............................................................................................................. Exploration of other culture conditions ..................................................................................................... Transcript and peptide degradation ............................................................................................................

  Thor-Ribo-Seq ..................................................................................................................................................... Expected insights ............................................................................................................................................... Accounting for missing heritability ........................................................................................................... REFERENCES ....................................................................................................................................................... APPENDIX ................................................................................................................................................... LIST OF PUBLICATIONS ..........................................................................................................................................

  .. ORAL COMMUNICATIONS ......................................................................................................................................

  .. TEACHING ..........................................................................................................................................................

Figure 1 :

 1 Figure 1: The origin of the phenotypic landscape.Example of a (A) bimodal or (B) normal distribution of a trait in a population. While traits that follow bimodal distributions are often related to a single gene, the complex or polygenic traits follow normal distributions. (C) Dissection of the phenotypic variation origins of a trait

Figure 2 :

 2 Figure 2: The different types of genetic variants. (A) Example of indels and SNPs, with the consequences that can be induced by a SNP on coding sequences. (B) Schematic of structural variants (>50 bp), some of which can induce copy number variations (CNVs).

Figure 3 :

 3 Figure 3: Minor allele frequency distribution of the SNPs in a natural population of 1,011 S. cerevisiae isolates.

  03,0.04] (0.04,0.05] (0.05,0.1] (0.1,0.15] (0.15,0.2] (0.2,0.25] (0.25,0.3] (0.3,0.35] (0.35,0.4] (0.4,0.45] (0.45,0.5]

Figure 4 :

 4 Figure 4: Gene expression is a complex process that underlies the genotype-phenotype relationship.

Figure 5 :

 5 Figure 5: Short-and long-read mRNA sequencing Short-read RNA sequencing is based on the extraction and isolation/enrichment of mRNA, which is then fragmented, ligated to adapters, and sequenced with Illumina. The reads are typically mapped to the genome, and the number of reads per gene is normalized, e.g., using transcript per million (TPM) normalization. Long-read sequencing is based on either PacBio or Nanopore technologies. Both methods are capable of sequencing full-length mRNA and are therefore very useful for detecting isoform changes in transcripts, de novo transcripts or base modification.

Figure 6 :

 6 Figure 6: Recent advances in LC-MS/MS have enabled for high-throughput proteomic exploration. Several advances have increased the throughput of proteomic experiments. Laboratory robots such as plate rotators (Step 1) or automated liquid handlers (Step 2) have greatly facilitated the preparation of proteomic samples. Also, new methods of LC-MS/MS have led to a decrease in time per sample of proteomic experiment (step 3), approximately 19min/sample(Muenzner et al., 2022). Since the resulting dataset of large-scale proteome exploration are usually computationally intensive, new

Figure 7 :

 7 Figure 7: Ribosome profiling is a powerful tool to study translation. (A) Ribosome-profiling consists of the sequencing of mRNA fragments that have been protected from a nuclease treatment by the translating ribosomes. The reads have a classical length of 28nt and are mostly mapped to the ORF sequence (the 5' leader and the 3' UTR and usually not covered). (B) The ribosome-profiling reads usually show a 3nt periodicity, reflecting the codon-wise movement of the ribosome on the mRNA. (C) The density of the ribosome along the mRNA is a powerful revelator of translations dynamics, such as the elongation speed or translation pauses.

Figure 8 :

 8 Figure 8: Overview of the linkage mapping studies.

Figure 9 :

 9 Figure 9: Overview of GWAS.(A) GWAS are performed by combining genotyping data, where the allelic state of each individual is known, and phenotypic data for each individual in the population. The association between the phenotype and each variant is performed using an LMM that considers confounding factors, such as population structure. The final result of a GWAS is typically presented using Manhattan plots where the association of each SNP along the genome has an association score(-log10(p-value)) that must exceed a threshold (represented by the horizontal line) to be considered as significant. (B) Overview of the human GWAS where each colored dot corresponds to an association between a locus and a

Figure 10 :

 10 Figure 10: Expression-based GWAS are powerful tools for detecting the mechanisms of gene expression regulation. (A) Expression-based GWAS are typically represented by plotting the eQTL or pQTL position against the target gene position. The diagonal line (red dots) represents the eQTL or pQTL affecting the abundance of nearby genes. The accumulation of points on a precise vertical line reveals the presence of QTL hotspots that affect numerous genes across the genome. (B) The QTL hotspots are easily observed by looking at the number of genes affected by the QTL in a genetic window (here 25kb).Both plots are adapted fromCaudal et al., 2023.

Figure 11 :

 11 Figure 11: The genetic variants associated with changes in gene expression can affect the expression either locally or distantly. (A) Local regulatory variants are located close to their target gene, usually in the promoter region and thus alter the ability of various proteins, such as transcription factors, to bind to DNA. (B) Distant

  Saccharomyces cerevisiae. Several S. cerevisiae strains isolated from winemaking environment enhanced tolerance to sulfite(Pérez-Ortín et al., 2002), a common chemical used to suppress the growth of various non-Saccharomyces yeasts or lactic bacteria (Ribéreau-

Figure 12 :

 12 Figure 12: The across-gene correlation shows a good match between transcript and protein levels within a sample.(A) The across-gene correlation quantifies the correlation between transcript and protein abundance within a sample, such as a cell line, a tissue, or a strain (as shown here). This correlation is calculated using the overlap between protein and transcript measurements from the sample. (B) Over the past two decades, several measurements or the across-gene correlation have been made across multiple species and the calculated coefficients typically fall between 0.4 and 0.8. Data gathered fromBuccitelli and Selbach, 2020, and from Aydin et al., 2023.

Figures 13 :

 13 Figures 13: The within-gene correlation highlights how the protein abundance variations match the mRNA abundance variations.

Figure 14 :

 14 Figure 14: Post-transcriptional buffering is a central phenomenon in dealing with transcriptional variation.(A) Post-transcriptional buffering allows the cell to counteract the effect of transcriptional variation and results in a reduced variation at the translatome and proteome levels. This suggests that these two levels are under stronger evolutionary constraints than transcriptome. (B) The mechanisms underlying post-transcriptional buffering are still under debate. Autoregulation and stoichiometry regulation are thought to be part of these mechanisms. Adapted fromBuccitelli and Selbach, 2020.

Figure 15 :

 15 Figure 15: Neighbor-joining tree of the collection of 1011 S. cerevisiae isolates.

Figure 16 :

 16 Figure 16: Environmental origin of the 1011-population.The colors correspond to domesticated (shades of red), wild (shade of green) or clinical (blue) ecological environments. The domesticated environments comprise very different trophic and high condition diversity. Adapted fromPeter et al., 2018.

Figure 1 :

 1 Figure 1: Summary of the PhD project.

  abundance variations are well-known sources of phenotypic diversity across individuals. Protein abundance is influenced by both transcriptional and posttranscriptional regulations, which ultimately affect the final phenotypes. There are several cellular mechanisms involved in the modulation of final protein abundance, including mRNA stability, translation initiation and protein degradation(Buccitelli and Selbach, 2020). In the last decades, various technologies have greatly facilitated the detailed exploration of all these

Figure 1 .

 1 Figure 1. Exploration of the transcriptional and translational variations. (A) Correlation matrixes (Spearman correlation test) of each RNA-seq and Ribo-seq isolate pairwise comparison (all the coefficients displayed significant Pvalues). (B) Difference between the RNA-seq and Ribo-seq correlation levels (Wilcoxon test p-value = 1.4x10 -9 ). (C, D) Number of occurrences in the variable genes using (C) RNA-seq and (D) Ribo-seq data.

Figure 2 .

 2 Figure 2. The transcription variations are buffered because of TE modulation. (A) Euclidean distance comparison between RNA-seq or Ribo-seq profiles. The distances were significantly higher in RNA-seq data (Wilcoxon test p-value = 0.0478). (B) Differences between the RNA-seq and Ribo-seq |log2(FC)| values obtained by comparing each gene's TPM in each pairwise comparison (Wilcoxon test p-value = 3.3x10 -193 ). (C) Comparison between the RNA-seq log2(FC) and the TE log2(FC) in the AMH vs BAN isolate pairwise comparison (Spearman correlation p-value = 2.8x10 -317 ). (D) Comparison between the RNA-seq log2(FC) and the Ribo-seq log2(FC) in the AMH vs BAN isolate pairwise comparison (Spearman correlation p-value < 4.94x10 -324 )

Figure 3 .

 3 Figure 3. Post-transcriptional buffering has a specific signature. (A) Detection of the buffered and unbuffered genes according to a 1.5-fold change threshold (corresponding to the lines on the plot) in the AMH vs BAN isolate pairwise comparison, blue points = unbuffered variation; red points = buffered variation. (B, C) Proportion of, respectively, essential gene and protein complex related genes among the two gene groups in the 28 pairwise comparison. The proportion are in both cases higher in the buffered group (Wilcoxon test). (D) RNA-seq and (E) tAI levels of the recurrently buffered or unbuffered genes (respective Wilcoxon test p-value: 2.6x10 -67 and 1.03x10 -40 ).

Figure 4 .

 4 Figure 4. Translation levels of the S. cerevisiae pangenome. (A, B) RNA-seq and Ribo-seq level of the ORFs acquired through introgression event in the CPI isolate and being homozygous for the S. paradoxus allele (n = 18) and their orthologs in the other isolates. No difference in term of transcription of translation were observed between the introgression related ORFs and their orthologs. (C) TE difference between the ORFs acquired through HGT in the BPL isolate and the other ORFs (Wilcoxon test p-value = 0.0002).

Figure S2 .

 S2 Figure S2. Detection of the genes displaying variable transcription and translation regulation. BPL (x axis) vs CMP (y axis) pairwise comparison of the log10 transformed Ribo-seq data. (A) The points highlighted in red are the variable genes detected using Mahalanobis distance only. (B) The points highlighted in red are the genes detected using the combination of Mahalanobis distance and linear regression residuals.

Figure S3 .

 S3 Figure S3. Total variable genes detected using both RNA-seq and Ribo-seq data Variable genes detection in the 28 (A) RNA-seq and (B) Ribo-seq pairwise comparisons. The variable genes for RNA-seq are highlighted in green and the ones for Ribo-seq in red. The black lines correspond to the linear regression obtained in each pairwise comparison.

Figure S4 .Figure S5 .

 S4S5 Figure S4. Gene expression variance is higher at the transcriptional level. Gene-wise variance using the log10 transformed data of the two datasets (Wilcoxon test p-value = 5.09x10 -38 ).

Figure

  Figure S6. Number of buffered genes in each isolate pairwise comparison. Number of genes in the buffered group for each isolate (each isolate has 7 values corresponding to the 7 pairwise comparisons against the other isolates). The number of genes is slightly different across the isolates (Kruskal-Wallis test).

Figure S10 .Figure S12 .Figure S13 .

 S10S12S13 Figure S10. Comparison between the tAI and the expression levels.Correlation between tAI and log10 value of: (A) RNA-seq, (B) Ribo-seq and (C) TE data. The spearman correlation coefficient is above each plot, all of the coefficient where significant (P-value<0.05)

Figure 1 :

 1 Figure 1: Generation of the proteomes of 8 isolates of S. cerevisiae.(A) We generated the accurate proteome of 8 S. cerevisiae isolates from different ecological origins from mid-log phase culture in synthetic complete medium. We combined these data with RNA-seq and ribo-seq data generated on the same set of strains(Teyssonnière et al., submitted). (B) Enclosed proteins tend to be more transcribed than the rest of the proteins. (C) Correlation matrix between the isolates. Overall, the proteomic profiles tend to be very similar.

Figure 2 :

 2 Figure 2: DEP survey reveals specific signature related to environmental adaptation.

  by LIMMA. (B) The genes annotated as lipid metabolic pathway in the GO annotation tend to be enriched in the over-expressed genes in the BTT isolate. The blue line corresponds to the enrichment of lipid metabolic pathway genes along the Log2(FC) ranks. The red line represents the maximum enrichment score (as reported by the R package fgsea). The vertical segments represent the position of each lipid metabolic pathway gene along the Log2(FC). (C) Normalized enrichment scores (NES, high values correspond to enrichment among over-expressed proteins, while low values correspond to enrichment among under-expressed proteins) of significant respiration-related GO annotations (from left to right: ATP biosynthetic process; proton motive forcedriven ATP synthesis; mitochondrial translation; proton-transporting ATP synthase activity, rotational mechanism; proton transmembrane transport. See Methods) for each strain.

Figure 3 :

 3 Figure 3: Gene expression variations across the gene expression levels. (A) Pairwise correlations (n=28) of the 8 isolates in each dataset show that the proteomic profiles are more similar than those obtained on the transcriptome and translatome. All p-values were obtained from Wilcoxon tests and were lower than 1x10 -5 . (B) Comparison of |log2(FC)| (n= 79,520) obtained in each dataset revealed that transcriptional variations are buffered during the gene expression process. All p-values were obtained from Wilcoxon tests and were lower than 1x10 -20 . (C) The acrossgene correlation of translatome vs. proteome is only slightly higher than the cross-gene correlation of transcriptome vs. proteome. All p-values were obtained from paired Wilcoxon tests and were all equal to 0.0078. (D) The within-gene correlation of the translatome vs. proteome is not higher than the

Figure 4 :

 4 Figure 4: Exploration of evolutionary constraints on gene expression in S. cerevisiae.

  (D) The categories related to metabolism are associated with faster gene expression evolution in the translatome and proteome. The p-values are all computed with Wilcoxon test and those annotated are below 0.01. (E) The categories related to interactions are often associated with stronger constraints on the evolution of gene expression. The p-values are all calculated using the Wilcoxon test and are all below 0.01.

Figure S2 :

 S2 Figure S2: Quantification and imputation of protein expression across strains.(A) The density area represents the distribution of protein expression determined from the average logarithmic intensity (log10) of detected peptides for each strain. A compact representation is also shown as a one-dimensional marginal distribution horizontally below each plot. The number of proteins quantified in each strain is written on the top left corner. A black density line indicates the distribution of protein expression after the imputation of missing values for proteins that were partially undetected among samples. The number of proteins partially detected is indicated in the top right corner. (B) The position of the circles relative to the x-axis is given by the average difference between imputed and quantified protein expression (Δ intensity) in each strain. For each strain, two arms extend in opposite directions to highlight the range of Δ intensity within one standard deviation from the mean represented by the circles.

Figure S4 :

 S4 Figure S4: Inter and intra protein abundance variation.(A) Heatmap of the Pearson correlation coefficients for pairwise comparisons of protein expression between proteomics samples. A hierarchical clustering based on the complete Euclidean distances between samples' expression profiles was applied to both rows and columns. The colors represent each isolate as shown in figure1. (B) Principal component analysis of samples expression profile. The expression profiles of all samples are plotted as points scatterplot against the first (xaxis) and second (y-axis) principal components, which capture 52% of the variability between them. Points are colored according to their matching strains. The shaded regions delimit the range of variability among samples of the same strain.

Figure S5 :

 S5 Figure S5: Distribution of the CVs across from each quantified protein.Histogram of the CV calculated on the median protein abundance of each of the 3,429 genes quantified in the proteomic data.

Figure S6 :

 S6 Figure S6: LOESS normalization of the 3 datasets. (A) Overview of the abundance values across the 8 isolates in the 3 data before normalization. (B) Overview of the abundance values across the 8 isolates in the 3 data after normalization.

Figure S8 :

 S8 Figure S8: Number of differentially expressed proteins for each isolate.Differential protein expression was calculated using LIMMA (see Methods) comparing the average values across samples for one strain versus the average across samples for the remaining strains (onevs-all) for every protein. We consider the decrease/increase of protein expression significant when the absolute fold-change reached at least 1.2 and if the adjusted p-value was greater than 10 -5 , at a false discovery rate of 5%. The number of differentially expressed proteins (DEP) is reported along the x-axis, with underexpressed and overexpressed proteins respectively shown as bars in the left and right directions. The total number of differentially expressed proteins for each strain is written on the right. The number of differentially expressed proteins using only quantified protein expression (i.e. non-imputed) is also shown in lighter colors and represented by dimmed colored bars.

Figure S9 :Figure S10 :

 S9S10 Figure S9: Overview of the BP GO features enriched in the DEPGraphical representation of the GO enrichment analysis on the DEP. The main GO categories (in white) were obtained using semantic similarities on the terms detected using the gprofiler2 package in R. The semantic similarity was performed using the rrvgo package in R.

Figure S11 :

 S11 Figure S11: Across and within gene correlation using the normalized data. (A) The across-gene correlation of translatome vs. proteome is only slightly higher than the crossgene correlation of transcriptome vs. proteome. The p-values were obtained from paired Wilcoxon tests and the ** correspond to p-values equal to 0.0078. (B) The within-gene correlation of the translatome vs. proteome is not higher than the across-gene correlation of the transcriptome vs. proteome. All p-values were obtained by paired Wilcoxon test. The *** correspond to p-values less than 1x10 -20 .

Figure S12 :

 S12 Figure S12: Effect of the gene number on the tree length in each dataset.For each gene number, 1,000 trees were constructed by randomly selecting the corresponding number of genes and the total length was computed.

Figure 1 .

 1 Figure 1. Quantitative proteomes and transcriptomes of a large S. cerevisiae population.

Figure 3 .

 3 Figure 3. Co-expression network is a major determinant of the proteome organization while the population structure is not. A-B-C. Comparison between the phylogenetic tree (A) obtained using the bi-allelic SNP (as in Peter et al. 2018) and the trees obtained from the Euclidean distances based on protein (B) or transcript (C) abundance. Colors correspond to the subpopulations. D-E. Cellular co-expression network computed with WGCNA using proteomic (D) or transcriptomic (E) data. Colors represent the cellular pathway detected for each co-expression module.

Figure 4 .

 4 Figure 4. DEPs reveal domestication-and subpopulation-specific metabolic adaptation. A. GSEA results on the DEPs (using 16 broad functional annotations from 44) of each subpopulation. Colors represent the normalized enrichment score (NES): Red -overexpression, blue -underexpression in subpopulation. B. Volcano plot of the comparison between wild and domesticated isolates. Colors highlight the genes belonging to two functional annotations related to chaperon mediated folding and respiration.

Figure 5 .

 5 Figure 5. SNP-and CNV-pQTL detection highlights strong differences in the genetic origin of transcript and protein abundance. A-B. Map of the SNP-(A) and CNV-(B) pQTL. The x-axis is the QTL position on the genome and the y-axis the position of the affected gene on the genome. The x and y-axis numbers represent the 16

Figure S1 .

 S1 Figure S1. Description of the characteristics and the normalization of the population proteome.

Figure S3 .Figure S5 .

 S3S5 Figure S2. Across-gene correlation. mRNA-protein correlation in each isolate (across-gene correlation). The blue line represents the median (0.53).

Figure S6 .

 S6 Figure S6. Functional exploration of the transcriptome WGCNA modules.Functional enrichment of each co-expression module detected using WGCNA on transcript abundance data. The enrichment was performed using the CEMiTool package. The dotted lines on each graph represent the significance threshold.

Figure S9 .

 S9 Figure S9. Location of the local SNP-pQTL.Distribution of the local SNP-pQTL around the start codon of their target gene. Downstream pQTL correspond to QTL located between the stop codon and 200 bp after the stop codon, upstream correspond to pQTL located between the start codon and 1,000 bp before the start codon.

Figure S10 .Figure S11 .Figure S12 .

 S10S11S12 Figure S10. Aneuploidy related CNV-pQTL have a higher effect-size than the other CNV-pQTL. Difference in effect size between the aneuploidy related CNV-pQTL and the other CNV-pQTL (*** = p-value < 2.2x10 -16 , Wilcoxon test).

  Figure S14. CV from the QCs and samples precursors.The CV was computing either the QCs set or the sample set.

Figure 1 :

 1 Figure 1: Thor-Ribo-Seq is a powerful approach to perform ribosome profiling on low input sample.

  have been developed to efficiently and accurately handle and process LC-MS/MS data (step 4). Figure adapted from Muenzner et al., 2022, images from Biorender and Singer websites.

mRNA abundance and tAI difference between the buffered gene the rest of the genes.

  Difference in the (A) mRNA levels and (B) tAI levels of the buffered genes against all the other genes. Respective Wilcoxon test p-values: 6.15x10 -33 and 2.42x10 -22 .

			Kruskal-Wallis, p = 0.032 Kruskal-Wallis, p = 0.093
	1300
	1000 1100 1200 900 Figure S8. Figure S9. % of Buffered genes Number of Buffered genes	BPL BPL	CMP CMP	CQC CQC	AMH AMH	CPI CPI	BTT BTT	BED BED	BAN BAN

S6. Number of buffered genes in each isolate pairwise comparison. Number of genes in the buffered group for each isolate (each isolate has 7 values corresponding to the 7 pairwise comparisons against the other isolates). The number of genes is slightly different across the isolates (Kruskal-Wallis test).

Figure S7. Proportion of buffered genes in each isolate pairwise comparison.

Percentage of gene having their variations buffered among the ones that displayed divergent transcription for each isolate (each isolate has 7 percentage values corresponding to the 7 pairwise comparisons against the other isolates). There is no difference in term of percentage among the isolates (Kruskal-Wallis test).

RNA-seq and Ribo-seq levels of the buffered and unbuffered genes in each pairwise comparison

	log 10 (RNA-seq TPM)	2.0 2.2 2.1	RNASeq RNASeq	50 500 10000 5 50 500 5000	Wilcoxon, p = 1e-15 BPL CMP Rho = 0.612 Rho = 0.464 50 500 10000 RNASeq 0.2 0.4 tAI 0.2 0.4 0.6 5 tAI CPI BTT Rho = 0.609 RNASeq 50 500 10000 Rho = 0.622	10000 10 100 log RNASeq 1 10 (Ribo-seq TPM) 50 500 5000 RNASeq	CQC Rho = 0.573 2.2 1.8 2.0 0.2 0.4 0.6 tAI BED Rho = 0.539	RNASeq RNASeq	Wilcoxon, p = 2.6e-16 AMH Rho = 0.62 10000 10 100 0.2 0.4 0.6 1 tAI BAN 50 500 5000 Rho = 0.532
				5					5				5
		1.9											5
					0.2	0.4	0.6			0.2	0.4	0.6	0.2	0.4	0.6	0.2	0.4	0.6
						tAI					tAI		tAI	tAI
						BPL					CMP		1.6	CQC	AMH
			RiboSeq	Buffered n = 28 1 10 100 10000 Rho = 0.686	Unbuffered n = 28 1 10 100 10000 RiboSeq Rho = 0.57	RiboSeq	1 10 100 10000	Buffered n = 28 1 10 RiboSeq 1000 Rho = 0.623	Unbuffered n = 28 Rho = 0.708
					0.2	0.4				0.2	0.4	0.6	0.2	0.4	0.6	0.2	0.4	0.6
						tAI					tAI		tAI	tAI
						CPI					BTT		BED	BAN
			RiboSeq	10 100 10000	Rho = 0.668	RiboSeq	10 100 10000	Rho = 0.67	RiboSeq	1 10 1000	Rho = 0.597	RiboSeq	5 50 500 10000	Rho = 0.624
									1
				1
					0.2	0.4	0.6			0.2	0.4	0.6	0.2	0.4	0.6	0.2	0.4	0.6
						tAI					tAI		tAI	tAI
						BPL					CMP		CQC	AMH
			TE	0.02 0.20 2.00 20.00	Rho = 0.373	TE	0.02 0.20 2.00 20.00	Rho = 0.368	TE	0.05 0.50 5.00	Rho = 0.352	TE	0.05 0.50 5.00	Rho = 0.408
					0.2	0.4				0.2	0.4	0.6	0.2	0.4	0.6	0.2	0.4	0.6
						tAI					tAI		tAI	tAI
						CPI					BTT		BED
			TE	0.02 0.20 2.00 20.00	Rho = 0.334	TE	0.05 0.50 5.00	Rho = 0.324	TE	0.02 0.20 2.00 20.00	Rho = 0.337	0.02 0.20 2.00 20.00
					0.2	0.4	0.6			0.2	0.4	0.6	0.2	0.4	0.6	0.2	0.4	0.6
						tAI					tAI		tAI

. Difference between the average expression level of each gene group (buffered and unbuffered) in each pairwise comparison using both RNA-seq and Ribo-seq data.

DEP detection in each isolate.

	A
	20
	15
	10
	5
	0
	CPI_Prot	CMP_Prot	AMH_Prot	CQC_Prot	BPL_Prot	BTT_Prot	BED_Prot	BAN_Prot	BPL_RNA	CMP_RNA	CQC_RNA	AMH_RNA	CPI_RNA	BTT_RNA	BED_RNA	BAN_RNA	BPL_Ribo	CMP_Ribo	CQC_Ribo	AMH_Ribo	CPI_Ribo	BTT_Ribo	BED_Ribo	BAN_Ribo
	B
	15
	10
	Figure S7:
	5
	0
	CPI_Prot	CMP_Prot	AMH_Prot	CQC_Prot	BPL_Prot	BTT_Prot	BED_Prot	BAN_Prot	BPL_RNA	CMP_RNA	CQC_RNA	AMH_RNA	CPI_RNA	BTT_RNA	BED_RNA	BAN_RNA	BPL_Ribo	CMP_Ribo	CQC_Ribo	AMH_Ribo	CPI_Ribo	BTT_Ribo	BED_Ribo	BAN_Ribo

  Precursors were filtered according to their detection rate in the samples, with a threshold set at 80% of detection rate across all the strains, while precursors with a coefficient of variation (CV) above 0.3 in the QC samples were excluded. The CVs of QCs and wild isolates samples were calculated and had a median CV of 15.15% and 34.21%, respectively (FigureS14; table

	Samples and entries with insufficient MS2 signal quality (< 1/3 of median MS2 signal) and
	with entries with Q.Value (> 0.01), PG.Q.Value (> 0.01), Global.Protein.Q.Value (> 0.01),
	Global.PG.Q.Value (> 0.01) were removed. A similar threshold was applied to Lib.PG.Q.Value
	and Lib.Q.Value to account for the MBR option used. Non-proteotypic precursors were also

excluded. Outlier samples were detected based on the total ion chromatograms (TIC) and number of identified precursors per sample (Z-Score > 2.5) and were excluded from further analysis.

phenotype. We also sought to detect QTL hotspots in both transcriptome and proteome GWAS.

We defined a hotspot as a concentration of at least 4 QTL in a 20 kb window.

We compared the protein turnover rate (Muenzner et al., 2022) of obtained on 619 proteins encompassed in our dataset to see whether turnover rate had an impact on the overlap between SNP-eQTL and SNP-pQTL. This data comprises protein degradation rates for 1,836 gene across 55 natural isolates. We computed an average turnover rate per gene and used this value to compare the level of protein degradation of the protein with or without an overlapping QTL.

Data availability

All sequencing reads are available in the Gene Expression Omnibus (GEO) under the accession number GSE173654. https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE173654

Overview of the project

Unraveling the genetic origins of the large phenotypic diversity observed in nature is a central goal of modern biology. In the case of complex traits such as cancer, height or autism, the association between genetic variants and the phenotype is often tedious to achieve, and overall, the mechanisms linking the variant to its cognate phenotype are still elusive. Gene expression is one key determinant of the genotype-phenotype relationship. More specifically, the regulation of gene expression plays a major role in translating genotypes into phenotypes, in particular when the genetic variants associated to complex traits are localized in non-coding or regulatory genomic regions. However, this regulation is truly puzzling because it involves regulatory processes that affect each level of gene expression. Over the past two decades, several technological and analytical advances have greatly facilitated the study of gene expression. Quantification of gene expression using RNA-seq, ribosome profiling as well as LC-MS/MS allowed for accurate quantification of each step of the process, while GWAS or linkage mapping helped to precisely map the origin of transcript or protein abundance variation on the genome. Nevertheless, large-scale studies of gene expression are sparse, and many aspects of gene expression remain to be explored and elucidated. For example, the similarities between the genetic origins of mRNA or protein abundance, and more globally, how protein and mRNA variation fit together, are still ongoing debates.

In this context, my project aims to take advantage of these technological advances and of the powerful S. cerevisiae model to study the variation of gene expression among individuals at the population level. A large collection of natural isolates of S. cerevisiae is available in our laboratory, which gathers more than a thousand strains for which genomes were completely sequenced using Illumina technology. The collection show very diverse ecological origins and includes both domesticated and wild strains, resulting in an accurate representation of the species diversity.

During these 4 years, I explored each level of gene expression in order to investigate their variation in a natural population. In a first chapter, I will describe the survey of gene expression at the transcriptional and translational levels, with the aim of characterizing a known determinant of variation between individuals, namely post-transcriptional buffering (Figure 1A). This phenomenon describes the fact that transcriptional variation tends to be buffered for as long as the expression process, suggesting increased evolutionary constraints on the later steps of gene expression. Although frequently observed, this phenomenon remains poorly understood, particularly at the level of translation. In a collaboration with the Riken Institute

CHAPTER I Translation variation across genetic backgrounds reveals a post-transcriptional buffering signature in yeast

Abstract

Gene expression is known to vary among individuals, and this variability can impact the phenotypic diversity observed in natural populations. While the transcriptome and proteome have been extensively studied, little is known about the translation process itself. Here, we therefore performed ribosome and transcriptomic profiling on a genetically and ecologically diverse subset of natural isolates of the Saccharomyces cerevisiae yeast. Interestingly, we found that the Euclidean distances between each profile and the expression fold changes in each pairwise isolate comparison were approximately 10% higher at the transcriptomic level. This observation clearly indicates that the transcriptional variation observed in the different isolates is buffered through a phenomenon known as post-transcriptional buffering at the translation level. Furthermore, this phenomenon seemed to have a specific signature by preferentially affecting essential genes as well as genes involved in complex-forming proteins, and low transcribed genes. We also explored the translation of the S. cerevisiae pangenome and found that the accessory genes related to introgression events displayed similar transcription and translation levels as the core genome. By contrast, genes acquired through horizontal gene transfer events tended to be less efficiently translated. Together, our results highlight both the extent and signature of the post-transcriptional buffering.

displaying 2.5% highest and lowest residuals. We kept the genes that were overlapping between the Mahalanobis detection method and the linear regression residual method. The detected genes are later mentioned as "variable genes".

Variable genes characteristics

Using genes descriptive data (Dowell et al., 2010;Giaever et al., 2002;Pu et al., 2009), we focused on describing the characteristics of the variable genes. Firstly, we questioned if the variable genes detected earlier displayed any enrichment or depletion of essential genes or genes part of protein complexes with Fisher's exact tests (FET). Gene ontology (GO) analysis (Ashburner et al., 2000) was performed on the geneontology.org website (Gene Ontology Consortium, 2021; Gene Ontology Consortium, 2019), using the subset of genes (N = 3755) encompassed by our Ribo-seq and RNA-seq experiments as the reference list. The p-values were corrected using Bonferroni correction. This was performed on RNA-seq or Ribo-seq variable genes, and on the overlapping variable genes (between the two datasets).

Detection of a post-transcriptional buffering phenomenon

In order to see if expression variation across the 8 isolates was stronger in RNA-seq or Riboseq data, we firstly computed the Euclidean distances between each strain using the log10transformed data from both datasets. We also obtained for each gene the expression variance using the log10-transformed data from both datasets. In addition, we generated log2 foldchange (log2(FC)) values in RNA-seq, Ribo-seq and TE datasets, where for each gene in each pairwise comparison:

Then, in each pairwise comparison, we used Spearman correlation test to compare the behavior of the 3 kinds of log2(FC) dataset (RNA-seq, Ribo-seq, and TE) against each other. We finally quantified the expression variation (in both RNA-seq and Ribo-seq data) based on the log2(FC):

The more this value increases, the more the difference between the TPM values is important.

Buffered and conserved regulation genes characteristics

In each isolate pairwise comparison and using the RNA-seq log2(FC) vs. TE log2(FC) comparison, we defined 2 groups of genes by applying a 1.5-foldchange threshold (log2(FC)≈0.58 or log2(FC)≈-0.58) for both RNA-seq log2(FC) vs. TE log2(FC) (see Figure 3A). This enables to capture genes with a buffered transcriptional variation (in red in Figure 3A) among the genes displaying at least a 1.5 foldchange transcriptional variation. The gene displaying transcriptional variation that were not capture among the buffered genes were considered as unbuffered genes (in blue in Figure 3A). We checked the percentage of genes concerned by post-transcriptional buffering among the genes that displayed transcriptional variation (using a minimal 1.5-foldchange threshold in the pairwise comparison). Genes were considered as recurrently buffered or unbuffered if they were detected in the corresponding group at least in more than half of the isolate pairwise comparison (i.e., detected 15 or more as buffered or unbuffered). We then performed Fisher's exact test with the two gene sets (buffered or unbuffered) to detect enrichment of essential or protein complex-related genes. We also check the proportion of essential genes and protein complex-related genes in the two groups in each pairwise comparison.

We then compare the transcription or translation level of the buffered and unbuffered groups by using for each gene the mean RNA-seq TPM value across the 8 isolates.

Codon usage bias influence

We used the tAI (tRNA adaptation index) index [START_REF] Dos Reis | Solving the riddle of codon usage preferences: a test for translational selection[END_REF][START_REF] Dos Reis | Unexpected correlations between gene expression and codon usage bias from microarray data for the whole Escherichia coli K-12 genome[END_REF] to estimate the codon bias usage of each gene. Briefly, tAI is an index showing how much a gene is adapted to the tRNA genome structure in terms of codon usage. In this perspective, we first calculated the tRNA copy number of our 8 strains with tRNAscan-SE with default parameters (Chan and [START_REF] Caudal | Pan-transcriptome reveals a large accessory genome contribution to gene expression variation in yeast[END_REF][START_REF] Lowe | tRNAscan-SE On-line: integrating search and context for analysis of transfer RNA genes[END_REF] using the assembled genome sequences from (Peter et al., 2018). Then, for each isolate, we used the tRNA copy number to compute the tAI using a Perl program available on https://github.com/mariodosreis/tai [START_REF] Dos Reis | Solving the riddle of codon usage preferences: a test for translational selection[END_REF][START_REF] Dos Reis | Unexpected correlations between gene expression and codon usage bias from microarray data for the whole Escherichia coli K-12 genome[END_REF] using default parameters. The resulting dataset was a tAI value for more than 98.8% of the 3,755 genes (some genes were discarded during the calculation) in each isolate. Ultimately, we could calculate an overall tAI (mean of the 8 or less tAI values) for 3,746 genes. For each isolate, we correlated the expression levels and TE of each gene to its tAI index.

Accessory ORF analysis

Using the mapping done on the S. cerevisiae pangenome, we selected the accessory ORFs that were previously detected in each of our 8 strains and selected the ones that had TPM values higher than 0 in both RNA-seq and Ribo-seq data. Then, we calculated the TE of each accessory ORF.

We first focused on the CPI isolate ORF acquired through introgression events (with the Saccharomyces paradoxus species). We selected the ORFs known to have an ortholog in S. cerevisiae genes. Then using homo/heterozygosity data (adapted from (Peter et al., 2018) gene presence/absence data), we selected the ORF that were homozygous for the S. paradoxus allele (n = 18) and we compared their expression with their orthologs in S. cerevisiae. We then compared the introgression TPM values vs. their orthologs mean TPM value (obtained from the 7 other strains) to see if we could observe over-or under-expression of these ORFs in comparison with their orthologs.

We explored then the expression levels of the BPL isolate accessory ORFs, especially the ones acquired through horizontal gene transfer (HGT) by comparing their expressions (RNA-seq and Ribo-seq) and TE with the other BPL gene values.

Supplementary Material

Supplementary tables available at: https://www.dropbox.com/scl/fi/782fpyfyarxkbwqas8494/Sup_Tables.xlsx?rlkey=2oztn6rye7 1qy3rqgotq91nwv&dl=0 Figure S1. Eight isolates were selected to ensure a maximum genetic diversity Neighbor-joining tree obtained from the biallelic SNPs of 1011 isolates (1). The colors of the points correspond to their clade, and the 8 isolates are highlighted in red.

CHAPTER II

Metabolism adaptation is a main driver of protein abundance evolution in the yeast Saccharomyces cerevisiae.

Abstract

Variation in gene expression among individuals is one of the major causes of the phenotypic diversity observed in natural populations. Expression variation can occur at any step of the gene expression process, and numerous human diseases have been linked to transcriptional or posttranscriptional variation. However, the mechanisms that influence the evolution of each level of gene expression remain poorly understood. Here, we monitored the proteome of 8 natural isolates of Saccharomyces cerevisiae and compared it with the transcriptome and translatome previously obtained from the same set of strains. We found that the proteome variations were mainly related to metabolism or respiration related genes. Interestingly, we found that the proteome variations differed from the transcriptome and proteome variations in part due to posttranscriptional buffering. In addition, we observed that translational variation (measured by ribosome profiling) was only slightly better at reflecting variation in protein abundance than transcriptional variation. Finally, we examined the factors that influence gene expression evolution at each gene expression level. We found that despite expression level specificities, similar evolutionary constraints affect all steps of gene expression. For example, genes encoding highly interacting proteins showed more conserved gene expression regulation, while metabolism-related genes showed faster gene expression evolution. Our results highlight that adaptation to different trophic conditions is a major driver of gene expression evolution across individuals.

Results

Proteome quantification of 8 S. cerevisiae isolates

We quantitatively profiled the proteomes of eight isolated strains of S. cerevisiae (Table S1)

enabling the detection of 28,800 (±3000) peptides per strain, including 93.8% of unique peptides (Figure S1). Our protein identification against the S. cerevisiae S288C reference proteome clustered these peptides in 3635 protein groups. Before analysis, we removed 20 hits corresponding to common lab contaminants, 52 hits corresponding to reverse sequences and 77 hits matching multiple-proteins (mostly duplicated proteins). For quantitative analysis, we only consider proteins identified by at least two unique peptides yielding 3,429 single-protein hits (Figure 1A, Table S2). The protein levels were quantified based on average peptide intensities determined by label-free quantification method (Cox et al., 2014). Therefore, the profiled intracellular proteome covers about half of the reference S288C S. cerevisiae proteome, and accounts for ~70% of cytoplasmic proteins in budding yeast. A significant proportion of the detected proteins (66%, n=2,280) were ubiquitously expressed across all eight strains while the remaining hits (34%, n=1,149) had partially missing protein levels within samples. Among those hits, a majority (73%, n=837) had been quantified in at least one replicate of every strain, and more than half (57%, n=660) in at least two replicates of every strain. To avoid discarding potentially valuable data, we decided to impute the missing protein intensities based on the distribution of quantified hits. Thus, 8% of protein intensities were imputed among all samples (ranging from 5% for AMH to 12% for CQC) for 10-30% of proteins per strain where at least one sample had a missing value (Figure S2A). For every strain, the variation of protein expression due to imputation should be negligible (below 2-fold between the 25th-75th percentiles) and at most comprised between 2 to 4-fold (5th-95th percentiles of Δ intensity), except for CQC where this difference rises to about 10-fold. Nevertheless, the difference of intensity (Δ intensity) between observed and imputed expression was centered on 0 for most strains, and slightly shifted for CQC (-0.12) and CMP (-0.04) towards lower expression after imputation (Figure S2B).

Overall, the 3,429 proteins encompassed genes with high transcripts level (Figure 1B) and a functional enrichment using GO annotation (Ashburner et al., 2000;Gene Ontology Consortium, 2021) revealed that several groups of genes tended to be more captured by our proteome exploration: genes related to protein transport, translation, ribosome, transcription and finally metabolism paths were overrepresented among the 3,429 proteins (Figure S3, table S3). Additionally, we observed an enrichment of essential genes and genes related to protein complexes among the captured proteins. Overall, the proteome profiles were very similar among the 8 isolates as indicated by the strong correlation of their average protein expression (Rho>0.9, Figure 1C), with BPL displaying a slightly different profile from the others.

Furthermore, protein levels are considerably stable within strains compared to their variation across different strains. Indeed, we note a very weak variation between biological and technical replicates among each strain (R>0.97), except for CMP and CQC (Figure S4A).

The CMP biological replicates are slightly less correlated than their technical replicates (R bio>0.95 vs. R tech>0.97). For CQC, the correlation between biological replicates is also not as strong (0.86>R bio>0.90) with one pair of technical replicates being quite reproducible (R tech.1=0.97) unlike the other pair (R tech.2=0.93). The principal component analysis captures 53% of the variability in expression among samples (Figure S4B) with noticeable higher biological variation for CMP and CQC strains as seen by the greater distance between their samples compared to the other strains. In the case of CQC, we can also observe a higher distance among the technical replicates, indicating their lack of reproducibility.

In addition, we also measured the variation of expression at the gene-level. We calculated the CV (measured as a percentage) for each gene and found that it ranged from 3 to 282%, with a median of 27% (Figure S5). Using gene set enrichment analysis (GSEA) with the R package fgsea (Korotkevich et al., 2021;Subramanian et al., 2005), we found that metabolism-related genes were enriched among the most variable genes (Table S4).

Recently, RNA sequencing and ribosome profiling were performed on the above 8 isolates (Teyssonnière et al., submitted). We sought to compare the transcriptomes, translatomes, and proteomes of these 8 isolates, which overlapped on 2,840 genes. Using the mean protein abundance for each gene across the 4 replicates and by combining it with the RNA-seq and ribo-seq data, we performed LOESS normalization of the 3 datasets together to accurately analyze the difference in abundance variation across expression levels (Figure S6, Table S5).

this phenomenon takes place in each post-transcriptional step and therefore, different mechanisms are certainly involved in its establishment.

Interestingly, our results suggest that post-transcriptional buffering affects the proteome more than the translatome. Overall, we observed that proteome variation was consistently lower than translatome variation. As ribosome profiling has been considered as a proxy for protein abundance (Brar and Weissman, 2015), we sought to question the similarities between the translatome and the proteome and the ability of ribosome profiling to reflect protein abundance.

In this context, we explore the correlations between each expression layer. To explore the relationship between the proteome and the other gene expression layers, two types of correlations are typically calculated: the across-gene correlation and the within-gene correlation (Buccitelli and Selbach, 2020;Liu et al., 2016). While the across-gene correlation explores the relationship between expression levels in a sample (here, an isolate), the within-gene correlation compares how each expression level varies in each gene and provides a better view of the similarities in variation between the transcriptome, translatome, and proteome. Computed with the non-normalized data, the across-gene correlation revealed that the mRNA-protein correlations were on average reaching 0.53 (Spearman correlation tests, Figure 3C) which is in line with previous explorations (Buccitelli and Selbach, 2020). Surprisingly, the translatomeproteome correlation was only a little higher than this: 0.59 (Figure 3C). The difference was nonetheless significant (paired Wilcoxon test p-value = 0.0078). Yet, this is greatly lower than the transcriptome-translatome correlation that reached 0.83. This suggests that although ribosome profiling is a better proxy for protein abundance than RNA sequencing, it is very limited as a predictor of protein abundance within an individual. Using the normalized data, this was even more apparent, as no difference was observed between the mRNA-protein correlation and the translatome proteome correlation (Figure S11A). For the within-gene correlation, even though the comparisons include 8 samples for each gene, the large number of genes compensates for this and allows us to compare variation differences between expression levels. Here, we found no difference between the predictability of protein abundance across the 8 isolates using either transcriptome or translatome data (Figure 3D) The transcriptomeproteome correlation and the translatome-proteome correlation had an average within-gene correlation of 0.11 and 0.11, respectively, while the transcriptome-translatome correlation reached 0.52. The results were identical when using the normalized data (Figure S11B). This highlights that ribosome profiling data may not be a reliable proxy for protein abundance across

Materials and methods

Sample preparation for proteomics profiling

In this study, we conducted a comprehensive proteomics profiling experiment using a subset of eight S. cerevisiae strains representing the diverse range of ecological, geographical, and genetic characteristics from a population of 1,011 natural isolates. The selected strains were cultured on synthetic defined media (SD), and their growth was closely monitored by measuring the optical density (OD). We specifically harvested cells when they reached the mid-log phase (OD ~0.5). Prior to sample processing, we performed two rounds of washing using Phosphate-Buffered Saline solution (PBS) to remove extraneous contaminants. The cell pellets were then flash-frozen in liquid nitrogen. In total, we prepared 32 samples, with OD values ranging from 0.4 to 0.8 units. Four replicates were prepared for each strain, consisting of two biological replicates derived from distinct colonies and two technical replicates that underwent identical sample preparation, only separating them before freezing them. Subsequently, all samples were sent to the proteomics facility for further analysis. The frozen cell pellets were then lysed and submitted to in-solution tryptic digestion using the S-Trap method (by Protifi). A solid phase extraction cleaning step using Oasis HLB was employed to purify the resulting peptides. The purified peptides were then subjected to nanoflow liquid chromatography (nanoAcquity) coupled with high-resolution, high-mass accuracy mass spectrometry (Thermo Exploris 480).

Proteomics identification and database searching

For data analysis, each sample was analyzed separately on the mass spectrometer in a randomized order during the discovery mode. The raw data acquired from the instrument were processed using MaxQuant v1.6.6.0. The Andromeda search engine was utilized to search the data against a database comprising protein sequences of Saccharomyces cerevisiae obtained from Uniprot.org. This database was supplemented with common lab protein contaminants.

During the search, we considered fixed modifications such as cysteine carbamidomethylation and variable modifications of methionine oxidation and/or protein N-terminal acetylation.

Quantitative comparisons were performed using Perseus v1.6.0.7. Decoy hits were filtered out, and only proteins detected in at least two replicates of at least one experimental group were retained for further analysis. This rigorous methodology ensured high-quality data for subsequent interpretation and downstream analysis. Ultimately, our proteomic dataset comprised 3,429 genes. We performed a GO analysis on this set of genes using the R package

Supplementary Material

Supplementary tables available at: https://www.dropbox.com/scl/fi/dljfmouqew8z7gkeb3dag/S_table.xlsx?rlkey=2v1slof41uecv gl21t5sguxkx&dl=0

CHAPTER III

Species-wide quantitative transcriptomes and proteomes reveal distinct genetic control of gene expression variation in yeast

Summary

Gene expression varies between individuals and corresponds to a key step linking genotypes to phenotypes. However, our knowledge regarding the species-wide genetic control of protein abundance, including its dependency on transcript levels, is very limited. Here, we have determined quantitative proteomes of a large population of 942 diverse natural Saccharomyces cerevisiae yeast isolates. We found that mRNA and protein abundances are weakly correlated at the population gene level. While the protein co-expression network recapitulates major biological functions, differential expression patterns reveal proteomic signatures related to specific populations. Comprehensive genetic association analyses highlight that genetic variants associated with variation in protein (pQTL) and transcript (eQTL) levels poorly overlap (3.6%). Our results demonstrate that transcriptome and proteome are governed by distinct genetic bases, likely explained by protein turnover. It also highlights the importance of integrating these different levels of gene expression to better understand the genotypephenotype relationship.

Results

Quantitative proteomes of a large collection of natural isolates

We generated a quantitative proteomic dataset for strains of the 1,011 strains collection (Peter et al., 2018) from cells cultivated in synthetic complete medium with amino acids in order to match the growth medium used for RNA sequencing (Caudal et al., 2023) (Figure 1A).We had previously acquired a proteome dataset of the 1,011 strains collection, measured with microflow chromatography and SWATH MS (Muenzner et al., 2022). For the acquisition of this new dataset we used a proteomic method that allows for an even higher throughput, using analytical flowrate chromatography and Scanning-SWATH MS with a 3 min gradient [START_REF] Messner | Ultra-fast proteomics with Scanning SWATH[END_REF]. After cultivation of the yeast isolates in 96 well plates, proteins were extracted, and subjected to reduction, alkylation, and trypsination in a semi-automated workflow using liquid handling robotics (Messner et al., 2020). Peptide preparations were separated using a 3-minute high-flow rate (800 µl/min) chromatographic gradient using an Infinity II chromatographic system (Agilent Technologies), coupled to a 6600 Triple TOF instrument (Sciex). Data was recorded using Scanning SWATH acquisition [START_REF] Messner | Ultra-fast proteomics with Scanning SWATH[END_REF] and the raw data was processed using the DIA-NN software (version 1.8), which was specifically developed for large scale proteomic exploration (Demichev et al., 2020). We applied several quality filters where poor-quality samples were removed from the analysis, and we excluded peptides that were not detected in more than 80% of the samples (see Methods). The generated dataset hence encompasses protein abundance quantification for 630 proteins among 942 isolates (Table S1 and Table S2). This dataset therefore covers the overall genetic diversity of the species and captures the subpopulations that were defined as part of the 1,011 yeast genomes project, including both domesticated and wild clades (Peter et al., 2018) (Figure S1A). We combined the proteomic dataset with transcriptomic data obtained from the 1,011 strains collection (Caudal et al., 2023), which gave access to the quantified expression of both levels for 629 genes across 889 isolates (Figure 1B, Table S1). To be able to properly compare these two datasets, we normalized them with quantile normalization after imputing the missing values using the KNN method (Table S3, Figure S1B-C).

To characterize the quantified proteins in our study, we first compared the level of transcription of both the identified and unidentified proteins. Low abundance transcripts are less likely to be quantified by proteomics as compared to high abundant transcripts (Figure S1D). Indeed, 489 ). The dotted lines correspond to the median correlation index for the proteomic (yellow) and transcriptomic (blue) data. C. Cellular functions that are preferentially affected by post-transcriptional buffering. Briefly, using either the proteome and the transcriptome abundances -1-we constructed expression-based neighbor joining trees -2-and compared the total sum of the branch lengths. We computed a ratio -3-defined by the proteome total branch lengths divided by the transcriptome total branch lengths. Using all the genes, this ratio was equal to 0.93 (overall, the expression evolution is more constrained at the proteome level). We performed the same procedure using subsets of genes corresponding to 101 biological process annotations. The biological processes

Materials and methods

Cultivation of library for proteomics

The yeast isolate collection was grown on agar containing synthetic complete medium (SC; 6.7 g/L yeast nitrogen base (MP Biomedicals, Cat#114027512-CF), 20 g/L glucose, 2 g/L synthetic complete amino acid mixture (MP Biomedicals, Cat#114400022)). After 48 h, colonies were inoculated in 200 µL SC liquid medium using a Singer Rotor and incubated at 30 °C overnight without shaking. These pre-cultures were then mixed by pipetting up and down, and diluted 20x by transferring 80 µL per culture to deep-well plates pre-filled with 1.55 mL SC liquid medium and one borosilicate glass bead per well. Plates were sealed with a permeable membrane and grown for 8 h at 1000 rpm, 30°C to exponential phase. The optical density at harvest was measured using an Infinite M Nano (Tecan). Per culture, 1.4 mL of cell suspension were harvested by transferring into a new deep-well plate and subsequent centrifugation [START_REF]Consistently with our previous findings on gene expression variation, the genes involved in different metabolism pathways were associated with fast gene expression evolution. As the 8 isolates came from different environments with very different conditions in terms of nutrient and resource availability, this highlights that trophic constraints play an important role in shaping gene expression. Conversely, genes that are involved in multiple protein interactions or that tend to play an essential role in cellular functions are likely to be associated with strong constraints on expression evolution[END_REF]220 x g, 5 min, 4°C). The supernatant was removed by inverting the plates. Cell pellets were immediately cooled on dry ice and stored at -80°C.

Sample preparation

Samples for proteomics were prepared as previously described (Messner et al., 2022a(Messner et al., , 2020;;Muenzner et al., 2022). In brief, samples were processed in 96-well format, with lysis being achieved by beat beating using a Spex Geno/Grinder and 200 µL of lysis buffer (100 mM ammonium bicarbonate, 7 M urea). Samples were reduced and alkylated using DTT (20 µL, 55 mM) and iodoacetamide (20 μL, 120 mM), respectively, diluted with 1 mL 100 mM ammonium bicarbonate, and 500 µL per sample were digested using 2 µg Trypsin/LysC (Promega, Cat#V5072). After 17 h of incubation at 37°C, 25 µL 20% formic acid were added to the samples, and peptides were purified using solid-phase extraction as described previously (Messner et al., 2020). Eluted samples were vacuum-dried and subsequently dissolved in 70 µL 0.1% formic acid. An equivoluminal pool of all samples was generated to be used as technical controls (QCs) during MS measurements. The peptide concentration of this pool was determined using a fluorimetric peptide assay kit (Thermo Scientific, Cat#23290). Peptide concentrations per sample were estimated by multiplying the optical density recorded at harvest with the ratio between pool peptide concentration and the median at-harvest optical density.

GO analysis looking at the functional enrichment present in the 630, the gene list reference was the genes encompassed in the transcriptomic data (Caudal et al., 2023). The others GO analyses used the 628 genes as the reference list. All the others GO analyses were performed using the same procedure, unless specified.

Expression variation exploration

We measured the strength of protein and transcript abundance variation using several methods.

We computed an absolute transformed Log2(fold change) value (|Log2(FC)|) where in each isolate pairwise comparison (ex: strain A vs strain B) and for each gene, we performed: Briefly, the more this value increases, the more different is protein abundance between two isolates for a specific gene. We also computed a pairwise spearman correlation between the isolates using the normalized proteomic and transcriptomic data. We also gathered the Euclidean distances between the expression profiles of each isolate, as well as the gene expression variance per gene.

We explored the post-transcriptional buffering phenomenon using an approach based on the computation of expression trees (Wang et al., 2020). First, on both protein and transcript normalized abundances, we constructed a neighbor-joining tree based on the Euclidean distance between each isolate. We computed the total branch length of these two trees and created a ratio of the proteome tree length on the transcriptome tree length. The ratio was equal to 0.93 which is line with the difference in Euclidean distance between the transcriptome and proteome. We performed 100 bootstrapping tests and used the resulting branch lengths to test the difference between the proteome and the transcriptome tree. We sought to check if some cellular pathways tended to be more affected by the post transcriptional buffering phenomenon. To do so, we gathered a reduced biological process GO annotation by computing the similarity between each GO term using the rrvgo R package and the 'Resnik" method (Sayols, 2022). We discarded terms that are at least 50% overlapping with another term and the terms encompassing no more than 5 genes, which resulted in a list of 101 terms. For each of these terms, we performed the same tree exploration, but this time with the genes encompassed by each term. We obtained therefore 101 tree length ratios. We selected the terms displaying a ratio lower than 0.93 or higher than 1, and for which the total branch length between the proteome and the transcriptome

Supplementary Material

Supplementary tables available at: https://www.dropbox.com/s/ukah0o3q2b8e4pw/Sup_tables.xlsx?dl=0 M3. 8. 25. 26. 3. 10. 5. 13. M2. 2. 6. 7. 12. 23. M1. 11. 24. 18. 21. 4. 9. 19. 17. 22. 14. 15. 20 M3. 8. 25. 26. 3. 10. 5. 13. M2. 2. 6. 7. 12. 23. M1. 11. 24. 18. 21. 4. 9. 19. 17. 22. 14. 15. 20.

Subpopulation

Signature genes & isolate number

CONCLUSION & PERSPECTIVES

In the context of unraveling the links between genetic variation and observable traits within a population or species, it is critical to better characterize intermediate traits (i.e. molecular traits).

In this perspective, we aimed at exploring gene expression through mRNA abundance, protein abundance, as well as translation itself in the budding yeast Saccharomyces cerevisiae. The study of their variation among individuals as well as their correlation and interaction were conducted to better understand how genetic diversity shapes such expression variation at the population level. We have used two complementary approaches to study gene expression: a deep gene coverage approach and a population-level approach. The combination of these two types of approaches is still essential today, as each approach emphasizes different aspects of gene expression variation while suffering from some specific limitations.

Complementary approaches for a better understanding of gene expression variation

Our investigation of gene expression variation using a large gene coverage on a limited number of individuals showed that despite functional similarities between the expression levels (transcriptome, translatome and proteome), there is an important and general trend of variation buffering: the transcriptional variations are buffered at the translatome and proteome levels. We observed that this phenomenon increases as long as the expression process progresses, suggesting higher constraints on the translatome and even more on the proteome. We showed that this buffering, usually referred to as post-transcriptional buffering (Blevins et al., 2019;Gonçalves et al., 2017;McManus et al., 2014;Wang et al., 2020), affects genes unevenly.

Indeed, genes such as essential genes or the ones involved in protein complexes are preferentially buffered, this trend being observed at both the translatome and the proteome level. This suggests that several mechanisms underlie this phenomenon. Some already described mechanisms are obviously good candidates, such as those preventing unassembled protein complex components (Chotewutmontri and Barkan, 2016;Jüschke et al., 2013;Lukoszek et al., 2016;Trösch et al., 2018) or, more generally, protein degradation (Gonçalves et al., 2017). However, this phenomenon remains largely elusive and should be further analyzed and explored to clarify it.

Overall, despite allowing for higher gene coverage and more precision, exploring gene expression in a limited number of strains is unlikely to reveal large-scale effects within a large population and poorly suited for the systematic exploration of the genetic origins of gene expression variation.

Population-level exploration of gene expression variation is a more appropriate strategy in this regard as some specific explorations, such as co-expression network or mRNA-protein variation correlation, can only be considered with large-scale gene expression surveys. We therefore examined mRNA and protein abundance at the population level and found that transcriptome and proteome variation were poorly correlated across individuals, in contrast to previous observations in smaller datasets (Albert et al., 2014;Aydin et al., 2023;Buccitelli and Selbach, 2020;Wang et al., 2019). We observed that genetic regulation of protein abundance is highly distinct from genetic regulation of mRNA abundance at the population level and that this is partly due to variation in protein degradation, which may play a central role in proteomespecific regulation. While the relationship and dependency between transcript and protein levels has been debated (Buccitelli and Selbach, 2020;Liu et al., 2016;Upadhya and Ryan, 2022), it remains poorly understood to date. Our dataset represents the first population-level, multiomics exploration and demonstrates that transcriptomes and proteomes are clearly two distinct layers of regulation, governed by different genetic bases in natural populations, highlighting the importance of integrating these different levels of gene expression to better understand the genotype-phenotype relationship. However, and despite these promising results, systematic studies at the population-scale level are still tedious and costly to perform. Indeed, the implementation of these strategies usually requires the reduction of experimental time and cost per individual. This represents a bottleneck at the proteomic side, for which gene coverage remains relatively low when applied at high-throughput. Consequently, large-scale proteome quantifications still require a sharp trade-off between the number of considered isolates and the proteome coverage.

Towards a larger view of gene expression

Despite representing one of the largest gene expression explorations to date, several crucial points were beyond the scope of our work, either because of technical limitations or the fact that gene expression encompasses tenths of different mechanisms.

A more exhaustive transcriptome and proteome exploration

A larger coverage of the proteome

As mentioned above, large-scale proteome research still suffers from gene coverage limitation. This is related to the methodology used for high-throughput proteomics where several technical

APPENDIX

List of publications

Submitted

Teyssonnière, E., Schichino, Y., Friedrich, A., Iwasaki, S., Schacherer, J., 2023. Translation variation across genetic backgrounds reveals a post-transcriptional buffering signature in yeast.

(In revision -Nucleic Acid Research)

Teyssonnière, E., Trébulle, P., Muenzner, J., Loegler, V., Ludwig, D., Amari, F., Mülleder, M., Friedrich, A., Hou, J., Ralser, M., Schacherer, J. Species-wide quantitative transcriptomes and proteomes reveal distinct genetic control of gene expression variation in yeast. (Submitted to Cell, September 2023)

Caudal, E., Loegler, V., Dutreux, F., Vakirlis, N., Teyssonnière, E., Caradec, C., Friedrich, A., Hou, J., Schacherer, J., 2023. Pan-transcriptome reveals a large accessory genome contribution to gene expression variation in yeast. https://doi.org/10.1101/2023.05.17.541122 (In revision -Nature Genetics)

In preparation

Teyssonnière, E., Dubreuil, B., Levy, E.D., Schacherer, J., Metabolism adaptation is a main driver of protein abundance evolution in the yeast Saccharomyces cerevisiae.

Oral Communications

Oral

International colloquium: iGénolevure "Yeast Biodiversity": Teyssonniere, E., Schichino, Y.,

A A A A A A A A

Abstract

An astonishing phenotypic diversity can be observed in natural populations. One of the major goals of modern biology is to unravel the genetic origins of this phenotypic landscape. Gene expression is known to be a main determinant of the relationship between genotypes and phenotypes. In recent decades, several analytical and technical advances have made it possible to study gene expression at every step of the expression process (e.g., transcriptome and proteome) and at very large scales. However, a complete exploration of gene expression across the entire process and at the population scale is still lacking. The goal of this dissertation is to get a more comprehensive view of how each layer of gene expression varies, influences each other, and is related to the natural genetic diversity observed within species. To this end, we analysed the transcriptomes and proteomes of a large natural population of S. cerevisiae (bringing together more than 1,000 individuals) and found unsuspected differences between mRNA and protein abundance regulation. Simultaneously, we studied the gene expression process at three different molecular levels (transcriptome, translatome and proteome) and found that important buffering mechanisms underlie the expression variation between individuals.

Elie Teyssonnière

Infraspecific variation and genome evolution