Reconstructing Bioinvasion Dynamics Through Micropaleontologic Analysis Highlights the Role of Temperature Change as a Driver of Alien Foraminifera Invasion

Invasive alien species threaten biodiversity and ecosystem structure and functioning, but incomplete assessments of their origins and temporal trends impair our ability to understand the relative importance of different factors driving invasion success. Continuous time-series are needed to assess invasion dynamics, but such data are usually difficult to obtain, especially in the case of small-sized taxa that may remain undetected for several decades. In this study, we show how micropaleontologic analysis of sedimentary cores coupled with radiometric dating can be used to date the first arrival and to reconstruct temporal trends of foraminiferal species, focusing on the alien Amphistegina lobifera and its cryptogenic congener A. lessonii in the Maltese Islands. Our results show that the two species had reached the Central Mediterranean Sea several decades earlier than reported in the literature, with considerable implications for all previous hypotheses of their spreading patterns and rates. By relating the population dynamics of the two foraminifera with trends in sea surface temperature, we document a strong relationship between sea warming and population outbreaks of both species. We conclude that the micropaleontologic approach is a reliable procedure for reconstructing the bioinvasion dynamics of taxa having mineralized remains, and can be added to the toolkit for studying invasions.


INTRODUCTION
Global trade and worldwide transport of people and goods have largely altered the natural distribution of species (e.g., Sardain et al., 2019), but large gaps of knowledge remain in assessing spatial and temporal patterns of invasions, especially in the marine environment (Seebens et al., 2017;Ojaveer et al., 2018). Incomplete assessments of temporal trends, origins and drivers of their spread may affect our ability to understand the ecology and history of communities and the relative importance of different factors driving invasion success (Carlton, 2009). Sporadic monitoring and surveillance, weak taxonomic knowledge, or elusive behavior of some species have left several marine invasions undetected for years, decades or even centuries (Carlton, 2009;Griffiths et al., 2010;Zenetos et al., 2019). Several approaches have been attempted to reconstruct the timing of first introduction events and subsequent stages of invasion: reexamination of old museum or herbarium collections (Ahnelt, 2016;Steen et al., 2017), analysis of published descriptions (Zullo, 1992;Galil et al., 2018), interviews to local fishermen (Bariche et al., 2014;Azzurro et al., 2019), molecular tools (Ordóñez et al., 2016;Deldicq et al., 2019) and radiometric dating (Petersen et al., 1992;Albano et al., 2018). In spite of different approaches utilized, the gap still persists because long and continuous timeseries are difficult to obtain, due to the often sporadic and intermittent availability of marine biodiversity data; hence several studies are restricted to inter-decadal comparisons (Parravicini et al., 2015) or are based on population model simulations (Clark et al., 2013;Walsh et al., 2016). Moreover, most of the current knowledge of invasions is from larger organisms, whereas very little is known about spatial and temporal patterns of bioinvasions by microscopic taxa such as foraminifera or other unicellular eukaryotes (Langer et al., 2012;Skarlato et al., 2018;Reavie and Cangelosi, 2020). These small-sized invaders can remain unnoticed for several years after their introduction and gain attention from beach visitors or scientists when population sizes reach a certain threshold and become visible to the naked eye (Guastella et al., 2019).
An example is provided by two algal symbiont-bearing benthic foraminifera: Amphistegina lobifera Larsen, 1976, an Indo-Pacific species that entered the Mediterranean via the Suez Canal (Prazeres et al., 2020), and its congener A. lessonii d'Orbigny, 1826, a species of uncertain native origin. Discussions on the non-indigenous status of A. lobifera and cryptogenic status (sensu Carlton, 2009) of A. lessonii are provided in Guastella et al. (2019). In particular, fossil specimens of A. lobifera have never been found in ancient Mediterranean records, while fossil remains of A. lessonii are common in Mediterranean shallow-water strata, late Pliocene to early Pleistocene in age (Di Bella et al., 2005). Both species are currently abundant in the Red Sea, but only A. lobifera is successfully spreading in the Mediterranean Sea, following a westward colonization process (Langer et al., 2012;Caruso and Cosentino, 2014;Langer and Mouanga, 2016;Guastella et al., 2019). On the contrary, A. lessonii is extremely rare in the Mediterranean Sea, and is mainly characterized by smaller-sized specimens compared to A. lobifera (Guastella et al., 2019). This sporadic occurrence and size range has been attributed to its different sensitivity to colder temperature that may limit test calcification (Titelboim et al., 2019).
Both species commonly proliferate in shallow-waters (<20 m depth) of warm tropical and subtropical areas, and their latitudinal distribution is strongly regulated by sea surface temperatures, with winter minimum temperatures around 18 • C (Hallock, 1984;Langer and Hottinger, 2000;Hohenegger, 2004;Renema, 2018). In the Mediterranean Sea, their geographic distribution and population abundance is believed to be primarily controlled by wintry sea surface temperature (winter isotherm 15 • C, Hollaus and Hottinger, 1997). However, little is known about Amphistegina spp. populations established in the Central Mediterranean region, which may have been selected for their higher tolerance to colder temperatures, compared to Red Sea and Eastern Mediterranean populations (Schmidt et al., 2016a;Titelboim et al., 2019). Nonetheless, a further spread of amphisteginids into the Western Mediterranean region could be expected in the near future, also promoted by the ongoing sea surface warming (Langer et al., 2013;Guastella et al., 2019). Profound ecological changes have been documented in the Eastern Mediterranean, where A. lobifera populations have become the dominant component of benthic foraminiferal assemblages. Massive deposits of amphisteginids produced an increased accumulation of carbonatic sands, which progressively altered native substrates. This, in turn, reduced biodiversity and caused the displacement of previously established foraminiferal biota (Langer et al., 2012;Mouanga and Langer, 2014).
In this study, we aim to reconstruct the invasion history of A. lobifera and A. lessonii using a multidisciplinary approach based on the analysis of sediment cores by techniques of micropaleontology, paleoecology, and radiometric dating. We document their population abundance along the timeline, by analyzing two sedimentary cores (CORE18 and CORE19) collected from Malta island (Central Mediterranean Sea). Additionally, we relate the temporal dynamics of the two populations with the trends of sea surface temperature (SST) measured in the area in the last 60 years, in order to assess whether sea warming has driven their invasion pattern. Finally, we elaborate hypotheses on the possible future spreading of both species in the Mediterranean Sea.

