Positive selection of senescence through increased evolvability: ageing is not a by-product of evolution

For the past century, scientists have debated whether or not ageing is directly selected by evolution. Because ageing occurs, by definition, late in life -that is, after the organism's development is complete -many people believe that it cannot be actively selected for as a process. Furthermore, because it reduces an individual's fitness, it is thought unlikely to be selected for. In agreement with this viewpoint, numerous theories have been proposed in the last 75 years to explain the observation of its widespread presence in the realm of life. The bd model of ageing is based on a simple life-history trait model that we recently introduced. Our model suggests that senescence can be positively selected through evolution due to the increased evolvability it confers on organisms, not through a specific mechanism but through a 'function ageing', limiting organismal maintenance and reproduction. It confirms the substrate for mutation accumulation and antagonistic pleiotropy theories while providing an elegant explanation for the apparent tradeoff between longevity and fertility that led to the disposable soma theory without requiring an energy tradeoff. Furthermore, it predicts that the Lansing effect will be present in organisms that exhibit rapid post-reproductive senescence. This formal and numerical modeling of the evolution of ageing also provides new hints for testing the validity of existing theories.

Introduction

Ageing is defined as the effect of elapsed time on an organism; it is also known as'senescence'. Depending on the organism studied, it manifests itself in a wide range of age-related patterns of mortality, ranging from negligible senescence to post-reproductive death via progressive or brutal age-dependent mortality increase [START_REF] Jones | Diversity of ageing across the tree of life[END_REF]. Numerous explanatory theories of ageing have been proposed in the past century (reviewed in [START_REF] Kirkwood | The evolution of ageing and longevity[END_REF]. Whether they are focusing on the molecular aspects that lead to the observed ageing phenomena or on the evolution of the process, one fundamental question about the selectability of ageing through evolution remains. Soon after Charles Darwin published his theory of evolution, August Weismann proposed evolutionary arguments to explain ageing [START_REF] Weismann | Pleiotropy, Natural Selection, and the Evolution of Senescence[END_REF]. "His initial idea was that there exists a specific death-mechanism designed by natural selection to eliminate the old, and therefore worn-out, members of a population" [START_REF] Gavrilov | Evolutionary theories of aging and longevity[END_REF]. Since then, and probably partially due to Weismann's later changes of mind on his early theory, it is mostly accepted that "ageing is not adaptive since it reduces reproductive potential" [START_REF] Kirkwood | The evolution of ageing and longevity[END_REF] and hence, fitness. Two other arguments are usually used against the selection of ageing. First, that this would have required a group selection effect stronger than selection at the individual level, which is very seldom the case [START_REF] Smith | What Determines the Rate of Evolution?[END_REF]. Second, the mortality rate is so high in early life that little senility is observed in the wild (reviewed in [START_REF] Nussey | Senescence in natural populations of animals: widespread evidence and its implications for bio-gerontology[END_REF]). There would thus be little opportunity for a removal mechanism to evolve (reviewed in [START_REF] Johnson | Revamping the evolutionary theories of aging[END_REF]). Since the 1950s onwards, evolutionary theories are mainly building upon the idea that ageing is a byproduct of natural selection [START_REF] Fabian | The Evolution of Aging[END_REF]. A first example is Peter Medawar's theory of mutation accumulation for which ageing is caused by the progressive accumulation of deleterious mutations with effects that manifest only late in life [START_REF] Medawar | A birth-death model of ageing: from individual-based dynamics to evolutive differential inclusions[END_REF]. Williams' antagonistic pleiotropy theory goes a step further than Medawar's by assuming the existence of antagonistic genes and mutations: They would not be selected against for negative effects manifesting later in life if they were beneficial at a young age. Finally, the "disposable soma" theory of ageing proposed by Thomas Kirkwood, is based on the idea that individuals have a limited amount of energy to be split between reproductive functions and (non-reproductive) maintenance of the organism, the "soma". According to Kirkwood's theory, increasing an organism's longevity would thus be associated with a decrease in growth and reproduction rates, delaying death [START_REF] Kirkwood | Evolution of ageing[END_REF]. Later, evolutionary conserved genes involved in both the regulation of longevity and organismal growth were discovered in the model organism C. elegans [START_REF] Kenyon | A C. elegans mutant that lives twice as long as wild type[END_REF], later shown to be conserved in flies [START_REF] Clancy | Extension of life-span by loss of CHICO, a Drosophila insulin receptor substrate protein[END_REF] , mice [START_REF] Bluher | Extended longevity in mice lacking the insulin receptor in adipose tissue[END_REF] and humans [START_REF] Van Heemst | Reduced insulin/IGF-1 signalling and human longevity[END_REF]. Thus, genetic modulators for longevity exist and manifest themselves through evolutionarily conserved physiological mechanisms. Nevertheless, although "the evolutionary theories of aging are closely related to the genetics of aging because biological evolution is possible only for heritable manifestations of aging" [START_REF] Gavrilov | Evolutionary theories of aging and longevity[END_REF] and we now have countless examples of evolutionarily conserved genes playing a role in ageing across species [START_REF] Partridge | Mechanisms of ageing: public or private?[END_REF], the subject is still a matter of vivid debate [START_REF] Kowald | Can aging be programmed? A critical literature review[END_REF]. In their 2016 review, Axel Kowald and Thomas Kirkwood state that the "idea that aging is a programmed trait that is beneficial for the species [...] is now generally accepted to be wrong". If ageing cannot be positively selected for through evolution, can it be, at least partially, programmed? This question raised no less vivid debates than the previous one in the past century [START_REF] Austad | Rebuttal to Bredesen: 'The non-existent aging program: how does it work?[END_REF]Bredesen, 2004aBredesen, , 2004b;;[START_REF] Gavrilov | Evolutionary theories of aging and longevity[END_REF][START_REF] Kirkwood | On the programmed/non-programmed nature of ageing within the life history[END_REF][START_REF] Longo | Programmed and altruistic ageing[END_REF][START_REF] Skulachev | Aging as a particular case of phenoptosis, the programmed death of an organism (A response to Kirkwood and Melov "On the programmed/non-programmed nature of ageing within the life history[END_REF]. Two good examples of these active debates can be found in [START_REF] Blagosklonny | Aging is not programmed: Genetic pseudo-program is a shadow of developmental growth[END_REF] and [START_REF] Longo | Programmed and altruistic ageing[END_REF], these two works showing dramatically opposed views. The first, affirming that ageing cannot be programmed, the second that it can and will occur as it brings a "kind of population-level selection" that can be explained by kin selection. One of the final arguments given by Kowald and Kirkwood is that if ageing were to be programmed, "it would be possible experimentally to identify the responsible genes and inhibit or block their action". Programmed or not, ageing in unicellular organisms is associated with mechanisms that discriminate new components from older ones as individuals replicate [START_REF] Henderson | Mother-daughter asymmetry of pH underlies aging and rejuvenation in yeast[END_REF]Lai et al., 2002, p. 2;[START_REF] Nyström | A Bacterial Kind of Aging[END_REF][START_REF] Sinclair | Extrachromosomal rDNA Circles-A Cause of Aging in Yeast[END_REF][START_REF] Steiner | Senescence in Bacteria and Its Underlying Mechanisms[END_REF]. In multicellular organisms, the Lansing effect is a good candidate for such a mechanism. It is the effect through which the "progeny of old parents do not live as long as those of young parents" in rotifers [START_REF] Lansing | A Nongenic Factor in the Longevity of Rotifers[END_REF][START_REF] Lansing | A Transmissible, Cumulative, and Reversible Factor in Aging[END_REF]. More recently, it has been shown that older drosophila females and to some extent males tend to produce shorter lived offspring [START_REF] Priest | The role of parental age effects on the evolution of aging[END_REF], zebra finch males give birth to offspring with shorter telomere lengths and reduced lifespans (Noguera et al., 2018) and finally in humans, "Older father's children have lower evolutionary fitness across four centuries and in four populations" [START_REF] Arslan | Older fathers' children have lower evolutionary fitness across four centuries and in four populations[END_REF] [START_REF] Arslan | Older fathers' children have lower evolutionary fitness across four centuries and in four populations[END_REF]. Despite the absence of consensus on the underlying mechanisms [START_REF] Monaghan | Intergenerational Transfer of Ageing: Parental Age and Offspring Lifespan[END_REF], the near-ubiquity of the Lansing effect is important for our understanding of the selective forces shaping the evolution of life histories. In a recent article [START_REF] Medawar | A birth-death model of ageing: from individual-based dynamics to evolutive differential inclusions[END_REF], we introduced an asexual and haploid age-structured population model implementing a strong Lansing effect, showing that this strong transgenerational effect of ageing could be maintained. Here, by extending this model to any system able to reproduce and maintain its homeostasis, we show that 1) ageing -i.e. non-infinite reproduction and homeostasis maintenance -is an adaptive force of evolution, 2) such a system evolves towards a configuration where fertility exceeds homeostasis capabilities, 3) propitious for pro-senescence mechanisms to appear -i.e. Lansing effect -within a few dozen generations. 4) Individuals carrying such a transgenerational effect of ageing have a non-null probability to be selected for when competing with individuals deprived of it. In addition, this model 5) allows to explain longevity-fertility tradeoffs mathematically without the need for energy investment strategies, while 6) suggesting that the selection of ageing relies on the maximization of a meta-characteristic named evolvability, in the form of fitness gradient. Finally, the landscape of this fitness gradient shows how mutation accumulation can accompany the evolution of ageing.

Results

A generalized bd model

The model and its population dynamics is similar to those described in [START_REF] Medawar | A birth-death model of ageing: from individual-based dynamics to evolutive differential inclusions[END_REF]. Briefly, the model describes an asexual and haploid population structured by a life-history trait that is defined by a pair of parametersor genes -(x b , x d ) where x b is defining the duration of fertility and x d the age at which mortality becomes non-null.

Here we generalized the model to any intensities of birth and death denoted (i b , i d ) as well as to populations without Lansing effect (Figure 1, see also Annex 1). The selective pressure is enforced by a logistic competition c mimicking a maximum carrying capacity of the environment, thus no explicit adaptive value is given to any particular trait and, for each reproduction event, a mutation (h) of probability p can affect both the genes x b and x d , independently, following a Gaussian distribution centered on the parental trait. In Figure 1, the different cases are explored, depending on the respective values of x b and x d . Individuals in the Figure 1b-c configuration (for x b ≤ x d ) will always give progeny with a genotype (x b , x d ) ∓ (h b or h d ). The case of individuals carrying a genotype with x d < x b (Figure 1a) is subtler, depending on the parental age a and whether the parent carries the possibility for a Lansing effect or not (Figure 1d-f). If a < x d or if the parent does not carry a Lansing effect, the genotype of the progeny will be as previously described. But if a > x d , and if the parent carries the Lansing effect, then the progeny inherits a dramatically reduced x d (here x d is set to 0), mimicking a strong Lansing effect. In the case of physiologically young parents (a < x d ), the progeny's genotype is that of its parent ∓ a Gaussian kernel of mutation centered on the parental gene. In the case of the reproduction event occurring after x d , for configuration (a) above, two cases are observed, (e) if the organism carries a Lansing effect ability, the x d of its progeny will be strongly decreased. (f) In the absence of the Lansing effect, the default rule applies.

In our previous work [START_REF] Medawar | A birth-death model of ageing: from individual-based dynamics to evolutive differential inclusions[END_REF], we formally and numerically showed the long-time evolution of the model to converge towards (x b -x d ) = 0 in the case of individuals carrying a Lansing effect. To test whether this convergence of (x b -x d ) still occurred without this strong transgenerational effect of ageing, we implemented a new version of the model devoid of Lansing effect and simulated its evolution for a viable -i.e. allowing the production of at least one progeny -trait (x b = 1.2, x d = 1.6). Surprisingly, we still observe a convergence of (x b -x d ). The dynamics of the trait (x b , x d ) is described by the canonical equation of adaptive dynamics depending on the Malthusian parameter and its gradient (Annexe 1). The latter can be interpreted as the age-specific strength of selection in the sense of Hamilton. We observe that the evolution speed of x b and x d decreases with time as the previous -less general -form of the model did [START_REF] Medawar | A birth-death model of ageing: from individual-based dynamics to evolutive differential inclusions[END_REF]. This allows us to recover the classical age-related decrease in the strength of selection [START_REF] Hamilton | The moulding of senescence by natural selection[END_REF][START_REF] Haldane | New Paths In Genetics[END_REF][START_REF] Medawar | A birth-death model of ageing: from individual-based dynamics to evolutive differential inclusions[END_REF]. Simulations of the generalized bd model presented here show that the x b -x d distance, i.e. the time separating the end of fertility from the increasing risk of death, converges -for any initial trait -towards a positive constant. Thus, it seems that the long term evolution of such a system is a configuration similar to Figure 2a (x d < x b ). The formal analysis of the generalized bd model confirmed that the long-time limit of the traits (x b -x d ) is the positive constant defined by the formula in Figure 2b (mathematical analysis Annexe 4), reached after a few dozen simulated generations, estimated by (time ≅ 75) / (average offspring) (Figure 2c). Although we formally demonstrate the long-time limit for any i b and i d , all our simulations are run using i b = i d = 1 in order to limit the number of conditions to assess and report. Surprisingly, the limit value of the trait is not affected by x b or x d values -the fertility and organismal maintenance durations per se -but only by their respective intensities i b and i d . These intensities can be interpreted as the instant mortality risk i d and the probability to give a progeny i b . Interestingly, the long-time limit values for any i b and i d shows a significantly stronger sensitivity to the increasing mortality risk i d than to reproduction by almost two orders of magnitude (Figure 2d). In addition, for extremely low values of i b and i d -i.e. below 0.01 -the apparent time correlation of the fertility and organismal integrity maintenance period is almost nonexistent since (x b -x d ) is large. A biological manifestation of this would be a loss of organismal maintenance occurring long before the exhaustion of reproductive capacity. Such an organism would be thus characterized as having no significant fertility decrease during the ageing process. On the other end, for individuals showing either a high instant mortality risk or a high probability to give a progeny, the (x b -x d ) trait is close to 0, meaning that fertility and organismal integrity maintenance are visibly correlated. Note that this mathematical study concerns individuals for which the mean number of descendants per individual is large enough, allowing us to define a viability set of traits (x b , x d ) (see Annexe 2.3). Because of these mathematical properties, a tradeoff emerges between i b , i d , x b and x d . Let's consider an organism -for both the Lansing and non-Lansing cases -with a low reproductive intensity i b = 0.01 and i d = 1. For it to be viable, the product i b * x b has to be strictly superior to 1, hence here x b ⩾ 100 (see Annexe 2.3). In this example, the long-time limit of the trait (x b -x d ) is equal to ln(2), thus x b ≅ x d . With the same reasoning, considering an organism significantly more fertile (with i b = 1, i d = 1), the long-time evolution lower limit of (x b -x d ) is 1/√3. This model thus allows an elegant explanation for the apparent negative correlation previously described between longevity and fertility without the need of implementing energy trade-offs or relative efficiency of energy allocation between maintenance and reproduction (see Annexe 2.3 -examples).

On the selection of Lansing effect

In our model, whatever the initial trait (x b , x d ) in the viability set, evolution leads to a configuration of the trait such that the risk of mortality starts to increase before the fertility period is exhausted. Similar to biochemical reactions involved in a given pathway being evolutionarily optimized through tunneled reactions and gated electron transfers, we hypothesize here that such a configuration, caused by simple mathematical constraints, creates the conditions for the apparition, selection and maintenance of a molecular mechanism coupling x b and x d . Such a coupling mechanism could thus be the Lansing effect, the only described age-related decline in progeny's quality that seems to affect numerous iteroparous species [START_REF] Lansing | A Transmissible, Cumulative, and Reversible Factor in Aging[END_REF][START_REF] Monaghan | Intergenerational Transfer of Ageing: Parental Age and Offspring Lifespan[END_REF]. More importantly, the vector of the Lansing effect could be as simple as a transgenerational transfer of maternal components known in yeast as asymmetrical inheritance, a mechanism of age-reset in Saccharomyces cerevisiae [START_REF] Hill | Asymmetric Inheritance of Aggregated Proteins and Age Reset in Yeast Are Regulated by Vac17-Dependent Vacuolar Functions[END_REF]. We calculated the chances of an organism with such a non-genetic pro-senescence mechanism surviving in competition with a population lacking such a mechanism. To do so, we looked at a population divided into two sub-populations: one made up of people who experienced the Lansing effect and the other made up of people who did not. As before, we assume that each individual is subject to the same competitive pressure. The Darwinian fitness of the two initial subpopulations is approximated by their Malthusian parameter (cf. Annexe 2, supplementary figure 2). Their traits are thus (1.5; 1.3) Lansing and (1.5; 0.83) non-Lansing . In order to simplify the analysis, both the birth and death intensities are as follows: i b =i d (the model is nevertheless generalized to any (i b ; i d ), see Annexe 5.1) and evolution was simulated for discrete pairs of mutation rate (p) and competition (c) parameters. Three indexes were calculated for each set of simulation: Table 1 a) the ratio of Lansing and non-Lansing populations that collapsed (a "-" indicates that all survived), Table 1 b) the ratio of total number of progenies produced during the simulation by each population and Table 1 c) the relative proportion of the Lansing population at the end of the simulation. Our 1200 simulations each of 2.10 5 birth-death events summarized in Table 1 show that the Lansing populations survive at least as well as non-Lansing ones (Table 1a) and show a significantly increased success for a moderate competition (c = 9.10 -4 ) and low (in our simulations) mutation rate (p = 0.1). With these parameters, Lansing populations show almost half the risk of disappearance of non-Lansing ones (Table 1a), producing nearly thrice as many descendants as non-Lansing populations (Table 1b) for up to a 20% faster growing population (Table 1c). Hence, although the Lansing effect leads to the production of a significant proportion of progeny with an extremely low fitness (x d = 0), it decreases the risk of population collapse for organisms carrying it and seems to allow a slightly better growth of the population, whatever the magnitude of the Lansing effect (Sup. Figure 1).

Table 1: Populations with Lansing effect are favorably selected under logistic competition when the mutation rate is non-zero. p is the mutation rate and c the intensity of the logistic competition. For each couple (p, c), 100 independent simulations were run with 500 individuals per population at t 0 of which traits are (1.5; 1.3) Lansing and (1.5; 0.83) non-Lansing so that their respective Malthusian parameters are equal. Each simulation corresponds to 2.10 5 events of birth or death. Table (a) shows the ratio of Lansing and non-Lansing populations (out of 100 simulations in each case) that did collapse by the end of the simulation. For the lowest competition, none of the populations collapsed within the timeframe of simulations (-). For an intermediate value of competition, approximately half less Lansing populations disappear relative to non-Lansing ones. Table (b) shows the ratio of the number of individuals generated between Lansing and non-Lansing populations. On average, Lansing populations generate approximately twice as many individuals as non-Lansing ones. (c) On average, Lansing populations grow 20% more than the non-Lansing. Values highlighted in green are discussed further below.

The Lansing effect increases fitness in ageing populations by increasing their evolvability.

We focused on the Malthusian parameter to understand the evolutionary success of a characteristic that appears to reduce an organism's fitness -the genotypic rate of increase [START_REF] Hairston | Natural Selection and the Parameters of Population Growth[END_REF] -in each population through time.

Here we present the results for an intermediate set of c and p -highlighted in green Table 1 -that we identified as associated with the highest success rate of Lansing bearing populations. First, we observe that, on average, Lansing populations (blue) grow while non-Lansing ones (red) have a decreasing size (Figure 3a -blue and red curves represent deciles 1, 5 and 9). Nonetheless, in simulations where both populations coexist continuously, the Lansing population's higher fitness is marginal, with populations growing 20% faster than the non-Lansing population (Figure 3b). This higher success rate appears to be due to the Lansing population's faster and broader exploration of the Malthusian parameter space (Figure 3c). This maximization of the Malthusian parameter is not associated with any significant difference in either population's lifespan (time of death -time of birth) distributions (Figure 3d). Indeed, despite having the same mutation rate p and being subject to the same competition c, the distribution of non-Lansing progenies is essentially of the parental trait in the first 5 generations, whereas Lansing progenies (unaffected by the Lansing effect -we excluded progeny with x d = 0 for the comparison) explore a broader part of the trait space (Figure 3e). Although their low fitness progeny (x d = 0) represent up to 10% of the population for a significant amount of time, the Lansing-bearing population has significantly higher success (Figure 3f). As a result, Lansing populations reach the equilibrium trait faster than non-Lansing populations (Figure 3g). Thus, the relatively higher success of Lansing bearing populations appears to be associated with increased genotypic diversity, which leads to a broader range of fitness on which natural selection can act, namely evolvability.

The relative success of Lansing-bearing populations with randomly distributed traits

We have proposed a simple model that demonstrates the mathematical basis driving the evolutionary pressure that connects organismal maintenance and reproductive mechanisms. Nonetheless, our numerical investigation of our model's behavior has so far been limited to initial conditions in which the competing populations had equal Malthusian parameters. The low number of generations involved suggest that the conditions for the development, selection and maintenance of mechanisms of ageing [START_REF] Lemoine | The Evolution of the Hallmarks of Aging[END_REF] would have occurred early on, in a population of mixed individuals. As such, we decided to test the evolution of the trait (x b -x d ) in Lansing and non-Lansing bearing individuals of traits uniformly distributed on [-10; +10] (Figure 4 -left panel). We chose to plot one (Figure 4 -central panel) of the hundred simulations we made, that is representative of the general results. Simulations show, over 110 million individuals, an early counterselection of extreme trait values, typically (x b -x d ) > 4. Interestingly, the whole space of (x b -x d ) trait is not explored evenly and the positive part of it represents approximately ⅔ of the individuals although the branched evolution process can lead a line on both the positive ('Too young to die' -Figure 1a) and negative ('Menopause' -Figure 1c) sides of the trait space. Both the Lansing and non-Lansing bearing populations manage to co-exist until the end of the simulation, each reaching a distribution centered on their respective theoretical solutions (Figure 4 -right panel), 0 for the Lansing (Méléard et al., 2019) and ln(3)/2 for the non-Lansing. When the initial condition does not restrict competition to populations with identical Malthusian parameters, the Lansing bearing population is significantly less successful than the non-Lansing one, accounting for only one-third of the final population size. As such, the evolution of a mixed population of individuals with a trait (x b -x d ) initially uniformly distributed on [-10; +10], bearing or not a strong inter-generational effect, will lead to a mixed solution of individuals carrying a trait that converged towards the theoretical solution such as x d ≲ x b thus allowing the maximization of fertility without cluttering the environment with individuals not producing descendants, very similar to Weismann's first intuition [START_REF] Weismann | Pleiotropy, Natural Selection, and the Evolution of Senescence[END_REF]. Nonetheless, this interpretation appears to be finalist and does not account for the relative success of each of our population types. We then investigate how ageing affects the evolvability of individuals who bear it, leading us to propose another interpretation of the fitness gain it brings. The fitness gradient, a mediator of evolvability

We focused our attention on the differential landscape of the Malthusian parameters as a function of the trait (x b , x d ) for both Lansing and non-Lansing populations to understand the origin of this relative success of individuals carrying the ability to transmit an ageing information to the next generation. We built this landscape numerically using the Newton method implemented in Annex 2. First of all, it is interesting to notice that we have derived from the equations that the maximum rate of increase for Malthusian parameters is 1/i d with a maximum fitness value capped by i b (Annex 2). Consistent with our previous characterization of the Trait Substitution Sequence in populations with Lansing effect (Méléard et al., 2019), Lansing individuals have a symmetrical fitness landscape (Figure 5, blue lines) centered on the diagonal x = x d (Figure 5, green diagonal). Along the latter, we can directly observe what is responsible for the "selection shadow". As (x b , x d ) increases, a mutation of same magnitude has smaller and smaller effects on the fitness, thus allowing for the accumulation of mutations (Figure 5, blue arrows). Non-Lansing individuals' cases are asymmetric, with symmetry breaking on the x b = x d diagonal. For x d > x b (Figure 5, upper diagonal), fitness isoclines fully overlap thus showing an equal response of both Lansing and non-Lansing fitness to mutations. In addition, as expected, the fitness of Lansing individuals is equal to that of non-Lansing ones for a given trait. On the lower part of the graph, corresponding to x d < x b , non-Lansing fitness isoclines separate from that of Lansing individuals, making the fitness of non-Lansing individuals higher to that of Lansing ones for a given trait. Nevertheless, the fitness gradient is significantly stronger for Lansing individuals as represented on Figure 5 by the yellow arrow and associated yellow area. For an individual of trait (x b = 2.45; x d = 1.05) a mutation making a non-Lansing individual gain 0.1 in fitness (isocline 0.7 to isocline 0.8) will make a Lansing individual increase its own by 0.42 (isocline 0.1 to above isocline 0.5). With a 4-fold difference, the Lansing population produces 4 times as many individuals as the non-Lansing ones for a given mutation probability. But this reasoning can be extended to any trait (x b , x d ) with or without Lansing effect. Organisms ageing rapidly -i.e. with low x b and x d -will see their fitness significantly more affected by a given mutation h than individuals with slower ageing affected by the same mutation. As such, ageing favors the emergence of genetic variants, it increases an organism's evolvability.

Figure 5: The Lansing effect is associated with an increased fitness gradient. We were able to derive Lansing and non-Lansing Malthusian parameters from the model's equations (see Annexe 1-2.3 and 1-5) and plot them as a function of the trait (x b , x d ). The diagonal x b = x d is drawn in light green. The corresponding isoclines are overlapping above the diagonal but significantly differ below, with non-Lansing fitness (red lines) being higher than that of Lansing's (light blue lines). In addition, the distance between two consecutive isoclines is significantly more important in the lower part of the graph for non-Lansing than Lansing bearing populations. As such, a mutation leading a non-Lansing individual's fitness going from 0.7 to 0.8 (yellow arrow) corresponds to a Lansing individual's fitness going from 0.1 to 0.52. Finally, Hamilton's decreasing force of selection with age can be observed along the diagonal with a growing distance between two consecutive fitness isoclines as x b and x d continue increasing.

Discussion

Ageing has a variety of effects on a wide range of organisms [START_REF] Jones | Diversity of ageing across the tree of life[END_REF]. However, its characteristics exhibit evolutionarily conserved so-called ageing hallmarks [START_REF] Lopez-Otin | The hallmarks of aging[END_REF], the most evolutionarily ancient of which is the loss of protein stability [START_REF] Lemoine | The Evolution of the Hallmarks of Aging[END_REF]. For nearly a century and a half, researchers interested in understanding the evolutive role of ageing have been arguing about why ageing exists at all. Although early theories [START_REF] Weismann | Pleiotropy, Natural Selection, and the Evolution of Senescence[END_REF] proposed a population-based adaptive role for it, the debate appears to have been settled in the twentieth century after seminal works demonstrating the time-dependent decline of selective pressure [START_REF] Hamilton | The moulding of senescence by natural selection[END_REF]Medawar, allowed for the existence of theories [START_REF] Kirkwood | Evolution of ageing[END_REF][START_REF] Medawar | A birth-death model of ageing: from individual-based dynamics to evolutive differential inclusions[END_REF][START_REF] Weismann | Pleiotropy, Natural Selection, and the Evolution of Senescence[END_REF] proposing that ageing is nothing more than a by-product of evolution. The model we presented here allows us to propose an alternative theory where ageing necessarily emerges for any system showing the two minimal properties of life, namely a) reproduction with variation (x b ) and b) organismal maintenance (x d ) [START_REF] Trifonov | Vocabulary of Definitions of Life Suggests a Definition[END_REF]. We formally show that a haploid and asexual organism with these two properties will rapidly evolve, within a few dozen generations, towards a solution such that (x b -x d ) is strictly positive. More importantly, the time separating both parameters is independent of their absolute values and is solely determined by their respective rates, respectively i b for x b and i d for x d . This property enables us to explain the observations that fueled the disposable soma theory of ageing, which is based on an apparent trade-off between an organism's fertility and lifespan. Indeed, in our model, the lower limit condition for an individual's viability is x * i b > 1. As such, an organism with a low fertility (i b << 1) will require a long fertility time (x b >> 1) to be viable. The formally demonstrated properties of our model imply that the duration for homeostasis maintenance (x d ) will rapidly evolve towards (x b -x d ) equal to a positive constant, implying that x d >> 1. A highly fertile organism, on the other hand, will evolve to its minimum viable condition, requiring only a small x d . The apparent trade-off between fertility and longevity is thus solely a consequence of x b i b > 1 and lim +∞ (x b -x d ) t , without the need to implement any constraint on resource allocations. By forcing xb and xd to converge, evolution creates the conditions for the emergence of a stage of life in which an individual's fertility declines while their risk of death increases. This time-coupling of the two properties would thus facilitate the selection of any molecular mechanism functionally coupling the two properties [START_REF] Echave | Evolutionary coupling range varies widely among enzymes depending on selection pressure[END_REF], and, contrary to what was proposed in [START_REF] Stearns | Trade-offs in life-history evolution[END_REF], we observe that two genes with no common genetic basis can be co-selected despite being unlinked by any direct tradeoff. By testing such a coupling mechanism, manifested as a strong Lansing effect, we demonstrated that not only is this seemingly fitness-decreasing mechanism capable of increasing fitness, but it also significantly increases the probability of a population bearing it to thrive when competing with a population of equal Malthusian parameter where the effect is absent. We found that a slight increase in fitness caused by the Lansing effect manifests itself numerically as an increase in genetic variability within the population. As a result, we propose that active mechanisms of ageing are chosen during evolution because of their ability to increase an organism's evolvability. The term refers to "an abstract, robust, dispositional property of populations, which captures the joint causal influence of their internal features upon the outcomes of evolution" [START_REF] Brown | What Evolvability Really Is[END_REF]. In other terms, it is "the capacity to generate heritable selectable phenotypic variation" [START_REF] Kirschner | Evolvability[END_REF]. It is an interesting concept as it allows for a character that has no direct effect on fitness -or even a negative one [START_REF] Smith | What use is sex?[END_REF] -to be under strong selection simply for its ability to bring the genetic-phenotypic variation that is the support of evolution. Furthermore, in a constantly changing environment, such a mechanism, triggered when age > x d , would be extremely useful. Indeed, when environmental conditions become less favorable, x d may be affected and individuals pushed to enter the [x d ; x b ] space earlier, increasing population evolvability. Our model predicts a high evolutionary conservation of the "function ageing" -ensuring that an individual has a limited time to propagate its genes -but not necessarily a conservation of the underlying mechanisms, by constraining the long-term evolution of the trait (x b -x d ) without any a priori on the underlying mechanism. This would lead to a layering of mechanisms through evolution as described recently in [START_REF] Lemoine | The Evolution of the Hallmarks of Aging[END_REF]. Being able to derive Hamilton's result [START_REF] Hamilton | The moulding of senescence by natural selection[END_REF] from the mathematical analysis of our model [START_REF] Medawar | A birth-death model of ageing: from individual-based dynamics to evolutive differential inclusions[END_REF] confirms such evolving organisms to be sensitive to the accumulation of mutations [START_REF] Medawar | A birth-death model of ageing: from individual-based dynamics to evolutive differential inclusions[END_REF] while allowing antagonistic pleiotropy to occur [START_REF] Weismann | Pleiotropy, Natural Selection, and the Evolution of Senescence[END_REF]. It is to note that this is, to our the second model showing that the time-dependent decreasing pressure of natural selection derives from ageing not the other way around [START_REF] Giaimo | The selection force weakens with age because ageing evolves and not vice versa[END_REF]. Although this simple model helps us to see ageing as an evolutionarily adaptive force for various ageing types, it is still a toy model. We are now developing more complex versions of it, notably to assess the interactions existing between i b , i d , x b and x d , but also to extend it to joint evolution of maturation, sex, ploidy and varying environmental conditions and ageing.

