Grassland expansions promoted global diversification of the Pardosa wolf spiders during the late Cenozoic (Araneae, Lycosidae)

The spiders in the genus Pardosa C.L. Koch, 1847, are a young lineage of the family Lycosidae Sundevall, 1833, that exhibit high species diversity and widespread distribution. Pardosa is abundant in open and disturbed environments. In fact, most of its species live in grasslands, and the few that live in forests switched habitats relatively recently. The genus markedly prefers grasslands with a broad range of climates. Thus, its origin and diversification were probably associated with grassland expansions during the late Cenozoic. To test this hypothesis, we developed a global phylogenetic hypothesis that helps reconstruct the biogeographic patterns of the genus Pardosa using three nuclear (18S, ITS2, and H3 ) and four mitochondrial (12S, 16S, NADH1 , and COI ) loci. Our phylo - genetic analyses cover 133 (125 described and 8 as yet undescribed) grassland species of Pardosa using Trochosa ruricola (De Geer, 1778) and Lycosa coelestis L. Koch, 1878, as outgroups. The results show that our selection of species in the genus is divided into four major clades: Clade I includes only P. crassipalpis Purcell, 1903, from South Africa; Clade II consists of a north-east African group (2 species) and a south-east Asian group (21 species); Clade III contains only P. sutherlandi (Gravely, 1924) from SE Asia; and Clade IV includes five species groups from Asia, Europe, and the Americas. The spiders of the genus probably originated in southern Africa or southern and eastern (SE) Asia at the Middle Miocene Climatic Optimum, about 19.40–14.18 Ma, and then ex - panded northwards to North America via the Bering Strait, as well as southwards to north-east Africa via the Arabian Peninsula, and westwards to Europe via western Asia between about 10.59 and 5.28 Ma. At least three exchanges occurred between North America and SE Asia, and at least two between Europe and North America. The biogeography of Pardosa in the past 14.18 Ma, associated with the evolution of grasses, suggested a late Cenozoic diversification of the genus as grasslands expanded.