Study Area, Selection of Sampling Site, and Sample Collection
The Maltese Archipelago is composed of three main islands (Malta, Gozo and Comino) and is located within the Sicily channel in the Central Mediterranean Sea, at a latitude of 35 • 48 28 -36 • 05 00 North and a longitude of 14 • 11 04 -14 • 34 37 East ( Figure 1A). The coasts of the archipelago, with a total perimeter of about 271 km, are mainly characterized by high rocky shores; only a few natural bays are present and most of them exhibit high levels of anthropogenic impact (Magri, 2006).
In order to identify suitable sites for collecting the sediment cores, where the sea-floor could reasonably reflect the pattern of sedimentation that occurred in the past decades and thus the history of colonization of Amphistegina spp. in the area, we complied with the following criteria: (I) occurrence of a well-established population of A. lobifera and subordinately of A. lessonii ( Figure 1B); (II) absence of fossil amphisteginids in neighboring outcropping rocks, in order to avoid their reworked occurrence in the studied cores, that may confuse the record of current foraminiferal assemblages of marine sediments. To satisfy requirement II, relative abundance of reworked foraminiferal FIGURE 1 | (A) Study area (Central Mediterranean Sea); (B) Current distribution and relative abundance of A. lobifera (in red) and A. lessonii (in green) along the coasts of the Maltese Islands. The pie charts in this panel derived from data available in Guastella et al. (2019) and are supplemented with new data regarding relative abundance of reworked fossil species (in light-blue); (C) Sampling site in Marsamxett Harbor. In the selected site, the relative abundance of A. lobifera and of reworked fossil species are about 40 and 3%, respectively. Additionally, the collection site is close to one site previously studied (Yokeş et al., 2007) (black star). Maps derived from Google Earth (https://www.google.it/intl/it/earth/; data providers: Google Landsat/Copernicus Data SIO, NOAA, U.S. Navy, NGA, GEBCO IBCAO) and modified with GIMP v. 2.10.20 (https://www.gimp.org/). species present in sediment samples was calculated through the comparison with fossil assemblages preserved in seven rocky samples collected in Malta from Miocene outcrops: three rocky samples were collected from the Globigerina Limestone Formation and prepared as washed residues, and four rocky samples were collected from the Upper Coralline Limestone Formation and prepared as thin rock sections ( Figure 1B and Supplementary Table 1 available on-line as Supplementary Data); (III) proximity to one of the four sites of the earliest record of A. lobifera in Malta (Yokeş et al., 2007), in order to have a temporal constraint for the reconstruction of the invasion dynamics; (IV) prevalence of fine-grained sediments, first requirement necessary for 210 Pb dating, a widely applied method to date recent sediments and to assess sediment accumulation rates (Tylmann et al., 2016;Andersen, 2017); (V) location within an enclosed bay, to ensure protection from storm waves and littoral currents that would alter the normal sedimentation process. In fact, the presence of an undisturbed sediment deposition is the second requirement to correctly apply 210 Pb-chronology (Tylmann et al., 2016); (VI) absence of dredging or beach nourishment operations in the area, for the same reason as above; (VII) depth >10 m to avoid sediment mixing by waves, as above, and <20 m, the maximum depth where A. lobifera has been observed to occur (Hallock, 1984;Hohenegger, 1995;Guastella et al., 2019).
Despite many efforts and surveys to search for multiple suitable sites along the coasts of the Maltese Archipelago, only a single locality satisfied all the stringent conditions listed above: Marsamxett Harbor (35 • 54 16.7 N; 14 • 30 27.5 E; Figure 1C). Located near station "4" of Yokeş et al. (2007) and known to have in 2018 a relative abundance of A. lobifera of ∼ 40% (Guastella et al., 2019), this natural bay displays the required condition of sediment grain size, shelter, depth and absence of human activities that could have altered sedimentation on the sea-floor. Furthermore, Amphistegina spp. populations are well established (Guastella et al., 2019) and the rocky outcrops around the selected site are exclusively characterized by rocks belonging to the Globigerina Limestone Formation (Aquitanian-early Langhian in age), formed in outer shelf to upper slope environments (Baldassini and Di Stefano, 2017). Thus, no fossil amphisteginids were recorded and the relative abundance of the other reworked fossil foraminifera (mainly deep-water taxa such as Cibicides pseudoungerianus, Heterolepa bellincionii, Neoeponides schreibersii, Reussella spinulosa, and Spiroplectammina carinata and the planktonic genera Globigerina, Globigerinoides, and Globoquadrina) is only 3% (Figure 1B  Within Marsamxett bay, two sediment cores were collected a few meters away from each other using a hand-corer operated by a scuba diver. The first core (CORE18) was sampled on 8th May 2018 at 16 m depth and the second one (CORE19) was collected on 4th September 2019 at 17 m depth. CORE18 was collected as a "pilot" sample in order to evaluate if all the stringent conditions listed above were respected and if the collected data had a high quality; CORE19 was sampled in order to have a replicate that could confirm the reconstruction of invasion dynamic provided by CORE18.
CORE18 is 41 cm long and CORE19 is 50 cm long, both of them mostly made up by fine-grained sediments (fine sands, silts and clays). After collection, both cores were longitudinally sectioned in two halves and then crosscut at each centimeter, obtaining 41 samples from the first one and 50 samples from the second one, respectively. Sediment samples were used for both grain-size and radiometric analyses and prepared for the micropaleontologic analyses.

Grain-Size and Radiometric Analyses
Grain-size analysis was carried out for CORE18 in order to verify its compliance with two prerequisites for radiometric dating: finegrained sediments and undisturbed record (not mixed vertically). Sediment samples were oven-dried at 40 • C for 1 day, weighed (up to a maximum of 17 g) and separated by wet sieving using five overlapped sieves (meshes, respectively of 1 mm, 500, 250, 125, and 63 µm) as indicated in Blott and Pye (2012). The fraction retained on each sieve was oven-dried and weighed again, obtaining weight percentage data in the form of cumulative curves (distinguishing between very coarse, coarse, medium, fine and very fine sands and indistinct mud fraction); the weight of the mud fraction (silts plus clays) was calculated by difference of the total weight.
Both cores were chronologically constrained by measuring activities of 210 Pb and 137 Cs isotopes along the sedimentary record. 210 Pb is a natural radionuclide belonging to the 238 U decay series and is characterized by a half-life of 22.3 years. The total 210 Pb activity in marine sediments has two components: supported 210 Pb activity derived from the decay of in situ 226 Ra, and unsupported 210 Pb activity derived from atmospheric fallout (Kosnik et al., 2015;Andersen, 2017). In marine environments, unsupported 210 Pb deposits associated with muddy particles accumulate at the sediment-water interface as excess 210 Pb. The dating is based upon the determination of the vertical distribution of unsupported 210 Pb, by subtracting supported 210 Pb activity from the total activity of 210 Pb (Hollins et al., 2011).
In the studied core samples, 210 Pb activity was measured at the CNR-ISMAR Institute (Bologna, Italy) via alpha counting of its daughter isotope 210 Po, assuming secular equilibrium between the two isotopes as described in Rizzo et al. (2009) and Incarbona et al. (2016). 210 Po was extracted from the sediment using hot HNO 3 and H 2 O 2 , and spiked with 209 Po (NIST standard SRM 4326, diluted to 0.43 Bq g −1 ) used as a yield monitor. After separation of the leachate from the residue, the solution was evaporated to near dryness and the nitric acid was eliminated using concentrated HCl, the residue was then dissolved in 1.5 N HCl and Iron was reduced using ascorbic acid. Finally, Po isotopes were plated onto a silver disk overnight, at room temperature. Repeated measurements of a certain number of sediment samples was carried out to estimate the analytical precision. Sedimentation rates for the last decades were calculated based on the decreasing concentration of excess 210 Pb, following the CF:CS model (Constant Flux: Constant Sedimentation) (Robbins, 1978), which assumes that both the atmospheric flux of excess 210 Pb to the sediment-water interface and the sediment supply remain constant over time (Abril and Brunskill, 2014). 137 Cs is an artificial radionuclide derived from nuclear fission and is commonly used as independent tracer for validation of the 210 Pb chronologies (Smith, 2001). Its main sources in the atmosphere were the nuclear weapons testing (peaked in 1963) and in some parts of the northern hemisphere the input from the accident that occurred in Chernobyl in 1986. Dried samples for 137 Cs measurements were placed in plastic containers and counted by gamma spectrometry. Gamma emissions of 137 Cs were counted at 661.7 keV photo-peak for 24-72 h using Ortec germanium detectors (OrtecHPGeGMX-20195P and GEM-20200) calibrated against a sediment spiked with the Amersham reference standard solution QCY48A by using the same counting geometry. The detectors were coupled to a multichannel analyzer and shielded by a 10 cm thick layer of lead.

Micropaleontologic Analyses
Core samples were oven-dried at 40 • C for 1 day, weighed (∼7 g), washed on a 63 µm sieve and then oven-dried again at 40 • C for 1 day. Washed residues were separated in discrete aliquots using a precision micro-splitter and finally analyzed at the stereomicroscope for the analysis of foraminiferal content. Since it is quite difficult to morphologically distinguish between the two target species when specimens are small in size (diameter <500 µm; Hohenegger et al., 1999), the counting was limited only to those specimens that clearly displayed distinctive morphological features (following Hottinger et al., 1993) and for which the specific attribution was unambiguous. Specimens were attributed to A. lobifera (Figure 2, images 1-3) when characterized by tests with a more rounded margin, if observed in profile (p), with lobulated sutures (lobes usually visible on both sides) that were bent backwards forming an unbroken arch in dorsal view (d), while they were mainly sigmoidal in ventral view (v). On the contrary, specimens were attributed to A. lessonii (Figure 2, images 4-6) when their tests were: (I) more flattened and always characterized by an angular peripheral margin if observed in profile (p); (II) not as rounded as in similar-sized A. lobifera specimens; (III) with sutures without any folding nor complications, which sharply bent backwards forming a typical falciform arch in dorsal view (d), while they were slightly depressed in the last chambers but always maintaining the typical sickle shape, in ventral view (v).
The absolute abundances of A. lobifera and A. lessonii were calculated along the cores as number of individuals recorded per gram of dry sediment (N g −1 ).
Additionally, the absolute abundance of calcareous nannoplankton was calculated in smear-slides as total number of individuals recorded per mm 2 and then used as proxy for hydrodynamism. Since nannoplankton deposition occurs only in low-energy systems (particle size <30 µm), this independent proxy was used to qualitatively estimate if the energy at the sampling site was low enough to permit deposition of finegrained particles and, consequently, of radionuclides used for the radiometric dating.

Sea Surface Temperature (SST) Dataset
In order to verify the relation between the abundance of Amphistegina spp. and the water temperature, Sea Surface Frontiers in Marine Science | www.frontiersin.org FIGURE 2 | Photomicrographs of adult and juvenile specimens of Amphistegina lobifera (images 1-3) and Amphistegina lessonii (4-6) collected from cm 2 to 3 bsf of CORE19. v) ventral view; p) profile; d) dorsal view. Note in juvenile specimens of A. lobifera (diameter size <300 µm, images 2-3) the presence of typical lobulated sutures (arrows) and the more rounded profile with respect to juvenile specimens of A. lessonii (5-6) characterized by more flattened tests with falciform sutures without any folding nor complication. Scale bar 300 µm.
Temperature (SST) of Sicily Channel was considered. SST was extracted from 3D temperature simulated by the reanalyses of Mediterranean Sea (ref. MEDSEA_REANALYSIS_ PHY_006_009, horizontal spatial resolution ∼ 6 km), available as monthly averages covering the period 1955-2015. The quality of the MEDSEA_REANALYSIS_PHY_006_009 was assessed for the entire period by comparing results with available observations, consolidated climatological products and current knowledge of the ocean circulation. The reanalysis data were simulated by means of the assimilation of in situ temperature profile and sea level anomaly in the Mediterranean basin.
In the framework of this study, the sea temperature was analyzed at 16 m depth in the grid point of the model closest to the sampling site of sediment cores [Lat. 35 • 56 15 N; Long. 14 • 30 0 E]. Starting from available data, curves of annual average SST and annual wintry average SST were obtained. The first one was elaborated by calculating the average SST for each year since 1955-2015, and the second one by calculating the average SST from values simulated during the winter season (January, February, and March). The average annual SST and wintry SST were plotted against the absolute abundances of A. lobifera and A. lessonii. In addition, the SST annual and wintry anomalies were computed taking as a reference the SST mean over the whole period and over the winter season, respectively. In particular, the annual anomalies were obtained by subtracting the SST averaged over the period 1955-2015 from each annual mean, while the wintry anomalies were obtained subtracting to each wintry mean the SST averaged over all January, February and March in the 1955-2015 period. All relationships were analyzed with Spearman correlation analysis using the software PAST v4.01 (Hammer et al., 2001). All graphs were generated in Microsoft Excel for Microsoft 365 MSO, and modified with GIMP v.2.10.20 1 .