Materials and Methods

See Annex 2 for code, packages and the software used.

Annexe 1: mathematical proofs

Figure 1 :

 1 Figure 1: Three typical configurations of the model with i b > i d and their effect on progeny's genotypes as a function of parental age. (upper panel) Each haploid individual is defined by a parameter x b defining its fertility period of intensity i b and a parameter x d defining the time during which it will maintain itself, with an intensity i d . These parameters can be positive or null. (a) 'Too young to die' : it corresponds to configurations satisfying x d < x b . (b) 'Now useless': it corresponds to configurations where x b = x d . (c) 'Menopause': it corresponds to configurations where x d > x b . (lower panel) Each individual may randomly produce a progeny during its fertility period [0; x b ]. (d)In the case of physiologically young parents (a < x d ), the progeny's genotype is that of its parent ∓ a Gaussian kernel of mutation centered on the parental gene. In the case of the reproduction event occurring after x d , for configuration (a) above, two cases are observed, (e) if the organism carries a Lansing effect ability, the x d of its progeny will be strongly decreased. (f) In the absence of the Lansing effect, the default rule applies.

Figure 2 :

 2 Figure 2: The bd model shows a convergence of x b -x d towards a positive value. Dynamics of the individual-based model shows a convergence of x b -x d towards a positive constant value in the absence of the Lansing effect. (a) The generalized b-d model shows a convergence of (x b -x d ) for any i b and i d towards a positive value given by (b) (cf. Annexe 4.3, figure 2). (c) Simulation of 1000 individuals with initial trait (x b = 1.2, x d = 1.6) of intensities i b = i d = 1, a competition c = 0.0009 and a mutation kernel (p = 0.1, σ = 0.05) show that the two parameters co-evolve and maintain x b -x d ⪆ 0 . (d) Landscape of solutions (x b -x d ) as a function of i b and i d (colors separate ranges of 50 units on the z-axis).