Introduction
Plants in the family Poaceae, usually called grasses, are composed of more than 12 thousand species (Govaerts et al. 2021).They originated during the late Cretaceous (about 100 million years ago, Ma; Gallaher et al. 2022) and currently cover over 25% of all land on Earth, com-prising 35% of the subtropics and tropics.Grasses also occur in all terrestrial habitats with a broad range of climates, from cold to hot and arid to wet.Climate, soils, fire, and herbivory shaped global grassy systems, including the C3 and C4 photosynthesis groups (Linder et al. 2018;Schubert et al. 2019;Strömberg and Staver 2022).The C3 grasses developed the tolerance needed to survive in colder climates (Schubert et al. 2019); the C4 grasses, evolutionarily derived from the C3 grasses, can prosper in hot and dry areas (Linder et al. 2018).Grassy biomes (>20 Ma; Strömberg and Staver 2022), associated with open-canopy habitats, harbor many spiders of the family Lycosidae Sundevall, 1833 (wolf spiders; Piacentini and Ramírez 2019).They interact ecologically with a diverse part of the lycosid fauna, both above and below ground; genera with representatives occurring in grasslands include Pardosa C. L. Koch, 1847;Lycosa Latreille, 1804, Trochosa C. L. Koch, 1847;and Hippasa Simon, 1885.These lycosids rely on native grasses as a substrate to live on and to forage on grass-eating insects.
Based on the remarkable abundance of wolf spiders in open habitats and in the fossil record (Wunderlich 2004a), the retention of the third tarsal claw, and the conservative morphology, Jocqué and Alderweireldt (2005) proposed that Lycosidae had a relatively recent origin and had co-evolved with grasslands and, by extension, all types of open habitats.Recently, a dated phylogenetic analysis (Piacentini and Ramírez 2019) inferred that the family appeared before the expansion of grasslands and diversified (about 50 Ma) with the reduction of tropical forests and the advance of open habitats.However, the linkages between the temporal-spatial diversification of Lycosidae, especially its most diverse genus Pardosa (based on the number of described species), and grassland expansions have not been specifically tested.
Pardosa wolf spiders inhabit nearly all terrestrial habitats worldwide and currently comprise 532 species (World Spider Catalog 2024).They are free-roaming predators that can disperse over short distances; the spiderlings of some species have been reported to disperse by ballooning on silken threads (e.g., Richter 1970;Greenstone 1982).They are abundant in grasslands and all types of open habitats, such as wetlands, stream banks, floodplains, glades, open deserts, farmlands, and human settlements (Jocqué and Alderweireldt 2005).The marked preference for grasslands contributed to their diversification and successful expansions when the arid and grasslands extended over the globe.Previous studies indicated that Pardosa diversified between about 14 and 10 Ma (Piacentini and Ramírez 2019), and most of its species are remarkably distributed in the four areas with different grassland environments and grassy (C3 and C4) evolutionary histories: (1) North America, NA; (2) eastern and southern Africa, EA; (3) central, eastern, and southern Asia, CA; and (4) the region surrounding the Mediterranean Sea, RM. (Fig. 1; Suppl.material 1: table S1; e.g., Vogel 1964Vogel , 1970Vogel , 2004;;Hänggi et al. 1995;Ivanov et al. 2023).Therefore, Pardosa is a particularly appropriate animal group to test the hypothesis of a diversification of specific lycosid lineages as grasslands expanded.However, the genus lacks a comprehensive phylogenetic hypothesis, and thus its lineage compositions, phylogenetic relationships, and geographical distribution remain unclear.
Spiders of the genus Pardosa have been proposed as effective biological control agents for pests in agricultural systems globally.As the biodiversity and ecological prominence of Pardosa spiders are increasingly recognized, there is a demand for illumination of their biogeographical patterns and diversification mechanisms for the purposes of protection and ecological management.The purpose of this study was to produce a robust phylogenetic hypothesis to elucidate the relationship of Pardosa grassland lineages based on a global sample of species using three nuclear (18S, ITS2, and H3) and four mitochondrial (12S, 16S, NADH1, and COI) loci.Furthermore, we explored the origin, diversification timeline, and global expansion history using a dated phylogenetic tree and tested the co-evolution between Pardosa spiders and grasslands.

Phylogenetic analyses
Phylogenetic relationships among the sampled Pardosa species were inferred using both the maximum likelihood and Bayesian inference approaches.The wolf spiders Trochosa ruricola (De Geer, 1778) and Lycosa coelestis L. Koch, 1878, were used as outgroups.Maximum likelihood analyses were implemented using the fast online phylogenetic tool W-IQ-TREE (Trifinopoulos et al. 2016).The optimal substitution model for each gene partition (TIM2+F+I+G4 for 12S, 16S, NADH1, and COI; K2P+I for 18S; K2P+R2 for H3; and K2P+I for ITS2) was estimated simultaneously using the greedy algorithm in Mod-elFinder (Kalyaanamoorthy et al. 2017) with the Bayesian information criterion and the FreeRate heterogeneity.We set the perturbation strength (p) and the number of iterations since the last best tree was found (c) to 0.3 and 1000, respectively.The SH-aLRT (Guindon et al. 2010) and the ultrafast bootstrap (UFBoot; Minh et al. 2013) were used to assess the support of the branching patterns estimated in the phylogeny with 0.99 of the minimum correlation coefficient and 1,000 of the maximum number of iterations.Bayesian analyses were performed with MrBayes 3.2.1.The best-fitting substitution models were selected by jModelTest under the Bayesian information criterion (Suppl.material 1: table S3; Posada 2008).The Markov Chain Monte Carlo (MCMC) chain was run for 150 million generations using parameters unlinked among partitions and sampled every 100 generations.We used Tracer v1.5 (Rambaut and Drummond 2009) to monitor the mixing of the MCMC chains.A burn-in sample of 375,000 trees was discarded, and a 50% majority rule consensus tree was computed with the remaining trees.

