Keywords: breeding system evolution, dioecy, evolutionary dead end, molecular dating, ancestral state reconstruction, fossil calibration, Ficus, Moraceae

Angiosperms are the most diversified clade of extant plants and are exceptionally species-rich in tropical regions. Flowers are the breeding organs of angiosperms. Flowers exhibit remarkable levels of both structural and functional diversity and for this reason have long thought to have a direct influence on the diversity and evolution of angiosperms.

In this thesis, I investigated breeding system evolution and biogeographic history in the family Moraceae (40 genera, ca. 1100 species, including its largest genus, Ficus, ca. 750 species), which I used as a model clade to understand the origin and evolution of diversity of angiosperms.

In Chapter I, I reconstructed and calibrated a new dated phylogenetic tree for Moraceae as a whole, based on a dataset of 320 species (including 272 species sampled from 36 genera in Moraceae, and 48 species representing 8 outgroup families in Rosales), eight molecular markers (three chloroplastic, five nuclear), and twelve fossil age constraints. I then used this phylogenetic tree to reconstruct ancestral states of breeding systems in Moraceae and Ficus using parsimony and model-based approaches, using and comparing six configurations of breeding system data in Moraceae, differing in the number of states recognized. The crown group ages of Moraceae and Ficus were estimated in the Cretaceous (73.2-84.7 Ma) and in the Eocene (40.6-55.9 Ma), respectively. Dioecy was inferred as the ancestral breeding systems of Moraceae with high support by all the approaches, models and configurations, followed by several subsequent transitions to monoecy, including in Ficus. This result suggests that dioecy is not necessarily an evolutionary dead end. While monoecy emerged as the most likely ancestral state of Ficus, this result remained uncertain and sensitive to model selection.

IV

In Chapter II, I reconstructed a dated phylogenetic tree for tribe Dorstenieae (15 genera, ca. 156 species, including Dorstenia, ca. 113 species), mainly distributed in tropical regions.

To do so, I produced a new dataset of nuclear genomic data (102 genes) generated with a target enrichment approach (Hyb-Seq), sampling all genera (15) and 83 (53%) species in the tribe and five outgroups. This phylogenomic tree allowed me to reconstruct the biogeographic history of this group using dispersal-extinction-cladogenesis models, including models with a founder-event speciation process. The crown group ages of Dorstenieae and Dorstenia were estimated in the Cretaceous (65.8-79.8 Ma) and in the Cretaceous/Paleocene to Eocene period (50.8-67.3 Ma), respectively. Two long-distance dispersal events from continental Africa to South America occurred in the Cenozoic, one in Dorstenia and one in Brosimum s.l. Species of Dorstenia further colonized Central America during the Oligocene to Miocene (12.0-34.7 Ma).

In Chapter III, I tested the climatic niche difference (temperature and precipitation) between the two breeding systems (monoecy and gynodioecy) in Ficus using a new dataset of cleaned spatial occurrence records and breeding systems for 183 species. To do so, I used two comparative approaches, differing in whether phylogeny is taken into account (generalized estimating equations, GEE) or not (generalized linear models, GLM). A positive relationship between precipitation and gynodioecy was supported by GLM, but not by GEE analyses, and no relationship between temperature and breeding systems was supported by either method. Higher dispersal ability and the potential for self-fertilization may explain why monoecious species of Ficus have been able to colonize and survive in drier environments. This thesis highlights the potential of phylogenetic comparative methods and V phylogenomic data to address questions of breeding system evolution and biogeography in Moraceae. More densely sampled phylogenies of Ficus and Moraceae as a whole will be required to confirm some of the results emerging from this thesis, such as ancestral monoecy in Ficus. Nevertheless, this thesis opens up several important new perspectives worth investigating in other plant clades, such as a relationship between breeding system and climatic niche.

on both research as well as my career have been invaluable. I still remember when I introduced my idea, which finally became chapter III of my thesis, to Hervé at the very beginning. This idea has been elaborated and finally become a PhD thesis with help from many people I will mention below.

I would also like to extend my deepest appreciation to my thesis committee, without whose help I will not improve so much in my studies. They are Dr. Thomas Couvreur, Dr. Elliot Gardner, Dr. Damien Hinsinger, Prof. Sophie Nadot, Prof. Jacqui Shykoff and Prof. Nina Rønsted.

I am extremely grateful to my jury of defense: Prof. Mathilde Dufaÿ, Dr. Jean-Yves Rasplus, Prof. Tatiana Giraud, Dr. Florian Jabbour, Dr. Jos Käfer and Dr. Hervé Sauquet. Thank them for their participation in my defense, helpful comments to my thesis and future studies.

I deeply indebted to my coauthors as well. I am grateful to Dr. Stefan Little for his profound knowledge in fossils. I would also like to extend my sincere thanks to Dr. Elliot Gardner and Prof. Nyree Zerega, for their insightful guides on experiments and analyses in Chapter II. I must thank them again for sharing their collections of Neotropical species of Dorstenieae with me.

I also wish to thank Véronique Normand and Alodie Snirc for their great help in my experiments.

II

Thanks should also go to all the members in team EVA (Evolution of Angiosperms), especially Prof. Sophie Nadot, who took care of me when Hervé was not in Orsay. I wish to thank Xavier Aubriot for sharing much knowledge in taxonomy and biogeography with me. Thanks also to Xiaolei Pang, who introduced the knowledge of ITCZ to me. I gratefully acknowledge Martine Fournier and Prof. Jacqui Shykoff's great help in the daily paper works. I would like to acknowledge Emmanuelle Jestin's excellent work in arranging the travel and accomodation for my Jury of thesis defense as well.

I had great pleasure of working with my three kind office mates Fantin Carpentier, Fanny Hartmann, Marine Duhamel and other members in team GEE (Evolutionary genetics and ecology), who brought much social life to me.

Finally, I cannot begin to express my thanks to my parents, without whose support I will never leave home so far away and have a chance to chase my dream. Thanks to Prof.

Xiaoyong Chen, who helped me to apply the CSC fellowship. Thanks to my three best friends: Xiaokun Zhang, Sen Lin, Dazheng Dong, for their company, no matter how physically far away we are. Thanks also to all those who helped and accompanied with me in these years. It is a pity that the space is too limmited to list all of you, but I remember you and all your help.

GENERAL INTRODUCTION

Angiosperms are the dominant plant lineage in most terrestrial ecosystems, especially in tropical regions. Breeding system evolution and biogeography are two of the keys to understand the extreme diversity of angiosperms. My PhD thesis uses the mulberry family, Moraceae, as a model clade to address these topics by answering three key questions: how did breeding systems evolve in Moraceae (Chapter I)? What is the biogeographic history of Dorstenieae, the second most diversified tribe in Moraceae (Chapter II)? Do the two breeding systems in Ficus occupy different climatic niches (Chapter III)? This introduction will present the background of three important aspects of my thesis: 1) broader biological background, introducing current knowledge of relationship of breeding system evolution and diversity in angiosperms, and the biogeographic history of the angiosperms in three main tropical regions; 2) the plant family Moraceae and the two large subclades studied in more detail in this thesis, genus Ficus and tribe Dorstenieae; 3) specific methods used in my thesis, including the target enrichment approach for genomic sequencing (Chapter II); fossil-calibrated molecular dating (Chapters I and II); parametric models for reconstruction of ancestral states (Chapter I) and biogeographic history (Chapter II).

Breeding system evolution and diversity in angiosperms

The term "breeding system" has not been used consistently throughout the literature and is sometimes treated as a synonym of "reproductive system" or "mating system" [START_REF] Neal | Are 'mating systems''breeding systems' of inconsistent and confusing terminology in plant reproductive biology? or is it the other way around?[END_REF][START_REF] Cardoso | Towards a unified terminology for angiosperm reproductive systems[END_REF]. In this thesis, to remain consistent with previous work on Ficus and Moraceae, I chose the term "breeding system" to describe the arrangement of flowers of different sex at different levels, including both individual and population levels.

Angiosperms are the most diversified lineage of extant plants, accounting for approximately 96% of vascular plants [START_REF] Christenhusz | The number of known plants species in the world and its annual increase[END_REF]. The diversity of their species and morphological attributes, including their flowers, have fascinated biologists for hundreds of years.

Because breeding systems directly impact genetic variation and reproductive success at both intraand interspecific levels, they are critical to understanding floral evolution and diversification.

Bisexual flowers are likely ancestral in angiosperms and subsequently evolved many times independently to unisexual flowers associated with a diversity of breeding systems (Sauquet et al. 2017). Dioecy is rare (ca. 6% species of angiosperms) but evolved independently thousands of times during the history of angiosperms (Renner 2014). This observation has prompted multiple lines of research on dioecy evolution, including resource allocation strategies, sexual selection, genetic determination of separate sexes in angiosperms [START_REF] Case | Environmental stress and the evolution of dioecy: Wurmbea dioica (Colchicaceae) in Western Australia[END_REF]Dufay et al. 2014;[START_REF] Charlesworth | Plant contributions to our understanding of sex chromosome evolution[END_REF]Käfer et al. 2017;[START_REF] Zemp | Has adaptation occurred in males and females since separate sexes evolved in the plant Silene latifolia?[END_REF]. Dioecy was once suggested to be an evolutionary dead end, following the assumption that loss of bisexuality appears to be easier than its gain and the observation that lower species richness was found in dioecious clades (Bull and Charnov 1985;Heilbuth 2000). However, this hypothesis may have been an artificial result from the sister-group comparison approach and has been challenged by two recent studies with new comparative approaches (Käfer et al. 2014;Sabath et al. 2016). Although bisexuality is predominant in angiosperms, no significant direction of breeding system transition was found in an angiosperm-wide metaanalysis, suggesting that different selective pressures and constraints have been at play in different clades for breeding system evolution in angiosperms (Goldberg et al. 2017).

Biogeographic history of angiosperms in the tropics

Tropical regions hold the highest terrestrial diversity on earth [START_REF] Mutke | Patterns of vascular plant diversity at continental to global scales[END_REF].

Species richness typically decreases from tropics to poles, a pattern commonly referred to as the latitudinal gradient and observed in most clades [START_REF] Wiens | Historical biogeography, ecology and species richness[END_REF], with a few exceptions (e.g., ray-finned fishes or Actinopterygii; [START_REF] Rabosky | An inverse latitudinal gradient in speciation rate for marine fishes[END_REF]. Before the Middle Miocene Climatic Transition (MMCT, ca. 14 Ma), tropical climates covered a much larger area than they do now [START_REF] Morley | Interplate dispersal paths for megathermal angiosperms[END_REF]. Thus, investigations on the origin, evolution and maintenance of diversity in tropical areas is critical not only to help us understand present patterns of diversity distribution, but also to predict their evolution in a changing environment. The origin of diversity of angiosperms in tropical regions, especially in the Neotropics and Africa, is one of the main questions in this thesis. Hence, I will briefly summarize below current knowledge of diversity of angiosperms in the Neotropics and Africa.

Biodiversity is not evenly distributed across tropical regions. Here, we consider the Neotropics as the area extending from central Mexico to Southern Brazil (including the Caribbean islands) (Antonelli and Sanmartín 2011). The Neotropics have been estimated to harbor the most number of vascular plant species, followed by the Indo-Pacific region [START_REF] Kreft | Global patterns and determinants of vascular plant diversity[END_REF][START_REF] Ulloa | An integrated assessment of the vascular plant species of the Americas[END_REF]. Both dispersal followed by in situ diversification and Andean uplift have been proposed as the main causes for extreme diversification in the Neotropical region [START_REF] Hoorn | Amazonia through time: Andean uplift, climate change, landscape evolution, and biodiversity[END_REF]). In addition, approximately ten biomes have been identified in the Neotropics (Antonelli and Sanmartín 2011;[START_REF] Hughes | Mating advantage for rare males in wild guppy populations[END_REF][START_REF] Antonelli | Amazonia is the primary source of Neotropical biodiversity[END_REF], suggesting high habitat diversity. On the one hand, biota exchange occurred between North and South America through the Central American land bridge during the Late Cretaceous to Eocene, and through the Panama Isthmus during the Oligocene-Miocene transition [START_REF] Morley | Interplate dispersal paths for megathermal angiosperms[END_REF]Jaramillo et al. 2006;Bacon et al. 2015).

Long-distance dispersal from continental Africa (e.g., Ficus section Americanae [Moraceae],

Begonia [Begoniaceae], Lobeliaceae), Asia-Oceania (e.g. Bocconia and Macleaya [Papaveraceae],

Osmorhiza [Apiaceae]) and Australia (e.g., Liliales) also contributed to the diversity of the Neotropics [START_REF] Moonlight | Continental-scale diversification patterns in a megadiverse genus: the biogeography of Neotropical Begonia[END_REF][START_REF] Yi | Chloroplast capture and intra-and inter-continental biogeographic diversification in the Asian -New World disjunct plant genus Osmorhiza (Apiaceae)[END_REF][START_REF] Givnish | Phylogenomics and historical biogeography of the monocot order Liliales: out of Australia and through Antarctica[END_REF][START_REF] Knox | The East Asian origin of the giant lobelias[END_REF]Pederneiras et al. 2018;[START_REF] Li | Differential diversifications of South American and Eastern Asian disjunct genera Bocconia and Macleaya (Papaveraceae)[END_REF]. On the other hand, analyses of biota interchange (including invertebrates, vertebrates, ferns, and angiosperms) among different Neotropical regions suggested Amazonia as the source of Neotropical diversity [START_REF] Antonelli | Amazonia is the primary source of Neotropical biodiversity[END_REF]. However, the drivers of Neotropical diversification still remain unclear at present. Asia-Oceania consists of a heterogenous assemblage of part of Laurasia and several fragments of Gondwana, and these fragments have reached their current positions in different times from the Cretaceous to the Miocene [START_REF] Hall | Southeast Asia's changing palaeogeography[END_REF]. This region has been proposed as the place of origin of angiosperms [START_REF] Buerki | Proto-South-East Asia as a trigger of early angiosperm diversification[END_REF][START_REF] Coiro | How deep is the conflict between molecular and fossil evidence on the age of angiosperms?[END_REF]. However, more evidence are needed to be found for this hypothesis. Borneo and continental South-East Asia have been suggested to be the major evolutionary hotspots for Southeast Asian Biodiversity during the latest 65.5 Ma in a study examining the assembly of local biota (including invertebrates, vertebrates, and plants) [START_REF] De Bruyn | Borneo and Indochina are major evolutionary hotspots for Southeast Asian biodiversity[END_REF]. Floristic exchange between Sunda and Sahul has been estimated to start approximately 33 Ma ago and proceed mainly eastwards (from Sunda including Malay Peninsula, Sumatra, Borneo to Sahul including New Guinea and other Pacific islands) [START_REF] Crayn | The Sahul-Sunda floristic exchange: dated molecular phylogenies document Cenozoic intercontinental dispersal dynamics[END_REF], suggesting the relationship of floras of South-East Asia and Pacific islands.

It remains obscure why continental Africa holds comparatively lower diversity than in the Neotropics and in South-East Asia. [START_REF] Couvreur | Odd man out: why are there fewer plant species in African rain forests?[END_REF] reviewed several hypotheses for the low diversity of tropical Africa (e.g., time-to-speciation effect, area-for-speciation effect), and suggested that lower speciation rates might be one potential explanation. Although continental Africa has been geographically isolated from the Late Cretaceous, when land bridges between continental Africa and Europe, South America, and the Indian Plate occurred, to the Middle Miocene [START_REF] Morley | Interplate dispersal paths for megathermal angiosperms[END_REF], evidence suggests that biota exchange has been continuing through longdistance dispersal. Six floristic regions have been defined in sub-Saharan Africa by geographical analyses, and these floristic regions show relationships with counterparts of other continents (Linder 2014), suggesting flora exchanges in the assembly of the African flora. Although the diversity in Africa as a whole is comparatively lower, clades such as Detarioideae (Fabaceae) harbor higher diversity in continental Africa than other tropical areas, suggesting clades of angiosperms have responded differently to historical climatic change in Africa (de la Estrella et al.

2017). Investigation of the biogeographic history of additional clades of angiosperms in Africa,

based on solid dated phylogenies, is critical to further understand the causes of low diversity in continental Africa [START_REF] Couvreur | Odd man out: why are there fewer plant species in African rain forests?[END_REF].

Moraceae: an ideal model clade

General presentation of family Moraceae

The angiosperm family Moraceae is an ideal model for investigating in breeding system evolution and the origin of diversity in tropical regions. Moraceae are now classified in order Rosales, together with eight other families (APG IV 2016). Urticaceae are well supported as the sister group of Moraceae (Soltis et al. 2011;Zhang et al. 2011). Moraceae consist of seven tribes, 40 genera and around 1100 species (Clement and Weiblen 2009;Zerega and Gardner 2019).

Moraceae have a pantropical geographic distribution with some species extending into temperate areas (Fig. 1). The family is characterized by milky latex, small, unisexual flowers (pistillode present or not), a uniseriate perianth, and compact but variable inflorescences [START_REF] Berg | Moraceae-genera other than Ficus[END_REF]Simpson 2010). The family exhibits high diversity in several key traits, including habit (tree, shrub, herb, woody liana), inflorescence architecture (from the open inflorescences of Morus to the closed ones of Ficus), pollination (by wind, e.g. Morus; or by insects e.g. Ficus) (Rohwer and Berg 1993) and breeding systems (monoecy, androdioecy, gynodioecy, dioecy) (Fig. 2) (Clement and Weiblen 2009). The seeds of Moraceae are mostly dispersed by vertebrates (Rohwer and Berg 1993;Shanahan et al. 2001). Dioecy has been proposed to be the ancestral breeding system of Moraceae (Datwyler and Weiblen 2004). The crown-group age of Moraceae has been estimated in the Cretaceous (72.6-110.0 Ma) (Zerega et al. 2005), while biogeographic analyses for Moraceae as a whole using model-based approaches have not been conducted so far. The ancestral area of Moraceae remains an enigma and may have been either in Gondwana or in Laurasia (Zerega et al. 2005). Except for Ficus (see below), biogeography analyses have been conducted so far in three genera (Artocarpus, Maclura and Dorstenia) in Moraceae to date. Artocarpus (ca. 70 spp.) has been estimated to originate in Borneo in the Eocene to and subsequently disperse to other Asia-Oceania areas ( Williams et al. 2017). Maclura (ca. 12 spp.) has been estimated to originate in South America in the Cretaceous to , followed by dispersal to Africa and the northern hemisphere (Gardner et al. 2017). The particular case of Dorstenia is developed further below.

Moraceae harbors several species of economic value and some of them are widely cultivated.

Leaves of mulberry (Morus) provide food for silk worms (Fig. 2 A, B, [START_REF] He | Draft genome sequence of the mulberry tree Morus notabilis[END_REF]. Artocarpus heterophyllus (Fig. 2 C) and Ficus carica (Fig. 2 I) are cultivated for their edible infructescences known as breadfruit and fig, respectively (Ghada et al. 2010;[START_REF] Williams | Out of Borneo: biogeography, phylogeny and divergence date estimates of Artocarpus (Moraceae)[END_REF]. Some species of Dorstenia have been used in folk medicine for the treatment of infection, snake bites, and as anti-inflammatory in Africa and South America [START_REF] Zapata-Sudo | Sedative and anticonvulsant activities of methanol extract of Dorstenia arifolia in mice[END_REF]. Lastly, the fibrous bark of Broussonetia papyrifera, known as paper mulberry, was used to make paper in ancient China [START_REF] Chung | Molecular recircumscription of Broussonetia (Moraceae) and the identity and taxonomic status of B. kaempferi var. australis[END_REF].

Genus Ficus

The most diversified genus in Moraceae is Ficus (ca. 750 spp.), which is known for its coevolution with its symbiotic pollinating wasps (Hymenoptera, Chalcidoidea, Agaonidae) (Cruaud et al. 2012). Ficus is characterized by very small unisexual flowers enclosed in urn-shape inflorescences called syconia (Fig. 2 I). There are three kinds of flowers in Ficus inflorescences: male flowers (producing pollen), female flowers (producing seeds) and gall flowers (develop wasps). Two breeding systems occur in Ficus, monoecy and gynodioecy, differing in the distribution of these three kinds of flowers among syconia (Fig. 3), each representing about half of the species of Ficus (Cook and Rasplus 2003). In monoecious Ficus, all three kinds of flowers co-occur in every syconium (Fig. 3). Pollinating wasps lay eggs and pollinate the host plants after entering the syconium, depending on the lengths of styles of female flowers in monoecious Ficus. Conversely, gynodioecious species of Ficus are characterized by two kinds of syconia: male and female (Fig. 3). On the one hand, functional female syconia contain only female flowers. On the other hand, functional male syconia contain both male and gall flowers, and as a result will not deliver seeds.

The male pollinating wasps are wingless, and usually never leave the syconium in which they developed. When the mature winged female wasps leave the male syconium in which they developed, they have to pass the cluster of male flowers near the exit of the syconium. These pollinating wasps then look for a new syconium and struggle into it (carrying pollen from another syconium). They may lose their wings and some legs during the process. Following entry into a male syconium, they lay eggs; if they enter into a female syconium, they pollinate the female flowers (Cook and Rasplus 2003). Thus, structural gynodioecy in Ficus (plants are structurally either bisexual, with monoecious inflorescences, or female only) is functionally equivalent to dioecy (plants are functionally either male or female only). Ficus has a broad, mostly pantropical geographic distribution across all the three major tropical regions, but gynodioecious Ficus are remarkably absent in the Neotropics (Cruaud et al. 2012). Monoecy has been proposed to be the ancestral state of breeding system of Ficus (Datwyler and Weiblen 2004).

Six subgenera were proposed for Ficus [START_REF] Berg | Flora Malesiana precursor for the treatment of Moraceae 1: the main subdivision of Ficus: the subgenera[END_REF]. Three of them have so far been confirmed to be monophyletic, whereas subgenera Ficus, Pharmacosycea, and Urostigma appear to be polyphyletic in the phylogeny of Cruaud et al. (2012) based on five nuclear markers and 200 species. Pederneiras et al. reconstructed the most densely sampled phylogenetic analysis for Ficus to date (249 species) with nine nuclear molecular markers (Pederneiras et al. 2018). The topology of Pederneiras et al. (2018) was similar to the one of Cruaud et al. (2012), but with higher support for the deeper nodes. With the development of high-throughput sequencing, initial phylogenomic studies of Ficus have also recently been conducted based on much larger gene datasets but more limited taxon sampling (Bruun-Lund et al. 2017;Rasplus et al. 2018). Phylogenetic relationships based on complete chloroplast genomes (Bruun-Lund et al. 2017) differ markedly from phylogenies based on a few nuclear markers (Cruaud et al. 2012;Pederneiras et al. 2018), suggesting hybridization in the early history of Ficus. However, phylogenetic relationships based on nuclear RAD-seq genomic data (Rasplus et al. 2018) also differ from previous nuclear phylogenies (Cruaud et al. 2012;Pederneiras et al. 2018), specifically in relationships at the base of Ficus.

The crown-group age of Ficus has been estimated in the Cretaceous to the Paleocene (60.0-101.9 Ma) or in the Eocene (34.9-50.6 Ma), depending on different calibration strategies (Cruaud et al. 2012;Pederneiras et al. 2018). The ancestral area of Ficus has been inferred to be Eurasia (Cruaud et al. 2012;Pederneiras et al. 2018), from which Ficus then dispersed from the Old World to the New World twice, leading to the two Neotropical lineages (sections Pharmacosycea and Americanae) (Pederneiras et al. 2018). Ficus subsect. Urostigma may have originated in the Paleocene to in Madagascar [START_REF] Chantarasuwan | Palaeotropical intercontinental disjunctions revisited using a dated phylogenetic hypothesis with nearly complete species level sampling of Ficus subsect. Urostigma (Moraceae)[END_REF].

Tribe Dorstenieae

Tribe Dorstenieae includes fifteen genera and 156 species, distributed on both sides of the Atlantic (Figure 1, Chapter II). Dorstenia, with ca. 113 species, is the most diversified genus in Dorstenieae and the second in Moraceae (Berg and Hijman 1999;Berg 2001).

New species have been discovered continuously in Africa and South America during the last decade (Mccoy and Massara 2008;dos Santos and Neto 2012;Machado and Marcelo Filho 2012;[START_REF] Chase | Dorstenia christenhuszii (Moraceae), a new species from the Taita Hills, Kenya[END_REF]dos Santos et al. 2013;Leal 2014;Machado et al. 2014;Rzepecky 2016). Dorstenia has very small unisexual flowers organized and partly embedded in disk-like monoecious inflorescences, except D. lavrani which is dioecious (Mccoy and Massara 2008). Dorstenia is further characterized by the marginal appendages around the inflorescence (i.e., extensions of the flattened inflorescence receptacle), and the length of these appendages varies among species (Fig. 2 H). Inflorescence shape in Dorstenia ranges from square to round. Species of Dorstenia are diverse in habit, ranging from treelets to shrubs to herbs and the genus contains almost exclusively all of the herbaceous species known in Moraceae (Berg and Hijman 1999;Berg 2001).

Most species of Dorstenia are distributed in either South America or Africa (incl.

Madagascar and Arabian Peninsula). Dorstenia has been estimated to originate in Africa in the Cretaceous (84.8-132.0 Ma) (Misiewicz and Zerega 2012). However, this estimate is older than the estimate of Moraceae as a whole (Zerega et al. 2005) and likely needs to be revised. Misiewicz and Zerega (2012) reconstructed the first phylogeny of Dorstenia using ITS sequences, but some of the deep nodes remained weakly supported. Interestingly, their results suggested a single dispersal into the Neotropics with a subsequent re-colonization of Africa. However, sampling of both ingroup (32 out of 113 species) and outgroups (seven species represented seven genera, including two in the same tribe) was limited and the phylogeny of tribe Dorstenieae as a whole remains incompletely understood.

Methodological approaches used in this thesis

During my PhD program, I took advantage of recent methodological improvements in phylogenetics and macroevolution, including both technical advances (high-throughput sequencing) and conceptual developments in phylogenetic tree reconstruction and dating, model-based biogeographic history reconstruction, and phylogenetic comparative methods.

Hyb-seq approaches

The rapid development of high-throughput sequencing (HTS) has revolutionized phylogenetic analysis in the last decade [START_REF] Drew | Another look at the root of the angiosperms reveals a familiar tale[END_REF][START_REF] Matasci | Data access for the 1,000 Plants (1KP) project[END_REF]Couvreur et al. 2019;[START_REF] Villaverde | Bridging the micro-and macroevolutionary levels in phylogenomics: Hyb-Seq solves relationships from populations to species and above[END_REF]. The five widely used HTS approaches at present are: 1) microfluidic PCR, based on PCR amplification of targeted regions; 2) restriction enzyme-based methods, using restriction enzymes to separate genomic DNA and sequence, for instance, RAD-seq; 3) genome skimming, consisting of total genomic DNA sequencing at low depth without enrichment (typically aimed at recovering organellar genomes); 4) target enrichment approach (incl. the strategy commonly referred to as "exon capture"), enriching shotgun sequencing libraries for target genes; 5) transcriptome sequencing (McKain et al. 2018). The Hyb-seq approach differs from other target enrichment strategies by not only targeting thousands of low-copy nuclear exons, but also their flanking regions, high-copy repeats, and organellar genomes simultaneously (Weitemier et al. 2014). One of the advantages of the target enrichment and Hyb-seq approaches is their suitability to work with degraded herbarium specimens (McKain et al. 2018). Herbarium collections are a treasure for biodiversity studies, when we consider the information about diversity and distribution from them, and the fact that some species are very difficult to be collected in the field (Hart et al. 2016).

