Exploring intraspecific variation, habitat differentiation, and evolutionary history in Arabidopsis thaliana: a biogeographical approach

Do˜nana Biological Station University of Sevilla EXPLORING INTRASPECIFIC VARIATION, HABITAT DIFFERENTIATION AND EVOLUTIONARY HISTORY OF ARABIDOPSIS THALIANA: A BIOGEOGRAPHICAL APPROACH PhD Thesis Bernardo Toledo Gonz´alez Sevilla 2023 Do˜nana Biological Station University of Seville EXPLORING INTRASPECIFIC VARIATION, HABITAT DIFFERENTIATION AND EVOLUTIONARY HISTORY OF ARABIDOPSIS THALIANA: A BIOGEOGRAPHICAL APPROACH Dissertation submitted by Bernardo Toledo Gonz´alez to obtain the PhD degree by the Universidad de Sevilla Sevilla, 2023 Doctoral Programme on Integrated Biology Supervisors: Dr. F. Xavier Pic´o Departamento de Biolog´ıa Evolutiva, Estaci´on Biol´ogica de Do˜nana, Consejo Superior de Investigaciones Cient´ıficas, Sevilla Dr. Arnald Marcer Ecological and Forestry Applications Research Centre (CREAF) Universitat Autonoma de Barcelona, Bellaterra Tutor: Dr. Montserrat Arista Departamento de Biolog´ıa Vegetal y Ecolog´ıa Facultad de Biolog´ıa Universidad de Sevilla Financiaci´on Los trabajos que aborda la presente tesis doctoral fueron financiados a trav´es de la Direcci´on General de Investigaci´on Cient´ıfica y T´ecnica del Ministerio Espa˜nol de Ciencia e Innovaci´on (MICINN). El autor, Bernardo Toledo Gonz´alez, recibi´o dicha financiaci´on de las becas para contratos predoctorales del fondo europeo de resiliencia (FEDER, CGL2016-77720-P). Los directores fueron financiados tambi´en por los fondos FEDER en el marco de proyecto nacional (CGL2016-77720-P) y la Agencia de Gesti´on de Ayudas Universitarias y de Investigaci´on de la Generalitat de Catalunya (2017-SGR-1006). Recommended citation Toledo, B. (2023). Exploring intraspecific variation, habitat differentiation and evolutionary history of Arabidopsis thaliana: a biogeographical approach. [Doctoral thesis, University of Seville]. Cover design and layout by: Ernesto Uriszar de Aldaca Document design and layout by: Xandra Campo Blanco Contents Agradecimientos i Abstract iii Resumen vii Publications xi 1. Publications................................... xi 2. Conferences ................................... xi 1. Introduction 1 1.1. Arabidopsis thaliana: a model system for ecology, genetics and evolution inplants..................................... 1 1.2. Intraspecific variation through transplant experiments and genetics analysis 5 1.3. Intraspecific variation: the key to better understand species’ distribution patterns ..................................... 7 1.4. Scopeofthethesis ............................... 9 2. Chapter I: An ecological history of the relict genetic lineage of A. thaliana 13 2.1. Introduction................................... 13 2.2. Materials and methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.2.1. Relict source populations . . . . . . . . . . . . . . . . . . . . . . . . 16 2.2.2. Species distribution models . . . . . . . . . . . . . . . . . . . . . . 18 2.2.3. Fossilpollendata............................ 21 2.2.4. Distribution maps of dominant Iberian tree species over time . . . . 23 2.2.5. Genetics of flowering time in relict A. thaliana ............ 24 2.3. Results...................................... 25 2.3.1. Potential distribution range of relict A. thaliana ........... 25 2.3.2. Vegetation changes from fossil pollen data . . . . . . . . . . . . . . 26 2.3.3. Niche overlap with tree species . . . . . . . . . . . . . . . . . . . . 28 List of Figures 3.5. Normalized effects of geography, genomics and environmental suitability on phenotypic differentiation in A. thaliana per region and trait group given by Maximum-Likelihood Population-Effects (MLPE) mixed models. See Table 3.3 for trait description. . . . . . . . . . . . . . . . . . . . . . . . . 53 4.1. Occurrences of A. thaliana for each genetic cluster (Central Asia, N = 53; Central Europe, N = 95; Germany, N = 47; Iberian Peninsula, N = 102; Italy-Balkan-Caucasus, N = 65; North Sweden, N = 27; Relicts, N = 63; South Sweden, N = 73; Western Europe, N = 77, Total N = 711). Each occurrence represents a unique 30 arcsecond grid cell. . . . . . . . . . . . 64 4.2. Variable importance and response curves of the predictor variables in the environmental suitability models for each genetic group. A) Relative contribution of each environmental variable expressed as the fraction of total variation accounted for (FTVA). B) Response curves of habitat suitability for each predictor variable expressed in normalized probability ratio output (PRO)values. ................................. 69 4.3. Habitat suitability maps for A. thaliana expressed in probability ratio output (PRO) values. Red lines represent the buffer area for the genetic clusters occurrences. Black lines represent the buffer area for the complete occurrence dataset. Shadowed areas are outside of the study area. . . . . 70 4.4. Representation of the environmental space as derived from the PCA analysis. A) Percent contributions of each variable to each of the first four principal components. The red line shows the 10% threshold we used to consider that a variable is important for a principal component. B) Kernel density estimations of the point distribution of each genetic cluster for the first two PCA components (arrows represent the magnitude and direction of the correlation of each environmental variable for each component). . . 72 4.5. Representation of the genetic distance between genetic clusters according to the first three dimensions of the MDS analysis. Shaded surfaces and isolines represent the kernel density estimations for the point distribution of each genetic cluster in genetic space. A) Kernel density estimations for the point coordinates for dimensions 1 and 3, B) dimensions 2 and 3, and C)dimensions1and2. ............................ 73 4.6. Overlap analyses between genetic clusters. We represent the Warren’s I (green) and the Schoener’s D (purple) statistics for the genetic overlap index (index), environmental (PCA) and suitability overlaps. . . . . . . . 74 List of Figures 4.7. MLPE model between genetic overlap index and suitability overlaps. Solid black line represents the MPLE regression line. Dashed line represents 95% Confidenceintervals. .............................. 75 A.1. Current potential distribution ranges of 17 tree species in the Iberian Peninsula. Data from (Benito-Garz´on et al.,2007). ................123 List of Tables 2.1. Variable percent permutation importance of the distribution model fit for the relict lineage of A. thaliana in southwest Mediterranean Basin (Figure 2.3). Mean and standard deviation (SD) values for each variable included in the model are also given for all populations (N = 63) as well as for Iberian (N = 39) and Maghreb (N = 24) populations. Variables with significant differences between the Iberian Peninsula and the Maghreb (P-value < 0.05; one-way ANOVA) are given in bold face. Dimensionless variables havenounits. ................................. 27 2.2. Niche overlap (Schoener’s D) between the potential distribution range of relict A. thaliana and each tree species in the present. Niche overlap values between the potential distribution range of each trees species between the present and the Mid-Holocene (MH), between the present and the Last Glacial Maximum (LGM), and between MH and LGM are also given. All D values were estimated for the Iberian Peninsula. Tree species are ranked in decreasing D values for relict A. thaliana. ................. 29 3.1. Sources for A. thaliana occurrences....................... 40 3.2. List of environmental predictors used in the environmental suitability model for A. thaliana. Columns: Variable (variable name), Description (description of the variable), MB (Model building, predictors selected after assessing collinearity), FM (Final model, predictors selected by MIAMaxent after performing a forward selection on the ’Model building’ set) . . . . . . . . 41 3.3. Environmental suitability model performance. AUC values for each fold in the cross-validation of the selected alpha were the following: fold 1 = 0.747, fold 2 = 0.807, fold 3 = 0.734, fold 4 = 0.739, fold 5 = 0.741. An alpha value of 0.001 (in bold) corresponded to the model with the highest AUC. ...................................... 44 List of Tables 3.4. Number of A. thaliana accessions included in the analyses per trait and region. Sources: [1] The 1001 Genome Consortium (2016); [2] Mart´ınez- Berdeja et al. (2020); [3] Davila et al. (2017); [4] Vasseur et al. (2018a); [5] Exposito-Alonso et al. (2019); [6] Katz et al. (2021). For variables of source [5], MH, ML, TH and TL letters should be read as: M - Madrid experiment site, T - T¨ubingen experiment site, H - high precipitation treatment, and L - low precipitation treatment. . . . . . . . . . . . . . . . . . . . . . . . 46 3.5. Summary of environmental suitability values per region. Summary values are in PRO units calculated from the areas enclosed by a buffer of about 50 km around the convex hull defined by the accessions for each region (Figure3.1b). ................................. 51 3.6. Correlations between geography, genomics and environmental suitability distance matrices in each region according to Mantel tests. . . . . . . . . 52 3.7. MLPE summary of model coefficients per phenotypic category and region. Columns mean, median sum and N show the mean value, the median value, the sum of all values, and the number of values for all coefficients different from zero for each region. The most important predictors per phenotypic category and region are given in pale blue. . . . . . . . . . . . . . . . . . 54 4.1. Summary table for the habitat suitability models of the genetic clusters. Columns: GROUP - Genetic cluster; NO. - number of occurrences; MODEL FORMULA - Best model; ALPHA - Regularisation parameter; and AUCc - Adjusted Area Under the Curve. . . . . . . . . . . . . . . . . . . . . . . 68 4.2. FTVA values for each variable across all genetic clusters in decreasing order. 68 4.3. Summary table for the MPLE models of the relation between genetic, environmental and suitability overlaps between genetic clusters. Models = Models’ formula; Int = Intercept of the model; PCA = PCA variable estimate; Suitability = Suitability variable estimate R2 N= R square of Nagelkerke; df = Degrees of freedom; LogLik = Log likelihood; AIC = Akaike’s Information Criterion; ∆AIC = Delta AIC; Wt = Model weight. Selected modelismarkedinbold. ........................... 75 A.1. Region, accession code, altitude (m.a.s.l.) and coordinates (degrees) of relict A. thaliana populations known to date (Figure 2.1). Unknown altitudes are indicated by n/a. Accessions sharing the same geographic coordinates are given once and the others indicated by dashes. . . . . . . . . . 117 List of Tables A.2. Location, country, maximum calibrated years before present (MCYBP), altitude (m.a.s.l.), coordinates (degrees) and source of the 80 paleobotanical sites from which fossil pollen data were available (Figure 2.2). Acronyms in column Source should be read as: EPD - European Pollen Database, RD - Rachid Cheddadi, PGS - Pen´elope Gonz´alez-Samp´eriz, YM - Yannick Miras, JALS - Jos´e Antonio L´opez-S´aez, GGR - Graciela Gil-Romera, WF - Will Fletcher, FB - Francesc Burjachs, LSF - Luisa Santos Fidalgo, SJ - SamanthaJones. ................................120 Agradecimientos Es dif´ıcil mencionar y hacer justicia en pocas l´ıneas a todas aquellas personas que con su apoyo, tanto laboral como personal, han hecho posible la finalizaci´on de esta tesis, y de este camino una aventura m´as agradable. Aun as´ı, voy a intentar no dejarme a nadie atr´as. A F. Xavier Pic´o, director de esta tesis, gracias por haberme dado la oportunidad de trabajar estos 4 a˜nos en la Estaci´on Biol´ogica de Do˜nana. Primero, agradecerte tu paciencia en los errores que he cometido durante este proceso de aprendizaje, sobre todo en mis continuos despistes administrativos y en los momentos m´as dif´ıciles de toma de decisiones, cuando has sido siempre paciente y resolutivo. Gracias por la confianza que has depositado en mi desde el principio y en los momentos m´as bajos donde no ve´ıa el camino. Siempre has estado ah´ı para contestar cualquier pregunta o resolver cualquier problema. Gracias por tu capacidad de escucha y tu entusiasmo por la ciencia, por esas charlas sobre la carrera cient´ıfica que me han ayudado a buscar mi propio camino, y por la confianza que me has dado para expresarte mis problemas o dudas. Tu labor como director ha ido m´as all´a de una labor exclusivamente laboral, siempre con optimismo y ganas. A Arnald Marcer, tambi´en director de esta tesis, por tus excelentes recomendaciones y tu visi´on t´ecnica y pragm´atica que me han llevado a no meterme en laberintos de an´alisis infinitos. Gracias por acogerme en Barcelona y ayudarme con mis primeros pasos anal´ıticos. A d´ıa de hoy sigo manteniendo tu orden y maneras de hacer en mis proyectos, lo que me ha ayudado a salir de mi tendencia al caos. Gracias por tu dedicaci´on y prontas correcciones, siempre al pie del ca˜n´on para que esto saliera adelante. Tus mensajes de apoyo y ´animo y tu actitud positiva me han sacado de mis momentos m´as bajos. Gracias por tu comprensi´on y tu paciencia en mis momentos m´as ca´oticos y por ayudarme a encontrar las respuestas que necesitaba. i Agradecimientos Tambi´en agradecer a mis compa˜neros del despacho 5.6 por acogerme como uno m´as del grupo. A Montse Vil`a por incluirme en sus reuniones y dejarme participar y nutrirme de sus discusiones cient´ıficas y sus cenas de grupo. A Javi Galan y Alejandro Trillo por aguantar mis d´ıas buenos y malos, y por hacer de un trabajo solitario una experiencia compartida. A Antonio Castilla, aunque compartimos brevemente despacho y grupo, gracias por tu energ´ıa, positivismo y experiencia, tu amor por la ciencia han acompa˜nado toda mi tesis desde el principio como ejemplo de perseverancia y buen hacer. Como no, agradecer a mis compa˜neros de la Estaci´on Biol´ogica de Do˜nana, muchos de ellos muy buenos amigos por todos estos a˜nos. A Pedro Romero, compa˜nero infatigable de viajes y aventuras. Sevilla no hubiera sido lo mismo sin ti, aunque a lo mejor mi tesis se hubiera depositado un poquito antes. Gracias por ense˜narme tu amor por el campo, por incluirme en el mundo de las expediciones y por las oportunidades de colaborar con tu grupo, que me han ayudado a mejorar como cient´ıfico y persona. A Pepe Tella y a Martina Carrete por acogerme en su grupo como uno m´as de la familia, donde he podido crecer como cient´ıfico y disfrutar de su maravillosa visi´on de la ciencia. Tambi´en tengo que agradecer a Dani Gracia su apoyo incondicional, por compartir tantos d´ıas de trabajo y convertirlos en un entretenimiento constante. A Mar por recibirme con un t´e cada ma˜nana y dejarme unos minutos de desahogo y risas antes de ponernos serios con el trabajo. A mis compa˜neros de piso adoptivos Eneko Redondo y Manu D´ıaz, gracias por acogerme en vuestra casa como si fuera la vuestra y aguantar mis locuras. A Vane y a Migue por las cenas y salidas juntos, para mi sois parte imprescindible de Sevilla. Gracias a los predoctorales en general por ser una familia con su apoyo, sus sue˜nos y sus locuras. Estoy especialmente agradecido a la persona que me introdujo a la Estaci´on Biol´ogica de Do˜nana y a Sevilla. Mi compa˜nero de vida durante los primeros a˜nos de tesis, Javi Manzano. Tras estudiar juntos varios a˜nos, tu recomendaci´on me llevo a esta grand´ısima ciudad, hemos compartido momentos inolvidables en Rayito de Luna, y los que nos quedan. A mis amigos de aventuras Noruegas Ana Sanz y Pablo De la Pe˜na sin los que no habr´ıa encontrado mi camino hasta el doctorado. Por ´ultimo, a mi familia, por su apoyo incondicional para perseguir mis metas y a mis padres en especial por el esfuerzo econ´omico para poder seguir esta dif´ıcil carrera. A Xandra, mi pareja, que me ha guiado y acompa˜nado estos ´ultimos tres a˜nos. Gracias por hacer esta tesis tambi´en tuya y estar vigilando que no me perdiera por el camino. No ha sido nada f´acil, con una pandemia de por medio, tres cambios de ciudades y muchas aventuras m´as. Gracias por estar ah´ı. ii Abstract Various evolutionary forces play a crucial role in shaping the genetic patterns observed between populations and contribute to the intraspecific variation of a species. These forces encompass a range of factors, including random changes, such as mutation and genetic drift, spatial dynamics determining migration and gene flow, and selection pressures driven by the environment. Over time, these forces lead to changes in allele frequencies and ultimately result in population divergences, some of which can facilitate adaptive divergences. As these evolutionary changes accumulate, they contribute to the genetic structure of populations and plasticity, giving rise to distinct lineages and, potentially, the emergence of new species. Therefore, the understanding of intraspecific variation, genomic and phenotypic variation, is crucial to disentangle the evolutionary history of the lineages and to reconstruct the natural events, such as climate fluctuations and anthropogenic factors, which alter the species demography and the balance between evolutionary forces, generating the patterns we observe. This doctoral thesis focuses on the population differentiation, intraspecific variation, and evolutionary history of Arabidopsis thaliana, a model plant for ecology, genetics, and evolution. Arabidopsis thaliana, a small herbaceous plant belonging to the mustard family, has been extensively studied due to its suitability as a model organism. Therefore, it is one of the best sampled species in the world and with more genetic information of wild variants publicly available. Here, I take a biogeographical approach to investigate the forces that shaped patterns of intraspecific variation in A. thaliana across its native range. This study combines the analysis of phenotypic characteristics observed in experimental studies, the genetic structure of the species determined through whole-genome sequences, and the environmental characteristics of the locations where the species’ accessions are found. The main focus of the thesis is to examine intraspecific variation at the levels of the relict lineage, the overall set of genetic lineages or clusters, and the species’ distribution across the Eurasian region. iii Publications 1. Publications Toledo, B., Marcer, A., M´endez-Vigo, B., Alonso-Blanco, C., Pic´o, F. X. (2020). An ecological history of the relict genetic lineage of Arabidopsis thaliana. Environmental and Experimental Botany, 170, 103800 https://doi.org/10.1016/j.envexpbot.2019. 103800 2. Conferences Toledo, B.; Marcer, A.; Pic´o, F. X. (2019, February). Biogegraphy of Arabidopsis thaliana across its native (relict and contemporary) and introduced ranges [Poster presentation]. First Meeting of the Iberian Ecological Society (SIBECOL), Barcelona, Spain. Toledo, B.; Marcer, A.; Pic´o, F. X. (2019 March). Biogeography of Arabidopsis thaliana relict linage across the Iberian Peninsula and Morocco [Poster presentation]. XVI Annual Meeting of the Spanish Group for Floral Ecology (EcoFlor), Granada, Spain Toledo, B.; Marcer, A.; Pic´o, F. X. (2021, November). Disentangling the current distribution of the annual plant Arabidopsis thaliana in its Native range [Poster presentation]. II Simposio Anual de Bot´anica Espa˜nola (SEBOT), Madrid, Spain xi 1. Introduction 1.1. Arabidopsis thaliana: a model system for ecology, genetics and evolution in plants In this thesis I focused on population differentiation, intraspecific variation and the drivers that govern the evolutionary history of A. thaliana in its native range. This small herbaceous plant was first published by Linnaeus in 1753 in the Species Plantarum II, previously described in 1577 by Johannes Thal and finally named as A. thaliana (L.) Heynh by Gustav Heynhold in 1842. A. thaliana is an annual plant that belongs to the mustard family (Brassicaceae). The mustard family is a wide and morphologically diverse plant group including economically important crops such as oilseed weeds, vegetables and spice plants. Although the mustard family has different important plant groups in the economic sphere, A. thaliana remains as the most studied genus in this family. The Arabidopsis genus has undergone multiple taxonomic changes, being at first a group with a variety of plants with convergent morphological traits. The latest taxonomic studies, incorporating an extensive molecular phylogeny, resulted in the reorganization of only nine species in the Arabidopsis genus (O’kane and Al-Shehbaz,1997). These phylogenetic and taxonomic studies also involved the exclusion of several other species that were historically assigned to the genus due to their close convergence, despite being phylogenetically unrelated (O’kane and Al-Shehbaz,2003). Among these nine species, six are considered as the closest wild relatives of A. thaliana:A. lyrata,A. croatica, A. cebennensis,A. arenosa and A. halleri.A. thaliana emerged as a model system for genetic and evolutionary studies more than 50 years ago (Koornneef and Meinke,2010). Recently, the research community has also adopted its closest relatives to understand the evolutionary history of all these species (Honjo and Kudoh,2019;Koch et al.,2008; Mitchell-Olds,2001). 1 1. Introduction Figure 1.1.: Life cycle of A. thaliana. (A) A. thaliana of the accession Columbia (Col) at different stages of its life cycle, from seed (bottom left) to seedling (11 days), to vegetative growth (39 days), and to reproductive growth (45 days). Photographs of (B) a flower, (C) a pollen grain (scanning electron micrograph), and (D) mature siliques (seed pods; left: closed; right: open with a few remaining unshattered seeds) at higher magnification. Image credits: B and C, Maria Bernal and Peter Huijser; other photographs, Ines Kubigsteltig and Klaus Hagemann. Figure and caption extracted from Kr¨amer (2015). 2 1.1. Arabidopsis thaliana: a model system for ecology, genetics and evolution in plants Figure 1.2.: Breakdown of 54 033 A. thaliana publications by research area, from 1965 to 2015 (partial data for 2015). BIOSIS concept codes were used to flag papers in the 54 033 set. If a paper had been assigned multiple codes, the least abundant code was used for this display, out of 46 that were chosen to provide a balance between not being too general (75% of the papers were tagged with ‘03504: Genetics – Plant’) and not being too specific (just one paper was tagged with ‘64500: Paleobiology’). Concept code terms ranged in prevalence from 11.6% for ‘51512: Plant physiology – Reproduction’ to 1.04% for ‘13014: Metabolism – Nucleic acids, purines and pyrimidines’. Selected categories that have waxed and waned as research areas over the years are highlighted. All categories and publications may be explored online at http://bar.utoronto.ca/50YearsOfArabidopsis/. Figure and caption extracted from Provart et al. (2016). Arabidopsis thaliana exhibits many ideal characteristics that position it as a prominent model organism for studying plant biology across multiple fields. Mainly, with its short generation time, A. thaliana can be manipulated in a wide range of experimental settings and easily grows on various artificial media, facilitating its practical handling (Figure 1.1,Meyerowitz 1987). However, what really makes A. thaliana an interesting model organism is its suitability for genetic studies (Meyer,2014;Provart et al.,2016; Weigel,2011). In this regard, A. thaliana has one of the smallest genomes of any flowering plant, with few repetitive sequences, little intergenic spacer DNA and packed into five chromosomes (Leonelli,2007). Furthermore, A. thaliana’s prolific seed production through self-pollination greatly facilitates molecular and quantitative genetic studies. All these features have reinforced A. thaliana as a model organism in plant research, becoming one of the most thoroughly studied of all known flowering plants with more than 54,000 published papers during the las 50 years (Figure 1.2;Provart et al. 2016). 3 1. Introduction Figure 1.3.: (a) General distribution map of A. thaliana (L.) Heynh. and the buffer area (knotted line). The shaded area indicates a, more or less, contiguous distribution of the species; points indicate single occurrences, and open circles imprecise indications of geographical occurrence. (b) Climatically modelled range of A. thaliana. The legend indicates the number of suitable monthly values for growth of the species. In the black area, all twenty-four monthly climatic values of temperature and precipitation support the occurrence of the plant. Figure and caption extracted from Hoffmann (2002). 4 1.2. Intraspecific variation through transplant experiments and genetics analysis Apart from being a valuable model plant in genetics, A. thaliana has also been gaining value in evolutionary and ecological studies (Weigel,2011). A. thaliana is a worldwide distributed plant (Figure 1.3;Hoffmann 2002) occurring in multiple types of habitats, from Mediterranean mixed forests to Scots pine, cedar and holm-oak forests, and to high-altitude grasslands. The origin of the species dates back to 1 Mya – 500 kya in Africa (Durvasula et al.,2017), an origin supported by the large predicted population sizes in this continent compared to those in Eurasia (Fulgione and Hancock,2018). Currently, it has been postulated that A. thaliana was split into its actual linages in two main events. First, the African expansion 120 – 90 kya ago split the populations into relict and nonrelict lineages that spread out to Eurasia. Second, a divergence between European and Asian sub-clusters took place ca. 45 ka ago around the Caucasus region (Durvasula et al., 2017;Zou et al.,2017). The species’ vast distribution, coupled with rapid contractions events attributed to rapid climate changes over glacial periods, has determined the great natural variation observed in both phenotypic traits and genetic composition (Beck et al., 2008;Pic´o et al.,2008). 1.2. Intraspecific variation: the importance of characterizing variation through transplant experiments and genetics analysis Variation in physiological traits among natural accessions was first noticed by Laibach (1943) when he pointed out the importance of A. thaliana as a model plant. In fact, the presence of natural variation enables us to study key aspects of evolutionary and population ecology at the molecular level (Koornneef et al.,2004;Shindo et al.,2007; Weigel,2011). For decades, researchers have focused on untangling the putative drivers, which underlie the natural variation observed in the phenotypic and genetic attributes among A. thaliana accessions (Horton et al.,2012;Koornneef et al.,2004;The 1001 Genome Consortium,2016;Weigel,2011). This interest in the study of natural variation has led to increased experimentation in natural settings, aimed at capturing as much variation as possible in populations of A. thaliana and its close relatives (Borevitz and Nordborg,2003). Recently, researchers have employed hundreds of accessions to establish a network of experiments aimed at investigating variations in morphological traits and fitness, with the goal of quantifying the extent of natural variation across the geographical distribution of A. thaliana (Figure 1.4;Horton et al. 2012;The 1001 Genome Consortium 5 1. Introduction Figure 1.4.: A catalogue of variation for A. thaliana. Two studies document genome sequences, gene expression and molecular modifications to DNA for more than 1,000 varieties of this plant species. Figure and caption extracted from Savolainen and Lascoux (2016). 2016). Intraspecific variation among populations across the species’ distribution arises from a combination of stochastic events (e.g.: genetic drift; Blanquart et al. 2012;Castilla et al. 2020), ecological constraints (e.g.: dispersal limitations; Castilla et al. 2020;Garant et al. 2005;Willi et al. 2007), and natural selection enhancing local adaptation (Exposito- Alonso et al.,2019;Kawecki and Ebert,2004;Leimu and Fischer,2008;M´endez-Vigo et al.,2011). Common garden and transplant experiments represent one approach to quantify intraspecific variation, aiming to identify phenotypic differences in plant traits throughout the developmental cycle. The translocation of populations across different environments has also revealed the rich natural intraspecific variation among A. thaliana accessions, such as their potential for adaptation to a variety of environments according to their phenological and morphological traits (˚ Agren and Schemske,2012;Exposito-Alonso et al., 2018;Fournier-Level et al.,2011). These studies have primarily focused on examining not only fitness-related traits (e.g., seed production and germination success) but also morphological traits (e.g., number of leaves, rosette and stalk size) in order to grasp 6 1.3. Intraspecific variation: the key to better understand species’ distribution patterns the extent of natural variation. Some interesting findings in A. thaliana include the correlation between natural variation in seed response to chilling for germination, flowering time and senescence which mark the seasonal life-history syndromes (Mart´ınez-Berdeja et al.,2020); the correlation between the latitudinal gradient of ecological strategies in A. thaliana as a result of within-species adaptation to climate (Vasseur et al.,2018b); the heritable genetic variation for responses to the different stresses such as drought, fungi and herbivory (Davila et al.,2017); and the changes in defense metabolite production in relation to the environment (Katz et al.,2021). Overall, common garden and transplant experiments that encompass a wide range of natural variation have been valuable tools for investigating how A. thaliana populations may respond to changes in the environment. In addition, these experiments have also contributed to identify key genes and regulatory pathways involved in important processes, such as biotic and abiotic stress responses and local adaptation to contrasting environmental conditions. Furthermore, these studies may also have important implications for society, for example by improving our understanding of the genetic and evolutionary mechanisms underlying crop and tree plantation resilience in a scenario of climate change (Yao et al.,2021). Despite their potential, this sort of experiments also presents some constraints when attempting to capture the complete general picture of relationships between environmental factors and trait variation. These limitations are mainly due to the limited number of simultaneous replicates available to characterize extensive collections of accessions or account for temporal variation (Savolainen et al.,2013). 1.3. Intraspecific variation: the key to better understand species’ distribution patterns Intraspecific variation is also characterized through genotypic analysis. In particular, accessions can be assigned to a certain genetic cluster based on their genetic associations for specific regions of the genome, and genetic clusters may vary in their degree of ecological and genetic differentiation due to various evolutionary processes (Marcer et al.,2016). Genetic analysis has helped identify a significant portion of the natural variation that may not be discernible through phenotypic characteristics alone, that in turn may be obscured by plasticity and trait convergence. As a result of the long-term collaborative efforts within the Arabidopsis community, A. thaliana has one of the most extensive population samplings, including whole-genome sequencing (Borevitz and Nordborg,2003;The 1001 7 2. Chapter I: An ecological history of the relict genetic lineage of A. thaliana pansion of a single lineage (The 1001 Genome Consortium,2016;Durvasula et al.,2017; Lee et al.,2017;Tsuchimatsu et al.,2017;Zeng et al.,2017;Fulgione et al.,2017;Hsu et al.,2019). Such divergent relict lineages clearly emerge in a phylogenomic tree using A. lyrata as the outgroup, the closest relative of A. thaliana Durvasula et al. (2017); Lee et al. (2017); Tsuchimatsu et al. (2017); Zeng et al. (2017); Fulgione et al. (2017); Hsu et al. (2019). Genomic, epigenomic and transcriptomic analyses have indicated that differentiation between relict and non-relict A. thaliana spans most of the genome, including a larger number of transposon elements-like methylated genes and approximately 20% of genes with differential expression between relict and non-relict accessions (The 1001 Genome Consortium,2016;Kawakatsu et al.,2016). In addition, some degree of genomic hybridization has also been detected in contact zones between relict and non-relict A. thaliana, such as the Iberian Peninsula (Durvasula et al.,2017) and China (Hsu et al., 2019), suggesting a complex demographic history that remains to be disentangled. Thus, relict lineages, in particular the African ones, probably represent the deepest history of A. thaliana. In short, the origin of the species is dated in 1 Mya – 500 kya and the split into the relict and non-relict lineages in 120 – 90 kya, which corresponds to the Abbassia Pluvial that was a warm, humid period in Africa (Durvasula et al.,2017). Interestingly, historical population sizes in African relict A. thaliana appear to have been larger than those in Eurasia (Fulgione and Hancock,2018). On top of that, African relict A. thaliana harbours genetic variation that was lost or never present in Eurasia, which overall confers the relict character to this genetic lineage. For example, genomic analysis of the self-incompatibility locus (S-locus) showed that only relict accessions from northern Africa carry all known major haplogroups, suggesting that the selfing ability of A. thaliana evolved in this region from multiple independent inactivation alleles (Tsuchimatsu et al., 2017;Durvasula et al.,2017). Furthermore, some flowering time genes have also revealed differentiation of relict accessions for the frequency of natural polymorphisms affecting gene function, such as a very low frequency of early loss-of-function alleles of FRIGIDA and multiple functional polymorphisms at FLC, FT, FLM, PHYB, PHYC and LHY genes (M´endez-Vigo et al.,2011;Brennan et al.,2014;Liu et al.,2014;Lutz et al.,2017; de Montaigu and Coupland,2017;James et al.,2018). In addition, other genes have been detected by phenotypic and environmental genome wide-association (GWA) analyses, like the flowering promoters TSF and SOC1, the transcription factors LUG and SLK1 involved in flower development, or the stress regulator EIN2 (The 1001 Genome Consortium,2016; Lee et al.,2017;Tabas-Madrid et al.,2018). 14 2.1. Introduction Overall, an out-of-Africa scenario as a plausible demographic model for A. thaliana history is replacing the classical model of the Central Asian origin of the species, which implied an expansion into Europe and Asia and leaving Africa out from any interpretation (Durvasula et al.,2017;Fulgione and Hancock,2018;Lee et al.,2017;Hsu et al.,2019). So far, the most important nuclei of relict A. thaliana can be found in the Iberian Peninsula and Morocco, which was detected first as a highly differentiated genetic cluster with genome-wide nuclear polymorphisms (Pic´o et al.,2008;Brennan et al.,2014;Marcer et al., 2016) and later confirmed with whole-genome sequencing (The 1001 Genome Consortium, 2016;Durvasula et al.,2017). Interestingly, other pools of relict A. thaliana populations have also been identified in geographic areas as distant as Madeira (Fulgione et al.,2017), western and southern Africa (Durvasula et al.,2017), Tibet (Zeng et al.,2017), the Chinese province of Yunnan (Hsu et al.,2019), and previously known accessions scattered across Atlantic oceanic islands (i.e. Cape Verde and the Canary Islands) and the Mediterranean Basin (i.e. Sicily, Algeria and Lebanon; (The 1001 Genome Consortium,2016;Durvasula et al.,2017;Hsu et al.,2019). However, an important piece of information is totally lacking to figure out how relict A. thaliana has been able to endure to the present day: the ecological history of its relict lineage. We only know that relict A. thaliana populations occur in a diverse array of natural environments with a limited impact of human activities (Marcer et al., 2016;The 1001 Genome Consortium,2016). This means that relict A. thaliana behaves more like an annual plant belonging to different natural plant communities (see further details on its biology in (Montesinos et al.,2009;Pic´o,2012;Manzano-Piedras et al.,2014; G´omez et al.,2018) rather than a weed able to colonise humanised environments from croplands to urban spots. Given the tight association of relict A. thaliana with natural environments and its genomic differentiation, we hypothesise that the relict lineage of A. thaliana in the Iberian Peninsula and Morocco has experienced its own evolutionary trajectory through long-term survival and adaptation to multiple environments across the region. The main reason we suggest this hypothesis has to do with the occurrence of Iberian and Moroccan relict A. thaliana in well-known glacial refugia for animals and plants in southern Europe and northern Africa (Cowling et al.,1996;Comes and Kadereit, 1998;Hewitt,1999;G´omez and Lunt,2007;Hughes et al.,2006;M´edail and Diadema, 2009). The region turned into a major biodiversity hotspot of the Mediterranean Basin (e.g. the Baetic–Rifan hotspot; (Medail and Quezel,1997), where several plant families have experienced intense local diversification and speciation (Rodr´ıguez-S´anchez et al., 2010;Lavergne et al.,2013;Molina-Venegas et al.,2016) and co-existing with old, pre- Mediterranean plant genera that still remain in the area (Herrera,1992). 15 2. Chapter I: An ecological history of the relict genetic lineage of A. thaliana Here, we sought for evidences from the past to disentangle the ecological history of relict A. thaliana across its major distribution range in southwest Mediterranean Basin. Given that small annual plants barely leave any trace in the fossil record, we addressed this issue in a different manner. First of all, we modelled the current potential distribution range of the relict lineage of A. thaliana with species distribution models (SDM) combining climate, remote sensing and soil data and using all known populations identified as relict so far in the Iberian Peninsula and Morocco. We used the resulting distribution range map to explore the association between current habitat suitability of relict A. thaliana and vegetation changes since the Last Glacial Maximum (LGM). To do that, we estimated species turnover and changes in arboreal vegetation from fossil pollen data as well as reconstructions of potential distribution ranges of the most important tree species occurring across the region over broad climatic periods since the LGM. Our hypothesis was that relict A. thaliana has long occurred in the Iberian Peninsula and Morocco. Given the current tight association between the occurrence of relict A. thaliana in natural environments and the inability of the species to behave like a weed (Marcer et al.,2016), we expected to find that current relict A. thaliana was more related to environments that experienced more stability over time. Taking advantage of the vast genomic resources available in A. thaliana, we also searched for genes with putative functional polymorphisms characteristic of the relict lineages and with a low frequency in non-relict lineages, which would also represent additional evidence of the long evolutionary history of relict A. thaliana (Tsuchimatsu et al.,2017;Durvasula et al.,2017). We discussed the results emphasising how relict A. thaliana endured past climatic changes and stressing the fact that relict A. thaliana is a lineage that continues its own evolutionary trajectory. 2.2. Materials and methods 2.2.1. Relict source populations We used our Iberian and Moroccan collection of relict A. thaliana populations (39 Iberian and 20 Moroccan accessions; Figure 2.1,2.2 and Table A.1) previously identified with sets of genome-wide neutral single nucleotide polymorphisms (Pic´o et al.,2008;Brennan et al.,2014;Marcer et al.,2016)) and later confirmed with genome-scale sequencing (The 1001 Genome Consortium,2016;Durvasula et al.,2017;Fulgione et al.,2017;Tabas- Madrid et al.,2018). Although the southwest Mediterranean Basin is the area with the most important distribution of relict A. thaliana, we also compiled the coordinates of all 16 2.2. Materials and methods Figure 2.1.: Geographic location of all known relict A. thaliana populations to date. The number of accessions and the number of populations with coordinates, respectively, are given in parenthesis. A. thaliana populations categorised as relicts from published articles conducting largescale genome sequencing of A. thaliana accessions (Figure 2.1 and Table A.1 in order to illustrate the geographic extension occupied by relict A. thaliana. Thus, there were nine additional populations from Africa (including four Moroccan accessions, four Tanzanian, four South African and one Algerian; (Durvasula et al.,2017), five from Madeira (including 11 accessions; (Fulgione et al.,2017), two from Italy (The 1001 Genome Consortium,2016; Durvasula et al.,2017), one from China (including two accessions; (Hsu et al.,2019), one from Lebanon (Durvasula et al.,2017), one from Tibet (Zeng et al.,2017), one from the Canary Islands (Hsu et al.,2019), and one from Cape Verde (Hsu et al.,2019). Given the scarcity of relict populations located out of the western Mediterranean distribution, final analyses were based on 63 southwest Mediterranean populations, including 39 from the Iberian Peninsula, 23 from Morocco and one from Algeria (Figure 2.2 and Table A.1. There is some confusion in the literature around the exact position of the Algerian population included in genome-scale sequencing studies. To deal with this shortcoming, we used the approximate geographical position of a population found in Algeria decades ago by (Qu´ezel and Barbero,1989). 17 2. Chapter I: An ecological history of the relict genetic lineage of A. thaliana Figure 2.2.: Geographic location of Iberian and African relict A. thaliana populations and paleobotanical sites used in this study. Photographs illustrate seven representative environments where relict A. thaliana occurs: montane Scots pine forests (Fun), Mediterranean mixed forests (Sca), high altitude grasslands (Vel), Mediterranean scrublands (Gra), Atlas cedar forests (Ket), and Mediterranean holm-oak forests (Ifr). The location of populations onto the map is indicated by numbers. 2.2.2. Species distribution models We modelled the current potential distribution range of relict A. thaliana to identify the environmental and ecological requirements of this genetic lineage in southwest Mediterranean Basin. The study area to model the potential distribution range of relict A. thaliana was defined by the known distribution of relict populations across the region accumulated in over 15 years of field campaigns and molecular work. Thus, the study area comprised the whole Iberian Peninsula, bordered by the Pyrenees to the north and by the central and western Maghreb to the south, including vast regions of Morocco, Algeria and Tunisia. Nevertheless, to avoid arbitrary decisions to delimit the study area in North Africa, somewhere between the Tell Atlas in Algeria, the Atlas mountains in Morocco and the Sahara desert, its southern limit was defined by the area affected by the so-called 18 2.2. Materials and methods greening of the Sahara during the Mid-Holocene, 6 kya (years before present), as reported by paleoclimatological reconstructions indicating the occurrence in North Africa of grassland and scrubland vegetation in wetter environmental conditions (De Noblet-Ducoudr´e et al.,2000;Compton,2011). In our opinion, this represents a reasonable distinction between potentially suitable habitat for A. thaliana and the inhospitable Sahara desert in an appropriate time frame for this study. Overall, the study area ranged between 28.0 ° N – 43.8 ° N and 11.2 ° E – 11.7 ° W, totalling 1,389,229 km2, of which 583,123 km2corresponded to the Iberian Peninsula and 806,106 km2to North Africa (Figure 2.2). We selected a set of climate, habitat, soil and ecophysiological predictors from different publicly available data sources at a resolution of 1 km2. It must be noted that we checked for co-linearity among all predictors to be below —0.8—, as measured by the Pearson’s correlation coefficient and using the whole set of data fed to the distribution model (i.e. background points plus relict occurrences). First of all, a total of six of 19 bioclimatic variables were selected from WorldClim v.2 (Fick and Hijmans,2017): BIO1 (annual mean temperature), BIO3 (isothermality), BIO4 (temperature seasonality), BIO12 (annual precipitation), BIO15 (precipitation seasonality) and BIO18 (precipitation of the warmest quarter). Vegetation data were represented by two variables: the percentage of crops and the percentage of non-arboreal vegetation. The percentage of crops referred to agricultural lands that reduce the occurrence of relict A. thaliana (Marcer et al.,2016). The percentage of crops was obtained from MODIS/Terra and Aqua Combined Land Cover (MCD12Q1) v.6 data product (Friedl and Sulla-Menashe,2015) combining croplands and croplands/natural vegetation mosaics using the global land cover classes of the International Geosphere-Biosphere Programme (IGBP). We obtained a layer of percentage of crops by assigning the most frequent value for each grid cell (1 km2) for the period 2003 – 2016. The percentage of non-arboreal vegetation was also obtained following the same procedure but using the MODIS/Terra Vegetation Continuous Fields (MOD44B) v.6 data product (DiMiceli et al.,2015). The percentage of non-arboreal vegetation was chosen because of the affinity of A. thaliana for forest openings and scrublands (Montesinos et al., 2009;G´omez et al.,2018). In this case, however, we used the mean value for each grid cell, instead of the most frequent value, to generate the single layer for this variable between 2003 and 2016, as this is a continuous variable. 19 2. Chapter I: An ecological history of the relict genetic lineage of A. thaliana Next, we obtained the Leaf Area Index (LAI), an ecophysiological trait, from the MODIS/Terra Leaf Area Index (MOD15A2H) v.6 data product (Myneni et al.,2015). Similarly, we used mean values for each grid cell for the period 2003 – 2016. LAI is broadly defined as the amount of leaf area (m2) in a canopy per unit ground area (m2). For this reason, it is an important structural parameter of the vegetation as it is directly and proportionally related to major parameters, such as photosynthetic primary production and evapotranspiration. LAI was highly correlated with other important ecophysiological variables, e.g. Net Photosynthesis, Gross Primary Production or Enhanced Vegetation Index (r >0.85 in all cases), and for this reason we only used LAI to represent the ecophysiological traits in our set of relict A. thaliana populations. When generating our layers from MODIS, we only considered those pixels labelled as highest quality in their respective quality control layers. Only 0.56% of pixels were excluded for this reason. Finally, a single topsoil pH layer was generated by merging data provided by the Geographical Database from Eurasia v.4 (JRC-ESDAC,2015) and the Harmonized World Soil Database v.1.2 (FAO, IIASA, ISRIC, ISSCAS, JRC,2012). We included this variable because it is an important factor accounting for the presence of A. thaliana in the Iberian Peninsula, particularly making A. thaliana increasingly scarce as topsoil pH becomes more basic (Marcer et al.,2016)). We used the software Maxent v.3.3.3k (Phillips et al.,2006;Elith et al.,2011) to model the potential distribution of the genetic lineage of A. thaliana in the Iberian Peninsula and the Maghreb. The model was fitted using the default parameters given by Maxent. We applied the 10-fold cross-validation area-under-the-curve (AUC) metric to evaluate prediction performance. Finally, we used all occurrence points to fit the final model and project it to the whole study area. All data preparation and modelling was carried out using the R computing environment (R Core Team,2017) using several packages (i.e. broom, corrgram, dismo, dplyr, ggplot2, jsonlite, magrittr, ncf, raster, rasterVis, readr, readxl, rgdal, rJava, sp, spdep, stringr, and tidyr) and the Geospatial Data Abstraction Library (GDAL and OGR,2015). Given that the Iberian Peninsula and the Maghreb are two separate sub-regions, we tested for differences between regions for all potential distribution range predictors with one-way ANOVA with SPSS v.25 statistical software (IBM Corp,2017). 20 2.2. Materials and methods 2.2.3. Fossil pollen data As A. thaliana does not appear in the fossil pollen record, we used indirect means to indentify the historical factors related to past vegetation dynamics accounting for the current potential distribution range of relict A. thaliana across the study region. To this end, we quantified the vegetation changes occurred over several millennia from paleobotanical records from all over the Iberian Peninsula and the Maghreb. Overall, we gathered data from 80 paleobotanical sites reporting fossil pollen data (Figure 2.2 and Table A.2). We obtained the database of paleobotanical records from two sources. First, we used pollen records publicly available in the European Pollen Database (EPD), resulting in 50 sites from Spain, five from Tunisia, three from Algeria, three from Portugal and two from Morocco. Second, we conducted a bibliographic search on recent published articles reporting fossil pollen data from the region that were not available in the EPD. We contacted the corresponding authors and requested for published raw data in the same format than that posted on the EPD (i.e. pollen counts per species, sediment sample depth and calibrated age). Eventually, we added to our database eight sites from Spain (Gil-Romera et al., 2008;Burjachs et al.,2016;Gonz´alez-Samp´eriz et al.,2017;Jones et al.,2018;L´opez-S´aez et al.,2015), four from Morocco (Campbell et al.,2017;Cheddadi et al.,2017), three from Andorra (Beaulieu et al.,2005;Miras et al.,2010,2015) and two from Portugal (Santos and Go˜ni,2003;Fletcher et al.,2007). Overall, the 80 paleobotanical sites were scattered across the Iberian Peninsula and the Maghreb, from sea side to montane environments (mean altitude ± SE = 884.0 ± 72.1 m.a.s.l., range = 0–2,406 m) and intermingled with relict A. thaliana popula tions (Figure 2.2; mean ± SE distance between relict populations and paleobotanical sites = 652.6 ± 390.7 km; range = 4.0 – 1686.7 km). Important for our study, data from paleobotanical sites included a calibrated chronology mostly expressed in calendar years before 1950 (cal yr BP). Our dataset contains sediment cores that were calibrated by different methods and age-depth models. However, the assessment of the value of these methodologies is clearly beyond our reach. The length of the age-depth series were highly variable among paleobotanical sites, ranging from a low of 228 cal yr BP to a high of 44,370 cal yr BP (mean ± SE = 10,061.2 ± 986.3 cal yr BP; Table A.2. Given the heterogeneity of paleobotanical data, in terms of depth, calibrated age, continuity of the chronology and number of paleobotanical sites, we divided the chronologies in three broad climatic periods (Blockley et al.,2012;Gonz´alez-Samp´eriz et al.,2017): present (up to 2 kya), Holocene (11.7 – 2 kya) and a final period combining the Younger 21 2. Chapter I: An ecological history of the relict genetic lineage of A. thaliana Dryas and the Allerød (14.0 – 11.7 kya). Each one of these periods is characterised by extensive climatic variability. The present includes the Medieval Climatic Anomaly (900–1300 CE), the Little Ice Age (1300–1850 CE) and the increasing warming since the Industrial Revolution. The Holocene (11.7 – 2 kya) is typically divided into various periods during which the vegetation experienced several important changes across the study region by virtue of climatic variation and human activities that modified the environment (e.g. fire and agriculture). Finally, the last period, including the Younger Dryas and the Allerød (14.0 – 11.7 kya), combined two distinctive periods, such as the Younger Dryas (12.9 – 11.7 kya) characterised by cold glacial conditions, and the Allerød (14.0 – 12.9 kya) characterised by warmer and moister environmental conditions. Beyond 14 kya, the number of paleobotanical sites was too low to attempt any solid statistical analysis. For each climatic period and paleobotanical site, we estimated the species turnover using the metric developed by MacArthur and Wilson (1963) and modified by Diamond (1969). Species turnover was calculated as the proportion of species that differs between two time points (i.e. the sum of species gained and lost divided by the total species observed in both time periods). To do that, we used time windows of 100 years within each climatic period from which we computed mean and SD species turnover. We tested narrower and wider time windows and results were quite consistent, although 100-year time windows captured better the high heterogeneity of fossil pollen data. Across paleobotanical sites, the mean ± SD number of 100-year time windows was 13.1 ± 5.9 (range = 1–23) for the present, 30.5 ± 22.1 (range = 1–93) for the Holocene, and 7.4 ± 5.0 (range = 1–17) for the Younger Dyras and the Allerød. For each climatic period, we also estimated the mean and SD of the percentage of arboreal pollen with and without pollen from the genus Pinus, which can be quite dominant across the region, using the same 100-year time windows. The percentage of arboreal pollen was estimated including pollen data from trees and shrubs. The drop or increase of arboreal pollen over time provides hints on the temporal change in vegetation structure and dynamics. We did not consider species abundance as the interpretation of this parameter is always controversial requiring mathematical approaches to reconstruct pollen diversity at each paleobotanical site (Weng et al.,2006). For each climatic period, we calculated the correlation between current habitat suitability of relict A. thaliana estimated with Maxent (mean and SD) and (i) the percentage of arboreal pollen with and without pollen from the genus Pinus (mean and SD across 100-year time windows), and (ii) species turnover (mean and SD across 100-year time windows). Habitat suitability was estimated for a circular area (5 km radius) around 22 2.2. Materials and methods the geographic position of each paleobotanical site. We sought to figure out whether the habitat suitability of relict A. thaliana in present times was related to the vegetation dynamics registered in paleobotanical sites at a local scale for the three broad climatic periods. We repeated the analyses with other radii (10, 20 and 50 km around each paleobotanical site) and results were quite consistent although larger areas obviously tended to blur the patterns. Given the spatially explicit nature of variables, we conducted Pearson’s correlation with and without correcting for the spatial autocorrelation of variables (Dutilleul’s modified t-test; (Dutilleul et al.,1993) using SAM software v.4.0 (Rangel et al., 2010). Spatial autocorrelation may result in several statistical flaws, such as violation of the assumption of independent error or bias in the estimation of correlation coefficients and variances (Dutilleul et al.,1993). The comparison of results with and without the correction for spatial autocorrelation allows the assessment of the extent of the problem caused by spatial autocorrelation of data. 2.2.4. Distribution maps of dominant Iberian tree species over time We explored the relationship between the potential distribution range of the most important tree species with that of relict A. thaliana. Tree species can be considered as representative taxa of their plant communities. Thus, we used them to pinpoint those plant communities more or less related to the current potential occurrence of relict A. thaliana. To this end, we compared the current potential distribution range of relict A. thaliana with the current potential distribution range of the most important tree species occurring in the Iberian Peninsula using data published elsewhere (Benito-Garz´on et al., 2007). The rasterized maps (1 km2grid format) were kindly provided by Marta Benito- Garz´on including the present potential distributions and two past potential distributions at two particular points in time: the Last Glacial Maximum (LGM; 21 kya) and the Mid- Holocene (6 kya). Past potential distributions were used to infer the dynamics of such plant communities over the last millennia. We excluded the Maghreb from these analyses as we did not have similar data from that region. In short, the present distributions of 17 tree species were determined from occurrences reported in the Spanish Forest Map and the Portuguese Invent´ario Florestal Nacional. Bioclimatic envelope models were used to estimate potential distributions using (i) Agroclimatic Characterization of Spain’s Provinces and the Portuguese AGRIBASE 23 2. Chapter I: An ecological history of the relict genetic lineage of A. thaliana Figure 2.5.: Relationship between niche overlap values for relict A. thaliana with 17 tree species in the present and stability in plant communities (see text for details) in two broad climatic periods: (A) since the Last Glacial Maximum (LGM) and (B), between the LGM and the mid-Holocene (MH). non-relict genomes. However, this analysis identified a new haplogroup (haplogroup C; Brennan et al. 2014), which was found at high frequency in Moroccan relicts (40%) as well as in the Algerian and Lebanon relicts (Figure 2.6A and B). This haplogroup was differentiated from the rest of accessions by three amino acids, and it appears at very low frequency in non-relict accessions from the rest of world (<1%). On the other hand, analysis of TSF sequences from the available genomes showed that most non-relict accessions from all over the world (97%) carried the major haplogroup (haplogroup A; Figure 2.6C and D). By contrast, 54% of the relict accessions had the socalled haplogroup B differentiated by multiple non-coding polymorphisms, which might affect TSF regulation. Although the functional effect of the missense and non-coding mutations differentiating CRY2-C and TSF-B need to be demonstrated, these haplogroups highlighted the potential usefulness of relict accessions to unveil novel functional genetic diversity. 30 2.4. Discussion Figure 2.6.: Neigbour-joining trees showing the genetic relationships among 1216 genomic sequences of CRY2 (A) and TSF (C). Bootstrap values larger than 50% are shown. Branches corresponding to partitions reproduced in less than 25% bootstrap replicates are collapsed. Geographical distribution of major CRY2 (B) and TSF (D) haplogroups. Relict and non-relict accessions carrying different haplogroups are shown in different colours as described in the legends of maps. 2.4. Discussion The collaborative efforts of a large and multidisciplinary research community merging extensive field sampling and whole-genome deep sequencing represent an outstanding asset to unravel the evolutionary history of the annual plant A. thaliana. The uninterrupted work of the community for decades has recently paid off by providing new insights into the early history of A. thaliana indicating the existence of relict lineages with a wide geographic distribution spanning up to three continents (The 1001 Genome Consortium,2016; Durvasula et al.,2017;Tsuchimatsu et al.,2017;Zeng et al.,2017;Fulgione et al.,2017; Hsu et al.,2019). However, the ecological history of relict A. thaliana is of paramount 31 2. Chapter I: An ecological history of the relict genetic lineage of A. thaliana importance to complete the picture provided by whole genome sequencing. This study attempted to reconstruct the ecological history of relict A. thaliana by seeking for indirect evidences of the recent past accounting for the occurrence of relict A. thaliana in southwest Mediterranean Basin, the most important nuclei of known relict populations. Our hypothesis was that relict A. thaliana has long survived and adapted across the region, as inferred for other annual plant species in the Iberian Peninsula (Fern´andez-Mazuecos and Vargas,2013;Bobo-Pinilla et al.,2018), probably in areas that experienced more stability over time. In fact, it is quite remarkable the occurrence of some relict A. thaliana populations growing in or very close to populations of wellknown relict plant species that long endured confined to mountain refugia in Spain and northern Africa (Carri´on et al., 2010), such as relict forests of silver birch (Betula alba) in central Spain, relict firs (Abies pinsapo,A. marocana and A. tazaotana) in southwest Spain and the Rif, the Sictus tree (Tetraclinis articulata) in southeast Spain, or the Atlas cedar (Cedrus atlantica) in the Atlas Mountains. Nevertheless, such coincidences, although stimulating and inspiring for this work, required of additional support to become more than mere observations. To test our hypothesis, we first modelled the current potential distribution range of relict A. thaliana in Morocco and the Iberian Peninsula, which can also be interpreted in terms of current habitat suitability across the region. Our results were consistent with a previous model for relict A. thaliana in the Iberian Peninsula (Marcer et al.,2016), even after the inclusion of some new elements, such as the extension of the study area to North Africa and the use of new vegetation parameters from remote sensing data. Overall, we confirm a known ecological attribute of relict A. thaliana, which is its preference for natural habitats and probably its difficulty to colonise humanised environments (Marcer et al., 2016;The 1001 Genome Consortium,2016). As a matter of fact, recent genomic analyses gave rise to the hypothesis that admixture of locally adapted relict A. thaliana with non-relict lineages helped the latter to colonise multiple natural and humanised habitats through rapid demographic expansion until dominating most of the species’ current distribution range (The 1001 Genome Consortium,2016;Lee et al.,2017;Hsu et al.,2019), which is what we observe today. Next, we used the current habitat suitability of relict A. thaliana in Morocco and the Iberian Peninsula, which showed small differences in some minor predictors of habitat suitability but not in the most relevant ones (Table 2.1), to explore its relationship with estimates of past vegetation dynamics over the last millennia. The analysis of fossil pollen data from sediment cores located across the study region indicated that paleobotanical sites with higher current habitat suitability for relict A. thaliana exhibited lower 32 2.4. Discussion temporal variation in vegetation dynamics for different broad climatic periods (Figure 2.4, as measured by species turnover and by the proportion of arboreal pollen. Despite the complexity of fossil pollen data, even in the interpretation of vegetation dynamics over millennia in the region (Carri´on et al.,2010), the existence of such relationships suggests that relict A. thaliana required of environments characterised by higher stability over time to make it through the present times. This conclusion lines up with the preference of known relict A. thaliana populations for natural environments in which they have long survived. It is worth emphasising that a previous analysis comparing historical climatic conditions between relict and non-relict A. thaliana populations suggested that temperature changed less for relict than for non-relict A. thaliana since the LGM (The 1001 Genome Consortium,2016). Although such analysis was carried out by comparing climatic stability of only 26 relict versus 1109 non-relict accessions worldwide, without taking the spatial autocorrelation of climatic data into account, and assuming that all populations (relict and non-relict) remained in the same location since the LGM, which is a huge assumption for an annual plant, our results might provide some support to that suggestion if we assume that vegetation changes have chiefly been mediated by variation in temperature. On top of the indirect evidence reported in this study on the relationship between current habitat suitability of relict A. thaliana and vegetation dynamics over the last millennia, we provided additional support to our original hypothesis by exploring the niche overlap between relict A. thaliana and some tree species in the Iberian Peninsula that represent whole plant communities. Our results allowed the identification of those plant communities with a higher (e.g. Quercus ilex subsp. ballota and Q. pyrenaica widely distributed in central Iberian Peninsula) and lower (Quercus ilex subsp. ilex and Pinus uncinata restricted to northeast Iberian Peninsula) affinity with relict A. thaliana in the present times. Interestingly, niche overlap between relict A. thaliana and tree species was significantly positively correlated with tree niche stability since the LGM (Figure 2.5), suggesting that relict A. thaliana in the present times is tightly associated to plant communities whose distribution experienced less variation over millennia. This result is pointing in the same direction than the previous one and reinforces the view that the persistence of relict A. thaliana through time depended to a large extent on the stability of the plant communities where it occurred. It must be emphasised that the dynamics of non-relict A. thaliana is thought to be completely different, as it is hypothesised that nonrelict lineages have experienced a demographic expansion colonising multiple environments (Durvasula et al.,2017;Fulgione and Hancock,2018;Lee et al.,2017;Hsu et al.,2019). Given the wide array of environments where relict A. thaliana can be found in Morocco 33 2. Chapter I: An ecological history of the relict genetic lineage of A. thaliana and the Iberian Peninsula in the present times, which is also supported by the important niche overlaps detected in our study between relict A. thaliana and several tree species (Figure A.1), it is plausible to suggest that the fate of relict populations was tightly determined by the fate of plant communities in which the plant was able to occur. This highlights the long evolutionary history of relict A. thaliana, which would have had time to adapt to very contrasting environments since the split of the species into relict and non-relict lineages. Finally, our survey of A. thaliana genome sequences available from world-wide accessions yielded additional support to the strong genomic differentiation between relict and non-relict A. thaliana at the functional level. This is important because relict A. thaliana is expected to harbour specific genetic variation resulting from its long evolutionary history (Tsuchimatsu et al.,2017;Durvasula et al.,2017). In particular, we identified highly differentiated haplogroups of CRY2 and TSF genes, which are present at high frequency among relict accessions but nearly absent in non-relict genomes (Figure 2.6). As shown with the two genes analysed here, it can be expected that genetic studies of relict accessions will provide natural functional polymorphisms that otherwise cannot be detected. Thus, we stress the need to step up our efforts to sample in depth the rest of African, Middle East and Asian regions where relict accessions have been detected to increase sample size in these regions and to further investigate the genetic differentiation of relict A. thaliana from a functional perspective. This will provide hints on the genetic mechanisms underlying the long evolutionary history of relict A. thaliana. In fact, the evolutionary history of relict A. thaliana still hides many fascinating secrets that remain to be discovered, such as the genetic relationships among relicts from distant geographic locations, the extent of admixture between relict and non-relict accessions and its evolutionary consequences, and the amount and structure of functional genetic variation that makes relict A. thaliana so particular. 34 3. Chapter II: Genomic relatedness chiefly drives multi-trait phenotypic differentiation in the annual A. thaliana across its wide native distribution 3.1. Introduction Widely distributed species represent important components of biodiversity across extensive geographic regions as well as relevant examples of adaptation to contrasting environmental conditions. In fact, the extent of phenotypic differentiation of widely distributed species reflects the phenotype-environment matching across the broad environmental gradients of their native distributions. In widely distributed plants, large-scale reciprocal transplant experiments revealed a trend for ecotypic variation, as individuals exhibited better performance and selection advantage at their home sites (Joshi et al.,2001;Becker et al.,2006;˚ Agren and Schemske,2012;Wilczek et al.,2014;Johnson et al.,2022). Furthermore, ecotypic variation on widely distributed plants may also emerge at the molecular level. In this case, alleles associated to higher fitness tend to be local alleles linked to local environmental factors or to environments more similar to those at the experimental settings (Fournier-Level et al.,2011;Lasky et al.,2018;Exposito-Alonso et al.,2019). Thus, adaptive variation (i.e. variation mediated by natural selection) appears as an important force in widely distributed plants, shaping genetic and phenotypic differentiation patterns across broad environmental gradients. Despite the potential of reciprocal transplant experiments to estimate individual performance from home and away locations, reciprocal transplant experiments may also 35 3. Chapter II: Genomic relatedness drives multi-trait phenotypic differentiation have some constraints due to the limited number of simultaneous experiments to characterize extensive collections of accessions or if the strength of selection varies over time (Savolainen et al.,2013). One way to overcome these issues is to estimate phenotypic differentiation in common environments between hundreds of accession pairs or subsets of populations spanning the entire distribution (Rehfeldt et al.,2002;Wilczek et al.,2009; Atwell et al.,2010;Neuffer,2011;Mart´ınez-Berdeja et al.,2020;Kahl et al.,2021;Villellas et al.,2021). In this varying in their genetic makeup, using genome-wide markers, and coming from different locations and environments across broad distributions. By using extensive collections of accessions, we can test the effect of three nonexclusive putative drivers of phenotypic differentiation across broad distributions: (i) geographic distance between accession pairs, (ii) genomic relatedness between accession pairs, and (iii) environmental distance between accession pairs. The first driver is geographic distance between accession pairs, i.e. isolation-by-distance (Slatkin,1993;Rousset,1997; Sexton et al.,2014), which provides insight into the effect of dispersal limitation on differentiation as effective dispersal is expected to decrease at increasing distances. Although isolation-by-distance acts primarily on the genetic differentiation of populations, as dispersal is not able to counteract the effects of random genetic drift, the same process may eventually contribute to phenotypic differentiation by intensifying local adaptation to increasingly varying environmental conditions with distance (Nosil et al.,2009). The second driver of phenotypic differentiation is the genomic relatedness between accession pairs, which deals with genomic similarity in phenotypic differentiation due the shared ancestry between accession, i.e. the probability that they share sets of genes that are identical-by-descent; (Weir et al.,2006). In this case, phenotypic differentiation between accessions pairs with greater genomic distances may increase due to overall lower trait variance transmitted among relatives over generations. This pattern should become stronger for traits with higher additive genetic variance, which is the fraction of the genetic variance transmitted from parents to offspring. In the case of adaptive variation in traits, selection may maintain beneficial allelic variants in populations, which is normally detected in genome-wide association studies (Atwell et al.,2010;Huang et al.,2010; Brachi et al.,2013;The 1001 Genome Consortium,2016) and in expression patterns at genes accounting for adaptive variation (Lopez-Arboleda et al.,2021). Thus, the effects of genomic distance on phenotypic differentiation may be particularly more intense in traits with high additive genetic variance under strong selection. The final driver of phenotypic differentiation has to do with environmental distance between accession pairs, i.e. isolation-by-environment (Sexton et al.,2014;Wang 36 3.1. Introduction and Bradburd,2014). In the case of accession pairs from similar environments, phenotypic differentiation may decrease due to the prevalent role of local adaptation to the same environmental attributes over the rest of forces influencing the genetic composition of populations. The most important problem of estimating environmental distance deals with the complexity of the environmental variation that matters for phenotypic differentiation in multiple traits, particularly when dealing with an extensive collection of accessions across a large, environmentally diverse area. One way to address this issue is to employ multiple environmental variables (i.e. climate, vegetation, soil) to estimate overall environmental suitability across the entire distribution with species distribution models. This is because lower and higher suitability areas across the distribution indicate lower and higher environmental adequacy for the study species, respectively. Importantly, species distribution models allow the identification of the most insightful variables accounting for environmental suitability, allowing us to pinpoint the specific environmental determinants of suitability. Here, we analyzed the main effects of geographic, genomic and environmental suitability distances on phenotypic differentiation in a wide array of quantitative traits of the annual plant A. thaliana across its native distribution. We estimated geographic, genomic and environmental suitability distances using geographic coordinates, whole-genome sequences and suitability from a species distribution model, respectively. For phenotypic differentiation, we used data from published studies characterizing variation in phenotypic traits in common environments with hundreds of accessions from all over the A. thaliana’s native distribution. In particular, phenotypic data encompassed various traits involving plant architecture and life history traits (The 1001 Genome Consortium,2016;Vasseur et al.,2018a;Exposito-Alonso et al.,2019;Mart´ınez-Berdeja et al.,2020), chemical defenses (glucosinates, Katz et al. 2021) and response traits to abiotic and biotic agents (Davila et al.,2017). Given the role of natural selection shaping up genetic variation in phenotypic traits in A. thaliana (Fournier-Level et al.,2011;˚ Agren and Schemske,2012; Wilczek et al.,2014;Exposito-Alonso et al.,2019), we expect a predominant effect of genomic distance on phenotypic differentiation. However, as we analyzed various traits varying in their additive genetic variance and presumably affected by selection in different ways and intensity across a broad distribution, we cannot completely rule out the potential contribution of geographic and environmental suitability distances to phenotypic differentiation. In addition, the A. thaliana’s native distribution is very diverse in climate (Hoffmann,2002), orography and land use, which affects several important plant attributes, 37 3. Chapter II: Genomic relatedness drives multi-trait phenotypic differentiation such as trait plasticity and competitive ability (Estarague et al.,2022). Furthermore, the extent of effective dispersal (natural or human-mediated dispersal), the balance between natural selection and random genetic drift affecting populations, and local adaptation to contrasting environments across the A. thaliana’s broad distribution, may also affect the relative contribution of geography, genomics and environmental suitability to multi-trait phenotypic differentiation across large spatial scales. To handle the heterogeneity of the A. thaliana’s native distribution, we also dissected the effects of geographic, genomic and environmental suitability distances on phenotypic differentiation into separate regions (i.e. Iberia, Great Britain, Fennoscandia and Eurasia). We discuss the results in terms of the difficulties and complexities posed by widely distributed plants to identify the drivers of phenotypic differentiation in multiple traits. 3.2. Methods 3.2.1. Study area The study area encompassed the best-sampled native distribution area of the annual plant A. thaliana (L.) Heyn. (Brassicaceae) in Eurasia and northern Africa (11 ° W–115 ° E, 30 ° N–80 ° N; Figure 3.1). We excluded areas of the native distribution in Africa (East and South Africa) and Asia (Himalayas and East Asia) due to the low number of accessions available and/or the lack of phenotypic data for accessions outside the selected area. 3.3. Environmental suitability modeling We compiled A. thaliana occurrences from the literature and several databases, such as GBIF, the 1001 Genomes Project, the Versailles Arabidopsis Stock Center and the authors’ accession collection from the Iberian Peninsula and Morocco (Table 3.1). For GBIF, we only used occurrences with a spatial uncertainty lower than 1,500 m, as coordinate uncertainty may considerably affect the outcomes of species distribution models (Marcer et al.,2022). We reduced occurrence bias by keeping a single occurrence per 30-arcsecond grid cell, matching the resolution of the environmental predictor layers. The final dataset used to generate the environmental suitability model contained 4,719 occurrences (Figure 3.1a). 38 3.3. Environmental suitability modeling Figure 3.1.: Occurrences and regions of A. thaliana used in this study. (a) Occurrences of A. thaliana used in the species distribution model (N = 7,419). Each occurrence represents a unique 30 arc-second grid cell. (b) Accessions with genomic and phenotypic data for each region (Iberia, N=188; Great Britain, N=57; Fennoscandia, N=228; Eurasia, N=483; Whole area, N=974). The whole range is the ensemble set of all accessions and regions. 39 3. Chapter II: Genomic relatedness drives multi-trait phenotypic differentiation Table 3.4.: Number of A. thaliana accessions included in the analyses per trait and region. Sources: [1] The 1001 Genome Consortium (2016); [2] Mart´ınez-Berdeja et al. (2020); [3] Davila et al. (2017); [4] Vasseur et al. (2018a); [5] Exposito-Alonso et al. (2019); [6] Katz et al. (2021). For variables of source [5], MH, ML, TH and TL letters should be read as: M - Madrid experiment site, T - T¨ubingen experiment site, H - high precipitation treatment, and L - low precipitation treatment. Variable description Ref. Great Britain Fennoscandia Iberia Eurasia Whole area Architecture Rosette leaf number [1] 49 202 179 444 891 Rosette diameter [1] 49 202 179 444 891 Rosette dry mass [4] 17 71 67 236 398 Rosette branch number [1] 49 202 179 444 891 Inflorescence size MH [5] 22 98 160 213 501 Inflorescence size ML [5] 22 98 160 213 501 Inflorescence size TH [5] 22 98 160 213 501 Inflorescence size TL [5] 22 98 160 213 501 Cauline leaf number [1] 49 202 179 444 891 Fruit production [4] 17 71 67 236 398 Seed production MH [5] 22 98 160 213 501 Seed production ML [5] 22 98 160 213 501 Seed production TH [5] 22 98 160 213 501 Seed production TL [5] 22 98 160 213 501 Development Germination time MH [5] 22 98 160 213 501 Germination time ML [5] 22 98 160 213 501 Germination time TH [5] 22 98 160 213 501 Germination time TL [5] 22 98 160 213 501 Base germination [2] 33 138 90 198 465 Germination (4 days, 4 ° C) [2] 33 138 90 198 465 Germination (4 days, 10 ° C) [2] 33 138 90 198 465 Germination (15 days, 4 ° C) [2] 33 138 90 198 465 Germination (15 days, 10 ° C) [2] 33 138 90 198 465 Germination (32 days, 4 ° C) [2] 33 138 90 198 465 Germination (32 days, 10 ° C) [2] 33 138 90 198 465 Rosette growth rate [4] 33 138 90 198 398 Survival to flowering MH [5] 22 98 160 213 501 46 3.3. Environmental suitability modeling Table 3.4.: (continued) Variable description Ref. Great Britain Fennoscandia Iberia Eurasia Whole area Survival to flowering ML [5] 22 98 160 213 501 Days to bolting [2] 33 138 90 198 465 Days to flowering [2] 33 138 90 198 465 Flowering time 1001 [1] 49 202 179 444 891 Flowering time MH [5] 22 98 160 213 501 Flowering time ML [5] 22 98 160 213 501 Flowering time TH [5] 22 98 160 213 501 Flowering time TL [5] 22 98 160 213 501 Days to fruiting [2] 33 138 90 198 465 Days to harvest [2] 33 138 90 198 465 Glucosinolates Allyl [6] 55 129 162 371 729 Branched chained [6] 55 129 162 371 729 Benzoyloxypropyl [6] 55 129 162 371 729 Indol-3-ylmethyl [6] 55 129 162 371 729 2S-OH-3-Butenyl [6] 55 129 162 371 729 3-OH-4-Pentenyl [6] 55 129 162 371 729 3-Butenyl [6] 55 129 162 371 729 3-Hydroxypropyl [6] 55 129 162 371 729 4-Hydroxybutyl [6] 55 129 162 371 729 4-Hydroxy-Indol-3-ylmethyl [6] 55 129 162 371 729 3-Methylsulfinylbutyl [6] 55 129 162 371 729 4-Methylsulfinylbutyl [6] 55 129 162 371 729 5-Methylsulfinylpentyl [6] 55 129 162 371 729 6-Methylsulfinylhexyl [6] 55 129 162 371 729 7-Methylsulfinylheptyl [6] 55 129 162 371 729 8-Methylsulfinyloctyl [6] 55 129 162 371 729 3-Methylthio [6] 55 129 162 371 729 4-Methylthiobutyl [6] 55 129 162 371 729 7-Methylthioheptyl [6] 55 129 162 371 729 8-Methylthiooctyl derivate [6] 55 129 162 371 729 4-Methoxy-Indol-3-ylmethyl [6] 55 129 162 371 729 8-Methylthiooctyl [6] 55 129 162 371 729 47 3. Chapter II: Genomic relatedness drives multi-trait phenotypic differentiation Table 3.4.: (continued) Variable description Ref. Great Britain Fennoscandia Iberia Eurasia Whole area 4-Pentenyl [6] 55 129 162 371 729 Response Pieris effects [3] 13 44 6 77 142 Pieris & Botrytis effects [3] 13 44 6 77 142 Pieris & Drought effects [3] 13 44 6 77 142 Thrip damage [3] 13 44 6 77 142 3.3.3. Analyses For each region, we evaluated the correlation between pairs of distance matrices for each driver of phenotypic differentiation (geography, genomics and environmental suitability) with a Mantel test using the ‘mantel’ function of the vegan R package Oksanen et al. (2020). We used Pearson correlations and tested significance with 999 permutations. To assess the relative contribution of geography, genomics and environmental suitability to phenotypic differentiation, we used Maximum-Likelihood Population-Effects (MLPE) mixed models (Clarke et al.,2002). We employed generalized least squares using the ‘gls’ function of the nlme R package (Pinheiro et al.,2020) and included a correlation structure with the ‘corMLPE’ function of the corMLPE R package (Pope,2020). For each phenotypic trait and region, we fitted a model using its pairwise distance matrix as response variable and the corresponding pairwise distance matrices for geography, genomics and environmental suitability as predictor variables. A pair of identifiers corresponding to each accession in each pair of the distance matrix provided the correlation structure of the datasets. Prior to analyses, we normalized all variables by subtracting the mean and dividing by the standard deviation. We generated a model selection table with all combinations of predictor variables, without interactions, using the ‘dredge’ function of the MuMIn v.1.46.0 R package (Barto´n,2020). We obtained the coefficient estimate for each predictor using model averaging with the ‘model.avg’ function of the MuMIn R package. Coefficients were averaged across all models if they had a maximum correlation coefficient of —0.65— between predictors and a maximum delta AIC of 2 with respect to the model with the lowest AIC. 48 3.4. Results 3.4. Results 3.4.1. Environmental suitability model We used MIAMaxent to estimate environmental suitability across A. thaliana’s native range based on climate, soil and vegetation variables (Figure 3.3). Based on the fraction of total variation explained by each variable (FTVA; Halvorsen et al. 2015), three variables accounted for 87.3% of the variation in our model: soil pH, the thermicity index, and the enhanced vegetation index (EVI) range (Figure 3.4a). These variables had interpretable patterns of variation to understand the distribution of environmental suitability across A. thaliana’s native range (3.4b). In particular, differences in environmental suitability showed two distinct general patterns across the A. thaliana’s native distribution. First, we found a latitudinal pattern with higher values concentrated at mid latitudes and with a sharp decrease at latitudes above 60 ° N (Figure 3.3b and 3.3c). Below 35 ° N and 45 ° N in the southwestern and southeastern part of the distribution, respectively, the pattern was more heterogeneous as environmental suitability exhibited sharp fluctuations along the southern margin (Figure 3.3b and 3.3c). Second, we detected a more homogeneous longitudinal pattern of environmental suitability in the western twothirds of the study area in comparison to the eastern one-third (Figure 3.3a and 3.3b). This higher homogeneity of environmental suitability in the western than in the eastern part of the study area matches the variation in occurrence density from West to East (Figure 3.1a). As far as regions were concerned, Great Britain showed low environmental suitability in northern Scotland and higher environmental suitability across the rest of the region (Figure 3.3b). In Iberia, unsuitable areas occurred along large areas of the eastern part of the peninsula, whereas the rest of the region exhibited higher environmental suitability (Figure 3.3b). In Fennoscandia, suitable areas were mostly concentrated along coastal areas with a sharp decrease in environmental suitability inland, except in some parts of central and southern Sweden (Figure 3.3b). Eurasia followed the same pattern described before for the Whole range (Figure 3.3b). Overall, Great Britain was the region with higher mean and median environmental suitability and Fennoscandia had the lowest values (Table 3.5). 49 3. Chapter II: Genomic relatedness drives multi-trait phenotypic differentiation Figure 3.3.: Environmental suitability (ES) maps for A. thaliana expressed in probability ratio output (PRO) values, which represent the suitability of an average training occurrence. Mean ES distribution (b) for the 30-arcsecond model aggregated to a 5-minute grid, (a) along the longitudinal axis and (c) along the latitudinal axis. ES’ standard deviation distribution (e) for the 30-arcsecond model aggregated to a 5-minute grid, (d) along the longitudinal axis and (f) along the latitudinal axis. 50 3.4. Results Figure 3.4.: Variable importance and response curves of the environmental suitability model for A. thaliana. (a) Relative contribution expressed as the fraction of total variation accounted for (FTVA) by each environmental variable included in the environmental suitability model. (b) Response curves of environmental suitability expressed in probability ratio output (PRO) values, which represent the suitability of an average training occurrence, along each environmental variable gradient. Table 3.5.: Summary of environmental suitability values per region. Summary values are in PRO units calculated from the areas enclosed by a buffer of about 50 km around the convex hull defined by the accessions for each region (Figure 3.1b). Region Minimum Maximum Mean SD Median Great Britain 0.105 9.162 2.177 1.192 1.991 Fennoscandia 0 5.644 0.404 0.49 0.196 Iberia 0.007 6.896 1.573 0.975 1.504 Eurasia 0 1749.581 2.003 23.547 0.879 Whole area 0 1749.581 1.731 20.88 0.785 51 3. Chapter II: Genomic relatedness drives multi-trait phenotypic differentiation 3.4.2. Drivers of phenotypic differentiation We conducted Mantel tests between driver pairs of phenotypic differentiation to evaluate the structure of our predictor variables in each region. The results indicated a positive correlation among geographic, genomic and environmental suitability distances in practically all regions (P <0.02 in all cases; Table 3.6). The exceptions were given by the Mantel tests between genomics and environmental suitability in Great Britain (P = 0.26) and between geography and environmental suitability in Eurasia (P = 0.72) and in the Whole range (P = 0.10). We fitted MLPE mixed models to discern the contributions of geographic, genomic and environmental suitability distance between accession pairs on phenotypic differentiation in A. thaliana. Overall, genomic distance was the driver with further significant effects on phenotypic differentiation, followed by geography and environmental suitability (Figure 3.5 and Table 3.7). The number of traits with genomics, geography or environmental suitability as the only significant driver of phenotypic differentiation was 75, 23 and 3, respectively, including all analyses from all regions (Figure 3.5). In Iberia, ge- Table 3.6.: Correlations between geography, genomics and environmental suitability distance matrices in each region according to Mantel tests. Region Predictor 1 Predictor 2 Statistic P-value Great Britain Geography Suitability 0.195 0.001 Great Britain Genomics Geography 0.091 0.018 Great Britain Genomics Suitability 0.032 0.257 Eurasia Genomics Geography 0.424 0.001 Eurasia Genomics Suitability 0.167 0.001 Eurasia Geography Suitability -0.016 0.721 Fennoscandia Genomics Geography 0.559 0.001 Fennoscandia Geography Suitability 0.479 0.001 Fennoscandia Genomics Suitability 0.393 0.001 Iberia Genomics Geography 0.291 0.001 Iberia Geography Suitability 0.166 0.001 Iberia Genomics Suitability 0.052 0.01 Whole range Genomics Geography 0.319 0.001 Whole range Genomics Suitability 0.101 0.001 Whole range Geography Suitability 0.022 0.104 52 3.4. Results Figure 3.5.: Normalized effects of geography, genomics and environmental suitability on phenotypic differentiation in A. thaliana per region and trait group given by Maximum-Likelihood Population-Effects (MLPE) mixed models. See Table 3.3 for trait description. 53 3. Chapter II: Genomic relatedness drives multi-trait phenotypic differentiation nomics and geography had similar contributions to phenotypic differentiation, whereas in Great Britain, genomics was clearly the driver with the highest relevance for phenotypic differentiation (Figure 3.5). We detected some negative significant coefficients for some drivers in some traits (Figure 3.5). We suspect, however, that they were the result of responses to correlated predictor variables rather than having a biological meaning, as negative coefficients mostly appeared in models combining drivers with positive and higher coefficients (Figure 3.5). As far as traits were concerned (Figure 3.5), glucosinolates and development traits accumulated the majority of significant effects of genomics, as the only driver, on their differentiation (35 and 28 significant coefficients including all regions, respectively), followed by architecture (8 significant coefficients) and response traits (4 significant coefficients). Geography, as the only driver, was the predictor with more significant coefficients for architecture (10 significant coefficients) and development traits (8 significant coefficients), and less important for response traits (3 significant coefficients) and glucosinolates (2 significant coefficients). Finally, environmental suitability was only the only driver in development traits (2 significant coefficients). When more than one driver was significant in the same analysis, they concentrated on glucosinolates (67 significant coefficients) and development (44 significant coefficients) traits and genomics was the driver with the greatest effects (Figure 3.5). Table 3.7.: MLPE summary of model coefficients per phenotypic category and region. Columns mean, median sum and N show the mean value, the median value, the sum of all values, and the number of values for all coefficients different from zero for each region. The most important predictors per phenotypic category and region are given in pale blue. Region Predictor Coefficient estimates summary Mean Median Sum N Architecture Great Britain Genomics 0.191 0.157 0.574 3 Iberia Genomics 0.062 0.047 0.309 5 Iberia Geography 0.057 0.042 0.284 5 Iberia Suitability -0.028 -0.028 -0.028 1 Fennoscandia Genomics 0.162 0.154 0.486 3 Fennoscandia Suitability 0.017 0.017 0.017 1 Fennoscandia Geography -0.056 -0.072 -0.224 4 Eurasia Geography 0.054 0.052 0.271 5 Eurasia Genomics 0.043 0.044 0.129 3 54 3.4. Results Table 3.7.: (continued) Region Predictor Coefficient estimates summary Mean Median Sum N Eurasia Suitability 0.021 0.021 0.021 1 Whole range Genomics 0.045 0.014 0.408 9 Whole range Geography 0.016 0.012 0.157 10 Whole range Suitability -0.007 -0.007 -0.015 2 Development Great Britain Genomics 0.255 0.227 3.565 14 Iberia Geography 0.127 0.119 1.906 15 Iberia Genomics 0.085 0.08 1.611 19 Iberia Suitability -0.037 -0.038 -0.519 14 Fennoscandia Geography 0.098 0.111 1.077 11 Fennoscandia Suitability 0.085 0.075 0.853 10 Fennoscandia Genomics 0.064 0.071 0.514 8 Eurasia Geography 0.203 0.243 0.813 4 Eurasia Suitability 0.033 0.033 0.033 1 Eurasia Genomics 0.019 0.016 0.113 6 Whole range Genomics 0.114 0.114 2.389 21 Whole range Geography 0.061 0.054 1.154 19 Whole range Suitability 0.026 0.016 0.289 11 Glucosinolates Great Britain Genomics 0.151 0.132 2.874 19 Iberia Genomics 0.068 0.063 1.163 17 Iberia Geography 0.02 0.018 0.219 11 Iberia Suitability -0.006 -0.015 -0.043 7 Fennoscandia Genomics 0.121 0.073 2.666 22 Fennoscandia Suitability -0.021 -0.021 -0.214 10 Fennoscandia Geography -0.092 -0.084 -1.192 13 Eurasia Genomics 0.063 0.055 1.381 22 Eurasia Suitability -0.006 -0.009 -0.09 15 Eurasia Geography -0.007 -0.009 -0.088 13 Whole range Genomics 0.082 0.079 1.813 22 Whole range Geography 0 0.004 0.008 18 Whole range Suitability -0.003 -0.003 -0.038 11 Response 55 4. Chapter III: The value of intraspecific variation to grasp niche differentiation brings new adaptive characteristics, we may observe niche divergence across a species’ distribution through adaptation (Warren et al.,2008) despite the observed species’ tendency towards niche conservatism (Wiens et al.,2010). Nevertheless, separating the relative effects of processes shaping the current patterns of the genetic structure and their relationship with the environment is a challenging task. The major difficulties are the lack of solid genomic data for many species along with a biased representation of their global distribution and limited tools for niche comparison. Thereby, the mechanisms by which heterogeneous natural environments influence species’ genetic variation and structure remain largely unknown. Some insights on the role of environment in niche divergence can be obtained from the ecological characteristics at the sites where populations occur. This is because genetic structure analysis provides the means to explore how populations can be split into more meaningful genetic units, which may also differ in ecological characteristics, to study niche differentiation. Furthermore, measures of niche overlap are commonly used as an indirect way to test niche environmental differentiation among genetic clusters (Broennimann et al.,2012). In particular, niche overlap calculations have been performed using both the environmental (E-space) and the geographical space (G-space) to characterize niche variation among groups, mostly among species (Herrando-Moraira et al.,2019), but also among genetic clusters (Marcer et al.,2016;Mart´ınez-Minaya et al.,2019). Most studies calculate the environmental space of the genetic clusters using ordination analysis by estimating the overlap among the areas occupied in multidimensional spaces. Recently, environmental niche models have also been used to test the influence of environmental meaningful characteristics in geographic space by considering the habitat suitability overlap predicted across the study area (Smith et al.,2019). Both approaches can help discern the level that is ecologically meaningful and present some evidence on niche differentiation as the common species-level approach to model habitat suitability has been questioned (Smith et al.,2019). Nevertheless, not all the species meet the characteristics to test niche divergence at the intraspecific level, such as climatic variability across the study range, good genetic characterization of populations and strong evidence for local adaptation. Model plant species, such as A. thaliana, are appropriate study systems to explore general ecological questions, owing to the extensive knowledge compiled over a few decades in relation to the intraspecific genetic, phenotypic and geographic ranges. Arabidopsis thaliana is a worldwide annual plant whose native distribution exhibits a great climatic variation across Eurasia (Hoffmann (2002); Chapter II). In addition, A. thaliana populations 62 4.1. Introduction are strongly structured genetically at different geographical scales (Pic´o et al.,2008; Platt et al.,2010;Brennan et al.,2014;The 1001 Genome Consortium,2016;Durvasula et al.,2017). The availability of genome-wide single-nucleotide polymorphism (SNP) data among a geographically diverse collection including hundreds of accessions (The 1001 Genome Consortium,2016) allows the analysis of A. thaliana’s genetic structure across its wide native distribution. Finally, several studies have shown evidence of local adaptation in A. thaliana populations at several geographic scales in common garden (M´endez-Vigo et al.,2011;Manzano-Piedras et al.,2014;Davila et al.,2017) and reciprocal transplant experiments (Fournier-Level et al.,2011;Exposito-Alonso et al.,2018). Importantly for this study, other works have shown evidence of niche adaptation to specific environments stressing the importance of accounting for intraspecific niche divergence in the characterization of A. thaliana’s habitat (Marcer et al.,2016). In particular, the relict genetic cluster seems relegated to some undisturbed environments (Marcer et al.,2016;The 1001 Genome Consortium,2016;Toledo et al.,2020), whereas non-relict lineages are those exhibiting a weedy behavior (Lee et al.,2017). Overall, the complex A. thaliana’s demographic and evolutionary history, including an African origin with a split in relict and non-relict lineages across multiple glacial periods, raises the chances to identify population niche divergence. In this study, we quantified the environmental differences among previously characterized genetic clusters of A. thaliana and the relationship between environmental and genetic distances among genetic clusters as a means to detect evidence for niche divergence and adaptation at the intraspecific level. For this purpose, we characterized environmental differences and the degree of environmental overlap among genetic clusters using two approaches. First, we used ordination analysis to describe the environmental space where each population occurs and the environmental overlap among genetic clusters. Second, we modeled the habitat suitability of each genetic cluster and estimated their overlap to describe the environmental characteristics across the geographical space. In addition, we used the genetic structure and the genetic distance analysis from previous studies (Chapter II) on A. thaliana across its native distribution to calculate the genetic overlap among genetic clusters. Finally, we estimated the relationship between genetic and environmental overlaps for both environmental space and geographic space approaches. Here, we addressed the following specific questions: (1) how did genetic clusters across the A. thaliana’s native distribution diverge in environmental conditions and habitat suitability? (2) Is the genetic overlap related to the environmental differences characterizing the distribution of genetic clusters? And (3), can we consider A. thaliana’s genetic clusters as ecological meaningful units? 63 4. Chapter III: The value of intraspecific variation to grasp niche differentiation 4.2. Materials and methods 4.2.1. Study area and source populations Our study area was the native range of A. thaliana in Eurasia and northern Africa (Figure 4.1) adjusted by accessions with available whole-genome sequence data. We used 739 A. thaliana accessions with genomic data structured in nine genetic clusters (The 1001 Genome Consortium,2016;Durvasula et al.,2017). Of these, 680 were available from The 1001 Genome Consortium (2016) and 59 from the Iberian and Moroccan collection including relict A. thaliana accessions (39 Iberian and 20 Moroccan accessions; Chapter I). Figure 4.1.: Occurrences of A. thaliana for each genetic cluster (Central Asia, N = 53; Central Europe, N = 95; Germany, N = 47; Iberian Peninsula, N = 102; Italy-Balkan-Caucasus, N = 65; North Sweden, N = 27; Relicts, N = 63; South Sweden, N = 73; Western Europe, N = 77, Total N = 711). Each occurrence represents a unique 30 arcsecond grid cell. 64 4.2. Materials and methods 4.2.2. Habitat suitability modeling We used the presence-only maximum entropy model as implemented by the MIAMaxent R package (Chapter II,Vollering et al. 2019). For occurrence records, we filtered the 739 available accessions to a single one per 30-arcsecond grid cell (approximately 1 km2), resulting in a final number of 711 occurrence points. For background points, we selected 10 000 random points within a 400-km buffer around species’ occurrences to better represent the environmental conditions of the species’ range. We spatially weighed the selection of these random background points with a sampling-effort proxy layer to further reduce the sampling bias. We created the bias layer as follows. First, we downloaded the available occurrences of the Brassicaceae family from GBIF (filtered by preserved specimen category and coordinate uncertainty below 2000 m; accessed 15 July 2019) to create a 2D kernel density layer using the ‘bkde2D’ function of the KernSmooth R package (Wand,2020). Second, we created a layer with distances to A. thaliana occurrences to increase the odds of selecting points that putatively represent absences. Finally, we subtracted the distance to occurrences layer from the sampling-effort proxy layer and scaled it between 0 and 1. For predictor variables, we used the same set of eight variables that were previously selected to model the habitat suitability of A. thaliana (Chapter II), namely mean temperature of the wettest quarter, precipitation of the coldest quarter, precipitation seasonality, thermicity index, potential evapotranspiration (PET) of the wettest quarter, enhanced vegetation index (EVI) range, percentage of soil sand and soil pH. We used the function ‘extract’ form the terra R package (Hijmans,2020) to assign each variable value to each occurrence and background point. We split the 711 occurrence points into nine separate datasets corresponding to their genetic clusters and modelled each one separately with the same 10 000 background points for each model. We used the same methodology as in Chapter II to calibrate the models for each genetic cluster. In this case, each spatial fold had the same number of occurrences and final models were also selected based on maximum Area Under the Curve (AUC). However, if a large gap existed between the highest and the second-highest AUCs, we selected the more parsimonious model to avoid overfitting. We built the final models again using the whole set of occurrences and background points to predict suitability across the whole study area. Finally, we normalized the Probability Ratio Output (PRO) predicted values as given by MIAMaxent (Vollering et al.,2019) for the matrix overlap comparison. 65 4. Chapter III: The value of intraspecific variation to grasp niche differentiation 4.2.3. Environmental space analysis We used Principal Components Analysis (PCA) to estimate the environmental distances between genetic cluster pairs. We used presence data from the genetic clusters and all non-presence grid cells available for the species within a 400-km buffer around cluster occurrences. We extracted the environmental values for each of the eight variables using the ’extract’ function and accordingly classified them as presence or background. We then run the PCA analysis with the function ’dudi.pca’ in the ade4 R package (Thioulouse et al., 2018). Finally, we calculated the kernel density values for each genetic cluster along the two principal components of the PCA analysis with the function ‘kde2d’ of the MASS R package (Ripley et al.,2013). 4.2.4. Genetic space analysis For the genetic space analysis, we used a genetic distance matrix between pairs of accessions as previously estimated (Chapter II). We assigned each accession to its corresponding genetic cluster and applied multidimensional scaling (MDS) to the pairwise genomic distance matrix using the function ’cmdscale’ of the stats R package (R Core Team, 2022). We then used the function ’kde2d’ of the MASS R package (Ripley et al., 2013) to calculate the kernel density values for each genetic cluster along the three MDS principal dimensions for each pair of dimensions. 4.3. Genetic cluster overlap We evaluated the degree of overlap between genetic cluster pairs based on genetic distances, environmental space (PCA environment) and habitat suitability (environmental habitat maps). We used the Schoener’s D statistic (Schoener,1968) and the Warren’s I (Warren et al.,2008) statistics to calculate the degree of overlap in each of the three cases. In the case of the genetic distance, we calculated the overlap for the three main axes. Finally, we summarized the three genetic overlap matrices values to obtain a genetic overlap index accounting for most of the variation in the genetic space. 66 4.4. Results 4.3.1. Determining factors of genetic cluster overlap We modeled how the genetic distance overlap between genetic clusters was determined by both environmental and habitat suitability overlaps using a Maximum-Likelihood Population-Effects model (MLPE; Clarke et al. 2002). We fitted each model using the genetic overlap index matrix as response variable and the environmental and habitat suitability overlap matrices as predictor variables. Models were fitted with generalized least squares using the function ‘gls’ of the nlme R package (Pinheiro et al.,2020) and including a correlation structure with the function ‘corMLPE’ of the corMLPE R package (Pope,2020). The correlation structure was provided to the model as a pair of identifiers corresponding to each genetic cluster. We used the Akaike’s Information Criterion (AIC) for model selection with a threshold value of 2 to consider a model better than the other. 4.4. Results 4.4.1. Genetic clusters habitat suitability We used presence-only maximum entropy models with climate, soil and vegetation predictor variables to estimate the habitat suitability of each of the nine genetic clusters across the Eurasian and northern African native ranges of A. thaliana. The nine genetic clusters were named as Central Asia, Central Europe, Germany, Iberian Peninsula, Italy-Balkans- Caucasus, North Sweden, Relicts, South Sweden, and Western Europe. We used spatial cross-validation and an optimized regularization (alpha) parameter for each genetic cluster model. Alpha values ranged from 0.0005 to 0.01, whereas the values of the adjusted AUC values ranged from 0.89 to 0.99 (Table 4.1). All variables were selected in the final models except for percentage of soil sand. Based on the fraction of total variation accounted for (FTVA) by each environmental variable, we observed large differences among genetic clusters. Each genetic cluster had a quite distinct set of variables contributing to their habitat suitability (Figure 4.2). Across all genetic clusters, the FTVA values for each variable in decreasing order are shown in Table 4.2 67 4. Chapter III: The value of intraspecific variation to grasp niche differentiation Table 4.1.: Summary table for the habitat suitability models of the genetic clusters. Columns: GROUP - Genetic cluster; NO. - number of occurrences; MODEL FORMULA - Best model; ALPHA - Regularisation parameter; and AUCc - Adjusted Area Under the Curve. Group NO Model formula alpha aucc Western Europe 77 Thermicity Index + pH + PET of wettest quarter + EVI range 0.001 0.96 South Sweden 73 Thermicity Index + PET of wettest quarter + Precipitation seasonality + Precipitation of coldest quarter 0.001 0.99 Relicts 63 Thermicity Index + Precipitation seasonality + pH + PET of wettest quarter 0.01 0.96 North Sweden 27 pH + mean temperature of wettest quarter + precipitation of Coldest quarter + precipitation of seasonality 0.0005 0.99 Italy-Balkan -Caucasus 65 pH+ PET of wettest quarter + mean temperature of wettest quarter + precipitation seasonality 0.001 0.93 Iberian Peninsula 102 Thermicity index + pH + PET of wettest quarter + mean temperature of wettest quarter + Precipitation seasonality 0.01 0.94 Germany 47 Thermicity index + EVI Range + Precipitation seasonality 0.001 0.93 Central Europe 95 Thermicity index + pH + EVI Range + mean temperature of wettest quarter 0.01 0.9 Central Asia 53 Thermicity index + pH + PET of wettest quarter + EVI Range + Precipitation seasonality 0.01 0.93 Table 4.2.: FTVA values for each variable across all genetic clusters in decreasing order. Variable No. models FTVA (mean ± SD) soil pH 7 34 ± 15 thermicity index 7 29 ± 18 mean temperature of the wettest quarter 4 25 ± 6 precipitation seasonality 6 25 ± 25 PET of the wettest quarter 6 21 ± 20 EVI range 5 13 ± 7 precipitation of the coldest quarter 2 9 ± 8 68 4.4. Results Figure 4.2.: Variable importance and response curves of the predictor variables in the environmental suitability models for each genetic group. A) Relative contribution of each environmental variable expressed as the fraction of total variation accounted for (FTVA). B) Response curves of habitat suitability for each predictor variable expressed in normalized probability ratio output (PRO) values. 69 4. Chapter III: The value of intraspecific variation to grasp niche differentiation Figure 4.3.: Habitat suitability maps for A. thaliana expressed in probability ratio output (PRO) values. Red lines represent the buffer area for the genetic clusters occurrences. Black lines represent the buffer area for the complete occurrence dataset. Shadowed areas are outside of the study area. The final predicted distributions also varied among genetic clusters (Figure 4.3). As expected, the highest suitability areas of the Italy-Balkans-Caucasus, Iberian Peninsula and Relicts spread out across the Mediterranean Basin. The Iberian Peninsula and Relicts showed a high suitability in the Iberian Peninsula and northern Africa, although the Relicts also extended to the Balkans, the Anatolian Peninsula, and to a lesser extent, East Asia. The habitat suitability of the Italy-Balkans-Caucasus extended across the northern Mediterranean Basin and throughout the Black and Caspian Basins. Western Europe, Central Europe and Germany occupied areas extending from western to central Europe. Western Europe also spread out across northern Iberian Peninsula, whilst Central European extended to the Caucasus. North Sweden and South Sweden occupied parts of northern Europe, although North Sweden also occupied areas in northeastern Asia. Finally, Central Asia extended across the mountain areas of East Asia and reached some western parts of Europe. 70 4.4. Results 4.4.2. Genetic clusters environmental space We performed Principal Component Analysis (PCA) to calculate the density areas of the environmental space (climate, soil and vegetation) of each genetic cluster across the study area (Figure 4.4). We retained the first four principal components that explained 77.96% of the total variation (PC 1: 28.42 %, PC 2: 24.89 %, PC 3: 12.49 %, and PC 4: 12.16 % of variation explained; Figure 4.4A). The variables that contributed at least 10% to each principal component were different for each PC (Figure 4.4A). For PC 1, the variables were PET of the wettest quarter (37.6%), mean temperature of the wettest quarter (30.5%) and precipitation of the coldest quarter (24%). For PC 2, the variables were soil pH (38.5%), thermicity index (22.6%), precipitation seasonality (19.3%) and EVI range (13.8%). For PC 3, the variables were thermicity index (38.3%), EVI range (18.6%), precipitation seasonality (16.4%) and mean temperature of the wettest quarter (13.3%). Finally, for PC 4, only percentage of soil sand contributed as much as 94.5%. The use of environmental space differed among genetic clusters (Figure 4.4B). Central Asian and Central Europe inhabited areas with higher PET of the wettest quarter and higher mean temperatures of the wettest quarter than the more coastal and Mediterranean clusters (e.g., Relicts and Iberian Peninsula). Genetic clusters from northern areas (e.g., North Sweden, South Sweden and Germany) inhabited locations with lower pH, higher EVI range and higher precipitation than the Mediterranean clusters (e.g., Relicts and Iberian Peninsula). Mediterranean clusters (Italian-Balkans-Caucasus, Relicts and Iberian Peninsula) occurred in areas with higher pH, higher precipitation seasonality, lower PET and lower mean temperatures in the wettest quarter. Finally, Western Europe exhibited a preference for high precipitation of the coldest quarter. 4.4.3. Genetic differentiation We applied multidimensional scaling (MDS) to the matrix of genomic distances to represent the genetic relationship between A. thaliana’s genetic cluster pairs (Figure 4.5). We retained the first three components of the MDS, which explained 49.65 %, 34.41 % and 15.94 % of variation, respectively. For the first and second dimensions, we found that Central Asia, North Sweden and Germany were strongly differentiated from the rest of the European clusters (Figure 4.5C). This strong differentiation also emerged between the first and the third axes for the Relicts and, to a lesser extent, for the Iberian Peninsula (Figure 4.5A). Finally, the relationship between the second and third axes, while 71 4. Chapter III: The value of intraspecific variation to grasp niche differentiation ture, we encourage the use of the whole species or genetic clusters, respectively, which in any way require extensive sampling and genotyping. In fact, extensive sampling across wide geographical ranges is mandatory if we pursue to grasp species’ niche divergence in depth. From a technical viewpoint, our study also provided some novel insightful aspects. For example, suitability comparisons may have two main advantages for overlap analysis. First, they simplify the environmental space by selecting variables according to the variation explained by the presence prediction. Second, they incorporate the geographical space into the projection while predicting suitability areas. This geographical component is important to establish its relationship with the genetic factors, as ecological factors, such as dispersal, may have an enormous weight shaping up the genetic relationships among accessions. Furthermore, comparisons in the environmental space tackles the environmental difference among accessions at the known locations instead of the potential ones, which highlights the real variation among genetic clusters and compares it with suitability projections. Thereby, combining both approaches allowed us to relate the intraspecific variation with the environmental and geographical spaces. On top of the indirect evidence reported in this study on niche divergence among genetic clusters of A. thaliana, we provided additional support to our original goals by exploring the relationship between genetic and environmental distances of the species at the intraspecific level. We found evidence on the relationship between the genetic and the environmental overlap for the genetic cluster level expressed as suitability. Nevertheless, we did not find significant relationships between the genetic and the environmental overlap in the environmental space. These results stressed the connection between the environmental characteristics and the genetic differentiation, as well as the importance of considering the environmental space and the geographical space by projecting the suitable environmental conditions for niche comparisons. Suitability predictions combine both the environmental space and the geographical space, that in the specific case of A. thaliana explain better the differences found among genetic clusters. Although the genetic-environmental overlap correlations were week, they set the way forward in determining further significant spatial modeling units using intraspecific variation on top of the species level. 78 5. Discussion Scientists have long endeavored to understand the driving factors behind the spatiotemporal patterns of biodiversity in order to minimize or prevent the effects of contemporary mass extinctions; it is estimated that over one million species are at risk of extinction in the Anthropocene era (Hannah,2012). One of the main concerns, given the current anthropogenic carbon emissions trajectory, is that biodiversity loss is predicted to accelerate with the increased warming by the rapid climate change (Lambers,2015). Historically, a large proportion of biodiversity studies has been devoted to the species level. However, in recent years, the research community has started paying attention to the role that variability within species plays in defining overall biodiversity patterns (Smith et al., 2019). The genetic and phenotypic diversity that exists within species, or intraspecific variation, has become an asset for conservation, being considered as an equally important component of biodiversity (Roches et al.,2021). The study of intraspecific variation will, not only help understand fundamental questions in ecology, such as natural selection and phenotypic diversity, but also the current crisis on several ecological processes (Bolnick et al.,2011;Des Roches et al., 2018). Intraspecific variation plays a critical role in species and ecosystem adaptation and resiliency, whilst providing important contributions to society in the form of food security and medicinal resources (Roches et al.,2021;Stange et al.,2020). Therefore, we urgently need to understand the importance of the forces shaping up the patterns of intraspecific variation. One approach is by looking into patterns of genetic and phenotypic differentiation and their major drivers of variation. This thesis takes a biogeographical approach to delve into the great intraspecific variation exhibited by Arabidopsis thaliana across its native range, with particular attention to the ecological and historical processes shaping the species’ distribution. Additionally, this thesis uses novel approaches to quantify the relative contributions that genetics, phenotypic and environmental elements play when defining the spectrum of intraspecific variation, which overall improve large-scale habitat distribution estimations. 79 5. Discussion Furthermore, in this thesis, I offered novel tools to study the distributional boundaries of an example of a widely distributed plant and its intraspecific units, which may have the potential to be applied to other vastly distributed organisms. 5.1. Ecological history and genetic differentiation of relict Arabidopsis thaliana populations in the Iberian Peninsula: insights into adaptation and persistence Genetic studies have emphasized the large intraspecific variation exhibited by A. thaliana populations, and have also attempted to reconstruct the species’ evolutionary history under the postulate of an early split in Africa followed by two consecutive colonizations in Europe (Fulgione et al.,2017;Fulgione and Hancock,2018). The current paradigm states that some populations have long remained as relicts in natural areas enduring glacial periods in Eurasia (Beck et al.,2008;Pic´o et al.,2008;Brennan et al.,2014;Fulgione et al.,2017), in contrast to the recent spread of non-relict populations associated with human disturbances that likely evolved after the out-of-Africa events (Lee et al.,2017). In this thesis, I aimed to investigate the ecological history of A. thaliana by specifically focusing on the genetically differentiated relict lineage found in the Iberian Peninsula. The hypothesis underlying this research suggested that relict populations in the Iberian Peninsula and North Africa have survived and adapted to diverse environments due to the region’s historical role as a glacial refuge. However, the lack of concrete ecological, historical and biogeographical evidence prompted me to delve deeper into the topic. This thesis provides significant insights into the ecological preferences and historical persistence of relict A. thaliana populations in the Iberian Peninsula (Chapter I). The findings of this chapter validate the species’ affinity for stable habitats and the challenges it faces in adapting to transformed environments, which is consistent with previous studies indicating the higher probability of finding relicts in non-disturbed environments (Marcer et al.,2016;The 1001 Genome Consortium,2016). The observed correlation between the current habitat suitability of relict populations and the long-term stability of vegetation dynamics suggests that temporal stability plays a pivotal role in the survival of this relict lineage. These findings support the hypothesis that the Iberian Peninsula and North Africa acted as glacial refuges for relict A. thaliana populations, as observed also in other plant species inhabiting the region (Agudo et al.,2023). This is in accordance with the patterns shown by close relatives of A. thaliana that also persist in multiple northern 80 5.2. Genomic distance as a key driver of phenotypic differentiation refugia in eastern Central Europe (ˇ Sr´amkov´a Fuxov´a et al.,2017). In addition to ecological factors, the analysis of genetic differentiation between relict and non-relict A. thaliana accessions, with a specific focus on flowering genes, revealed distinct haplogroups within the relict A. thaliana. These haplogroups were predominantly found among relict accessions while they were scarce in non-relict genomes. This genetic differentiation underscores the unique functional genetic variations that have evolved in relict populations over time. These findings, along with the great difference observed in niche dimensions for different genetic clusters in A. thaliana (Chapter III), supports the theory of a complex history with numerous gaps that still need to be bridged. Further steps will be necessary to complete this picture, such as the replication of this study to disentangle the relationship between current presence and historical environmental and vegetation stability for other relict and non-relict A. thaliana lineages in other regions of its native distribution. A complete reconstruction of the evolutionary history of the species will help understand its current distribution and the genetic structure of its extant populations. 5.2. Genomic distance as a key driver of phenotypic differentiation in Arabidopsis thaliana across diverse environments Phenotypic variation can explain many of the differences observed in individual fitness (Nolting et al.,2021), which may have important consequences for population dynamics (Forsman,2014;Schreiber et al.,2011). Some of these phenotypic differences emerge as a result of evolutionary responses manifested as changes in allele frequencies across successive generations (Zhao et al.,2020). However, other observed phenotypic variations occur as plastic responses without changes in the genetic structure between and within generations (Sultan,1987). In Chapter II, I explored the factors that underlie the largescale phenotypic variation observed across the diverse environmental gradients present in the species’ native distribution. The comprehensive analysis of phenotypic traits obtained from published studies, using hundreds of accessions spanning the species’ entire native distribution, shed light into the genetic and evolutionary aspects underlying the observed phenotypic differentiation in the species. 81 5. Discussion This thesis presents evidence of the greater role played by genomic distance in explaining phenotypic differentiation in A. thaliana, surpassing the influence of geographic and environmental suitability distances (Chapter II). I found that genomic distance is the primary driver for differentiation of many traits and has the greatest effect when multiple drivers contribute to phenotypic differentiation. Furthermore, the notable correlations observed between geographic, genomic, and environmental suitability distances, indicate their interconnectedness as drivers of phenotypic differentiation. This suggests that natural selection, rather than random genetic drift, may be playing a prominent role in shaping the phenotypic and genetic structure of the species across its diverse environments. Moreover, it is reasonable to interpret that phenotypic differentiation stems more from a genetic basis rather than a mere plastic phenotypic expression in response to heterogeneous environmental conditions. This observation aligns with the significant genetic structure observed in A. thaliana, which can be attributed to limited dispersal, high self-fertilization, and significant local adaptation in fitness-related traits (Castilla et al.,2020). There are evident regional differences in the strength of correlations between driver pairs, indicating the relative influence of historical, geographic and environmental backgrounds on phenotypic differentiation. Selective pressures change across the distribution range and over time, with greater significance in widely distributed organisms such as A. thaliana. This thesis emphasizes the importance of identifying the appropriate geographical units to accurately determine the putative drivers of phenotypic variation. Furthermore, it may also be necessary to take into account regional differences when setting up new transplant experiments in order to improve our understanding of the patterns of variation that these drivers undergo across the species distribution. I suggest that regional differences can vary for the different traits under investigation, making it worthwhile to explore distinct regional units for each phenotypic trait. Identifying these regional units can lead to accurately identify the drivers for different groups of populations that can experience different evolutionary pressures. Interestingly, I show in this thesis that different groups of traits exhibit distinct patterns of significant contributions in relation to the drivers of phenotypic differentiation (Chapter II). In particular, glucosinolates (compounds produced for defense against herbivores) and development traits show more significant contributions, likely due to their stronger impact on overall fitness as they involve defense abilities and life-cycle completion. In contrast, architecture traits and traits associated with responses to abiotic and biotic agents have fewer significant contributions, possibly due to their secondary role as 82 5.2. Genomic distance as a key driver of phenotypic differentiation targets of selection or to their complex biological responses. These results underscore the importance of considering multiple traits simultaneously at a global scale in order to comprehend general patterns of phenotypic diversification in organisms with wide geographic distributions. In this study, the fundamental role of fitness-related traits indicates that selection may be acting as a major force shaping up intraspecific phenotypic differentiation. In widely distributed organisms, pinpointing the drivers of phenotypic differentiation across their native distributions requires the standardization of phenotypic measurements in multiple common garden experiments and the creation of a unified global trait differentiation database. This titanic task can only be achieved through an extensive collaboration among the plant community members. In the case of A. thaliana, there are valuable antecedents of such world-wide collaborations within the community; i.e., the creation of a detailed map of variation in 1 135 high-quality sequenced natural inbred lines (The 1001 Genome Consortium,2016) along with their metabolome and transcriptome (Kawakatsu et al.,2016). This approach is crucial in order to effectively uncover general patterns in phenotypic differentiation. However, it is important to acknowledge the challenges associated with homogenizing trait databases. Another conclusion of this thesis is that large-scale studies on widely distributed organisms should prioritize comprehensive sampling schemes across the distribution range as well as the use of whole-genome sequences or extensive sets of genome-wide markers. However, given that studies are independent and usually focused on small regions, and that funding is rarely available for such large-scale sampling, I suggest that A. thaliana researchers take into account the bias in the data and prioritize the sampling of understudied areas. Additionally, it is important to include multiple traits to characterize phenomes rather than a limited number of related traits when the objective is to understand phenotypic differentiation. Finally, the environmental heterogeneity across the distribution of widely distributed organisms should be tackled in depth, including climate, soil type and vegetation, to understand its biogeographical patterns. For example, I demonstrate in this thesis that factors, such as soil pH or vegetation index (NDVI) are variables highly correlated with, and putatively drivers of, the native distribution of A. thaliana. 83 5. Discussion 5.3. Intraspecific genetic variation and its importance in modeling the distribution and environmental suitability of A. thaliana across its native range. In this thesis, I stress the importance that intraspecific genetic variation plays in determining the distributional boundaries and ecological limits of A. thaliana across its native area. I quantified the differences in both the environmental space and suitability areas of the nine genetic clusters detected across A. thaliana’s native distribution (Chapter III). My results reinforce the idea of niche divergence among intraspecific units within the species’ native distribution, and hence, emphasize the importance of considering intraspecific variation for estimating the potential distribution of A. thaliana. Overall, the results of this chapter highlight the importance of grouping populations into genetically-informed ecologically meaningful units, which will eventually improve our suitability models of genetically structured organisms (Gotelli and Stanton-Geddes,2015;Marcer et al.,2016; Pauls et al.,2013). Previous studies provided evidence of local adaptation in A. thaliana populations, supporting the notion that environmental differences observed among genetic clusters may be the result from different environmental pressures and demographic behavior. In particular, Marcer et al. (2016) used the four genetic clusters of A. thaliana detected in the Iberian peninsula to explore their distinctive sets of drivers governing their distribution. In the same line, I quantify the divergence between the nine genomic clusters proposed by the 1001 genomes project (The 1001 Genome Consortium,2016). However, the existing spatial bias in the number of occurrences across the species’ distribution and between genetic clusters is a matter of concern for conducting such analyses. Therefore, as suggested previously, I stress the need for filling these gaps in spatial and genetic data availability if we want to improve the estimation of the environmental suitability across the species’ range. Moreover, it may be worth placing a special emphasis in surveying those areas known to represent transitions between genetic clusters. For this purpose, this thesis represents a contribution to the identification of those areas where further sampling efforts should be placed. Environmental variables influence each genetic cluster differently, suggesting the effect of different selection forces acting on A. thaliana populations across its wide distribution range. In general, temperature and precipitation are known to determine the environmental suitability and distributional boundaries of A. thaliana (Hoffmann,2002; 84 5.3. Intraspecific genetic variation for modeling distribution and environmental suitability Zou et al.,2017). Although this study reinforces the significant role of temperature in the distribution of A. thaliana, it also highlights the importance of soil pH and vegetation index (e.g. NDVI or EVI), which is often overlooked in environmental modeling analyses (Chapter I,II and III). Accurate calculations of species suitability require consideration of multiple factors. One of these crucial factors is the ecological community that hosts the species in consideration. This may include aspects such as canopy type, biome, species composition, and community stability over time. As I show in Chapter I, the relict A. thaliana accessions are highly correlated with the historical stability of vegetation over millennia. This thesis suggests that the observed niche divergence between genetic clusters of A. thaliana likely emerges from the environmentally and geographically diverse areas where the species occur. Thus, any attempt to investigate niche divergence and to disentangle the distribution boundaries of widely distributed organisms, while preventing arbitrary spatial correlations between populations pertaining to different evolving units, should first conduct preliminary molecular analysis to reveal the hidden genetic structure of the species. In addition, the use of environmental suitability modeling provides certain advantages over a multivariate environmental approach when conducting overlap analysis, as it simplifies the environmental space and incorporates the geographical space. This approach allows the assessment of the relationship between intraspecific variation and the environmental and geographical spaces at the same time. Furthermore, this thesis also provides evidence of a relationship between environmental suitability and the species’ genetic cluster structure. However, I found no significant relationship between environment and genetic distance in the environmental multivariate space. This result highlights the premise that suitability models may prove adequate to embody both environmental and geographical factors when comparing niches. Despite the fact that correlations between genetic and environmental overlaps may be weak, this thesis sets the groundwork for determining meaningful spatial modeling units by considering intraspecific variation. The use of matrix distance analysis enables the comparison of multivariate spaces of different nature, (e.g. genetic, environmental and phenotypic), which allows the analysis of the correlation of complex factors or multiple variables at the same time. It represents a promising and powerful tool for conducting muti-trait comparisons when examining complex relationships. This thesis emphasizes the importance of carrying out systematic transplant experiments across the entire distributional range of a species in order to compare phenotypic, genetic and environmental 85 5. Discussion aspects and understand the intricate relationships within a species. 5.4. A way forward Overall, this thesis contributes to our knowledge of the ecological and genetic factors that shape the evolutionary history of the relict A. thaliana lineage and provides a solid foundation for future investigations in this respect. Moving forward, collaboration among research groups, particularly involving countries with higher socio-economic instability, is of utmost importance in order to effectively address the challenges associated with quantifying phenotypic differentiation in widely distributed organisms. In particular, northern and central European countries have invested great efforts in intensive sampling in their territories. In contrast, little is known about populations in eastern European and Asian countries, which are essential to determine the distributional limits of the species and to uncover the whole intraspecific structure in order to ultimately reconstruct the evolutionary history of A. thaliana. In addition, further research should be encouraged to explore the uncharted territories where relict A. thaliana accessions have been identified, such as African, Middle Eastern, and Asian regions. These areas are essential to reconstruct both the provenance and the subsequent expansions of the species and to understand its phenotypic and genetic divergence since the early history of the species. Furthermore, I expect that genetic and phenotypic variation in these areas will be equal or even greater than the hotspots of intraspecific variation like the Iberian Peninsula, which to date is the region with the largest environmental, phenotypic and genomic variation detected so far. By expanding the sample size and conducting more comprehensive genetic studies, we will be able to improve our understanding of both functional genetic variation (i.e., polymorphisms in genes affecting its function, such as stop codons or non-synonym substitutions) within the relict lineage and intraspecific variation of A. thaliana. These efforts will not only contribute to the advancements in A. thaliana eco-evolutionary research, but will also aid in the development of new tools to identify the overall range of intraspecific variation and reconstruct the species’ historical development. This is particularly important since these methods (species distribution models, matrix models, and genetic cluster analysis) typically require large datasets to obtain accurate results and draw reliable conclusions. 86 5.4. A way forward To conclude, I would like to emphasize that the rapid advancement in the study of phenotypic variation, historical factors, and the open access to large environmental datasets, are fundamental contributors to the rapid development of consolidated biological disciplines such as biogeography and evolutionary ecology. I hope that these advancements will help fill data gaps and further demonstrate the potential of these tools to guide future research on the ecology and evolution of A. thaliana and potentially on other widely distributed organisms too. 87 Bibliography Barber, R. A., Ball, S. G., Morris, R. K., and Gilbert, F. (2022). Target-group backgrounds prove effective at correcting sampling bias in maxent models. Diversity and Distributions, 28(1):128–141. https://doi.org/10.1111/ddi.13442. Barto´n, K. (2020). MuMIn: Multi-Model Inference. R package version 1.43.17. https: //CRAN.R-project.org/package=MuMIn. Beaulieu, J.-L. D., Miras, Y., Andrieu-Ponel, V., and Guiter, F. (2005). Vegetation dynamics in north-western mediterranean regions: Instability of the mediterranean bioclimate. Plant Biosystems - An International Journal Dealing with all Aspects of Plant Biology, 139(2):114–126. https://doi.org/10.1080/11263500500197858. Beck, J. B., Schmuths, H., and Schaal, B. A. (2008). Native range genetic variation in Arabidopsis thaliana is strongly geographically structured and reflects Pleistocene glacial dynamics. Molecular Ecology, 17(3):902–915. https://doi.org/10.1111/j. 1365-294X.2007.03615.x. Becker, U., Colling, G., Dostal, P., Jakobsson, A., and Matthies, D. (2006). Local adaptation in the monocarpic perennial Carlina vulgaris at different spatial scales across Europe. Oecologia, 150:506–518. https://doi.org/10.1007/s00442-006-0534-9. Benito-Garz´on, M., Dios, R. S. D., and Ollero, H. S. (2007). Predictive modelling of tree species distributions on the iberian peninsula during the last glacial maximum and mid-holocene. Ecography, 30(1):120–134. https://doi.org/10.1111/j.0906-7590. 2007.04813.x. Benito-Garz´on, M., Robson, T. M., and Hampe, A. (2019). ∆TraitSDMs: species distribution models that account for local adaptation and phenotypic plasticity. New Phytologist, 222(4):1757–1765. https://doi.org/10.1111/nph.15716. Blanquart, F., Gandon, S., and Nuismer, S. (2012). The effects of migration and drift on local adaptation to a heterogeneous environment. Journal of evolutionary biology, 25(7):1351–1363. https://doi.org/10.1111/j.1420-9101.2012.02524.x. Blockley, S. P., Lane, C. S., Hardiman, M., Rasmussen, S. O., Seierstad, I. K., Steffensen, J. P., Svensson, A., Lotter, A. F., Turney, C. S., and Bronk Ramsey, C. (2012). Synchronisation of palaeoenvironmental records over the last 60,000 years, and an extended intimate1 event stratigraphy to 48,000 b2k. Quaternary Science Reviews, 36:2–10. https://doi.org/10.1016/j.quascirev.2011.09.017. Bobo-Pinilla, J., Pe˜nas de Giles, J., L´opez-Gonz´alez, N., Mediavilla, S., and Mart´ınez- Ortega, M. M. (2018). Phylogeography of an endangered disjunct herb: long-distance 94 Bibliography dispersal, refugia and colonization routes. AoB PLANTS, 10(5). https://doi.org/ 10.1093/aobpla/ply047. Bolnick, D. I., Amarasekare, P., Ara´ujo, M. S., B¨urger, R., Levine, J. M., Novak, M., Rudolf, V. H., Schreiber, S. J., Urban, M. C., and Vasseur, D. A. (2011). Why intraspecific trait variation matters in community ecology. Trends in Ecology & Evolution, 26(4):183–192. https://doi.org/10.1016/J.TREE.2011.01.009. Borevitz, J. O. and Nordborg, M. (2003). The Impact of Genomics on the Study of Natural Variation in Arabidopsis.Plant Physiology, 132(2):718–725. https://doi. org/10.1104/pp.103.023549. Brachi, B., Villoutreix, R., Faure, N., Hautek`eete, N., Piquot, Y., Pauwels, M., Roby, D., Cuguen, J., Bergelson, J., and Roux, F. (2013). Investigation of the geographical scale of adaptive phenological variation and its underlying genetics in Arabidopsis thaliana. Molecular Ecology, 22(16):4222–4240. https://doi.org/10.1111/mec.12396. Brennan, A. C., M´endez-Vigo, B., Haddioui, A., Mart´ınez-Zapater, J. M., Pic´o, F. X., and Alonso-Blanco, C. (2014). The genetic structure of Arabidopsis thaliana in the south-western Mediterranean range reveals a shared history between North Africa and southern Europe. BMC Plant Biology, 14(1):17. https://doi.org/10.1186/ 1471-2229-14-17. Broennimann, O., Fitzpatrick, M. C., Pearman, P. B., Petitpierre, B., Pellissier, L., Yoccoz, N. G., Thuiller, W., Fortin, M.-J., Randin, C., Zimmermann, N. E., et al. (2012). Measuring ecological niche overlap from occurrence and spatial environmental data. Global ecology and biogeography, 21(4):481–497. https://doi.org/10.1111/j. 1466-8238.2011.00698.x. Burjachs, F., Jones, S. E., Giralt, S., and Fern´andez-L´opez de Pablo, J. (2016). Lateglacial to Early Holocene recursive aridity events in the SE Mediterranean Iberian Peninsula: The Salines playa lake case study. Quaternary International, 403:187–200. https: //doi.org/10.1016/j.quaint.2015.10.117. Campbell, J. F. E., Fletcher, W. J., Joannin, S., Hughes, P. D., Rhanem, M., and Zielhofer, C. (2017). Environmental drivers of Holocene dorest development in the Middle Atlas, Morocco. Frontiers in Ecology and Evolution, 5(September):1–22. https://doi.org/10.3389/fevo.2017.00113. Carri´on, J. S., Fern´andez, S., Gonz´alez-Samp´eriz, P., Gil-Romera, G., Badal, E., Carri´on- Marco, Y., L´opez-Merino, L., L´opez-S´aez, J. A., Fierro, E., and Burjachs, F. (2010). 95 Bibliography Expected trends and surprises in the Lateglacial and Holocene vegetation history of the Iberian Peninsula and Balearic Islands. Review of Palaeobotany and Palynology, 162(3):458–475. https://doi.org/10.1016/j.revpalbo.2009.12.007. Castilla, A. R., M´endez-Vigo, B., Marcer, A., Mart´ınez-Minaya, J., Conesa, D., Pic´o, F. X., and Alonso-Blanco, C. (2020). Ecological, genetic and evolutionary drivers of regional genetic differentiation in Arabidopsis thaliana.BMC Evolutionary Biology, 20(1):71. https://doi.org/10.1186/s12862-020-01635-2. Chang, C. C., Chow, C. C., Tellier, L. C., Vattikuti, S., Purcell, S. M., and Lee, J. J. (2015). Second-generation PLINK: rising to the challenge of larger and richer datasets. Gigascience, 4(1):s13742–015. https://doi.org/10.1186/s13742-015-0047-8. Cheddadi, R., Henrot, A.-J., Fran¸cois, L., Boyer, F., Bush, M., Carr´e, M., Coissac, E., De Oliveira, P. E., Ficetola, F., Hambuckers, A., Huang, K., L´ezine, A.-M., Nourelbait, M., Rhoujjati, A., Taberlet, P., Sarmiento, F., Abel-Schaad, D., Alba-S´anchez, F., and Zheng, Z. (2017). Microrefugia, climate change, and conservation of Cedrus atlantica in the Rif mountains, Morocco. Frontiers in Ecology and Evolution, 5:114. https: //doi.org/10.3389/fevo.2017.00114. Clarke, R. T., Rothery, P., and Raybould, A. F. (2002). Confidence limits for regression relationships between distance matrices: Estimating gene flow with distance. Journal of Agricultural, Biological, and Environmental Statistics, 7(3):361–372. https://doi. org/10.1198/108571102320. Collart, F., Heden¨as, L., Broennimann, O., Guisan, A., and Vanderpoorten, A. (2021). Intraspecific differentiation: Implications for niche and distribution modelling. Journal of Biogeography, 48(2):415–426. https://doi.org/10.1111/jbi.14009. Comes, H. P. and Kadereit, J. W. (1998). The effect of quaternary climatic changes on plant distribution and evolution. Trends in Plant Science, 3(11):432–438. https: //doi.org/10.1016/S1360-1385(98)01327-2. Compton, J. S. (2011). Pleistocene sea-level fluctuations and human evolution on the southern coastal plain of South Africa. Quaternary Science Reviews, 30(5):506–527. https://doi.org/10.1016/j.quascirev.2010.12.012. Cowling, R. M., Rundel, P. W., Lamont, B. B., Arroyo, M. K., and Arianoutsou, M. (1996). Plant diversity in mediterranean-climate regions. Trends in Ecology and Evolution, 11(9):362–366. https://doi.org/10.1016/0169-5347(96)10044-6. 96 Bibliography Davila, N. H., Frago, E., Thoen, M. P. M., Kloth, K. J., Becker, F. F. M., van Loon, J. J. A., Gort, G., Keurentjes, J. J. B., van Heerwaarden, J., and Dicke, M. (2017). Natural variation in life history strategy of Arabidopsis thaliana determines stress responses to drought and insects of different feeding guilds. Molecular Ecology, 26(11):2959–2977. https://doi.org/10.1111/mec.14100. de Montaigu, A. and Coupland, G. (2017). The timing of gigantea expression during day/night cycles varies with the geographical origin of Arabidopsis accessions. Plant Signaling & Behavior, 12(7):e1342026. https://doi.org/10.1080/15592324.2017. 1342026. De Noblet-Ducoudr´e, N., Claussen, M., and Prentice, C. (2000). Mid-holocene greening of the sahara: first results of the gaim 6000 year bp experiment with two asynchronously coupled atmosphere/biome models. Climate Dynamics, 16(9):643–659. https://doi. org/10.1007/s003820000074. DeLeo, V. L., Menge, D. N., Hanks, E. M., Juenger, T. E., and Lasky, J. R. (2020). Effects of two centuries of global environmental variation on phenology and physiology of Arabidopsis thaliana.Global Change Biology, 26(2):523–538. https://doi.org/10. 1111/gcb.14880. Des Roches, S., Post, D. M., Turley, N. E., Bailey, J. K., Hendry, A. P., Kinnison, M. T., Schweitzer, J. A., and Palkovacs, E. P. (2018). The ecological importance of intraspecific variation. Nature ecology & evolution, 2(1):57–64. https://doi.org/10. 1038/s41559-017-0402-5. Diamond, J. M. (1969). Avifaunal equilibria and species turnover rates on the channel islands of california. Proceedings of the National Academy of Sciences, 64(1):57–63. https://doi.org/10.1073/pnas.64.1.57. DiMiceli, C., Carroll, M., Sohlberg, R., Kim, D., Kelly, M., and Townshend, J. (2015). Mod44b modis/terra vegetation continuous fields yearly l3 global 250m sin grid v006 [data set]. https://doi.org/10.5067/MODIS/MOD44B.006. Donoghue, M. J. and Edwards, E. J. (2014). Biome shifts and niche evolution in plants. Annual Review of Ecology, Evolution, and Systematics, 45(1):547–572. https://doi. org/10.1146/annurev-ecolsys-120213-091905. Donohue, K. (2005). Niche construction through phenological plasticity: life history dynamics and ecological consequences. New Phytologist, 166(1):83–92. https://doi. org/10.1111/j.1469-8137.2005.01357.x. 97 Bibliography Durvasula, A., Fulgione, A., Gutaker, R. M., Alacakaptan, S. I., Flood, P. J., Neto, C., Tsuchimatsu, T., Burbano, H. A., Pic´o, F. X., Alonso-Blanco, C., and Hancock, A. M. (2017). African genomes illuminate the early history and transition to selfing in Arabidopsis thaliana.Proceedings of the National Academy of Sciences, 114(20):5213– 5218. https://doi.org/10.1073/pnas.1616736114. Dutilleul, P., Clifford, P., Richardson, S., and Hemon, D. (1993). Modifying the t test for assessing the correlation between two spatial processes. Biometrics, 49(1):305–314. https://doi.org/10.2307/2532625. El-Din El-Assal, S., Alonso-Blanco, C., Peeters, A. J., Raz, V., and Koornneef, M. (2001). A QTL for flowering time in arabidopsis reveals a novel allele of CRY2. Nature genetics, 29(4):435–440. https://doi.org/10.1038/ng767. Elith, J. and Leathwick, J. R. (2009). Species distribution models: Ecological explanation and prediction across space and time. Annual Review of Ecology, Evolution, and Systematics, 40(1):677–697. https://doi.org/10.1146/annurev.ecolsys.110308. 120159. Elith, J., Phillips, S. J., Hastie, T., Dud´ık, M., Chee, Y. E., and Yates, C. J. (2011). A statistical explanation of maxent for ecologists. Diversity and Distributions, 17(1):43– 57. https://doi.org/10.1111/j.1472-4642.2010.00725.x. Estarague, A., Vasseur, F., Sartori, K., Bastias, C. C., Cornet, D., Rouan, L., Beurier, G., Exposito-Alonso, M., Herbette, S., Bresson, J., et al. (2022). Into the range: a latitudinal gradient or a center-margins differentiation of ecological strategies in Arabidopsis thaliana?Annals of Botany, 129(3):343–356. https://doi.org/10.1093/ aob/mcab149. Exposito-Alonso, M., Brennan, A. C., Alonso-Blanco, C., and Pic´o, F. X. (2018). Spatiotemporal variation in fitness responses to contrasting environments in Arabidopsis thaliana.Evolution, 72(8):1570–1586. https://doi.org/10.1111/evo.13508. Exposito-Alonso, M., G´omez Rodr´ıguez, R., Barrag´an, C., Capovilla, G., Chae, E., Devos, J., Dogan, E. S., Friedemann, C., Gross, C., Lang, P., Lundberg, D., Middendorf, V., Kageyama, J., Karasov, T., Kersten, S., Petersen, S., Rabbani, L., Regalado, J., Reinelt, L., Rowan, B., Seymour, D. K., Symeonidi, E., Schwab, R., Tran, D. N., Venkataramani, K., de Weyer, A. L. V., Vasseur, F., Wang, G., Wedeg¨artner, R., Weiss, F., Wu, R., Xi, W., Zaidem, M., Zhu, W., Garc´ıa-Arenal, F., Burbano, H. A., Bossdorf, O., Weigel, D., Burbano, H. A., Bossdorf, O., Nielsen, R., and Weigel, D. (2019). Natural selection on 98 Bibliography the Arabidopsis thaliana genome in present and future climates. Nature, 573(7772):126– 129. https://doi.org/10.1038/s41586-019-1520-9. FAO, IIASA, ISRIC, ISSCAS, JRC (2012). Harmonized world soil database (version 1.2) [data set]. https://doi.org/10.5067/MODIS/MCD12Q1.006. Fern´andez-Mazuecos, M. and Vargas, P. (2013). Congruence between distribution modelling and phylogeographical analyses reveals quaternary survival of a toadflax species (Linaria elegans) in oceanic climate areas of a mountain ring range. New Phytologist, 198(4):1274–1289. https://doi.org/10.1111/nph.12220. Fick, S. E. and Hijmans, R. J. (2017). Worldclim 2: new 1-km spatial resolution climate surfaces for global land areas. International Journal of Climatology, 37(12):4302–4315. https://doi.org/10.1002/joc.5086. Fletcher, W. J., Boski, T., and Moura, D. (2007). Palynological evidence for environmental and climatic change in the lower guadiana valley, portugal, during the last 13 000 years. The Holocene, 17(4):481–494. https://doi.org/10.1177/0959683607077027. Forsman, A. (2014). Effects of genotypic and phenotypic variation on establishment are important for conservation, invasion, and infection biology. Proceedings of the National Academy of Sciences of the United States of America, 111(1):302–307. https: //doi.org/10.1111/jbi.12562. Fournier-Level, A., Korte, A., Cooper, M. D., Nordborg, M., Schmitt, J., and a. M. Wilczek (2011). A map of local adaptation in Arabidopsis thaliana.Science, 334(6052):86–89. https://doi.org/10.1126/science.1209271. Friedl, M. and Sulla-Menashe, D. (2015). Mcd12q1 modis/terra+aqua land cover type yearly l3 global 500m sin grid v006 [data set]. https://doi.org/MODIS/MCD12Q1.006. Fulgione, A. and Hancock, A. M. (2018). Archaic lineages broaden our view on the history of Arabidopsis thaliana.New Phytologist, 334(6052):86–89. https://doi.org/ 10.1111/nph.15244. Fulgione, A., Koornneef, M., Roux, F., Hermisson, J., and Hancock, A. M. (2017). Madeiran Arabidopsis thaliana reveals ancient long-range colonization and clarifies demography in Eurasia. Molecular Biology and Evolution, 35(3):564–574. https: //doi.org/10.1093/molbev/msx300. Garant, D., Kruuk, L. E., Wilkin, T. A., McCleery, R. H., and Sheldon, B. C. (2005). Evolution driven by differential dispersal within a wild bird population. Nature, 433(7021):60–65. https://doi.org/10.1038/nature03051. 99 Bibliography GDAL and OGR (2015). GDAL/OGR Geospatial Data Abstraction Software Library [Data set]. http://gdal.org. Gil-Romera, G., Ant´on, M., and Calleja, J. (2008). The late holocene palaeoecological sequence of serran´ıa de las villuercas (southern meseta, western spain). Vegetation History and Archaeobotany, 17:653–666. https://doi.org/10.1007/s00334-008-0146-x. Gonz´alez-Samp´eriz, P., Aranbarri, J., P´erez-Sanz, A., Gil-Romera, G., Moreno, A., Leunda, M., Sevilla-Callejo, M., Corella, J. P., Morell´on, M., Oliva, B., and Valero-Garc´es, B. (2017). Environmental and climate change in the southern Central Pyrenees since the Last Glacial Maximum: A view from the lake records. CATENA, 149:668–688. https://doi.org/10.1016/j.catena.2016.07.041. Goslee, S. C., Urban, D., and Urban, D. L. (2007). The ecodist package for dissimilaritybased analysis of ecological data journal of statistical software the ecodist package for dissimilarity-based analysis of ecological data. Article in Journal of Statistical Software. https://doi.org/10.18637/jss.v022.i07. Gotelli, N. J. and Stanton-Geddes, J. (2015). Climate change, genetic markers and species distribution modelling. Journal of Biogeography, 42(9):1577–1585. https://doi.org/ 10.1111/jbi.12562. G´omez, A. and Lunt, D. H. (2007). Refugia within refugia: Patterns of phylogeographic concordance in the Iberian Peninsula, pages 155–188. Springer Netherlands, Dordrecht. https://doi.org/10.1007/1-4020-4904-8_5. G´omez, R., M´endez-Vigo, B., Marcer, A., Alonso-Blanco, C., and Pic´o, F. X. (2018). Quantifying temporal change in plant population attributes: insights from a resurrection approach. AoB PLANTS, 10(5). https://doi.org/10.1093/aobpla/ply063. Halvorsen, R., Mazzoni, S., Bryn, A., and Bakkestuen, V. (2015). Opportunities for improved distribution modelling practice via a strict maximum likelihood interpretation of maxent. Ecography, 38(2):172–183. https://doi.org/10.1111/ecog.00565. Hancock, A. M., Brachi, B., Faure, N., Horton, M. W., Jarymowycz, L. B., Sperone, F. G., Toomajian, C., Roux, F., and Bergelson, J. (2011). Adaptation to climate across the Arabidopsis thaliana genome. Science, 334(6052):83–86. https://doi.org/10. 1126/SCIENCE.1209244. Hannah, L. (2012). Saving a Million Species, pages 389–394. Island Press/Center for Resource Economics, Washington, DC. https://doi.org/10.5822/ 978-1-61091-182-5_20. 100 Bibliography Herrando-Moraira, S., Nualart, N., Herrando-Moraira, A., Chung, M. Y., Chung, M. G., and L´opez-Pujol, J. (2019). Climatic niche characteristics of native and invasive lilium lancifolium. Scientific Reports, 9(1):1–16. https://doi.org/10.1038/ng.695. Herrera, C. M. (1992). Historical effects and sorting processes as explanations for contemporary ecological patterns: Character syndromes in Mediterranean woody plants. The American Naturalist, 140(3):421–446. https://doi.org/10.1086/285420. Hewitt, G. M. (1999). Post-glacial re-colonization of european biota. Biological Journal of the Linnean Society, 68(1):87–112. https://doi.org/10.1006/bijl.1999.0332. Hijmans, R. J. (2012). Cross-validation of species distribution models: removing spatial sorting bias and calibration with a null model. Ecology, 93(3):679–688. https://doi. org/10.1038/ng.695. Hijmans, R. J. (2019). geosphere: Spherical Trigonometry. R package version 1.5-10. https://CRAN.R-project.org/package=geosphere. Hijmans, R. J. (2020). terra: Spatial Data Analysis. R package version 0.9-8. https: //CRAN.R-project.org/package=terra. Hijmans, R. J., Phillips, S., Leathwick, J., and Elith, J. (2017). Package dismo: Species Distribution Modeling.https://CRAN.R-project.org/package=dismo. Hoffmann, M. H. (2002). Biogeography of Arabidopsis thaliana (l.) heynh. (Brassicaceae). Journal of Biogeography, 29(1):125–134. https://doi.org/10.1046/j.1365-2699. 2002.00647.x. Honjo, M. N. and Kudoh, H. (2019). Arabidopsis halleri: a perennial model system for studying population differentiation and local adaptation. AoB PLANTS, 11(6):plz076. https://doi.org/10.1093/aobpla/plz076. Horton, M. W., Hancock, A. M., Huang, Y. S., Toomajian, C., Atwell, S., Auton, A., Muliyati, N. W., Platt, A., Sperone, F. G., Vilhj´almsson, B. J., Nordborg, M., Borevitz, J. O., and Bergelson, J. (2012). Genome-wide patterns of genetic variation in worldwide Arabidopsis thaliana accessions from the RegMap panel. Nature Genetics, 44(2):212– 216. https://doi.org/10.1038/ng.1042. Hsu, C. W., Lo, C. Y., and Lee, C. R. (2019). On the postglacial spread of human commensal Arabidopsis thaliana: journey to the East. New Phytologist, 222(3):1447– 1457. https://doi.org/10.1111/nph.15682. 101 Bibliography Hu, Z. M., Zhang, Q. S., Zhang, J., Kass, J. M., Mammola, S., Fresia, P., Draisma, S. G., Assis, J., Jueterbock, A., Yokota, M., and Zhang, Z. (2021). Intraspecific genetic variation matters when predicting seagrass distribution under climate change. Molecular Ecology, 30(15):3840–3855. https://doi.org/10.1111/MEC.15996. Huang, X., Wei, X., Sang, T., Zhao, Q., Feng, Q., Zhao, Y., Li, C., Zhu, C., Lu, T., Zhang, Z., et al. (2010). Genome-wide association studies of 14 agronomic traits in rice landraces. Nature genetics, 42(11):961–967. https://doi.org/10.1038/ng.695. Hughes, P. D., Woodward, J. C., and Gibbard, P. L. (2006). Quaternary glacial history of the mediterranean mountains. Progress in Physical Geography: Earth and Environment, 30(3):334–364. https://doi.org/10.1191/0309133306pp481ra. IBM Corp (2017). Ibm spss statistics for windows, version 25.0. James, A. B., Sullivan, S., and Nimmo, H. G. (2018). Global spatial analysis of Arabidopsis natural variants implicates 5’utr splicing of late elongated hypocotyl in responses to temperature. Plant Cell and Environment, 41(7):1524–1538. https://doi.org/10. 1111/pce.13188. Johnson, L. C., Galliart, M. B., Alsdurf, J. D., Maricle, B. R., Baer, S. G., Bello, N. M., Gibson, D. J., and Smith, A. B. (2022). Reciprocal transplant gardens as gold standard to detect local adaptation in grassland species: New opportunities moving into the 21st century. Journal of Ecology, 110(5):1054–1071. https: //doi.org/10.1111/1365-2745.13695. Jones, S. E., Burjachs, F., Ferrer-Garc´ıa, C., Giralt, S., Schulte, L., and Fern´andez-L´opez de Pablo, J. (2018). A multi-proxy approach to understanding complex responses of salt-lake catchments to climate variability and human pressure: A Late Quaternary case study from south-eastern, Spain. Quaternary Science Reviews, 184:201–223. https: //doi.org/10.1016/j.quascirev.2017.12.015. Joshi, J., Schmid, B., Caldeira, M. C., Dimitrakopoulos, P. G., Good, J., Harris, R., Hector, A., Huss-Danell, K., Jumpponen, A., Minns, A., Mulder, C. P., Pereira, J. S., Prinz, A., Scherer-Lorenzen, M., Siamantziouras, A. S., Terry, A. C., Troumbis, A. Y., and Lawton, J. H. (2001). Local adaptation enhances performance of common plant species. Ecology Letters, 4(6):536–544. https://doi.org/10.1046/J.1461-0248. 2001.00262.X. JRC-ESDAC (2015). Soil geographical database for eurasia and the mediterranean (version 4.0) [data set]. https://esdac.jrc.ec.europa.eu/resource-type/ european-soil-database-soil-properties. 102 Bibliography Jump, A. S., Marchant, R., and Pe˜nuelas, J. (2009). Environmental change and the option value of genetic diversity. Trends in plant science, 14(1):51–58. https://doi.org/10. 1016/j.tplants.2008.10.002. Kahl, S. M., Kappel, C., Joshi, J., and Lenhard, M. (2021). Phylogeography of a widely distributed plant species reveals cryptic genetic lineages with parallel phenotypic responses to warming and drought conditions. Ecology and Evolution, 11(20):13986– 14002. https://doi.org/10.1002/ece3.8103. Katz, E., Li, J.-J., Jaegle, B., Ashkenazy, H., Abrahams, S. R., Bagaza, C., Holden, S., Pires, C. J., Angelovici, R., and Kliebenstein, D. J. (2021). Genetic variation, environment and demography intersect to shape Arabidopsis defense metabolite variation across europe. eLife, 10. https://doi.org/10.7554/ELIFE.67784. Kawakatsu, T., shan Carol Huang, S., Jupe, F., Sasaki, E., Schmitz, R. J., Urich, M. A., Castanon, R., Nery, J. R., Barragan, C., He, Y., Chen, H., Dubin, M., Lee, C. R., Wang, C., Bemm, F., Becker, C., O’Neil, R., O’Malley, R. C., Quarless, D. X., Alonso- Blanco, C., Andrade, J., Bemm, F., Bergelson, J., Borgwardt, K., Chae, E., Dezwaan, T., Ding, W., Ecker, J. R., Exp´osito-Alonso, M., Farlow, A., Fitz, J., Gan, X., Grimm, D. G., Hancock, A., Henz, S. R., Holm, S., Horton, M., Jarsulic, M., Kerstetter, R. A., Korte, A., Korte, P., Lanz, C., Lee, C. R., Meng, D., Michael, T. P., Mott, R., Muliyati, N. W., N¨agele, T., Nagler, M., Nizhynska, V., Nordborg, M., Novikova, P., Pic´o, F. X., Platzer, A., Rabanal, F. A., Rodriguez, A., Rowan, B. A., Salom´e, P. A., Schmid, K., Schmitz, R. J., ¨ Umit Seren, Sperone, F. G., Sudkamp, M., Svardal, H., Tanzer, M. M., Todd, D., Volchenboum, S. L., Wang, C., Wang, G., Wang, X., Weckwerth, W., Weigel, D., Zhou, X., Schork, N. J., Weigel, D., Nordborg, M., and Ecker, J. R. (2016). Epigenomic diversity in a global collection of Arabidopsis thaliana accessions. Cell, 166(2):492–505. https://doi.org/10.1016/j.cell.2016.06.044. Kawecki, T. J. and Ebert, D. (2004). Conceptual issues in local adaptation. Ecology Letters, 7(12):1225–1241. https://doi.org/10.1111/J.1461-0248.2004.00684.X. Kerwin, R. E., Feusier, J., Muok, A., Lin, C., Larson, B., Copeland, D., Corwin, J. A., Rubin, M. J., Francisco, M., Li, B., et al. (2017). Epistasis×environment interactions among Arabidopsis thaliana glucosinolate genes impact complex traits and fitness in the field. New Phytologist, 215(3):1249–1263. https://doi.org/10.1111/nph.14646. Koch, M. A., Wernisch, M., Schmickl, R., Url, S., and Koch, A. (2008). Arabidopsis thaliana’s wild relatives: An updated overview on systematics, taxonomy and evolution. Taxon, 57(3):933–943. https://doi.org/10.1002/tax.573021. 103 Bibliography J. D., Frommer, W. B., Grotewold, E., Meyerowitz, E., Nemhauser, J., Nordborg, M., Pikaard, C., Shanklin, J., Somerville, C., Stitt, M., Torii, K. U., Waese, J., Wagner, D., and Mccourt, P. (2016). 50 years of Arabidopsis research: Highlights and future directions. New Phytologist, 209(3):921–944. https://doi.org/10.1111/nph.13687. Qu´ezel, P. and Barbero, M. (1989). Les formations `a gen´evriers rampants du djurdjura (alg´erie). leur signification ecologique, dynamique et syntaxinomique dans une approche globale des cedraies kabyles. Lazaroa, 11:85–100. R Core Team (2017). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/. Rangel, T. F., Diniz-Filho, J. A. F., and Bini, L. M. (2010). Sam: a comprehensive application for spatial analysis in macroecology. Ecography, 33(1):46–50. https:// doi.org/10.1111/j.1600-0587.2009.06299.x. Rehfeldt, G. E., Tchebakova, N. M., Parfenova, Y. I., Wykoff, W. R., Kuzmina, N. A., and Milyutin, L. I. (2002). Intraspecific responses to climate in Pinus sylvestris. Global Change Biology, 8(9):912–929. https://doi.org/10.1046/j.1365-2486. 2002.00516.x. Ripley, B., Venables, B., Bates, D. M., Hornik, K., Gebhardt, A., Firth, D., and Ripley, M. B. (2013). Package ‘mass’. Cran r, 538:113–120. Roches, S. D., Pendleton, L. H., Shapiro, B., and Palkovacs, E. P. (2021). Conserving intraspecific variation for nature’s contributions to people. Nature Ecology and Evolution 2021 5:5, 5(5):574–582. https://doi.org/10.1038/s41559-021-01403-5. Rodr´ıguez-Rodr´ıguez, E. J., Beltr´an, J. F., Tejedo, M., Nicieza, A. G., Llusia, D., M´arquez, R., and Arag´on, P. (2020). Niche models at inter- and intraspecific levels reveal hierarchical niche differentiation in midwife toads. Scientific Reports, 10:10942. https://doi.org/10.1038/s41598-020-67992-6. Rodr´ıguez-S´anchez, F., Hampe, A., Jordano, P., and Arroyo, J. (2010). Past tree range dynamics in the iberian peninsula inferred through phylogeography and palaeodistribution modelling: A review. Review of Palaeobotany and Palynology, 162(3):507–521. https://doi.org/10.1016/j.revpalbo.2010.03.008. Rousset, F. (1997). Genetic differentiation and estimation of gene flow from F-Statistics under isolation by distance. Genetics, 145(4):1219–1228. https://doi.org/10.1093/ genetics/145.4.1219. 110 Bibliography Santos, L. and Go˜ni, M. F. S. (2003). Lateglacial and Holocene environmental changes in Portuguese coastal lagoons 3: vegetation history of the Santo Andre coastal area. The Holocene, 13(3):459–464. https://doi.org/10.1191/0959683603hl638rp. Savolainen, O. and Lascoux, M. (2016). Geography matters for Arabidopsis.Nature 2016 537:7620, 537(7620):314–315. https://doi.org/10.1038/nature19466. Savolainen, O., Lascoux, M., and Meril¨a, J. (2013). Ecological genomics of local adaptation. Nature Publishing Group, 14:807. https://doi.org/10.1038/nrg3522. Schoener, T. W. (1968). The Anolis lizards of Bimini: Resource partitioning in a complex fauna. Ecology, 49(4):704–726. https://doi.org/10.2307/1935534. Schreiber, S. J., B¨urger, R., and Bolnick, D. I. (2011). The community effects of phenotypic and genetic variation within a predator population. Ecology, 92(8):1582–1593. https://doi.org/10.1890/10-2071.1. Sexton, J. P., Hangartner, S. B., and Hoffmann, A. A. (2014). Genetic isolation by environment or distance: Which pattern of gene flow is most common? Evolution, 68(1):1–15. https://doi.org/10.1111/EVO.12258. Shindo, C., Bernasconi, G., and Hardtke, C. S. (2007). Natural genetic variation in Arabidopsis: Tools, traits and prospects for evolutionary ecology. Annals of Botany, 99(6):1043. https://doi.org/10.1093/AOB/MCL281. Sides, C. B., Enquist, B. J., Ebersole, J. J., Smith, M. N., Henderson, A. N., and Sloat, L. L. (2014). Revisiting darwin’s hypothesis: Does greater intraspecific variability increase species’ ecological breadth? American journal of botany, 101(1):56–62. https: //doi.org/10.3732/ajb.1300284. Slatkin, M. (1993). Isolation by distance in equilibrium and non-equilibrium populations. Evolution, 47(1):264–279. https://doi.org/10.1111/j.1558-5646.1993.tb01215. x. Smith, A. B., Godsoe, W., Rodr´ıguez-S´anchez, F., Wang, H.-H., and Warren, D. (2019). Niche estimation above and below the species level. Trends in Ecology & Evolution, 34(3):260–273. https://doi.org/10.1016/j.tree.2018.10.012. Stange, M., Barrett, R. D., and Hendry, A. P. (2020). The importance of genomic variation for biodiversity, ecosystems and people. Nature Reviews Genetics, 22(2):89–105. https: //doi.org/10.1038/s41576-020-00288-7. 111 Bibliography Sultan, S. E. (1987). Evolutionary Implications of Phenotypic Plasticity in Plants, pages 127–178. Springer US, Boston, MA. https://doi.org/10.1007/ 978-1-4615-6986-2_7. Tabas-Madrid, D., M´endez-Vigo, B., Arteaga, N., Marcer, A., Pascual-Montano, A., Weigel, D., Xavier Pic´o, F., and Alonso-Blanco, C. (2018). Genome-wide signatures of flowering adaptation to climate temperature: Regional analyses in a highly diverse native range of Arabidopsis thaliana.Plant, Cell and Environment, 41(8):1806–1820. https://doi.org/10.1111/pce.13189. Takou, M., Wieters, B., Kopriva, S., Coupland, G., Linst¨adter, A., and Meaux, J. D. (2019). Linking genes with ecological strategies in Arabidopsis thaliana.Journal of Experimental Botany, 70(4):1141–1151. https://doi.org/10.1093/JXB/ERY447. Tamura, K., Peterson, D., Peterson, N., Stecher, G., Nei, M., and Kumar, S. (2011). MEGA5: molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Molecular biology and evolution, 28(10):2731–2739. https://doi.org/10.1093/molbev/msr121. The 1001 Genome Consortium (2016). 1,135 genomes reveal the global pattern of polymorphism in Arabidopsis thaliana.Cell, 166(2):481–491. https://doi.org/10.1016/ j.cell.2016.05.063. Thioulouse, J., Dray, S., Dufour, A.-B., Siberchicot, A., Jombart, T., and Pavoine, S. (2018). Multivariate Analysis of Ecological Data with ade4. Springer. https://doi. org/10.1007/978-1-4939-8850-1. Title, P. O. and Bemmels, J. B. (2018). Envirem: an expanded set of bioclimatic and topographic variables increases flexibility and improves performance of ecological niche modeling. Ecography, 41(2):291–307. https://doi.org/10.1111/ecog.02880. Toledo, B., Marcer, A., M´endez-Vigo, B., Alonso-Blanco, C., and Pic´o, F. X. (2020). An ecological history of the relict genetic lineage of Arabidopsis thaliana.Environmental and Experimental Botany, 170:103800. https://doi.org/10.1016/j.envexpbot. 2019.103800. Tsuchimatsu, T., Goubet, P. M., Gallina, S., Holl, A.-C., Fobis-Loisy, I., Berg`es, H., Marande, W., Prat, E., Meng, D., Long, Q., Platzer, A., Nordborg, M., Vekemans, X., and Castric, V. (2017). Patterns of polymorphism at the self-incompatibility locus in 1,083 Arabidopsis thaliana genomes. Molecular Biology and Evolution, 34(8):1878–1889. https://doi.org/10.1093/molbev/msx122. 112 Bibliography Tuanmu, M.-N. and Jetz, W. (2015). A global, remote sensing-based characterization of terrestrial habitat heterogeneity for biodiversity and ecosystem modelling. Global Ecology and Biogeography, 24(11):1329–1339. https://doi.org/10.1111/geb.12365. Vaissi, S. and Rezaei, S. (2022). Niche divergence at intraspecific level in the Hyrcanian wood frog, Rana pseudodalmatina: A phylogenetic, climatic, and environmental survey. Frontiers in Ecology and Evolution., 10:774481. https://doi.org/10.3389/ fevo.2022.774481. Valavi, R., Elith, J., Lahoz-Monfort, J. J., and Guillera-Arroita, G. (2018). blockCV: An R package for generating spatially or environmentally separated folds for k-fold crossvalidation of species distribution models. Methods in Ecology and Evolution, 10(2):225– 232. https://doi.org/10.1111/2041-210X.13107. Vasseur, F., Exposito-Alonso, M., Ayala-Garay, O. J., Wang, G., Enquist, B. J., Vile, D., Violle, C., and Weigel, D. (2018a). Adaptive diversification of growth allometry in the plant Arabidopsis thaliana.Proceedings of the National Academy of Sciences of the United States of America, 115(13):3416–3421. https://doi.org/10.1073/pnas. 1709141115. Vasseur, F., Sartori, K., Baron, E., Fort, F., Kazakou, E., Segrestin, J., Garnier, E., Vile, D., and Violle, C. (2018b). Climate as a driver of adaptive variations in ecological strategies in Arabidopsis thaliana.Annals of Botany, 122(6):935–945. https://doi. org/10.1093/aob/mcy165. Vidigal, D. S., Marques, A. C., Willems, L. A., Buijs, G., M´endez-Vigo, B., Hilhorst, H. W., Bentsink, L., Pic´o, F. X., and Alonso-Blanco, C. (2016). Altitudinal and climatic associations of seed dormancy and flowering traits evidence adaptation of annual life cycle timing in Arabidopsis thaliana.Plant Cell and Environment, 39(8):1737–1748. https://doi.org/10.1111/pce.12734. Villellas, J., Ehrl´en, J., Crone, E. E., Cserg˝o, A. M., Garcia, M. B., Laine, A.-L., Roach, D. A., Salguero-G´omez, R., Wardle, G. M., Childs, D. Z., et al. (2021). Phenotypic plasticity masks range-wide genetic differentiation for vegetative but not reproductive traits in a short-lived plant. Ecology Letters, 24(11):2378–2393. https://doi.org/10. 1111/ele.13858. Vollering, J., Halvorsen, R., and Mazzoni, S. (2019). The MIAmaxent R package: Variable transformation and model selection for species distribution models. Ecology and Evolution, 9(21):12051–12068. https://doi.org/10.1002/ece3.5654. 113 Bibliography Wand, M. (2020). KernSmooth: Functions for Kernel Smoothing Supporting Wand and Jones (1995). R package version 2.23-18. https://CRAN.R-project.org/package= KernSmooth. Wang, I. J. and Bradburd, G. S. (2014). Isolation by environment. Molecular Ecology, 23(23):5649–5662. https://doi.org/10.1111/MEC.12938. Warren, D. L., Glor, R. E., and Turelli, M. (2008). Environmental niche equivalency versus conservatism: quantitative approaches to niche evolution. Evolution, 62(11):2868–2883. https://doi.org/10.1111/j.1558-5646.2008.00482.x. Weigel, D. (2011). Natural variation in Arabidopsis: From molecular genetics to ecological genomics. Plant Physiology, 158(1):2–22. https://doi.org/10.1104/pp.111.189845. Weir, B. S., Anderson, A. D., and Hepler, A. B. (2006). Genetic relatedness analysis: modern data and new challenges. Nature Reviews Genetics, 7(10):771–780. https: //doi.org/10.1038/nrg1960. Weiß, C., Schuenemann, V., Devos, J., Shirsekar, G., Reiter, E., Gould, B., Stinchcombe, J., Krause, J., and Burbano, H. (2016). Temporal patterns of damage and decay kinetics of DNA retrieved from plant herbarium specimens. Royal Society Open Science, 3:160239. https://doi.org/10.1098/rsos.160239. Weng, C., Hooghiemstra, H., and Duivenvoorden, J. F. (2006). Challenges in estimating past plant diversity from fossil pollen data: statistical assessment, problems, and possible solutions. Diversity and Distributions, 12(3):310–318. https: //doi.org/10.1111/j.1366-9516.2006.00230.x. Wiens, J. J., Ackerly, D. D., Allen, A. P., Anacker, B. L., Buckley, L. B., Cornell, H. V., Damschen, E. I., Jonathan Davies, T., Grytnes, J.-A., Harrison, S. P., et al. (2010). Niche conservatism as an emerging principle in ecology and conservation biology. Ecology letters, 13(10):1310–1324. https://doi.org/10.1111/j.1461-0248.2010.01515. x. Wilczek, A. M., Cooper, M. D., Korves, T. M., and Schmitt, J. (2014). Lagging adaptation to warming climate in Arabidopsis thaliana.Proceedings of the National Academy of Sciences of the United States of America, 111(22):7906–7913. https://doi.org/10. 1073/pnas.1406314111. Wilczek, A. M., Roe, J. L., Knapp, M. C., Cooper, M. D., Lopez-Gallego, C., Martin, L. J., Muir, C. D., Sim, S., Walker, A., Anderson, J., et al. (2009). Effects of genetic 114 Bibliography perturbation on seasonal life history plasticity. Science, 323(5916):930–934. https: //doi.org/10.1126/science.1165826. Willi, Y., Van Buskirk, J., Schmid, B., and Fischer, M. (2007). Genetic isolation of fragmented populations is exacerbated by drift and selection. Journal of evolutionary biology, 20(2):534–542. https://doi.org/10.1111/j.1420-9101.2006.01263.x. Wood, Z. T., Wiegardt, A. K., Barton, K. L., Clark, J. D., Homola, J. J., Olsen, B. J., King, B. L., Kovach, A. I., and Kinnison, M. T. (2021). Meta-analysis: Congruence of genomic and phenotypic differentiation across diverse natural study systems. Evolutionary applications, 14(9):2189–2205. https://doi.org/10.1111/eva.13264. Wordley, C. F., Sankaran, M., Mudappa, D., and Altringham, J. D. (2015). Landscape scale habitat suitability modelling of bats in the western ghats of india: Bats like something in their tea. Biological Conservation, 191:529–536. https://doi.org/10. 1016/J.BIOCON.2015.08.005. Yao, T., Zhang, J., Xie, M., Yuan, G., Tschaplinski, T. J., Muchero, W., and Chen, J. G. (2021). Transcriptional regulation of drought response in Arabidopsis and woody plants. Frontiers in Plant Science, 11:572137. https://doi.org/10.3389/fpls.2020. 572137. Zeng, L., Gu, Z., Xu, M., Zhao, N., Zhu, W., Yonezawa, T., Liu, T., Qiong, L., Tersing, T., Xu, L., Zhang, Y., Xu, R., Sun, N., Huang, Y., Lei, J., Zhang, L., Xie, F., Zhang, F., Gu, H., Geng, Y., Hasegawa, M., Yang, Z., Crabbe, M. J. C., Chen, F., and Zhong, Y. (2017). Discovery of a high-altitude ecotype and ancient lineage of Arabidopsis thaliana from Tibet. Science Bulletin, 62(24):1628–1630. https://doi.org/10.1016/j.scib. 2017.10.007. Zhao, W., Sun, Y.-Q., Pan, J., Sullivan, A. R., Arnold, M. L., Mao, J.-F., and Wang, X.- R. (2020). Effects of landscapes and range expansion on population structure and local adaptation. New Phytologist, 228:330–343. https://doi.org/10.1111/nph.16619. Zou, Y.-P., Hou, X.-H., Wu, Q., Chen, J.-F., Li, Z.-W., Han, T.-S., Niu, X.-M., Yang, L., Xu, Y.-C., Zhang, J., Zhang, F.-M., Tan, D., Tian, Z., Gu, H., and Guo, Y.- L. (2017). Adaptation of Arabidopsis thaliana to the Yangtze river basin. Genome Biology, 18(1):239. https://doi.org/10.1186/s13059-017-1378-9. ˚ Agren, J. and Schemske, D. W. (2012). Reciprocal transplants demonstrate strong adaptive differentiation of the model organism Arabidopsis thaliana in its native range. New Phytol, 194(4):1112–1122. https://doi.org/10.1111/j.1469-8137.2012.04112.x. 115 Bibliography ˇ Sr´amkov´a Fuxov´a, G., Z´avesk´a, E., Kol´aˇr, F., Luˇcanov´a, M., ˇ Spaniel, S., and Marhold, K. (2017). Range-wide genetic structure of Arabidopsis halleri (Brassicaceae): glacial persistence in multiple refugia and origin of the Northern Hemisphere disjunction. Botanical Journal of the Linnean Society, 185:321–342. https://doi.org/10.1093/ botlinnean/box064. 116 A. Appendix: Supplementary data Table A.1.: Region, accession code, altitude (m.a.s.l.) and coordinates (degrees) of relict A. thaliana populations known to date (Figure 2.1). Unknown altitudes are indicated by n/a. Accessions sharing the same geographic coordinates are given once and the others indicated by dashes. Region Accession Altitude Long(x) Lat(y) Algeria Alg-h1 1950 4.1500 36.4500 Canary Islands Can-0 n/a -13.4811 29.2144 Cape Verde Cvi-0 n/a -23.6167 15.1111 China Tibet-0 4200 90.9338 29.6903 China SRR2204316 n/a 98.6300 28.0200 China SRR2204703 – – – Iberian Peninsula Bad-1 720 -7.5548 41.5693 Iberian Peninsula Bdm-0 2124 -3.9915 40.7842 Iberian Peninsula Bon-0 1 -6.2898 36.8768 Iberian Peninsula Can-0 879 -4.5141 39.3884 Iberian Peninsula Cat-0 733 -3.6870 40.5421 Iberian Peninsula Cem-0 898 -4.3173 41.1484 Iberian Peninsula Con-0 756 -5.5991 37.9382 Iberian Peninsula Don-0 23 -6.3598 36.8323 Iberian Peninsula Fun-0 1804 -4.0503 40.7872 Iberian Peninsula Gra-0 1284 -5.3949 36.7678 Iberian Peninsula Hee-0 787 -3.4268 39.3048 Iberian Peninsula Her-12 584 -5.7828 39.3948 Iberian Peninsula Hum-2 924 -3.6889 42.2272 Iberian Peninsula Las-0 1071 -4.1527 40.5680 Iberian Peninsula Lso-0 769 -3.1629 38.8577 Iberian Peninsula 01/03/01 1030 -3.9316 39.5813 Iberian Peninsula Mdc-0 725 -3.5296 38.8743 Iberian Peninsula Med-0 771 -6.7657 42.4552 117 A. Appendix: Supplementary data Table A.1.: (continued) Region Accession Altitude Long(x) Lat(y) Iberian Peninsula Moj-0 705 -5.2814 36.7567 Iberian Peninsula Mor-0 803 -6.0732 40.8705 Iberian Peninsula Nac-0 1420 -3.9939 40.7510 Iberian Peninsula Nav-0 735 -4.6466 40.4177 Iberian Peninsula Ped-0 1096 -3.8978 40.7385 Iberian Peninsula Per-0 472 -1.1216 37.5952 Iberian Peninsula Pol-1 520 -4.2512 39.8163 Iberian Peninsula Pub-0 906 -6.6334 42.0567 Iberian Peninsula Pun-0 826 -4.7757 40.4033 Iberian Peninsula Rdc-0 800 -4.2866 40.4517 Iberian Peninsula Sca-0 501 -5.3700 37.9622 Iberian Peninsula Slc-3 945 -3.7443 38.4706 Iberian Peninsula Sln-22 646 -3.9948 40.3569 Iberian Peninsula Sne-0 1650 -3.2625 37.0746 Iberian Peninsula Toe-0 929 -3.9151 40.5925 Iberian Peninsula Vdc-0 593 -5.0302 38.0498 Iberian Peninsula Vel-0 2662 -3.3818 37.0854 Iberian Peninsula Ven-0 1553 -4.0141 40.7585 Iberian Peninsula Vim-0 652 -6.5101 41.8771 Iberian Peninsula Vis-0 244 -6.0455 39.8460 Iberian Peninsula Zar-0 1250 -4.1948 40.5518 Italy Etna-2 1800 14.9800 37.6900 Italy Mr-0 1000-1500 9.6500 44.1500 Lebanon Qar-8a 1160-1213 35.8400 34.1000 Madeira A1 1812 -16.9300 32.7400 Madeira A2 – – – Madeira A3 – – – Madeira A4 – – – Madeira A5 – – – Madeira A6 1600 -16.9370 32.7450 Madeira R1 1200 -17.1300 32.7540 Madeira R2 – – – Madeira R3 – – – Madeira R4 1200 -17.1320 32.7580 Madeira U1 1800 -16.9430 32.7590 118 Table A.1.: (continued) Region Accession Altitude Long(x) Lat(y) Morocco Agl-0 1530 -5.4486 32.9724 Morocco Ait-1 2075 -7.8129 31.2366 Morocco Arb-0 1070 -7.5262 31.4199 Morocco Azr-0 1726 -5.1791 33.4236 Morocco Bab-0 1186 -5.0237 35.0426 Morocco Bba-2 1320 -4.0839 34.0193 Morocco Bbe-0 1340 -4.8314 34.9952 Morocco Elh-2 1510 -7.4064 31.4720 Morocco Elk-1 1230 -6.0150 32.5352 Morocco Ifr-2 1605 -5.1747 33.5501 Morocco Ket-10 1614 -4.6661 34.9608 Morocco Khe-0 1470 -5.5103 32.9274 Morocco Meh-0 1336 -4.0515 33.9561 Morocco Oua-0 1825 -6.2753 32.0785 Morocco Set-0 1430 -7.6736 31.2266 Morocco Tah-0 1224 -4.2225 34.0529 Morocco Taz-0 1386 -4.1026 34.0917 Morocco Til-2 1420 -6.2296 32.0421 Morocco Tiz-0 1433 -4.0265 33.8723 Morocco Zin-9 116 -5.4270 35.4528 Morocco Toufl-1 n/a -7.4217 31.4687 Morocco Aitba-1 n/a -7.4499 31.4840 Morocco Aitba-2 – – – Morocco Ita-0 2000 -4.1989 34.0787 South Africa SA-h1 n/a 19.3750 -34.1250 South Africa SA-h2 – – – South Africa SA-h3 – – – South Africa SA-h4 n/a 19.2820 -33.3990 Tanzania Tanz-1 2775 36.1200 -2.8739 Tanzania Tanz-2 – – – Tanzania TZ- Kilamanjaro n/a 37.1596 -3.0216 Tanzania TZ-Mt. Meru n/a 36.7167 -3.2333 119