Divergence time estimation
An uncorrelated lognormal relaxed molecular clock model was used to estimate divergence time in BEAST v1.8.1 (Drummond and Rambaut 2007).The birth-death speciation process was chosen as the tree prior.Partitioned strategies (Brandley et al. 2005) were incorporated in BEAST analyses, and each gene was used as a separate partition.For each partition, the specific model of evolution was recommended by jModelTest (GTR+I+G for 12S, 16S, NADH1, and COI; HKY+G for 18S, ITS2, and H3).We ran the MCMC for 50 million generations and sampled every 1000 generations.The maximum clade credibility tree was computed using TreeAnnotator v1.8.0 based on the remaining trees, after discarding the first 25% of the yielded trees as burn-in.Tracer v1.5 was used to determine convergence and measure the effective sample size (>200) for all parameters.We dated the tree of all globally sampled Pardosa species.For the molecular clock analysis, we used the minimum ages based on fossils of Lycosidae The Agelenidae, Thomisidae, Oxyopidae, Psechridae, Trechaleidae, Selenopidae, and some Lycosidae species were used as outgroups.Their gene sequences were available from our study and GenBank (Suppl.material 1: table S1).

Biogeographical reconstruction
The biogeographical history of Pardosa was reconstructed in RASP v3.0 (Yu et al. 2010(Yu et al. , 2015) ) using the Bayesian binary MCMC analysis (BBM; Sanmartín et al. 2001), the statistical dispersal-vicariance analysis (S-DIVA; Yu et al. 2010), and the dispersal-extinction-cladogenesis analysis (DEC; Ree and Smith 2008).We set two areas for the maximum number of ancestral states at each node.The fossil-calibrated trees obtained from BEAST analysis, from which the outgroups were trimmed, were used for biogeographic reconstruction.The distributional data of Pardosa species were available from previous literature, GenBank, and information from our samples.Based on the current distribution and the phylogeny of Pardosa as well as geographic divisions and climate, we defined five geographic areas occupied by the genus: southern Africa, north-east Africa, West Asia+Europe/Palearctic, southern and eastern (SE) Asia/Orient, and North America/Nearctic.diuturnal group (DIG) comprise 7 and 2 Pardosa species from North America, respectively; the paludicola group (PAG) is comprised of 18 species from Europe and North America; the alacris group (ALG) includes 11 species from SE Asia, Europe, and North America; and the brevivulva group (BRG) consists of 95 species from SE Asia, West Asia, Europe, and North America.Bayesian analyses suggest that all sampled African species (P.crassipalpis, P. injucunda, and P. naevia) constitute Clade I with a low posterior probability (only 0.62), and the phylogenetic relationships among most species/groups within Clade IV are unclear (Suppl.material 1: fig.S2).

Divergence time
The fossil-calibrated phylogeny is shown in Fig. 3.The recent ancestor of Pardosa probably appeared during the early Miocene, approximately 19.40 Ma (95% credibility interval, CI: 24.72-14.47Ma).The divergence of the South African P. crassipalpis (Clade I) from the species of the clades (II-IV) was predicted to have occurred around 14.18 Ma (95% CI: 18.19-10.74Ma).Thus, Pardosa wolf spiders originated probably at about 19.40-14.18Ma (the Middle Miocene Climatic Optimum, MMCO).Within Clade II, the Pardosa species from north-east Africa diverged from those from SE Asia around 9.64 Ma.The split between the SE Asian P. sutherlandi (Clade III) and the taxa of Clade IV from Asia, Europe, and the Americas started about 12.59 Ma.Clade IV first diversified during the late Miocene, approximately 8.82 Ma (95% CI: 10.77-8.73Ma), and then split into distinct species groups during the Miocene.The North American diuturna group (DIG) diverged from the nigriceps group (PABG) about 8.17 Ma.Diversification of the widespread species group (PABG) started around 7.77 Ma.

