Boraginales are a member of Lamiidae, a subclass that includes eight orders (Solanales, Lamiales, Vahliales, Gentianales, Boraginales, Garryales, Metteniusales and Icacinales) and 46 families (APG IV, 2016). The classification of this order has changed substantially over time. In the broad sense, Boraginaceae (Boraginaceae s.l. or Boraginales) (i.e. Candolle, 1845; G€ urke, 1897) have included five subfamilies, namely Boraginoideae, Cordioideae, Ehretioideae, Heliotropioideae and, recently, Wellstedioideae. Currently, these subfamilies are largely recognized at the family level. Recently, Luebert et al. (2016) suggested the recognition of 11 monophyletic families: Boraginaceae s.s., Codonaceae, Coldeniaceae, Cordiaceae, Ehretiaceae, Heliotropiaceae, Hoplestigmataceae, Hydrophyllaceae, Lennoaceae, Namaceae and Wellstediaceae.
In the present study, Boraginaceae are used in the strict sense (Boraginaceae s.s.; see Luebert et al., 2016), and the family currently includes c. 90 genera and c. 1600 species distributed across three subfamilies: Echiochiloideae, Boraginoideae and Cynoglossoideae (Chac on et al., 2016). Boraginoideae and Cynoglossoideae are further divided into two and eight tribes, respectively.
Boraginaceae have a wide geographical distribution, occurring in temperate, Mediterranean and alpine climates (Kadereit and Bittrich, 2016). More than half of the genera of the subfamilies and tribes have their centre of diversity in Eurasia. The centres of diversity for Echiochiloideae, Boragineae and Lithospermeae are in Western Asia and the Mediterranean region (Kadereit and Bittrich, 2016). The centre of diversity of Cynoglossoideae is East and Central Asia. The Americas, Africa and Australia have representatives of the tribes, but Africa and the Americas have only three (Afrotysonia Rauschert, Echiostachys Levyns and Lobostemon Lehmn.) and seven (Amsinckia Lehm., Antiphytum DC. ex Meisn., Cryptantha Lehm. ex G. Don, Dasynotus I.M. Johnst., Moritzia DC. ex Meisn., Pectocarya DC. ex Meisn. and Thaumatocaryon Baill.) endemic genera, respectively. Australia includes only one endemic genus (Embadium J.M. Black). Some genera are cosmopolitan and present on five continents (Hackelia Opiz, Lappula Moench) or even six continents (Cynoglossum L. and Myosotis L.) (Kadereit and Bittrich, 2016).
Boraginaceae, along with Gentianaceae (Ferreira de Sousa et al., 2019), Saxifragaceae (Perveen and Qaiser, 2009) and other families, are eurypalynous (i.e. characterized by substantial variation in pollen morphology) (Halbritter et al., 2018). Palynological studies in Boraginaceae have a long history. Indeed, the first studies on the pollen morphology of some members of Boraginaceae were conducted by Fritzsche (1832) and Mohl (1835), who characterized the structure and number of apertures in many genera (Sahay, 1979). The first comprehensive studies on the pollen morphology of Boraginales were by Erdtman (1952Erdtman ( , 1966) ) and Avetisian (1956). Erdtman (1966) demonstrated similarities between pollen morphology of Ehertiaceae and some species of Hydrophyllaceae, Lennoaceae, Verbenaceae and Hoplestigmataceae. Avetisian (1956) studied pollen morphology of 216 species occurring in the former Soviet Union, and he recognized pollen traits of Cordiaceae and Ehertiaceae. Clarke (1977) employed electron microscopy to study pollen of Boraginales from northwest Europe, and using five qualitative and six quantitative features recognized 17 pollen types. Sahay (1979) used scanning electron microscopy (SEM) and transmission electron microscopy (TEM) to compare pollen of Boraginales with allied families and demonstrated that Boraginaceae, with 28 pollen morphoforms, are a highly eurypalynous family. Pollen morphology of Boraginales from Korea and Pakistan was investigated by Ahn and Lee (1986) and Perveen et al. (1995), respectively. The studies emphasize the role of pollen morphology in the classification of tribes and subfamilies, and intergeneric relationships within the order.
Starting in the early 1990s, multiple surveys using electron and light microscopy were conducted on tribes and subtribes of the family. One of the most comprehensive of these was undertaken by Bigazzi and Selvi (1998) who recognized 15 pollen types based on wide variation in pollen shape, exine ornamentation and number of apertures. Liu et al. (2010) investigated the pollen morphology of five native genera of Lithospermeae in China, and proposed an identification key for genera of the tribe. Studies of Cynoglosseae (Mathez, 1974;Liu et al., 2001;Fukuda and Ikeda, 2012;Attar et al., 2018) and Eritichieae (D ıez and Vald es, 1991) provided evidence that the pollen of these two tribes is heterocolpate, whereas in Lithospermeae and Boragineae, pollen is colporate. In addition, studies on single genera or species also were undertaken using SEM and TEM (i.e. Hargrove and Simpson, 2003;Bigazzi et al., 2006;Coutinho et al., 2012). These studies were valuable in delimiting tribes and genera, and also in evolutionary studies across the family. However, studies have focused mostly on specific clades, and despite the large amount of pollen study throughout the family, no comprehensive examination of variation across the entire family has been undertaken.
In recent decades, some researchers have conducted phylogenetic studies focused on relationships within genera and tribes, or among closely related genera and tribes (i.e. Selvi et al., 2006;Cohen, 2011;Simpson et al., 2017;Mabry and Simpson, 2018). Recently, Nazaire and Hufford (2012), Cohen (2014) and Chac on et al. ( 2016) reconstructed more comprehensive phylogenies of the family. At the same time, there have been minimal efforts to employ these molecular phylogenies to examine pollen evolution.
Recent analytical advancements have facilitated our abilities to collectively explore phylogenetics, character evolution and rates of evolutions (Kriebel et al., 2017). The evolutionary changes and adaptive significance of pollen characters have been debated for decades throughout the angiosperms. General evolutionary trends in angiosperm pollen include: increasing pollen grain size, shifting aperture position from equatorial to global, changing aperture shape from colpate to porate, and increasing complexity of tectum ornamentation (Zhang et al., 2017).
Several hypotheses concerning adaptive evolution of angiosperm pollen characters have been proposed, and possible relationships have been observed between pollination strategies and pollen features (Muller, 1979;Linder, 2000;Zhang et al., 2017). Regarding tectum ornamentation, for example, entomophilous flowers often appear to produce pollen grains with an ornamented exine surface, whereas psilate pollen may be characteristic of anemophilous and hydrophilous plants (Hu et al., 2008). The relationship of pollen size and pollination system also has been repeatedly discussed. For instance, Muller (1979) suggested that wind-pollinated pollen grains range between 20 and 40 µm in size because larger pollen is transported mostly by animal pollinators and very small pollen often is associated with highly specialized pollination systems. Recently Zhang et al. (2017) confirmed Muller's results and suggested a possible hypothesis that different-sized pollen grains correspond to different pollination strategies.
The evolutionary trends of pollen morphology in some taxa of the Lamiids have been discussed in some studies: an evolutionary shift in exine sculpture, transitions from retipilate sculpturing to retirugulate and reticulate sculpturing, has been described for the tribe Rhinantheae (Orobanchaceae) (Lu et al., 2007); an increase in aperture number, transitions from zonoapertures to pantoapertures, and from colpi to pores have been inferred for American species of Borreria G. Mey. (Rubiaceae) (Pire, 1996); and evolutionary trends including an increase in aperture number, and a transition from imperforate to perforate to microreticulate tectum sculpture has been found in Cuscuta L. (Convolvulaceae) (Welsh et al., 2010).
Boraginaceae exhibit a great array of specialized pollen types and, thus, are an interesting plant group for investigating the evolution of the pollen characters and their possible role in the diversification of the group. Recent, increasingly robust phylogenies of Boraginaceae provide a foundation for investigating the distribution of pollen features. Indeed, investigating pollen characters in a molecular phylogenetic framework may provide some clues to understanding the distribution and evolutionary patterns of pollen features, and also identify the possible morphological characters associated with shifts in diversification rates of the family.
The present study has two main goals. The first is to examine pollen diversity across Boraginaceae using data from published papers and our own work, incorporating data from >500 species. Special attention was paid to Lithospermeae; previous studies (e.g. Liu et al., 2010;Cohen, 2014) have pointed to high diversity of pollen grain in Lithospermeae, yet the palynological knowledge of the tribe is still fragmentary and pollen of some genera of the tribe, such as Ancistrocarya Maxim., remain poorly examined. The second goal is to resolve patterns of evolution of multiple features of pollen, such as size, shape, aperture number and rate of diversification, and to examine the potential association between diversification of species and pollen size.
Forty-three flowers belonging to 43 species representing 25 genera of Lithospermeae were obtained from herbarium specimens in the Missouri Botanical Garden (MO), Makino Botanical Garden (MBK), Haettepe University (HUB) and Conservatoire et Jardin botaniques de la Ville de Gen eve (G). A list of species, voucher numbers and locations is presented in Table S1.
Staminate flowers were taken from herbarium specimens and immersed in Pohl's solution for 24 h to soften and hydrate tissues. One to a few anthers from each sample were transferred to small sheets of folded filter paper, placed in sample bottles, then dehydrated using a graded ethanol series consisting of 80%, 90% and 95% EtOH, and subsequently soaked in acidified 2,2dimethoxypropane (DMP, MilliporeSigma, St Louis, Missouri, USA) to remove residual water (De Moura et al., 2018). The packets were then placed in acetone and critical point-dried in a K850 critical point dryer (Electron Microscopy Services, Hatfield, Pennsylvania, USA). Pollen samples were placed on stubs and sputtercoated with gold (Denton Vacuum, Moorestown, California, USA) on a Denton Desk V Sputter Coater (Denton Vacuum), at 35 mA for 120 s, yielding a thickness of c. 71.4 nm. The pollen grains were photographed with a NeoScope JCM-5000 scanning electron microscope (JEOL USA Inc., Peabody, Massachussets, USA) using 15 kV accelerating voltage and operating under high vacuum. Detailed qualitative features, including pollen shape in polar and equatorial view and aperture type, and quantitative features, including pollen size (polar and equatorial diameters) and polar:equatorial ratio, were determined (Table S2). For quantitative characters, ≥15 pollen grains were measured for each sample. The terminology followed the "Illustrated Pollen Terminology" (Halbritter et al., 2018) and Punt et al. (2007).
Based on the variation observed across the family from data from the literature and our own survey, a set of seven palynological characters were coded for 561 taxa belonging to 538 species from three subfamilies of Boraginaceae and relatives in closely related families, including Codonaceae, Wellstediaceae, Hydrophyllaceae, Heliotropiaceae, Ehretiaceae and Cordiaceae. Data were collected for polar (Paxis) and equatorial (E-axis) diameters, aperture type (Heteroaperture or Iso-aperture), aperture number, pollen shape and pollen outline. P-and E-axes data were gathered from previous and current studies and in cases where these values were not mentioned in the papers, photos of pollen grains from the literature were measured using ImageJ2 software (Rueden et al., 2017). Prism v.8.0 (Graph-Pad Software, San Diego, California USA; www.graphpad.com) was employed for statistical analyses among tribes for the P-and E-axes and P:E ratio, using a one-way ANOVA followed by Tukey's multiple comparisons test. Violin and box plots were constructed from the P-and E-axes, P:E ratio, pollen shape, pollen outline, pollen type and aperture number.
A total of 303 species from 71 genera of Boraginaceae plus seven genera from related families of Boraginales, including Hydrophyllaceae, Heliotropiaceae, Ehretiaceae, Cordiaceae, Wellstediaceae and Codonaceae (in total 1390 sequences), were gathered using the NCBI GenBank portal (accessed 3 April 2019) (Table S3). Seventeen chloroplast DNA regions,including trnL,rps16,ndhF,rpl16,trnK,rbcL, and atpB, and the nuclear ribosomal (nrDNA) internal transcribed spacer (ITS) and external transcribed spacer (ETS) were used for phylogenetic analysis. Sequences were aligned separately for each gene using MAFFT v.7 (Katoh et al., 2019), under default parameters (FFT-NS-1 strategy), then were adjusted manually, and concatenated for subsequent phylogenetic analyses. For the molecular matrix, c. 70% of sequences were missing. This high number is the result of sequences obtained from available data deposited in GenBank.
Phylogenetic analyses were conducted with the total evidence molecular matrix that includes all sequence data (cpDNA and nrDNA). Parsimony (MP) phylogenetic analyses were conducted in TNT v.1.5 using the following parameters: 1000 ratchet iterations, with a 10% probability of up-and down-weighting, 100 rounds of tree fusing, 1000 drift cycles and random sectorial searches (Goloboff and Catalano, 2016). Ten thousand bootstrap replicates were implemented in TNT using the same search parameters as above and a 37% deletion probability.
Maximum-likelihood (ML) analyses were conducted using RAxML-HPC on XSEDE (v.8.2.12) at CIPRES Science Gateway (Miller et al., 2010) using the GTR + I + G substitution model as suggested by jModelTest2 on XSEDE (v.2.1.6) according to the Akaike Information Criterion (AIC) at the CIPRES portal. Rapid bootstrap analyses (1000 replicates) and searching for the bestscoring ML tree occurred simultaneously (the "-f a" option).
Bayesian inference (BI) analyses were undertaken with BEAST 2 v.2.6.1 (Bouckaert et al., 2019) using the Kettering University High-Performance Computing Cluster (KUHPC). The Yule model, the relaxed clock model, and the GTR + I + G substitution model were employed. Four Boraginaceae fossils were used that were applied to stem nodes of (i) Cryptantha (Cynoglosseae) and (ii) Echiochilon Desf. (Echiochileae), and the crown nodes (iii) Lappula (Rochelieae) and (iv) Lithospermum L. (Lithospermeae) which were all set with uniform priors. This calibration scheme is based on Chac on et al. (2017), with details on the exact priors available in Table S4. Parameters were estimated using 200 million generations with trees sampled every 20 000 generation. Tracer v.1.7.1 (Rambaut et al., 2018) was employed to examine the sampling adequacy and convergence of the chains to a stationary distribution. Finally, TreeAnnotator v.2.6.4 was used (Drummond and Rambaut, 2007) to summarize the postburn-in trees (burn-in was set at 80%) and produce a maximum clade credibility (MCC) tree.
In order to estimate the ancestral value of the P-axis, E-axis and P:E ratio in Boraginaceae and related families, three independent analyses were conducted with four different coding schemes by using two approaches, MP and ML, to infer the patterns of pollen evolution. The minimum (MIN) and maximum (MAX) values of characters for each species as well as the mean (MEAN) and range (RANGE) were calculated as four different coding schemes. Parsimony ancestral state reconstruction was performed in two different software programs, TNT v.1.5 (Goloboff andCatalano, 2016) andMESQUITE v.3.6 (Maddison andMaddison, 2019). These analyses were conducted with the combined matrix that included all molecular and morphological data. In MESQUITE, ancestral values for MIN, MAX and MEAN of continuous characters (P-axis, E-axis and P:E ratio) were reconstructed using squared-change parsimony (Maddison, 1991). MIN, MAX, MEAN and RANGE values were conducted with TNT, with the aforementioned strategy. In addition, to evaluate the detailed evolutionary trend of pollen size variation along the phylogeny, we used a node-by-node approach in the framework proposed by Gould and Macfadden (2004). We interpreted nodal location and magnitude of change for each node.
The ML ancestral state reconstruction of continuous characters for MIN, MAX, MEAN and RANGE was investigated with the fast Anc function in the R package Phytools (Revell, 2012). Using these data, the magnitude of change for each node (increase or decrease) was interpreted. In addition, the ancestral state reconstruction for continuous traits was visualized with the contMap function in Phytools.
Ancestral state reconstructions were estimated for pollen shape, outline-in-equatorial-view, and aperture number and type at each node using multiple software packages, including TNT (Goloboff and Catalano, 2016), MESQUITE (Maddison and Maddison, 2019) and the R package corHMM (Beaulieu and Donoghue, 2013), in MP and ML frameworks. Parsimony ancestral state reconstructions of quantitative characters were performed using TNT with the same strategy as described previously, and in MESQUITE, an unordered parsimony model was applied. In corHMM, the rayDISC function, which accommodates polymorphisms, missing data and multi-state characters, was employed. Character optimizations using all rates different (ARD) model and marginal reconstruction of ancestral states were explored for these data, and all remaining parameters were set to their default values.
In order to model patterns of diversification in the family and to identify shifts in P-axis, E-axis and P:E ratio, Bayesian Analysis of Macroevolutionary Mixtures (BAMM v.2.5.0) (Rabosky, 2014) was used, which models complex dynamics of speciation, extinction, and trait evolution on an ultrametric phylogenetic tree. BAMM model shifts in macroevolutionary regimes across a MCC tree were performed using reversible-jump Markov chain Monte Carlo (rjMCMC) sampling. Initial prior settings were calculated with the setBAMMpriors function in the BAMMtools package in R (Rabosky et al., 2014;Rabosky, 2014), and the expectedNumberOfShifts parameter was set at 1.0. Because of incomplete sampling, we provided the program with a file containing sample fractions (one for genera, and one for tribes) using the sampleProbsFilename argument (Table S5). BAMM was run for 100 million generations on the MCC tree, sampling parameters every 1000 generations. The BAMM outputs were used to perform multiple analyses with R/BAMMtools. For both analyses (diversification and pollen size evolution), the 95% credible set of rate shift configurations was obtained with the function credibleShiftSet. Single best shifts configuration and the cumulative shift probability tree were conducted using the getBestShiftConfiguration and cumulativeShiftProbsTree functions, respectively. Finally, rates through time were computed for each tribe using the getRateThroughTimeMatrix function and visualized by plotRateThroughTime.
The correlation between the three continuous characters (P-and E-axes, and P:E ratio) and macroevolutionary parameters (Speciation (k), Extinction (l) and Net diversification (r)) based on Pearson, Spearman and Kruskal methods was estimated using Structured Rate Permutations on Phylogenies (STRAPP) (Rabosky and Huang, 2016) in BAMMtools. Results from BAMM were compared with those from ParSplit, an MP approach to evolutionary shifts, which is useful given controversy surrounding BAMM (Moore et al., 2016;Rabosky et al., 2017;Meyer and Wiens, 2018;Meyer et al., 2018;Rabosky, 2019). ParSplit v.3 (Chabrol and Didier, 2017) was used for performing the parsimonious split of a phylogenetic tree with regard to tip values of quantitative characters (P-and E-axes, and P: E ratio).
Nineteen DNA regions were included in the current analyses (Table S6). The 17 cpDNA regions comprised a total of 22 747 aligned nucleotides, 3887 of which were parsimony-informative. The nrDNA ITS and ETS yielded 1796 aligned nucleotides, and 751 were parsimony-informative. In total, the concatenated molecular matrix included 24 543 aligned nucleotides, with 4638 parsimony informative sites.
In all analyses, Boraginaceae were resolved as monophyletic. The tribes Boragineae and Lithospermeae were recovered as sister group with 99% bootstrap (BS) and posterior probability (PP) of 1.0. According to the MP tree, in Lithospermeae, there was weak support for the clade comprising Echium L. + Lobostemon + Echiostachys + Pontechium B€ ohle & Hilger + Onosma L. + Maharanga DC. (< 50% MP-BS) (Fig. S1); however, both ML and BI analyses provided greater support for this clade (ML-BS = 91, PP = 1.0). Apart from different support values, the consensus tree of each analysis resolved distinct relationships within the closely related species. For example, in ML analyses Echiostachys incanus (Thunb.) Levyns is sister to Lobostemon trigonus H. Buek + Lobostemon montanus H. Buek, but in the BI analysis, E. incanus is only sister to Lobostemon montanus. Based on the resulting trees from the ML and BI analyses, most clades along the backbone of Cynoglossoideae were well-supported (> 96% BS, PP = 1.0), but in MP analyses, some clades, such as Myosotideae, received low support (< 50% BS). In Cynoglossoideae, Trichodesmeae is sister to the rest of the members of the subfamily (ML-BS = 100%, MP-BS = 64%, PP = 1.0). The remaining Cynoglossoideae members were grouped into five larger clades representing tribes. Omphalodeae was sister to Asperugeae (ML-BS = 96%, PP = 1.0, MP-BS = 90%), and these two were sister to the rest of the clades. In Rochelieae, based on ML and BI trees, Lappula was nonmonophyletic and sister to Hackelia. These two genera were sister to a clade consisting of Eritrichium Schrad. + Rochelia Rchb.; however, in MP analyses, Rochelia was sister to Lappula, but this relationship received weak support (MP-BS < 50%). In Myosotideae, Brachybotrys Maxim. ex Oliv. was sister to Trigonotis Steven (ML-BS = 100, MP-BS = 94, PP = 1.0), and these two genera were sister to a clade that includes Myosotis + Decalepidanthus Riedl (ML-BS = 97, MP-BS = 69, PP = 1.0). Cynoglosseae was well-supported (ML-BS = 99, MP-BS = 99, PP = 1.0) and divided into four main clades. In all analyses, these four clades were well-resolved, but the clade composed of Cynoglossinae has little resolution and varying degrees of support (Fig. S1).
Pollen size. In Boraginaceae, the P-and E-axes were in the ranges 4.84-56 µm and 1.68-44 µm, respectively (Table S2). The smallest and largest pollen grains were found in Cynoglossoideae and Boraginoideae, respectively. In Cynoglossoideae, the ranges of the Pand E-axes were estimated as 3.5-41.26 µm and 1.68-25.60 µm, respectively. Among members of Cynoglossoideae, the smallest and largest values for the P-and E-axes were found in Cynoglosseae. In Boraginoideae, the ranges of the P-and E-axes were estimated at 9-52 µm and 4.63-42 µm, respectively. Among members of Boraginoideae, the smallest values for the P-and E-axes were observed in Lithospermeae, and the largest values for these axes were found in the Boragineae. For both axes, greater variation was observed in Boragineae than in Lithospermeae. In Echiochiloideae, the ranges of of the P-and E-axes were estimated at 10.5-17.5 µm and 8.7-13.3 µm, respectively (Fig. 1). The P:E ratio within members of Boraginaceae ranged between 0.63 and 3.8, and the greatest variability was present in Cynoglossoideae. The P:E ratio ranged between 0.99 and 2.45 in Boraginoideae. Among the tribes of Boraginoideae, Lithospermeae had greater variation than Boragineae in P:E ratio. In Echiochiloideae, the P:E ratio ranged between 0.86 and 1.40 (Fig. 1).
Results of the ANOVA-Tukey's multiple comparisons test demonstrated that within Boraginoideae, Boragineae (n = 164) and Lithospermeae (n = 350) were statistically significant (P < 0.0001) and the mean of the P-axis of Boragineae (mean = 33.76 µm) was larger than that in Lithospermeae (mean = 19.07 µm). However, within Cynoglossoideae, the primary statistically significant differences were between Trichodesmeae (n = 27) and other members of Cynoglossoideae (P < 0.0001), with the mean of the Paxis in Trichodesmeae larger than those in rest of the members of the subfamily. Additionally, the P-axis of Asperugeae (n = 124) was larger than that of either Cynoglosseae (n = 235) or Omphalodeae (n = 74) (P < 0.0001) (Table S7). The results for the E-axis were the same as for the P-axis, with few exceptions (Table S7). No significant differences were identified within the P:E ratio of Boraginoideae (i.e. between Boragineae and Lithospermeae) (P = 0.9); however, the P:E ratio of Boraginoideae was smaller than that of all tribes of Cynoglossoideae (P < 0.0001), except Trichodesmeae (P = 0.99) (Table S7).
Pollen shape. The shape of the pollen grains, based on the P:E ratios, was divided into various shape classes. In Boraginaceae, four pollen shapes were identified: spheroidal, subprolate, prolate and perprolate.
In both Cynoglossoideae and Boraginoideae, all four shapes were observed. The pollen grains of Echiochiloideae were subprolate and spheroidal for the five species examined (Fig. 2). Overall, the largest and the smallest percentages of prolate and perprolate shapes were observed in Asperugeae and Omphalodeae, respectively. In these two tribes, both spheroidal and subprolate shapes were absent. The prolate shape was absent only in Echiochiloideae. Perprolate pollen was not observed in Echiochiloideae, Boragineae or Trichodesmeae. The largest and the smallest percentages of subprolate shapes were identified in Echiochiloideae and Rochelieae, respectively. The largest and the smallest percentages of spheroidal pollen were observed in Trichodesmeae and Myosotideae, respectively (Fig. 2). Spheroidal pollen was absent in three tribes: Asperugeae, Omphalodeae and Rochelieae.
Pollen outline. According to the pollen outline-inequatorial-view, seven shapes--elliptical, elliptical with a slightly constricted equator, elliptical-rectangular, rectangular, spheroidal, dumbbell-shaped and ovoid--were recognized in Boraginaceae (Fig. 2). Greater variation in pollen outline-in-equatorial-view was found in Boraginoideae compared with Cynoglossoideae. Within Cynoglossoideae, the largest percentage of pollen outlines was dumbbell-shaped, which was observed in all tribes except Trichodesmeae. Within Boraginoideae, the greatest diversity was recognized in Lithospermeae (Fig. 3), the only tribe that includes all seven pollen outlines shapes. Spheroidal pollen was identified only in Trichodesmeae and Lithospermeae. Rectangular pollen was known from four tribes: Echiochileae, Boragineae, Lithospermeae and Cynoglosseae (Fig. 2).
Pollen aperture. In Boraginaceae, two main aperture types were recognized: isoapertures and heteroapertures. Heteroaperturate pollen was observed in all tribes of Cynoglossoideae except Trichodesmeae, and isoaperturate pollen was identified in the other two subfamilies, Boraginoideae and Echiochiloideae, as well as in Trichodesmeae (Fig. 4). In the family, the most frequently observed aperture number was three. Within all members of Cynoglossoideae, only three apertures and three pseudoapertures were observed. Greater diversity in pollen aperture number was observed in Boraginoideae than in Cynoglossoideae (Fig. 4).
The MIN, MAX, MEAN and RANGE values for quantitative characters (P-axis, E-axis and P:E ratio) were reconstructed in TNT, Mesquite and Phytools using MP and ML approaches. According to the MP and ML results, multiple shifts in length, for both the P-and E-axes, were observed on the tree. One of the main shifts occurred in Boraginoideae with the ancestral value of Boragineae increasing compared to that of Lithospermeae (Table 1). Within Boragineae, at least five increases and two decreases in pollen size were observed. Within Lithospermeae, at least four increases and two decreases were resolved (Fig. 5).
Another shift was in Cynoglossoideae with the ancestral value of Trichodesmeae increasing relative to the rest of the subfamily (Table 1). Other shifts occurred primarily at the origin of various taxa, such as the largest increase at the crown of Amsinckia (Figs 5 and S3). The ancestral values of the P:E ratio revealed shifts between the subfamilies. However, no notable changes occurred between Boragineae and Lithospermeae. Unlike Boraginoideae, the range of ancestral values was greater in Cynoglossoideae and showed remarkable changes among the tribes (Table 1). For example, a shift was reconstructed at the crown of Trichodesmeae, and this shift differs from other members of the Cynoglossoideae (Fig. S3).
The evolution of quantitative characters, pollen shape, outline-in-equatorial-view and aperture number, was reconstructed using three approaches (i.e. in TNT, corHMM and MESQUITE).
Pollen shape. According to all analyses, prolate and/ or subprolate pollen shapes were the ancestral states for Boraginaceae (Figs S2 and S4,and Table 2). The prolate shape was resolved as ancestral for Boraginoideae and both tribes of the subfamily, and this pollen shape has been lost multiple times within the taxon. Based on the ML methodology, the perprolate shape had the greatest probability of being the ancestral condition for Cynoglossoideae; however, MP methods revealed that the prolate pollen shape was the ancestral condition for the subfamily (Table 2). The perprolate pollen shape evolved independently twice across Boraginoideae and more than five times within Cynoglossoideae. Spheroidal and subprolate pollen shapes evolved independently at least five and eight times within the family, respectively. The ancestral state of pollen shape in Trichodesmeae was more similar to that in Lithospermeae and Boragineae than to those from Cynoglossoideae (Figs 6 and S2, and Table 2). According to the MP methods, the ancestor of Echiochiloideae was subprolate (Figs 6 and S2, and Table 2).
Outline-in-equatorial-view. According to results from the MP methods, the ancestor of Boraginaceae had pollen with elliptical and/or elliptical-rectangular outlines (Table 2). The elliptical-rectangular and/or elliptical outlines were resolved as ancestral for Boraginoideae as well as both tribes in it. The rectangular outline was the ancestral state for Echiochiloideae and evolved independently at least four times within the ingroup (Figs 6 and S2). The results from both MP and ML methods revealed that, at the origin of sampled species of the clades of Omphalodeae and Rochelieae, dumbbell-shaped pollen was the ancestral condition, and this type of outline characterizes almost all species in each clade (Figs 6,S2 and S4). This type of pollen also evolved once in Boraginoideae. The ancestral state of the outline was elliptical in Trichodesmeae (Figs 6 and S2, and Table 2). Ovoid pollen originated twice across Lithospermeae--once at the origin of Echiinae and the other in the clade that includes Alkanna Tausch (Figs 6 and S2). Spheroidal pollen arose independently at least five times within the ingroup (Figs 6 and S2). Pollen aperture. According to results from both the MP and ML methods, aperture number was resolved as two to three at the origin of Boraginaceae (Fig. S2 and Table 2). Three apertures originated twice in Lithospermeae--once in Alkanna and once in Echiinae (Figs 6 and S4). The presence of two to three apertures was reconstructed as the ancestral state of Cynoglossoideae and characterizes most species in the subfamily. Four to six apertures was resolved as ancestral for Boragineae, and this character state evolved independently twice within Lithospermeae (Fig. 6). Pollen with more than ten apertures arose four times across the Boraginoideae: once each in Glandora D.C. Thomas, Weigend & Hilger, Moltkia Lehm., Symphytum L. and Nonea Medik. This character was restricted to the members of Boraginoideae. Seven to nine apertures evolved multiple times across Lithospermeae and at least twice within Myosotideae (Figs 6 and S4).
Diversification analyses. Based on the results of BAMM, the best shift configuration had three rate shifts: the origin of Boraginaceae, in Plagiobothrys Fisch. & C.A. Mey., and in Cynoglossinae (Fig. S5). This configuration was also supported by marginal shifts. The 95% credible set of shift configurations also showed a shift in diversification close to the common ancestor of Cynoglosseae and one accelerated shift within Plagiobothrys (Fig. S6). According to BAMM estimates of the speciation rate, these rates tended to increase in Boraginaceae, especially in Plagiobothrys. Speciation-through-time plots for each tribe also provided evidence that speciation rates were mostly constant in all tribes. Only in Cynoglosseae did rates tend to increase, starting c. 10 Ma (Fig. S7).
Pollen size diversification. Results of the best shift configuration analyses in BAMM recognized nine shifts in diversification rates for the P-axis in Boraginaceae: the stem nodes of Boragineae and of Omphalodeae, and multiple shifts within Lithospermeae and Cynoglosseae, (Fig. S8). BAMM analyses for the E-axis indicated ten shifts across the family. The position of these shifts were almost the same as for the P-axis, with only minor differences (Fig. S8). For example, shifts in Buglossoides Moench and the stem node of Trichodesmeae were restricted to the E-axis (Fig. S8). Finally, according to the P:E ratio, 11 shifts were recovered. Most of the shifts occurred in small clades or genera, such as shifts in Microula Benth. and in Plagiobothrys. Shifts in Asperugeae and Boragineae, respectively, were the notable differences between the results of the P:E ratio and both the P-and E-axes (Fig. S8).
The pollen size rate-through-time plots were generated for each tribe of Boraginaceae (Fig. 7). The results of the P-axis revealed that, among tribes of Cynoglossoideae, Cynoglosseae exhibited higher rates of change than other tribes in the subfamily. In contrast, Boragineae had a low rate of change from 30 Ma to the present; however, Lithospermeae had an increase in rate from c. 25 Ma to the present. Finally, a constant rate of change was observed in Echiochileae starting c. 40 Mya (Fig. 7). The results of the E-axis rate-through-time were almost the same as for the Paxis, with a few exceptions, such as Myosotideae having an increase after c. 10 Ma (Fig. S9). Rates for the P:E ratio increased gradually in all tribes of Cynoglossoideae after 10 Ma, except Trichodesmeae and Rochelieae, which had a fairly consistent rate. The rate of Lithospermeae increased continuously, and rates for Echiochileae and Boragineae were constant (Fig. S9).
In addition, parsimonious splits of the tree for the P-axis, E-axis and P:E ratios were identified by ParSplit. All of the parsimonious splits are A-splits (see: Chabrol and Didier, 2017). Results for both P-and Eaxes demonstrated an average of 60 changes per tree, but for P:E ratio, fewer changes were identified (54 changes). For both the P-and E-axes, the split with the lowest cost occurred in the common ancestor of the Amsinckia (Fig. S10). Based on the P-axis, the next shift was in Arnebia Forssk. and multiple other changes were in branches within Boragineae, near the crown of Nonea and Anchusa L. According to the E- axis, the second significant shift within the in-group was in the clade with Nonea, and the next shift occurred in Arnebia. Additionally, based on results of both the P-and E-axes, none of the parsimonious changes were located on the stem of Cynoglossoideae. The split with lowest cost, according to the P:E ratio, was observed in a small clade in Omphalodeae that includes Mimophytum Greenm. (Fig. S10). The next several shifts occurred in branches within Plagiobothrys and near the stem of Cynoglosseae.
In the results obtained from BAMM for both the Pand E-axes, there were multiple accelerated shifts across the family. Prominent shifts occurred in the stem of Boragineae, the stem of Amsinckia, and the crown of Arnebia. Interestingly, the parsimonious splits obtained from ParSplit reconstructed an almost identical pattern, with clades of Amsinckia and Arnebia being two splits with low costs. Additionally, multiple splits ranked three, four, six and eight occurred in the crown of Boragineae. Accelerated shifts for the P:E ratio, such as at the crown of Omphalodeae and the clade of Plagiobothrys were identified at the same nodes of ParSplit as two splits ranked one and four (Figs S8 and S10).
The results of the STRAPP analyses are presented in Table S8. Pearson, Spearman and Kruskal correlations were calculated between k, l and r for each of the P-axis, E-axis and P:E ratio. Moderate correlations between the P-axis and k (Pearson's q = 0.6542121, P = 0.015), l (Pearson's q = 0.6543182, P = 0.017) and r (Pearson's q = 0.6554462, P = 0.012) were identified based on both Pearson and Spearman correlations. Additionally, correlations between the Eaxis and k (Pearson's q = 0.6619504, P = 0.030), l (Pearson's q = 0.66144114, P = 0.027) and r (Pearson's q = 0.6606201, P = 0.030) were identified based on both Pearson and Spearman correlations. However, no significant associations for P:E ratio were identified for the macroevolutionary parameters.
The present study corroborates the findings of others (Clarke, 1977;Ahn and Lee, 1986;Bigazzi and Selvi, 1998), but with much broader sampling of nearly a third of the 1600 species of Boraginaceae.
Pollen size is one of the most variable pollen features in the family. The smallest and largest pollen grains were found in Plagiobothrys hispidulus (Greene) I.M. Johnst. (5.60 9 1.68 µm) and Anchusa arvensis (L.) M. Bieb. (52 9 35 µm), respectively. The pollen grains of these two species represent a larger pattern in Boraginaceae: the pollen of Cynoglossoideae is significantly smaller than that of Boraginoideae and Echoichloideae (Fig. 1). However, reports from the literature demonstrate a few exceptions to this pattern. Rochelia Rchb., Amsinckia (Ahn and Lee, 1986;D ıez and Vald es, 1991;Hargrove and Simpson, 2003), and Trichodesma R. Br. (Perveen et al., 1995) (a taxon that is unusual in its subfamily) have larger pollen than the rest of the genera of the Cynoglossideae and some genera of Boraginoideae (Fig. 1). Pollen size also is useful to delimit the two tribes of Boraginoideae: Boragineae (P-axis mean = 33.76 µm, E-axis mean = 24.91 µm) (Table S7) has larger pollen than Lithospermeae (Paxis mean = 19.07 µm, E-axis mean = 13.54 µm).
Although pollen size is a diagnostic feature within Cynoglossoideae, in combination with other pollen features it is even more useful. Several pollen shapes and outlines are observed across Boraginaceae, confirming and extending previous observations (Clarke, 1977;Perveen et al., 1995). The members of the family can be distinguished by pollen shape and shape of the equatorial outline. The P:E ratio of Cynoglossoideae, except Trichodesmeae, is larger than in other subfamilies. The most elongated pollen grains are observed in Omphalodeae and Rochelieae, with pollen in these tribes having a P:E ratio greater than that of other tribes and subtribes (except a few species of subtribe Amsinckiinae) of Boraginaceae (Fig. 1 and Table S7), and the most spherical pollen grains are observed in Trichodesmeae and Echiochileae (Fig. 2 and Table S7). The most common pollen shape outlines in Cynoglossoideae are dumbbell-shaped and elliptical with a slightly restricted equator, shapes that exist in all tribes of the subfamily except Trichodesmeae, and this shape is present in >180 of the 314 examined species of Cynoglossoideae (Fig. 2 and Table S2). Indeed, these two types of outlines vary more or less in the extent of equatorial constriction, and this variation may have been affected by pollen hydration (Hargrove and Simpson, 2003). Alternatively, the variation in equatorial constriction could also be a real biological phenomenon and further studies could help confirm the extent of constriction in the family. Pollen with these outlines was not found in either Boragineae or Echiochileae, and only a few taxa in Lithospermeae bear pollen with these outlines, which originated independently (e.g. Lithospermum, Ancistrocarya) (Table S2). Therefore, it can be concluded that pollen with a dumbbell-shaped outline (or/ and elliptical with slightly restricted equator) is diagnostic of Cynoglossoideae, even though other types of pollen are present in the subfamily, and particularly in combination with other palynological features. Within Cynoglossoideae, Asperugeae + Omphalodeae, which are sister groups (Fig. 6), share the same pollen shape and outline (dumbbell and elliptical with a slightly constricted equator outlines, and prolate and perprolate shapes), and Cynoglosseae + Myosotideae (Fig. 6) share similar pollen shape and outlines (elliptical and elliptical-rectangular outlines and spheroidal shape) that are absent in the rest of Cynoglossoideae (Fig. 2). Trichodesmeae bears pollen with elliptical and circular outlines and spheroidal and subprolate shapes, and this type of pollen differs from other members of the subfamily, which is consistent with it being sister to the rest of Cynoglossoideae. However, these results are based on nine investigated (of c. 50) species of the group, so further studies can provide a greater understanding of these patterns.
Aperture type and number are also diverse in the family. Boraginoideae has isoaperturate pollen, and the number of apertures ranges from three to than ten (Fig. 4 and Table S2). All Cynoglossoideae, except Trichodesmeae, bear heteroaperturate pollen. The members of the subfamily usually have three pseudoapertures and three apertures; however, reports from the literature provided evidence that a few species of Myosotideae have four to six apertures and four to six pseudoapertures (D ıez and Vald es, 1991;Hargrove and Simpson, 2003) (Fig. 4 and Table S2). Additionally, some members of Amsinckiinae, such as some species of Cryptantha, Plagiobothrys and Pectocarya, also have polar pseudoapertures (Spaeth, 2014). The presence of heteroaperturate pollen grains is a shared trait that exists in Cynoglossoideae and not in Boraginoideae or in Echiochiloideae, which are both isoaperturate (Fig. 4) confirming previous observations (D ıez and Vald es, 1991;Perveen et al., 1995;Bigazzi and Selvi, 1998).
Given the most recent taxonomy of Boraginaceae, the present study provides evidence that pollen characteristics are useful for diagnosing various subfamilies and tribes, such as Cynoglossoideae and Boraginoideae (Table S2). The pollen grains of tribes of Cynoglossoideae are small, primarily perprolate in shape, and dumbbell-shaped in outline, and include six apertures (three apertures and three pseudoapertures), making the pollen easy to distinguish from that of Boraginoideae and Echiochiloideae. Boragineae and Lithospermeae, the two tribes of Boraginoideae, are similar in aperture type and number, but they can be distinguished based on pollen size and outlines. As mentioned earlier, pollen size in Boragineae is significantly larger than in Lithospermeae (Fig. 1), so pollen size could be considered, with the other features, diagnostic of the tribe. In addition, ovoid, dumbbellshaped, spheroidal and elliptical with a slightly constricted equator outlines are four pollen outlines that do not exist in Boragineae, with the other three shapes (elliptical, elliptical-rectangular and rectangular) present in members of Boragineae. Interestingly, the pollen grains of Trichodesmeae, which is sister to the rest of the tribes of Cynoglossoideae, is more similar to those of Boraginoideae than members of Cynoglossoideae. Thus, the inclusion of Trichodesmeae in Cynoglossoideae is not supported by the pollen features, and the pollen may actually provide evidence for the recognition of another subfamily--Trichodesmoideae--should the pollen differences with Cynoglossoideae be considered sufficiently distinct (along with other possible features) to elevate the tribe to subfamilial status.
Although pollen characters lend support to the higher-level systematics in Boraginaceae, these features also could be helpful for lower taxonomic levels, such as subtribes and genera. For example: Ancistrocarya appears to be most closely related to Lithospermum (Johnston, 1954), and Chac on et al. ( 2019) provide molecular phylogenetic evidence that Ancistrocarya is a member of a clade that includes Lithospermum, Buglossoides, Aegonychon Gray and Glandora. The present study is the first to incorporate an analysis of pollen morphology of Ancistrocarya, and the results support the phylogenetic placement of the genus in the position resolved by Chac on et al. ( 2019). Additionally, we observed that the pollen size (c. 11 9 8 µm), perprolate shape and dumbbell outline of Ancistrocarya japonica Maxim. are very similar to these features in the pollen of Lithospermum officinale L. Given that pollen features are diagnostic at various taxonomic levels, it could be quite useful to identify pollen fossils of Boraginaceae and reconstruct changes in past vegetation and determine flexibility of taxa during past climate and environmental changes.
Pollen size. The detailed reconstruction of the evolution of pollen size obtained from the TNT, Phytools and Mesquite analyses are shown in Table 1. The values of the P-and E-axes reconstructed from the ML method in Phytools were more similar to those of squared-change parsimony than to those from TNT; however, minimal differences were observed among all three methods.
Several hypotheses about the relationship between pollen grain size and genome size (Beaulieu et al., 2008;Knight et al., 2010;Srisuwan et al., 2019), pollinators (Harder, 1998), floral morphology, origins of heterostyly (Cohen, 2014), pollination mechanisms and climatic conditions (Ejsmond et al., 2011) have been proposed to describe evolutionary changes in pollen size. Kobrlov a and Hrone s (2019) studied variation of genome size in Boraginaceae by examining all native taxa of the family in the Czech Republic. Their efforts revealed that members of Boraginoideae had significantly larger genome sizes than members of Cynoglossoideae, and our observations demonstrate that pollen grains are significantly larger in Boraginoideae than Cynoglossoideae. These results are, in general, consistent with data on Boraginaceae from the Plant DNA C-values Database (Leitch et al., 2019). For instance, Anchusa arvensis, a species with one of the largest pollen grains (52 9 35 µm) in the family, has a genome more than four-fold larger (2C = 9.10 pg; see Leitch et al., 2019) than those in other members of the family, such as Cynoglossum officinale L. (2C = 1.82 pg; see: Leitch et al., 2019), which has one of the smallest pollen grains (13 9 9 µm). There are a few exceptions, for example Amsinckia furcata Suksd. has larger pollen grains (34 9 24 µm) than Halacsya sendtneri D€ orfl. (21 9 20 µm), and yet both have the same genome size (2C = 2.60 pg; see: Leitch et al., 2019). In addition, Onosma stellulata Waldst. & Kit. has a larger genome (2C = 5.56 pg; see: Leitch et al., 2019) than Amsinckia furcata, despite having smaller pollen grains (23 9 20 µm). Regardless of these exceptions, comparisons of these studies with the findings of the present study (Fig. S11) provides evidence for a correlation between pollen size and genome size. Although the results may be a weak trend (Fig. S11), the different evolutionary trajectories of pollen size in Boraginoideae and Cynoglossoideae would be an interesting area to explore in the future. Recently, the evolution of genome size and its relationship with pollen grain size was investigated in Ipomoea batatas (L.) Lam. (Srisuwan et al., 2019). The authors suggested that genome size and ploidy levels were positively correlated with pollen size. Juli a and Villodre (1994) also showed that pollen size increases with ploidy level in Puccinellia Parl. (Poaceae). Nevertheless, there are cases, such as Tarasa Phil. (Malvaceae), in which pollen size and ploidy level are not correlated (Tate and Simpson, 2004). In general, larger genome sizes are not usually associated with larger plants. However, just because the plants themselves are not larger does not mean that their cells are not larger, which is often the case with polyploid species (Tsukaya, 2013). Therefore, genome size may not be associated with larger plants overall but could be related to larger pollen as pollen includes only two to three cells.
Pollen shape. Patterns of evolution of pollen shape and outline-in-equatorial-view revealed that the ancestral condition for the family was prolate pollen with an elliptical-rectangular outline. These two shapes evolved independently multiple times across Boraginoideae. However, there are trends toward perprolate and dumbbell-shaped pollen within Cynoglossoideae (except Trichodesmeae), and similar patterns are resolved in Boraginoideae. In addition, results indicated that pollen shapes and outlines are more diverse in Boraginoideae, especially in Lithospermeae, than in Cynoglossoideae.
As mentioned earlier, pollen outlines are more diverse in Lithospermeae than in other tribes in the family. It seems that the pattern of evolution of pollen outlines of Lithospermeae proceeded from ellipticalrectangular to ovoid and dumbbell-shaped. However, as a consequence of some reversals to the ancestral condition and some independent origins, it is not possible to specify a particular pattern for the evolution of pollen outline throughout the entire tribe. The ovoid outline occurs only in Lithospermeae and originated twice in separate clades--once in Alkanna and once in a large clade that includes Echium L. + Lobostemon + Pontechium B€ ohle & Hilger + Onosma L.--and both of these clades originated in the early to middle Miocene. During the Miocene, the crown node of Lithospermum, which have elliptical with slightly constricted equator and dumbbell-shaped outlines, diverged. Throughout the Miocene, an increase in global temperature and reduced rainfall led to expanding open vegetation systems and decreasing closed vegetation (Polly et al., 2011). Two potential patterns may have occurred during the evolution of these outline shapes. On the one hand, the ovoid shape arose before the pear shape (Fig. 8), then elliptical pollen with a slightly restricted at equator (this state looks like pear shaped but has ends of equal size) and finally, the dumbbell-shaped outline (this state looks like elliptical with slightly restricted equator but is extremely restricted at the equator) (Fig. 8). On the other hand, there is no linear pattern between these four states of outline shapes, and each likely evolved independently as a result of a particular advantage, such as adaptation to a new environment. For instance pollen shape in Amsinckia, which grows in dry conditions (Spaeth, 2014), tends to be mostly prolate and has an elliptical outline, whereas pollen of Plagiobothrys, which grows in wet conditions (Spaeth, 2014), tends to have mostly perprolate shape and outline-in-equatorial-view is diverse between rectangular to dumbbell-shaped. Recently, Yang et al. (2020) recognized a significant correlation between pollen shape and habitat moisture in lamiid members, and the researchers confirmed that habitat moisture could be a selective force on pollen morphology in angiosperms. However, clarifying the effect of new conditions on pollen shape in Boraginaceae will require further study.
The advantages of pollen shape diversity have been investigated by previous authors (Ghorbel and Nabli, 1998;Bigazzi and Selvi, 2000;Cohen, 2010Cohen, , 2014)). These authors hypothesized that pollen shape and stigma papillae were associated functionally, such as restricting entrance and survival of alien pollen grains from adhering and germinating on the stigma. Consequently, the current data can be used for further research, in order to clarify this association of stigmatic papillae or climatic condition and pollen shape.
Aperture number. The analysis of pollen aperture number revealed that the presence of three apertures is ancestral in the family. The ancestor of Cynoglossoideae is three-aperturate, which is constant across the subfamily. However, three apertures originated twice in Lithospermeae: once in the clade of species of Alkanna and the other one in a large clade that includes Cystostemon Balf.f. + Maharanga DC. + Onosma and Echium + Ponthechium + Lobostemon + Echiostachys.
Interestingly, this character state is associated with ovoid pollen in these two groups of the tribe. There are multiple origins of various aperture numbers, and patterns differ between the tribes. The direction of aperture number evolution suggests that aperture number is quite stable in Cynoglossoideae, but increased in Boraginoideae, with some reversals to three apertures in particular clades. This increase in aperture number is associated with greater diversity in pollen shape and outline as well. Many aperture patterns exist in flowering plants, and pollen with three apertures is dominant across eudicots (Albert et al., 2018). An increase in pollen aperture number in Boraginaceae evolved multiple times. This trend is observed in other families as well, such as Plantaginaceae (Al-Quran, 2004), Rosaceae (Chung et al., 2010) and Lamiaceae (Doaigey et al., 2018). This provides a possible benefit by increasing the number of potential germination sites, which increases the probability of connection between an aperture and the stigmatic surface resulting in an increase in the fertilization rate (Furness and Rudall, 2004). Dajoz et al. (1991) found that an increase in aperture number is associated with faster germination and lower life expectancy, but Albert et al. (2018) demonstrated that three aperturate pollen had a significantly higher rate of germination than pollen with more apertures. Furthermore, aperture number has an influence on pollen viability under osmotic stress, and pollen with fewer apertures survives better than pollen with more apertures (Prieu et al., 2016). Thus, the climate of species of Boraginoideae might be more stable in areas inhabited by species of Cynoglossoideae, like Mediterranean-type ecosystems vs. alpine regions, and this can explain the advantage of three apertures resulting in the prevalence of the trait throughout Boraginaceae, and especially in Cynoglossoideae and other species of eudicots (and see Aperture type section below). Although speculative, pollen aperture number could be one possible reason that Cynoglossoideae is quite species-rich and has a larger geographical range compared to the other subfamilies.
Like other pollen grain traits, the type of aperture varies across Boraginaceae. Both heteroaperturate and isoaperturate pollen are found in the family. Isoaperturate pollen is observed at the origin of Boraginaceae (Fig. 6), and heteroaperturate pollen originated c. 32 Ma. In heteroaperturate pollen, two different types of apertures are present with only one type assumed to function as a potential germination site (Halbritter et al., 2018). Heteroaperturatate pollen is the main type in Cynoglossoideae, excepting Trichodesmeae. The effect of pseudoapertures in harmomegathy, a characteristic folding of pollen grains to accommodate the decrease in cellular volume resulting from water loss (Volkova et al., 2013), has been examined only rarely in previous studies; however, there is evidence that pseudoapertures have a role in harmomegathy (Payne, 1972;Blackmore and Barnes, 1986;Scotland et al., 1990). Volkova et al. (2013) studied different structures involved in the harmomegathic mechanisms in five species in Boraginaceae. Their studies demonstrated that although apertures and pseudoapertures seem to function as the primary elements of harmomegathy, the remainder of the pollen-wall-like surface features have harmomegathic functions as well. For example, the pollen of Cryptantha coryi I.M. Johnst. and C. celosioides Payson includes pseudoapertures, a transverse groove and a perforated tectum. The pseudoapertures and transverse groove decrease the folding ability of the pollen, so in order to provide greater flexibility, the pollen grains possess a perforated tectum on the poles (Volkova et al., 2013). The pollen grains of Myosotis palustris (L.) Lam. have pseudoapertures, a perforated tectum and two well-defined triangular polar poroid areas that are involved in harmomegathy. These polar poroid areas are folded inward in a dehydrated state and structurally resemble meridional pseudoapertures (Volkova et al., 2013). This evidence provides a clarification for the tendency of pseudoapertures across most of the members of Cynoglosoideae. Other than Boraginaceae, pseudoapertures are known in several families of angiosperms, such as Acanthaceae (Wood and Scotland, 2009), Verbenaceae (El Ghazali and Krzywinski, 1989) and Fabaceae (Caccavari, 2002). The role of pseudoapertures is not clearly understood in these families, but it seems that it has harmomegathic functions as well. Consequently, in heteroaperturate pollen, harmomegathy and germination functions are performed by different parts of the pollen. Apertures associated with function of germination and pseudoaperture also provide additional flexibility (Volkova et al., 2013). Thus, this feature of heteroaperturate pollen (transference of function to different structures; Volkova et al., 2013) could be considered as another advantage for Cynoglossoideae over other subfamilies.
Examination of evolutionary patterns of pollen size, shape and outline, and type and number of apertures demonstrates that after the divergence of Trichodesmeae, the sister tribe of the rest of Cynoglossoideae, the pollen size of Cynoglossoideae became significantly smaller than in other subfamilies. The number of apertures was reduced to three, and a new type of aperture, the pseudoaperture, evolved throughout the subfamily. In addition, pollen shape and outline in Cynoglossoideae are less diverse than in other subfamilies and are limited to the prolate and perprolate shapes, and elliptical with slightly constricted equator and dumbbell-shaped outlines. The four other tribes in the family (excepting Trichodesmeae)--Echiochileae, Boragineae, Lithospermeae, and Trichodesmeae--started to diversify during the late Eocene to early Oligocene. During this time, lower temperatures and increased seasonality led to a shift towards a reduction in forest and increasingly savanna-like vegetation (Polly et al., 2011). This was the period when pollen grain variation was increasing in Boraginaceae leading to high diversity in pollen shapes and outlines, larger pollen and more apertures across the family, especially in these four tribes. Except Trichodesmeae, all tribes of Cynoglossoideae diversified during the Oligocene to early Miocene. The pattern of pollen shape, size and aperture in the subfamily could be correlated with a cooling trend throughout the Oligocene to part of the Miocene epoch. Thus, all of these trends, especially the origin of pseudoapertures and decreasing pollen size, in Cynoglossoideae might have been possible adaptive responses to these new environmental and ecological conditions. Additionally, this pollen may have provided suitable conditions for colonization following long-distance dispersal and migration to new habitats, which may have helped facilitate the large geographical distribution for this group. On the one hand, the high number of widespread taxa in the subfamily (e.g. Lappula, Myosotis and Mertensia Roth) supports this hypothesis. On the other, high diversity in pollen shape in Boraginoideae could be associated with diversification of flower morphology and pollination systems, which are more diverse in Boraginoideae than other subfamilies.
Macroevolutionary analyses provide evidence for a shift in Boraginaceae where diversification rates slow toward the present. Additionally, the biogeographical history of the crown of Boraginaceae dates to the Eocene (or perhaps earlier--Palaeocene c. 55-72 Ma; see Chac on et al., 2017). Global warming during the Palaeocene-Eocene (Polly et al., 2011) may have had a remarkable effect on patterns of diversification. Phylorate plots (Figs S5 and S6) demonstrate that the rate of speciation for most of the family is faster than that for the early diverging subfamily Echiochiloideae. This taxon is sister to other members of the family (which originated c. 45 (41.5-50.1) Ma; see Chac on et al., 2017), and it includes the smallest number of species of the three subfamilies of Boraginaceae. Additionally, Chac on et al. ( 2017) suggested that the extinction of lineages occurred near the backbone of Boraginaceae during the Palaeogene that could be related to changes in the climate in regions inhabited by Boraginaceae, such as increased aridification in some parts of Eurasia. Thus, these explanations provide reasons that Echiochiloideae had a slower rate of diversification than the rest of the family, which is much more diverse. The results obtained in the present study demonstrate that the highest rates of diversification occurred within Plagiobothrys and Cynoglossinae. The plots of speciation through time (Fig. S7) revealed an increase in the rates for the clade of Cynoglosseae at c. 10 Ma, during the middle Miocene, when several events occurred, such as increasing global temperatures, decreasing precipitation, expanding open vegetation systems and decreasing closed vegetation (Polly et al., 2011). Aridification increased in Eurasia during this epoch, and steppe vegetation started to appear. In southern Asia, grassland ranges extended, producing greater diversity in habitats. Additionally, aridity increased in southern Europe because the Tethys Sea connection between the Mediterranean and Indian Ocean was separated in the mid-Miocene (Polly et al., 2011). These palaeoclimatic and tectonic events could have affected this pattern of diversification and geographical distribution. The biogeographical events that were inferred by Chac on et al. ( 2017) for the Cynoglossoideae involved many Oligocene and Miocene dispersals, including from Western Europe to the Irano-Turanian region (e.g. in the Cynoglossineae), long-distance dispersals from eastern Asia to North America, and more dispersal events from North America to South America (Amsinckiinae) (Guilliams et al., 2017). The deformation of the Iranian and Anatolian plateaus and the Afghan block during the late Miocene and Pliocene may have caused the formation of a land corridor between the western and eastern areas of the western Irano-Turanian region, which may have prompted diversification rate of Cynoglossinae through increasing aridity and topographic heterogeneity in this region (Sherafati et al., 2021). Thus, the increased rate of diversification for species-rich groups, like Cynoglossum and relatives and Plagiobothrys and relatives, seems consistent with radiations in the Mediterranean region (s.l.) and western North America, respectively.
Evolutionary rates and patterns of shifts in pollen size were estimated using two different methods, ParSplit and BAMM. The results of the novel parsimonybased method, ParSplit, were quite similar to those from BAMM. Indeed, splits with low costs occurred at the same (or very near) nodes as those with accelerated shifts in BAMM. However, some of the high-cost splits (such as splits in the Mertensia clade that ranked 50 th and 55 th ) were not detected by BAMM. The most parsimonious split occurred in Amsinckia, which has the largest pollen grains among members of Cynoglossoideae. This lowest cost split may indicate a significant selective pressure of the environment on the evolution of pollen size in this genus. In both ParSplit and BAMM analyses, multiple shifts were observed within or toward the base of each of the tribes. Therefore, these clade-specific patterns may correlate with specific biogeographical and fairly recent ecological conditions. BAMM results revealed that most of the clades were nearly stable or had low rates of change. However, there were multiple accelerated shifts across the family. The most notable of them occurred in the stem of Boragineae, which also is confirmed by ParSplit analyses. The accelerated rate in Boragineae and a trend to greater pollen size may be linked to species adaptations to novel habitats during the late Eocene. Alternatively, according to the punctuated equilibrium theory (Eldredge and Gould, 1972), these shifts and evolutionary trends may be due to the connection between speciation and morphological evolution. Although morphological changes appear also to shift with an accelerated rate of speciation, we did not find an association between an accelerated rate of evolution for pollen size in Boragineae and speciation. Also, the increased rate of diversification for species rich groups, like Cynoglossum and relatives and Plagiobothrys and relatives, is not correlated with pollen size evolution. Thus, there is likely a nonpollen-related explanation for diversification of these groups. There could also be a feature related to pollen that was not tested in this study, and perhaps the presence of pyrrolizidine alkaloids (El-Shazly and Wink, 2014;Ahmad et al., 2018) could provide another explanation for the success of particular groups. However, results of the STRAPP analyses revealed that the P-and E-axes had a significant correlation with macroevolutionary parameters based on Pearson and Spearman methods, resulting in moderate evidence demonstrating and association between pollen size and patterns of macroevolution. Moreover, it is notable that the P:E ratio is not correlated with patterns of macroevolution, which could suggest that it is pollen size, not shape (i.e. P:E ratio), that results in the shifts in diversification. Lithospermum and relatives are resolved as one of the few clades with an identified shift in P:E ratio, but not P-or E-axes, a case that differs from the usual pattern. Chac on et al. ( 2019) indicated four shifts in diversification rates of tribe Lithospermeae, one in each Alkanna, Onosma, Echium and Lithospermum. Different sampling may be the reason that the results from the present study differ from theirs. Nevertheless, comparisons of shifts in highly diverse genera of Lithospermeae with those for pollen size demonstrated that there might be a correlation between pollen size and shifts in diversification rates of Echium and relatives.
According to our results, pollen size was likely a factor in diversification, and it is likely one of multiple factors. Perhaps other features, such as stigma papillae, pollinators and pyrrolizidine alkaloids, could also have played a role, with pollen being just one of the features involved in the resulting diversification. Collectively, these results suggest that pollen size alone is not correlated with diversification of species of Boraginaceae. Consistent with previous analyses (Vamosi et al., 2018), it seems more likely that a combination of features, rather than one trait alone, may be associated with innovation in the family, and the diversification rate shifts represented here cannot be assigned to single factor, such as pollen size, or to changes in one organ or stage in the life cycle.
This is the first comprehensive investigation on pollen diversification and evolution across Boraginaceae. Distinctive and unique combinations of pollen features characterize most of the tribes in the family. Across tribes, pollen size, shape and aperture number are most variable in Boragineae and Lithospermeae. Indeed, variation in floral morphology, such as zygomorphy, heterostyly (Cohen, 2014) and specific pollination systems (Bigazzi and Selvi, 1998), is accompanied by greater pollen grain diversity in Boraginoideae than in Cynoglossoideae and Echiochiloideae. Generally, pollen characteristics differ among the three subfamilies (i.e. Echiochiloideae, Boraginoideae and Cynoglossoideae), and there are pollen size, shape and aperture type differences between Trichodesmeae and the other tribes in Cynoglossoideae. The ancestral character reconstructions of pollen morphological traits identify diagnostic characters to many of the clades and lend additional support to the recent classification of the family (Chac on et al., 2016). Further examination is suggested for pollen exine ornamentation as it is diverse in Boraginaceae and can play an important role in the evolution, ecology, biogeography, diversification and classification of the family.