Molecular dating

Molecular dating, which transforms the relative time from molecular branch lengths of phylogenetic trees to absolute time, using calibrations typically derived from the fossil record. Molecular dating is essential to connect the evolutionary history of a studied group of organisms with extrinsic geographic and climatic events (Sauquet 2013). Originally developed under the assumption of a strict clock (constant molecular rate) model, current molecular dating methods now allow heterogenous rates (relaxed clock models). Fossil calibrations have traditionally been applied as minimum age constraints. More recently, following the development of Bayesian approaches to molecular dating, various prior distributions (e.g., lognormal, exponential) have been implemented to describe the expected time between the fossil and calibrated node age (Ho and Phillips 2009). Further, to mitigate the impact of potential errors in calibration and the conflicts among calibration nodes, soft boundaries (allowing non-zero probability beyond both minimum and maximum age constraints) were introduced in prior distributions for calibrations (Yang and Rannala 2005). However, even more recent molecular dating approaches (not used in this thesis), including the fossilized birth-death process and the so-called total-evidence (or tipdating) methods treat fossils more explicitly as tips (and/or ancestors) of the phylogenetic tree, thereby avoiding the arbitrariness of prior distribution parameterization [START_REF] Pyron | Divergence time estimation using fossils as terminal taxa and the origins of Lissamphibia[END_REF]Ronquist et al. 2012;[START_REF] Heath | The fossilized birth-death process for coherent calibration of divergence-time estimates[END_REF].

The development of molecular dating methods in the genomic era has mainly focused on handling the much larger datasets than in the Sanger sequencing era, including the challenge of adequately partitioning such datasets [START_REF] Foster | Estimating the number and assignment of clock models in analyses of multigene datasets[END_REF]Foster and Ho 2017).

Although widely used software such as BEAST 2 [START_REF] Bouckaert | BEAST 2: A Software Platform for Bayesian Evolutionary Analysis[END_REF]) may still be used in theory to analyze genomic data, computational times are in practice a limiting factor.

New software such as MCMCTree in package PAML (Yang 2007) use approximate likelihood to solve this problem. Contrary to some early expectations, phylogenomic data so far appear to have a limited impact on divergence times, whereas fossil calibrations remain the most important factor [START_REF] Dos Reis | Bayesian molecular clock dating of species divergences in the genomics era[END_REF]Foster and Ho 2017). Thus, robust and well justified fossil calibration remains as critical with phylogenomic datasets as in the Sanger sequencing era (Sauquet et al. 2012;[START_REF] Ksepka | The Fossil Calibration Database-A new resource for divergence dating[END_REF].

Parametric models for ancestral state reconstruction

Ancestral state reconstruction (ASR) helps evolutionary biologists to understand the transition of traits in the evolutionary history of a clade. Parsimony approaches, which find the states that minimize the number of transition events given the states on tips of the phylogeny, were first used in ASR [START_REF] Maddison | Outgroup analysis and parsimony[END_REF]). However, model-based approaches are now preferred by many, as they consider the probability of transition events according to the branch lengths and allow to measure the relative confidence of each state at internal nodes of the phylogeny, expressed as proportional marginal likelihoods [START_REF] Pagel | The maximum likelihood approach to reconstructing ancestral character states of discrete characters on phylogenies[END_REF]. Asymmetrical (directional) transition can also be estimated with ML approaches, allowing biologists to test more hypotheses in studies of trait evolution. Originally developed under a maximum likelihood (ML) framework, model-based approaches to ASR are also now available in a Bayesian framework. Bayesian approaches to ASR take into account uncertainty of both phylogenies (topologies and branch lengths) and model parameters (Pagel et al. 2004;[START_REF] Ronquist | Bayesian inference of character evolution[END_REF]). The reversible-jump Markov Chain Monte Carlo approach takes further advantage of the Bayesian framework to explore and visit multiple models of morphological evolution in proportion to their posterior probabilities (Pagel and Meade 2006).

Parametric models for historical biogeography

The reconstruction of ancestral distribution areas is similar to ancestral state reconstruction, and both fields have similar histories of methodological improvement.

Historical biogeographic reconstruction started with parsimony-based method, for instance event-based approaches [START_REF] Ronquist | Parsimony analysis of coevolving species associations In: Tangled trees: phylogeny, cospeciation and coevolution[END_REF]. Inspired by the methods of ancestral state reconstruction and phylogenetic model selection, parametric models were introduced to historical biogeographic reconstruction about ten years ago, for instance the now widely used dispersal-extinction-cladogenesis (DEC) model [START_REF] Ree | Maximum likelihood inference of geographic range evolution by dispersal, local extinction, and cladogenesis[END_REF]Ree and Sanmartín 2009). Parametric biogeographic methods not only consider the topology but also the relative or absolute time, which parsimony methods do not. As in ASR models, the amount of biogeographic change (range expansion or local expansion) is proportional to branch lengths (relative time). Absolute time links the historical biogeographic reconstruction with extrinsic geographic events (e.g., land bridges, fluctuation of sea level, climate in certain area). When multiple parametric models are considered, model selection metrics such as Akaike Information Criterion (AIC) and Bayes Factor (BF) may be used to determine the most appropriate model to explain the data.

Objectives of this thesis

This thesis focuses on breeding system evolution and the origin of diversity of angiosperm in tropical regions by using Moraceae as a model clade. My thesis is composed of three chapters, each written as a standalone manuscript for publication. Chapter I mainly focuses on breeding system evolution in Moraceae as a whole. To tackle this question, a new dated phylogenetic tree and ancestral states of breeding systems were reconstructed in Moraceae and Ficus. This chapter was recently published in Annals of Botany (Zhang et al. 2019a). In Chapter II, I investigated the biogeographic history of Dorstenieae, the second most species-rich tribe of Moraceae. To do so, a new dated phylogenomic tree of tribe Dorstenieae, sampled from herbarium specimens, was reconstructed and used to address the biogeographic history of the group. This chapter was recently made available as a preprint in bioRxiv and is currently being considered for recommendation by Peer Community in Evolutionary Biology [START_REF] Zhang | Long-distance dispersal shaped the diversity of tribe Dorstenieae (Moraceae)[END_REF]

ABSTRACT

 Background and Aims Although dioecy, which characterizes only 6% of angiosperm species, has been considered an evolutionary dead end, recent studies have demonstrated that this is not necessarily the case. Moraceae (40 genera, 1100 spp., incl.

Ficus, 750 spp.) are particularly diverse in breeding systems (incl. monoecy, gynodioecy, androdioecy, and dioecy) and thus represent a model clade to study macroevolution of dioecy.

 Methods Ancestral breeding systems of Ficus and Moraceae were inferred. To do so, a new dated phylogenetic trees of Ficus and Moraceae was first reconstructed by combining a revised twelve-fossil calibration set and a densely sampled molecular dataset of eight markers and 320 species. Breeding system evolution was then reconstructed using both parsimony and model-based (maximum likelihood and Bayesian) approaches with this new timescale.

INTRODUCTION

Flowers are the reproductive structures of angiosperms. The ca. 350,000 species of angiosperms are highly diverse in floral morphological traits, including breeding systems, ranging from bisexuality (hermaphroditism) to separate sex on distinct plants (dioecy) through several presumed intermediate states between these two ends (monoecy, andromonoecy, gynomonoecy, androdioecy, gynodioecy) (Renner 2014). Although a recent study demonstrated that bisexual flowers are ancestral in angiosperms and evolved many times independently to unisexual flowers (Sauquet et al. 2017), the exact number and the context of these transitions remains to be characterized. Dioecy is rare in angiosperms (only ca. 6%) (Renner 2014) and has been suggested to be an evolutionary dead end (Bull and Charnov 1985;Heilbuth 2000). Transitions from hermaphroditism to dioecy have been suggested to occur through three alternative pathways: dimorphic pathway (through gynodioecy, androdioecy, or polygamodioecy) (Dufay et al. 2014) ; monomorphic pathway (through monoecy, andromonoecy, gynomonoecy, or polygamomonoecy) [START_REF] Renner | Dioecy and its correlates in the flowering plants[END_REF]; and direct pathway [START_REF] Barrett | The evolution of plant sexual diversity[END_REF]Goldberg et al. 2017). However, the view of dioecy as an evolutionary dead end has gradually changed over the last decade after it was found that transitions from dioecy to other breeding systems are possible and that the flexibility of breeding systems transitions may be caused by different selective pressures and constraints [START_REF] Barrett | The evolution of plant reproductive systems: how often are transitions irreversible?[END_REF]Goldberg et al. 2017). Thus, it remains unclear why dioecy is in fact so rare in angiosperms (Käfer et al. 2017).

The mulberry family (Moraceae), consists of approximately 40 genera and 1100 species, and represents a good model system to study dioecy transitions (Clement and Weiblen 2009). Four breeding systems are observed in the family: monoecy, androdioecy, gynodioecy, and dioecy (Datwyler and Weiblen 2004). Ficus, the largest genus in Moraceae, contains almost 75% of the species in the family (750 spp.). About half of the species of Ficus are monoecious, while the other half are gynodioecious (Cook and Rasplus 2003).

In monoecious species of Ficus, both functionally male and female flowers coexist in the same inflorescence (syconium); female flowers vary in style length. Fig wasp pollinators (Agaonidae) lay eggs in part of these female flowers and pollinate the others (Cook and Rasplus 2003). Gynodioecious species of Ficus are characterized by two kinds of plants and inflorescences: female individuals bear syconia that contain only functionally female flowers (i.e., flowers that can develop seeds), and functionally male individuals bear syconia that contain both functionally male and gall flowers (female flowers that develop only wasps) (Cook and Rasplus 2003). Therefore, structurally gynodioecious species of figs are functionally dioecious.

A previous study, based on parsimony reconstruction sampling 46 Ficus species (but no outgroup), suggested that monoecy is ancestral in Ficus, and that gynodioecy originated at least once or twice within the genus [START_REF] Weiblen | Phylogenetic relationships of functionally dioecious Ficus (Moraceae) based on ribosomal DNA sequences and morphology[END_REF]. A subsequent family-level parsimony study, sampling 83 Moraceae species (incl. 11 species of Ficus) suggested monoecy and dioecy to be ancestral in Ficus and Moraceae, respectively (Datwyler and Weiblen 2004). However, the accuracy and precision of ancestral state reconstruction depends on the reliability of the phylogenetic tree and sampling density, and important progress has been made to reconstruct phylogenetic relationships in Ficus and Moraceae since these two studies were published (Zerega et al. 2005;Clement and Weiblen 2009;Zerega et al. 2010;Cruaud et al. 2012;Williams et al. 2017). In addition, new probabilistic approaches for ancestral state reconstruction, taking into account phylogenetic uncertainty and divergence times, are now available (Pagel et al. 2004). Parsimony, maximum likelihood (ML) and Bayesian approaches today are the most commonly used methods for reconstructing ancestral states. The accuracy of branch length estimation may be important in model-based approaches (ML, Bayesian), which rely on them for computing the likelihood of evolutionary change along the phylogeny (Pagel and Meade 2006). However, to date only few attempts have been made to estimate divergence times in Moraceae (Datwyler and Weiblen 2004;Zerega et al. 2005). Although more studies have addressed the timescale of Ficus evolution [START_REF] Rønsted | 60 million years of co-divergence in the fig-wasp symbiosis[END_REF]Zerega et al. 2005;Cruaud et al. 2012), partly inconsistent results have been obtained, for instance with crown-group age estimates for Ficus ranging from 40.1 to 101.9 Mya. In addition to providing a framework to reconstruct ancestral states, estimating divergence times in Ficus and Moraceae is also important for other evolutionary questions, including the evaluation of the codiversification and biogeographic history of the genus (Cruaud et al. 2012).

Here, we investigate two key questions: 1) what are the crown-group ages of Ficus and Moraceae?; and 2) what are the ancestral breeding systems of Ficus and Moraceae, and how many times did monoecy and dioecy evolve in Ficus and Moraceae? To do so, we reconstruct a new phylogenetic tree for Moraceae using a densely sampled molecular dataset, estimate divergence times using a relaxed molecular clock calibrated with a revised set of 12 fossil age constraints, and reconstruct ancestral breeding systems in Ficus and Moraceae with state-of-the-art model-based approaches. We estimate the age of crown group of Ficus and Moraceae in the Eocene and Late Cretaceous, respectively. Our results suggest that dioecy is not an evolutionary dead end in Moraceae and the transitions from dioecy to the other breeding systems (androdioecy, gynodioecy and monoecy) occurred several times during the evolutionary history of the family. This study sheds light on the evolution of dioecy into other breeding systems and more generally improves our understanding of breeding system evolution in angiosperms.

MATERIALS AND METHODS

Taxonomic sampling and molecular dataset assembly

We selected GenBank sequences from 320 species belonging to 65 genera and eight families of Rosales, including 272 species and 36 genera (3/4 of the circa 40 genera) of Moraceae (Table S1). Our efforts include a comprehensive sample of outgroups in order to estimate accurately divergence times in the family while taking advantage of the rich fossil record of the order (see below). We used eight markers from two genomes: three chloroplast genes (matK, rbcL, ndhF) and five nuclear markers including two noncoding regions (ITS, ETS) and three genes (G3pdh, ncpGS, the waxy region). This combination of coding and non-coding markers was selected to resolve both deep-and shallow-level relationships. MatK and rbcL have been proposed as standard barcoding regions for land plants [START_REF] Hollingsworth | A DNA barcode for land plants[END_REF] and ndhF has been proved to be very useful in previous phylogenetic studies of Moraceae (Datwyler and Weiblen 2004;Zerega et al. 2005;Clement and Weiblen 2009). ITS, ETS, G3pdh, ncpGS, and the waxy region were used in the latest and most densely sampled phylogenetic study of Ficus (Cruaud et al. 2012).

Considering the comparatively fast molecular rates of ITS and ETS, these markers may be difficult to align and may introduce more noise than signal at the family level. Therefore, ITS and ETS sequences were here only used for species of Ficus. All sequences were aligned using Muscle v3.7 [START_REF] Edgar | MUSCLE: multiple sequence alignment with high accuracy and high throughput[END_REF]) as implemented on CIPRES (Miller et al. 2010); alignments were then checked and adjusted by hand. In the ITS alignment we deleted three short regions (ca. 10 bp) because of the uncertainty of gap length and position of base pairs in these regions. Combined multi-marker alignments were assembled with Mesquite v3.04 [START_REF] Maddison | Mesquite: a modular system for evolutionary analysis[END_REF].

Phylogenetic reconstruction

To identify and exclude problematic GenBank sequences (e.g., incorrectly identified), separate gene analyses were conducted first with RAxML v8.2.9 (Stamatakis 2014) on CIPRES (Miller et al. 2010), using the GTRCAT model and 100 bootstrap replicates. We then compared single-marker phylogenetic trees with each other and with the latest published phylogenies of Ficus (Cruaud et al. 2012), Moraceae (Clement and Weiblen 2009), and other families of Rosales [START_REF] Wiegrefe | The Ulmaceae, one family or two? Evidence from chloroplast DNA restriction site mapping[END_REF][START_REF] Potter | Phylogeny and classification of Rosaceae[END_REF]Zhang et al. 2011;[START_REF] Wu | Molecular phylogeny of the nettle family (Urticaceae) inferred from multiple loci of three genomes and extensive generic sampling[END_REF][START_REF] Yang | Molecular phylogenetics and character evolution of Cannabaceae[END_REF][START_REF] Onstein | Do Mediterraneantype ecosystems have a common history?-Insights from the buckthorn family (Rhamnaceae)[END_REF][START_REF] Hauenschild | Phylogenetic relationships within the cosmopolitan buckthorn family (Rhamnaceae) support the resurrection of Sarcomphalus and the description of Pseudoziziphus gen[END_REF]. A sequence was excluded when its position in one tree conflicted strongly (i.e., with high support) with those in other trees or published phylogenies. After excluding problematic sequences (Celtis: L12638, AY263941, AY263961, AY263925, AY263899; Morus: L01933) and confirming that no supported conflict remained among single-marker trees, combined phylogenetic analyses of chloroplast genes (matK, rbcL, ndhF), nuclear markers (ITS, ETS, G3pdh, ncpGS, the waxy region), and of all eight markers were conducted with maximum likelihood and Bayesian approaches using RAxML and BEAST on CIPRES.

For maximum likelihood analyses, the dataset was divided into eight partitions according to marker, each with the GTRCAT model and 1000 bootstrap replicates. All trees are presented rooted on the most external outgroup, Rosaceae, which have been shown to be the sister group of all remaining Rosales in all recent higher-level phylogenetic analyses (e.g., [START_REF] Wang | Rosid radiation and the rapid rise of angiospermdominated forests[END_REF]Soltis et al. 2011). The phylogenetic reconstruction by Bayesian approach was done simultaneously with divergence time estimation (see below).

Fossil calibration

The fossil record of Moraceae is comparatively poor and in need of critical revision (Collinson 1989), whereas unambiguous fossils exist for other families of Rosales (Burge and Manchester 2008;Benedict et al. 2011;Friis et al. 2011). When reliable fossils are absent or scarce for the ingroup, outgroup calibration may provide more accurate estimates than secondary calibration (Sauquet et al. 2012;[START_REF] Hipsley | Beyond fossil calibrations: realities of molecular clock practices in evolutionary biology[END_REF]. Therefore, we specifically designed the taxon sample of this study to include sufficient outgroup nodes to take advantage of the fossil record of Rosales. Our set of fossil age constraints includes four in Moraceae and eight distributed among the remaining families of Rosales (Table S2).

To revise fossil calibrations in Rosales, we proceeded as follows [START_REF] Parham | Best practices for justifying fossil calibrations[END_REF]Sauquet et al. 2012): 1) we started from lists of calibrations used in previous molecular dating studies (Zerega et al. 2005;Cruaud et al. 2012;Magallón et al. 2015) and completed this list with specific reviews (Collinson 1989;Friis et al. 2011) and recently published fossil taxa [START_REF] Manchester | Biogeographical relationships of North American tertiary floras[END_REF]Calvillo-Canadell and Cevallos-Ferriz 2007;[START_REF] Manos | Phylogeny of extant and fossil Juglandaceae inferred from the integration of molecular and morphological data sets[END_REF]Benedict et al. 2011); 2) for each fossil, we critically assessed the phylogenetic assignment based on original descriptions and subsequent reviews, and using the latest reference phylogeny for each family [START_REF] Manchester | Biogeographical relationships of North American tertiary floras[END_REF]Calvillo-Canadell and Cevallos-Ferriz 2007;[START_REF] Manos | Phylogeny of extant and fossil Juglandaceae inferred from the integration of molecular and morphological data sets[END_REF]Benedict et al. 2011;Friis et al. 2011;[START_REF] Jud | Flowering after disaster: Early Danian buckthorn (Rhamnaceae) flowers and leaves from Patagonia[END_REF]; 3) for each fossil, we also critically revised the absolute age or age range of the fossil using the latest stratigraphy and geological time scale from the International Commission on Stratigraphy [START_REF] Cohen | The ICS International Chronostratigraphic Chart[END_REF]. Because none of these fossil taxa have been included in total evidence phylogenetic analyses, our assignments here are at best 'apomorphy-based' and therefore we have been particularly conservative in both selecting our final calibrations and assigning them to clades. This implies that some taxa previously used as calibrations are typically used here to constrain the age of a more inclusive node than in previous studies (e.g., the fossil achenes we used to calibrate the stem group of Ficus were used to calibrated the crown group of Ficus in former studies; Datwyler and Weiblen 2004;[START_REF] Rønsted | 60 million years of co-divergence in the fig-wasp symbiosis[END_REF]Zerega et al. 2005;Cruaud et al. 2012).

Molecular dating analyses

We used BEAST v1.8.0 [START_REF] Drummond | Bayesian phylogenetics with BEAUti and the BEAST 1.7[END_REF]) implemented on CIPRES (Miller et al. 2010) to estimate divergence times and topology simultaneously. To reduce computational burden, we used empirical base frequencies and divided the dataset into three partitions: chloroplast DNA (matK, rbcL, ndhF), noncoding nuclear markers (ITS, ETS), and coding nuclear genes (G3pdh, ncpGS, the waxy region). The substitution model for each partition was selected using the Akaike Information Criterion (AIC) with jModelTest v2.1.6 [START_REF] Darriba | jModelTest 2: more models, new heuristics and parallel computing[END_REF] as implemented on CIPRES. The best substitution model was TVM+I+G for the chloroplast (matK, rbcL, ndhF) and noncoding nuclear (ITS, ETS) partitions, and TPM3uf+I+G for the coding nuclear (G3pdh, ncpGS, the waxy region) partition.

For each calibrated node, we chose the oldest reliable fossil as a minimum age constraint. Because most fossils typically provide only minimum ages for the clade to which they can be safely attributed, a uniform prior distribution was used for each fossil calibration, using the upper (younger) boundary of the fossil oldest stratigraphic age range as the minimum bound and the maximum constraint for the root (see below) as the maximum bound. Molecular dating analyses usually require at least one maximum age constraint (Ho and Phillips 2009;Sauquet 2013). Here, we chose to set the maximum age bound for the root (i.e., the crown-group node of Rosales) to 107 Ma based on the crowngroup age estimate of in the latest, large-scale molecular dating study of angiosperms (Magallón et al. 2015). Although the phylogenetic relationships among sampled members of Rhamnaceae are inconsistent with previous work [START_REF] Onstein | Do Mediterraneantype ecosystems have a common history?-Insights from the buckthorn family (Rhamnaceae)[END_REF][START_REF] Hauenschild | Phylogenetic relationships within the cosmopolitan buckthorn family (Rhamnaceae) support the resurrection of Sarcomphalus and the description of Pseudoziziphus gen[END_REF]; see Discussion), they are unlikely to have affected divergence times because we chose conservative fossil assignments in Rhamnaceae.

Paliurus clarnensis has been proposed to be more closely related to extant species of Paliurus than any other member of Rhamnaceae (Burge and Manchester 2008), thus providing a minimum age for the stem node of Paliurus. Because of conflict in the position of the genus, we decided to calibrate the crown group node of Rhamnaceae instead.

Coahuilanthus belindae has been proposed to be a member of Rhamnaceae, but its position in the phylogeny remains unclear (Calvillo-Canadell and Cevallos-Ferriz 2007), therefore we used it as a minimum age constraint for the stem node of the family.

The tree prior in BEAST was set as a Birth-Death process. To produce a starting tree that conforms with the hard age constraints enforced in this analysis, we transformed the best-known RAxML tree into an ultrametric tree using the bladj function of phylocom v4.2 [START_REF] Webb | Phylocom: software for the analysis of phylogenetic community structure and trait evolution[END_REF]. Four separate runs were conducted, each with 100 million generations, sampling trees and parameters every ten thousand generations. Chain convergence was checked in Tracer v1.5 (Rambaut and Drummond 2009), with the first 10% of chain length excluded as burn-in. After confirming that chains of these four runs had converged, we combined them using LogCombiner v 1.8.0, resampling every 100 thousand generations and discarding the first 10% (i.e., 10 million) generations as burn-in. Tree statistics were summarized on the Maximum Clade Credibility (MCC) tree by TreeAnnotator v1.8.0, using median ages as node heights.

In addition to our main dating analysis, we also conducted a series of sensitivity experiments. First, we tested the impact of calibrations on the estimated topology by conducting another run keeping the maximum age constraint on the root but without internal age constraints (Sauquet et al. 2012). Second, we produced another run by excluding the age constraint on the stem node of Ficus, but maintaining all other 11 age constraints to test the influence of the former on the estimated crown-group age of Ficus.

Last, we tested the impact of our root maximum age constraint on divergence time estimates for the ingroup. Indeed, there remains considerable uncertainty on the crowngroup age of angiosperms [START_REF] Doyle | Molecular and fossil evidence on the origin of angiosperms[END_REF]Magallón et al. 2015;[START_REF] Herendeen | Palaeobotanical redux: revisiting the age of the angiosperms[END_REF]) and recent analyses that do not constrain this age consistently estimate it to be much older than the first accepted crown-group fossils (Foster et al. 2017). Therefore, we also re-ran our analysis using an older maximum age constraint on the root, set to 112 Ma, based on the crown-group age of estimated in the unconstrained analysis of Magallón et al. (2015), in which the crown-group age of angiosperms was estimated to be 160-256 Ma rather than 136-140 Ma (S. Magallón, pers. comm.).

Ancestral state reconstruction

The breeding system state of each species was determined from the literature, and then recorded in the PROTEUS database [START_REF] Sauquet | PROTEUS: a database for recording morphological data and creating NEXUS matrices[END_REF]. A list of all data records (each linked to an explicit reference) and the matrix are provided as Table S3. For the main analysis, we reconstructed the ancestral state with all 320 species in the dataset (i.e., including all outgroups). To characterize a reliable evolutionary scenario supported by different classification methods of breeding systems, (e.g., [START_REF] Anger | Dioecy in Amborella trichopoda: evidence for genetically based sex determination and its consequences for inferences of the breeding system in early angiosperms[END_REF], we analyzed six configurations of the same character and data, differing in the number of character states considered (Table 1; Fig. 1). Configuration A distinguishes among five breeding system states (as recorded and listed in Table S3). Configuration B (four states) treats androdioecy as missing data. In our dataset, there are only three androdioecious species, accounting for less than 1% of all species sampled. Configuration C (three states) is similar to Configuration B, but with gynodioecy and dioecy pooled as a single state (dioecy). In Ficus, gynodioecious species are functionally dioecious (Cook and Rasplus 2003), and all gynodioecious species but one (Dryas octopetala) in this study belong to Ficus.

Configuration D (three states) is similar to Configuration B, but treats bisexual flowers as missing data. In Rosales, families Cannabaceae, Moraceae and Urticaceae all have unisexual flowers (Simpson 2010), thus maintaining the bisexual state is not essential to this study and could be detrimental by unnecessarily increasing the number of parameters to estimate. Configuration E (two states) combines the two transformations of Configurations C and D: gynodioecy and dioecy pooled, bisexuality treated as missing. In Configuration F (two states), we contrast gender monomorphic and gender dimorphic species. Gender monomorphic and gender dimorphic means that in one population there are one or two functional classes of sex [START_REF] Lloyd | Sexual strategies in plants III. A quantitative method for describing the gender of plants[END_REF]. Therefore bisexuality and monoecy are gender monomorphic and androdioecy, gynodioecy, and dioecy are gender dimorphic.

For each configuration, we used and compared parsimony, maximum likelihood (ML) and Bayesian approaches to reconstruct ancestral states [START_REF] Sauquet | Zygomorphy evolved from disymmetry in Fumarioideae (Papaveraceae, Ranunculales): new evidence from an expanded molecular phylogenetic framework[END_REF]Sauquet et al. 2017). Parsimony analyses were conducted in Mesquite v3.04 [START_REF] Maddison | Mesquite: a modular system for evolutionary analysis[END_REF], ML and Bayesian analyses were conducted in BayesTraits V2 [START_REF] Pagel | BayesTraits V2[END_REF]. Parsimony and ML analyses were conducted with the Maximum Clade Credibility (MCC) tree produced from the BEAST analysis, whereas Bayesian analyses were conducted with 3600 trees randomly sampled from the posterior of the BEAST analysis.

