\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[T2A]{fontenc}
\usepackage[ngerman,polish,english]{babel}



\usepackage{float}








\begin{document}

\title{The holocentric chromosome microevolution: from phylogeographic patterns
to genomic associations with environmental gradients}



\author[1]{José Ignacio Márquez-Corro}%
\author[2]{Santiago Martín-Bravo}%
\author[3]{José Luis Blanco-Pastor}%
\author[4]{Modesto Luceño}%
\author[3]{Marcial Escudero}%
\affil[1]{Royal Botanic Gardens Kew}%
\affil[2]{Affiliation not available}%
\affil[3]{University of Seville}%
\affil[4]{Pablo de Olavide University}%


\vspace{-1em}



  
  \date{June 15, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
Geographic isolation and chromosome evolution are two of the major
drivers of diversification in eukaryotes in general, and specifically,
in plants. On one hand, range shifts induced by Pleistocene glacial
oscillations deeply shaped the evolutionary trajectories of species in
the Northern Hemisphere. On the other hand, karyotype variability within
species or species complexes may have adaptive potential as different
karyotypes may represent different recombination rates and linkage
groups that may be associated with locally adapted genes or supergenes.
Organisms with holocentric chromosomes are ideal to study the link
between local adaptation and chromosome evolution, due to their high
cytogenetic variability, especially when it seems to be related to
environmental variation. Here, we integrate the study of the
phylogeography, chromosomal evolution and ecological requirements of a
plant species complex distributed in the Western Euro-Mediterranean
region (Carex gr. laevigata, Cyperaceae). We aim to clarify the relative
influence of these factors on population differentiation and ultimately
on speciation. We obtained a well-resolved RADseq phylogeny that sheds
light on the phylogeographic patterns of molecular and chromosome number
variation, which are compatible with south-to-north postglacial
migration. In addition, landscape genomics analyses identified candidate
loci for local adaptation, and also strong significant associations
between the karyotype and the environment. We conclude that karyotype
distribution in C. gr. laevigata has been constrained by both range
shift dynamics and local adaptation. Our study demonstrates that
chromosome evolution may be responsible, at least partially, for
microevolutionary patterns of population differentiation and adaptation
in Carex.%
\end{abstract}%



\sloppy


\emph{Original article}

\textbf{The holocentric chromosome microevolution: from phylogeographic
patterns to genomic associations with environmental gradients}

\selectlanguage{ngerman}\textbf{José Ignacio Márquez-Corro\textsuperscript{1,2}, Santiago
Martín-Bravo\textsuperscript{1}, José Luis
Blanco-Pastor\textsuperscript{3}, Modesto Luceño\textsuperscript{1} and
Marcial Escudero\textsuperscript{3}}

\textsuperscript{1}Departamento de Biología Molecular e Ingeniería
Bioquímica, Universidad Pablo de Olavide, Carretera de Utrera km 1,
41013, Seville, Spain

\textsuperscript{2}Jodrell Laboratory, Royal Botanic Gardens, Kew, TW9
3AB, Richmond, United Kingdom

\textsuperscript{3}Departamento de Biología Vegetal y Ecología,
Universidad de Sevilla, 41012, Seville, Spain

\textbf{Author for correspondence}

Marcial Escudero. Email:\emph{amesclir@gmail.com}

\textbf{Short title:} Landscape genomics and chromosome evolution

\section*{Abstract}

{\label{abstract}}

Geographic isolation and chromosome evolution are two of the major
drivers of diversification in eukaryotes in general, and specifically,
in plants. On one hand, range shifts induced by Pleistocene glacial
oscillations deeply shaped the evolutionary trajectories of species in
the Northern Hemisphere. On the other hand, karyotype variability within
species or species complexes may have adaptive potential as different
karyotypes may represent different recombination rates and linkage
groups that may be associated with locally adapted genes or supergenes.
Organisms with holocentric chromosomes are ideal to study the link
between local adaptation and chromosome evolution, due to their high
cytogenetic variability, especially when it seems to be related to
environmental variation. Here, we integrate the study of the
phylogeography, chromosomal evolution and ecological requirements of a
plant species complex distributed in the Western Euro-Mediterranean
region (\emph{Carex} gr. \emph{laevigata} , Cyperaceae). We aim to
clarify the relative influence of these factors on population
differentiation and ultimately on speciation. We obtained a
well-resolved RADseq phylogeny that sheds light on the phylogeographic
patterns of molecular and chromosome number variation, which are
compatible with south-to-north postglacial migration. In addition,
landscape genomics analyses identified candidate loci for local
adaptation, and also strong significant associations between the
karyotype and the environment. We conclude that karyotype distribution
in \emph{C.} gr. \emph{laevigata} has been constrained by both range
shift dynamics and local adaptation. Our study demonstrates that
chromosome evolution may be responsible, at least partially, for
microevolutionary patterns of population differentiation and adaptation
in \emph{Carex} .

\textbf{Keywords:} \emph{Carex} , Cyperaceae, chromosomal speciation,
hybrid dysfunction, landscape genomics, suppression recombination

\section*{Introduction}

{\label{introduction}}

Chromosome evolution is one of the most important drivers of
biodiversity (Stebbins, 1950; Butlin, 2005; Ayala \& Coluzzi, 2005;
Hoffmann \& Rieseberg, 2008). Chromosomal rearrangements are often
associated with speciation (Rieseberg, 2001), as they have the potential
to reduce gene flow between diverging populations (Grant, 1963). In
angiosperms, karyotypic diversity is significantly associated with both
clade species richness and diversification rates, which suggests
chromosomal rearrangements promote or reinforce the speciation process
(Carta \& Escudero, 2023). Much research into chromosome evolution in
plants has focused on the physiological and ecological implications of
polyploidy (e.g. Otto \& Whitton, 2000; Otto, 2007), while the
evolutionary consequences of dysploidy have received comparatively much
less attention (Escudero et al., 2014).

Speciation mediated by chromosome evolution has been explained by two
models. First, the hybrid-dysfunction model of chromosomal speciation
presumes reduced fitness of hybrids between chromosome races. This model
is undermined by the fact that strong selection against structural
heterozygotes (underdominance) would make new cytotypes extremely
difficult to succeed, meaning that underdominant mutations are not
likely to cause speciation (Coyne \& Orr, 2004). The alternative model,
the suppressed recombination model of chromosomal speciation is better
supported theoretically as it does not involve underdominance. In this
case, chromosomal rearrangements promote divergence between populations
by preventing recombination in clusters of locally coadapted gene
complexes (supergenes; Schwander et al., 2014; Black \& Shuker, 2019).
These rearrangements may facilitate species divergence even in the
presence of ongoing gene flow between chromosomatically differentiated
populations (Navarro \& Barton, 2003; Ortiz-Barrientos et al., 2002;
Noor et al., 2001; Lowry \& Willis, 2010).

In holocentric chromosomes, centromeric regions are distributed along
the entire length of the chromosome, which may therefore attach to
microtubules during mitosis and meiosis (Márquez-Corro et al., 2019a).
In contrast, monocentric chromosomes have microtubule attachments
localized exclusively in a single centromeric region. Holocentric
chromosome organization has been described for ca. 20 lineages in three
of the six kingdoms in the domain Eukarya (the Eukaryotes): plants
(angiosperms and green algae); animals (at least six species-rich
arthropod clades, plus velvet worms and nematodes); and Rhizaria
(Escudero et al., 2016; Márquez-Corro et al., 2018). As a consequence,
approximately, ca. 20\% of all Eukaryotes have holocentric chromosomes
(Márquez-Corro et al., 2018). The evolutionary implications of
holocentricity are potentially profound but largely misunderstood.
Because of these special features of holocentric chromosomes, they have
been suggested as ideal systems to test the proposed models of
chromosomal speciation (Lucek et al., 2022).

Chromosome fragments from fissions that would be lost in monocentric
chromosomes may be inherited and become fixed in organisms with
holocentric chromosomes (Márquez-Corro et al., 2019a). Likewise,
enlarged chromosomes from fusion events can align and segregate
correctly in holocentric chromosomes. Conversely, in organisms with
monocentric chromosomes fusions usually result in the formation of
dicentric chromosomes that fail to segregate properly (Márquez-Corro et
al., 2019a). Consequently, large series of chromosome numbers and high
rates of karyotype evolution are found in holocentric lineages
(\emph{i.e.} , butterflies and moths --Lepidoptera-- 2\emph{n} = 10 to
250, sedges --Cyperaceae-- 2\emph{n} = 4 to 226; de Vos et al. 2020;
Márquez-Corro et al., 2019b, 2021). In Lepidoptera, phylogenetic
comparative evidence suggests that chromosome fission and fusion drive
cladogenesis (de Vos et al. 2020; Augustijnen et al., 2023). In the
sedge genus \emph{Carex} (Cyperaceae), chromosome rearrangements
contribute to genetic diversity within species (Hipp et al., 2010;
Escudero et al., 2013a) and the karyotype diversity is positively
associated with the time of coalescence of the species (Escudero et al.,
2010). Finally, chromosome number changes from fission and fusion are
related to linkage groups and, in holocentrics, they also determine
recombination rates as the number of chiasmata during meiosis is
directly proportional to chromosome number (Nokkala et al., 2004). As a
consequence, linkage groups and recombination rates (chromosome numbers)
have been positively associated with adaptive speciation towards
different environmental conditions in holocentric sedges (Escudero et
al., 2012, 2013b; Spalink et al., 2018; Márquez-Corro et al., 2021).
Some of the patterns of chromosome number association with environmental
conditions are rather weak at the macroevolutionary level (Escudero et
al., 2012; Spalink et al., 2018; Márquez-Corro et al., 2021), probably
based on the fact that chromosome evolution is much faster than
speciation rates in sedges (Márquez-Corro et al., 2021).

We hypothesize that holocentric chromosomes may be playing an important
role in species' ability to adapt and colonize new environments.
Therefore, having a direct impact on niche expansion for some of the
most widespread species of the study group. Karyotype rearrangements may
leave a traceable signature in the different individuals' genomes,
probably linked to specific loci in the linkage groups. Landscape
genomics constitutes a currently developing approach that may help to
identify specific genetic markers, including chromosome number
variation, associated with environmental variables and therefore
involved in local adaptation processes (Rellstab et al., 2015; Hoban et
al, 2016; Ahrens et al., 2018). To our knowledge, no study using
landscape genomics has addressed this topic in organisms with
holocentric chromosomes.

The goals of this study are: (i) to infer patterns of holocentric
chromosome evolution at phylogeography scale in a lineage of four
closely related and recently diversified \emph{Carex} species, (ii) to
identify genomic signatures of adaptive evolution and investigate the
associations between chromosome numbers and local environmental
conditions, (iii) to elucidate the role of cytogenetic evolution on the
diversification of the group and specifically the model of chromosomal
speciation.

\section*{Materials and Methods}

{\label{materials-and-methods}}

\emph{Study group}

We selected a lineage within \emph{Carex} sect.
\emph{Spirostachyae}(hereafter \emph{C.} gr. \emph{laevigata} ) which is
composed of four species (\emph{C. laevigata} , \emph{C. binervis} ,
\emph{C. camposii} and\emph{C. paulo-vargasii} ) endemic to western
Mediterranean-Atlantic Europe (Fig. 1), due to several characteristics
that make it ideal for the study of the possible association of
chromosome number variation and specific genomic loci with environmental
gradients. In the first place, their well-known cytogenetic variability,
with a representative sampling of chromosome counts across their
relatively restricted distribution range (Luceño \& Castroviejo, 1991).
In the second place, their wide chromosome number variation at the
inter-, intraspecific (\emph{C. binervis} , 2\emph{n} = 72--74; \emph{C.
camposii} , 2\emph{n} = 72;\emph{C. paulo-vargasii} , 2\emph{n} =
74--75; \emph{C. laevigata} , 2\emph{n} = 69--84), and even
intrapopulation level, especially for the more widely distributed
species (\emph{C. binervis} , \emph{C. laevigata} ; Fig. 1). Moreover,
chromosome number variation in \emph{C. laevigata} appears to be
influenced by population latitude (Fig. 1), which suggests a role of
chromosome evolution in the adaptation to the environment (Luceño \&
Castroviejo, 1991; Escudero et al., 2013b).\emph{Carex} gr.
\emph{laevigata} has a relatively recent origin (Pleistocene crown
diversification; ca. 2 mya, Martín-Bravo et al., 2019), and complex,
previously unresolved evolutionary relationships, in which rapid
chromosome speciation have been suggested to play an important role
(Escudero et al., 2013b).

\emph{DNA extraction, sequencing and data clustering}

A total of 189 individuals from 33 populations (Fig. 1; Table S1) of the
four species of \emph{C.} gr. \emph{laevigata} (plus one sample each
of\emph{C. mairei} --sister group--, \emph{C. distans} , \emph{C.
boryana}and \emph{C. borbonica} as outgroup) were sampled for
Restriction-site-Associated DNA sequencing (RAD-seq). This technique has
been proven as a useful High Throughput Sequencing (HTS) approach for
investigating recent evolutionary processes in non-model organisms,
including different scenarios of speciation among closely related
species (Otero et al., 2022 and references therein), in which
traditional Sanger-sequenced nuclear and plastid markers have been
unable to solve phylogenetic relationships (Escudero et al., 2013b). A
special effort was conducted for \emph{C. laevigata} in order to obtain
a wide intraspecific genomic sampling representative of chromosome
number variation within and across populations. Thus, 130 individuals
from 20 populations of \emph{C. laevigata} , most of which with
available chromosome counts, were included. DNA was extracted using a
DNeasy Plant Mini Kit (Qiagen, Valencia, CA, USA). RAD-seq libraries
preparation using restriction enzyme PstI were followed by sonication
and barcoding performed by Floragenex (Eugene, Portland, OR, USA)
following Baird\emph{et al.} (2008).

Raw data were demultiplexed based on the code for each sample and
adapters were removed from the reads using ipyrad v0.7.24 (Eaton \&
Overcast, 2020) using a reference genome (\emph{Carex scoparia} ; Planta
et al., 2022). We used the matrix resulting from a clustering threshold
of 90\% similarity and a minimum sequencing depth of 6 sequences per
locus for within-sample clustering (c90m6). We compared this matrix to
others with a clustering of 95\% and sequencing depth of 6 and 11, but
no significant difference was found. Then we discarded those samples
with higher missing data. Finally, we reanalyzed the data according to
the c90m6 parameters and a minimum number of samples with data at a
given locus of 90 out of 156 (ca. 60\%). We discarded those samples with
\textgreater{}95\% of missing data for the phylogeny construction, so a
final matrix of 156 accessions was retained (three outgroup taxa, one
sister taxa and 152 \emph{Carex} gr. \emph{laevigata} samples). The
phylogenetic relationships of the retained samples were reconstructed
using IQTREE v1.6.11 software (Nguyen \emph{et al.} , 2015;
Hoang\emph{et al} ., 2018).

\emph{Environment-associated loci (RDA analysis)}

We used redundancy analysis (RDA) to detect candidate genes under
environmental selection in \emph{C.} gr. \emph{laevigata} . RDA can be
used as a genotype-environment association (GEA) method to detect loci
under selection based on multivariate ordination (Forester et al.,
2016). RDA determines how groups of loci (here also chromosome number)
covary in response to the multivariate environment and can detect
processes that result in weak, multilocus molecular signatures (Rellstab
et al., 2015; Forester et al., 2018). Compared to other methods, RDA has
shown a superior combination of low false-positive and high
true-positive rates across a variety of selection scenarios (Forester et
al., 2018). Another advantage of RDA is that it can be used to analyze
many loci and environmental predictors simultaneously. For RDA we
considered only SNP sites for which there was sequence coverage in at
least 100 individuals, and for which the less-common allele was present
in at least 10\% of sampled individuals {[}minor allele frequency (MAF)
filter \selectlanguage{english}[?]0.10{]}. The chromosome number was added as an additional column
in the SNPs database. RDA is a regression-based method, and so it is
subject to problems when using highly correlated predictors. Hence, we
removed correlated predictors with a correlation value of
r\textgreater{}0.7. We used bioclimatic variables from the Worldclim 2
database (Fick \& Hijmans, 2017). The variable reduction was guided by
an ecological interpretation of the relevance of possible predictors. We
implemented a variable reduction protocol: First, we performed cluster
analyses of factors according to a matrix of absolute correlation values
\textbar{}r\textbar{}. For that, we used the complete linkage clustering
method of the `hclust' function in R (R Core Team, 2021). After
subsequent cluster analyses, we retained one variable in clusters with a
distance among variables lower than 0.3 (correlation higher than 0.7).
We favored quarterly over monthly variables.

We also checked for multicollinearity using variance inflation factors
(VIF) and confirmed that VIF of selected variables was \textless{}10. We
also performed a permutation test using 1000 permutations to assess the
statistical significance of environmental variables and RDA axes used in
the models. We used the loadings of the SNPs (stored as species in the
RDA object from the R environment) in the ordination space to determine
which SNPs are candidates for local adaptation. We extracted the SNP
loadings from the two significant constrained axes. Outlier SNPs are
those that load in the tails of SNP loading distributions. Outliers were
identified as SNPs with the greatest loadings along the significant RDA
axes (i.e., those in the 2.5\% upper and lower tails; Capblancq et al.,
2018), and were putatively considered as extremely associated with
environmental variables. Then we investigated the correlations between
environmental predictors and outlier SNPs. Finally, we represented SNPs
in the ordination space and color-coded them based on the predictor
variable with which they are most strongly correlated. Outlier SNPs with
the same signal as chromosome number were selected for investigating
their functional annotation, thus providing additional evidence of their
role in environmental adaptation. We used their genomic position to
identify candidate genes potentially under environmental selection
located \textless{}50 Kbp upstream or downstream in the
annotated\emph{Carex scoparia} genome (Planta et al., 2022).

\textbf{Results}

\emph{RAD-seq output}

After filtering the data through the ipyrad pipeline, we obtained two
datasets: i) one composed of 156 samples, with 1761 loci and 10138 SNPs
filtered, that was subsequently used for phylogenetic analysis, and ii)
another retaining 153 samples (three outgroup samples excluded) with 698
SNPs filtered, that was used for RDA analyses.

\emph{Phylogenetic relationships}

The phylogenetic relationships among \emph{C.} gr.
\emph{laevigata}species established by the c90m6 matrix of RAD-seq SNPs
(Fig. 2, Fig. S1) are congruent with the previous studies but resolved
the intermingled relationships that have been previously obtained
for\emph{C. laevigata-C. binervis} (Escudero et al., 2009; Escudero \&
Luceno, 2010). \emph{Carex paulo-vargasii} and \emph{C. camposii} were
retrieved as subsequent sisters to the remaining species. \emph{Carex
laevigata} constitutes a monophyletic species, while \emph{C.
binervis}was retrieved as polyphyletic, with two distinct lineages, an
early-diverging one (sister to all \emph{C.} gr. \emph{laevigata} ) that
is composed of populations from France and the British Islands
(Eurosiberian part of the species range), and an Iberian (Mediterranean
part of the species range) \emph{C. binervis} lineage that was retrieved
as sister to \emph{C. laevigata} . In turn, \emph{C. laevigata} was
composed of a main clade restricted to the southern Iberian Peninsula
(2\emph{n} = 72--76), sister to another main clade with a wider
distribution, spanning the species range and chromosome number
variability (2\emph{n} = 69--84).

\emph{RDA analysis}

Uncorrelated variables (r\textless{}0.7) retained for the RDA analysis
were BIO4 (temperature seasonality --standard deviation \selectlanguage{english} 100--), BIO6
(minimum temperature of the coldest month), BIO16 (precipitation of
wettest quarter) and BIO17 (precipitation of driest quarter).
Permutation tests showed statistical significance for the first two axes
and all four environmental variables (p-value\textless{}0.01).
Bioclimatic variables included in the RDA explained 17.2\% of the
genetic variation (inertia of constrained axes), and the remaining
82.8\% of the variation is unexplained and not accounted for by the
included bioclimatic variables (inertia of unconstrained axes).
Eigenvalues for the two significant constrained axes were 84.65
(proportion of variance explained (PVE): 0.70) for RDA1 and 19.41 (PVE:
0.16) for RDA2. RDA analysis (Fig. 3a) showed a clear association
between \emph{C. paulo-vargasii} genotypic information and BIO4, and
between Iberian \emph{C. binervis} populations and BIO17. On the other
hand, Eurosiberian \emph{C. binervis} populations, in congruence with
the phylogenetic relationships, formed a distinct cluster that
overlapped with \emph{C. camposii} . Finally, \emph{C.
laevigata}populations showed a wider spread across RDA axes, although
two groups of populations mostly fitted BIO16 and BIO17 values, whereas
another group placed more along BIO6.

We obtained a set of 72 outlier loci highly putatively associated with
the set of retained environmental predictors (Fig. 3b). Thus, BIO6,
BIO17, BIO16 and BIO4 showed the strongest correlation with 43, 14, 8
and 7 outlier SNPs, respectively. Remarkably, the chromosome number of
individuals used as a response variable was retrieved as an ``outlier
locus'' in the RDA analysis and showed the strongest correlation with
BIO6 (Fig. 3b).

There are four loci that are significantly associated with BIO6 and
highly correlated with chromosome number (Table S2). One of the four
loci is linked to genome regions with high amounts of repetitive DNA.
The other three loci are linked to important metabolic genes. One locus
was found within an urease subunit (related to urea metabolism; although
four other important genes were very close: lipoxygenase protein,
aspartyl aminopeptidase, adenylate kinase and homologous to rice gene
LOC). A second locus was found within a diacylglycerol kinase (DGK,
related to diacylglycerol metabolism, a lipid from cell membrane; a
second protein kinase domain was found very close). Finally, the last
locus was quite close to a formyl transferase.

\section*{Discussion}

{\label{discussion}}

\emph{Phylogeography constrains chromosome number distribution in the
western Mediterranean Carex gr. laevigata}

Phylogeographic and phylogenetic inference in \emph{C.}
gr.\emph{laevigata} has remained obscure in previous studies based on
Sanger sequencing of a few DNA regions, especially for the \emph{C.
laevigata} -\emph{C. binervis} species pair (Escudero et al., 2008,
2010; Escudero \& Luceno, 2009). On one hand, systematic studies focused
on the Carex sect. \emph{Spirostachyae} based on ITS and 5'\emph{trn} K
intron (Escudero et al., 2008, 2010; Escudero \& Luceno, 2009) retrieved
non-monophyly and intermingled phylogenetic relationships for \emph{C.
laevigata} -\emph{C. binervis} . On the other hand, the only previous
attempt to elucidate the phylogeographic structure of the group, based
on two plastid DNA regions (5\selectlanguage{english}'\emph{trn} K and the intergenic\emph{trn}
V-\emph{ndh} C) and a wide population sampling, obtained a highly
intricate haplotype network (Escudero et al., 2013b). Previously
inferred ribotype additivities and haplotype sharing between \emph{C.
laevigata} -\emph{C. binervis} have been attributed to hybridisation
(Escudero et al., 2008, 2013b; Escudero \& Luceno, 2009) or incomplete
lineage sorting (Escudero et al., 2013b). The phylogenetic resolution
provided by the highly increased amount of loci obtained across the
entire genomes of the study species with the RAD-seq approach (1761 loci
here vs. 1--2 in previous studies) allowed us to shed light on the main
phylogeographic patterns of \emph{C.} gr. \emph{laevigata} . While the
phylogenetic relationships for \emph{C. camposii} and \emph{C.
paulo-vargasii} are congruent with previous studies, we uncovered a
previously unknown genetic lineage for \emph{C. binervis} from the
Eurosiberian part of its range, and distinct from the clade
including\emph{C. binervis} samples from the Iberian Peninsula. Overall,
the southern Iberian peninsula-north Africa displayed higher
phylogenetic, taxonomic (all four species of \emph{C.} gr.
\emph{laevigata} present) and cytogenetic (most chromosome number range)
diversities (Figs. 1--2). Besides, endemics from that region (\emph{C.
paulo-vargasii} from Morocco and \emph{C. camposii} from Sierra Nevada),
as well as southern populations of \emph{C. binervis} and \emph{C.
laevigata} appeared as successive sisters to the rest of the lineages
(Fig. 2). On the other hand, only \emph{C. laevigata} and \emph{C.
binervis} are distributed in the northern Iberian Peninsula and the rest
of western Europe, where a reduced subset of their phylogenetic
diversity (Fig. 2) and chromosome numbers are found (2\emph{n} = 70--74,
Fig. 1). Therefore, the taxonomic, molecular and cytogenetic evidence
points to southern Iberia-north Africa (especially the region around the
Strait of Gibraltar) as the evolutionary cradle for \emph{C.}
gr.\emph{laevigata} . These phylogeographic and karyotypic patterns are
compatible with the signature of Pleistocene glaciations, including
southern refugia and northwards postglacial recolonization (e.g. Hewitt,
2011). This temporal-geographic scenario, which frequently caused
speciation (e.g. Kadereit \& Abbott, 2021), has been previously reviewed
for many plant groups (Rodriguez-Sanchez et al., 2008; Lavergne et al.,
2013; Nieto Feliner, 2014; Molina-Venegas et al., 2017).

\emph{Landscape genomics meets karyotype evolution}

Despite alleged limitations of RAD-seq genomic scanning to detect loci
related to local adaptation (Lowry et al., 2017), our landscape genomic
analyses have shown that 74 outlier loci are significantly correlated
with climatic variables related to temperature (BIO4 and BIO6) and
precipitation (BIO16 and BIO17). Interestingly, chromosome number was
among these outliers and was also correlated with BIO6. Together with
chromosome number we also identified four loci significantly correlated
with BIO6 and with similar associations with RDA1 and RDA2 axes as
chromosome number (Fig. 3b).

One of the four loci is linked to genome regions with high amounts of
repetitive DNA. This is unsurprising as highly repetitive sequences can
be found at centromeres and telomeres, which play critical roles in
genome integrity maintenance throughout the cell cycle (Onishi-Seebacher
\& Korbel, 2011). Repetitive DNA sequences were previously referred to
as ``junk DNA'' because few discernible functions could be assigned to
these regions, but recent studies have shown that they have important
functional roles in genome organization and evolution (von Sternberg \&
Shapiro, 2005; Biscotti et al., 2015). Although the precise functions of
repetitive regions in the genome are still poorly understood, in plants
they have been associated with abiotic stress resistance or adaptation
to the environment. For example, repetitive non-coding DNA has proven to
be susceptible to environmental stress, which can induce specific
quantitative modifications on the affected genome fraction (Bassi et
al., 2015). Repetitive DNA has been associated with ectopic
recombination (Gonzalez \& Petrov, 2012), also in holocentric organisms
(Hook et al., 2023), which could result in chromosomal fission and
fusion rearrangements and explain, at least partially, the resulting
patterns of karyotype variation. In this way, very recently, the genomic
architecture of fission and fusion has been deciphered for holocentric
sedges (Hofstaetter et al., 2022; Escudero et al., 2023) and butterflies
(Hook et al., 2023), which is characterized by a high density of
repetitive DNA.

The other three significantly correlated loci were associated with genes
involved in important metabolic routes, being kinases specially relevant
(three hits) and associated to two of the three loci. Among these,
diacylglycerol kinase (DGK) genes are of special relevance given that
they are involved in the modulation of plant growth, development, and
adaptation in both biotic and abiotic stress conditions (Kue Foka et
al., 2020). Additionally, DGK plays an important role in generating
membrane-derived oligosaccharides that protect cells against osmotic
stress conditions (Jefferson et al., 2013). Furthermore,
glycerophospholipids (GPL), which are synthesized by DGK, have been
linked with adaptation to high-altitude cold environments (Wei et al.,
2022). Kinase domains seem to be associated with stress response and
adaptation to the environment in plants. Receptor-like kinases (RLKs)
and histidine kinases have been found to play a role in plant response
to abiotic stresses (Osakabe et al., 2013). Additionally,
calcium-dependent protein kinases (CDPKs) are essential for plant
development and stress responses (Alves et al., 2021). SNF1-related
protein kinases have been found to have a unique function in plant
glucose and stress signaling (Lumbreras et al., 2001). Furthermore,
receptor-like kinases have been found to play significant roles in plant
growth regulation and responses to stresses (Fenglian et al., 2012).
Finally, kinase domains have also been found to function in plant growth
and salt-stress responses (Zhou et al., 2006).

Remarkably, different chromosome numbers (which are related to linkage
groups and recombination rates) are significantly associated with
allelic variants at these loci (with important genes), which in turn are
also significantly associated with environmental climatic variables.
This study is a step forward in supporting the hypothesis that
cytogenetic variation (chromosome number, linkage groups and
recombination rates) in species with holocentric chromosomes are
selected towards different optima of climatic regimes (Escudero et al.,
2012, 2013b; Spalink et al., 2018; Marquez-Corro et al., 2021).

\emph{Implications for speciation}

Karyotype diversity is directly associated with diversification rates
and species richness in angiosperms (Carta \& Escudero, 2023) which
suggests a key role of chromosome evolution in the plant diversification
process. Chromosome evolution has been hypothesized as a major driver of
diversification in the Mediterranean Basin (Escudero et al., 2018;
Thompson, 2020). In addition, chromosome evolution has been suggested to
play a major role in species diversification in holocentrics (de Vos et
al., 2020; Escudero et al., 2010, 2013b; Hipp et al., 2010). Finally,
chromosome evolution may also have an indirect role reinforcing
geographical or ecological speciation (Coyne \& Orr, 2004).

There are two overarching models for chromosomal speciation, hybrid
dysfunction and recombination suppression (Ayala \& Coluzzi, 2005). They
are not mutually exclusive and they can both play a partial role in the
process of chromosomal speciation, especially in holocentric chromosomes
where hybrid dysfunction is better supported theoretically than in
monocentrics (Lucek et al., 2022). The suggested adaptive chromosomal
speciation here could be helpful to support both models of chromosomal
speciation. On one hand, strong selection towards new chromosomal
variants may help them to become established in the populations which
can ameliorate one of the strongest critics against the hybrid
dysfunction model (the minority cytotype exclusion; Levin, 1975). On the
other hand, recombination suppression models already suggest locally
adapted genes in which recombination is protected by chromosome
rearrangements (Ayala \& Coluzzi, 2005). Accordingly, under this
chromosomal speciation model, we certainly expect that different
karyotypes are locally adapted. Our landscape genomics analyses suggest
that karyotypes are locally adapted in our study group. In addition, the
phylogenetic placement of sampled populations suggests multiple origins
of different chromosome numbers (e.g. 2\emph{n} = 72, 76), which is
congruent with recurrent chromosome fission and fusion processes in
relatively short time periods (Escudero et al., under review) and
further support its adaptive value. We hypothesize here that clusters of
locally adapted genes in \emph{C.} gr. \emph{laevigata} may have been
protected for recombination by the fission and fusion rearrangements. In
fact, fissions and fusions in holocentric have been observed to suppress
recombination in holocentrics in a study of comparative genomics (whole
sequence genome vs. linkage map, Escudero et al., 2023). For instance,
in \emph{C. laevigata} , this would have allowed different karyotypes
(and, consequently, gene combinations) to either disperse northwards
after the last glacial period or remain in the progressively warmer
habitats in the southern Iberian Peninsula, further contributing to
lineage divergence.

\section*{Acknowledgements}

{\label{acknowledgements}}

This project was funded by the Andalusian regional government (Spain)
through project ProyExcel\_00125 towards SMB, and by the Spanish
Ministry of Economy and Competitiveness through the BES-2017-079621
scholarship towards JIM-C. JIM-C was also granted by the `Next
Generation EU' funding, the Recovery Plan, Transformation and Resilience
and the Ministry of Universities, under the grants `Margarita Salas' for
the requalification of the Spanish university system 2021-2023 called by
the Universidad Pablo de Olavide, Seville. ME was supported by
FEDER/MICINN -- AEI (PID2021-122715NB-I00). JLB-P was supported by the
``Ministerio de Ciencia e Innovacion'' (Spain) and by the ``Agencia
Estatal de Investigacion'' through the CHRONOS project
(PID2020119327RJ-I00 MCIN/AEI/10.13039/501100011033). The authors thank
the Andalusian Scientific Information Technology Center (CICA, Spain)
for providing computational resources, and J. Viruel (RBG Kew, United
Kingdom) for the suggestions and comments at the early stages of the
analyses.

\textbf{Conflict of interest}

The authors have no conflicts of interest to declare.

\textbf{Authors contributions}

ME, SMB and ML designed the research; JIM-C, ME, SMB and ML authors
collected the material samples; ME and JIM-C prepared and sent the DNA
samples to Floragenex; JIM-C, ME, JLB-P analyzed the data; JIM-C, ME and
SMB wrote the manuscript; and ML and JLB-P contributed to the manuscript
writing.

\textbf{Data accessibility}

Obtained RAD-seq data are available on the NCBI under BioProject
accession PRJNA982410.

\textbf{ORCID}

Jose Ignacio Marquez-Corro\emph{https://orcid.org/0000-0003-4277-2933}

Santiago Martin-Bravo\emph{https://orcid.org/0000-0003-0626-0770}

Jose Luis Blanco-Pastor\emph{https://orcid.org/0000-0002-7708-1342}

Modesto Luceno\emph{https://orcid.org/0000-0002-0336-4083}

Marcial Escudero\emph{https://orcid.org/0000-0002-2541-5427}

\section*{References}

{\label{references}}

Ahrens, C. W., Rymer, P. D., Stow, A., Bragg, J., Dillon, S., Umbers, K.
D., \& Dudaniec, R. Y. (2018). The search for loci under selection:
trends, biases and progress. \emph{Molecular Ecology} , \emph{27} (6),
1342--1356. \emph{https://doi.org/10.1111/mec.14549}

Alves, H. L. S., Matiolli, C. C., Soares, R. C., Almadanim, M. C.,
Oliveira, M. M, \& Abreu, I. A. (2021). Carbon/nitrogen metabolism and
stress response networks -- calcium-dependent protein kinases as the
missing link?. \emph{Journal of Experimental Botany} , \emph{12} ,
4190--4201.\emph{https://doi.org/10.1093/jxb/erab136}

Augustijnen, H., Baetscher, L., Cesanek, M., Chkhartishvili, T., Dinca,
V., Iankoshvili, G., Ogawa, K., Vila, R., Klopfstein,S., de Vos, J. M.,
\& Lucek, K. (2023). A macroevolutionary role for chromosomal fusion and
fission in \emph{Erebia} butterflies. bioRxiv,
2023.01.16.524200.\emph{https://doi.org/10.1101/2023.01.16.524200}

Ayala, F. J., \& Coluzzi, M. (2005). Chromosome speciation:
humans,\emph{Drosophila} , and mosquitoes. \emph{Proceedings of the
National Academy of Sciences} , \emph{102} (suppl\_1),
6535--6542.\emph{https://doi.org/10.1073/pnas.0501847102}

Baird, N. A., Etter, P. D., Atwood, T. S., Currey, M. C., Shiver, A. L.,
Lewis, Z. A., Selker, E. U., Cresko, W. A., \& Johnson, E. A. (2008).
Rapid SNP discovery and genetic mapping using sequenced RAD
markers.\emph{PLoS One} , \emph{3} , 1--7.
https://doi.org/10.1371/journal.pone.0003376

Bassi, P., Basile, A., \& Sorbo, S. (2015) Behaviour of repetitive
non-coding DNA in response to heavy metal stress in the protonemata
of\emph{Funaria hygrometrica} . \emph{Plant Biosystems - An
International Journal Dealing with all Aspects of Plant Biology} ,
\emph{149} , 315--321, https://doi.org/10.1080/11263504.2013.824040

Biscotti, M. A., Olmo, E., \& Heslop-Harrison, J. S. (2015). Repetitive
DNA in eukaryotic genomes. \emph{Chromosome Research,} \emph{23} ,
415--420. https://doi.org/10.1007/s10577-015-9499-z

Black, D., \& Suker, D. M. (2019). Supergenes. \emph{Current Biology}
,\emph{29} , R615--R617. https://doi.org/10.1016/j.cub.2019.05.024

Butlin, R.K., (2005). Recombination and speciation. \emph{Molecular
Ecology, 14} , 2621--2635.

\emph{https://doi.org/10.1111/j.1365-294X.2005.02617.x}

Capblancq, T., Luu, K., Blum, M. G. B., \& Bazin, E. (2018). Evaluation
of redundancy analysis to identify signatures of local
adaptation.\emph{Molecular Ecology Resources, 18} (6), 1223--1233.

\emph{https://doi.org/10.1111/1755-0998.12906}

Carta, A., \& Escudero, M. (2023). Karyotypic diversity: a neglected
trait to explain angiosperm diversification?. \emph{Evolution, 77} (4),
1158--1164.\emph{https://doi.org/10.1093/evolut/qpad014}

Coyne, J. A., \& Orr, H. A. (2004). Speciation Sinauer Associates.
Sunderland, MA, 276, 281.

de Vos, J. M., Augustijnen, H., Batscher, L., \& Lucek, K. (2020).
Speciation through chromosomal fusion and fission in
Lepidoptera.\emph{Philosophical Transactions of the Royal Society B} ,
\emph{375} , 20190539.\emph{https://doi.org/10.1098/rstb.2019.0539}

Eaton, D. A. R., \& Overcast, I. (2020). ipyrad: Interactive assembly
and analysis of RADseq datasets. \emph{Bioinformatics, 36} , 2592--2594.

\emph{https://doi.org/10.1093/bioinformatics/btz966}

Escudero, M., Valcarcel, V., Vargas, P., \& Luceno, M. (2008). Evolution
in \emph{Carex} L. sect. \emph{Spirostachyae} (Cyperaceae): a molecular
and cytogenetic approach. \emph{Organisms Diversity \& Evolution, 7}
(4), 271--291.\emph{https://doi.org/10.1016/j.ode.2006.08.006}

Escudero, M., \& Luceno, M. (2009). Systematics and evolution
of\emph{Carex} sects. \emph{Spirostachyae} and \emph{Elatae}
(Cyperaceae).\emph{Plant Systematics and Evolution, 279} ,
163--189.\emph{https://doi.org/10.1007/s00606-009-0156-x}

Escudero, M., Hipp, A. L., \& Luceno, M. (2010). Karyotype stability and
predictors of chromosome number variation in sedges: a study
in\emph{Carex} section \emph{Spirostachyae} (Cyperaceae).
\emph{Molecular Phylogenetics and Evolution, 57} (1),
353--363.\emph{https://doi.org/10.1016/j.ympev.2010.07.009}

Escudero, M., Hipp, A. L., Hansen, T. F., Voje, K. L., \& Luceno, M.
(2012). Selection and inertia in the evolution of holocentric
chromosomes in sedges (\emph{Carex} , Cyperaceae). \emph{New
Phytologist, 195} (1),
237--247.\emph{https://doi.org/10.1111/j.1469-8137.2012.04137.x}

Escudero, M., Weber, J. A., \& Hipp, A. L. (2013a). Species coherence in
the face of karyotype diversification in holocentric organisms: the case
of a cytogenetically variable sedge (\emph{Carex scoparia,}
Cyperaceae).\emph{Annals of Botany, 112} (3),
515--526.\emph{https://doi.org/10.1093/aob/mct119}

Escudero, M., Maguilla, E., \& Luceno, M. (2013b). Selection by climatic
regime and neutral evolutionary processes in holocentric chromosomes
(\emph{Carex} gr. \emph{laevigata} : Cyperaceae): A microevolutionary
approach. \emph{Perspectives in Plant Ecology, Evolution and
Systematics, 15} (2),
118--129.\emph{https://doi.org/10.1016/j.ppees.2013.01.001}

Escudero, M., Martin-Bravo, S., Mayrose, I., Fernandez-Mazuecos, M.,
Fiz-Palacios, O., Hipp, A. L., Pimentel, M., Jimenez-Mejias, P.,
Valcarcel, V., Vargas, P., \& Luceno, M. (2014). Karyotypic changes
through dysploidy persist longer over evolutionary time than polyploid
changes. \emph{PLoS One, 9} (1), e85266.

\emph{https://doi.org/10.1371/journal.pone.0085266}

Escudero, M., Marquez-Corro, J. I., \& Hipp, A. L. (2016). The
phylogenetic origins and evolutionary history of holocentric
chromosomes. \emph{Systematic Botany, 41} (3), 580--585.

\emph{https://doi.org/10.1600/036364416X692442}

Escudero, M., Balao, F., Martin\selectlanguage{english}-Bravo, S., Valente, L., \& Valcarcel, V.
(2018). Is the diversification of Mediterranean Basin plant lineages
coupled to karyotypic changes?. \emph{Plant Biology, 20} ,
166--175.\emph{https://doi.org/10.1111/plb.12563}

Escudero, M., Marques, A., Lucek, K., \& Hipp, A. (2023). Highly
conserved synteny despite massive chromosome fusion and fission suggest
fragile sites in holocentric plants. \emph{Authorea}
.\emph{https://doi.org/10.22541/au.167593554.47685772/v1}

Fenglian, X., Jiaojiao, L., Chunhui, D., Na, W., Suqiao, Z. (2012).
OsSRL, a Receptor-like Kinase with S-locus Site, Is Involved in Drought
Tolerance in Rice. \emph{Chinese Bulletin of Botany} , \emph{47} ,
474--482.\emph{https://doi.org/10.3724/SP.J.1259.2012.00474}

Fick, S. E., \& Hijmans, R. J. (2017). WorldClim 2: New 1-km spatial
resolution climate surfaces for global land areas. \emph{International
Journal of Climatology, 37} (12),
4302--4315.\emph{https://doi.org/10.1002/joc.5086}

Forester, B. R., Jones, M. R., Joost, S., Landguth, E. L., \& Lasky, J.
R. (2016). Detecting spatial genetic signatures of local adaptation in
heterogeneous landscapes. \emph{Molecular Ecology, 25} (1),
104--120.\emph{https://doi.org/10.1111/mec.13476}

Forester, B. R., Lasky, J. R., Wagner, H. H., \& Urban, D. L. (2018).
Comparing methods for detecting multilocus adaptation with multivariate
genotype--environment associations. \emph{Molecular Ecology, 27} (9),
2215--2233. \emph{https://doi.org/10.1111/mec.14584}

Gonzalez, J., \& Petrov, D. A. (2012). Evolution of genome content:
population dynamics of transposable elements in flies and
humans.\emph{Methods in Molecular Biology, 855} ,
361--383.\emph{https://doi.org/10.1007/978-1-61779-582-4\_13}

Grant, V. (1963). The origin of adaptations. Columbia University Press:
New York. 704 pp.

Hewitt, G. M. (2011). Mediterranean peninsulas: the evolution of
hotspots. In: Zachos, F., \& Habel, J. (eds.). \emph{Biodiversity
hotspots: Distribution and protection of conservation priority areas} ,
123--147. Springer:
Berlin.\emph{https://doi.org/10.1007/978-3-642-20992-5\_7}

Hipp, A. L., Rothrock, P. E., Whitkus, R., \& Weber, J. A. (2010).
Chromosomes tell half of the story: the correlation between karyotype
rearrangements and genetic diversity in sedges, a group with holocentric
chromosomes. \emph{Molecular Ecology, 19} (15),
3124--3138.\emph{https://doi.org/10.1111/j.1365-294x.2010.04741.x}

Hoang, D. T., Chernomor, O., von Haeseler, A., Minh, B. Q., Vinh, L. S.
(2018). UFBoot2: Improving the ultrafast bootstrap
approximation.\emph{Molecular Biology and Evolution,} \emph{35} ,
518--522.\emph{https://doi.org/10.1093/molbev/msx281}

Hoban, S., Kelley, J. L., Lotterhos, K. E., Antolin, M. F., Bradburd,
G., Lowry, D. B., Poss, M. L., Storfer, A., \& Whitlock, M. C. (2016).
Finding the genomic basis of local adaptation: pitfalls, practical
solutions, and future directions. \emph{The American Naturalist}
,\emph{188} (4), 379--397.\emph{https://doi.org/10.1086\%2F688018}

Hoffmann, A. A., \& Rieseberg, L. H. (2008). Revisiting the impact of
inversions in evolution: from population genetic markers to drivers of
adaptive shifts and speciation?. \emph{Annual Review of Ecology,
Evolution, and Systematics, 39} ,
21--42.\emph{https://doi.org/10.1146/annurev.ecolsys.39.110707.173532}

Hofstatter, P. G., Thangavel, G., Lux, T., Neumann, P., Vondrak, T.,
Novak, P., Zhang, M., Costa, L., Castellani, M., Scott, A., Toegeloa,
H., Fuchs, J., Mata-Sucre, Y., Dias, Y., Vanzela, A. L. L., Huettel, B.,
Almeida, C. C. S., \selectlanguage{polish}Š\selectlanguage{english}imkov\selectlanguage{ngerman}á, H., Souza, G., Pedrosa-Harand, A., MAcas,
J., Mayer, K. F. X., Hoube, A., \& Marques, A. (2022). Repeat-based
holocentromeres influence genome architecture and karyotype
evolution.\emph{Cell} , \emph{185} (17),
3153--3168.\emph{https://doi.org/10.1016/j.cell.2022.06.045}

Höök, L., Näsvall, K., Vila, R., Wiklund, C., \& Backström, N. (2023).
High-density linkage maps and chromosome level genome assemblies unveil
direction and frequency of extensive structural rearrangements in wood
white butterflies (\emph{Leptidea} spp.). \emph{Chromosome Research, 31}
(1), 2.\emph{https://doi.org/10.1007/s10577-023-09713-z}

Jefferson, R. E., Blois, R. M., Bowie, J. U. (2013). Membrane proteins
can have high kinetic stability. \emph{Journal of the American Chemical
Society, 40} , 15183--15190.\emph{https://doi.org/10.1021/ja407232b}

Kadereit, J. W., \& Abbott, R. J. (2021). Plant speciation in the
Quaternary. \emph{Plant Ecology \& Diversity}, \emph{14},
105--142.\emph{https://doi.org/10.1080/17550874.2021.2012849}

Kue Foka, I. C., Ketehouli, T., Zhou, Y., Li, X.-W., Wang, F.-W., \& Li,
H. (2020). The emerging roles of Diacylglycerol Kinase (DGK) in plant
stress tolerance, growth, and development. \emph{Agronomy} , \emph{10} ,
1375. \emph{https://doi.org/10.3390/agronomy10091375}

Lavergne, S., Hampe, A., \& Arroyo, J. (2013). In and out of Africa: how
did the Strait of Gibraltar affect plant species migration and local
diversification?. \emph{Journal of Biogeography}, \emph{40}(1),
24--36.\emph{https://doi.org/10.1111/j.1365-2699.2012.02769.x}Levin, D.
A. (1975). Minority cytotype exclusion in local plant populations.
\emph{Taxon}, \emph{24}(1),
35--43.\emph{https://doi.org/10.2307/1218997}Lowry, D. B., \& Willis, J.
H. (2010). A widespread chromosomal inversion polymorphism contributes
to a major life-history transition, local adaptation, and reproductive
isolation. \emph{PLoS Biology, 8}(9),
e1000500.\emph{https://doi.org/10.1371/journal.pbio.1000500}Lowry, D.
B., Hoban, S., Kelley, J. L., Lotterhos, K. E., Reed, L. K., Antolin, M.
F., \& Storfer, A. (2017). Breaking RAD: An evaluation of the utility of
restriction site\selectlanguage{english}-associated DNA sequencing for genome scans of
adaptation. \emph{Molecular Ecology Resources, 17}(2), 142--152.
\emph{https://doi.org/10.1111/1755-0998.12635}Lucek, K., Augustijnen,
H., \& Escudero, M. (2022). A holocentric twist to chromosomal
speciation?. \emph{Trends in Ecology \& Evolution, 37}(8),
655--662.\emph{https://doi.org/10.1016/j.tree.2022.04.002}Luceno, M., \&
Castroviejo, S. (1991). Agmatoploidy in \emph{Carex laevigata}
(Cyperaceae). Fusion and fission of chromosomes as the mechanism of
cytogenetic evolution in Iberian populations. \emph{Plant Systematics
and Evolution, 177}, 149--159.\emph{https://doi.org/10.1007/BF00937952}

Lumbreras, V., Alba, M. M., Kleinow, T., Koncz, C., \& Pages, M. (2001).
Domain fusion between SNF1\selectlanguage{english}-related kinase subunits during plant
evolution. \emph{EMBO Reports, 1} , 55--60.

\emph{https://doi.org/10.1093/embo-reports/kve001}

Marquez-Corro, J. I., Escudero, M., \& Luceno, M. (2018). Do holocentric
chromosomes represent an evolutionary advantage? A study of paired
analyses of diversification rates of lineages with holocentric
chromosomes and their monocentric closest relatives. \emph{Chromosome
Research, 26} , 139--152.

\emph{https://doi.org/10.1007/s10577-017-9566-8}

Marquez-Corro, J. I., Martin-Bravo, S., Pedrosa-Harand, A., Hipp, A. L.,
Luceno, M., \& Escudero, M. (2019a). Karyotype evolution in holocentric
organisms. In: \emph{Encyclopedia of Life Sciences} , John Wiley \&
Sons, Ltd:
Chichester.\emph{https://doi.org/10.1002/9780470015902.a0028758}

Marquez-Corro, J. I., Martin-Bravo, S., Spalink, D., Luceno, M., \&
Escudero, M. (2019b). Inferring hypothesis-based transitions in
clade-specific models of chromosome number evolution in sedges
(Cyperaceae). \emph{Molecular Phylogenetics and Evolution, 135} (4):
203--209. \emph{https://doi.org/10.1016/j.ympev.2019.03.006}

Marquez-Corro, J. I., Martin-Bravo, S., Jimenez-Mejias, P., Hipp, A. L.,
Spalink, D., Naczi, R. F. C., Roalson, E. H., Luceno, M., \& Escudero,
M. (2021). Macroevolutionary insights in sedges (\emph{Carex}:
Cyperaceae): The effects of rapid chromosome number evolution on lineage
diversification. \emph{Journal of Systematic and Evolution, 59}(4):
776--790.\emph{https://doi.org/10.1111/jse.12730}Martin-Bravo, S.,
Jimenez-Mejias, P., Villaverde, T., Escudero, M., Hahn, M., Spalink, D.,
Roalson, E.H., Hipp, A.L., Benitez-Benitez, B., Bruederle, L.P., Fitzek,
E., Ford, B.A., Ford, K.A., Garner, M., Gebauer, S., Hoffmann, M.A.,
Xiao-Feng, J., Larridon, I., Leveille-Bourret, E., Lu, Y.-F., Luceno,
M., Maguilla, E., Marquez-Corro, J. I., Miguez, M., Naczi, R., Reznicek,
A. A., \& Starr, J. R. (2020). A tale of worldwide success: Behind the
scenes of\emph{Carex} (Cyperaceae) biogeography and
diversification.\emph{Journal of Systematics and Evolution, 57}(6),
695--718.\emph{https://doi.org/10.1111/jse.12549}Molina-Venegas, R.,
Aparicio, A., Lavergne, S., \& Arroyo, J. (2017). Climatic and
topographical correlates of plant palaeo-and neoendemism in a
Mediterranean biodiversity hotspot. \emph{Annals of
Botany},\emph{119}(2),
229--238.\emph{https://doi.org/10.1093/aob/mcw093}

Navarro, A., \& Barton, N. H. (2003). Chromosomal speciation and
molecular divergence--accelerated evolution in rearranged
chromosomes.\emph{Science, 300} (5617),
321--324.\emph{https://doi.org/10.1126/science.1080600}

Nguyen, L.-T., Schmidt, H. A., von Haeseler, A., \& Minh, B.Q. (2015).
IQ-TREE: A fast and effective stochastic algorithm for estimating
maximum likelihood phylogenies. \emph{Molecular Biology and
Evolution,32} , 268--274.\emph{https://doi.org/10.1093/molbev/msu300}

Nieto Feliner, G. (2014). Patterns and processes in plant phylogeography
in the Mediterranean Basin. A review. \emph{Perspectives in Plant
Ecology, Evolution and Systematics} , \emph{16} (5),
265--278.\emph{https://doi.org/10.1016/j.ppees.2014.07.002}

Nokkala, S., Kuznetsova, V. G., Maryanska-Nadachowska, A., \& Nokkala,
C. (2004). Holocentric chromosomes in meiosis. I. Restriction of the
number of chiasmata in bivalents. \emph{Chromosome Research, 12} ,
733--739.\emph{https://doi.org/10.1023/B:CHRO.0000045797.74375.70}

Noor, M. A., Grams, K. L., Bertucci, L. A., \& Reiland, J. (2001).
Chromosomal inversions and the reproductive isolation of
species.\emph{Proceedings of the National Academy of Sciences, 98} (21),
12084--12088.\emph{https://doi.org/10.1073/pnas.221274498}

Onishi-Seebacher, M., \& Korbel, J. O. (2011). Challenges in studying
genomic structural variant formation mechanisms: The short-read dilemma
and beyond. \emph{BioEssays} , \emph{33} (11),
840--850.\emph{https://doi.org/10.1002/bies.201100075}

Ortiz-Barrientos, D., Reiland, J., Hey, J., \& Noor, M. A. (2002).
Recombination and the divergence of hybridizing species. \emph{Genetica,
116} , 167--178.\emph{https://doi.org/10.1023/A:1021296829109}

Osakabe, Y., Yamaguchi-Shinozaki, K., Shinozaki, K., \& Tran, L. S. P.
(2013). Sensing the environment: key roles of membrane-localized kinases
in plant perception and response to abiotic stress. \emph{Journal of
Experimental Botany} , \emph{2} , 445--458.

\emph{https://doi.org/10.1093/jxb/ers354}

Otero, A., Vargas, P., Fernandez-Mazuecos, M., Jimenez-Mejias, P.,
Valcarcel, V., Villa-Machio, I., \& Hipp, A. L. (2022). A snapshot of
progenitor--derivative speciation in Iberodes
(Boraginaceae).\emph{Molecular Ecology, 31} ,
3192--3209.\emph{https://doi.org/10.1111/mec.16459}

Otto, S. P., \& Whitton, J. (2000). Polyploid incidence and
evolution.\emph{Annual Review of Genetics, 34} (1),
401--437.\emph{https://doi.org/10.1146/annurev.genet.34.1.401}

Otto, S. P. (2007). The evolutionary consequences of
polyploidy.\emph{Cell, 131} (3),
452--462.\emph{https://doi.org/10.1016/j.cell.2007.10.022}

Planta, J., Liang, Y.-Y., Xin, H., Chansler, M. T., Prather, L. A.,
Jiang, N., Jiang, J., \& Childs, K. L. (2022). Chromosome-scale genome
assemblies and annotations for Poales species \emph{Carex cristatella}
,\emph{Carex scoparia} , \emph{Juncus effusus} , and \emph{Juncus
inflexus} . \emph{G3 (Bethesda), 12} ,
jkac211.\emph{https://doi.org/10.1093/g3journal/jkac211}

R Core Team. (2021). R: A language and environment for statistical
computing. \emph{R Foundation for Statistical Computing} , Vienna,
Austria.\emph{https://Www.R-Project.Org}

Rellstab, C., Gugerli, F., Eckert, A. J., Hancock, A. M., \&
Holderegger, R. (2015). A practical guide to environmental association
analysis in landscape genomics. \emph{Molecular Ecology}, \emph{24}(17),
4348--4370.\emph{https://doi.org/10.1111/mec.13322}

Rieseberg, L. H. (2001). Chromosomal rearrangements and
speciation.\emph{Trends in Ecology \& Evolution, 16} (7),
351--358.\emph{https://doi.org/10.1016/S0169-5347(01)02187-5}

Rodriguez-Sanchez, F., Perez-Barrales, R., Ojeda, F., Vargas, P., \&
Arroyo, J. (2008). The Strait of Gibraltar as a melting pot for plant
biodiversity. \emph{Quaternary Science Reviews}, \emph{27},
2100--2117.\emph{https://doi.org/10.1016/j.quascirev.2008.08.006}

Schwander, T., Libbrecht, R., \& Keller, L. (2014). Supergenes and
complex phenotypes. \emph{Current Biology} , \emph{24} ,
R288--R294.\emph{https://doi.org/10.1016/j.cub.2014.01.056}

Spalink, D., Pender, J., Escudero, M., Hipp, A. L., Roalson, E. H.,
Starr, J. R., Waterway, M. J., Bohs, L., \& Sytsma, K. J. (2018). The
spatial structure of phylogenetic and functional diversity in the United
States and Canada: An example using the sedge family
(Cyperaceae).\emph{Journal of Systematics and Evolution, 56} (5),
449--465.\emph{https://doi.org/10.1111/jse.12423}

Stebbins, G. L. (1950). Variation and evolution in plants. In Variation
and evolution in plants. Columbia University
Press.\emph{https://doi.org/10.7312/steb94536}

Thompson, J. D. (2020). Plant evolution in the Mediterranean: insights
for conservation. Oxford University Press,
USA.\emph{https://doi.org/10.1093/oso/9780198835141.001.0001}

von Sternberg, R., \& Shapiro, J. A. (2005) How repeated retroelements
format genome function. \emph{Cytogenetic and Genome Research, 110} ,
108--116.\emph{https://doi.org/10.1159/000084942}

Wei, P., Yu, X.-L., Yang, Y.-J., Chen, Z.-Y., Zhao, S.-Q., Li, X.-Z.,
Zhang, W.-C., Liu, C.-L., Li, X.-Y., \& Liu, X. (2022). Biased gene
expression reveals the contribution of subgenome to altitude adaptation
in allopolyploid Isoetes sinensis. \emph{Ecology and Evolution}
,\emph{12} , e9677.\emph{https://doi.org/10.1002/ece3.9677}

Zhou, H.-L., Cao, W.-H., Cao, Y.-R., Liu, J, Hao, Y.-J., Zhang, J.-S.,
\& Chen, S.-Yi. (2006). Roles of ethylene receptor NTHK1 domains in
plant growth, stress response and protein phosphorylation. \emph{FEBS
Letters, 580} ,
239--1250.\emph{https://doi.org/10.1016/j.febslet.2006.01.037}

\section*{Figures and tables}

{\label{figures-and-tables}}

Fig. 1. \emph{Carex} group \emph{laevigata} species distribution and
chromosome number reports for the specimens included in the study. The
dashed lines represent the tentative distributions of the two lineages
of \emph{C. binervis} hypothesized by the phylogenetic analysis
(Eurosiberian vs. strictly Mediterranean). Distribution of the main
lineages is shown next to the tree.

Fig. 2. RAD-seq phylogeny of 152 samples of \emph{C.}
gr.\emph{laevigata} and one sample of \emph{C. mairei} (the three
outgroup samples were dropped) based on 1761 loci and 10138 SNPs
filtered from IQTREE analysis.

Fig. 3. The RDA analysis. (a) Dots represent the coordinates of
populations in the biplot of the first two RDA axes and dot colors
represent \emph{Carex} species. (b) Dots represent the coordinates of
SNPs in the biplot of the first two RDA axes and dot colors represent
environmental predictors showing the highest correlation values with
outlier SNPs. The projection of environmental predictors used in the RDA
analysis is also displayed.

\textbf{Supporting Information}

Fig. S1. RAD-seq phylogeny of 156 samples of \emph{Carex}
gr.\emph{laevigata} plus the outgroup.

Table S1. SRA accession, voucher and chromosome data for the material
of\emph{Carex} gr. \emph{laevigata} , related species and outgroups used
in the study.

Table S2. The four loci associated with chromosome number and BIO6 are
shown indicating locus number, chromosome, location in base pairs (bp),
the closest genome feature and its location and other genomic features
found within 50 Kpb range upwards and downwards.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-2/Figure-2}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-3/Figure-3}

\end{center}
\end{figure}

\selectlanguage{english}
\FloatBarrier
\end{document}