Figure 3 :

 3 Figure3: The Lansing effect maximizes populational survival by increasing its evolvability. 100 independent simulations were run with a competition intensity of 9.10 -4 and a mutation rate p = 0.1 on a mixed population made of 500 non-Lansing individuals and 500 individuals subjected to such effect. At t 0 , the population size exceeds the maximum load of the medium thus leading to a population decline at start. At t 0 , all individuals are of age 0. Here, we plotted a subset of the 100.10 6 plus individuals generated during the simulations. Each individual is represented by a segment between its time of birth and its time of death. In each graph, blue and red curves represent deciles 1, 5 and 9 of the distribution at any time for each population type. (a) The higher success rate of Lansing bearing populations does not seem to be associated with a significantly faster population growth but with a lower risk of collapse. (b) For cohabitating populations, the Lansing bearing population (blue) is overgrowing by only 10% the non-Lansing one (red). (c) This higher success rate is associated with a faster and broader exploration of the Malthusian parameter -surrogate for fitness -space in Lansing bearing populations (d) that is not associated with significant changes in the lifespan distribution (e) but a faster increase in genotypic variability within the [0; 10] time interval. (f) This occurs although progeny from physiologically old parents can represent up to 10% of the Lansing bearing population and leads to it reaching the theoretical optimum within the timeframe of simulation (g) with the exception of Lansing progenies. (e-g) horizontal lines represent the theoretical limits for (x b -x d ) in Lansing (blue) and non-Lansing (red) populations.