To test the influence of topological uncertainty on ancestral state reconstruction, we also conducted additional Bayesian analyses with a fixed (MCC) tree.

Maximum likelihood analyses presented here explored two models for each configuration: an equal-rates model (ER, or Mk1), assuming equal transition rates among all character states [START_REF] Lewis | A likelihood approach to estimating phylogeny from discrete morphological character data[END_REF], and an all-rates-different model (ARD), allowing distinct (asymmetric) transition rates among character states (Pagel 1994). The best model (equal or unequal rate) for each configuration was selected according to the Akaike Information Criterion (AIC). We calculated ΔAIC between two models for each configuration, using ΔAIC of 2 or more as a criterion for positive support of the best-fit model [START_REF] Posada | Model selection and model averaging in phylogenetics: advantages of Akaike information criterion and Bayesian approaches over likelihood ratio tests[END_REF]. Recognition of bisexuality and monoecy as gender monomorphic, and gynodioecy, androdioecy and dioecy as gender dimorphic [START_REF] Lloyd | Sexual strategies in plants III. A quantitative method for describing the gender of plants[END_REF] We also reconstructed ancestral breeding system states using both a "common" (i.e., fixed-model) Bayesian approach and a reversible-jump Markov chain Monte Carlo (rjMCMC) strategy. While both allowed us to take parameter and phylogenetic (incl.

molecular dating) uncertainty into account in ancestral state reconstruction (Pagel et al. 2004), the rjMCMC approach allowed us to explore and visit multiple models of morphological evolution in proportion to their posterior probabilities (Pagel and Meade 2006). Both the equal-rates (ER) and the all-rates-different (ARD) model were tested in the common Bayesian approach, and their relative fit was compared with the Bayes Factor [START_REF] Kass | Bayes Factors[END_REF] according to the criteria of [START_REF] Lodewyckx | A tutorial on Bayes factor estimation with the product space method[END_REF]. Chain lengths were set to ten million generations (or two million generations for fixed-tree analyses) and parameters and ancestral states were sampled every 1000 generations. Chain convergence was checked in Tracer v1.5 (Rambaut and Drummond 2009), with the first 10% generations excluded as burn-in.

To investigate the impact of different taxonomic sampling strategies (with or without outgroups of Ficus) on ancestral state reconstruction of breeding systems in Ficus, we also conducted analyses with only Ficus species (same topology). Here, after excluding non-Ficus species in configuration A, we used Mesquite instead of BayesTraits for ML reconstructions because we wished to test the impact of different root states priors on ancestral state reconstruction. We applied both the equal-rates (Mk1) and the asymmetric 2-rate (AsymmMk, equivalent to ARD for a binary character) models to reconstruct the ancestral state of Ficus. In addition, two different root states priors (equal and equilibrium frequencies) were applied for the asymmetric 2-rate (AsymmMk) model in Mesquite (whereas BayesTraits assumes an equal frequency root state prior for all states).

RESULTS

Phylogenetic reconstruction

Chloroplast markers were insufficient to resolve relationships within Ficus, as most internal nodes had low support (Fig. S1A-D). We found three instances of conflict (F.

trigonata clustered with F. tinctoria; F. religiosa clustered with F. benghalensis; F. pumila clustered with F. hirta, bootstrap probability 79, 80, and 57) between topologies reconstructed from nuclear and chloroplast sequences (Fig. S1 D andE). Except these, we found no well supported conflict between chloroplast and nuclear trees (Fig. S1). We thus combined the eight markers and focus here on the results from this combined analysis (Fig.

2 andS2). Hereafter we only discuss the trees reconstructed from the whole dataset. The topologies reconstructed with maximum likelihood (ML) and Bayesian approaches were consistent (Fig. S1 M and Fig. S2). Differences were only observed in some weakly supported nodes. All families, all tribes of Moraceae except Moreae, and all genera were supported as monophyletic. Tribe Moreae was reconstructed as polyphyletic because

Streblus smithii clustered with Maclureae (the other two species of Streblus sampled were instead found to be nested in Moreae; Fig. S2). In Ficus, subgenera Phamacosycea, Urostigma, and Ficus were found to be paraphyletic, and subgenus Synoecia nested in subgenus Ficus. In addition, some deep nodes (e.g., the most recent common ancestor of F. carica and Urostigma) were poorly supported in the Bayesian analysis (Fig. S2). With or without fossil age constraints, the topologies obtained with BEAST were identical and the posterior probability values of each nodes were similar as well (result not shown).

Excluding the crown-group calibration of Ficus or using an older maximum root age constraint did not influence the topology either (Table S4). (marked with orange star). For full details of this tree, see Fig. S2.

Molecular dating

After combination of the four separate runs, the effective sample size (ESS) of each parameter was over 100, most of them over 200. The crown-group ages of Ficus and Moraceae were estimated asEocene (40.6-55.9 Ma) and Late Cretaceous (73.2-84.7 Ma), respectively (Table 2; Fig. 2). The oldest tribe was Artocarpeae (64.0-68.6 Ma) and

Maclureae originated most recently (8.9-41.1 Ma). In Ficus, subgenus Sycomorus was the oldest (28.0-41.1Ma) (Table 2). Excluding the Ficus stem-group calibration or using an older maximum root age constraint did not have strong impact on estimated ages (Table S4).

Ancestral state reconstruction

Parsimony, maximum likelihood (ML) and Bayesian analyses of the full dataset (320 species) reconstruct dioecy as ancestral for Moraceae with strong support in all six character state configurations (Table 3, Fig. 3, Fig. S3 andS4). Furthermore, we find that dioecy has evolved into monoecy at least five times in Moraceae. However, several alternative scenarios are possible due to uncertainty in the most recent common ancestor of Dorstenieae, Castilleae, and Ficeae (Ficus) (Fig. 3). In one possible scenario, monoecy evolved independently once in Artocarpeae (Artocarpus); once in Moreae (Morus); once in Dorstenieae (with at least two subsequent reversals); twice in Castilleae (Perebea humilis; and the clade of Antiaris toxicaria and Mesogyne insignis); and at least once in Ficeae (Ficus). In another possible scenario, monoecy evolved independently in Artocarpeae and Moreae (as above), but shares a common origin in the clade of Dorstenieae, Castilleae, and Ficus, with at least two reversals (to dioecy) in Dorstenieae, and one in the ancestor of Castilleae (followed with two subsequent gains again in the tribe as above).

Within Ficus, no matter whether gynodioecy (functional dioecy) is ancestral (and therefore intermediate) or not (see below), at least three transitions from gynodioecy to monoecy are reconstructed (within section Sycomorus). Some of the transitions such as the direct transition between gynodioecy and dioecy in configurations A, B or D were estimated to be unlikely (i.e., transition rates were estimated to be zero) in the rjMCMC analyses (Fig 2). The 95% HPD of transition rates in each configuration broadly overlapped. In the Bayesian analysis with the unequal rate model in configuration A (Table S5), the chain did not converge. This remained true when fixing the tree (i.e. using the MCC tree) instead of using a sample of trees from the BEAST posterior. Reconstructions with only species of Ficus also showed an uncertain result with respect to the ancestral state of the genus (Fig. S5). Both parsimony and ML with the equalrate model strongly supported monoecy to be ancestral in Ficus, whereas the ML approach with the unequal rate model supported gynodioecy as ancestral in the genus. The support for ancestral gynodioecy differed with root state prior: with the equal frequencies root state prior, gynodioecy was highly supported (0.99 out of 1) whereas with the equilibrium frequencies root state prior, the probability for gynodioecy decreased sharply to 0.56, suggesting high uncertainty.

DISCUSSION

Phylogenetic relationships in Moraceae

The general topology of phylogenetic trees reconstructed in this study (both maximum likelihood and Bayesian approach) is consistent with previous studies in Ficus (Cruaud et al. 2012), Moraceae (Clement and Weiblen 2009) and Rosales (Zhang et al. 2011), except for relationships within Rhamnaceae, for which our results partly differ from those of two recent phylogenetic studies [START_REF] Onstein | Do Mediterraneantype ecosystems have a common history?-Insights from the buckthorn family (Rhamnaceae)[END_REF][START_REF] Hauenschild | Phylogenetic relationships within the cosmopolitan buckthorn family (Rhamnaceae) support the resurrection of Sarcomphalus and the description of Pseudoziziphus gen[END_REF]. This conflict may be due to the use of different molecular markers and to our sampling of this outgroup family with low density. Here we used matK, rbcL, and ndhF, while Hauenschild et al. (including all the markers we used and trnL-F, psbA, psbA-trnH, rpl16 and ITS). However, these differences are unlikely to affect our ancestral state reconstructions, given the distant positions of Rhamnaceae and Moraceae in the phylogeny of Rosales, and the fact that the breeding system state of most of the species in Rhamnaceae has been scored as missing data in this study. In our phylogenetic reconstruction, outgroup relationships are consistent with previous work [START_REF] Wang | Rosid radiation and the rapid rise of angiospermdominated forests[END_REF]Soltis et al. 2011;Zhang et al. 2011).

Within Moraceae, phylogenetic relationships among genera are consistent with recent studies (Zerega et al. 2005;Clement and Weiblen 2009;Zerega et al. 2010;Williams et al. 2017) except for Moreae (Fig. 2 and Fig. S2). Tribe Moreae was found to be paraphyletic due to the position of Streblus smithii as sister to Maclureae. This relationship is not strongly supported (posterior probability 0.76) and may be caused by too few informative sites for this species. In our dataset, S. smithii is only represented by the ndhF sequence and our separate analysis of the ndhF dataset, in which all three species of Streblus were sampled, supports the same result as our combined analysis (Fig. S1C). In the original paper where the sequence came from (Datwyler and Weiblen 2004), the species was found in an uncertain phylogenetic position. However, in a later reconstruction (Clement and Weiblen 2009) where ndhF, 26S, and morphological data were combined, S. smithii clustered with the other two species in the genus, possibly because the nuclear (26S) and morphological signal overcame a divergent ndhF sequence in their combined analysis.

Within Ficus, subgenera Phamacosycea, Urostigma and Ficus were found to be paraphyletic, and subgenus Synoecia is nested within subgenus Ficus and some deep nodes remain poorly supported. These results are consistent with the most recent comprehensive phylogenetic study of the genus (Cruaud et al. 2012), while they are in conflict with a recent phylogenomic study of Ficus based on full chloroplast genomes and this conflict could be caused by potential cyto-nuclear discordance (Bruun-Lund et al. 2017).

Chloroplast markers performed poorly in Ficus as the informative sites they provide are too few (Cruaud et al. 2012). Almost no branch length can be observed in the chloroplast tree for the shallow nodes (Fig. S1). In our analysis, Ficus species were represented not only by chloroplast markers but also by nuclear ones which were found informative enough in former studies [START_REF] Rønsted | Reconstructing the phylogeny of figs (Ficus, Moraceae) to reveal the history of the fig pollination mutualism[END_REF]Cruaud et al. 2012). In the Bayesian analyses, 3600 trees were taken into account, therefore the uncertainty of topology was considered.

A new time scale for Moraceae diversification

The Misiewicz and Zerega (2012) for Dorstenia (84.8-132.0 Ma).

BEAST analyses without internal calibrations showed identical topology and similar posterior support compared to our reference analysis using all calibrations, suggesting that here calibration did not influence the topology. Here, we used the maximum age estimated for crown-group Rosales by Magallón et al. (2015) as the maximum age constraint on the root in our analyses, but also tested the influence of older root constraints. Consistent with previous work [START_REF] Sauquet | Using fossils and molecular data to reveal the origins of the Cape proteas (subfamily Proteoideae)[END_REF]Massoni et al.2014;Foster et al. 2017), we found that alternative maximum root constraints had very little impact on estimated divergence times.

In addition, the crown-group age of Ficus was similar with or without the stem-group Ficus calibration (Table S4). These results suggest that our estimates are not an artifact of the calibration points near the nodes of interest.

Breeding system transitions in Moraceae and Ficus

It was found that increased taxon sampling density can be helpful in ancestral state reconstruction [START_REF] Salisbury | Ancestral state estimation and taxon sampling density[END_REF]. Ancestral states of breeding systems in Moraceae have been previously reconstructed with parsimony and a phylogenetic tree of 46 species [START_REF] Weiblen | Phylogenetic relationships of functionally dioecious Ficus (Moraceae) based on ribosomal DNA sequences and morphology[END_REF], 83 species (Datwyler and Weiblen 2004) and 73 species (Clement and Weiblen 2009). In this study, we reconstructed ancestral states with a more densely sampled phylogenetic tree (incl. 200 species of Ficus and 72 species from other genera of Moraceae), compared different model-based approaches, and took into account phylogenetic uncertainty in our Bayesian analyses. Here, we found similar results for two focal nodes (Moraceae and Ficus) whether using a single tree or thousands of trees sampled from the posterior of BEAST analyses, suggesting that, in our analyses, topology did not affect the estimation. In all the approaches and models, dioecy was strongly supported to be ancestral for Moraceae, while the ancestral state for Ficus remains unclear.

Dioecy has been reported to correlate with several ecological traits, including abiotic pollination, fleshy fruits, and woody growth form [START_REF] Freeman | Ecology of plant dioecy in the intermountain region of Western North America and California[END_REF][START_REF] Renner | Dioecy and its correlates in the flowering plants[END_REF][START_REF] Vamosi | Phylogenetic analysis of the ecological correlates of dioecy in angiosperms[END_REF]. During the history of Moraceae, dioecy evolved into monoecy at least five times (see above; Fig. 3, Fig. S3 andS4). Although a previous study suggested no strong statistical support for a relationship between breeding system and pollination syndrome in Moraceae (Clement and Weiblen 2009), our results reinforce the general pattern that would be consistent with such a relationship: dioecious taxa tend to be wind-pollinated, whereas monoecious taxa tend to be insect-pollinated. For instance, we infer Artocarpeae as ancestrally dioecious (Fig. 3, Fig. S3 and S4) and Datwyler and Weiblen (2004) reconstructed the clade as ancestrally anemophilous. In our reconstruction, dioecy has evolved to monoecy in Artocarpus, which is characterized by inflorescences with an insect pollination syndrome. This relationship may also apply to tribe Dorstenieae as a whole, most species of which are monoecious and were hypothesized to be insect pollinated according to their inflorescence structure (Berg and Hijman 1999). However, the relationship becomes less clear in Castilleae and Ficus (Datwyler & Weiblen 2004), partly because all species of Ficus are insect-pollinated (Cook and Rasplus 2003), yet only half of them are monoecious. Unfortunately, the current lack of sufficient data on actual pollination modes in Moraceae (outside Ficus) precludes us from further testing this potential correlation at the family level.

Although the ancestral state of Ficus proved to be particularly difficult to reconstruct in this study due to our exploration of various models and character state configurations, most of the evidence supports monoecy rather than gynodioecy as the ancestral state in the genus. As highlighted in recent studies of floral traits using the same methods [START_REF] Sauquet | Zygomorphy evolved from disymmetry in Fumarioideae (Papaveraceae, Ranunculales): new evidence from an expanded molecular phylogenetic framework[END_REF](Sauquet et al. , 2017)), our results suggest that great caution should be exercised when interpreting results from ML analyses exploring a limited set of models. Indeed, the models with highest posterior probability identified through our reversible-jump MCMC analyses corresponded to neither the equal rates or the unequal rates models exlpored with ML. These best-fit models typically involved only one free parameter, but excluded some transitions (Fig. 1).

In addition, Bayesian analyses present the advantage of taking into account phylogenetic and branch length uncertainty (Pagel et al. 2004). All our Bayesian analyses that took into account phylogenetic and branch length uncertainty suggested monoecy as ancestral in Ficus, whether models were allowed to vary (Table 4) or were fixed (Table S5), except for configuration F. In addition, the model-averaged rates from the rjMCMC analyses of the three configurations where dioecy and gynodioecy were treated as separate states (A, B, F)

were estimated to zero, suggesting that direct transitions between the two states do not occur in Moraceae. Therefore, we argue it is more probable that monoecy is ancestral in Ficus and represents an intermediate state between dioecy and gynodioecy in the genus (Fig. 3). From a functional point of view, it also seems more likely that the highly specialized pollination association between Ficus and fig wasps started in monoecious figs.

Assuming that monoecy is ancestral in Ficus, the phylogenetic distribution of monoecy and gynodioecy in the genus and the uncertainty remaining in the ancestral states of several deep nodes allow for various scenarios (Figs. 3, S3-S5). It is possible that gynodioecy evolved only once (after divergence of the Pharmacosycea 1 and Urostigma 1 lineages), followed by two main reversals to monoecy (Pharmacosycea 2 and Urostigma 2).

Alternatively, gynodioecy may have evolved twice independently (once in Sycomorus and once in the ancestral lineage of Sycidium, Ficus 1, and Synoecia). In all scenarios, gynodioecy reverted to monoecy at least three times within Sycomorus (Fig. 3). The biological reasons for these fluctuations between monoecy and gynodioecy in Ficus remain unclear, but several hypotheses have been proposed. Gynodioecy in Ficus may be linked with an adaptation to seasonality [START_REF] Kjellberg | The stability of the symbiosis between dioecious figs and their pollinators: a study of Ficus carica L. and Blastophaga psenes L[END_REF]; reduction of non-pollinating wasps [START_REF] Kerdelhué | The evolution of dioecy among Ficus (Moraceae): an alternative hypothesis involving non-pollinating fig wasp pressure on the fig-pollinator mutualism[END_REF]; persistence of pollinator populations when the host fig populations are small [START_REF] Kameyama | Persistence of a fig wasp population and evolution of dioecy in figs: a simulation study[END_REF]; seed protection [START_REF] Greeff | Can seed protection lead to dioecy in Ficus?[END_REF]; disperser selection; de-coupling of wasp and seed size; survival from predation;

non-pollinators and predator satiation; and chronic pollinator shortages, crop asynchrony, and inbreeding depression (Harrison and Yamamura 2003). In addition, breeding system evolution in Ficus may have a relationship with biogeographic distribution. For instance, there appear to be no gynodioecious species of figs in South America (Cruaud et al. 2012).

Different climates and habitats in different continents may have accelerated breeding system evolution in Ficus. The genus has been inferred to have originated in Eurasia, then migrated and diversified in South America during the Miocene (Cruaud et al. 2012).

Migration to a new continent may have led to breeding system reversal from gynodioecy to monoecy in the clade Urostigma 2 (Fig. 3, Fig. S3 and S4) under the scenario where gynodioecy originated only once in Ficus.

Dioecy not an evolutionary dead end in Moraceae

Our results provide additional strong support for dioecy to be ancestral in Moraceae as a whole, regardless of the approach or model used for reconstructing ancestral states (Fig3, Table 3 and4), consistent with previous work [START_REF] Weiblen | Phylogenetic relationships of functionally dioecious Ficus (Moraceae) based on ribosomal DNA sequences and morphology[END_REF]Datwyler and Weiblen 2004;Clement and Weiblen 2009). Furthermore, our extensive sample of outgroups shows that dioecy in Moraceae is ancestral to a larger clade including at least Cannabaceae, Urticaceae, and Moraceae (Fig. 3).

These results add to a growing number of studies that have challenged the notion that dioecy is an irreversible trait in flowering plants. Dioecy was once suggested to be evolutionary dead end (Bull and Charnov 1985). It is found in only approximately 6% of angiosperms and distributed widely in around 43% families of flowering plants (Renner 2014). However, "evolutionary dead end" has been used to describe different patterns of macroevolution [START_REF] Bromham | Detecting macroevolutionary self-destruction from phylogenies[END_REF], including: 1) irreversible evolution (Bull and Charnov 1985;[START_REF] Glémin | Mating systems and selection efficacy: a test using chloroplastic sequence data in angiosperms[END_REF]; or 2) a state that presents more long-term disadvantages than advantages compared to other states, ultimately leading to probable extinction (Käfer et al. 2014;Käfer et al. 2017). Dioecy used to be thought as irreversible (Bull and Charnov 1985) and high extinction rates were found on dioecious clades by comparative approach (Heilbuth 2000). However, it has been known for a long time that although dioecious angiosperms are exposed to a higher risk of failure to find a sex partner, they are also advantageous in having obligate cross-fertilization (Darwin 1876). Various studies have shown that dioecious species may include individuals with other breeding systems and that dioecy in general may be more labile than previously thought [START_REF] Pannell | Widespread functional androdioecy in Mercurialis annua L. (Euphorbiaceae)[END_REF][START_REF] Korpelainen | Labile sex expression in plants[END_REF][START_REF] Ainsworth | Boys and girls come out to play: The molecular biology of dioecious plants[END_REF]Käfer et al. 2017). These transitions could be linked to ecological processes including changing plant-pollinator relationships and environments [START_REF] Case | A phylogenetic study of evolutionary transitions in sexual systems in Australasian Wurmbea (Colchicaceae)[END_REF]; long distance dispersal [START_REF] Schaefer | A three-genome phylogeny of Momordica (Cucurbitaceae) suggests seven returns from dioecy to monoecy and recent long-distance dispersal to Asia[END_REF]; hybridization and polyploidy [START_REF] Obbard | Hybridization, polyploidy, and the evolution of sexual systems in Mercurialis (Euphorbiaceae)[END_REF]). In addition, various studies have documented clades where dioecy is the ancestral rather than the derived state, challenging

the assumption of irreversibility of dioecy evolution. For instance, dioecy has been reconstructed as ancestral in Myristicaceae, with at least four subsequent shifts to monoecy in the family [START_REF] Sauquet | Androecium diversity and evolution in Myristicaceae (Magnoliales), with a description of a new Malagasy genus, Doyleanthus gen. nov[END_REF]) and appears to be a transition state between monoecy and polygamy in Arecaceae [START_REF] Nadot | The palm family (Arecaceae): a microcosm of sexual system evolution[END_REF]. A recent analysis of macroevolutionary dynamics of breeding system in 40 angiosperm genera supported transitions away from dioecy in Bursera (Burseraceae) and Dodonaea (Sapindaceae; Goldberg et al. 2017).

Furthermore, two recent angiosperm-wide analyses showed that dioecy does not have a negative relationship with diversity: dioecy was shown to correlate with an increased diversification rate (Käfer et al. 2014) or not to have a significant relationship with the latter (Sabath et al. 2016).

Conclusion

Here we reconstructed divergence times and ancestral breeding systems in Moraceae.

With 12 fossil calibrations, the crown group ages of Ficus and Moraceae were estimated as 40.6-55.9 Ma (Eocene) and 73.2-84.7 Ma (Late Cretaceous). Dioecy was supported as the ancestral state of Moraceae and we showed that it evolved into monoecy at least five times in the family, and subsequently reversed in some clades (e.g. Artocarpeae, Dorstenieae). In Ficus, which represents 75% of species in Moraceae, the ancestral state was estimated to be monoecy with moderate support and at least one transition from monoecy to gynodioecy followed by at least three reversals were estimated. In future work, to investigate the exact breeding system evolutionary scenario of Ficus, it would be important to sample basal lineages of Ficus more densely and improve phylogenetic resolution of the backbone of the genus. The new time scale of Moraceae and Ficus we provide here will also be useful for future analyses of biogeography and co-diversification in the family. Finally, our results lend further support to the growing idea that dioecy does not always represent an evolutionary dead end, shedding light on the understanding of breeding system evolution in angiosperms more generally. group age of Dorstenia has been estimated using a range of taxonomic sampling density to be 3.5-18.4 Ma (Zerega et al., 2005, with The pantropical distribution, almost equal diversity on both sides of the Atlantic, and the diverse traits of Dorstenieae make it a good model for understanding the origin and evolution of Neotropical biodiversity and the relationships between African and Neotropical floras.

Genomic targeted enrichment approaches have been shown to be more efficient and economic than Sanger sequencing [START_REF] Lemmon | High-throughput genomic data in systematics and phylogenetics[END_REF]McKain et al. 2018) and have been widely used in phylogenetic studies in recent years [START_REF] Xi | Coalescent versus concatenation methods and the placement of Amborella as sister to water lilies[END_REF]Fisher et al. 2016;Hart et al. 2016;[START_REF] Mitchell | Anchored phylogenomics improves the resolution of evolutionary relationships in the rapid radiation of Protea L[END_REF]Couvreur et al. 2019). In this study, we targeted nuclear genes as they have been suggested to hold the greatest potential for investigating the evolutionary history of angiosperms for several reasons.

Firstly, nuclear genes have worked well in reconstructing more strongly supported phylogenetic relationships than organellar markers in both deep and shallow time scales [START_REF] Xi | Coalescent versus concatenation methods and the placement of Amborella as sister to water lilies[END_REF][START_REF] Mitchell | Anchored phylogenomics improves the resolution of evolutionary relationships in the rapid radiation of Protea L[END_REF]. Secondly, nuclear genes are assumed to be unlinked, decreasing the probability of misleading phylogeny by part of the genes used (Fisher et al. 2016). Lastly, nuclear genes present multiple lineage histories, contrary to plastid genes, which are usually considered to represent a single locus.

In this study, we included samples representing all fifteen genera of Dorstenieae (Table 1) to reconstruct a comprehensive dated phylogenetic tree of Dorstenieae with two primary goals: 1) to test phylogenetic relationships and monophyly of Dorstenieae genera, and relationships within Dorstenia; and 2) to investigate divergence times and the biogeographic history of tribe Dorstenieae and genus Dorstenia. This study also provided an opportunity to test the potential of the targeted enrichment strategy to resolve species-level phylogenetic relationships using herbarium material, taking advantage of a recently developed nuclear bait set for Moraceae [START_REF] Gardner | Low-coverage, whole-genome sequencing of Artocarpus camansi (Moraceae) for phylogenetic marker development and gene discovery[END_REF], Johnson et al., 2016).

MATERIALS AND METHODS

Specimen and sample collection

In this study, we follow the classification of Clement & Weiblen (2009) 1). We included 83 species (93 taxa) representing all currently recognized Dorstenieae genera and 53% of the species within tribe Dorstenieae (Table S1, we also collected samples of 22 more species in Dorstenieae and they were finally excluded in the main analyses, see below). Additionally, five outgroup taxa in Moraceae were included. They represent five of the six other Moraceae tribes: Artocarpeae, Castilleae, Ficeae, Moreae, and the newly created Parartocarpeae (Zerega and Gardner 2019). Most taxa were extracted from herbarium material sampled from the Field Museum (F), the Missouri Botanical Garden (MO), the New York Botanical Garden (NY), and the Muséum national d'Histoire naturelle (P) (Table S1). Specimens collected within the last twenty years, and reliably identified material with inflorescences or infructescences were preferred. Samples include 55 (49%) species of Dorstenia, representing eight out of the nine sections proposed by Berg and Hijman (1999). Section Bazzemia, which contains only one species from Mozambique (Berg & Hijman, 1999), was not sampled.

DNA extraction and sequencing