Biogeography
The BBM analysis (Fig. 4) inferred the two possible ancestral ranges occurring during the middle Miocene for Pardosa: (1) southern Africa (the marginal probability, MP: 0.6616); and (2) SE Asia/Orient (MP: 0.2391).The inference from S-DIVA (Suppl.material 1: fig.S3a) favored southern Africa+SE Asia (MP: 0.6534), SE Asia+north-east Africa (MP: 0.2066), and southern Af-rica+North America/Nearctic (MP: 0.0586) as the most likely ancestral areas of the genus, whereas southern Af-rica+SE Asia was preferred under DEC (Suppl.material 1: fig.S3b).The evolutionary routes of Pardosa inferred from the BBM (Fig. 4), S-DIVA (Suppl.material 1: fig.S3a), and DEC (Suppl.material 1: fig.S3b) were generally concordant.The ancestral reconstructions suggest that SE Asia is an important dispersal center for Pardosa, with dispersal from SE Asia to North America, northeast Africa, and Europe at around 10.59, 9.64, and 5.28 Ma, respectively.The three intercolonial evens between South-east Asia and North America were observed during the period of about 10.59-5.31Ma.The inter-colony dispersal between Europe and North America occurred between 5.48 and 2.80 Ma.

Discussion
A detailed exploration of the evolutionary history of Pardosa requires a robust phylogenetic framework (Graybeal 1998;Zwickl and Hillis 2002).Current knowledge of the evolutionary relationships among Pardosa species is based on assessments of morphological similarity or phylogenetic analyses of lycosids, including very few species of Pardosa (e.g., Zehethofer and Sturmbauer 1998;Vink et al. 2002;Muster and Berendonk 2006;Murphy et al. 2006;Piacentini and Ramírez 2019;Roslin et al. 2022).This study first comprehensively addresses the phylogeny of global grassland species of Pardosa using multiple nuclear and mitochondrial markers.The topology of the maximum likelihood tree (Fig. 2) is similar to that of the dated tree (Fig. 3) and provides a well-supported hypothesis for the relationships among the major clades/groups of Pardosa within our dataset.The biogeographic histories of Pardosa suggest that the evolution of this genus was strongly affected by grassland expansions resulting from historic climatic and environmental shifts.These findings mirror other studies highlighting the vital role of the evolution of grasslands in lycosid biogeography (Jocqué and Alderweireldt 2005;Piacentini and Ramírez 2019).
Lycosidae appeared about 33.80 Ma, after the Eocene-Oligocene extinction event, but well before the grassland expansions.Pardosa (Pardosinae) diverged from the clade (Lycosa+Trochosa) around 19.40 Ma, which is consistent with the age (25-16 Ma) obtained in the phylogenetic analysis of Piacentini and Ramírez (2019).The initial diversification of Pardosa spiders was about 14.18 Ma, coinciding with the Middle Miocene Climatic Optimum (17-14 Ma).The time  when Pardosa appeared is much later than the original age (55 Ma; Strömberg 2011) of C4 grasses with the hot and dry tolerance (Linder et al. 2018), and the age (30 Ma) when in colder climates, C3 open-habitat grasses developed the tolerance needed to survive frosts (Schubert et al. 2019), but coincides with the time (the late Cenozoic) when open-habitat grasses began to become ecologically dominant (Strömberg and Staver 2022).Four major clades of Pardosa emerged during the middle to late Miocene, between 14.81 and 8.82 Ma.This pattern coincides with the rapid expansion of open-habitat grasses and the retreat of tropical forests due to the cooling and aridification climate and the low CO 2 conditions, as well as fire and herbivory activities after the MMCO (Zachos et al. 2001;Strömberg and Staver 2022).Our results support the idea that Pardosa is a young lineage of wolf spiders that became vagrant and diversified during the Miocene, when the grasslands expanded.In the past 8 Ma, including the  pattern is likely associated with the continent-specific trajectories of grasses (Karp et al. 2021;Kukla et al. 2022).For example, the subtropical C3 open-habitat grasses first spread to the colder regions (Kukla et al. 2022).In the time (about 10-6 Ma), the tropical open-habitat C4 grasses expanded to form grasslands and savannas at low to mid-latitudes (Karp et al. 2021;Lu et al. 2020), and the frost-tolerant grasses spread to higher latitudes.
During around 10.59 to 5.31 Ma, multiple exchanges (at least three times) for the spiders occurred between North America and South-east Asia.The Pardosa exchanges also appeared between Europe and North America during around 5.48 to 2.80 Ma.These results indicated that intercolony dispersals were associated with grassland extensions and retreats resulting from climatic shifts such as historical glaciations.The strong dispersal capacity and adaptability to complex climates and disturbed environments (Richter 1970;Greenstone 1982;Samu and Szinetár 2002;Jocqué and Alderweireldt 2005;Woolley et al. 2016) led to a world-wide distribution of Pardosa spiders that were dominant in open habitats within only about 15 Ma.Frequent ballooning and cursorial dispersal confer their high mobility and area expansion (e.g., Richter 1970;Greenstone 1982).Our findings indicate that the long-distance overwater dispersal and subsequent range expansions of Pardosa spiders occurred between South-east Asia and southern Africa and between North America and Europe.Furthermore, the marked preference of Pardosa wolf spiders for disturbed habitats, such as clearings in grasslands or forested areas (Jocqué and Alderweireldt 2005), made them especially successful when arid and open habitats, including grasslands, extended over the globe during the late Cenozoic.