Radiometric Dating
Granulometric curves obtained from CORE18 (see Figure 1 for map) show that the sample meets both prerequisites for radiometric dating through 210 Pb: the mud fraction (particle size <63 µm) is continuously present with discrete abundances between 8 and 40%, and high-energy episodes (e.g., storm waves) that could have altered the normal deposition can be reasonably excluded, due to the absence of abrupt variations down core and the continuous occurrence of nannoplankton (Figure 3).
The curves of 210 Pb decay are very similar in the two cores, indicating a good replicability of the collected data. Both curves show the typical activity-depth profile, with higher activities at the core top that rapidly decrease down core (Figure 4 and Supplementary Table 3 available on-line as Supplementary Data) halving within the first 25 cm below sea floor (bsf). This interval was used to calculate a constant Sediment Accumulation Rate (SAR) of about 0.22 cm yr −1 (Figure 4A); the derived age model provides an estimated time interval of about 4.5 years for each centimeter of sediment. Unfortunately, no 137 Cs was recorded in either core, thus an independent validation through this method was not possible. While the absence of the Chernobyl peak of 137 Cs is common in sediments collected in the southern part of the Mediterranean Sea due to the dispersion pattern of 137 Cs fallout that followed the accident, it is surprising to have no signal of nuclear bomb experiments. Nevertheless, the absence of 137 Cs in both cores supports the finding that in this area the 137 Cs supply is negligible, as also reported by other works (e.g., Hassen et al., 2019). 1 https://www.gimp.org/ Amphistegina Content Down-Core Along the studied cores, the absolute abundance of both A. lobifera and A. lessonii is characterized by a decreasing trend starting from the core top to the lowest occurrence (Figure 4 and Supplementary Table 4 available on-line as Supplementary Data). All curves show the same distribution pattern down core indicating that the two records are well replicated and have not suffered episodes of mixing by high energy events. In CORE18, A. lobifera is continuously present from cm 0-1 to cm 14-15 bsf, while in CORE19 its lowest occurrence is slightly deeper and corresponds to cm 16-17 bsf. The highest abundances were recorded in the upper portion of sediment cores. In both cores, the absolute abundances of A. lobifera decrease abruptly starting from cm 6-7 bsf down to the lowest occurrences, where the minimum values were recorded. Similarly, A. lessonii shows the highest abundances in the upper portion of the cores and an abrupt decrease starting from cm 6-7 bsf down to the lowest occurrences, which correspond to cm 14-15 bsf in CORE18 and to cm 17-18 bsf in CORE19.