Whole genome DNA was extracted using a modified CTAB method (Doyle & Doyle, 1987). Extracted DNA was re-suspended in 50 μl light TE. The DNA concentration of each sample were measured using Qubit 2.0 fluorometer (Life Technologies, Carlsbad, CA, USA) following the standard protocol. For library preparation, we used 200 ng DNA when possible, but never less than 20 ng. All samples were run on an agarose gel to assess fragment size. Samples with DNA fragment lengths longer than 500 base pairs (bp) were sonicated on a Covaris M220 (Covaris, Wobum, MA, USA) for 45 seconds at 50 W peak power and a duty factor of 20%, which typically produces average fragment sizes of 550 bp. We prepared dual-indexed sequencing libraries using the KAPA Hyper Prep Library Construction Protocol (KAPA biosystems, Wilmington, MA, USA), generally following the manufacturer's protocol except that most steps were performed at one-quarter volume to save costs. Low-input or very degraded samples were not size selected. Libraries were amplified using 12 cycles of PCR, but half of the unamplified template was retained in case PCR needed to be repeated. Products were cleaned using Solid Phase Reversible Immobilization (SPRI)

and quantified using a high-sensitivity dsDNA assay on a Qubit 2.0. For samples with a concentration less than 5 μg/ml, we repeated PCR amplification with 14 cycles.

Successful libraries were combined into seven pools of 13 to 14 libraries each. We hybridized the libraries to custom Moraceae probes [START_REF] Gardner | Low-coverage, whole-genome sequencing of Artocarpus camansi (Moraceae) for phylogenetic marker development and gene discovery[END_REF] manufacturered by Arbor Biosciences (Ann Arbor, MI, USA) as a MYbaits kit.

Hybridization for 20 hours followed the manufacturer's protocol, and products were reamplified 14 cycles of PCR. Amplified products were quantified on a Qubit 2.0, and fragment sizes were determined using a High-Sensitivity DNA assay on a BioAnalyzer 2100 (Agilent Technologies, Palo Alto, Califonia, USA). When required, adapter dimer was removed using 0.7x SPRI beads. All libraries were then sequenced on a single lane of an Illumina HiSeq 2000 (2 × 100 bp, paired-end) by Genewiz (Genewiz, South Plainfield, NJ, USA).

To the 96 libraries prepared for this study, we added sequences from 14 samples Gardner (Zerega and Gardner 2019). These 20 samples were all sequenced on an Illumina MiSeq with somewhat longer reads (2x300bp, v3). Finally, we used transcriptomic reads for Broussonetia papyrifera obtained from GenBank.

Sequence cleaning, assembly, and filtering

Sequences were assembled using HybPiper [START_REF] Johnson | HybPiper: extracting coding sequence and introns for phylogenetics from high-throughput sequencing reads using target enrichment[END_REF] are predicted using Exonerate [START_REF] Slater | Automated generation of heuristics for biological sequence comparison[END_REF]. When a gene is assembled into several disconnected contigs-common in degraded samples where the fragments (and therefore effective read length) are very short-HybPiper scaffolds these short contigs in the correct order based on the reference. In the event that multiple genes are assembled for a single target, HybPiper distinguishes orthologs from paralogs using a combination of alignment length and identity relative to the reference. HybPiper outputs include in-frame CDS sequences as well as "supercontigs," which contain CDS as well as any flanking non-coding sequences. For our analyses, we used only the supercontig sequences. We then filtered the sequences within loci to remove those less than 150 bp long or shorter 25% of the average length for the locus.

After filtering, the number of genes recovered for each taxon varied widely among taxa (from zero for D. scaphigera and 515 for Ficus macrophylla) and the number of taxa retrieved for each gene differed sharply among genes (from two to 101). These extreme differences would result in dataset with a high proportion of missing data, which could impact the accuracy of both phylogenetic reconstruction and divergence time estimate. Therefore, we filtered the dataset further by selecting 102 genes with high taxon occupancy and excluding taxa for which less than 30 of these genes were recovered (Table S1).

Phylogenomic reconstruction and molecular dating

Sequences were aligned with MACSE [START_REF] Ranwez | MACSE: multiple alignment of coding sequences accounting for frameshifts and stop codons[END_REF], gene by gene (-fs 30 -stop 50). We reconstructed phylogenetic relationships using both a concatenated maximum likelihood (ML) and a coalescent approach. Milicia exceisa and Artocarpus venenosus were specified as the most external outgroups based on recent phylogenetic analyses of Moraceae (Clement and Weiblen 2009;Zerega et al. 2010;[START_REF] Chung | Molecular recircumscription of Broussonetia (Moraceae) and the identity and taxonomic status of B. kaempferi var. australis[END_REF]Gardner et al. 2017;Zhang et al. 2019). Sequences were concatenated into a supermatrix and partitioned according to genes using the fasta_merge.py script from HybPiper; a ML tree was reconstructed with RAxML v8.2.10 (Stamatakis 2014) as implemented on CIPRES (Miller et al. 2010) 1000 rapid bootstrap replicates.

GTRGAMMA was chosen as the substitution model. For the coalescent approach, we first reconstructed ML trees with RAxML v8.2.10 for each gene, using with GTRGAMMA model and 200 rapid bootstrap replicates. These ML trees for each gene were used to infer a species tree using the summary coalescent approach implemented in ASTRAL-II [START_REF] Mirarab | ASTRAL-II: coalescent-based species tree estimation with many hundreds of taxa and thousands of genes[END_REF]. Node support was calculated based on 200 bootstrap replicates, with resampling within loci.

Because the topologies of the ML and species trees were broadly consistent, we used the ML tree to estimate divergence times with a relaxed clock model and four and then used them to estimate the divergence times. We set the prior in both steps according to the estimates of baseml. In both steps, we ran the process for 38.5 million generations, with the first 10% of the chain length discarded as burnin, sampling a total of 10,000 generations at a frequency of once every 3,500 generations. Two independent runs with the same settings were conducted to confirm the convergence of the MCMC.

To check the influence of the prior on the estimation, we used another prior setting (program defaults), followed the same steps as above, then ran the chain for 22 million generations, sampling a total of 10,000 generations at a frequency of once every 2,000 generations. After checking the convergence of the runs, we combined the results of them for each prior setting. Two additional independent runs with the two different prior settings, but without data, were also conducted to test the impact of the priors on the results. The convergence of Markov chain Monte Carlo (MCMC) was checked by reading the result log files in Tracer v1.7 [START_REF] Rambaut | Posterior summarisation in Bayesian phylogenetics using Tracer 1.7[END_REF]. For the Bayesian approach, divergence was checked by confirming that the effective sample size (ESS) of parameters in the two independent runs was over 100 after removing the first 10% of chain length as burnin.

Reconstruction of biogeographic history

Distribution data for the taxa sampled were collected from monographs and revisions (Rohwer and Berg 1993;Berg and Hijman 1999;Berg 2001;[START_REF] Chung | Molecular recircumscription of Broussonetia (Moraceae) and the identity and taxonomic status of B. kaempferi var. australis[END_REF]. We separated the distribution area of extant species of Dorstenieae into six areas (Fig. 1) following a recent study of Annonaceae, which share a similar global distribution (Couvreur et al. 2011): A, South America; B, North/Central America; C, Africa; D, Madagascar; E, India and Sri Lanka; F, Southeast Asia and Oceania. We decided to combine the two areas west and east of Wallace's line (areas F and G of Couvreur et al., 2011) into a single area (area F) because all species of Dorstenieae distributed east of Wallace's line were also observed west of Wallace's line in our dataset (e.g., Allaeanthus luzonicus, Fatoua villosa).

Only Dorstenieae (93 taxa, 83 species) were kept in the biogeographic reconstruction to avoid the bias of incompletely sampled outgroups. BioGeoBEARS (Matzke 2013) in R v3.5.1 (R Core Team, 2018) was used to reconstruct the biogeographic history of the clade with the DEC and DEC+J models. Prior studies have suggested that the DEC+J model can exacerbate the bias of preferring cladogenetic events (i.e., sympatry, vicariance and founder event speciation) over anagenetic processes in the DEC model, implying that the higher likelihoods typically obtained with DEC+J do not necessarily mean a better fit to the data than DEC (Ree and Sanmartín 2018). Therefore, we employed both the DEC and the DEC+J models. In addition, we ran the analyses with the DIVA-like and BAYAREA-like models including or not including founder-event speciation (+J). All analyses were run both with a simple dispersal matrix and a time-stratified model, hereafter referred to as model 0 and model 1, respectively. Both models were inspired by Couvreur et al. (2011) and are constrained by unequal dispersal relative rates aimed at reflecting the connectivity of biogeographic areas (Fig. 1). In the time-stratified model, the constraints on dispersal probabilities vary through five time periods (0-5, 5-30, 30-45, 45-65, and 65-75 Ma) according to physical distance between areas (Fig. 1). For instance, the dispersal rates from continental Africa (C) and Madagascar (D) to India and Sri Lanka (E) were constrained to be very low during the rafting of the Indian plate 30-65 Ma. Because all sampled species of Dorstenieae currently occupy no more than two areas, the maximum number of co-occurrence areas was set as three in all analyses.

RESULTS

Sequencing

Sequencing and assembly statistics appear in Table S1. Enrichment ranged from 0.34% to 84% reads on target, with the enrichment found in Artocarpeae and Moreae (58% and 84% on target, respectively) and the least efficient found in the genus Dorstenia (0.34% on target). The final dataset consisted of 98 taxa (89 species, spanning all of the currently recognized genera) and 102 genes. The width of the aligned supermatrix was 132,753 bp with 29.51% gaps and missing data.

Phylogenetic relationships

The maximum likelihood (ML) and coalescent analyses produced broadly identical topologies, with some weakly supported differences at shallow phylogenetic depths (Fig. 2, Figure S1a,b). Most of the nodes on the tree were strongly supported (over 90% bootstrap support , Fig. 2). Fatoua villosa was found to be sister to the remaining of Dorstenieae. Malaisia scandens and Broussonetia papyrifera together formed a clade that is sister to Allaeanthus. All genera of Dorstenieae sampled more than once were found to be monophyletic with two notable exceptions. Trymatococcus oligandrus and

Helianthostylis sprucei appear to be nested in Brosimum, and Scyphosyce (two species)

and Utsetela gabonensis were found to be nested among early-diverging lineages of Dorstenia. Hereafter we refer to the clade of Brosimum, Trymatococcus and

Divergence times of Dorstenieae and Dorstenia

The S2a,b). Divergence times estimated with the PL approach were all compatible with those from the Bayesian approach (i.e., falling within the 95% credibility intervals) except the age of crowngroup Broussonetia s.l., which was significantly older with PL (Table 2, Figure S2c).

Biogeographic history of Dorstenieae

We analyzed the dataset with both the dispersal-extinction-cladogenesis (DEC) and DEC+J (DEC with founder-event speciation) models. The time-stratified model fit the data significantly better than the time-constant model in all analyses conducted (Table S2), therefore we focus mainly on the results from the time-stratified analyses here, unless otherwise mentioned. In addition, the DEC+J model had a lower value of Akaike information criterion (AIC) than that of DEC (Table S2). The ancestral distribution area of Dorstenieae was estimated to be the combined area of continental Africa, Madagascar, and Southeast Asia and Oceania (CDF) with both the DEC and DEC+J models (Fig. 3, Fig. S3c). The ancestral area for both the stem and crown-group nodes of Dorstenia s.l. (incl. Scyphosyce and Utsetela) were estimated in continental Africa (C) with both models. South America (A) was found as the ancestral area of both the stem and crown-group nodes of the Neotropical Dorstenia clade with the DEC+J model of alternative ancestral areas (for full details, see Figure S3d).

(Fig. 3), while the DEC model reconstructed the combined area of South America and continental Africa (AC) as ancestral for the stem-group node of this clade (Figure S3c).

Dorstenia indica, endemic to India and Sri Lanka, was found to be nested in an African clade and diverged from its sister group in the Eocene to .

The Central/North American Dorstenia species clustered into two clades with strong support (Fig. 2), suggesting two independent colonizations from South America to Central/North America during the Oligocene to Miocene (Table 2). The ancestral states of the crown-group nodes of both Central/North American clades were estimated as

Central/North American. The stem-group nodes of clade 1 (from D. erythrandra to D. caimitensis) and clade 2 (from D. excentrica to D. contrajerva) were estimated to be Central/North American and South American, respectively (Fig. 3). The stem and crown-group nodes of Brosimum s.l. were estimated in the joint area of South America and continental Africa (AC), and in South America (A), respectively, with both the DEC and DEC+J models.

The additional biogeographic models produced similar results to those of the DEC and DEC+J models, with a few exceptions (Figure S3, Table S3). Reconstruction with time-constant models showed similar results as with time-stratified models. The ancestral distribution area for the crown-group node of Dorstenieae was estimated to be Southeast Asia and Oceania (F) by all the models based on the DIVA-like model (Figure S3i-l). In the time-stratified BayArea-like models, with or without founderevent speciation, the ancestral area of the crown-group node of Dorstenieae was estimated as the joint area of continental Africa and Southeast Asia and Oceania (CF), or as Southeast Asia and Oceania (F), respectively. The same pattern was found in reconstructions with time-constant BayArea-like models (Figure S3e-h). The stem and crown-group nodes of Dorstenia were estimated in continental Africa (C) in all the models. The crown-group node of Neotropical Dorstenia was estimated in South America (A) in all analyses, while the stem-group node was different among models.

All the BayArea-like models estimated this node in continental Africa (C) (Figure S3eh). The time-stratified DIVA-like models estimated the joint area of South America and continental Africa (AC), or South America (A) alone as the ancestral area of this node, with or without founder-event speciation respectively. The same results were found in time-constant DIVA-like models (Figure S3i-l). Lastly, we also ran all biogeographic analyses with the chronograms reconstructed with the penalized likelihood approach.

The results were similar with some exceptions (results not shown). For instance, Southeast Asia and Oceania (F) was estimated as the ancestral area of the crown-group node of Dorstenieae with the time-stratified DEC+J model. The stem-group node of Neotropical Dorstenia was estimated as continental Africa (C) by the same model.

DISCUSSION

Success of the targeted enrichment strategy with herbarium specimens

Most (99%) of the samples in this study were from herbarium specimens (Table S1). Some of them were collected more than 40 years ago and the amount of sample collections from the herbarium was typically limited (around 3 to 20 mg) due to destructive sampling policies. Three samples (Dorstenia aristeguietae, D. choconiana, and D. prorepens) were filtered by HybPiper because of the low matching of reads to the reference. We excluded another fifteen taxa for the low number of genes recovered (less than 30, Table S1). The lowest amount of DNA used among the final 98 taxa retained in our dataset was 25.1 ng (Dorstenia brasiliensis). Samples that were excluded ranged from being over 100 years old to 11 years old, while samples that were included were collected as long ago as 1923. These results suggest that the degradation of DNA in old herbarium specimens appears to have had little influence on this study.

This may be because of the short DNA fragments (on average less than 500 bp long) needed in the library preparation, those longer than this size requiring sonication. Thus, our results suggest that the targeted enrichment strategy and HybPiper pipeline [START_REF] Johnson | HybPiper: extracting coding sequence and introns for phylogenetics from high-throughput sequencing reads using target enrichment[END_REF] worked well with a broad range of ages of herbarium specimens.

Sequencing herbarium samples is a valuable approach for phylogenetics, as many species can be difficult to collect. Herbarium specimens often represent reliable and accurate vouchers of species identification, and some species may be rare or have even gone extinct in the wild [START_REF] Särkinen | How to open the treasure chest? Optimising DNA extraction from herbarium specimens[END_REF][START_REF] Staats | Genomic treasure troves: complete genome sequencing of herbarium and insect museum specimens[END_REF]. The success of similar targeted enrichment strategies with historical specimens for phylogenetic studies has previously been highlighted in other lineages of angiosperms at various scales, including Arabidopsis thaliana (Brassicaceae, [START_REF] Staats | Genomic treasure troves: complete genome sequencing of herbarium and insect museum specimens[END_REF], Inga (Fabaceae, Hart et al., 2016), and Annonaceae.

Our results also suggest that the baits, which were originally designed for tribes Artocarpeae and Moreae in the same family [START_REF] Gardner | Low-coverage, whole-genome sequencing of Artocarpus camansi (Moraceae) for phylogenetic marker development and gene discovery[END_REF], worked well in Dorstenieae. In addition to the phylogenetic markers developed from 333 inferred single-copy exons for Moraceae, we retrieved approximately another 200 genes in this study for Dorstenieae. Those assembled untargeted genes had a mean identity to at least one target gene of 84%, and a mean alignment length to at least one target gene of 87%

(as a percentage of the untargeted genes). These Moraceae specific baits worked well throughout the entire Moraceae family (Zerega and Gardner 2019). Five untargeted genes were included in our final 102-gene dataset. The most genes retrieved for one taxon was 515 from the outgroup sample Ficus macrophylla, suggesting a high probability that these baits would also work well to explore relationships within Ficeae (Ficus), the largest tribe in Moraceae (Couvreur et al. 2019).

Phylogenetic relationships in Dorstenieae and Dorstenia

With all the extant genera of Dorstenieae included, this is the most densely sampled phylogenetic study of Dorstenieae to date. The reconstructed relationships (Fig. 2) are generally consistent with previous work (Zerega et al. 2005;Clement and Weiblen 2009;Misiewicz and Zerega 2012;[START_REF] Chung | Molecular recircumscription of Broussonetia (Moraceae) and the identity and taxonomic status of B. kaempferi var. australis[END_REF]Zhang et al. 2019), but with stronger support, especially in Brosimum s.l. and Dorstenia s.l..

Prior to this study, the most densely sampled molecular phylogenetic analysis of Dorstenia was provided by Misiewicz and Zerega (2012), based on ITS sequences of 35 taxa (32 species) of Dorstenia and seven outgroup species. Our results are similar to those of Misiewicz and Zerega (2012) with respect to shallow-level relationships.

However, they differ markedly at a deeper level in that the authors had found three African species (D. variifolia, D. tayloriana var. tayloriana, and D. cuspidata) to be nested in the Neotropical clade of Dorstenia. In our analyses, D. variifolia and D.

tayloriana were also sampled and they were sister species as in Misiewicz & Zerega (2012), but we found these species nested in African clades with strong support. The difference may be explained by the root setting methods, the variation in each study of both number of genes and density of taxon sampling for both Dorstenia species as well as of non-Dorstenia species within the tribe Dorstenieae. Most differences between the ML and coalescent trees reconstructed in this study concentrated in the Neotropical Dorstenia clade (Fig. 2). This result and the short branches observed in this clade suggest incomplete lineage sorting in the diversification of Neotropical Dorstenia species [START_REF] Pamilo | Relationships between gene trees and species trees[END_REF].

Divergence times of Dorstenieae and its genera

Runs with or without data showed different results suggesting the data were informative. Runs with different prior settings showed similar estimates, suggesting that our results are robust to various assumptions on rate variation. The crown-group age of Dorstenieae was estimated in the Upper Cretaceous (65.8-79.8 Ma; Table 2), which overlaps with former studies that used fewer genes (Zerega et al. 2005;Gardner et al. 2017;Zhang et al. 2019). The crown-group age of Dorstenia was estimated in the Paleocene (50.8-67.3 Ma), which is younger than in Misiewicz and Zerega (2012).

Although the estimated stem-group node of Neotropical Dorstenia clade in our study was younger than that in Misiewicz and Zerega (2012), the crown-group node of the same clade fell within a similar range in both studies. The difference in the stem-group node may be caused by different topologies and number of genes used to estimate the divergence time. The crown-group age of Brosimum s.l. and Broussonetia s.l. were estimated in the period from Oligocene to Miocene and the Eocene (14.9-31.1 respectively), during which time the whole earth cooled down from Mid-Eocene Climatic Optimum and was warmer than the current climate [START_REF] Zachos | An early Cenozoic perspective on greenhouse warming and carbon-cycle dynamics[END_REF].

Biogeographic history of Dorstenieae

Our results suggest that the most recent common ancestor of Dorstenieae was widely distributed in the joint area of continental Africa, Madagascar and Southeast result is plausible or an example of pathological behavior of DEC models (Ree and Sanmartín 2018) will require further investigation.

Two independent dispersals from South America to Central/North America in Dorstenia were estimated from the Oligocene to Miocene (period from Node II to stemgroup node of Clade 1, 16.7-30.6 Ma and period from stem to crown group of Clade 2, 12.0-34.7 Ma), during which time waves of dispersal from South to North America have been found in other lineages (Bacon et al. 2015). This result supports the dispersal of plants from South America to Central/North America in the Neogene [START_REF] Bagley | Phylogeography and biogeography of the lower Central American Neotropics: diversification between two continents and between two seas[END_REF]). The branch length or time between the stem and crown group of Neotropical Dorstenia was very short (less than 10 Ma, Table 2), suggesting the rapid divergence of the clade. The branch length between the stem and crown group of the other Neotropical Dorstenieae lineage (Brosimum s.l.) was also short (ca. 12 Ma, Table 2). Therefore, LDD followed by rapid diversification would explain the extant distribution pattern of Neotropical Dorstenieae species. A similar pathway was found in the pantropically distributed tribe Annoneae (Annonaceae) [START_REF] Thomas | Historical biogeography of Goniothalamus and Annonaceae tribe Annoneae: dispersal-vicariance patterns in tropical Asia and intercontinental tropical disjunctions revisited[END_REF]Williams et al. 2017).

A diversity of seed dispersal modes has been reported in Dorstenieae, including autochory by expulsion or ejection of endocarp in Dorstenia, Bleekrodea, Fatoua, and zoochory in Brosimum lactescens and Trymatococcus amazonicus (Berg, 2001).

Dorstenia was suggested to be poorly adapted for LDD (Berg and Hijman 1999). In addition, small Dorstenia seeds of forest undergrowth species often germinate shortly after maturity, further reducing the chances of LDD (Berg 2001). Low probability of LDD may explain the single origin of Neotropical Dorstenia and the monophyly of the two Central/North America clades. Furthermore, our results suggest a comparatively faster succession of speciation events following establishment of Dorstenia in the Neotropics, suggesting rapid speciation after LDD has been an important process in shaping the origin of Neotropical diversity.

Ree and Sanmartín (2018) recently raised methodological concerns with models including both anagenetic and cladogenetic processes, especially the DEC and DEC+J models. In the main analyses presented here, we focused on the results of DEC-based models (i.e., DEC m0, DEC m1, DEC+J m0, DEC+J m1). We also did reconstructions with DIVA-like and BayArea-like based models. The difference among DEC, DIVAlike and BayArea-like models lie mainly in the cladogenetic process they assume: DEC and DIVA models explore both sympatric speciation and vicariance processes while BayArea only explores sympatric speciation [START_REF] Nicholas | Probabilistic historical biogeography: new models for founder-event speciation, imperfect detection, and fossils allow improved accuracy and model-testing[END_REF]. We compared the models with the Akaike information criterion (AIC). Whether or not founder-event speciation was allowed, led the BayArea model to rank from the best to the worst model (Table S2) in both time-stratified and time-constant analyses. AIC has been argued not to be a good criterion for models including both cladogenetic and anagenetic events due to a bias of the likelihood to favor time-independent cladogenetic processes, a problem exacerbated with the introduction of founder event-speciation in the model (Ree and Sanmartín 2018). Zero-estimate for dispersal and strong counterintuitive unparsimonious reconstruction may be signals for this bias. These two phenomena were not observed in our reconstruction (Table S2, Figure S3). We did not rely on model selection in this study. Instead, we emphasize that all of the models we used (incl. DEC, DIVA-like and BayArea-like based models) led to very similar reconstructions (Fig. 3, Fig. S3).

Taxonomic implications

Some of the currently recognized genera of Dorstenieae and all sections of Dorstenia may need to be modified based on the results from our analyses. Using 102 genes (132,753 bp), we obtained the same relationships among Malaisia, Broussonetia

and Allaeanthus as a previous study based on one chloroplast and one nuclear gene [START_REF] Chung | Molecular recircumscription of Broussonetia (Moraceae) and the identity and taxonomic status of B. kaempferi var. australis[END_REF], providing additional support for the recognition of Allaeanthus as a separate genus.

Both Scyphosyce and Utsetela were found to be nested in Dorstenia. Scyphosyce is a genus of two species from western Africa (Berg 1977). We sampled both species of Scyphosyce in our analyses and found them sister to each other, suggesting their placement within Dorstenia was unlikely to be caused by misidentification. The basal grade of the Dorstenia s.l. clade was formed by species of sections Nothodorstenia and Xylodorstenia of the genus Dorstenia and by the genera Scyphosyce and Utsetela, all of which share woody habit and larger seeds, which are referred to by Berg and Hijman (1999) as macrosperms. This basal grade shares other characteristics as well. The inflorescences of most Dorstenieae genera are bisexual (some species of Broussonetia, Allaeanthus, and Fatoua have unisexual inflorescences). The macrospermous species commonly have only one to a few pistillate flowers per inflorescence which produce one to few large seeds per infructescence). The remaining species of Dorstenia are herbaceous and have several to numerous pistillate flowers per inflorescence which can produce numerous smaller seeds (Berg and Hijman 1999). Utsetela is a genus of two species from western Africa (Berg 1977;[START_REF] Jongkind | Novitates gabonenses (24) A New Species in Utsetela (Moraceae)[END_REF]. Only one species of Utsetela was sampled in this study. Considering the comparatively long branches of U.

gabonensis, D. alta, and D. mannii, this relationship could be caused by a long-branch attraction artefact. To test this, we excluded the three species of Dorstenia which taxonomic sampling and use of many more genes, and our reconstructions from both ML and coalescent approaches strongly supported the non-monophyletic status of these sections. An alternative taxonomic solution to sinking Scyphosyce and Utsetela into Dorstenia is to include some Dorstenia species (most of section Nothodorstenia and at least one member of section Xylodorstenia -D. angusticornis) into the genus Scyphosyce. Dorstenia africana, D. kameruniana, D. oligogyna, D. djettii and D.

dorstenioides (the former four all included in section Nothodorstenia) were once classified as genus Craterogyne [START_REF] Lanjouw | Studies in Moraceae I. The genera Trymatococus Poepp. et Endl. and Graterogyne Lanj[END_REF]. Dorstenia dorstenioides, which has been proposed as the link between sections Xylodorstenia and Nothodorstenia (Berg and Hijman 1999), was excluded in the present study due to the low number of genes represented in our main analyses (Table S1). As Scyphyosyce [START_REF] Baillon | Stirpes exoticae novae[END_REF])is an older name than Craterogyne, Scyphosyce would take priority for the name of a new genus.

Regarding the clade containing Utsetela, some of the Dorstenia species in that clade were recently transferred to a new genus (Maria) established by [START_REF] Vianna | Maria,a new genus of Moraceae[END_REF], comprising four species in Dorstenia section Xylodorstenia (D. alta, D.

angusticornis, D. scaphigera, D. turbinata). Maria was later found to be a homonym and renamed Hijmania [START_REF] Vianna Filho | Hijmania, a replacement name for Maria (Moraceae)[END_REF]. It may be necessary to include those species into the same genus under the name Utsetela [START_REF] Pellegrin | Utsetela Pellegrin, genre nouveau d'Urticacées-Artocarpées du Gabon[END_REF]), which has priority, but until more complete taxon sampling is completed, we do not presently propose any taxonomic changes. One of the woody African macrospermous species that warrants further attention is D. elliptica. It was included in section Nothodorstenia by Berg (1978) because it had bracts resembling other members of that section. However both Misiewicz & Zerega (2012) and the present study, found D. elliptica to be sister to all Neotropical Dorstenia.