Figure 4 :

 4 Figure 4: Mixed populations lead to (x b -x d ) theoretical limit in a limited time and cohabitation of Lansing and non-Lansing populations. Starting with a homogenous population of 5000 Lansing bearing and 5000 non-Lansing individuals with traits uniformly distributed from -10 to +10 (left panel), we ran 100 independent simulations on time in [0; 1000]. (center panel) Plotting the trait (x b -x d ) as a function of time for one simulation shows a rapid elimination of extreme traits and branching evolution. (right panel) The final distribution of traits in each population type is centered on the theoretical convergence limit for each. N total ≅ 110 million individuals, c = 9.10 -4 , p = 0.1

Acknowledgments

We would like to thank Dr. Sarah Kaakai for her help in transposing our initial simulation Python codes into the IBMPopSim framework, Dr. Allon Weiner for a few hours of "naive" discussion that helped explore the interpretations of this model's impact on our perception of ageing, Dr. André Klarsfeld for his numerous useful comments on the manuscript and Dr. Bastian Greshake Tsovaras for helping with setting up the mybinder.org implementations of the simulations' codes. This work was granted access to the GENCI-sponsored HPC resources of TGCC@CEA under allocation A0090607519. This work has been supported by the Chair "Modélisation Mathématique et Biodiversité" of Veolia Environnement-Ecole Polytechnique-Museum National d'Histoire Naturelle-Fondation X.