Sequential Stages of Invasion and Comparison With SST Trends
Based on the age model, where 1 cm of sediment corresponds to a time interval of 4.5 years, the first occurrence of A. lobifera in Marsamxett Harbor can be dated to the mid-1940s (Figure 5). During the first decades of its colonization, up to the end of the 1980s, A. lobifera is present but with very low abundances, probably in response to environmental conditions still not completely favorable for the growth of dense populations. Starting from the 1990s, A. lobifera increases more rapidly, reaching the maximum peak of abundance between 2005 and 2010, when it was recorded in the Maltese Islands by Yokeş et al. (2007)  The first occurrence of A. lessonii in Marsamxett Harbor can be dated at the end of 1930s. Although with lower abundances than A. lobifera, A. lessonii displays a similar invasion dynamic: its population is restricted to a few specimens during the first decades, then rapidly increases in abundance starting from the 1990s, up to the first record made in 2018 (Guastella et al., 2019). Thus, the first record of A. lessonii in Malta was reported with a lag time of 80 years, while the establishment of permanent populations with a lag time of 30 years (Figure 5).
From the comparison of the increasing abundances of A. lobifera and A. lessonii recorded along the sediment cores with the annual average SST measured during the last 60 years, a strong similarity emerges (Figure 6 and Supplementary Table 5 available on-line as Supplementary Data). Between the 1980s and the 1990s, a progressive rise of the annual average SST is recorded, with an average overall increase of ∼ 1 • C; in the same timeframe, A. lobifera and A. lessonii abundances start to increase as well. The same pattern results when considering the annual FIGURE 3 | From left to right: photo of CORE18, granulometric curves according to Blott and Pye (2012), and curve of absolute abundance of calcareous nannoplankton down core. The content of finest particles continuously recorded without abrupt variations along the core is mandatory in order to proceed with radiometric dating.
average wintry SST, for both A. lobifera and A. lessonii. Since the 1980s, the average SST increased rapidly during the winter, exceeding 15 • C at the beginning of the 1990s and leading to an average overall rise of ∼ 0.7 • C (Figure 6). The aforementioned relation is further supported by the time series of annual SST anomalies evaluated against the 1955-2015 average, at 16 m depth in the grid point of the model closest to the sampling site of sediment cores in Malta (Figure 7 and Supplementary Table 5).
In fact, Figure 7 shows that both A. lobifera and A. lessonii increase their abundance only when the annual SST anomaly curve exhibits positive values, corresponding to the warmer phase (with temperature higher than 1955-2015 average) starting from 1985 and still in progress (Marullo et al., 2011), and finally, the pattern is similar when considering also the annual wintry SST anomaly. All correlations between absolute abundance of both foraminiferal species and annual average SST, wintry SST, FIGURE 4 | Curves reporting from left to right: 210 Pb decay, showing the typical activity profile decreasing with depth, and absolute abundances of A. lobifera (in red) and A. lessonii (in green) recorded in both cores with a characteristic decreasing trend moving down core. The boxes on the bottom report a constant SAR, respectively of 0.20 cm yr −1 for CORE18 (A) and 0.22 cm yr −1 for CORE19 (B); this last value was utilized in the applied age model, which leads to an estimated time interval of about 4.5 years for each cm of sediment (see the text for further explanations).
SST anomaly and wintry SST anomaly resulted positive and significant in both core samples, with an average correlation coefficient of 0.74 (Figures 6, 7 and Supplementary Table 6 available on-line as Supplementary Data).

