Jeanne Bovet
email: jeanne.bovet@iast.fr

Mélissa Barkat-Defradas

Valérie Durand

Charlotte Faurie

Michel Raymond

WOMEN'S ATTRACTIVENESS IS LINKED TO EXPECTED AGE AT MENOPAUSE (WOMEN'S ATTRACTIVENESS AND AGE AT MENOPAUSE)

Keywords: humans, sexual selection, mate choice, female attractiveness, face, menopause, residual reproductive value, fertility

A great number of studies have shown that features linked to immediate fertility explain a large part of the variance in female attractiveness. This is consistent with an evolutionary perspective, as men are expected to prefer females at the age at which fertility peaks (at least for short-term relationships) in order to increase their reproductive success. However, for long-term relationships, a high residual reproductive value (the expected future reproductive output, linked to age at menopause) becomes relevant as well. In that case, young age AND late menopause are expected to be preferred by men. However, the extent to which facial features provide cues to the likely age at menopause has never been investigated so far. Here, we show that expected age at menopause is linked to facial attractiveness of young women.

As age at menopause is heritable, we used the mother's age at menopause as a proxy for her daughter's expected age of menopause. We found that men judged faces of women with a later expected age at menopause as more attractive than those of women with an earlier expected age at menopause. This result holds when age, cues of immediate fertility and facial ageing were controlled for. Additionally, we found that the expected age at menopause was not correlated with any of the other variables considered (including immediate fertility cues and facial ageing). Our results show the existence of a new correlate of women's facial attractiveness, expected age at menopause, which is independent from immediate fertility cues and facial ageing.

Introduction

Why some women are more attractive than others has been the focus of many studies, and it has been repeatedly shown that features linked to immediate fertility (i.e. fertility at a given time) explain a large part of the variance in attractiveness (see [START_REF] Rhodes | The evolutionary psychology of facial beauty[END_REF] for a review).

Cues of immediate fertility can be found in body shape, voice pitch, or facial features; in fact, all traits linked to sexual maturity, age and parity or substantially influenced by sex hormones (which are factors influencing immediate fertility) have the potential to influence physical attractiveness (see, for example, [START_REF] Symons | Beauty is in the adaptations of the beholder: The evolutionary psychology of human female sexual attractiveness[END_REF][START_REF] Henss | Waist-to-hip ratio and female attractiveness. Evidence from photographic stimuli and methodological considerations[END_REF][START_REF] Singh | Female mate value at a glance: Relationship of waist-to-hip ratio to health, fecundity and attractiveness[END_REF][START_REF] Jasienska | Large breasts and narrow waists indicate high reproductive potential in women[END_REF][START_REF] Feinberg | The voice and face of woman: One ornament that signals quality?[END_REF][START_REF] Smith | Facial appearance is a cue to oestrogen levels in women[END_REF][START_REF] Pipitone | Women's voice attractiveness varies across the menstrual cycle[END_REF][START_REF] Singh | Cross-cultural consensus for waist-hip ratio and women's attractiveness[END_REF][START_REF] Little | Facial attractiveness: evolutionary based research[END_REF][START_REF] Puts | Sexual Selection on Human Faces and Voices[END_REF][START_REF] Puts | Women's attractiveness changes with estradiol and progesterone across the ovulatory cycle[END_REF][START_REF] Pfluger | Cues to fertility: perceived attractiveness and facial shape predict reproductive success[END_REF][START_REF] Jones | Agreement and Individual Differences in Men's Preferences for Women's Facial Characteristics[END_REF][START_REF] Mondragón-Ceballos | Waist-to-Hip Ratio, but Not Body Mass Index, Is Associated with Testosterone and Estradiol Concentrations in Young Women[END_REF][START_REF] Sugiyama | Physical Attractiveness: An Adaptationist Perspective[END_REF][START_REF] Butovskaya | Waist-to-hip ratio, body-mass index, age and number of children in seven traditional societies[END_REF] but see [START_REF] Marcinkowska | Cross-cultural variation in men's preference for sexual dimorphism in women's faces[END_REF][START_REF] Marcinkowska | Men's Preferences for Female Facial Femininity Decline With Age[END_REF] for some modulations of these preferences). To increase their reproductive success, males should prefer females of the age at which age-specific fertility peaks, at least for short-term mates (Maestripieri et al. 2014). For long-term relationships such as marriage, however, residual reproductive value -the expected future reproductive output -becomes pertinent too. In humans, this trend is strengthened by the existence of reproductive senescence or menopause (i.e., the permanent cessation of menstruation), associated with the ultimate cessation of child-bearing potential, long before the somatic senescence.

Consequently, for long-term relationships, the number of offspring produced by the couple will depend in part on the number of years before the woman reach menopause, particularly in natural fertility populations where reproduction could theoretically extend until menopause. Thus, men can increase their reproductive success by choosing a long-term mate with a longer reproductive window (i.e. higher residual reproductive value). The first criterion of a long reproductive window is of course age. And indeed, youth is one of the most important factor of women's attractiveness [START_REF] Buunk | Age preferences for mates as related to gender, own age, and involvement level[END_REF]. Then, for a given age, the temporal reproductive window will vary according to the age at menopause, which is highly variable both across and within populations [START_REF] Avis | Is there a menopausal syndrome? Menopausal status and symptoms across racial/ethnic groups[END_REF][START_REF] Thomas | International variability of ages at menarche and menopause: Patterns and main determinants[END_REF][START_REF] Velde | The variability of female reproductive ageing[END_REF][START_REF] Dratva | Is age at menopause increasing across Europe? Results on age at menopause and determinants from two population-based studies[END_REF][START_REF] Morris | Familial concordance for age at natural menopause: results from the Breakthrough Generations Study[END_REF][START_REF] Stepaniak | Age at natural menopause in three Central and Eastern European urban populations: The HAPIEE study[END_REF]. If the future age at menopause is somehow detectable in young adults, it could contribute substantially to mate choice and thus influence female attractiveness. It is possible that previously identified components of immediate fertility are also informative of female residual reproductive value: reproductive value at a young age and later in life could be positively correlated due to some common underlying factors, or negatively correlated as a result of a biological trade-off [START_REF] Hamilton | The moulding of senescence by natural selection[END_REF][START_REF] Wood | The evolution of menopause by antagonistic pleiotropy[END_REF][START_REF] Carter | Antagonistic pleiotropy as a widespread mechanism for the maintenance of polymorphic disease alleles[END_REF]. Alternatively, future age at menopause could be independent of immediate fertility cues, and thus a correlate of attractiveness that has not yet been identified.