Among herbaceous species of Dorstenia, the presence of bracts on receptacles has traditionally been used to distinguish among sections Emygdioa, Dorstenia, Lecanium on the one hand and sections Acauloma, Bazzemia, Lomatophora, Kosaria on the other hand (Berg and Hijman 1999). Although none of these sections were found to be monophyletic, the Neotropical clade contains all of the species with bracteate receptacles except D. picta (Figure S1). Our results suggest that traditional morphological characters for sectional delimitation within Dorstenia do not hold up to molecular phylogenetic scrutiny, that a close examination of alternative characters is needed, and a new intrageneric classification is warranted.

Conclusion

The targeted enrichment sequencing strategy, paired with the HybPiper pipeline, proved to be an effective approach at reconstructing phylogenetic relationships in Dorstenieae using herbarium specimens. Further molecular and morphological work will be required before solving some of the taxonomic issues highlighted in this study, such as sinking Utsetela and Scyphosyce into Dorstenia or separating some Dorstenia species into either the genera of Utsetela or Scyphosyce. Dorstenieae as a whole may have originated in the joint area of continental Africa, Madagascar and Asia-Oceania area, followed by at least two independent colonizations of South America (Brosimum s.l. and Dorstenia s.l.). Some species in these two clades further dispersed to Central/North America. The mechanical processes for the long-distance dispersal of species of Dorstenieae remains an enigma. More studies on pollination and dispersal in this tribe will be required to further elucidate the biogeographic history of this group.

The development of new biogeographic models and new model selection procedures will also be essential to help to clear the biogeographic history of Dorstenieae and other pantropically distributed lineages. The robust and most densely sampled Dorstenieae phylogeny presented here will be a valuable resource for further studies on character evolution in this fascinating tribe and will assist with future taxonomic revisionary work.

INTRODUCTION

The genus Ficus (ca. 750 spp., Moraceae) has long intrigued researchers given its taxonomic and ecological diversity and the coevolution with its wasp pollinators (Agaonidae, Cruaud et al. 2012). Two different breeding systems have been observed in Ficus so far: monoecy and gynodioecy, the latter is functionally similar to dioecy in this genus (hereafter referred to as dioecy, Cook and Rasplus 2003). Because of the physiological and ecological differences between breeding systems, Ficus species of the two breeding systems differ not only in their inflorescences (syconia), but are also expected to differ in dispersal ability, distribution and other traits (Harrison and Yamamura 2003;[START_REF] Nazareno | Dioecy, more than monoecy, affects plant spatial genetic structure: the case study of Ficus[END_REF][START_REF] Yang | The incidence and pattern of copollinator diversification in dioecious and monoecious figs[END_REF].

Ficus is a typical tropical clade, but its global distribution is uneven across tropical areas (Fig. 1). Approximately 120, 105 and 510 species are known from the Neotropics, Africa (including Madagascar and adjacent Indian Ocean islands, and the Arabian Peninsula) and Asia-Oceania, respectively (Berg et al. 2005). Broad-scale distribution patterns also differ between breeding systems (Fig. 1). Dioecious Ficus does not occur in the Neotropics (Cruaud et al. 2012) and is mainly distributed in Asia-Oceania. Several dioecious species of Ficus (e.g., F. sarmentosa, F. heteromorpha, F. erecta, F. coronata, F. tikoua and F. pumila) have extended their distribution into more temperate areas in Asia (Berg et al. 2005;[START_REF] Zhao | Phenological adaptations in Ficus tikoua exhibit convergence with unrelated extra-tropical fig trees[END_REF]. Monoecious Ficus occurs across all tropical realms and also occurs in dryer areas, for example in dry parts of Atlantic Costal regions or (semi-)desert of northern Mexico (Berg 2001). However, the underlying physiological and ecological causes for these distinct distribution patterns remain debated.

Dispersal distances of Ficus species and their pollen are strongly influenced by dispersal abilities of their wasp pollinators. Pollinators of monoecious Ficus were observed in the canopy layer in both Bornean and Neotropical forests, whereas the pollinators of dioecious Ficus occur most often in lower canopy layers [START_REF] Hespenheide | Selective predation by two swifts and a swallow in Central America[END_REF][START_REF] Harrison | Fig wasp dispersal and the stability of a keystone plant resource in Borneo[END_REF]. Considering the higher wind speed in the canopy, pollinating wasps of monoecious hosts would have longer dispersal distances than the counterpart of dioecious Ficus, suggesting that the pollen of monoecious Ficus, on average, disperse farther [START_REF] Compton | The flight heights of chalcid wasps (Hymenoptera, Chalcidoidea) in a lowland Bornean rain forest: fig wasps are the high fliers[END_REF]. The pollinators of Ficus sycomorus (monoecious) were found to disperse tracking the seasonal nighttime wind and the average dispersal distance was 88.6 km in Africa (160 km maximum, [START_REF] Ahmed | Wind-borne insects mediate directional pollen transfer between desert fig trees 160 kilometers apart[END_REF]). As a result of more frequent long-distance dispersal of pollinators (pollen), monoecious Ficus populations showed lower levels of genetic structure than dioecious Ficus populations [START_REF] Nazareno | Dioecy, more than monoecy, affects plant spatial genetic structure: the case study of Ficus[END_REF].

Besides pollen dispersal, the figs (syconia) are mainly eaten and dispersed by a wide range of birds, arboreal mammals and fruit bats, and while considering the diverse habits of Ficus, it is difficult to address differences in seed dispersal distances between breeding systems (Shanahan et al. 2001).

In a previous study, ancestral breeding systems were reconstructed in Ficus and Moraceae (Zhang et al. 2019). Although some uncertainty remained, depending on models and approaches, monoecy was shown to be the likely ancestral state in Ficus, and later evolved into dioecy once or twice, with at least three subsequent reversals to monoecy. However, the conditions in which these transitions occurred remain unclear. Understanding the ecological difference of the two breeding systems in Ficus will help us better understand the environmental conditions of breeding system transitions not only in Ficus but also in angiosperms. Phylogenetic comparative methods provide the tools to investigate the relationship between breeding systems and environmental niche (Table 1). Phylogenetic independent contrasts [START_REF] Felsenstein | Phylogenies and the Comparative Method[END_REF] and directional comparative methods [START_REF] Lauder | Form and function: structural analysis in evolutionary morphology[END_REF] were first used in comparative biology to consider the phylogenetic relationship of the data when testing the correlations among traits. These two methods were later found to be very similar, as both of them estimate the evolutionary regression coefficient [START_REF] Pagel | Seeking the evolutionary regression coefficient: an analysis of what comparative methods measure[END_REF]. The directional comparative approach tests the correlation by first estimating the ancestral state for each internal node. Considering the uncertainty of topology and ancestral state reconstructions, Pagel (1994) proposed a maximum likelihood approach which estimates the correlation of continuous traits without reconstructing ancestral states. Another approach based on general linear models called phylogenetic generalized least squares (PGLS) was developed for investigating the relationships among continuous traits while taking the entire phylogeny into account [START_REF] Grafen | The phylogenetic regression[END_REF][START_REF] Martins | Phylogenies and the comparative method: a general approach to incorporating phylogenetic information into the analysis of interspecific data[END_REF][START_REF] Pagel | Inferring evolutionary processes from phylogenies[END_REF][START_REF] Mundry | Statistical issues and assumptions of phylogenetic generalized least squares In: Modern phylogenetic comparative methods and their application in evolutionary biology[END_REF]). However, constrained by the algorithm, the PGLS approach can only deal with continuous variables. This deficiency was solved by the development of the generalized estimating equations (GEE) approach, which can analyze both continuous and discrete data and more complex models [START_REF] Paradis | Analysis of comparative data using generalized estimating equations[END_REF].

Our study aims to understand how niche differences (temperature and precipitation) between breeding systems may have influenced the global distribution of Ficus. Using novel approaches to study correlated evolution, we are now able to address four key questions: 1) are dioecious species of Ficus distributed across lower temperatures (of the coldest month) than monoecious species of Ficus? 2) do monoecious species of Ficus occur more often in drier environments than species of dioecious Ficus? 3) is the evolution of dioecy in Ficus associated with the colonization of cooler niches? 4) is the evolution of monoecy in Ficus associated with the colonization of dry climatic niches)? To investigate these four key questions, we assembled and cleaned a comprehensive database of Ficus occurrence records, extracted climatic niche variables for the species, and conducted phylogenetic comparative analyses with the GEE approach for climatic variables against the breeding systems in Ficus. This study sheds light on the understanding of distribution patterns and niche differences in Ficus according to breeding systems, and helps to understand the relationship of breeding systems and distribution in Moraceae and angiosperms more generally in a climate changing world.

MATERIALS AND METHODS

Occurrence records download and cleaning

We chose 200 species of Ficus (Cruaud et al. 2012;Zhang et al. 2019) to make the following tests. A total of 102,486 occurrence records were downloaded from the Global Biodiversity Information Facility (GBIF) with R package rgbif [START_REF] Chamberlain | R Python, and Ruby clients for GBIF species occurrence data[END_REF][START_REF] Chamberlain | rgbif: Interface to the Global Biodiversity Information Facility API[END_REF]) by searching by species name. We then flagged and excluded potentially erroneous coordinates (e.g. both latitude and longitude are zero; same latitude and longitude; coordinate on the sea; coordinate as the center or capital city of one country and so on) using R package WorldClim), we did not include them in the following analyses. The fourteen included variables and their names used in the analyses are listed in Table 1. The median of each variable of each species was calculated to conduct the principal component analysis (PCA). Considering the value of precipitation data are one order of magnitude larger than those of temperature data, we normalized climatic variables related with precipitation (bio12-19) by taking their square root before PCA. Principal components 1 and 2 (PC1 and PC2) were calculated with R package ade4 [START_REF] Chessel | The ade4 package-I-One-table methods[END_REF]Dray and Dufour 2007;Dray et al. 2007;[START_REF] Bougeard | Supervised multiblock analysis in R with the ade4 package[END_REF].

We then conducted two types of regression, differing fundamentally in whether or not phylogeny was taken into account. First, we used generalized linear models (GLM) to understand how the species' climatic niches influence breeding system variation, without taking phylogeny into account. To this end, breeding system was regressed against PC1, PC2, PC1 and 2, and for each individual climatic variable using logistic regressions by function compar.gee in R package ape [START_REF] Paradis | ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R[END_REF]. However, the correlation between climatic niches and breeding systems may simply result from phylogenetic dependence (i.e. inherited from the ancestral species [START_REF] Felsenstein | Phylogenies and the Comparative Method[END_REF]. Thus, we used the generalized estimating equations (GEE) approach with logistic distribution (breeding system is a discrete tratit, [START_REF] Paradis | Analysis of comparative data using generalized estimating equations[END_REF] to test the correlation of PC1, PC2, PC 1 and 2, and each climatic variables versus breeding systems, while correcting for the phylogenetic effect. For these analyses, we used the ultrametric (dated)

Maximum Clade Credibility tree from Zhang et al. (2019) but only kept the tips of the 183 Ficus species in the geographic dataset. The results from the two regression methods were then compared.

climate variables of the geographic coordinates where each species has the lowest precipitation of the warmest quarter (bio18), rather than the media value of the entire set of occurrence records of each species used in the main dataset. The rationale for this analysis stems from the fact that the median value of climatic variables may not adequately represent the climatic tolerance of the species to extreme values. Hereafter, we refer to this dataset as bio18min. For both the subsetAsia and bio18min datasets, we conducted principal component analysis, regression (GLM, GEE) and ancestral state reconstruction in the same way as described above.

RESULTS

Global distribution pattern of Ficus according to breeding systems

There were 60,432 occurrence records representing 183 species (105 monoecious, 78 gynodioecious, Table S1) in the final dataset. Most of the distribution area of monoecious Ficus overlapped with that of dioecious Ficus in Africa (incl. Madagascar) and Asia-Oceania (Fig. 1).

Monoecious species of Ficus were also observed in comparatively drier areas of Africa and Australia (Fig. 1). In Asia, dioecious Ficus extends further North than monoecious Ficus, reaching the temperate zone, including Japan and central China (Fig. 1). Continental Southeast Asia (including southern China), northern Borneo and New Guinea hold the highest concentration of dioecious species of Ficus (Fig. 1A). Central America, the tropical Andes, West and Central Africa, eastern South Africa, and the Australian Wet Tropics are the hotspots of monoecious Ficus (Fig. 1B).

Principal component analysis and regressions

In the main analysis (median of each climatic variable), principal component 1 and 2 (PC1 and PC2) explained 45.6% and 34.4%, respectively, of the variation in climatic niche variables across all Ficus (Fig. 2). PC1, PC2, their combination were significantly negatively correlated with dioecy. In addition, all of the single variables related with precipitation (bio 12-19) were significantly positively correlated with dioecy in regressions using generalized linear models (GLM) analyses (Table 1), and thus negatively associated with monoecy. Conversely, most variables related to temperature were not significantly associated with breeding systems, except the maximum temperature of warmest month (bio5), which was negatively correlated with dioecy.

This means that monoecious species occur in drier (and sometimes warmer) climates than dioecious Ficus. In the regression using the generalized estimating equations (GEE) approach, which considers the phylogenetic non-independence of the data, none of the principle components and their combination, nor any single climatic variable was significantly correlated with breeding system variation. Results were almost identical in the two additional analyses (subsetAsia and bio18min), with a few exceptions (Table S2). In GLM analyses of both the subsetAsia and bio18min dataset, the same direction of significance was found for all variables except bio19 with the bio18min dataset. However, the strength of significance was weaker for all variables of the two additional datasets than in the analyses of the main dataset. No significance was found for any variable in GEE analyses of either additional dataset, consistent with the GEE analysis of the main dataset.

Temperature may not relate with the distribution of dioecious Ficus

Although temperature was suggested to be an important factor for the development of Ficus at its northern boundary in previous studies of several species [START_REF] Peng | The reproductive success of Ficus altissima and its pollinator in a strongly seasonal environment: Xishuangbanna, Southwestern China[END_REF][START_REF] Chen | Latitudinal effects on phenology near the northern limit of figs in China[END_REF], our results do not support any relationship between temperature and breeding systems in Ficus as a whole. Maximum temperature of warmest month (bio5) is the only temperature variable found to have a (negative) significant relationship with dioecy in GLM analyses, a result not supported by the GEE approach. In the GLM analysis with the subsetAsia dataset, bio5 was again the only variable showing (negatively) a significant relationship with dioecy.

Monoecious Ficus occur in drier environments more often than dioecious Ficus

Our richness maps of Ficus with two breeding systems imply that precipitation has shaped the distribution pattern within the climatic boundaries of the genus. The overlapping areas of the two breeding systems are areas with higher humidity (Fig. 1). On the African, Asian and Australian continents, this distribution pattern coincides with the summer position in both the northern and southern hemispheres of the intertropical convergence zone (ITCZ), which is identified as the maximum in time-mean precipitation [START_REF] Schneider | Migrations and dynamics of the intertropical convergence zone[END_REF]) (Fig. 1). Ficus is mainly distributed in peninsulas, islands, archipelagos and coastal areas in Asia-Oceania, where climates are humid.

This region is also one of the distribution centers for both breeding systems. In addition, most of monoecious Ficus in the Neotropics is distributed in the ITCZ area (Fig. 1B).

In addition to the overlapping areas of both breeding systems, monoecious species also occupy areas with drier climates in continental Africa, Madagascar and Australia. In continental Africa, monoecious species occur in Namibia, northern South Africa, Zimbabwe and Ethiopia, which have been suggested to be part of an arid "corridor", spreading from southwest to northeast of continental Africa [START_REF] Jürgens | Floristic biodiversity and history of African arid regions[END_REF] and documented in other groups of flowering plants [START_REF] Thiv | Old-New World and trans-African disjunctions of Thamnosma (Rutaceae): Intercontinental long-distance dispersal and local differentiation in the succulent biome[END_REF][START_REF] Freitas | Phylogenetic and morphological investigation of the Mochlus afer-sundevallii species complex (Squamata: Scincidae) across the arid corridor of sub-Saharan Africa[END_REF]. In Madagascar, most of dioecious Ficus occurs in the humid forest along the east coast [START_REF] Evans | Insights on the evolution of plant succulence from a remarkable radiation in Madagascar (Euphorbia)[END_REF], while monoecious Ficus is distributed throughout the country (Fig. 1B) including the dry west edge. In Australia, with some exceptions, dioecious Ficus is distributed in the Monsoonal Tropics and aseasonal-wet areas [START_REF] Greenwood | Eocene monsoon forests in central Australia?[END_REF][START_REF] Crisp | Radiation of the Australian flora: what can comparisons of molecular phylogenies across multiple taxa tell us about the evolution of diversity in present-day communities?[END_REF]).

However, monoecious species are also distributed in the arid Eremean zone [START_REF] Crisp | Radiation of the Australian flora: what can comparisons of molecular phylogenies across multiple taxa tell us about the evolution of diversity in present-day communities?[END_REF]) of central and western Australia. The distribution pattern of Ficus in Asia is different. Dioecious species of Ficus are distributed farther North than monoecious species (Ficus sarmentosa extends to 37°N in our dataset, Fig. 1). The northern boundary of dioecious Ficus in continental Asia coincides with the summer (July) position of ITCZ [START_REF] Schneider | Migrations and dynamics of the intertropical convergence zone[END_REF], while monoecious species seem to be constrained in tropical areas with few records in subtropical areas (Fig. 1B). In addition to its presence in typically year-round arid regions, monoecious Ficus also occupies areas with temporary dry periods, where dioecious Ficus is absent (Fig. 1). For instance, monoecious

Ficus occurs in the Cape Floristic Region of South Africa, the most species-rich Mediterraneantype ecosystem given its size (Linder 2014) We hereby propose several explanations for the colonization of drier environments by monoecious Ficus. Firstly, the dispersal ability of dioecious Ficus is limited compared to that of monoecious species. Dioecious species of Ficus contribute more than monoecious ones to spatial genetic structure [START_REF] Nazareno | Dioecy, more than monoecy, affects plant spatial genetic structure: the case study of Ficus[END_REF]. Mean parent-offspring distance has been estimated to be 200 m in a population genetic structure investigation of two dioecious Ficus [START_REF] Dev | Fine-scale population genetic structure of two dioecious Indian keystone species, Ficus hispida and Ficus exasperata (Moraceae)[END_REF], while pollen of a monoecious species (F. sycomorus) was reported to be delivered as far as 160 km [START_REF] Ahmed | Wind-borne insects mediate directional pollen transfer between desert fig trees 160 kilometers apart[END_REF]. Secondly, disadvantage of the separation of sex in dioecious Ficus is another possible explanation. With both functional male and female flowers on the same individual, monoecious Ficus may self-fertilize (wasps from one syconium may pollinate another syconium of the same plant individual) in a highly seasonal environment [START_REF] Cook | Effects of within-tree flowering asynchrony on the dynamics of seed and wasp production in an Australian fig species[END_REF], while dioecious Ficus cannot. Third, seasonal climate, for instance dry summers, may impact the production of seeds in dioecious Ficus. Summer is an important period for the development of both host plants and larva of pollinators [START_REF] Peng | The reproductive success of Ficus altissima and its pollinator in a strongly seasonal environment: Xishuangbanna, Southwestern China[END_REF]. At least one crop of wasps released from male syconia of dioecious Ficus will serve mainly the purpose of pollination of female syconia in a period (year), rather than laying eggs in other male syconia. In dioecious Ficus, this crop tends to be in the summer [START_REF] Zhao | Phenological adaptations in Ficus tikoua exhibit convergence with unrelated extra-tropical fig trees[END_REF]; but see [START_REF] Chen | Winter cropping in Ficus tinctoria: an alternative strategy[END_REF], when higher evapotranspiration implies additional constraints on water use for seed development. Therefore, dry summers may lead to reduction of seeds in the next phase of the hosts' life cycle. These three hypotheses may explain the absence of dioecious Ficus in Mediterranean climates, and dry climates more generally.

Conclusion

Ficus is a pan-tropically distributed genus occupying several bioregions [START_REF] Crisp | Radiation of the Australian flora: what can comparisons of molecular phylogenies across multiple taxa tell us about the evolution of diversity in present-day communities?[END_REF]Hughes et al. 2012;Linder 2014). Different climatic variables may played main roles in different bioregions [START_REF] González-Orozco | A biogeographical regionalization of Australian Acacia species[END_REF]. Thus, analyses combining species from different bioregions may suffer from confounding factors. We sampled 183 species (ca. 24% of all Ficus species) in this study. Future more densely sampled phylogenetic trees will be required to start disentangling the main environmental factors in different bioregions. However, phylogenetic relationships in Ficus remain themselves incompletely understood (Cruaud et al., 2012;Bruun-Lund et al. 2017), potentially limiting our understanding of breeding system evolution in the genus (Zhang et al. 2019). A new, highly supported phylogenetic backbone of Ficus was proposed in a recent study based on 600 RAD-seq loci, in which subgenus Sycomorus was found to be sister to the remaining of Ficus (Rasplus et al. 2018), contrary to previous work (Cruaud et al. 2012;Pederneiras et al. 2018) . If confirmed, this significant change of topology of phylogeny may affect the results of comparative analyses such as those presented in this study.

The overlapping distribution areas of monoecious and dioecious Ficus correspond to the July and January positions of the intertropical convergence zone (ITCZ) on continents, consistent with the idea that precipitation may have played a role in determining the general distribution pattern of Ficus within the climatic boundaries of the genus. A plausible positive relationship between precipitation and dioecy was found in Ficus, and conversely between tolerance to drier environments and monoecy in Ficus. The longer dispersal ability and potential for self-fertilization in monoecious Ficus are two hypotheses we propose to explain this pattern. Contrary to previous ideas, no significant relationship of temperature and breeding system was detected. A more densely sampled phylogenetic tree of Ficus and further improvements of phylogenetic comparative methods will be necessary for confirming the results obtained in this study. I investigated breeding system evolution, biogeographic history and their impacts on the diversity of angiosperms in this thesis, using the pantropically distributed family Moraceae as a model clade. My work has led to several important new results: 1) dioecy was probably the ancestral state of breeding system in Moraceae and subsequently evolved into other three breeding systems, suggesting that dioecy is not an evolutionary dead end; 2) the origin of tribe Dorstenieae was estimated to be in the Cretaceous (65.8-79.8 Ma); two independent dispersal events from Africa to South America occurred in Dorstenieae, leading to Brosimum s.l. and the Neotropical clade of Dorstenia; 3) precipitation appears to have been a factor in shaping the biogeographic distribution pattern of Ficus according to breeding systems.

In Chapter I, I reconstructed the ancestral state of breeding systems in Moraceae and Ficus with a new dated tree built for the purpose and six morphological models differing in the number of breeding system categories. While this study highlighted that dioecy was the starting point, rather than a dead end of breeding system evolution in the family, much uncertainty remained with respect to the ancestral state of Ficus. This uncertainty is partly the consequence of a combination of both monoecious and dioecious lineages close to the base of the Ficus phylogeny. However, basal relationships in Ficus themselves remain uncertain and it is possible that future phylogenetic work in the genus will lead to relationships that more clearly support either monoecy or gynodioecy as the ancestral state of Ficus. A new phylogenomic tree of 40 species of Ficus based on RAD-seq data has recently been proposed (Rasplus et al. 2018), in which, all the phylogenetic relationships of subgenera or sections were highly supported and subgenus Sycomorus was sister to the remaining of Ficus. This topology is different from ours and two other previous phylogenies of Ficus reconstructed from nuclear markers (Cruaud et al. 2012;Pederneiras et al. 2018), in which these results imply two long-distance dispersal (LDD) events from continental Africa to South America lead to explain the origin of Brosimum s.l. and Neotropical Dorstenia. The mechanism for the LDD of Dorstenia and Brosimum remain unclear, especially given that Dorstenia has been suggested to be poorly adapted for LDD (Berg and Hijman 1999). Further studies of dispersal related traits of Dorstenia and Brosimum will improve the understanding of the mechanism of LDD for species unlikely to do so [START_REF] Onstein | To adapt or go extinct? The fate of megafaunal palm fruits under past global change[END_REF], a general and difficult question that has been raised in numerous other plant clades [START_REF] Doyle | Phylogeny, molecular and fossil dating, and biogeographic history of Annonaceae and Myristicaceae (Magnoliales)[END_REF][START_REF] Barker | Molecular dating of the 'Gondwanan'plant family Proteaceae is only partially congruent with the timing of the break-up of Gondwana[END_REF]).

The latest biogeographic study of Moraceae as a whole was provided by Zerega et al. (2005), processes may be biased towards cladogenetic processes, and as a result over-reliance on the founder-event speciation process to explain biogeographic distributions (Ree and Sanmartín 2018).

Most of the genera in Moraceae consist of fewer than 10 species (Rohwer and Berg 1993), and some species may be hard to collect in the field. As demonstrated in Chapter II and several recent studies (Hart et al. 2016;Couvreur et al. 2019) Species have preferentially remained in the same or similar niche/biome while colonizing new environments [START_REF] Crisp | Phylogenetic biome conservatism on a global scale[END_REF]), which may have played an important role in shaping the current distribution pattern of plant diversity. In Chapter III, I combined knowledge of breeding system evolution and detailed spatial data to characterize and understand differences in the geographic distribution of monoecious and gynodioecious species of Ficus. At the beginning of this PhD and chapter, I hypothesized that temperature was a key factor influencing the distribution of Ficus, based on the tendency of gynodioecious Ficus to occupy higher latitudes and colder environments in Asia than monoecious Ficus (Berg et al. 2005;[START_REF] Bain | Plasticity and diversity of the phenology of dioecious Ficus species in Taiwan[END_REF]. However, none of my analyses supported a correlation of temperature with breeding system in Ficus, even with a subset of the data restricted to the Asia-Oceania region. Precipitation was found to be a plausible factor in shaping the current distribution pattern of different breeding systems in Ficus: both breeding systems of Ficus are largely overlapping in tropical climates, but monoecious Ficus extends to drier environments than gynodioecious species. To our knowledge, this is the first time that such a correlation has been suggested and discussed.

In Chapter I, I estimated the age of crown group Ficus in the Eocene (40.6-55.9 Ma), during when the Eurasia region, which has been proposed as the ancestral area of Ficus (Pederneiras et al. 2018), was still separated by the Siberian and Turan Seas as Europe and Asia [START_REF] Akhmetiev | Paleogene floral assemblages around epicontinental seas and straits in Northern Central Eurasia: proxies for climatic and paleogeographic evolution[END_REF]. In addition, the average global temperature was higher and tropical regions were much more widespread than they are today [START_REF] Zachos | Trends, rhythms, and aberrations in global climate 65 Ma to present[END_REF][START_REF] Morley | Cretaceous and Tertiary climate change and the past distribution of megathermal rainforests In: Tropical rainforest responses to climatic change[END_REF]). The distribution areas of Ficus may have extended and retreated several times to track tropical climates during the history of Ficus [START_REF] Costa | The phylogeography of two disjunct Neotropical Ficus (Moraceae) species reveals contrasted histories between the Amazon and the Atlantic Forests[END_REF]. During these processes, some monoecious Ficus dispersed into the current drier areas, which may have been previously wetter. However, the cooccurring gynodioecious species would have been swept away during past aridification events due to their lower tolerance of dry climate. In addition, monoecious Ficus probably also successfully dispersed into drier areas where the genus was previously absent. To our knowledge, ours is the first study to test climatic niche difference of congeneric species with different breeding systems.