The Micropaleontologic Approach as a Reliable Procedure in Studies of Invasion Dynamics
This study is one of the few works (e.g., Ribeiro et al., 2012;Lishawa et al., 2013;Albano et al., 2018;Deldicq et al., 2019;Holman et al., 2019) applying analyses and methods commonly used in Micropaleontology and Paleoecology for the investigation of a relatively recent phenomenon, such as marine bioinvasions. Here we unified the recent fossil record of two alien foraminiferal species preserved in sediment cores with field observations, radiometric dating and environmental variables, such as SST trends. We also documented a possible causality link between temperature increase and population outbreaks.
Advantages of this approach include: (I) possibility of analyzing long historical records (up to about one century) on a continuous base; (II) possibility of analyzing early stages of invasions (e.g., see Walsh et al., 2016), which are very often unknown because several alien species tend to remain overlooked for a long time after their first arrival (Carlton, 2009); (III) possibility of directly exploring the response of the target species to environmental parameters, provided that long time-series are available. More in general, the analysis of sediment cores has proven useful in detecting anthropogenic change in a variety of coastal and marine communities, by reconstructing ecological scenarios before human intervention (e.g., Yamamuro and Kanai, 2005;de Boer et al., 2013;Lin et al., 2019;Handley et al., 2020;O'Dea et al., 2020).
Thanks to the availability of continuous, direct measurements of species abundance, the first three points above are here treated with a quantitative, replicable approach, differently from several other studies where bioinvasion patterns and dynamics are reconstructed from qualitative or non-structured observations (e.g., Delaney et al., 2008;Azzurro et al., 2019) or from modeling approaches (Salihoglu et al., 2011;Kanary et al., 2014;Mellin et al., 2016), to compensate for the lack of continuous and direct measurements.
However, using the analysis of sediment core samples to unravel invasion histories has some stringent restrictions. The first and most obvious restriction is that this method applies only to taxa with mineralized remains that can persist in thanatocoenoses of sediments (e.g., Albano et al., 2018;Deldicq et al., 2019). Another challenging requirement is related to the nature of the sediment containing the target species: it should reflect both its habitat (in terms of substratum, grain-size, depth, etc.) and its real distribution through time, should be sufficiently sheltered from wave/current motion to avoid vertical mixing, free from bioturbation and free from reworked fossil specimens that, in the case of living species having fossil counterparts, could misrepresent the original first occurrence along the core record. This may present a problem when one of these conditions is not met. Moreover, we must keep in mind that faunal FIGURE 6 | Curves reporting from left to right: absolute abundance of A. lobifera (A) and A. lessonii (B) plotted against the annual SST and wintry SST simulated since 1955-2015. Beside: scatter-plots with Spearman's coefficients (r s ) relative to the absolute abundances of A. lobifera and A. lessonii recorded in CORE18 and CORE19 compared to the annual average and wintry SST (p-values of correlation indicated in brackets as follows: *p < 0.05; **p < 0.01; ***p < 0.001).
variations observed down core strongly depend on the sediment accumulation rate (SAR): a single centimeter of sediment can correspond to multiple years, thus only in the case of very high SAR we can reconstruct annual or seasonal time series.
In spite of the restrictions mentioned above, this approach deserves to be included in the toolkit for studying invasion dynamics by these types of marine organisms, for which literature data are poorest and which are most likely to meet all the required conditions.

Sea Warming Trends and Range Expansions
Sea surface warming in the semi-enclosed Mediterranean basin has been largely documented and predicted, based on satellite measurements coupled to complex ocean models (Pastor et al., 2017;Sakalli, 2017;Macias et al., 2018). Analysis of decadal SST trends has documented a consistent warming increase starting from 1980 to 1983, which could be part of a 70-year variation linked to the Atlantic multi-decadal oscillation (Marullo et al., 2011;Pastor et al., 2017). Despite being uneven across time and space, the warming trend of the Mediterranean Sea surface has notably accelerated during the last two decades (Pastor et al., 2017;Sakalli, 2017). Since a global increase in sea temperatures has been reported to promote the redistribution of marine biodiversity, supporting the spread of thermophilic invaders (Occhipinti-Ambrogi, 2007;Marras et al., 2015;Molinos et al., 2016;Walsh et al., 2016), an accelerated warming rate should be taken into serious account when studying relationships between invasion processes and climate change, especially in the semienclosed Mediterranean Sea where the climatic signals can be amplified (Pastor et al., 2017). In fact, there is evidence that sea warming has been progressively shifting the horizontal and vertical distribution of Mediterranean marine species (Lejeusne et al., 2010;Marbà et al., 2015), and favoring tropical alien species to the detriment of temperate native species (Raitsos et al., 2010;Azzurro et al., 2019).
Our results, derived from historical records preserved in radiometrically-dated sediment cores, provide a unique opportunity to unravel the invasion history of cryptic invaders FIGURE 7 | Curves reporting from left to right: absolute abundance of A. lobifera (A) and A. lessonii (B) plotted against the annual SST anomaly and wintry SST anomaly evaluated over the period 1955-2015. Beside: scatter-plots with Spearman's coefficients (r s ) relative to the absolute abundances of A. lobifera and A. lessonii recorded in CORE18 and CORE19 compared to the annual and wintry SST anomalies (p-values of correlation indicated in brackets as follows: *p < 0.05; **p < 0.01l; ***p < 0.001). The beginning of the exponential growth of both populations starts around the 1990s, when the warmer phases (in orange) abruptly increase in intensity and frequency. and to follow the 70 year-long dynamics of two tropical species that have been exposed to the progressive heating of the central Mediterranean Sea. We document that these two target species, the Indo-Pacific A. lobifera and the cryptogenic A. lessonii, reached Malta at the beginning of the 1940s, several decades earlier than their first records (Yokeş et al., 2007;Guastella et al., 2019). However, they did not develop dense populations for at least 50 years. Protracted lag times between initial introduction and population explosion have been reported for other marine alien species, such as the mussel Brachidontes pharaonis (Fischer, 1870) in the Mediterranean Sea (Rilov et al., 2004), the barnacle Austrominius modestus (Darwin, 1854) in the North Sea (Witte et al., 2010) and the cladoceran (water flea) Bythotrephes longimanus (Leydig, 1860) in Lake Mendoka, United States (Spear et al., 2020). Such delayed population outbreaks may be due to multiple independent introduction events occurring along time, that may enrich the population gene pool and provide genotypes more fit for the local conditions (Dlugosch and Parker, 2008;Handley et al., 2011), or to changes of the receiving environment (e.g., increased vulnerability, relaxation of biotic pressure, more favorable abiotic environment, etc.), resulting in population growth of the alien species (Sakai et al., 2001;Crooks, 2005;Walsh et al., 2016).
In the case of A. lobifera and A. lessonii, a possible explanation for the delayed population outbreak could be the altered seasonal patterns of sea temperature triggered by global climate change. For both populations, in fact, the beginning of their exponential growth started only after A.D. 1990, and more pronouncedly after A.D. 2003, when the wintry SST repeatedly exceeded 15 • C (Figure 6) and the warmer phase abruptly increased its intensity and frequency along the timeline (Figure 7). Indeed, laboratory experiments and direct observations have demonstrated that 15 • C represents a thermal threshold for both species, which are unable to calcify the tests at temperatures below that value (Titelboim et al., 2019). In particular, A. lessonii displays higher sensitivity than A. lobifera, and this fact could explain why the latter has so far achieved a higher invasion success in the Mediterranean basin (Titelboim et al., 2019). Additionally, A. lobifera collected from the eastern Mediterranean has a similar thermal response to high temperatures than A. lobifera from the Red Sea, showing that the thermal tolerance is retained during its invasion (Schmidt et al., 2016a). On the other hand, laboratory experiments have shown that higher temperatures (>32 • C) are better tolerated by A. lessonii than A. lobifera (Schmidt et al., 2016b;Titelboim et al., 2019). This implies that in the future, according to the predicted Mediterranean warming (Sakalli, 2017;Macias et al., 2018) and the different tolerance ranges of the two species, A. lessonii could be favored, overcoming A. lobifera in relative abundance.
Historical data suggest that such switch is not unlikely. The southern portion of the Sicily Channel, an area characterized today by mean annual SST varying between 19 and 21 • C (Drago et al., 2010), during Late Pliocene (3.1-2.51 Ma), was warmer on average 6-7 • C than modern conditions Plancq et al., 2015;Tzanova and Herbert, 2015). In that epoch, Mediterranean amphisteginid populations were strongly dominated by A. lessonii, with minor abundances of the congeneric species A. targionii and A. gibbosa (Di Bella et al., 2005). Their test remains were so abundant on the sea floor after death to form biogenic deposits named Amphistegina-rich beds, over 1 m thick and cropping out all over the Mediterranean basin (Di Bella et al., 2005;Caruso and Cosentino, 2014;Cau et al., 2019). Amphistegina lessonii probably developed its acclimation ability earlier, thanks to a longer colonization history during the warmer Pliocene Mediterranean, at least since ∼ 5.3 Ma, when the species recolonized the basin after the abrupt foraminiferal extinction following the Messinian salinity crisis (Hayward et al., 2009). This ability could advantage the species even in the future, taking over from the congeneric Indo-Pacific A. lobifera. On the other hand, A. lobifera is characterized by a high dispersal capacity, and can modify its geographic range in response to global warming without requiring further adaptation (Prazeres et al., 2020).
This fact, however, might be considered realistic only if invaders spread and acclimate to new recipient environments fast enough to keep pace with climate change (Hiddink et al., 2012;Marras et al., 2015;Molinos et al., 2016). Based on a literature search and meta-analyses, Sorte et al. (2010) estimated that the poleward shift of marine species in response to climate change happens at an average rate of ∼ 19 km yr −1 , which is an order of magnitude faster than the estimated range shift of terrestrial species. According to a recent species distribution model (Guastella et al., 2019), A. lobifera is spreading in the Mediterranean Sea at a rate of 13.2 km yr −1 . However, this estimate was based on the known sequence of first records of the species along the eastern and central Mediterranean coasts. Here, we demonstrate that a serious bias exists in the assessment of the "first arrival" of A. lobifera, showing that in Malta the species went overlooked for as much as 60 years. Additionally, it could be assumed that this bias also exists in other localities around the Mediterranean Sea, because our results precede the first record of A. lobifera in the entire Mediterranean, making all first records of this alien species presumably incorrect, as well as questioning hypothesis of spreading due to ichtyochory by lessepsian rabbitfish (Guy-Haim et al., 2017). On the other hand, we show that the species arrived in the central Mediterranean at least in the mid-1940s, that is only 70 years after the opening of the Suez Canal occurred in 1869. Malta is about 3,200 km far from Port Said, the outlet of the Suez Canal (as calculated along the North-African coasts), a distance that both A. lobifera and A. lessonii have covered in about 70 years with a spread rate of ∼ 45 km yr −1 . This value is consistent with the spread rates estimated for other Erythraean invaders that move in response to climate change (referred to interquartile range of 25-75%: 17.0-49.8 km yr −1 for continental shelf dispersal and 12.6-35.1 km yr −1 for straight-line dispersal; Hiddink et al., 2012).