Authors' Contribution

RT wrote the Python code presented in Annexe 2 and developed the mathematical analysis presented in Annexe 1, JP translated the Python code into C and ran the broad trait simulations. MS verified the mathematical analysis from Annexe 1, developed the analysis presented in Annexe 2 and wrote the manuscript. RM designed the study, designed the figures and wrote the manuscript.


MD Maria Betânia Ferreira

MD Tomás Fonseca

MD Rita Costa

MD António Marinho

MD José Carlos Oliveira

MD, PhD Faiez Zannad

MD, PhD Patrick Rossignol

MD Patrícia Rodrigues

PhD António S Barros

; João

MD, PhD Pedro Ferreira

Dra Maria

Betânia Ferreira
email: betania_ferreira82@hotmail.com

Sex differences in circulating proteins of patients with Rheumatoid Arthritis: a cohort study

Keywords: rheumatoid arthritis, sex differences, proteomic biomarkers, pathophysiology

Objetives: Differences in proteomic profiles between men and women may provide insights into the biological pathways that contribute to known sex differences in Rheumatoid Arthritis (RA). Studies focusing on sex differences in circulating proteins in RA patients are scarce. Our objective was to investigate the sex differences in circulating proteins of RA patients.

Methods: Cohort study enrolling 399 RA patients. 94 circulating proteinbiomarkers (92CVDIIOlink®+troponin-T+c-reactive protein) were measured.

Clinical, demographic, and echocardiographic characteristics were compared between men and women. The sex differences in the biomarker's expression were assessed using regression modelling.

Results: 306 (76.7%) patients were women. Compared with men, women had less visceral fat, smoked less, had diabetes and COPD less frequently and expressed more fatigue, anxiety and depression. The association with cardiovascular outcomes did not differ between sexes. After adjusting for potential confounders, women expressed higher levels of circulating proteins related to adipokine signaling and vascular function (e.g., leptin and vascularendothelial-growth factor), whereas men expressed higher levels of circulating proteins related to extracellular matrix organization and inflammation (e.g., matrix metalloproteinase-2 and c-reactive protein). These results were not found in patients without RA.

Conclusion: Sex differences in circulating proteins reflect distinct pathways implicated in pathogenesis of RA, including inflammation, adiposity, angiogenesis, and extracellular matrix organization. These findings may help further investigations into factors underlying sex-based differences and allow future studies focused on sex specific personalized treatment approaches in RA.

ClinicalTrials.

Introduction

Rheumatoid arthritis (RA) incidence is three times higher among women compared with men and has increased over the last decades [START_REF] Wolfe | The epidemiology of rheumatoid arthritis: a review. II. Incidence and diagnostic criteria[END_REF][START_REF] Myasoedova | Is the incidence of rheumatoid arthritis rising?: results from Olmsted County, Minnesota, 1955-2007[END_REF]. In women, RA incidence increases linearly with age; however, men tend to have a later RA onset, with a higher risk after the sixth decade of life [START_REF] Symmons | The incidence of rheumatoid arthritis in the United Kingdom: results from the Norfolk Arthritis Register[END_REF]. Genetic, hormonal, and environmental/behavioral factors may explain, at least partly, why women have higher susceptibility and prevalence of RA [START_REF] Favalli | Sex and Management of Rheumatoid Arthritis[END_REF].

Furthermore, RA presentation may also differ between women and men, with many studies suggesting that women have higher disease activity, lower proportion of achieved remission and poorer treatment response [START_REF] Iikuni | The influence of sex on patients with rheumatoid arthritis in a large observational cohort[END_REF][START_REF] Sokka | Women, men, and rheumatoid arthritis: analyses of disease activity, disease characteristics, and treatments in the QUEST-RA study[END_REF][START_REF] Kvien | A Norwegian DMARD register: prescriptions of DMARDs and biological agents to patients with inflammatory rheumatic diseases[END_REF][START_REF] Forslind | Sex: a major predictor of remission in early rheumatoid arthritis?[END_REF][START_REF] Symmons | Epidemiology of rheumatoid arthritis: determinants of onset, persistence and outcome[END_REF][START_REF] Flouri | Comparative effectiveness and survival of infliximab, adalimumab, and etanercept for rheumatoid arthritis patients in the Hellenic Registry of Biologics: Low rates of remission and 5-year drug survival[END_REF][START_REF] Tengstrand | The influence of sex on rheumatoid arthritis: a prospective study of onset and outcome after 2 years[END_REF][START_REF] Kuiper | Influence of sex, age, and menopausal state on the course of early rheumatoid arthritis[END_REF]. Some of the reported differences may depend on how RA disease activity is measured rather than the RA disease activity itself, because women generally report more symptoms and poorer scores on most questionnaires [START_REF] Van Vollenhoven | Sex differences in rheumatoid arthritis: more than meets the eye[END_REF].