Future directions
This study generated hypotheses regarding the origin and dispersal of Pardosa grassland lineages and suggested that grassland expansions drove its global diversification during the late Cenozoic using the nuclear 18S, ITS2, and H3 and mitochondrial 12S, 16S, COI, and NADH1 loci.Sampling is the process of choosing a subset of a target lineage that will serve as its representative.In this paper, the total sampling specimens for grassland species and some regions are on the low side.However, we hope our study can aid in strategic resampling, reflecting known lineage divergences from grasslands.Moreover, the taxonomic revisions of Pardosa need to be made in advance.

Figure 1 .
Figure 1.The major distributional areas of Pardosa spiders and the main biomes across the globe.White spots indicate the probable type localities of all known Pardosa species; red spots indicate the sample localities of Pardosa spiders sequenced in this study.

Figure 2 .
Figure 2. Phylogenetic tree of 133 Pardosa species reconstructed using the maximum likelihood method.The numbers at the nodes represent bootstrap support values from the maximum likelihood analyses.
global Ice Age (2.6 Ma onward), climatic cooling further accelerated, and grasslands extended as the forests retreated and the open and arid environments expanded.At this time, Pardosa spiders rapidly diversified, currently representing the most abundant lycosid genus and being dominant in grasslands and some open environments.Our ancestral reconstructions suggested that the Pardosa spiders originated in southern Africa or SE Asia/ Orient during the middle Miocene, likely evolving concurrently with the spread of grasses and the diversification of herbivores (~17 Ma; Charles-Dominique et al. 2016; Strömberg and Staver 2022), followed by the geographic expansions from middle altitudes of SE Asia to high altitudes, and then to North America via the Bering Strait, as well as to north-east Africa southwards via the Arabian Peninsula, and to Europe westwards via West Asia between about 10.59 and 5.28 Ma.This dispersal

Figure 4 .
Figure 4. Biogeography of Pardosa and potential global dispersion routes (arrows).Reconstruction using Bayesian binary MCMC (BBM) in RASP v3.0.The colors of pie wedges at each node represent geographical areas inferred to have been occupied by ancestral taxa.Pink circles in the maps indicate the possible ancestral ranges occurring during the middle Miocene for Pardosa.The numbers at the nodes represent support.