CONCLUDING REMARKS
Accurate data on the arrival and temporal dynamics of alien taxa are essential to understand invasion patterns and their drivers, and hence to design and implement effective measures. Such data are often lacking, especially for small-sized taxa that often remain undetected or are only sporadically recorded until their populations reach outstanding densities. The present work demonstrates how, in the case of taxa such as foraminifera with mineralized remains that persist in sediments after their death, micropaleontologic analysis of radiometrically-dated sediment cores can be used to reconstruct invasion histories, accurately determining the true date of introduction and subsequent temporal changes in abundance. In turn, the population abundance time-series can be linked to changes in environmental parameters. Using this approach, we unraveled the invasion history of two foraminiferal species in the Central Mediterranean Sea, reporting a considerable lag-time between their arrival and their first record, and highlighting how inadequate knowledge can lead to misleading or incorrect conclusions (e.g., expansion rates and vectors).
Finally, by relating the temporal dynamics of the two Amphistegina populations with trends in sea surface temperature, we also show that a link exists between sea warming and population density of both species. This suggests that temperature could represent an important driver of invasion patterns, with sea surface warming having the potential to trigger population outbreaks of these thermophilic alien species. We therefore propose the micropaleontologic approach adopted in the present work as an important addition to the toolkit for studying bioinvasions, serving to elucidate invasion histories, identify the main drivers of invasion success, and ultimately better inform management decisions.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
NM and AM developed the idea. AC, JE, and CC collected and sampled the cores. RG and MC performed the micropaleontologic analyses. LL carried out the radiometric dating calculations. RL provided and elaborated sea surface temperature data. All authors contributed to interpretation of results and discussions and wrote parts of the manuscript.

FUNDING
This work was part of a Ph.D. project of the University of Pavia (RG) financially supported by the MIUR (Italian Ministry of Education, University and Research) and by FRG-2018 and FFABR funds (to NM and AM, University of Pavia) and by R1D14-PLHA2010_MARGINE (to AC, University of Palermo). JE received financial support from the University of Malta's Research Fund. Two reviewers have contributed to improve an earlier version of the manuscript.

ACKNOWLEDGMENTS
We acknowledge Sonia Albertazzi (CNR-ISMAR) for help with 210 Pb and 137 Cs measurements, Maria Pia Riccardi (University of Pavia) for photos at the stereomicroscope, and Nadia Pinardi (University of Bologna) for facilitating access to SST data. The reviewers are gratefully acknowledged for their constructive comments.