Whether niche difference is a general pattern in other diversified genera with at least two breeding systems such as Acer (Sapindaceae, [START_REF] Renner | The evolution of dioecy, heterodichogamy, and labile sex expression in Acer[END_REF] remains to be tested in more clades.

Furthermore, the loci responsible for sex determination in Moraceae are still unclear. Did any sex chromosomes develop in the Moraceae genome? Are they the same in other model clades such as Carica (Caricaceae) [START_REF] Wang | Sequencing papaya X and Y h chromosomes reveals molecular basis of incipient sex chromosome evolution[END_REF], Cucurbitaceae [START_REF] Boualem | A cucurbit androecy gene reveals how unisexual flowers develop and dioecy emerges[END_REF] and Silene Codes for distribution areas were the same as in Figure 1 (A, South America; B, Central/North America; C, continental Africa; D, Madagascar; E. India and Sri Lanka; F, Southeast Asia and Oceania.)

Table S2. Estimates for interaction effects of climatic variables and breeding systems in Ficus from the regressions by generalized linear models (GLM) and generalized estimating equations (GEE) for three datasets: 1) full dataset (same as Table 1); 2) subsetAsia; and 3) bio18min. In this thesis, I investigated breeding system evolution and biogeographic history in the family Moraceae, which I used as a model clade to understand the origin and evolution of diversity of angiosperms. In Chapter I, I reconstructed and calibrated a new dated phylogenetic tree for Moraceae as a whole. I then used this tree to reconstruct ancestral states of breeding systems in Moraceae and Ficus. The crown group ages of Moraceae and Ficus were estimated in the Cretaceous and in the Eocene, respectively. Dioecy was inferred as the ancestral breeding systems of Moraceae, with several subsequent transitions to monoecy, including in Ficus. This result suggests that dioecy is not necessarily an evolutionary dead end. In Chapter II, I reconstructed a dated phylogenetic tree for tribe Dorstenieae, mainly distributed in tropical regions, with a new data set of nuclear genomic data generated with a Hyb-Seq approach. Biogeographic history was then reconstructed using dispersalextinction-cladogenesis models. The crown group ages of Dorstenieae and Dorstenia were estimated in the Cretaceous and in the Cretaceous/Paleocene period, respectively. Two long-distance dispersal events from continental Africa to South America occurred in the Cenozoic (Dorstenia and Brosimum s.l.). In Chapter III, I tested the climatic niche difference (temperature and precipitation) between the two breeding systems (monoecy and gynodioecy) in Ficus using a new dataset of cleaned spatial occurrence records and breeding systems for 183 species. I used two comparative approaches: generalized estimating equations (GEE) and generalized linear models (GLM). A positive relationship between precipitation and gynodioecy was supported by GLM, but not GEE analyses, and no relationship between temperature and breeding systems was supported by either method. Higher dispersal ability and the potential for self-fertilization may explain why monoecious species of Ficus have been able to colonize and survive in drier environments. This thesis highlights the potential of phylogenetic comparative methods and phylogenomic data to address questions of breeding system evolution and biogeography in Moraceae, and opens up several important new perspectives worth investigating in other plant clades, such as a relationship between breeding system and climatic niche.

  le clade le plus diversifié des plantes actuelles et sont exceptionnellement riches en espèces dans les régions tropicales. La fleur est l'organe reproducteur des Angiospermes. Il existe une diversité remarquable de fleurs tant sur le plan structurel que fonctionnel et pour cette raison, l'on pense depuis longtemps que les fleurs ont une influence directe sur la diversité et l'évolution des Angiospermes. Dans cette thèse, j'ai étudié l'évolution des systèmes sexuels et l'histoire biogéographique de la famille des Moraceae (40 genres, env. 1100 espèces, dont le genre le plus diversifié, Ficus, env. 750 espèces), clade modèle utilisé pour comprendre l'origine et l'évolution de la diversité chez les Angiospermes. Dans le Chapitre I, j'ai reconstruit et calibré un nouvel arbre phylogénétique daté pour les Moraceae à partir d'un jeu de données de 320 espèces (dont 272 issues de 36 genres chez les Moraceae, et 48 espèces représentant 8 familles de groupes externes chez les Rosales), 8 marqueurs moléculaires (3 chloroplastiques, 5 nucléaires) et 12 contraintes d'âge fondées sur le registre fossile. J'ai ensuite utilisé cet arbre phylogénétique pour reconstruire les états ancestraux des systèmes sexuels chez les Moraceae et Ficus avec les méthodes de parcimonie et de maximum de vraisemblance, en utilisant et comparant six configurations de données de systèmes sexuels chez les Moraceae différant par le nombre d'états distingués. Les âges des groupes-couronne des Moraceae et du genre Ficus sont estimés au Crétacé (73.2-84.7 Ma) et à l'Eocène (40.6-55.9 Ma), respectivement. La dioécie est inférée comme l'état ancestral des systèmes sexuels chez les Moraceae, un résultat très bien soutenu par toutes les approches, modèles et configurations. Plusieurs transitions ultérieures vers la monoécie se sont ensuite produites, y compris chez Ficus. Ce VII résultat suggère que la dioécie ne représente pas nécessairement un cul-de-sac évolutif. La monoécie semble être l'état ancestral le plus vraisemblable chez Ficus, mais ce résultat est peu robuste et sensible au choix du modèle. Dans le Chapitre II, j'ai reconstruit un arbre phylogénétique daté pour la tribu des Dorstenieae (15 genres, env. 156 espèces, dont Dorstenia, env. 113 espèces), distribuée principalement dans les régions tropicales. Pour ce faire, j'ai produit un nouveau jeu de données génomiques nucléaires (102 gènes) à partir d'une approche d'enrichissement ciblé (Hyb-Seq) échantillonnant tous les genres (15) et 83 espèces (53%), ainsi que 5 groupes externes. Cet arbre phylogénomique m'a permis de reconstruire l'histoire biogéographique du groupe en utilisant les modèles de dispersion-extinction-cladogenèse, dont des modèles avec un processus d'événement-fondateur de dispersion. Les âges des groupes-couronne des Dorstenieae et du genre Dorstenia sont estimés au Crétacé (65.8-79.8 Ma) et dans la période du Crétacé/Paléocène à l'Eocène (50.8-67.3 Ma), respectivement. Deux évènements de dispersion à longue distance depuis l'Afrique continentale vers l'Amérique du Sud ont eu lieu au Cénozoïque, l'un chez Dorstenia et l'autre chez Brosimum s.l. Le genre Dorstenia a ensuite colonisé l'Amérique Centrale entre l'Oligocène et le Miocène (12.0-34.7 Ma). Dans le Chapitre III, j'ai testé les différences de niche climatique (température et précipitation) entre les deux systèmes sexuels (monoécie et gynodioécie) chez Ficus avec un nouveau jeu de données fiables d'occurrences spatiales et de systèmes sexuels chez 183 espèces. À cette fin, j'ai utilisé deux approches comparatives, différant dans la prise en compte de la phylogénie (équations d'estimation généralisées, GEE) ou non (modèles linéaires généralisés, GLM). Une relation positive entre précipitation et gynodioécie est VIII soutenue par les analyses GLM, et aucune méthode ne soutient une relation entre température et système sexuel. Une meilleure capacité à se disperser et le potentiel d'autopollinisation sont deux explications possibles pour la colonisation et la survie des espèces monoïques dans des environnements plus secs. Cette thèse démontre le potentiel des méthodes phylogénétiques comparatives et des données phylogénomiques pour répondre aux questions d'évolution des systèmes sexuels et de biogéographie chez les Moraceae. Des phylogénies plus densément échantillonnées pour Ficus et les Moraceae sont requises pour confirmer certains des résultats émergeant de cette thèse, tels que le caractère ancestral de la monoécie chez Ficus. Cependant, cette thèse ouvre plusieurs nouvelles perspectives importantes méritant d'être approfondies chez d'autres clades de plantes, telles que la relation entre système sexuel et niche climatique.

Figure 1 .

 1 Figure 1. Geographic distribution and species richness of Moraceae. This map was compiled using a dataset of cleaned occurrence records from the Global Biodiversity Information Facility (GBIF), following the same protocol as described in Chapter III. To decrease the impact from cultivated individuals, records of widely cultivated species in Moraceae (Artocarpus heterophyllus, A. altilis, Broussonetia papyrifera, Castilla elastica, Ficus benjamina, F. benghalensis, F. carica, F. elastica, F. lyrata, F. microcarpa, F. pumila, F. religiosa, Morus alba, M. nigra) were excluded from this map.

Figure 2 .

 2 Figure 2. Representative species of the seven tribes of Moraceae. A-B. Morus alba female and male inflorescence (Moreae); C. Artocarpus heterophyllus (Artocarpeae); D-E. Parartocarpus venenosus female and male inflorescence (Parartocarpeae); F. Castilla elastica (Castilleae); G. Maclura pomiera syncarp (Moreae); H. Dorstenia barteri (Dorstenieae); I. Ficus carica female inflorescence (Ficeae). Photo: A,I by Qian Zhang; H by Hervé Sauquet; C from Williams et al. (2017); D, E from Zerega and Gardner (2019); G from (Gardner et al. 2017), reproduced with authors' permission; B, F from Wikimedia commons by Fastily and Dick Culbert, respectively.

Figure 3 .

 3 Figure 3. Diagram of syconia (enclosed, urn-shaped inflorescences) in the two breeding systems of Ficus, monoecy and gynodioecy.

Figure 1 .

 1 Figure 1. Schematic representation of the six configurations for the breeding system character analyzed in this study. The arrows denote transition rate parameters, as estimated in the reversible-jump Markov Chain Monte Carlo analyses, with mean rates reported. The absence of an arrow indicates a near-zero estimate for the corresponding parameter, suggesting the model does not support direct transition between the two states.

Figure 2 .

 2 Figure 2. Dated phylogenetic tree of 320 species with family names in Rosales, tribe names in Moraceae, and section names in Ficus. This is the Maximum Clade Credibility (MCC) from the BEAST analysis of eight molecular markers (see text). Fossil-calibrated nodes are indicated in red. The posterior distribution of estimated ages is shown for crown-group Moraceae and Ficus

Figure 3 .

 3 Figure 3. Ancestral state reconstruction of breeding systems in Moraceae. Here we show results from the maximum likelihood analysis of configuration A (five character states) with the equalrate model (for other results, see Figs. S3-S5). Clade labels: family names in Rosales, tribe names in Moraceae, and section names in Ficus.

(

  2016) used ITS and trnL-F, and Onstein et al. (2015) used eight molecular markers

FUNDING

  This work was supported by a China Scholarship Council PhD grant [grant numberINTRODUCTIONThe Neotropical ecozone has been defined as the region from central Mexico to southern Brazil[START_REF] Morrone | Biogeographical regionalisation of the Neotropical region[END_REF]. The Neotropics hold the highest terrestrial biodiversity on earth(Antonelli and Sanmartín 2011) and harbor all major tropical biomes: lowland rain forests, seasonally dry forests, mid-elevation montane forests, savannas, high elevation grasslands and deserts(Hughes et al. 2012). A recent estimate of tree species based on a pantropical tree inventory database suggested that the number of tree species in the Neotropics was as many as in the Indo-Pacific region and almost triple the counterparts in continental Africa[START_REF] Slik | An estimate of the number of tropical tree species[END_REF]. Several hypotheses have been proposed for the origins and evolution of Neotropical biodiversity during the past two decades. These hypotheses can be coarsely classified as biotic (e.g., dispersal ability, niche conservatism) and abiotic (time, climate, mountain uplift,Antonelli & Sanmartín, 2011).After completely splitting from continental Africa in the Cretaceous (ca. 105 Ma) (McLoughlin 2001), South America was isolated until the uplift of the Panama Isthmus (ca.15 Ma), which connected North and South America[START_REF] Montes | Arc-continent collision and orocline formation: Closing of the Central American seaway[END_REF] Bacon et al. 2015). Evidence of long-distance dispersals (LDD) among all of these three continents has been found with the development of molecular dating and phylogenetic approaches(Christenhusz and Chase 2013). The relationships between the floras of the Neotropics and of other continents has intrigued researchers, as investigating these relationships not only can shed light on the origin and evolution of Neotropical biodiversity, but also help to understand disjunct distributions and long-distance dispersal. Taxonomically and genetically densely sampled phylogenetic analyses represent an ideal approach to improve the understanding of the origin and evolution of Neotropical diversity(Antonelli and Sanmartín 2011; Hughes et al. 2012).The generic composition of the angiosperm tribe Dorstenieae in the Moraceae family has been under recent scrutiny. Based on taxonomic treatments and recent molecular phylogenetic studies, Dorstenieae are currently thought to consist of fifteen genera (

  prepared for other projects, including Brosimum (10), Trymatoccocus oligandrus, Helianthostylis sprucei, Allaeanthus luzonicus, and Malaisia scandens, as well as five samples (Artocarpus heterophyllus, Milicia excelsa, Parartocarpus venenosus, Ficus macrophylla, and Antiaropsis decipiens) from Johnson et al. (2016) and Zerega and

Figure 1 .

 1 Figure 1. Delimitation of distribution areas of Dorstenieae and relative dispersal matrices for the time-constant (model 0) and time-stratified (model 1) models.A, South America; B, North/Central America; C, Africa; D, Madagascar; E, India and Sri Lanka; F, Southeast Asia, and Oceania. Five levels of dispersal probability, 0.01, 0.25, 0.5, 0.75 and 1, representing the probability from low to high.

Figure 2 .

 2 Figure 2. Maximum Likelihood phylogenomic tree (a) and ASTRAL tree (b) of Dorstenieae. Fossil-calibrated nodes are marked as red stars. Nodes with bootstrap support value (less than 90% are indicated with an asterisk. Tip names are colored by general distribution area.

Figure 3 .

 3 Figure 3. Biogeographic history reconstruction of Dorstenieae based on the time-stratified DEC+J model. Inferred ancestral distribution areas prior to speciation are indicated on the nodes. Pie charts for selected nodes represent the relative probability (proportional likelihoods)

Figure 1 .

 1 Figure 1. Distribution map and species richness of Ficus according to breeding system, based on cleaned occurrence records from GBIF: A) dioecy; B) monoecy.

Figure 2 .

 2 Figure 2. Principal component analysis (PCA) of climatic variables in Ficus. Each dot represents the median value of a species. Bioclimatic variables are described inTable 1.

Figure 3 .

 3 Figure 3. Ancestral state reconstruction of precipitation of the warmest quarter (bio18) in Ficus. Breeding systems are highlighted next to the tips. The ancestral breeding system of Ficus and the number of transitions remain uncertain (Zhang et al. 2019).

  who discussed the potential biogeographic history of Moraceae based on their phylogeny and divergence time estimates. However, parametric biogeographic approaches, such as those used in Chapter II for Dorstenieae andPederneiras et al. (2018) for Ficus, have not yet been applied to Moraceae as a whole. Thus, it remains unclear whether Gondwana or Laurasia was the ancestral area of Moraceae. Reconstructing the biogeographic history of a pantropically distributed family such as Moraceae will improve the understanding of the origin of plant diversity in tropical regions.However, biogeographic model-based approaches are still relatively young and present some limitations that are currently being debated in the community. The widely used dispersalextinction-cladogenesis (DEC) model does not explore the scenario in which speciation occurs with dispersal events. This process might be common in scenarios of long-distance dispersal to a new environment, especially in island systems(Ree and Sanmartín 2009;[START_REF] Matzke | Model selection in historical biogeography reveals that founder-event speciation is a crucial process in island clades[END_REF]. The introduction of the founder-event speciation process (jump process) in the DEC model, leading to the now widely used DEC+J model, offered a solution to this problem(Matzke 2013). However, a recent study demonstrated that models which includes both anagenetic and cladogenetic

Figure 1 .

 1 Figure 1. A. Divergence time estimates for Brosimum by MCMCTree with default prior setting (time unit in 100 Ma).

Figure 1 .

 1 Figure 1. B. Divergence time estimates for Brosimum by MCMCTree with prior setting referring to the results from baseml (time unit in 100 Ma).

Figure 1 .

 1 Figure 1. C. Divergence time estimation for Brosimum by r8s (time unit in 1 Ma).

Figure S4 .

 S4 Figure S4. Ancestral state reconstruction with 320-species dataset by equal-rate maximum likelihood with configurations A to F (configuration B presented as Fig.3).

Figure S5 .

 S5 Figure S5. Ancestral state reconstruction with only Ficus species in the dataset with tip names by different approaches: A) parsimony; B) equal rate maximum likelihood; C) unequal rate maximum likelihood with equal root state prior; D) unequal rate maximum likelihood with equilibrium root state prior.

Titre:

  Phylogénie, biogéographie et évolution des systèmes sexuels chez les Moraceae Mots clés : Ficus, Dorstenieae, phylogénomique, datation moléculaire, reconstruction des états ancestraux, niche climatique Résumé : Les Angiospermes sont le clade le plus diversifié des plantes actuelles et sont exceptionnellement riches en espèces dans les régions tropicales. Dans cette thèse, j'ai étudié l'évolution des systèmes sexuels et l'histoire biogéographique de la famille des Moraceae, clade modèle utilisé pour comprendre l'origine et l'évolution de la diversité chez les Angiospermes. Dans le Chapitre I, j'ai reconstruit et calibré un nouvel arbre phylogénétique daté pour les Moraceae. J'ai ensuite utilisé cet arbre pour reconstruire les états ancestraux des systèmes sexuels chez les Moraceae et Ficus. Les âges des groupes-couronne des Moraceae et du genre Ficus sont estimés au Crétacé et à l'Eocène, respectivement. La dioécie est inférée comme l'état ancestral des systèmes sexuels chez les Moraceae, avec plusieurs transitions ultérieures vers la monoécie, y compris chez Ficus. Ce résultat suggère que la dioécie ne représente pas nécessairement un cul-de-sac évolutif. Dans le Chapitre II, j'ai reconstruit un arbre phylogénétique daté pour la tribu des Dorstenieae, distribuée principalement dans les régions tropicales, à partir d'un nouveau jeu de données génomiques nucléaires produit avec une approche Hyb-Seq. L'histoire biogéographique du groupe a ensuite été reconstruite en utilisant les modèles de dispersion-extinction-cladogenèse. Les âges des groupescouronne des Dorstenieae et du genre Dorstenia sont estimés au Crétacé et dans la période du Crétacé au Paléocène, respectivement. Deux évènements de dispersion à longue distance depuis l'Afrique continentale vers l'Amérique du Sud ont eu lieu au Cénozoïque (Dorstenia et Brosimum s.l.). Dans le Chapitre III, j'ai testé les différences de niche climatique (température et précipitation) entre les deux systèmes sexuels (monoécie et gynodioécie) chez Ficus avec un nouveau jeu de données fiables d'occurrences spatiales et de systèmes sexuels chez 183 espèces. À cette fin, j'ai utilisé deux approches comparatives : équations d'estimation généralisées (GEE) et modèles linéaires généralisés (GLM). Une relation positive entre précipitation et gynodioécie est soutenue par les analyses GLM, et aucune méthode ne soutient une relation entre température et système sexuel. Une meilleure capacité à se disperser et le potentiel d'autopollinisation sont deux explications possibles pour la colonisation et la survie des espèces monoïques dans des environnements plus secs. Cette thèse démontre le potentiel des méthodes phylogénétiques comparatives et des données phylogénomiques pour répondre aux questions d'évolution des systèmes sexuels et de biogéographie chez les Moraceae et ouvre plusieurs nouvelles perspectives importantes méritant d'être approfondies chez d'autres clades de plantes, telles que la relation entre système sexuel et niche climatique. Title: Phylogeny, biogeography, and breeding system evolution in Moraceae Keywords: Ficus, Dorstenieae, phylogenomics, molecular dating, ancestral state reconstruction, climatic niche Abstract: Angiosperms are the most diversified clade of extant plants and are exceptionally species-rich in tropical regions.

CHAPTER I. ESTIMATING DIVERGENCE TIMES AND ANCESTRAL BREEDING SYSTEMS IN FICUS AND MORACEAE

	This chapter was recently published in Annals of Botany in January 2019
	(https://academic.oup.com/aob/article/123/1/191/5092035), in collaboration with Stefan
	Little, Renske Onstein and Hervé Sauquet.
	. In addition, this work led to a
	collaboration on a smaller phylogenomic study focused on the genus Brosimum. While this
	has not yet reached the stage of a draft manuscript, I have outlined my contribution in
	Supplementary Data part of this thesis. Chapter III focused on characterizing niche
	differences of the two breeding systems in Ficus. To do so, I built a spatial occurrence

Table 1 .

 1 Summary of the six character state configurations used to reconstruct ancestral breeding systems in Moraceae.

	Configuration	States	Justification
	A	Bisexuality, gynodioecy,	Original states recorded in species
		androdioecy, monoecy,
		dioecy
	B	Bisexuality, gynodioecy,	Excluding state androdioecy for its low frequency in
		monoecy, dioecy	the dataset of this study
	C	Bisexuality, monoecy,
		dioecy + gynodioecy

Based on configuration B, gynodioecy combined with dioecy: all gynodioecious species (but one) in this study belong to Ficus, where they are functionally dioecious D Gynodioecy, monoecy, dioecy Based on configuration B, excluding state bisexuality: only distant outgroups of Moraceae include species with bisexual flowers (e.g., Rosaceae, Elaeagnaceae, Rhamnaceae) E Monoecy, dioecy + gynodioecy Based on configuration B, gynodioecy combined with dioecy (as in C) and bisexuality excluded (as in D) F Gender monomorphic, gender dimorphic

Table 2 .

 2 Estimated divergence times for key nodes of Moraceae. Here, we follow the tribal classification of Clement and Weiblen (2009), as shown in Fig 1.

	Nodes	Breeding systems	Support values	Mean age (95%
			(BP/PP) 1	HPD) 2 (Ma)
	CG 3 Rosales	Bisexuality, monoecy, androdioecy,	/	105.5 (102.8-107.0)
		gynodioecy, dioecy
	SG 4 Moraceae	/	100/1.00	87.5 (81.7-93.3)
	CG Moraceae	Monoecy, androdioecy, gynodioecy,	99/1.00	79.0 (73.2-84.7)
		dioecy
	CG Artocarpeae	Monoecy, dioecy	99/1.00	65.6 (64.0-68.6)
	CG Moreae 5	Monoecy, dioecy	91/1.00	40.3 (34.9-46.8)
	CG Maclureae	Dioecy	98/1.00	24.7 (8.9-41.1)
	CG Dorstenieae	Monoecy, dioecy	98/1.00	60.5 (51.4-70.2)
	CG Castilleae	Monoecy, androdioecy, dioecy	100/1.00	31.2 (18.7-47.0)
	SG Ficeae	/	98/1.00	62.2 (56.0-68.6)
	CG Ficeae	Monoecy, gynodioecy	99/1.00	48.5 (40.6-55.9)
	CG Pharmacosycea	Monoecy	58/1.00	11.6 (5.5-19.4)
	1
	CG Pharmacosycea	Monoecy	100/1.00	24.2 (13.2-36.7)
	2
	CG Sycomorus	Monoecy, gynodioecy	77/1.00	38.7 (28.0-41.1)
	CG Sycidium	Gynodioecy	100/1.00	29.2 (22.6-36.4)
	CG Synoecia	Gynodioecy	100/1.00	20.1 (13.7-26.8)
	CG Ficus	Gynodioecy	82/1.00	31.0 (24.7-38.0)
	CG Urostigma 1	Monoecy	100/1.00	21.5 (12.0-31.2)
	CG Urostigma 2	Monoecy	99/1.00	35.5 (28.0-42.3)
	1 PP: posterior probability; BP: bootstrap probability; 2 95% HPD: 95% highest posterior density; 3 CG:
	crown group;

4 SG: stem group; 5 Streblus smithii was here excluded from crown-group Moreae (see text for details).

Table 3 .

 3 Summary results of ancestral state reconstruction for the complete (320 species) dataset by maximum likelihood approach.

	CG Moraceae 1	Prob 2	CG Ficus 3	prob -Lh 4 Nr	AIC
					of
					par 5
	Equal rate
	A Dioecy	0.91	Monoecy	0.92 121.6 10	263.2
	B Dioecy	0.90	Monoecy	0.91 102.7 6	217.4
	C Gynodioecy + dioecy	0.95	Monoecy	0.87 91.5	3	189
	D Dioecy	0.88	Monoecy	0.89 78.9	3	163.8
	E Gynodioecy + dioecy	0.92	Monoecy	0.80 63.1	1	128.2
	F Gender dimorphic	0.89	Gender	0.76 76.1	1	154.2
			monomorphic
	Unequal rate
	A Androdieocy; dioecy;	Each ca.	Gynodioecy	0.95 124.7 20	289.4
	gynodioecy	0.33
	B Dioecy	0.99	Gynodioecy	0.98 90.8	12	205.6
	C Gynodioecy + dioecy	0.99	Gynodioecy +	0.84 76.9	6	165.8
			dioecy
	D Dioecy	1.00	Gynodioecy	1.00 70.0	6	152
	E Gynodioecy + dioecy	1.00	Gynodioecy +	0.83 61.6	2	127.2
			dioecy
	F Gender dimorphic	1.00	Gender dimorphic 0.84 74.6	2	153.2

1 CG Moraceae: breeding system state of crown-group Moraceae; 2 Prob: probability of support estimated by maximum likelihood; 3 CG Ficus: breeding system state of crown-group Ficus; 4 Lh: minus loglikelihood; 5 Nr of par: number of free parameters; The best-fit model identified by AIC for each configuration is indicated in boldface.

  These results suggest that the unequal rate model is overparameterized for this configuration with five character states. Indeed, twenty free parameters (transition rates) are estimated in this model. In contrast, a four-state unequal rate model requires 12 parameters, a three-state model requires six parameters, and a binary model requires only two parameters. The potential overparameterization of configuration A was also suggested by the results of reversible-jump Markov Chain Monte Carlo (rjMCMC) analyses. For each configuration, the number of free parameters in the model was estimated between one and three.

Table 4 .

 4 Summary

			One MCC tree		3600 posterior trees
	Configuration CG	95%	CG	95%	CG	95%	CG Ficus 95%
		Moraceae 1	HPD 2	Ficus 3	HPD	Moraceae	HPD		HPD
	A	Dioecy	0.82-	Monoecy 0.00	Dioecy	0.80-	Monoecy	0-
			1.00		-0.88		1.00		0.99
	B	Dioecy	0.84-	Monoecy 0.71-	Dioecy	0.80-	Monoecy	0.71-
			1.00		0.90		1.00		1.00
	C	Gynodioecy	0.86-	Monoecy 0.00-	Gynodioecy	0.74-	Monoecy	0.00-
		+ dioecy	1.00		0.84	+ dioecy	1.00		0.97
	D	Dioecy	0.85-	Monoecy 0.66-	Dioecy	0.82-	Monoecy	0.65-
			1.00		0.90		1.00		1.00
	E	Gynodioecy	0.87-	Monoecy 0.00-	Gynodioecy	0.75-	Monoecy	0.00-
		+ dioecy	1.00		0.85	+ dioecy	1.00		0.98
	F	Gender	0.85-	Gender	0.19-	Gender	0.76-	Gender	0.03-
		dimorphic	1.00	dimorphic	1.00	dimorphic	1.00	Dimorphic	1.00
	1 CG Moraceae: breeding system state of crown-group Moraceae; 2 95% HPD: 95% highest posterior density;

results of ancestral state reconstruction for the complete (320 species) dataset obtained with the Bayesian approach (reversible-jump MCMC).

3 CG Ficus: breeding system state of crown-group Ficus.

  stem node of the clade with which they share apomorphies(Sauquet 2013) (TableS2). For instance, the fossil achenes attributed to Ficus(Collinson 1989) were here used to calibrate the stem node of Ficus, whereas[START_REF] Rønsted | 60 million years of co-divergence in the fig-wasp symbiosis[END_REF] and Cruaud et al.(2012) used the same fossils to calibrate the crown node of Ficus. We also used different prior distributions for fossil calibrations. Although some authors have suggested that priors such as lognormal or exponential distributions may help to improve the accuracy of divergence time estimation(Yang and Rannala 2005; Ho and Phillips 2009), the prior

new time scale presented here for Moraceae was estimated with 12 fossil calibrations (four ingroup, eight outgroup), more than in any previous study of the family so far. Here the crown-group age of Moraceae was estimated as 73.2-84.7 Ma. This estimate is similar to that of

Zerega et al. (2005) (72.6-110.0 Ma)

, but with a narrower confidence interval. Our stem-group age estimate for the family

is older than that reported by

Magallón et al. (2015) (57.7-77.8 Ma)

, most likely because our increased sampling of the family and its outgroups allowed us to use more fossil calibrations in Rosales. The crown-group age of Ficus was here estimated as 40.6-55.9 Ma, which is similar to the age reported by

Zerega et al. (2005)

(40.1-51.0 Ma), but younger than the age found by Rønsted et al. (2005) (51.4-78.0 Ma) and Cruaud et al. (2012) (60.0-

101.9 Ma). These differences may be explained by different calibration strategies. When fossil species have similar morphological traits with extant species of a clade, these fossils are often used optimistically to calibrate the crown group of this clade. This practice is problematic because it rules out the possibility that such fossils are stem relatives of the clade

(Ronquist et al. 2012; Sauquet 2013

). Here we used an apomorphy-based approach, using fossils as minimum age constraints modeled with uniform prior distributions to calibrate the setting of parameters such as mean or standard deviation for these distributions is often arbitrary and the shape of these distributions implicitly assumes that the fossil diverged close to the node calibrated. Our strategy here was instead to use uniform distribution with fossil ages as the minimum bound and maximum root age as the maximum bound. The crown-group ages of tribes Artocarpeae and Dorstenieae were estimated as 64.0-68.6 Ma and 51.4-70.2 Ma, respectively, both of which are younger than the ages reported by

Williams et al. (2017)

for Artocarpeae

and

Table 1 .

 1 Table1) and approximately 156 species(Berg et al. 2001; Clement & Weiblen 2009; Zerega et al. 2010;[START_REF] Chung | Molecular recircumscription of Broussonetia (Moraceae) and the identity and taxonomic status of B. kaempferi var. australis[END_REF]. The most diverse genus in Dorstenieae is Dorstenia, which includes 113 species(Berg and Hijman 1999; Mccoy and Massara 2008; dos Santos and Neto 2012; Machado and Marcelo Filho 2012;[START_REF] Chase | Dorstenia christenhuszii (Moraceae), a new species from the Taita Hills, Kenya[END_REF]; Classification of Dorstenieae and the number of species for each genus.

	Estimating absolute divergence times as accurately as possible is essential to
	biogeographic reconstruction to connect the evolutionary history of target taxa with
	past climate change and geographic events (Sauquet 2013). Previous studies have
	estimated the crown-group age of Dorstenieae to be at least 71 Ma (Zerega et al. 2005),

dos

Santos et al. 2013; Leal 2014; Machado et al. 2014; Rzepecky 2016)

. The genera in Dorstenieae are restricted to either side of the Atlantic, except Dorstenia, which has almost the same number of species in South America and continental Africa (50 in the Neotropics; 62 in continental Africa, Madagascar and Arabian Peninsula; 1 in India and Sri Lanka)

(Berg and Hijman 1999)

. In a recent phylogenetic study of Moraceae

(Zhang et al. 2019)

, the relationships among the genera of Dorstenieae (except Bosqueiopsis and Scyphosyce) were reconstructed. Most of them were strongly supported, but the relationships of Trilepisium, Brosimum and Treculia remained unclear. The most densely sampled phylogenetic study of Dorstenia to date sampled 32 species (28%) and found Dorstenia to have originated in Africa, with three African species nested inside the Neotropical clade

(Misiewicz and Zerega 2012)

. 50.6-72.5 Ma

(Gardner et al. 2017), or 51.4-70.2 Ma (Zhang et al. 2019)

. While these estimates overlap, the crown-group age of Dorstenia has remained unclear. The crown-

  two Neotropical species included), 12.7-31.7 Ma[START_REF] Chen | Latitudinal effects on phenology near the northern limit of figs in China[END_REF] with two neotropical and one African species included), or 84.8-132.0 Ma(Misiewicz & Zerega, 2012, with 15 neotropical and 14 African species included).Approximately 10% of the species in Moraceae are herbaceous and all of them belong to Dorstenieae (more specifically all in Dorstenia and Fatoua;Berg, 2001).Species of Dorstenieae show diversity in pollination modes, dispersal mechanisms, and habit(Berg 2001). The species of Dorstenia are found in a wide variety of habitats (e.g., tropical rain forest, savannas, or crevices of cliffs) and life forms (e.g., tree, shrubs,

caulescent, herbaceous;

Berg & Hijman, 1999; Berg, 2001; Misiewicz & Zerega, 2012)