Given the different incidence, presentation, and prognosis of RA in men and women (including cardiovascular events), it is likely that the mechanisms involved in the pathogenesis of RA are sex-specific [START_REF] Elemam | Role of Chemokines and Chemokine Receptors in Rheumatoid Arthritis[END_REF][START_REF] Burska | Cytokines as biomarkers in rheumatoid arthritis[END_REF][START_REF] Crowson | Impact of risk factors associated with cardiovascular outcomes in patients with rheumatoid arthritis[END_REF][START_REF] Lew | Sex-Based Differences in Cardiometabolic Biomarkers[END_REF][START_REF] Despres | Body fat distribution and risk of cardiovascular disease: an update[END_REF]. Therefore, a better understanding of the sex differences in RA regarding pattern and burden of disease and cardiovascular outcomes is needed. A multimarker approach may potentially provide insight on the biological processes underlying the sex differences in patients with RA and help ascertaining if the observed sex differences are due to distinct underlying mechanisms or to different clinical expression of the disease. Another important question is whether the sex differences in the biological processes are specific of RA or also common to patients without RA.

Using a prospective cohort with detailed phenotyping and multiple circulating proteins, this study aims to investigate the sex differences in RA patients from a clinical and biological perspective, and also to assess whether the differences in protein expression (and underlying biological pathways) are specific to RA.

Methods

Study design

We conducted a prospective single-center cohort study with patients attending the Autoimmune Disease Outpatient Unit of Centro Hospitalar do Porto, Portugal. The study was registered in clinicaltrials.gov under the number NCT03960515 [START_REF] Ferreira | Prevalence, risk factors and proteomic bioprofiles associated with heart failure in rheumatoid arthritis: The RA-HF study[END_REF].

This study was conducted following the principles of the Declaration of Helsinki and approved by the Hospital de Santo António, Porto, Portugal ethics committee under number #2016-023 (020-DEFI/020-CES). All participants provided the written informed consent before enrollment in the study.

The data used for this study are available upon reasonable request to the corresponding author.

Patient selection

We included all adult patients (18 years or older) with a diagnosis of RA and that were followed-up in our specialized RA clinic. All RA patients fulfilled the 2010 ACR/EULAR Classification Criteria [START_REF] Villeneuve | ACR-EULAR classification criteria for rheumatoid arthritis[END_REF]. Patients with active neoplasm or other severe co-morbid conditions potentially associated with life expectancy shorter than six months, as well as severe dementia, inability to walk or totally dependent on a third person were excluded from our study.

The screening period for patient enrollment occurred from June 2016 to June 2018. All enrolled patients maintained their usual care under the guidance of their treating physician.

Patients were followed-up for a median of 1.5 (0.7 to 2.3) years. The cause of death and hospitalization were independently adjudicated.

Patient evaluation, echocardiography and routine laboratory tests

We collected demographic data, cardiovascular history and risk factors, comorbid conditions, RA history, daily medication intake and a cardiovascular symptom questionnaire. Anxiety and depression were assessed using the Portuguese version of the Hospital Anxiety and Depression Scale [START_REF] Pais-Ribeiro | Validation study of a Portuguese version of the Hospital Anxiety and Depression Scale[END_REF]. Physical examination included vital signs, pulmonary and cardiac auscultation, and anthropometric measures. Visceral fat was assessed using a bioelectrical impedance analysis (BIA) with Tanita® BC545 instrument. Six-minute walking test was performed according to the American Thoracic Society protocol [START_REF] Atscopsfcpf | ATS statement: guidelines for the six-minute walk test[END_REF].

The transthoracic echocardiogram (TTE) was performed by an experienced echocardiographer and included the M-mode, 2D and Doppler measurements acquired according to standard recommendations [START_REF] Lang | Recommendations for cardiac chamber quantification by echocardiography in adults: an update from the American Society of Echocardiography and the European Association of Cardiovascular Imaging[END_REF][START_REF] Lang | Recommendations for cardiac chamber quantification by echocardiography in adults: an update from the American Society of Echocardiography and the European Association of Cardiovascular Imaging[END_REF].

Routine laboratory tests (blood and urine samples) included a complete blood count, c-reactive protein (CRP), glucose, blood electrolytes, lipid profile, creatinine, high sensitivity troponin T (hs-TnT), and N-Terminal-pro Brain Natriuretic Peptide (NT-pro BNP).

Plasma Proteomic Biomarkers