In this study, we investigated whether expected age at menopause is related to the attractiveness of young women's faces. As age at menopause is heritable (heritability between 0.42 and 0.72 depending on the study: [START_REF] Snieder | Genes Control the Cessation of a Woman's Reproductive Life: A Twin Study of Hysterectomy and Age at Menopause[END_REF][START_REF] De Bruin | The role of genetic factors in age at natural menopause[END_REF][START_REF] Van Asselt | Heritability of menopausal age in mothers and daughters[END_REF][START_REF] Murabito | Heritability of age at natural menopause in the Framingham Heart Study[END_REF][START_REF] Morris | Familial concordance for age at natural menopause: results from the Breakthrough Generations Study[END_REF]), mother's age at menopause was used as a proxy for the daughter's expected age of menopause. To better understand how the residual reproductive value could be facially detectable, immediate fertility cues were included in the model as possible explanatory variables. It could also be speculated that women who reach menopause later have an overall slower life history trajectory, with a slower general physical ageing, and may look younger, at the same age, than women who will experience menopause at an earlier age (in line with this idea, a later age at menopause is associated with a longer life expectancy, for a review see [START_REF] Gold | The Timing of the Age at Which Natural Menopause Occurs[END_REF]. To see if general physical ageing could mediate the relationship between facial attractiveness and expected age at menopause, we included facial ageing (measured as the difference between actual age and perceived age) as another explanatory variable.

Materials and Methods

Ethics Statement.

The protocol used to recruit participants and collect data was approved (#1226659) by the French National Committee of Information and Liberty (CNIL). For each participant, the general purpose of the study was explained ("a study on the determinants of mate choice") and a written voluntary agreement was requested for a statistical use of data (private information and photographs). Data were analysed anonymously.

Stimuli.

A total of 97 women between 25 and 35 years of age were recruited by social networks and advertising in Montpellier, France. Sixty-eight Caucasian women whose mother had a natural and known age at menopause constituted our final stimuli sample (mean age = 28.4, age range: 25-35, see table S4). Volunteers were instructed to come to the lab after collecting information about their mother's and (when possible) grandmothers' menopause (menopause was defined as the first full year without any menstruation) and without wearing any make-up. For each woman, the following information was collected: date and place of birth for themselves, their parents and grandparents; monthly income (divided into ten classes from less than 760€ to more than 4705€) for themselves and their parents; education level; and age at menopause for their mother and grandmothers. A facial photograph was taken with the same digital still camera (Canon EOS 20D) at a distance of 1 m using the same general settings. Each woman was asked to have a neutral facial expression, to remove any glasses or earrings, and to wear a hairband (to make sure that all the face was visible). All photographs were electronically processed using Adobe Photoshop CS3 to normalise size (photographs were aligned on eyes position, with a fixed distance between eyes and chin), colour balance, contrast and luminosity (using the Photoshop auto-corrections tools). Hairstyle was cropped, and the background was replaced by a uniform grey colour (see figure 1). A compensation of 20€ was provided for the subjects' participation.

Procedure.

A Delphi-based computer program was generated to present randomly drawn pairs of photographs of the 68 women (see figure 1). For each pair, the rater was instructed to click on the photograph depicting the woman he found the most attractive (the outcome measure of our study). The position of the photograph on the screen (left or right) was randomly ascribed. Each rater had 30 distinct pairs of photographs to assess, corresponding to 60 different women. If the rater knew one of the women he had to judge, the trial was removed. Also, the first pair seen by each participant was not used for the analyses, because the task could require some habituation. Three pairs, randomly chosen from among those previously viewed, were presented again at the end to estimate judgement reliability.

Raters.

A total of 156 male raters assessed the relative attractiveness of these women.

Volunteer raters were recruited in public places in Montpellier (France) and were unaware of the purpose of the study when assessing the pairs of pictures. For each rater, the following information was collected: date and place of birth, grandparents' origins, monthly income, occupation, house ownership, taxability, education level, and sexual orientation. Only data from Caucasian and heterosexual raters were used for the analyses. Assessments of unreliable raters (i.e., with more than one incorrect answer during the test of judgement reliability) were removed. A total of 119 raters were retained in the final sample, with a mean age of 36.2 (age range: 17-72, see table S4). Each woman was observed, on average, by 101.1 raters (range: 93-108).

Immediate fertility cues.

We collected 3 physical features hypothesised to be linked to women's immediate fertility: the waist-to-hip ratio (the ratio between body circumference at the waist and the hips or WHR, [START_REF] Singh | Body Shape and Women's Attractiveness -the Critical Role of Waist-to-Hip Ratio[END_REF][START_REF] Jasienska | Large breasts and narrow waists indicate high reproductive potential in women[END_REF][START_REF] Singh | Beauty is in the eye of the plastic surgeon: Waist-hip ratio (WHR) and women's attractiveness[END_REF][START_REF] Lassek | Menarche is related to fat distribution[END_REF][START_REF] Mondragón-Ceballos | Waist-to-Hip Ratio, but Not Body Mass Index, Is Associated with Testosterone and Estradiol Concentrations in Young Women[END_REF][START_REF] Butovskaya | Waist-to-hip ratio, body-mass index, age and number of children in seven traditional societies[END_REF]; but see [START_REF] Nenko | Fertility, body size, and shape: An empirical test of the covert maternal depletion hypothesis[END_REF], the fundamental frequency of the voice (F0, [START_REF] Harries | Changes in the male voice at puberty: vocal fold length and its relationship to the fundamental frequency of the voice[END_REF][START_REF] Abitbol | Sex hormones and the female voice[END_REF][START_REF] Feinberg | Are human faces and voices ornaments signaling common underlying cues to mate value?[END_REF][START_REF] Evans | The relationship between testosterone and vocal frequencies in human males[END_REF][START_REF] Abend | The sound of female shape: a redundant signal of vocal and facial attractiveness[END_REF]) and the facial femininity [START_REF] Farkas | Age-and sex-related changes in facial proportions[END_REF][START_REF] Johnston | Is Beauty in the Eye of the Beholder[END_REF][START_REF] Perrett | Effects of sexual dimorphism on facial attractiveness[END_REF][START_REF] Fink | The biology of facial beauty[END_REF][START_REF] Smith | Facial appearance is a cue to oestrogen levels in women[END_REF][START_REF] Little | Facial attractiveness: evolutionary based research[END_REF][START_REF] Pfluger | Cues to fertility: perceived attractiveness and facial shape predict reproductive success[END_REF][START_REF] Jones | Agreement and Individual Differences in Men's Preferences for Women's Facial Characteristics[END_REF]. The WHR was measured in the lab by the investigator. To measure the F0, women were recorded reading the French version of a standard text ("La bise et le soleil") using a Tascam DR-07 MKII digital recorder. Across each recording, the fundamental frequency (F0, the acoustic correlate of pitch) was measured using Praat software [START_REF] Boersma | Praat: doing phonetics by computer[END_REF]. To generate morphological facial femininity scores, a geometric morphometric analysis of the faces was used following methods described in [START_REF] Scott | Does Masculinity Matter? The Contribution of Masculine Face Shape to Male Attractiveness in Humans[END_REF][START_REF] Lee | Genetic Factors That Increase Male Facial Masculinity Decrease Facial Attractiveness of Female Relatives[END_REF][START_REF] Dixson | Beneath the beard: do facial morphometrics influence the strength of judgments of men's beardedness?[END_REF] First, the coordinates of 142 landmarks (anatomical points present in all individuals, e.g. the corners of the lips) and semi-landmarks (sliding points positioned along some anatomical curves, such as the bow of the eyebrow) were delineated for each female face, as well as for 26 male facial photographs retrieved from another database. These 26 additional men were recruited according to the same criteria as for the women's recruitment: heterosexuality, Caucasian origin, 25-35 years old. The delineation of the landmarks and semi-landmarks were done using Psychomorph [START_REF] Tiddeman | Towards Realism in Facial Image Transformation: Results of a Wavelet MRF Method[END_REF]. The R package Geomorph (version 3.0.3) was used to carry out a Procrustes superimposition of the landmark and semi-landmark data, which removes non-shape information such as translation, size, and rotational effects [START_REF] Zelditch | Geometric Morphometrics for Biologists: A Primer[END_REF][START_REF] Zelditch | A Practical Companion to Geometric Morphometrics for Biologists: Running analyses in freely-available software[END_REF]. A Principal Component Analysis (PCA) was conducted on the Procrustes-registered landmarks and semi-landmarks data of the 26 male faces and 26 female faces randomly drawn from this study dataset, matching the age distribution of the 26 men. This PCA produced shape variables which are a decomposition of the landmark coordinates of the male and female faces (see figure S1). The values on the factors of the PCA were computed for the remaining 42 female faces (not used to create the PCA). Then, a linear discriminant analysis (LDA) incorporating the two first components of the PCA was used to discriminate between male and female faces. The resulting analysis provided correct sex classification for 92.6% of faces. Discriminant function scores were therefore used as an index of facial femininity, with high scores indicating a more feminine facial morphology scores [START_REF] Scott | Does Masculinity Matter? The Contribution of Masculine Face Shape to Male Attractiveness in Humans[END_REF][START_REF] Lee | Genetic Factors That Increase Male Facial Masculinity Decrease Facial Attractiveness of Female Relatives[END_REF][START_REF] Dixson | Beneath the beard: do facial morphometrics influence the strength of judgments of men's beardedness?[END_REF]. All analyses were performed using R software, version 3.4.2.

Estimation of facial ageing.

It could be speculated that women who reach menopause later have an overall slower life history trajectory, with a slower general physical ageing, and may look younger, at the same age, than women who will experience menopause at an earlier age. To control for this possibility, a second Delphi-based computer program was generated to present the facial photographs of the 68 women to 136 raters. For each photograph, the rater was instructed to estimate the age of the woman (see figure S2). Each rater had 20 distinct photographs to assess, randomly drawn among the 68 pictures. Three photographs, randomly chosen among those previously viewed, were presented again at the end to estimate judgement reliability. For each rater, information about sex and age was collected. The reliability of raters was assessed by computing the sum of absolute differences between the 3 first estimations and the corresponding repetitions. Raters with values higher than 15 years were removed (a more stringent threshold of 9 years, or no threshold, did not change qualitatively the results). A total of 107 facial ageing raters (including 61 women) were retained in the final sample, with a mean age of 35.1 (age range: 16-65). Each woman was observed, on average, by 31.2 raters (range: 24-43). For each woman, the difference between their real age and the mean age given by the raters was calculated. This variable represents the facial ageing of the women, a critical explanatory variable for the analysis of attractiveness.

Statistical analyses.

Logistical regressions were used to analyse raters' attractiveness preferences. The binary response variable corresponded to being chosen or not for the focal woman (arbitrarily the woman presented at the left position) during the presentation of each pair. Women and attractiveness raters were considered random samples from a larger population of interest and were thus random-effect variables. Therefore, generalised linear mixed models with a binomial error structure were used. For each choice made by a rater, the difference between the ages at menopause of the focal and the non-focal woman's mothers was calculated. The value of this difference was integrated into the model as the main variable of interest. To control for potential confounding effects, the differences between the focal and the non-focal woman's ages and socio-economic status (SES, a PCA combination of education level and the woman's and her parents' monthly incomes) were introduced into the model. Because the subjects displayed a perceptible smile in some photographs, a qualitative variable describing this aspect was also introduced (this binary variable was coded by three independent raters, blind to the other characteristics of the women). Variables concerning the raters' characteristics were also included in the model as potential confounding effects. These variables were the rater's age and SES (a PCA combination of the variables "monthly income", "occupation", "house ownership", "taxability", and "education level"). All variables were standardized for the analysis. In a second model, the difference between the two women for cues of immediate fertility (facial femininity, WHR and F0), and facial ageing were included as explanatory variables. Indeed, these variables, which may be linked to expected age at menopause, could potentially mediate the relationship between expected age at menopause and facial attractiveness.

Linear regression was used to analyse the mother's age at menopause according to the maternal grandmother's age at menopause. The "one-parent -one child" regression coefficient represents half the heritability (h 2 ) of the trait [START_REF] Lynch | Genetics and analysis of quantitative traits[END_REF]. Pearson correlations were used to analyse relationships between each woman's measured traits: expected age at menopause, immediate fertility cues, facial ageing, age, and SES. All statistical analyses were performed using R software, version 3.4.2.

Results

Attractiveness, age at menopause and immediate fertility.

The variable expected age at menopause (estimated by mother's age at menopause) had a significant positive effect on the probability of a woman being chosen as the most attractive: Men tend to prefer women who are likely to reach menopause later (β = 0.22, SE = 0.053, P < 0.001 see model 1 in table 2). This result holds when immediate fertility cues and facial ageing are controlled for (β = 0.24, SE = 0.053, P < 0.001, see model 2 in table 1 and figure 2), showing that the effect of expected age at menopause on attractiveness is not due to know cues of immediate fertility cues or facial ageing. Additionally, the 3 cues of immediate fertility had a significant effect on the probability of being chosen: men tend to prefer more feminine faces (β = 0.4, SE = 0.057, P < 0.001, see model 2 in table 1), and the faces of women who have a higher (i.e. more feminine) F0 (β = 0.2, SE = 0.058, P = 0.001). In contradiction with our prediction, men also tend to prefer the faces of women who had a higher (i.e. more masculine) WHR (β = 0.13, SE = 0.057, P = 0.02). Facial ageing had a negative effect on the probability to be chosen (β = -0.56, SE = 0.06, P < 0.001): Independently of their actual age, women who look younger than their actual age were preferred by men. Women's age and smile demonstrated a significant effect on the probability of being chosen: men preferred women who were younger (β < -0.2, P < 0.001, see models 1 and 2 in table 1) and smiled more (β > 0.27, P < 0.001, see models 1 and 2 in table 1). The socio-economic status of the woman had a negative effect on the probability to be chosen as the more attractive, but this result was significant only in model 2 (β = -0.17, SE = 0.057, P = 0.003). The raters' age and socio-economic level had no significant effect on their choices (all P > 0.5, see models 1 and 2 in table 1).

Menopause and other women's features.

There was no significant correlation between expected age at menopause and immediate fertility cues (facial femininity, F0 and WHR), facial ageing, woman's age or socio-economic status (all P > 0.1, see table 2), showing that expected age at menopause is capturing a distinct correlate of facial attractiveness. There was no significant correlation between the other variables, except a negative correlation between facial ageing and SES (r = -0.29, P = 0.021, see table 2).

Heritability of age at menopause.

To help establish the validity of our sample, we measured the heritability of age at menopause, for comparison with previous studies done with various samples and methods [START_REF] Snieder | Genes Control the Cessation of a Woman's Reproductive Life: A Twin Study of Hysterectomy and Age at Menopause[END_REF][START_REF] De Bruin | The role of genetic factors in age at natural menopause[END_REF][START_REF] Van Asselt | Heritability of menopausal age in mothers and daughters[END_REF][START_REF] Murabito | Heritability of age at natural menopause in the Framingham Heart Study[END_REF][START_REF] Morris | Familial concordance for age at natural menopause: results from the Breakthrough Generations Study[END_REF]). Here we used the regression between grandmother and mother, without controlling for possible shared environments (this is not crucial for the present study, as the main purpose here is to establish the validity of the mothers' age at menopause as a proxy for the daughters' expected age at menopause). Among the 97 women, 42 completed information about both their mother's and maternal grandmother's menopause. The mean age at menopause was 51.0 (range: 42-58) for the mothers and 50.3 (range: 39-60) for the grandmothers. The maternal grandmother's age at menopause had a significant effect on the mother's age at menopause (β = 0.275, SE = 0.11, P = 0.016, see figure 3). This implies a heritability of age at menopause from the maternal side of 27.5% and thus an overall heritability h² ~ 55% (SE = 0.22).

Discussion

Here, we show that expected age at menopause is significantly related to the facial attractiveness of young women. This result holds when variables potentially linked to immediate fertility (age, facial femininity, voice pitch and WHR), facial ageing and socioeconomic status are controlled for. Additionally, immediate fertility cues and facial ageing were not correlated with the expected age at menopause. Thus, the expected age at menopause seems to be linked to independent facial information, which is not deducible from previously documented components of facial attractiveness linked to immediate fertility or facial ageing.

Results classically found in the literature were replicated in this study, suggesting that this sample of young females is not different from those described elsewhere. First, our sample provided an heritability estimate for age at menopause of 0.55, consistent with previous values ranging from 0.42 to 0.72 [START_REF] Snieder | Genes Control the Cessation of a Woman's Reproductive Life: A Twin Study of Hysterectomy and Age at Menopause[END_REF][START_REF] De Bruin | The role of genetic factors in age at natural menopause[END_REF][START_REF] Van Asselt | Heritability of menopausal age in mothers and daughters[END_REF][START_REF] Murabito | Heritability of age at natural menopause in the Framingham Heart Study[END_REF][START_REF] Morris | Familial concordance for age at natural menopause: results from the Breakthrough Generations Study[END_REF]. Moreover, immediate fertility cues were positively correlated with facial attractiveness, as in previous studies: men in our sample preferred more feminine faces as in [START_REF] Cunningham | Their Ideas of Beauty Are, on the Whole, the Same as Ours -Consistency and Variability in the Cross-Cultural Perception of Female Physical Attractiveness[END_REF][START_REF] Perrett | Effects of sexual dimorphism on facial attractiveness[END_REF][START_REF] Rhodes | Does sexual dimorphism in human faces signal health?[END_REF][START_REF] Koehler | The relationship between sexual dimorphism in human faces and fluctuating asymmetry[END_REF][START_REF] Little | Facial attractiveness: evolutionary based research[END_REF], and facial attractiveness was positively correlated to pitch voice (F0), consistent with Collins and Missing 2003;[START_REF] Feinberg | The voice and face of woman: One ornament that signals quality?[END_REF][START_REF] Wheatley | Women's faces and voices are cues to reproductive potential in industrial and forager societies[END_REF][START_REF] Smith | Concordant Cues in Faces and Voices: Testing the Backup Signal Hypothesis[END_REF]. Also, attractive faces were rated younger than their true age, as in the study of [START_REF] Kwart | Age and beauty are in the eye of the beholder[END_REF]. However, two variables had significant effects on attractiveness, but in the opposite direction than expected. First, men in our sample had a preference for faces of women who have a higher (i.e. more masculine) WHR, which is in contradiction with the idea that face and body would signal one same quality. However, our results go in the same direction than results of [START_REF] Thornhill | The Body and Face of Woman: One Ornament that Signals Quality?[END_REF], who found a positive but insignificant correlation between facial attractiveness and WHR. Further studies are needed to investigate if face and WHR are signaling different aspects of female mate quality. Secondly, we found that socio-economic status of the women had a negative effect on facial attractiveness. We were expecting a positive correlation, as a higher socio-economic status is related to less stress during development, better nutrition, less unhealthy behavior, etc. [START_REF] Adler | Socioeconomic status and health: The challenge of the gradient[END_REF][START_REF] Kalick | Does Human Facial Attractiveness Honestly Advertise Health? Longitudinal Data on an Evolutionary Question[END_REF]. But here again, even if counterintuitive, our results are in the same direction than a previous study showing a negative correlation between facial attractiveness and SES for female [START_REF] Hume | Facial attractiveness signals different aspects of "quality" in women and men[END_REF]. We cannot speculate on this result, as the variance of SES in our female sample is very narrow (women were all students at the university). However, we can suggest that it deserves more study (a lot of research focused on the effect of attractiveness on SES, but less on the reverse relationship).

Spurious significant results may sometimes arise following model simplification [START_REF] Whittingham | Why do we still use stepwise modelling in ecology and behaviour?[END_REF]. This statistical bias is unlikely in this study, as all terms were kept and no model simplification was performed. However, it is still possible that a confounding variable, which remains to be identified, explains the link between expected age at menopause and facial attractiveness. Moreover, further study is needed to specify the facial cues used by men to detect residual reproductive value in young women (in this respect, figure S3 shows the differences in face shape between women of high and low expected age at menopause).

Here, we investigated the mediating role of immediate fertility cues and facial ageing, as residual reproductive value could potentially be linked to these two traits [START_REF] Wood | The evolution of menopause by antagonistic pleiotropy[END_REF][START_REF] Gold | The Timing of the Age at Which Natural Menopause Occurs[END_REF]). The next step would be to consider a large range of facial features such as facial symmetry, averageness, adiposity, skin colour, skin homogeneity, hair, etc., in an exploratory attempt to detect the facial cues linked to expected age at menopause. Also, age at menopause could additionally be assessed with more objective methods then self-report, for example by using the serum antimüllerian hormone concentration [START_REF] Disseldorp J Van, Faddy | Relationship of Serum Antimüllerian Hormone Concentration to Age at Menopause[END_REF], or through genetic analysis [START_REF] He | Genome-wide association studies identify loci associated with age at menarche and age at natural menopause[END_REF][START_REF] He | A large-scale candidate gene association study of age at menarche and age at natural menopause[END_REF]. Assuming that preferences reflect actual mate choice (which remains to be established for this trait), men's preference for women signalling a late age at menopause -and thus a longer temporal reproductive window -could suggests a current selection towards a later menopause, at least in Western societies. Interestingly, a secular trend of increased menopausal age has recently been observed [START_REF] Vannoord | Age at natural menopause in a populationbased screening cohort: The role of menarche, fecundity, and lifestyle factors[END_REF][START_REF] Rodstrom | Evidence for a secular trend in menopausal age: a population study of women in Gothenburg[END_REF][START_REF] Dratva | Is age at menopause increasing across Europe? Results on age at menopause and determinants from two population-based studies[END_REF]. This suggests that the social, familial and cultural conditions found today in Western societies are promoting current selection towards a later menopause in women. As a consequence, studies carried out in these societies are probably not adequate to empirically test evolutionary hypotheses on the ancestral selection on the emergence or maintenance of an extensive post-reproductive life-span in the human lineage [START_REF] Thouzeau | Emergence and maintenance of menopause in humans: A game theory model[END_REF]. Our results also suggest that sexual selection should be included in studies investigating the evolution of menopause in humans.

Somatic senescence, or ageing, is readily detectable in a face through wrinkles, skin texture and colouration, etc. (for a review, see [START_REF] Rhodes | Age Estimation of Faces: A Review[END_REF]. Those traits are specifically targeted by cosmetic changes when a younger appearance is desirable. Facial traits correlated with reproductive senescence have yet to be identified, and it remains to be determined whether they can be, or are already, manipulated to increase attractiveness by signalling a later menopause.

Figure S2.

Example of screenshot during the evaluation of women's age by the raters. The rater was instructed to estimate the women's age. The age appearing by default for each picture was randomized. For each woman, the difference between their real age and the mean age given by the raters was calculated. This variable represents the facial ageing of the women. Photograph reproduced with permission.

Figure S3

. Top: Shape differences of a low expected reproductive value face (average face of the 7 women with the earlier expected age at menopause in our sample) relative to the average face (average of all the women's face in our study), on a Thin-Plate Spline (TPS) grid.

Bottom: shape differences of a high expected reproductive value face (average face of the 7 women with the later expected age at menopause) relative to the average face. For both, the shape deformation has been magnified by 2. Figure created using R package Geomorph 3.0.3.

Figure 1 .

 1 Figure1. Example of screenshot during the evaluation of women's facial attractiveness by the raters. For each pair of women, which was randomly chosen, the rater was instructed to click on the photograph of the woman that he found the most attractive. Photographs reproduced with permission.

Figure 2 .

 2 Figure 2. Predicted probability to be chosen during the test of attractiveness, according to the difference in expected age at menopause (unstandardized values) between the focal and the non-focal women, controlling for all the other variables of model 2. Circle areas are proportional to the number of choices made by the raters: chosen (1) or not (0) as the most attractive woman. 95% confidence interval is represented by the dotted lines.

Figure 3 .Figure S1 .

 3S1 Figure 3. Relationship between age at menopause of mothers and maternal grandmothers.Plain line represents the linear regression with all points (β = 0.275, SE = 0.11, P = 0.016, corresponding to a heritability of age at menopause h 2 = 0.55). 95% confidence interval in grey.

Acknowledgements

The authors would like to thank Mona Joly for assisting in data collection, as well as the City Hall of Montpellier, Luc Gomel and all the staff from the Serre Amazonienne for providing places to recruit raters. We also thank two anonymous reviewers for their valuable comments on the manuscript. This study was supported by the Agence Nationale pour la Recherche 'HUMANWAY' project (ANR-12-BSV7-0008-01). Funding from the ANR through the Labex IAST to J.B. is also gratefully acknowledged. This is contribution 2017-283 of the Institut des Sciences de l'Evolution de Montpellier (UMR CNRS 5554).

Tables and figures captions 510

Table 1. Effects of the different variables on the probability to be chosen during the test of 511 attractiveness (male raters had to choose the woman found to be the most attractive between 512 two facial photographs). N = 3439 observations (119 male raters and 68 female faces).


Laurent Duchemin

Christophe Eloy

Eric Badel

Bruno Moulia

Tree crowns grow into self-similar shapes controlled by gravity and light sensing

Plants have developed di erent tropisms: in particular, they re-orient the growth of their branches towards light (phototropism) or upwards (gravitropism). How these tropisms a ect the shape of a tree crown remains unanswered. We address this question by developing a propagating front model of tree growth. This model being length-free, it leads to self-similar solutions, independent of the initial conditions, at long time. Varying the intensities of each tropism, di erent self-similar shapes emerge, including singular ones. Interestingly, these shapes bear similarities with existing tree species. It is concluded that the core of speci c crown shapes in trees relies on the balance between tropisms.

I. INTRODUCTION

In many tree species the outer shape of the crown is a criterion for species identi cation, especially when trees have grown in isolation. The Mediterranean cypress has a distinctive elongated shape, while birch and oaks show more spherical shapes, and the Norway spruce has a conical prole. The growth process that leads to these di erent silhouettes is not well understood. It is associated with genetically xed species-dependent factors like the branching angles, the phyllotaxy, and more generally the developmental characteristics of the branching architecture known as the "architectural model" of the species [1]. But crown shapes are also affected by the environmental growth conditions (in a geneticdependent manner): trees grow di erently if they are isolated or in a forest [2], if subject to climatic stresses like wind or snow cover [3]. So far, the few studies about crown shaping have investigated the hypothesis that it may result from an interplay between a genetically-de ned architectural development and branch growth or shedding linked to light competition [START_REF] Měch | Proceedings of SIGGRAPH[END_REF].

However tropisms in the outer growing branches can be viewed as another possible candidate for the control of crown shape. Tropisms are de ned as the re-orientation of a growing branch following a vectorial cue from the environment. Two tropisms are shared by most plants: phototropism is the process that leads to growth in the direction of light, while gravitropism is driven by the direction of earth gravity. Genetic di erences in tropic sensitivities have been reported [START_REF] Sierra-De Grado | [END_REF]6]. The study of these tropisms is still an active eld today, ranging from shoot scale [6][7][8][9], to the molecular networks that regulate them [10]. Yet, a question remains largely open: how do these tropisms and their balance in uence the shape of a tree crown?

In this work we investigate the simple hypothesis that the two tropisms of the peripheral primary shoots at the tip of the branches can be major players in driving the shaping of the crown, accounting for the interplay between genetic control and responses to the local environment. To do so we developed a simple model of the crown growth. In the literature, popular growth models generally consider branches and leaves as the elementary bricks of a numerical simulation [11][12][13][14][15]. However these models usually involve a large number of empirical parameters, and thus remain too complex to study the speci c roles of phototropism and gravitropism in selecting the tree shapes. Recently, Beyer [16] proposed a simpler and more parsimonious model that considers the whole tree as a continuous medium. This approach seems very promising, but has not yet been fully exploited to address the combined role of phototropism and gravitropism. In this Letter, we propose a new idea: considering the growth of a tree crown as a continuous front propagation process. Despite some similarities with other front propagation processes, such as crystal growth [17], premixed combustion [START_REF] Clavin | Combustion Waves and Fronts in Flows: Flames, Shocks, Detonations, Ablation Fronts and Explosion of Stars[END_REF], or Sa man-Taylor ngering [START_REF] Pelcé | Dynamics of curved fronts[END_REF], tree growth has some distinctive features.

In particular, since growth involves photosynthesis and is controlled through light sensing (photomorphogenesis), our model includes both photosensitivity (the foliage receiving more light will grow faster), phototropism (growth is preferentially oriented in the direction of light), and also gravitropism (growth is preferentially oriented against gravity). As we shall see below, the proposed model requires only two parameters to describe the shape of the crown (irrespective of its size): the intensities of phototropic and gravitropic growing responses. With this model, we intend to address the following questions: What is the family of shapes generated when these two parameters vary? Do these shapes converge towards steady solutions? How these steady solutions can be described analytically? And nally do these solutions t the major classes of crown shapes observed in undisturbed isolated trees in nature?

II. CROWN GROWTH MODEL

We rst assume that the average sunlight is uniformly distributed in the upper hemisphere, but we shall come back to this strong assumption in the discussion section. When the crown envelope is convex and axisymmetric, the amount of light received daily at each point of the front is thus assumed to be proportional to ψ , the angle between the local tangent and the horizontal (Fig. 1). At each point on the front, n rep-resents the outward normal unit vector, is the unit vector that represents the average direction of light, v is the upward vertical unit vector, and κ is the in-plane curvature, positive when the surface is locally concave (Fig. 1). According to this nomenclature, the velocity of the front is written as

U = ψ n + α v + α p n + α v + α p + γκn, (1)
where α and α p are the intensities of the crown gravitropic and phototropic responses respectively, and γ is similar to a "surface tension", but acting on the growth velocity.

n v ϕ ϕ ψ h t θ R U n δt rC δt π/2 -θ ϕ β Figure 1
. Sketch of a growing tree crown. h and v are unit horizontal and vertical vectors respectively, n and t are the unit vectors respectively normal and tangent to the front and the unit vector points towards the mean direction of light, which is the rst bisector of (h, t). The angle ψ = φ + π /2 represents the local amount of sunlight intercepted for this axisymmetric shape. The inset shows a zoom around the front to highlight the conditions for self-similarity of the growing shape.

Our model is based on two basic assumptions: [A1] the crown foliage is dense enough all across the crown so that it may be described as a continuous medium with a continuous boundary line (i.e., a growth front); [A2] all the growing apices sitting on the front have similar dependency on light and gravity. Assumption [A1] requires that the development of the tree crown is su cient and the crown displays a steady characteristic shape, a condition that is usually observed in most species as soon as they have reached the "architectural unit" (AU) stage, i.e. that they have di erentiated all their branch types [START_REF] Edelin | [END_REF]. [A2] requires that "rogue branches" do not occur. It is well known that traumatic reiteration resulting from branch accidental breakage or from large herbivory damages produces re-juvenilised branches that behave differently than regular branches (e.g., suckers). These traumatic reiterations should be distinguished from regular reiterations involved in the process of crown morphogenesis (a process known as branch metamorphosis [1,[START_REF] Edelin | [END_REF]), as these reiterated branches still complie with a regular behaviour and the crown remains dense.

The rst term on the right-hand-side of Eq. ( 1) includes two important biological mechanisms. First, the front velocity is assumed to be proportional to ψ , the light intercepted, to account for photosensitive growth. Second, the gravitropism and the phototropism are modelled through a re-orientation of the growth in a direction computed as the weighted average of the three unit vectors n, v, and . This modelling of tropisms is in accordance with the experimental observations made at the level of a single shoot, and with the ArC model of shoot gravi-photo-tropism, that was assessed against experiments [21]. Additionally, it has been demonstrated experimentally that the sensitivity to the gravity is not dependent on the angle made by the shoot versus gravity. The mechanism subtending this behaviour lies in the functioning of specialized cells called statocytes (see [6,7,22,23] for more explanations on the physics and biology behind this mechanism). For phototropism, the equation and the de nition of phototropic sensitivity are mainly similar (see [21] and references therein).

The last term in Eq. ( 1) stabilises the front dynamics and smoothes out the front shape by damping the velocity uctuations, and hence has a similar role as surface tension in Sa man-Taylor dynamics for example. There is obviously no surface tension acting at the boundary of the tree crown. However recent ecophysiological investigations have revealed that two biological mechanisms are acting as to reduce the di erences in growth velocities of neighbouring shoots, and hence resulting in a attening tendency of the canopy boundary, which is qualitatively similar to a surface tension. The rst mechanism is due to lateral sensing of the spectral signature of the light re ected by neighbouring plants, through the phytochrome pigment [24]. This sensing results in a photomorphogenetic synchronisation of the growth in length of the neighbouring stems, keeping the top of the canopy at. The second mechanism involves the sensing of wind-induced strains [25]. Whenever a shoot overreach its neighbours, it is not sheltered by the canopy anymore and its growth speed is reduced until the front is at again. It has been shown nicely in herbaceous populations of shoots that these two mechanisms are responsible for the attening tendency of canopy tops, as dramatically illustrated by the atness of the top surface of crops such as wheat or corn [26]. On a practical point of view, when simulating the growth of a crown, we rescale its shape at each time-step to keep its axi-symmetric volume V constant. As explained in the next section, this is equivalent to a coe cient γ growing in time in Eq. (1). Although this rescaling may seem unphysiological, we perform it for two reasons: it ensures that the code is stable at long times, and it allows to reach selfsimilarity (when surface tension is rescaled, Eq. (1) has no length scale and self-similar solutions emerge).

Finally it may be noted that in our simpli ed model, the curvature κ is calculated in a vertical plane, and not as the three-dimensional mean curvature of the axisymmetric interface. The latter re nement would be possible, but does not a ect the main features of the solutions.

III. NUMERICAL SOLUTIONS

Equation ( 1) is solved numerically by placing N = 200 marker points on the front. The points are then advected along the normal to the front. Furthermore, at each time-step, the coordinates of the front are rescaled by the volume to the power 1/3, such that the crown volume remains constant (V = 1). This rescaling makes the problem scale-invariant and allows the "surface tension" to remain of the same order during the growth process. As a consequence, the scaling of the true curvature in the physical plane increases like V 1/3 . Since the "surface tension" coe cient in the numerical space γ * is kept constant, it means that :

γ = γ * V 1/3 . (2)
Since V = 1, we can drop the star and use γ in place of the constant γ * . The surface points are advected in time using an explicit second-order scheme [27], with a time-step δt ∼ δs 2 /γ , where δs is the typical distance between two successive marker points, chosen to ensure the stability of this explicit scheme. The markers are redistributed on the front at each time-step, in order to conserve a regular spacing. Note that the "surface tension" becomes signi cant when γκ ∼ O(1), as seen in Eq. (1).

Figure 2 shows two examples of such a computation for two distinct initial conditions. In this example, there is no phototropism and no gravitropism (α = α p = 0) and surface tension is set to γ = 0.01, small compared to the global length scale: γ 1. We rst observe that, although initial conditions significantly di er, the successive shapes converge towards a unique self-similar shape (Fig. 2a). This self-similarity was expected because the problem has been made scale-invariant.

To further test this idea, we performed extensive numerical simulations with di erent parameter sets (α , α p ), at low surface tension (γ = 0.01). The results are presented in gure 3: green shapes correspond to computations that have converged towards self-similar shapes, meaning that a typical distance between two successive rescaled shapes (Fréchet distance) has reached 10 -6 , whereas grey shapes correspond to non-convergent computations (in these cases found in the lower-left part of the diagram , the shape attens out indefinitely and convergence is never reached). The grey curve which segregates converged from non-converged shapes will be discussed in section V.Figure 3 displays a large diversity of shapes as a function of the parameter sets (α , α p ), showing that the (genetic) variation in the sensitivities to phototropism and gravitropism can indeed produce di erent crown shapes in our model. Moreover, we performed numerical simulations with decreasing surface tension γ (Fig. 4). The results showed that all the nal self-similar shapes converge towards a universal shape in the limit γ → 0. In the next section, we focus on this limit and show that this universal shape can be derived analytically, providing insight into the control of the steadystate crown shape, when it exists, by the parameter sets (α , α p ).

IV. ANALYTICAL SOLUTIONS

In the limit of vanishing surface tension (i.e. γ → 0), Eq. ( 1) becomes scale invariant. We look for a self-similar solution of this equation by imposing that the front is described by the homothetic surface:

R(θ, ϕ, t) = r (θ )C(t), (3)
in spherical coordinates (Fig. 1).

Between two successive times, t and t + δt, the front has been increased by r (θ )C (t)δt in the radial direction. This condition can be expressed as

r (θ )C (t)δt cos β = U • n δt = ψ B(φ) δt, (4)
where

ψ = φ + π /2, β = π /2 -θ -φ and B(φ) = n + α v + α p n + α v + α p • n. (5)
To solve this problem, we rst note that cones are trivial self-similar solutions of Eq. ( 4), with φ the constant angle between the cone surface and the vertical direction (Fig. 1). Moreover, Eq. ( 4) implies that C (t) does not depend on time. The cones having the same velocity C (that we choose equal to 1 without loss of generality) intersect the vertical axis at the elevation

z c (φ) = r (0) = (φ + π /2)B(φ)/sin φ. (6)
Figure 4 shows the cones described by Eq. ( 6) for two different sets of parameters. The universal shape we found numerically in the limit of vanishing surface tension appears as the inside envelope of all those cones. Interestingly, this inside envelope can be obtained analytically.

For convenience, we switch to cartesian coordinates. The cones of Eq. ( 6) intersect the plane = 0 along the straight line described by the equation: z = z c (φ)x/tan φ. Consequently, two cones corresponding to angles φ and φ + δφ will intersect at the point

x(φ) = -z c (φ) sin 2 φ, (7)
z(φ) = z c (φ) + z c (φ) sin φ cos φ. (8)
When -π /2 ≤ φ ≤ π /2, Eqs. (7-8) describe the cone envelope as a parametric curve (red curves in Fig. 4). Because this parametric curve is everywhere tangent to a cone with the same time evolution C(t), it is also a self-similar solution of Eq. ( 1) for γ = 0. This approach is similar to the Wul construction used in the physics of crystal growth [START_REF] Einstein | Handbook of crystal growth: Fundamentals[END_REF][START_REF] Pimpinelli | Physics of crystal growth[END_REF].

Another interesting feature of these self-similar shapes arises from following the trajectories of surface points backward in time. These trajectories can be thought of as traces of the underlying structure, i.e. the tree branches. Figure 6 shows these trajectories for two cases: α = α p = 0 and α = α p = 2. These trajectories are numerically calculated by advecting backward in time 31 points along the vector -U . At each time step, the points -lying on a smaller selfsimilar shape-are associated to an updated angle ψ using cubic splines. For α = α p = 0, there is a region without branches close to the bottom. This region is a consequence of the at bottom, where ψ = 0 and the front velocity is zero. Note that the curvature of the "branches" is only due to the variation of ψ and along these trajectories, and not to the bending under self-weight of the branches or any other global re-orientating mechanisms such as those related to reaction wood [START_REF] Moulia | [END_REF], which are not taken into account in our model.

We will discuss below why numerical simulations converge towards the parametric solution given by Eqs. (7)(8). Meanwhile, we describe how the respective intensities of gravitropic and phototropic responses modify the crown propagation and the nal shape.

V. COMPARISON BETWEEN NUMERICAL AND ANALYTICAL SOLUTIONS

We note that these attracting shapes do not exist for any choice of the parameters α and α p , as already observed in gure 3. Figure 5a shows with grey levels how fast the numerical solution converges towards a self-similar solution. It also shows that in the lower part of this diagram, there is no convergence. Equation (6) indicates the cases when selfsimilar shapes exist. The convergence occurs when the cone with a vertical angle (φ = 0) exists (i.e. x(0) > 0 in Eq. 7). This condition simpli es to: 6) is the expected self-similar shape. When γ → 0, the time-evolving shapes converge towards this self-similar shape.

α p > - √ 2 (9) γ γ γ γ z c (a) ϕ (b)
Moreover, the self-similar shapes cease to exist when their top (given by z c (φ) = 0, φ ∈ [0, π /2]) reaches their bottom (z c (φ) = 0, φ ∈ [-π /2, 0]). This condition can be computed numerically and corresponds to the light green curve in gure 5a. A good approximation of this curve for α ∈ [-2, √ 2 -1] can be obtained by writing that the selfsimilar shapes exist roughly when the two cones with horizontal angles (φ = ±π /2) are one above the other (i.e. z c (π /2) > z c (-π /2) in Eq. 6), a condition that simpli es to:

α p + α > -1 (10)
These existence conditions correspond to the light green curve in the phase diagram (Fig. 5), below which there is no analytical self-similar solution. The curve is also reproduced in gure 3 and is perfectly consistent with the convergence region observed for the time-evolving shapes.

When they exist, self-similar solutions exhibit di erent characteristics that correspond to di erent hatchings in Fig. 5a. First, we note that the self-similar solutions always feature a pointed top. It happens for z c (π /2) > 0, with z c given by Eq. ( 6). This condition is equivalent to α p +α > -1; which is guaranteed by Eq. (10).

Second, the self-similar shape can exhibit a pointed bottom when z c (-π /2) > 0, or equivalently when x(-π /2) < 0. This con guration is equivalent to α > 1 (dark green hatching on the right of Fig. 5a). Examples of such shapes are displayed in Figs. 5c,i,j. Note that, contrary to at-bottomed shapes (Figs. 5b,e,g, for instance), these shapes have a non zero growth velocity at their pointed bottom.

Third, the bottom of the self-similar shape can present a cusp when x (-π /2) ≤ 0; which is equivalent to α p (2α 2 -2α + α 2 p ) ≤ 0 (red hatching in Fig. 5a). This cusp can clearly be seen in the self-similar solutions displayed in Figs 5i,j where the value of x at the bottom (for φ = -π /2) initially decreases. The same cusp is also perceptible in Figs 5d,g,h, although it is less pronounced.

Finally, the self-similar shape described by the parametric curve given in Eqs. (7)(8) can exhibit a loop, as seen in Figs. 5h,i (it is also perceptible in Fig. 5f, but less visible). This loop corresponds to the zone with the yellow hatching in Fig. 5a, whose boundary has been determined numerically.

As seen in Figs. 4 and 5b-i, the numerical self-similar solutions obtained at long time always converge, in the limit of vanishing surface tension, towards the inward envelope of the cone solutions given by the parametric curve of Eqs. (7)(8). However, when surface tension is zero, any shape constructed as an assembly of pieces of cones could be a valid self-similar solution. These solutions do not occur in the numerical simulations for two reasons: rst, any connection between two pieces of cones tend to "recess" towards the inward envelope when the surface tension is non zero. Second, the front propagation is described as the advection of markers in the direction normal to the front, whereas a special numerical treatment of the corners would be necessary to maintain a solution made of several assembled cones. The grey level is proportional to the number of time-steps necessary for convergence towards a self-similar shape (black corresponds to no convergence). The hatchings correspond to di erent features of the self-similar shapes: light green for no self-similar solutions, red when a cusp is present at the bottom, yellow for a loop, and dark green for a pointed bottom (more details in the main text). The white square corresponds to gure 3. (b-j) Characteristic shapes found numerically for γ = 0.01 (green curve and ll) and analytically in the limit of vanishing surface tension with Eqs. (7-8) (red curve). The thin black lines correspond to the vertical symmetry axis (Oz) and to the horizontal plane z = 0. Shapes (i) and (j) have been rescaled by a factor of 2 in both directions, whereas all others shapes have the same volume.

VI. COMPARISON WITH OBSERVED CROWN SHAPES

Care should be taken when comparing the predictions of the model with real tree crowns observed in nature. Indeed the impinging of light rays on the crown should not be intercepted by a neighbouring object (wall, rock, other trees). Moreover tree shoot should not be able to sense the re ected light from neighbouring crowns even in the absence of direct shading of the incident light. The few studies about the distance of neighbouring sensing is limited to approximately 5 m [31]. Finally no external action such as pruning, grazing or wind-break should have occurred (at least at a scale that could a ect the growth and shaping of the crown).

Another important aspect is that the external growing shoots should have a similar physiological and developmental status. Indeed some reiterations (especially traumatic ones) can result in obvious disruptions of the regularity of the crown shape. Therefore, only two developmental stages were retained: (i) young trees approaching the "architectural unit stage" and (ii) mature trees having developed their scaffold limbs through the process of developmental metamorphosis [1,[START_REF] Edelin | [END_REF].

We used two sources of tree crown pictures. The rst source is an arboretum (E. Badel's park, sitting in southern France in a sub-Mediterranean climate at 44.904°N, 4.830°E) where no pruning or grazing has ever been achieved and trees were planted at large distances from each other. The second source results from a search on botanical websites and Wikimedia Commons (https://commons.wikimedia.org), retaining only photographs for which the species and the location were known and the criteria concerning the absence of signi cant neighbours or traumatic reiteration could be assessed.

For each picture, the best t is found using the following method. First, the right half part of the tree crown is described as a series of hand-picked points (typically about 20 points). Then, we assume that the shape is axisymmetric and we rescale and shift this shape such that its volume is equal to unity and its center of mass is at the origin. After that, we look for the parametric curve described by Eqs. (7-8) that minimises the distance with the real tree crown, where the distance d is de ned as

d 2 = 1 π ∫ π 0 (R model (θ ) -R crown (θ )) 2 dθ, (11)
with R model and R crown the polar representations of the model and tree crown respectively. This procedure allows to nd the best parameters α and α p . The same procedure is repeated with the additional constraint α p > 0. This tting method has been applied to 36 di erent tree crown pictures (see Supplementary Materials), out of which four representative cases have been extracted ( gure 7 For some trees, the shape of the crown could not be well tted with positive values of α and α p . However, allowing a negative value for α ( gure 7b) or α p ( gure 7c,d), leads to a good t. Negative values for α p would mean a negative shoot phototropism, a condition that is usually not observed in studies about shoot phototropism. We analyse in the discussion section the biological insights that such a situation may provide. Figure 7(c) and (d) provide two sets of parameters: the rst one (red pro le) for which we allowed any values for α p , the second one (black pro le) for which α p is constrained to be positive. These two comparisons clearly show that in some cases, like (c), the crown shape can be well tted using several parameter sets (possibly even a continuum set of parameters), whereas others show a poor comparison with the real shape when α p is constrained to be positive. On this matter, the diagonal described by equation (10) and seen in gure 5 plays a particular role: close to this line, shapes look similar, allowing for a family of solutions.

VII. DISCUSSION

Many tree species display characteristic and heritable crown shapes, at least at some stages and when trees have grown in isolation. Yet crown shape is remarkably variable with the environment (what biologists call "platicity"): trees grow di erently if they are isolated in full sunlight or in a forest [2], if they are submitted to wind or protected from it [3]... Reconciling these two features has remained a challenge. Indeed there is a large number of biological mechanisms that can a ect the intensity and orientation of shoot growth, e.g. phyllotaxis, apical dominance [1,[START_REF] Edelin | [END_REF], response of buds break to light [24], changes in growth directions at di erent developmental stages, tropisms, internal correlations between organs [1], as well as repulsive e ects between branches mediated though phytochromes [24], biomechanical processes at the branch level, such as the bending of the branches under self-weight, its xing through secondary growth and global re-orientating mechanisms related to reaction wood [START_REF] Moulia | [END_REF]. Moreover, the development of the crown is modifying the environment of the growing apices, creating a feedback of the actual crown shape over its development (shading, wind sheltering...). All the challenge is then to identify the leading physical and physiological mechanisms which can explain the interaction between genetic control and sensitivity to the environment in the development of crown shape. Two standpoints can be adopted.

One is to concentrate on the development of the branched architecture through bud breaks, shoot growth, and its remodelling through environmental factors interacting with the crown size and spread such as light competition or wind e ects [START_REF] Měch | Proceedings of SIGGRAPH[END_REF]32]. The second stand-point is to concentrate on what is occurring at the boundary of the crown, and more speci cally on the apical buds that are sitting on its vicinity, and which are responsible for the increment in crown spread through the annual ush of shoot growth. This standpoint can be visualised when considering for example a coniferous tree in the springtime, on which the crown of the last year has kept its dark-green needles and the new sprouts undergoing primary growth (elongation) shows up with a much lighter green colour. The descriptive parameters of this crown increment are the distribution of the buds, the amount of their elongation, and their direction. And the interaction between the genetic speci cation and the local environment at the boundary of the crown is speci ed through the physiological processes that may in uence these descriptive parameters. This second standpoint is the one that has been investigated in this work. Three di erent processes were thus retained: the photomorphogenetic photosensitivity of shoot elongation (referred to as photosensitivity) and the two major tropisms known to a ect the direction of the primary growth: photo-and gravi-tropism [6,21]. The initial orientation of the bud was considered to have negligible e ect as growing shoot can reach their photo-gravitpropic set-point angle irrespective of the original orientation in an interval of time that is very short compared to the morphogenesis of the crown shape (typically less than a day) [6,21].

In this article, we model this process of incremental primary growth at the boundary of the crown with a front propagation equation taking into account the e ects of photosensitivity, phototropism and gravitropism, and the feedback of the actual crown shape through its shading e ect of the growing apices, as a function of their position. We also included a synchronisation of the growth velocity among neighbouring shoots [24][25][26] resulting in a local attening tendency of the canopy boundary, which is qualitatively similar to a surface tension. Through numerical simulations, we show that this front equation generally yields self-similar solutions at long time, independently of the initial shape. In the limit of vanishing surface tension, we show that these selfsimilar solutions can be described analytically by parametric curves, which, for particular values of the phototropic and gravitropic responses, can exhibit singularities.

The outputs of the model were then compared to real tree crowns. As the dynamics of the crown growth is at the core of the model, a direct comparison between the proposed front velocity, given by Eq. ( 1), and empirical observation of growth directions and amount would have been the best assessment. However, no such data was available, and their collection would require long term experiments. The assessment of the predictions of the model was thus achieved by comparing its self-similar solutions with real tree crowns at the special stages at which the crown displays a typical and stationary shape (namely the "architectural unit stage" and mature stages [1,[START_REF] Edelin | [END_REF]). This was achieved on 36 crown images of di erent tree species sampling many taxonomic clades from di erent climates. Using the distance d de ned by Eq. ( 11) measuring how the shapes extracted from the pictures di er from the analytical parametric curves, we could nd, for each picture, the set of phototropic and gravitropic intensities that best t the tree crown. It appears that in 97% of the cases, we nd a distance d < 0.05, which we consider as good (see Supplementary Materials).

These results however raise two major questions: one about the values of the photo-and gravisensitivity parameters; and the second one about the signi cance of the two important assumptions of the model (aside from the choice of the three mechanisms): the linear photosensitivity curve and the uniform lighting. Regarding the photo-and gravisensitivity parameters, positive values of α p and α allow for a satisfactory comparison with real crown shapes in only 33% of the cases. For 75% of the cases, a reasonable t can be obtained with the constraint α p ≥ 0. However negative values of α p are required to t the remaining 25% ( gure 7d and Supplementary Materials). Negative shoot phototropism has never been documented in wild-type plants living in natural conditions and can be considered non-physiological. By the same token, negative shoot gravitropism is restrained to rare mutants, and is unlikely to be involved in the studied trees. However, these negative values may re ect the in uence of another tropism, namely auto-tropism [START_REF] Moulia | [END_REF]. Auto-tropism is the tendency for growing plant shoots to reduce their curvature, through a proprioceptive sensing [6,33]. It was found that the proprioceptive auto-tropism is a major player of the orientation of shoot growth [6,21]. Auto-tropism was not included in our model and negative phototropism could compensate for that. But it could also compensate for other processes acting at the whole branch level (e.g, branch biomechanical reorientation or shedding) that were not taken into account in the model. Some explicit representation of the growing shoot or of the whole branch would be required to asses if these additional mechanisms may account for the crown shapes that cannot be tted with our model. Therefore a new research avenue emerging from the present work is to extend the model to account for the average orientation of branches along the crown boundary, the e ect of, for example, autotropism, or self-weight [START_REF] Moulia | [END_REF] as well as some mechanisms of spatial competition [34] or of branch repulsion [24]. In this work we have assessed whether following the trajectories of surface points backward in time gives an idea of the underlying branches structure. But e ects such as auto-tropism or bending under self-weight were not taken into account in our model and the "inferred branch patterns" were clearly not very realistic.

The two other hypotheses of the model: i) linear photosensitivity of growth (the front velocity is assumed to be proportional to the average light intercepted) and ii) uniform lighting (light, on average, is assumed to be uniformly distributed in the upper hemisphere) remain also to be discussed. Indeed, whereas the basis for the phototropic and gravitropic response lays on a wealth of published reports, as shown in the introduction section, the situation is less favourable here. Both hypotheses can be advocated to be a good starting point that allows to keep the model tractable and understandable (parsimony argument); yet they should not be ad hoc. The rst hypothesis is grounded in the fact that plants have light sensors called phytochromes that allow them to be sensitive to the ratio between the Red light and the Far Red light (R:FR) [24]. FR light is almost fully re ected by foliage elements, whereas R light is largely absorbed. Then, for a nite front area A, the R:FR ratio can be well approximated by the ratio of visible sky and of visible surrounding foliage. When A → 0, assuming that the density of surrounding foliage is uniform, this ratio remains proportional to the angle of visible sky, and this corresponds to our model. In general it has been found that there is a linear relationship between the elongation rate of shoots and the R:FR ratio [24], which is usually negative for many open-habitat species (leading to a shade avoiding) and positive for some tree species (a behavior known as shade tolerance) [35]. The second hypothesis seems much less straightforward. It is only realistic for fully overcast skies on cloudy days. On sunny days, the di use light from the sky displays a fairly uniform radiance, but obviously the direct sunlight is highly non-uniform, with a bias towards the South (resp. North) in the Northern (resp. Southern) hemisphere [36]. And including a direct sunlight component in our model would clearly break the axi-symmetry of the crown shape due to a faster growth in the sun-facing side. On the contrary, tree crowns in Nature do not seem to be biased towards the sun [37] (with the exception of Araucaria columnaris [38]). Is this a serious drawback for our model? Not necessarily. Indeed the phytochrome sensing of the R:FR ratio of the light is compensated for light intensity [24]. Besides the elongation growth is known to be impaired by water stress, that is related to the direct sunlight, so that higher growth may actually happen on overcast cloudy or rainy days, and hence under uniform lighting [39,40]. Yet, to our knowledge, no explanation has ever been published to this puzzling question of crown axi-symmetry in isolated trees despite the asymmetry of many environmental factors (sunlight of course, but also evaporation, temperature...). But now that we have shown that the two previous assumptions coupled with phototropism and gravitropism led to realistically looking tree crown shapes, it is clear that investigating more thoroughly light-compensated sensing mechanisms and their relation to shoot growth and to crown shaping in trees is likely to be very rewarding. And including more mechanistic light responses in our model can be a tool to do so. Despite this set of open points, our work clearly shows that the hypothesis that the observed crown shapes could be the result of a front propagation process driven by the lightand gravity-sensing at the crown boundary is a valid starting point. Moreover, it provides an insight on how to reconcile the genetic speci city and the environmentally-driven variability (plasticity) of crown shapes. Indeed if a tree has enough time to reach a stationary shape, then there should be no more in uence of the initial state and/or the pathway to this form. This stationary shape emerges from the recurrent e ects of elongation photosensitivity and photo-and gravitropic sensitivities, and the feedback on these three processes of the interaction between the lighting environment and the actual crown shape. But the values of the three sensitivities are genetically-xed. This may explain how a process driven by three responses to the environment (and hence allowing plasticity) yields a steady-state form that is heritable and probably genetically controlled. The situation is similar with what has been found for the shape of single shoots by [6,21]. But for this to be achieved there should be no other player in the environment*crown shape interaction, a case only met when trees are grown in isolation and other factors -e.g. shading by neighbours, wind, frost-are negligible. It is likely not to be the whole story though. Indeed there are indications that the steady shapes of the same tree at the (young) "architectural unit stage" and at the mature stage are di erent [1,[START_REF] Edelin | [END_REF] (although this does not rely on observations of the same tree at the two stages, grown in isolation). One may think that the shape control through photomorphogenetic and tropic mechanisms captured in our model only stabilises each of the two forms, whereas the shift between the two stages could be driven by changes in the branching development when trees developed their sca old limbs through the process of developmental metamorphosis.

To provide di erent attractors, though, non linear behaviour is required in our model; but there are many candidates for that (eg. interaction between di erent light sensors, phytochromes, cryptochromes [24]. . . ). However, before that, our work also points at the need for data following the development of isolated trees from di erent species during decades. Only this type of data may allow to assess whether self-similarity is reached. Additionally the e ect of minor interventions, such as partial pruning, could be studied experimentally by looking at the shape dynamics and relaxation towards a self-similar shape. Such detailed data are required to assess our model more accurately, as well as to assess the alternative models. It would be also very informative to have a full 3D shape (or at least view at right angles) to assess the axisymmetry of the crown and its robustness to environmental gradients across the crown. This should be a priority now. But hopefully new tools such as terrestrial LIDARs or high-speed photogrametric devices are now available for that [41]. More detailed experimental studies of the diversity of photomorphogenetic and tropic responses for di erent species or under di erent environmental conditions would also be extremely useful (so that the sensitivities in our model could be measured independently rather then tted ).

From a broader perspective, our approach also opens the way to interdisciplinary research to whole-plant morphogenesis. A rst avenue is to change scale in our understanding of plant morphogenesis, form the scale of the meristem or single shoot [6,33] to that of the whole tree. Advances on these aspects would help to further elucidate the mechanisms driving the natural variability of crown shapes [6,9] and how plants read their own shape [33]. Finally, in a di erent context, the approach proposed in the present article could be useful to study the growth of other biological systems that display crown growth, such as stony corals [42], for which local growth velocities and direction depend both on internal regulation processes (nutrient exchange, for instance) and external cues ( ow velocity, nutrient capture, etc. . . ).

Figure 3 .

 3 Figure 3. A set of crown shapes obtained from a numerical computation, varying α and α p , and with γ = 0.01. Green shapes have converged towards a self-similar-one, whereas grey shapes have not converged. For each shape, the parameter values correspond to the grid intersection.

Figure 4 .

 4 Figure 4. Attracting shapes obtained for di erent values of γ : γ = 0.04, 0.02, 0.01, 0.005 and 0.0025, with α = α p = 0 (a) and α = 0.8, α p = -1 (b). The envelope of the grey lines given by equation (6) is the expected self-similar shape. When γ → 0, the time-evolving shapes converge towards this self-similar shape.

Figure 5 .

 5 Figure 5. (a) Phase diagram in the α -α p plane.The grey level is proportional to the number of time-steps necessary for convergence towards a self-similar shape (black corresponds to no convergence). The hatchings correspond to di erent features of the self-similar shapes: light green for no self-similar solutions, red when a cusp is present at the bottom, yellow for a loop, and dark green for a pointed bottom (more details in the main text). The white square corresponds to gure 3. (b-j) Characteristic shapes found numerically for γ = 0.01 (green curve and ll) and analytically in the limit of vanishing surface tension with Eqs. (7-8) (red curve). The thin black lines correspond to the vertical symmetry axis (Oz) and to the horizontal plane z = 0. Shapes (i) and (j) have been rescaled by a factor of 2 in both directions, whereas all others shapes have the same volume.

Figure 6 .

 6 Figure 6. Trajectories of points along the front, for the cases α = α p = 0 (left) and α = α p = 2 (right). Starting from the self-similar pro le corresponding to t = 1, 31 points are advected backward in time using Eq. (1) until t = 0.1.

  Comment citer ce document : Duchemin, L., Eloy, C., Badel, E., Moulia, B. (2018). Tree crowns grow into self-similar shapes controlled by gravity and light sensing. Journal of The Royal Society Interface, 15 (142), 10 p. , DOI : 10.1098/rsif.2017.0976

Figure 7 .

 7 Figure 7. Comparison between real tree crowns and self-similar shapes of the model: (a): Betula pubescens, α = 0.927, α p = 0.267, d = 0.022; (b): Quercus castaneifolia, α = -1.00, α p = 0.353, d = 0.019; (c): Enterolobium cyclocarpum, α = -0.490, α p = -0.488, d = 0.023 (black curve: α = -1.00, α p = 0.034, d = 0.020), (d): Thuja occidentalis, α = 0.379, α p = -1.23, d = 0.018 (black curve: α = 5.00, α p = 2.69, d = 0.079). The red curve represents the best t and the black curve the best t with α p > 0. See Supplementary materials for comparisons with other tree species.

  Comment citer ce document : Duchemin, L., Eloy, C., Badel, E., Moulia, B. (2018). Tree crowns grow into self-similar shapes controlled by gravity and light sensing. Journal of The Royal Society Interface, 15 (142), 10 p. , DOI : 10.1098/rsif.2017.0976

  Comment citer ce document : Duchemin, L., Eloy, C., Badel, E., Moulia, B. (2018). Tree crowns grow into self-similar shapes controlled by gravity and light sensing. Journal of The Royal Society Interface, 15 (142), 10 p. , DOI : 10.1098/rsif.2017.0976

  both positive, as expected from what is known about shoot tropisms, namely upward gravitropism and positive phototropism.

). In gure 7a, the t is fair with parametric values of α and Comment citer ce document : Duchemin, L., Eloy, C., Badel, E., Moulia, B. (2018). Tree crowns grow into self-similar shapes controlled by gravity and light sensing. Journal of The Royal Society Interface, 15 (142), 10 p. , DOI : 10.1098/rsif.2017.0976 α p

ACKNOWLEDGMENTS

We warmly thank C. Coutand, A. Lacointe, X. Leoncini and B. Roman for stimulating discussions.

AUTHOR CONTRIBUTIONS

All authors conceived the study. L.D. and C.E. developed the numerical and analytical model. E.B. contributed to tree pictures. All authors wrote the paper.

DATA ACCESSIBILITY

The datasets analysed during the study are available from the corresponding author on reasonable request.

FUNDING STATEMENT

The project leading to this publication has received funding from Excellence Initiative of Aix-Marseille University -A*MIDEX, a French "Investissements d'Avenir" programme. It has been carried out in the framework of the Labex MEC. We also acknowledge support from the CNRS (Mission pour l'interdisciplinarité, project ARBRE).