.

  , which uses reference sequences to guide local de novo assemblies of each target gene. Because Dorstenieae are phylogenetically distant from the baits sequences, which come from

	Artocarpeae and Moreae (Gardner et al. 2016), a new HybPiper reference was
	generated using six samples sequenced here, each with at least eight million read pairs
	(Dorstenia bahiensis, D. cayapia, D. erythrandra, D. kameruniana, Treculia africana,
	and Fatoua villosa). Reads were trimmed using Trimmomatic (Bolger et al., 2014)
	(LEADING:20 TRAILING:20 SLIDINGWINDOW:4:20 MINLEN:20) and assembled
	with SPAdes 3.10.1 (Bankevich et al. 2012) using default parameters. Coding
	Assembly of all reads then proceeded as follows. We first trimmed low quality
	bases and adapter sequences using Trimmomatic (Bolger et al., 2014) with the
	following parameters: ILLUMINACLIP:TruSeq3-PE-2.fa:2:30:10 LEADING:20
	TRAILING:20 SLIDINGWINDOW:4:20 MINLEN:40. To ensure that quality
	trimming worked as expected, we examined a subset of reads were both before and
	after	trimming	with	FastQC	(https://www.bioinformatics.babraham.ac.uk/
	projects/fastqc/). Reference-guided assembly then proceeded with HybPiper (Johnson
	et al. 2016), using the new reference described above. Briefly, the program works as

sequences were predicted with Augustus

[START_REF] Stanke | AUGUSTUS: a web server for gene finding in eukaryotes[END_REF]

, using Arabidopsis thaliana genes as a reference. A seven-way orthology search was carried out with ProteinOrtho5

[START_REF] Lechner | Proteinortho: Detection of (Co-)orthologs in large-scale analysis[END_REF]

) using all CDS over 200 bp from the de novo assemblies in addition to the Artocarpus HybPiper

[START_REF] Johnson | HybPiper: extracting coding sequence and introns for phylogenetics from high-throughput sequencing reads using target enrichment[END_REF]

) reference from Kates et al.

[START_REF] Kates | Allele phasing has minimal impact on phylogenetic reconstruction from targeted nuclear gene sequences in a case study of Artocarpus[END_REF]

. Orthologs present in at least three taxa were included in a new seven-taxon HybPiper reference consisting of in-frame CDS. This expanded reference contained approximately 500 genes.

follows: reads are sorted by gene based on the reference. Local de novo assemblies are carried out using SPAdes

[START_REF] Bankevich | Evaluating the impact of genomic data and priors on Bayesian estimates of the angiosperm evolutionary timescale[END_REF]

, and coding DNA sequences (CDS)

  crown-group age of tribe Dorstenieae was estimated in the and that of Dorstenia in the Paleocene (50.8-67.3 Ma) (Table2, FigureS2). The stem and crown-group ages of the Neotropical Dorstenia clade were dated in the Eocene to early Oligocene (34.6-51.8 Ma and 29.8-44.7 Ma, respectively). Brosimum s.l., the other Neotropical clade in Dorstenieae, was estimated to date from the late Eocene to early Miocene (stem node: 19.6-41.5 Ma; crown node: 14.9-31.1 Ma). Runs with the prior only showed different results to those with data, indicating the data had a significant impact on the posterior (results not shown). Results of the two different prior settings were similar to each other (Table2, Figure

Table 2 .

 2 Divergence time estimates with penalized likelihood (PL with r8s) and Bayesian (MCMCTree) approaches for key nodes of Dorstenieae.

	Node	r8s
		(Ma)

MCMCTree set1 (Ma) MCMCTree set2 (Ma)

  Dorstenia species distributed in Central and North America, from D.

	SG Dorstenia MAm clade1	23.1	16.5-28.2	16.7-28.0
	CG Dorstenia MAm clade1	20.2	12.6-24.9	12.8-24.6
	SG Dorstenia MAm clade2	25.4	23.9-35.1	24.1-34.7
	CG Dorstenia MAm clade2	13.3	12.1-28.9	12.0-28.4
	CG: crown group; SG: stem group
	MCMCTree set1: results of estimate in MCMCTree with the prior setting referring to the results of
	baseml
	MCMCTree set2: results of estimate with the default prior setting in MCMCTree
	Brosimum s.l. includes Brosimum, Trymatococus and Helianthostylis
	Broussonetia s.l. includes Broussonetia, Malaisia, and Allaeanthus
	Dorstenia Neo: Neotropical Dorstenia species
	Dorstenia MAm Clade1:
	SG Dorstenieae	77.3	68.7-82.1	68.9-82.1
	CG Dorstenieae	75.4	65.2-79.8	65.8-79.8
	CG Treculia	35.3	10.9-32.7	11.2-32.4
	SG Brosimum s.l.	40.7	19.4-42.9	19.6-41.5
	CG Brosimum s.l.	27.9	14.7-31.7	14.9-31.1
	SG Broussonetia s.l.	70.4	60.2-76.4	61.1-76.2
	CG Broussonetia s.l.	52	30.5-41.8	30.2-40.7
	CG Scyphosyce	34.7	10.8-31.5	10.7-30.9
	SG Dorstenia	61.9	51.2-68.5	52.3-68.7
	CG Dorstenia	61.3	49.8-67.0	50.8-67.3
	SG Dorstenia Neo	43.4	34.3-52.1	34.6-51.8
	CG Dorstenia Neo	33.4	29.7-44.7	29.8-44.7

erythrandra to D. caimitensis in Figure

3

Dorstenia MAm Clade2: Dorstenia species distributed in Central and North America, from D. excentrica to D. contrajerva in Figure

3

Table 1 .

 1 Estimates for interaction effects of climatic variables and breeding systems in Ficus from the regressions by generalized linear models (GLM) and generalized estimating equations (GEE).

	Code	Variable	GLM	GLM coefficient	GEE
	bio1	Annual mean temperature	not sig.	/		not sig.
	bio5	Max temperature of warmest month		-0.022318
			**			not sig.
	bio6	Min temperature of coldest month	not sig.	/		not sig.
	bio8	Mean temperature of wettest quarter	not sig.	/		not sig.
	bio9	Mean temperature of driest quarter	not sig.	/		not sig.
	bio10	Mean temperature of warmest quarter	not sig.	/		not sig.
	bio11	Mean temperature of coldest quarter	not sig.	/		not sig.
	bio12	Annual precipitation	***	0.10051		not sig.
	bio13	Precipitation of wettest month	***	0.2675		not sig.
	bio14	Precipitation of driest month	***	0.23656		not sig.
	bio16	Precipitation of wettest quarter	***	0.17727		not sig.
	bio17	Precipitation of driest quarter	***	0.13509		not sig.
	bio18	Precipitation of warmest quarter	***	0.28233		not sig.
	bio19	Precipitation of coldest quarter	*	0.04934		not sig.
	PC1	Principal component 1	***	-0.2242		not sig.
	PC2	Principal component 1	***	-0.29111		not sig.
	PC1+2	Combination of principal component 1 and 2		-0.25431	(PC1),	-
			***	0.34284 (PC2)	not sig.

Significance codes: '***' P<0.0001; '**' P<0.001; '*' P<0.01; '.' P<0.05; 'not sig.' not significant.

  , herbarium specimens are now a promising resource for phylogenomic work using approaches such as target enrichment, where DNA fragmentation (typical of herbarium material) is much less of an issue than traditional Sanger sequencing of PCR-amplified markers. The rapid adoption of this and other genomic approaches by the plant systematic community is likely to lead to considerable improvement of phylogeny at all levels[START_REF] Johnson | A universal probe set for targeted sequencing of 353 nuclear genes from any flowering plant designed using k-medoids clustering[END_REF] Couvreur et al. 2019), including in Moraceae(Rasplus et al. 2018;

Zerega and Gardner 2019)

. These developments may in turn result in further necessary adjustments of the current phylogenetic classification of Moraceae

(Clement and Weiblen 2009; Zerega and Gardner 2019

) and the circumscription of several genera. For instance, in Chapter II, I found the genera Scyphosyce and Utsetela to be nested in Dorstenia and Trymatococcus and Helianthostylis in Brosimum s.l. The latter result was also confirmed in the more densely sampled side study focusing on Brosimum (see supplementary data).

Table S2 .

 S2 Detailed information for the fossil calibrations used in this study.

	SG Dorstenia angusticornis Morac Artocarpoxylon deccanensis Dorstenie Xylodorstenia wood Dorstenia drakeana Dorstenie Dorstenia Dorstenia ramosa subsp. Dorstenie Dorstenia	India MO25064 / V018140	9,459,4 6,289,0 6,955,0	Mehrotra et al., 1984 42.64% 45.96% 0.93%	Williams et al., 2017 102 08/06/1995 102 09/10/1987 84 s.d.	64-67 Ma BO Daramola 665(MO) 64.0 SD Koch et al 87195(NY) M.D.M. Vianna Filho 2020(F) Hooper et al., 2010,
	Artocarp dolichocaula eae	Mehrotra, Prakash, and Bande ae ae ae		15 1F		37 82 51				Selena Smith, pers
	us sl Dorstenia appendiculata Dorstenia elliptica Dorstenia roigii	Dorstenie Dorstenie Dorstenie	Lecanium Nothodorsteni Emygdioa	NY MO25064 /	11,461, 1,984,6 331,405	0.34% 5.33% 35.23%	57 25/01/1996 95 28/01/1997 30 07/11/1923	com Andre MA Amorim et al. 1941(NY) J. Nning 20(MO) Ekman 17973(F)
	SG Rhamna ceae Dorstenia arifolia Rham nacea e Dorstenia erythrandra Dorstenia soerensenii CG Rhamna ceae Rham nacea Dorstenia arifolia-1 Dorstenia excentrica Dorstenia tayloriana var. Coahuilanthus belindae Calvillo-Canadell & Cevallos-ae ae a ae Dorstenie Dorstenia Dorstenie Emygdioa Dorstenie Lomathophora flowers ae ae ae Ferriz Paliurus clarnensis Burge & Manchester fruits Dorstenie Dorstenia Dorstenie Emygdioa Dorstenie Lomathophora ae ae tayloriana ae Dorstenia arifolia-2 Dorstenie Dorstenia Dorstenia excentrica-1 Dorstenie Emygdioa Dorstenia tenera Dorstenie Lomathophora e CG Rosacea Rosac eae Prunus cathybrownae Benedict, DeVore & Pigg flowers+ 565961 Coahuila, Mexico 715 00 35 / 10,830, / 14,844, MO25064 2,189,3 558 492 03 29 Red Gap, Oregon, USA V020231 2,191,7 / 10,039, / 8,308,2 8F 07 982 20 V020231 1,350,1 / 7,735,0 P0682282 3,333,6 fruits Boot Hill locality, ae 9F 71 ae 69 ae 7 99 Dorstenia astyanactis Dorstenie Kosaria P0682277 10,658, Dorstenia fawcetti Dorstenie Emygdioa / 5,545,8 Dorstenia tenuis Dorstenie Emygdioa V022835 1,542,6 Washington, USA ae 9 870 ae 92 ae 1F 31 e SG Ulmacea Ulmac eae Ulmites ulmifolius (Schloemer-Jäger) Kvacek leaves Spitsbergen, Antiaropsis descipiens Castilleae 1,151,7 Dorstenia bahiensis Dorstenie Dorstenia P0004758 14,590, Dorstenia flagellifera Dorstenie Emygdioa / 10,711, Dorstenia turbinata Dorstenie Xylodorstenia MO25064 3,233,7 23 ae 7 733 ae 580 ae 06 05 Norway Allaeanthus greveana Dorstenie / 3,425,7 Dorstenia barnimiana Dorstenie Acauloma MO25064 9,790,2 Dorstenia foetida Dorstenie Kosaria / 13,237, Dorstenia umbricola Dorstenie Lecanium NY 4,586,2	Calvillo-Canadell and Cevallos-Ferriz, 2007 33.04% 22.60% 13.48% Burge and Manchester, 2008 19.26% 5.70% 29.65% 23.70% 46.05% 4.96% Benedict et al., 2011 57.88% 9.95% 14.55% 21.87% Kvacek et al., 1994 5.37% 17.61% 39.13% 4.00% 47.68% 19.93% 1.32% 21.23%	Calvillo-Canadell and Cevallos-Ferriz, 2007 102 04/04/2002 102 1927August 92 21/12/1965 Burge and Manchester, 2008 98 04/03/1980 102 25/10/1981 94 28/01/2002 66 12/04/1987 102 10/08/1991 97 14/04/1990 Benedict et al., 2011 97 15/05/2014 102 23/06/2010 102 31/01/1980 102 14/12/1984 Manchester, 1999 99 102 11/07/1993 101 1927 July 96 14/10/2002 101 84 1981 spring 102 25/05/1963 84 24/02/1989	Late Campanian AM de Carvalho 7169(NY) 72.1 WJ Eyerdam 298(F) WJJO de Wilde 9392(MO) 44 Ma (middle Eocene) 44.0 J.P.P. Carauta 3423(F) M Nee 22374(F) MA Mwangoka & A Kalage 2673(MO) Calvillo-Canadell and Cevallos-Ferriz, 2007 Burge and Manchester, 2008 Amorim, A.M. 7(F) Liborio Lopez 6(F) F.J. Breteler 9984(P) 49.42 ± 0.54 Ma 49.42 Benedict et al., 2011 E. Gardner 98 (SAN) Haba, P.K. 677(P) CC Berg 990(NY) S.G. Tressens, E. Cabral, S. Cáceres, M. Urbani & C. Zamudio 2911(F) Paleocene 56.0 N. Zerega 281 (NY) H174(P) W.J. Eyerdam 196(F) M Cheek 11086(MO) Manchester, 1999 O. Pascal and F. Hallé 677 (NY) J Lavranos & MB Bleck 19480(MO) William Burge 2844(F) M. Alexiades, V. Pesha(NY)
	e			ae ae ae ae			14 196401	25 93 026 73
	SG Allaeanthus kurzii Urtica Dorstenia benguellensis Boehmeria sibirica Dorofeev Dorstenie Dorstenie Kosaria Dorstenia goetzei Dorstenie Kosaria Dorstenia urceolata Dorstenie Lecanium	achenes / MO25064 Kireevskoe, USSR 913,251 4,419,7 / 4,806,8 V022835 8,823,4	Collinson, 1989 29.55% 38.60% 17.79% 20.96%	Collinson, 1989 100 101 05/12/1998 98 06/08/1997 102 16/10/1977	Oligocene and Al Gentry and Chawalit Niyomdham 66546 (MO) 23.03 Collinson, 1989 AS Mkeya et al 1109(MO) PB Phillipson 4801(MO) P.J.M. Maas & P. Carauta 3233(F)
	Boehme	ceae		ae ae ae ae			22 4F		56 38 36			Miocene
	ria Allaeanthus luzonicus Dorstenia bonijesu Dorstenia hildebrandtii var. Dorstenia variifolia	Dorstenie Dorstenie Dorstenie Dorstenie	Dorstenia Kosaria Lomathophora	P0084063 MO25064 MO25064	76,606 4,914,4 5,145,9 1,620,3	71.24% 10.77% 5.19% 32.25%	98 102 10/05/1981 100 12/04/2000 47 10/05/1987	boundary Kuo-fang Chung 2016 (HAST) V.F. Ferreirov 1729(P) L Festo 648(MO) J&J Lovett et al. 2135(MO)
	SG Pilea schlechteri Bleekrodea Dorstenia bowmaniana Urtica ceae Pilea lithuanica Dorofeev ae ae ae ae Dorstenie Dorstenie Lecanium Dorstenia hirta Dorstenie Lecanium Dorstenia warneckei Dorstenie Kosaria madagascariensis ae ae ae ae Bosqueiopsis gilletii Dorstenie Dorstenia brasiliensis Dorstenie Emygdioa Dorstenia holstii Dorstenie Lomathophora fruits Dorstenia yambuyaensis Dorstenie Lomathophora	USSR P0686294 0 28 25 V020230 NY P0681971 1 6F 777756 9 P0682285 / / P0681971	77 77 09 10,190, 8,048,9 4,311,8 2,582,2 180 50 69 28 8,006,4 3,508,8 8,494,8 1,644,6	Collinson, 1989 42.29% 27.98% 35.57% 1.17% 21.83% 21.49% 11.47% 4.22%	Collinson, 1989 102 24/11/2009 102 22/01/1971 101 17/07/2000 50 08/09/1984 102 15/11/2009 89 13/07/1995 95 08/12/2000 57 22/09/1983	Oligocene and MYA 375(P) P. Carauta 1277(F) 23.03 FO Souza et al 21(NY) D.W. Thomas 3660(P) Miocene boundary J.R. Timberlake, T. Müller & F. Crawford 5767(P) Collinson, 1989 MJ Jansen-Jacobs et al 4436(NY) A Ntemi Sallu 645(MO) Mikio Kaji 24(P)
				ae ae ae ae			1 8		93 61 98 28
	Brosimum acutifolium Dorstenia brasiliensis-1 Dorstenia indica Dorstenia zanzibarica	Dorstenie Dorstenie Dorstenie Dorstenie	Emygdioa Kosaria Kosaria		V022835 / MO25064	377,383 5,416,0 2,795,9 1,249,3	26.38% 28.00% 23.85% 19.54%	102 102 24/11/1974 95 15/12/1975 73 04/08/1986	T. Lileen 4457 L.C. Gurken 37(F) L Bernardi 16022(MO) J Lovett 895(MO)
				ae ae ae ae			2F 23		62 44 45
	Brosimum alicastrum Dorstenia brasiliensis-3 Dorstenia kameruniana Fatoua villosa	Dorstenie Dorstenie Dorstenie Dorstenie	Emygdioa Nothodorsteni	V020232 /	235,089 4,612,3 15,998, 13,261,	64.90% 1.99% 9.27% 38.91%	102 7/07/2013 102 11/02/1993 102 21/08/2003 100 10/07/2013	E. Gardner 23 (CHIC) G. Hatschbac 58895(F) OA Kibure & H Bofu 1045(MO) E. Gardner 27 (CHIC)
				ae ae ae ae	a		9F		29 183 729
	Brosimum costaricanum Dorstenia brownii Dorstenia lindeniana Helianthostylis sprucei	Dorstenie Dorstenie Dorstenie Dorstenie	Lomathophora Lecanium	MO25064 /	584,789 3,049,3 7,149,7 201,836	50.02% 5.53% 37.24% 39.65%	102 31/01/1988 92 10/11/2004 102 26/04/1966 101	C. Kernan 27 (F) MA Mwangoka et al 3606(MO) Gayle C Jone et al 3130(F) F16142
	MRCA Brosimum gaudichaudii Family Taxon Dorstenia buchananii Dorstenia lindeniana-1 Malaisia scandens Brosimum guianense Dorstenia caimitensis Dorstenia mannii Scyphosyce manniana	ae ae ae ae Dorstenie Dorstenie Dorstenie Dorstenie ae ae ae ae Dorstenie Dorstenie Dorstenie Dorstenie	Fossil Lomathophora Kosaria Lecanium Emygdioa	Locality MO95910 20 / 7 / P0682277 P0682500	65 62 62,019 833,457 4,738,6 499,782 96 313,169 2,090,3 3,093,5 531,584	Reference (description) 58.54% 41.67% 16.39% 64.65% 37.28% 2.24% 2.12% 17.92%	Reference (node assignment) 81 37 04/12/1954 102 28/11/2002 102 25/05/2014 102 99 23/08/1927 62 07/12/1999 84 13/01/1987	Oldest fossil age S. Mendes & al. 332 (F) Time NC Chase 5466(F) D Alvarez 2721(MO) E. Gardner 122 (SAN) set (Ma) J.F. Morales 5258 WJ Eyerdam(F) Simons, E.L.A.N. & R. Westerduijn 300(P) Age reference H.F. 1220(P)
	CG	Canna	Aphananthe cretacea ae ae ae ae		fruits	3 5	Walbeck,	62 59	Knobloch and Mai,	Friis et al., 2011	Maastrichtian	66.0	Magallón et al., 2015
	Cannaba Brosimum multinervium bacea Knoblock & Mai Dorstenie Dorstenia carautae Dorstenie Dorstenia milaneziana Dorstenie Scyphosyce pandurata Dorstenie	Lecanium Lecanium	(endoca MO14913 Germany NY P0682497	269,526 2,804,1 7,146,0 10,995,	1986 40.27% 24.69% 23.36% 29.27%	101 102 05/01/1990 102 17/03/1996 102 1987	Jaramillo 574 JM Silva 782(MO) WW Thomas, A Amorin, J Jardim 11078(NY) D.W. Thomas 6869(P)
	ceae	e		ae ae ae ae		rp)	75 95156 7	60 35 555
	SG Brosimum parinarioides Canna Humulus rolundatus Dorofeev Dorstenie Dorstenia cayapia subsp. Dorstenie Emygdioa Dorstenia oligogyna Dorstenie Nothodorsteni fruits Sloetia elongata Dorstenie	Isakovka, USSR 360,342 20,545, P0682275 NY 6,466,6 P0688953 6,790,7	Takhtajan, 1982 43.65% 31.03% 12.44% 37.13%	Collinson, 1989 102 21/01/1983 102 27/02/1999 102 31/10/1994 102 30/12/1992	Eocene and Vasquez and Jaramillo 3817 (F) 33.9 Collinson, 1989 Andre Ma Carvalho et al 6710(NY) J.J. Wieringa 3005(P) P. Thomas & L.E. Teo no number(P)
	Humulus cayapia Brosimum potabile bacea e Dorstenia cayapia subsp. Dorstenia peruviana Treculia acuminata SG Ficus Morac eae Ficus lucidus Chandler ae ae ae ae Dorstenie Dorstenie Dorstenie Dorstenie ae vitifolia ae ae ae Brosimum rubescens Dorstenie Dorstenia cayapia subsp. Dorstenie Dorstenia petraea Dorstenie Treculia africana Dorstenie SG Brousso netia Morac eae Broussonetia rugosa Chandler a Emygdioa Lecanium Emygdioa Emygdioa ae vitifolia-1 ae ae ae Brosimum utile Dorstenie Dorstenia ciliata Dorstenie Lomathophora (endoca rp) 585795 8 9 / MO14087 P0676332 achenes Southern England 986 87 11 466,812 2,952,9 2,200,4 210,985 11 84 76 5 624,585 V020230 1,860,8 / 3,210,6 / 13,473, fruits (endoca rp) Southern England 8F 88 33 575 159,850 P0681957 1,219,0 Dorstenia picta Dorstenie Lecanium / 6,315,5 Treculia obovoidea Dorstenie P0687715 10,640, ae ae 0 24 ae 11 ae 7 372 SG Morus Morac eae Morus tymensis Dorofeev fruits USSR Broussonetia papyrifera Dorstenie 51,643, Dorstenia contrajerva Dorstenie Dorstenia / 3,208,9 Dorstenia poinsettiifolia Dorstenie Lomathophora P0682278 3,057,4 Trilepisium Dorstenie P0106014 536,683 ae 997 ae 13 ae 9 80 madagascariense ae 3 ae a 71 ae 93 scabra ae 4F 47 ae Dorstenia africana Dorstenie Nothodorsteni / 4,238,8 Dorstenia contrajerva-1 Dorstenie Dorstenia / 5,296,8 Dorstenia psilurus var. Dorstenie Lomathophora V009638 3,595,8 Trymatococcus oligandrus Dorstenie 319,596	30.89% 32.63% 23.38% 18.89% Chandler, 1962 42.95% 38.60% 40.91% 38.00% Chandler, 1961 67.41% 13.59% 2.67% 27.93% Takhtajan, 1982 0.90% 3.16% 3.24% 13.89% 34.22% 42.67% 26.95% 37.93%	102 18/01/1989 101 09/10/1990 97 12/04/2002 33 21/09/1983 Collinson, 1989 102 68 19/04/1977 102 16/08/1951 102 26/06/2013 Chandler, 1961; Collinson, 1989 100 69 16/01/2000 94 05/03/1993 102 04/10/2009 Collinson, 1989 89 102 24/05/1993 74 08/12/1999 39 16/09/2013 and 1911 102 27/05/1993 102 16/10/1997 79 between 1893 101	Oligocene Maciel & Rosario 1552 Roberto Fontes Vieira et al 554(NY) A Fuentes 4353(MO) J.J. Floret, A.M. Louis 1351(P) boundary Paleocene and 56.0 Collinson, 1989 McPherson 20944 (F) A. Krapovickas & A. Schinini 31524(F) Grady L Webster 4073(NY) N. Zerega 909 (SAN) Eocene boundary Eocene and 33.9 Collinson, 1989 C. Galdames 4403 (F) X.M. van der Burgt 576(P) RE Gereau et al. 5189(MO) Cheek, M. 15081(P) Oligocene boundary Eocene and 33.9 GenBank SRA reads: accession no. SRR1477753 Ulises Chavarria 809(F) Simons, E.L.A.N. & R. Westerduijn 319(P) Z. Ranktonirina & Georges Be 197(P) Collinson, 1989 Oligocene boundary Carvalho 5308(NY) E. Martinez S. et al 28913(NY) G.A. Zenker s.n. W. Hahn 3649
	Dorstenia alta Dorstenia djettii Dorstenia psilurus-1 Utsetela gabonensis		Dorstenie Dorstenie Dorstenie Dorstenie	Xylodorstenia Nothodorsteni Lomathophora	MO25064 MO25064 / P0103782	2,711,2 4,776,1 2,712,0 5,988,5	3.60% 18.56% 7.00% 22.98%	87 24/01/2001 102 12/03/1996 86 27/04/1999 102 27/09/1997	A Ntemi Sallu 638(MO) J Amponsah et al 1413(MO) Tim Flynn et al 6522(NY) F.J. Breteler 14096(P)
				ae ae ae ae	a		16 27 9		47 55 56 93