We measured a large protein biomarker panel (Olink® CVDII panel) that comprised 92 biomarkers associated with cardiovascular and inflammatory of pathophysiological domains. Results are expressed in normalized NPX units (https://www.olink.com/). The overview of the 92 circulating proteins is depicted in the Supplement Table 1, including their full names and Uniprot ID. In the Supplement Table 2 the proteins are described according to their main biological functions. The full description and technical details of the analytical assay are described in the Supplemental Addenda.

Plasma concentrations of high-sensitivity troponin T (hsTnT) and CRP were also measured using the Elecsys Troponin T Gen 5 STAT test by Roche Diagnostics® and the enzyme-linked immunosorbent assay (Olympus CRP Latex Calibrator Normal Set®), respectively. The CVDII Olink® panel displayed a mean intra-assay and inter-assay coefficient of variation between 9.1% and 11.7%. The hs-TnT and CRP showed an intra-assay and inter-assay coefficient of variation was 10% at a concentration lower than the 99 th percentile cut-off.

External assessment in patients without rheumatoid arthritis

To assess if our findings were unique to RA or whether if the sex differences in the biological processes are common to patients without RA we externally assessed our findings in two independent case-control cohorts of patients with cardiovascular risk factors, but without RA, both part of the European granted program Heart Omics and AGEing (HOMAGE). These cohorts have been previously described and the full results published [START_REF] Ferreira | Proteomic Bioprofiles and Mechanistic Pathways of Progression to Heart Failure[END_REF]. Briefly, a nestedmatched case-control design was used with cases (patients subsequently progressing to incident HF), and controls (patients without progression to HF). Controls were matched on cohort, follow-up time, age, and sex. Two independent sample sets (Ia with 286 cases and 591 controls and Ib with 276 cases and 280 controls) were used. The analysis presented in this report are fully adjusted on the matching variables age, heart failure case, follow-up time, smoking, diabetes mellitus, history of coronary artery disease, body mass index, systolic blood pressure, use of antihypertensive medication, and heart rate.

Statistical and bioinformatical analysis

Baseline clinical, demographic, and echocardiographic characteristics were compared between male and female patients using chi-square, t test, or Mann-Whitney tests, as appropriate. The level of significance for all hypothesis tests (P value) was set at 0.05. The normality of the distribution of continuous variables was evaluated by Kolmogorov-Smirnov's test and Shapiro-Wilk's test. Continuous variables are presented as mean and standard deviation (if normally distributed) or median and 25-75% quartiles (if non-normally distributed). We used the student's t-test to examine the differences of the independent continuous variables between male and female groups when normally distributed; the Mann-Whitney U was used for variable without normal distribution. Categorical variables are presented as absolute (n) and relative (%) frequencies. To study whether biomarker patterns differ between male and female patients, we first identified which clinical variables better differentiated sex using stepwise (forward) logistic regression analysis including all variables with a p-value <0.1 from the Table 1 ("clinical model"). To identify the circulating proteins that better differentiated sex, we performed logistic regression analysis fully adjusted for the "clinical model" and corrected for test multiplicity using a false discovery rate (FDR) of 5% or less. As means of internal validation we Network analyses according to sex were performed using the STRING database (Search Tool for the Retrieval of Interacting Genes/Proteins) and the involved pathways were over-represented using the KEGG (Kyoto Encyclopedia of Genes and Genomes) tool.

Results

Patient selection and baseline characteristics

A total of 408 patients with RA were screened between June 2016 and June 2018. Seven patients had exclusion criteria and 2 did not perform biomarker assessment. The remaining 399 patients fulfilled the eligibility criteria and performed the biomarker study. The baseline characteristics comparing men and women are depicted in Table 1. Of the total 399 RA patients, 306 (76.7%) were women. Compared with men, women had less visceral fat, smoked less, and had diabetes and COPD less often. Women expressed more fatigue, anxiety and depression (p <0.05 for all). Men and women displayed similar history of RA and receive similar treatment for either disease-modifying antirheumatic drugs (biologic or synthetic DMARDs) and non-DMARDs. In general, men expressed more features associated with a higher cardiovascular risk such as high blood pressure, lower HDL cholesterol, higher CRP levels, lower eGFR, and lower left ventricular ejection fraction. Table 1.

Association between proteomic biomarkers and male sex

On top of the clinical model presented in the Supplemental Table 3, that included age, waist circumference, smoking, diabetes, diastolic blood pressure, hemoglobin, and eGFR as the variables with strongest association with sex (AUC =0.89), of the 94 proteomic biomarkers, 13 were retained, as differentially expressed according to sex, with an FDR <5%. Table 2. Of these 13 proteins, 6 had higher concentrations expressed in women: lipoprotein lipase (LPL), growth/differentiation factor 2 (GDF2), leptin (LEP), vascular endothelial growth factor D (VEGFD), interleukin-1 receptor-like 2 (IL1RL2) and a disintegrin and metalloproteinase with thrombospondin motifs 13 (ADAMTS13), and 7 had higher concentrations expressed in men: CRP, troponin-T, lactoylglutathione lyase (GLO1), placenta growth factor (PGF), C-C motif chemokine 3 (CCL3), interleukin-18 (IL18) and matrix metalloproteinase-12 (MMP12). The network analysis representation of the findings and the interplay between these proteins among men and women is displayed in the Figure 1. The full biomarker regression analysis is displayed in the Supplemental Table 4.

The LASSO model largely confirmed the same biomarkers that better differentiated sexes, indicating excellent internal validity of the data.

Supplemental Table 5.

Proteomic external assessment in patients without RA

The external assessment in patients without RA was generally not consistent with the finding of RA patients. In women, only VEGF and ADAMTS13 were directionally and consistently replicated in HOMAGE phases 1a & 1b. In men, no biomarkers were replicated in both phases (CRP and hs-TnT were not measured in HOMAGE). Supplemental Table 6.

Cardiovascular outcomes

There was no difference in the rates of the composite endpoint of cardiovascular death or cardiovascular hospitalization between women and men: events 37/306 (12.1%) in women vs. 19/93 (20.4%) in men, adjusted HR 1.13, 95% CI 0.56-2.31, p-value =0.73. Figure 2.

Discussion

The present study, including 399 RA patients of whom 77% women, shows that women had less visceral fat, were less frequently smokers, and had less comorbid conditions and cardiovascular risk factors (including than men.

However, women experienced more fatigue, anxiety and depression than men.

Women expressed higher levels of circulating proteins related to adipokine signaling and vascular function (e.g., LPL, LEP, GDF2, VEGFD, ADAMTS13), whereas men expressed higher levels of circulating proteins related to ischemia, extracellular matrix organization and inflammation (e.g., hs-TnT, MMP12, CCL3, CRP, IL18, GLO1), suggesting that men and women may have sex-specific biological pathways involved in RA pathogenesis. Furthermore, most of the sex differences in the expression of investigated circulating proteins were not observed in control patients with CV risk factors suggesting that these might be specific to RA.

This work aimed at substantiating potential differences in mechanistic pathways which ay underly the many differences in incidence, presentation, and prognosis of RA in men and women. We examined whether a multiple biomarker approach may detect sex-specific biological pathways involved in RA, providing basis for a better understanding of the disease processes which may help conceptualizing individualized treatments.

To the best of our knowledge, this is the first report assessing a large panel of circulating proteins to study the sex differences in protein expression in RA. For the majority of the circulating proteins found in this study, we could not find any previous study in RA populations allowing to compare our sex-specific findings. However, data from non-RA populations and animal models may serve as basis to support and reinforce our findings. Adipokines usually have higher expression in women [START_REF] Lau | Sex Differences in Circulating Biomarkers of Cardiovascular Disease[END_REF]. In RA, leptin has pro-inflammatory effects and correlates with measures of disease activity [START_REF] Toussirot | The role of leptin in the pathophysiology of rheumatoid arthritis[END_REF][START_REF] Fatel | Adipokines in rheumatoid arthritis[END_REF][START_REF] Skalska | Adiponectin Isoforms and Leptin Impact on Rheumatoid Adipose Mesenchymal Stem Cells Function[END_REF]. Our study suggests that leptin may be involved in RA pathogenesis specifically in women, independently of abdominal obesity and diabetes. Lipoprotein lipase (LPL), also with higher expression among women in our study, has the dual function of hydrolyzing triglycerides and serve as a ligand factor for receptor-mediated lipoprotein uptake. A study with RA patients has shown that the higher level of LPL activity in this population might result from the inhibition of LPL degradation by corticosteroids [START_REF] Kidani | Active-dimeric form of lipoprotein lipase increases in the adipose tissue of patients with rheumatoid arthritis treated with prednisolone[END_REF]. However, this study does not provide sex-specific associations and in our study the use of corticosteroids was similar between men and women. Therefore, our data suggest that LPL may have a sex-specific role in RA, deserving further investigation. Women also expressed higher levels of angiogenesis-related markers (GDF2, VEGFD and ADAMTS13). Growth differentiation factor 2 (or bone morphogenetic protein 9, BMP9) is a member of the TGF-β family and a potent circulating inhibitor of angiogenesis [START_REF] Morrell | Targeting BMP signalling in cardiovascular disease and anaemia[END_REF]. Furthermore, GDF2 signaling increases vascular inflammation through direct effects on the endothelium [START_REF] Helbing | Bone Morphogenetic Protein-Modulator BMPER Regulates Endothelial Barrier Function[END_REF]. Therefore, blocking the GDF2/BMP9 pathways may be a target for reducing vascular inflammation [START_REF] Yao | Inhibition of bone morphogenetic proteins protects against atherosclerosis and vascular calcification[END_REF]. The VEGFD also plays an important role in angiogenesis by stimulating the proliferation and migration of endothelial cells and is upregulated by proinflammatory cytokines and hypoxia. Circulating VEGF concentrations are elevated in RA and correlate with disease activity and progression [START_REF] Ramirez | Patients with rheumatoid arthritis in clinical remission and ultrasound-defined active synovitis exhibit higher disease activity and increased serum levels of angiogenic biomarkers[END_REF][START_REF] Taylor | VEGF and imaging of vessels in rheumatoid arthritis[END_REF]. ADAMTS13 controls the von Willebrand Factor (vWF) -mediated platelet thrombus formation by cleaving the vWF multimers into smaller forms; therefore, preventing microvascular thrombosis and inflammation [START_REF] Lu | ADAMTS13 ameliorates inflammatory responses in experimental autoimmune encephalomyelitis[END_REF]. The IL1RL2 is a receptor for interleukin-36 mediating local (including joint) inflammatory responses [START_REF] Boutet | Interleukin-36 family dysregulation drives joint inflammation and therapy response in psoriatic arthritis[END_REF]. Compared with women, men with RA exhibit higher concentrations of systemic inflammation markers and circulating proinflammatory cytokines/chemokines (CRP, IL-18, CCL3 and GLO1) [START_REF] Zhang | CCL3 participates in the development of rheumatoid arthritis by activating AKT[END_REF][START_REF] Joosten | Association of interleukin-18 expression with enhanced levels of both interleukin-1beta and tumor necrosis factor alpha in knee synovial tissue of patients with rheumatoid arthritis[END_REF][START_REF] Imas | Optical Biosensors for the Detection of Rheumatoid Arthritis (RA) Biomarkers: A Comprehensive Review[END_REF], tissue injury, fibrosis and remodeling (MMP12) [START_REF] Wang | Overexpression of human matrix metalloproteinase-12 enhances the development of inflammatory arthritis in transgenic rabbits[END_REF], cardiac injury (hs-TnT) [START_REF] Ferreira | Prevalence, risk factors and proteomic bioprofiles associated with heart failure in rheumatoid arthritis: The RA-HF study[END_REF], and synoviocyte proliferation (PGF) [START_REF] Kumkumian | Platelet-derived growth factor and IL-1 interactions in rheumatoid arthritis. Regulation of synoviocyte proliferation, prostaglandin production, and collagenase transcription[END_REF]. Together, our data suggest that adipokines and endothelial/vascular inflammation play a central role in the pathogenesis of RA in women, whereas in men the pathogenesis of RA may be more related to systemic inflammation and tissue injury [START_REF] Forslind | Sex: a major predictor of remission in early rheumatoid arthritis?[END_REF].

The association with cardiovascular outcomes did not differ between sexes, but the number of events was overall low. Nonetheless, men experienced nominally more cardiovascular events than women, which may reflect a greater cardiac involvement in men with RA (also supported by the higher concentration of hs-TnT in men).

Limitations

Several limitations should be acknowledged in our study. This is an observational study and thus causality cannot be ascertained. Our cohort came from a single center and some of the findings may reflect specific sociodemographic characteristics and local practice patterns which may limit the generalization of the findings. The lower representativeness of men (with a ratio of 1:3 vs. women, compatible with the epidemiology of the disease) and a small number of events, decreased the statistical power for comparing outcomes between sexes. .We have studied the Olink® CVDII panel, hs-TnT and CRP.

Whether a larger protein screening could provide more information might worth exploring. Our study included ambulatory patients only, and our findings cannot be generalized to hospitalized RA patients. Finally, studies in larger populations to validate our results are needed.

Conclusion

Sex differences in circulating proteins reflect distinct pathways implicated in pathogenesis of RA, adipokines and endothelial/vascular inflammation play a central role in the pathogenesis of RA in women, whereas in men the pathogenesis of RA may be more related to systemic inflammation and tissue injury. Our study provides novel hypotheses and highlights the need for a better understanding of the sex differences in RA. VEGFD, Vascular endothelial growth factor D; GDF2, Growth/differentiation factor 2; LEP, Leptin; IL1RL2, Interleukin-1 receptor-like 2; CRP, c-reactive protein; CCL3, C-C motif chemokine 3; MMP12, Matrix metalloproteinase-12; IL18, Interleukin-18; GLO1, Lactoylglutathione lyase; hs-TnT, high sensitivity troponin T; PGF, Placenta growth factor. *All the biomarkers had a FDR lower than 5% after adjustment (see the methods section for details).

The biomarkers are ordered by strength of the association (p-value).

To facilitate the interpretation of the data, we present all the biomarkers in terms of "higher association" with women and men i.e., the inverse odds ratio (1/OR) is presented for women by changing the indicator variable. Model adjusted on age, waist circumference, smoking, diabetes, diastolic blood pressure, hemoglobin, and estimated glomerular filtration rate.

  have also performed a 10-fold cross-validation least absolute shrinkage and selection operator (LASSO) regression analysis. The association between sex and the composite endpoint of cardiovascular death or cardiovascular hospitalizations was performed with a Cox model. All statistical analyses were performed with Stata® (version 16; StataCorp 2019, Stata Statistical Software: Release 16; College Station, TX: StataCorp LP).

Figure 1 .

 1 Figure 1. Network analysis of protein expression in men and women

Figure 2 .

 2 Figure 2. Composite endpoint of CV death or CV hospitalization between women and men

Table 1 .

 1 Baseline characteristics of men and women.

	Baseline characteristics	Women (N=306)	Men (N=93)	P
				value
	Biometric and demographic
	Age (years)	58.0 ± 13.4	61.0 ± 12.3	0.056
	BMI (Kg/m 2 )	26.5 ± 4.5	26.6 ± 4.3	0.93
	Waist circumference category*	157 (51.3%)	28 (30.1%)	<0.001
	Visceral fat	7.8 ± 3.0	11.9 ± 4.1	<0.001
	Smoking status
	Never	235 (76.8%)	24 (25.8%)
	Past	39 (12.7%)	45 (48.4%)	<0.001
	Current	32 (10.5%)	24 (25.8%)
	Comorbidities
	Diabetes mellitus	33 (10.8%)	20 (21.5%)	0.008
	Dyslipidemia	152 (49.7%)	38 (40.9%)	0.14
	Myocardial infarction	3 (1.0%)	3 (3.2%)	0.12
	Angina pectoris	22 (7.2%)	3 (3.2%)	0.17
	Ischemic stroke or transient	10 (3.3%)	4 (4.3%)	0.64
	ischemic attack
	Valvular heart disease	9 (2.9%)	1 (1.1%)	0.31
	Atrial fibrillation	11 (3.6%)	3 (3.2%)	0.87
	COPD	9 (2.9%)	9 (9.7%)	0.006
	Prior heart failure	18 (5.9%)	6 (6.5%)	0.84
	Hypertension	133 (43.5%)	51 (54.8%)	0.054
	Anemia	59 (19.3%)	11 (11.8%)	0.098
	Chronic kidney disease	12 (3.9%)	1 (1.1%)	0.18
	RA history
	Time since RA diagnosis	12.5 ± 10.8	10.8 ± 9.0	0.18
	(years)
	RF or anti-CCP positive	241 (78.8%)	75 (80.6%)	0.69
	Articular erosions	105 (37.6%)	34 (41.5%)	0.53
	DAS28 value with ESR	2.9 ± 1.2	2.6 ± 1.3	0.13
	Symptoms and physical exam
	Dyspnea	91 (29.7%)	22 (23.7%)	0.25
	Fatigue	196 (64.1%)	46 (49.5%)	0.012
	Systolic blood pressure (mmHg)	132.5 ± 19.6	137.2 ± 18.4	0.042
	Diastolic blood pressure	75.2 ± 10.0	79.7 ± 10.4	<0.001
	(mmHg)
	Heart rate (bpm)	80.2 ± 13.2	82.2 ± 18.0	0.24
	Echocardiogram
	LVEF (biplane)	61.7 ± 6.9	58.7 ± 7.3	<0.001
	LVEF <50%	12 (4.4%)	8 (10.0%)	0.059
	Left ventricular hypertrophy	28 (10.3%)	8 (9.5%)	0.85
	Left atrium volume (biplane	32.6 ± 11.2	31.8 ± 9.7	0.57
	indexed)

Table 2 .

 2 Circulating proteins that better differentiate between sexes. OR. Odds ratio; CI, confidence interval. ADAMTS13, Adisintegrin and metalloproteinase with thrombospondin motifs 13; LPL, Lipoprotein lipase;

	Biomarker	OR (95%CI)	P value*
	Higher expression in women
	LPL	8.33 (3.70-20.0)	0.0001
	GDF2	4.76 (2.50-9.09)	0.0001
	LEP	4.00 (2.63-5.88)	0.0001
	VEGFD	7.14 (2.43-20.0)	0.0003
	IL1RL2	3.84 (1.81-8.33)	0.0004
	ADAMTS13	33.3 (3.03-99.0)	0.0037
	Higher expression in men
	CRP	1.68 (1.32-2.14)	0.0001
	hs-TnT	3.14 (1.76-5.61)	0.0001
	GLO1	2.64 (1.55-4.50)	0.0003
	PGF	5.27 (2.05-13.6)	0.0006
	CCL3	1.82 (1.28-2.58)	0.0008
	IL18	2.26 (1.29-3.98)	0.0044
	MMP12	1.88 (1.20-2.95)	0.0055
	Legend:

Acknowledgements

JPF, PR, and FZ are supported by the RHU Fight-HF, a public grant overseen by the French National Research Agency (ANR) as part of the second "Investissements d'Avenir" program (reference: ANR-15-RHUS-0004), by the French PIA project "Lorraine Université d'Excellence" (reference: ANR-15-IDEX-04-LUE), Contrat de Plan Etat Lorraine IT2MP and FEDER Lorraine.

Disclosures

None.

Funding

None.

Authors' Contributions

M.B.F., drafted the manuscript; J.P.F., performed the statistical analysis and drafted the manuscript; A.S.B., performed the bioinformatical analysis; all authors read, provided critical input and approved the final version of the manuscript.

Dr. João Pedro Ferreira is the guarantor of this work and, as such, had full access to all the data in the study and takes responsibility for the integrity of the data and the accuracy of the data analysis.