Table S3 .

 S3 List of ancestral distribution area estimated for several nodes by BioGeoBEARS. a) time-constant models (model0); b) time-stratified

	models (model1)
	a)
	model	CG	SG	CG	SG	CG	SG Neo	CG Neo	Node
		Dorstenieae	Brosimum	Brosimum	Dorstenia	Dorstenia	Dorstenia	Dorstenia	I
			s.l.	s.l.
	DEC+Jm0	CDF	C	A	C	C	C	A	C
	DECm0	CDF	AC	A	C	C	AC	A	C
	BayArea+Jm0 F	C	A	C	C	C	A	C
	BayAream0	CF	C	A	C	C	C	A	C
	DIVA+Jm0	F	AC	A	C	C	C	A	C
	DIVAm0	F	AC	A	C	C	AC	A	C
	b)
	model	CG	SG	CG	SG	CG	SG Neo	CG Neo	Node
		Dorstenieae	Brosimum	Brosimum	Dorstenia	Dorstenia	Dorstenia	Dorstenia	I
			s.l.	s.l.
	DEC+Jm1	CDF	AC	A	C	C	A	A	AC
	DECm1	CDF	AC	A	C	C	AC	A	AC
	BayArea+Jm1 F	C	A	C	C	C	A	C
	BayAream1	CF	AC	A	C	C	C	A	C
	DIVA+Jm1	F	AC	A	C	C	A	A	A
	DIVAm1	F	AC	A	C	C	AC	A	AC

  , generalized linear models; GEE, generalized estimating equations; M, use median for the tip value; FD, full dataset; Asia, subset dataset with Asia-Oceania Ficus only; bio18min, dataset assembeld from the corresponding values of each variables from coordinate of the minimum value of variable "bio18" of each species

	var	GLM_M_	GLM_M_FD	GEE_M_	GLM_M_	GLM_M_Asia	GEE_M_	GLM_bio18mi	GLM_bio18min_F	GEE_bio18mi
		FD	_coef	FD	Asia	_coef	Asia	n_FD	D_coef	n_FD
	PC1	***	-0.2242 not sig.	.	-0.13621 not sig.	*	-0.16617 not sig.
	PC2	***	-0.29111 not sig.	*	-0.23273 not sig.	***	-0.33875 not sig.
	PC1+2	***,***	-0.25431, -	not sig.	.,*	-0.14615,-	not sig.	***,***	-0.17919,-0.34749	not sig.
			0.34284			0.24375
	bio1 = Annual Mean Temperature	not sig.	/	not sig.	not sig.	/	not sig.	not sig.	/	not sig.
	bio5 = Max Temperature of	**	-0.022318 not sig.	*	-0.020294 not sig.	**	-0.012796 not sig.
	Warmest Month
	bio6 = Min Temperature of Coldest	not sig.	/	not sig.	not sig.	/	not sig.	not sig.	/	not sig.
	Month
	bio8 = Mean Temperature of Wettest	not sig.	/	not sig.	not sig.	/	not sig.	not sig.	/	not sig.
	Quarter
	bio9 = Mean Temperature of Driest	not sig.	/	not sig.	not sig.	/	not sig.	not sig.	/	not sig.
	Quarter
	bio10 = Mean Temperature of	not sig.	/	not sig.	not sig.	/	not sig.	not sig.	/	not sig.
	Warmest Quarter
	bio11 = Mean Temperature of	not sig.	/	not sig.	not sig.	/	not sig.	not sig.	/	not sig.
	Coldest Quarter
	bio12 = Annual Precipitation	***	0.10051 not sig.	**	0.07439 not sig.	***	0.06578 not sig.
	bio13 = Precipitation of Wettest	***	0.2675 not sig.	*	0.2823 not sig.	.	0.07163 not sig.
	Month
	bio14 = Precipitation of Driest	***	0.23656 not sig.	*	0.11984 not sig.	***	0.23511 not sig.
	Month
	bio16 = Precipitation of Wettest	***	0.17727 not sig.	*	0.15527 not sig.	*	0.05061 not sig.
	Quarter
	bio17 = Precipitation of Driest	***	0.13509 not sig.	*	0.06915 not sig.	***	0.13624 not sig.
	Quarter
	bio18 = Precipitation of Warmest	***	0.28233 not sig.	***	0.21337 not sig.	***	0.19395 not sig.
	Quarter
	bio19 = Precipitation of Coldest	*	0.04934 not sig.	.	0.04571 not sig.	not sig.	/	not sig.
	Quarter
	Significance codes: '***' P<0.0001; '**' P<0.001; '*' P<0.01; '.' P<0.05; 'not sig.' not significant.

GLM(see text).

Université Paris-Saclay

Espace Technologique / Immeuble Discovery Route de l'Orme aux Merisiers RD 128 / 91190 Saint-Aubin, France

ACKNOWLEDGEMENTS

I would like to express my deepest gratitude to my supervisors Dr. Hervé Sauquet and Dr. Renske Onstein, for their patience and help throughout my PhD journey. Their advice

ACKNOWLEDGEMENTS

We thank Jos Käfer for his constructive review of this paper; Sophie Nadot for comments on an earlier draft; Selena Smith for comments on fossils of Artocarpus; Susana Magallón for discussion of divergence time estimation; and Nyree Zerega and Elliot Gardner for discussion of breeding system evolution in Moraceae.

ACKNOWLEDGEMENTS

We thank the Field Museum (F), the Missouri Botanical Garden (MO), the New York Botanical Garden (NY), and the Muséum national d'Histoire naturelle (P) for sharing herbarium material; Tracy Misiewicz and Caroline Storer for assistance with DNA extraction; Charles Foster for providing helpful advice on running MCMCTree; and Santiago Ramirez-Barahona for information on fossil calibrations. This work was supported by a China Scholarship Council (CSC) PhD grant (grant No. 201506140077) and a research grant from the International Association for Plant Taxonomy (IAPT) to Q.Z. Sequencing costs were covered by Agence Nationale de la Recherche grant ANR-12-JVS7-0015-01 to H.S.

ACKNOWLEDGEMENTS

This work was supported by a China Scholarship Council PhD grant [grant No. 201506140077] to Q.Z.

its applications to single-cell sequencing. Journal of computational biology : a journal of computational molecular cell biology 19: 455-477. Berg CC. 1977. Revisions of African Moraceae (excluding Dorstenia, Ficus, Musanga and Myrianthus). Bulletin du Jardin botanique national de Belgique/Bulletin van de Nationale Plantentuin van Belgie 47: 267-407. Berg CC. 2001. Moreae, Artocarpeae, and Dorstenia (Moraceae): with introductions to the family and Ficus and with additions and corrections to Flora Neotropica Monograph 7. Flora Neotropica 83: 1-346. Berg CC, Hijman MEE. 1999. The genus Dorstenia (Moraceae). Ilicifolia 2: 1-211. Bouckaert R, Heled J, Kühnert D, et al. 2014. BEAST 2: A Software Platform for Bayesian Evolutionary Analysis. PLOS Computational Biology 10: e1003537. Chandler MEJ. 1961. Flora of the Lower Headon beds of Hampshire and the Isle of Wight. British Museum (Natural History). Chandler MEJ. 1962. The Lower Tertiary Floras of Southern England: Flora of the pipe-clay series of Dorset (lower Bagshot). order of the Trustees of the British Museum. Chase MW, Thijs KW, Kamau P, Fay MF. 2013. Dorstenia christenhuszii (Moraceae), a new species from the Taita Hills, Kenya. Phytotaxa 81: 45-48. Christenhusz MJM, Chase MW. 2013. Biogeographical patterns of plants in the Neotropicsdispersal rather than plate tectonics is most explanatory. Botanical Journal of the Linnean Society 171: 277-286. Chung K-F, Kuo W-H, Hsu Y-H, Li Y-H, Rubite RR, Xu W-B. 2017. Molecular recircumscription of Broussonetia (Moraceae) and the identity and taxonomic status of B. kaempferi var. australis. Botanical studies 58: 11. Clement WL, Weiblen GD. 2009. Morphological evolution in the mulberry family (Moraceae). Systematic Botany 34: 530-552. Collinson ME. 1989. The fossil history of the Moraceae, Urticaceae (including Cecropiaceae), and Cannabaceae. Crane, P, R,, Blackmore, S. (Ed s.), Evolution, systematics, and fossil history of the Hamamelidae 2: 319-339. Couvreur TLP, Helmstetter AJ, Koenen EJM, et al. 2019. Phylogenomics of the Major Tropical Plant Family Annonaceae Using Targeted Enrichment of Nuclear Genes. Frontiers in Plant Science 9: 1941. Couvreur TLP, Pirie MD, Chatrou LW, et al. 2011. Early evolutionary history of the flowering plant family Annonaceae: steady diversification and boreotropical geodispersal. Journal of Biogeography 38: 664-680. Fisher AE, Hasenstab KM, Bell HL, Blaine E, Ingram AL, Columbus JT. 2016. Evolutionary history of chloridoid grasses estimated from 122 nuclear loci. Molecular phylogenetics and evolution 105: 1-14.

as ancestral in Ficus, whereas the unequal-rate model supported gynodioecy (or gynodioecy + dioecy or gender dimorphic) as ancestral in the genus (Table 3). The strength of support differed among different configurations but the ancestral state did not change.

In the ML approach, we compared the models in each configuration using the Akaike Information Criterion (AIC). Only with configuration A was the equal-rate model selected, whereas the unequal rate model better fit the data with configurations B-E. As a result, considering best-fit models only, gynodioecy (or gynodioecy + dioecy, or gender dimorphic) was found as ancestral in Ficus with five out of six configurations (Table 3). Third, the Bayesian approach with reversible-jump Markov chain Monte Carlo (rjMCMC) weakly supported monoecy (or gender monomorphic) as ancestral in Ficus, but with a broad 95% highest posterior density (HPD; 0.00 to 1.00), suggesting high uncertainty in the estimation (Table 4). Bayesian analyses with fixed model selected the unequal rate model in 4 out of 6 configurations, two with weak support), and monoecy as ancestral in Ficus in 4 out of 6 configurations.

SUPPLEMENTARY INFORMATION

Figure S3

. Ancestral state reconstruction with 320-species dataset by parsimony approach with tip names with all six configurations (see text).

Figure S4.

Ancestral state reconstruction with 320-species dataset by equal-rate maximum likelihood with configurations A to F (configuration B presented as Fig. 3).

Figure S5

. Ancestral state reconstruction with only Ficus species in the dataset with tip names by different approaches.

Table S1. List of GenBank accession numbers for the sequences used in this study.

Table S2. detailed information for the fossil calibrations used in this study.

Table S3. Extraction of the PROTEUS database including a list of all data records and linked references, lists of characters and character states used in this study, and the matrix in ancestral state reconstruction analyses.

CHAPTER II. LONG-DISTANCE DISPERSAL SHAPED THE DIVERSITY OF TRIBE DORSTENIEAE (MORACEAE)

This study was submitted to the preprint server bioRxiv in January 2019, in collaboration with Elliot Gardner, Nyree Zerega and Hervé Sauquet. Both penalized likelihood (PL) and Bayesian relaxed clock approaches were used to estimate divergence times in Dorstenieae. PL was implemented in r8s v1.7 [START_REF] Sanderson | r8s: inferring absolute rates of molecular evolution and divergence times in the absence of a molecular clock[END_REF] with strict minimum and maximum age constraints as described above. For the PL approach, the best smoothing value was first determined using cross validation by testing 21 values of smoothing parameter scaling from 0.1 to 1000. The optimum value (i.e., with lowest chi-square) of 1.6 was then used as the smoothing parameter in divergence time estimation.

ABSTRACT

For the Bayesian approach, we used MCMCTree as implemented in the PAML v4.9 package (Yang 2007). MCMCTree has been used to estimate divergence times with other similar phylogenomic datasets (e.g., Foster et al., 2017), where other programs such as BEAST2 [START_REF] Bouckaert | BEAST 2: A Software Platform for Bayesian Evolutionary Analysis[END_REF] would take too long to converge. The tree topology was fixed, using the best-scoring tree from the RAxML analysis, and all calibrations were implemented as uniform priors with soft boundaries (2.5% on both sides) (Ho and Phillips 2009). We set the minimum age of the fossil and as the minimum boundary and the younger bound of the root (73 Ma) as the maximum boundary of the uniform distribution for each fossil calibration. The birth-death process was used as the tree prior. For the MCMCTree analysis, we used the supermatrix, analyzed as a single partition under the GTR substitution model.

We first ran baseml in PAML v4.9 with a strict clock model to estimate the rough mean of parameters such as the shape parameter for the overall rate and the transition/transversion rate ratio. brasiliensis-1 and D. arifolia were not (but still lie in the Neotropical clade). Although in similar positions on the phylogenetic trees reconstructed by two approaches, the support for the splits of these two species were low in the coalescent tree (bootstrap support less than 50%) but high in the maximum likelihood tree (bootstrap support over 90%). Several nodes were comparatively weakly supported in Dorstenia (Fig. 2). Our samples represented eight out of nine recognized sections in Dorstenia. None of them was found to be monophyletic in our analyses (Figure S1). All of the Neotropical species of Dorstenieae outside of Dorstenia also formed a clade (i.e., Brosimum s.l., incl. Trymatococcus and Helianthostylis).

Asia and Oceania (CDF, Fig. 3 (Jaramillo et al. 2006).

During this period, global temperature increased, leading to the middle Eocene climatic optimum [START_REF] Zachos | An early Cenozoic perspective on greenhouse warming and carbon-cycle dynamics[END_REF]). The increasing temperature has been suggested as one of the factors for the extension of the Neotropical region in the Eocene (Hughes et al. 2012). A larger Neotropical area at that time would have increased the probability of the successful colonization of species from Africa to a suitable habitat in the Neotropics.

Curiously, in these reconstructions, D. elliptica (an African species sister to the Neotropical Dorstenia clade) was inferred to be the result of a back dispersal from South America to continental Africa (Fig. 3,S3c). Whether this somewhat unexpected clustered with U. gabonensis (D. alta, and D. mannii and D. turbinata) and reran the phylogenetic analyses with both methods. Utsetela gabonensis was still in the same position (nested in Dorstenia) after excluding these three species (results not shown).

Sampling the other species, U. neglecta [START_REF] Jongkind | Novitates gabonenses (24) A New Species in Utsetela (Moraceae)[END_REF], would be necessary to confirm this relationship and draw any taxonomic conclusions. Some differences among this basal grade of taxa include: tepals of individual flowers in the inflorescences of Dorstenia are connate, while they are free in Scyphosyce and Utsetela (Berg 1977).

While all Scyphosyce, Utsetela and Dorstenia have drupe(let) of fruit, their receptacles are different in shape and the filaments are far more elongated in Utsetela then in other two genera (Berg 1977;Berg and Hijman 1999). Despite some differences, similarities in woody habit, fruit type, macrospermy, and the phylogenetic reconstruction presented here suggest that merging Scyphosyce and Utsetela into Dorstenia may be a reasonable taxonomic outcome to preserve the monophyly of Dorstenia. An alternative option would be to separate those clades into two separate genera as elaborated on below.

Nine sections were proposed by Berg and Hijman (1999) in Dorstenia based on inflorescence, habit and life form (i.e. geophytes, phanerophytes, hemicryptophytes) characters. Our sampled species represented all the sections except Bazzemia, which consists of a single species in Mozambique. None of the eight sampled sections were found to be monophyletic in this study (Figure S1) and the phylogeny can help inform future subgeneric classification of Dorstenia. Of particular interest to consider taxonomically are the two most basal Dorstenia clades that include members of the sections Xylodorstenia and Nothodorstenia as well as species from two other genera (Scyphosyce and Utsetela). Misiewicz and Zerega (2012) did not include any species of Scyphosyce and Utsetela, and they found sections Xylodorstenia monophyletic. This is not the case in our reconstruction (Figure S1), and is likely due to our increased Table S1. List of specimens collected in this study.

SUPPLEMENTARY INFORMATION

Table S2. AIC of biogeographic reconstruction with time-constant (model0, a) or timestratified (model1, b) DEC, DIVA-like and BayArea-like based models.

Table S3. List of ancestral distribution area estimated for several nodes by BioGeoBEARS. a) time-constant models (model0); b) time-stratified models (model1) .

CHAPTER III. BREEDING SYSTEM EVOLUTION AND CLIMATE IN FICUS

This study has not been submitted to any journal yet and was conducted in collaboration with Renske Onstein and Hervé Sauquet.

ABSTRACT

 Background and Aims Physiological and ecological differences of breeding systems may lead to different distribution patterns. Two breeding systems, gynodioecy and monoecy, occur in Ficus (ca. 750 spp., Moraceae), an ecologically important genus distributed in all major tropical regions. We used Ficus as a model to investigate the niche difference of different breeding systems.

 Methods Richness maps of 183 species of Ficus (105 monoecy, 78 gynodioecy) were plotted.

To investigate the niche difference of the two breeding systems in Ficus, regression was conducted with fourteen climatic variables against breeding systems with generalized linear models (GLM) and generalized estimating equations (GEE). Ancestral states of precipitation niche were reconstructed to test its relationship with the transitions of breeding systems in the history of Ficus.

 Key Results Although gynodioecious and monoecious Ficus overlap in their distribution, gynodioecious Ficus does not occur in the Neotropics and monoecious Ficus is distributed in drier environments. The distribution patterns of both breeding systems fit the summer positions of the intertropical convergence zone on continents. A significantly positive relationship of precipitation and gynodioecy was supported by GLM but not GEE analyses.

 Conclusions Niche differences were found in Ficus with different breeding systems but may be the result of phylogenetic effect rather than historical correlation. This study sheds light on understanding of the origin, diversification, and distribution pattern of the two breeding systems in Ficus.

Key words: breeding systems, Ficus, monoecy, gynodioecy, generalized estimating equations (GEE), generalized linear models (GLM), climatic variables, intertropical convergence zone (ITCZ)

CoordinateCleaner [START_REF] Antonelli | Amazonia is the primary source of Neotropical biodiversity[END_REF]. In addition, we chose records later than the year 1945; of which precision was higher than 10 km; of which sources were "human observation", "observation"

or "preserved specimen". Given both the presence of invasive and cultivated records in GBIF, we finally skimmed the records manually to exclude records outside the reported natural distribution area of each species by searching local floras and additional literature (Table S1). For instance, records in New Zealand were excluded as Ficus is naturally absent there [START_REF] Gardner | The naturalisation of banyan figs (Ficus spp., Moraceae) and their pollinating wasps (Hymenoptera: Agaonidae) in New Zealand[END_REF].

Considering the long history and wide range of the cultivated plant Ficus carica, it is hard to determine the original native distribution of this species (Ghada et al. 2010), which was therefore excluded from our analyses. Sixteen additional species were excluded due to either the lack of GBIF data (Ficus arfakensis, Ficus bullenei, Ficus dugandii, and Ficus palmeri) or synonymy (e.g., two species included as distinct tips in our dataset, but treated as synonyms in GBIF). As a result, our final dataset consists of 183 species of Ficus with GBIF data, including 23 Neotropical species (all monoecious), 42 African species (38 monoecious, 4 dioecious), and 118 species in Asia-Oceania (44 monoecious, 74 dioecious). After all of the above steps of data cleaning, we plotted separate richness maps for each breeding system in a spatial resolution of one minute using package SpeciesGeoCodeR [START_REF] Töpel | SpeciesGeoCoder: fast categorization of species occurrences for analyses of biodiversity, biogeography, ecology, and evolution[END_REF]. Species assignment to breeding systems follows Zhang et al. (2019). The complete dataset of cleaned occurrence records is available from supplementary data 1.

Principal component analysis and regression

Climatic data were retrieved from the WorldClim database [START_REF] Fick | WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas[END_REF] at the spatial resolution of 2.5 minutes for the cleaned occurrence records obtained in the above steps.

Because some of these variables were calculated from others (e.g. mean diurnal range, BIO2 in

Ancestral state reconstruction of climatic variables

Results from the regression analyses suggested a possible relationship of precipitation with breeding system. To further investigate and visualize which transitions in breeding systems were associated with change of precipitation in the history of Ficus, we reconstructed the ancestral value of precipitation of the warmest quarter (bio18) with a Brownian motion model using the maximum likelihood approach implemented in the fastAnc function of R package phytools [START_REF] Revell | phytools: an R package for phylogenetic comparative biology (and other things)[END_REF].

Here, we aimed to test the influence of precipitation on the distribution of Ficus. The environmental data in our dataset were retrieved from the coordinates of the occurrence records, suggesting they are all in the tolerant ranges of Ficus. Both temperature and precipitation are important to the development to both the Ficus and their pollinating wasps. We selected bio18 to control one factor and observe the performance of the other. When analyzing precipitation of the warmest quarter (bio18), temperature is the warmest in the period (year), hence we control the temperature as the best and are testing the relationship of precipitation with breeding system.

Additional analyses

Two additional analyses were conducted to further explore and characterize potential signal for niche difference between breeding systems in the dataset. Firstly, considering the fact that dioecious Ficus extended more poleward than monoecious Ficus only in Asia, we conducted all the analyses above with a subset of the dataset including only Asia-Oceania Ficus (118 species).

These analyses were prompted by the observation that dioecious species of Ficus tend to extend to higher latitudes than monoecious species only in this region. Hereafter, we refer to this dataset was named as subsetAsia. Secondly, we built a separate species-climate dataset by extracting the

Ancestral state reconstruction for precipitation

The ancestral state reconstruction of precipitation of the warmest quarter (bio18) reveals that nearly all of the drier niches in Ficus are occupied by monoecious species (Fig. 3), providing a visual confirmation of the relationship between breeding systems and precipitation outlined by the GLM analyses. The estimated mean value of root of the tree was 612 mm (95% confidence interval:

308-916 mm). Sampled species in the monoecious section Galoglychia, which is distributed from Africa to Arabian Peninsula including Madagascar and surrounding islands [START_REF] Berg | Flora Malesiana precursor for the treatment of Moraceae 7: Ficus subgenus Urostigma[END_REF]), showed comparatively lower values of bio18.

DISCUSSION

Different results in GLM and GEE

All the climatic variables related with precipitation (bio 12-19) were significantly correlated with dioecy in Ficus in regressions with generalized linear models (GLM), but not those with the generalized estimating equations (GEE) approach [START_REF] Paradis | Analysis of comparative data using generalized estimating equations[END_REF], which takes phylogenetic relationships into account. This difference could be a type I error for GLM to find significance with higher numbers of degrees of freedom compared with the counterpart which adapted with phylogenetic relationship in GEE [START_REF] Paradis | Analysis of comparative data using generalized estimating equations[END_REF].

SUPPLEMENTARY INFORMATION

Table S1. Species list and distribution information of the dataset in this study.

Table S2. Estimates for interaction effects of climatic variables and breeding systems in Ficus from the regressions by generalized linear models (GLM) and generalized estimating equations (GEE) for three datasets: 1) full dataset (same as Table 1); 2) subsetAsia; and 3) bio18min.

Supplementary data 1. The cleaned occurrence records used in the analysis of this study.

SUPPLEMENTARY DATA

All the supplementary data attached to this thesis are available from this online folder: https://www.dropbox.com/sh/r09q6qnpxlezpji/AACth8ARb0xr670UrMyEEtKJa?dl=0. As some of the supplementary information are too long to print in a thesis, I have selected some proper supplementary figures and tables to print here.

Summary of my contribution to the

Ficus hispidioides

Rhamnus utilis