To assess the proposed concepts of Solum, Soror, and Moliones islands and their effect on marine evolutionary and biogeographical dynamics of shallow-water organisms, we selected a total of 324 islands for study (see Fig. 6). A key condition we imposed is that the selected islands were not severely impacted by ice during glacial times. Consequently, they are located in temperate and tropical climates that range in latitude from 47.74 °N (Ostrov Rasshua Island, Kuril Archipelago, Okhotsk Sea) to 20.19 °S (Aneityum Island, Vanuatu Archipelago). Further details can be found in Suppl. material 1: tables S1–S3. Within archipelagos, the selection of islands was random, as we did not use all the islands from the archipelagos sampled in this study. However, islands were not randomly picked from a table containing all islands on our planet, as this would not be feasible, considering the estimated 20,000 islands existing today with an area of over 1 km². Moreover, island selection fulfilled three important criteria. Firstly, they should be representative of a broad geographic distribution, with islands/archipelagos located in the Pacific and Atlantic Oceans, the Mediterranean Sea, the Japan Sea, the Okhotsk Sea, the Banda Sea and the Solomon Sea. Secondly, the dataset should include different types of islands, comprising ocean island volcanoes, island arc volcanoes, uplifted coral platforms, and continental shelf islands. Lastly, and for statistical reasons, after analysing the bathymetric data (see details below), we did our best to have at least 30 islands/groups of islands from each of the proposed categories (i.e., Solum, Soror and Moliones islands; see Suppl. material 1: tables S1–S3). This condition was fulfilled for Solum (50 islands) and Moliones (32 groups of islands) but not for Soror (20 groups of islands), despite our best efforts to do so.

The insular littoral area (ILA) of Solum islands or groups of islands (in the case of Soror and Moliones) was calculated for the last 150 ka (thousands of years). Following Ávila et al. (2018)m we define ILA as the area corresponding to the insular shelf surface located between the mean sea level (MSL) at a given time and the contemporary 50-m isobath. Thus, this 50-m spread in bathymetric range encompasses the intertidal down to the infralittoral biological zones. Calculations were performed using QGIS 3.28.8 (QGIS 2023). For bathymetric data, we used the GEBCO 2023 Grid, which is a global terrain model for ocean and land, providing elevation data, in meters, on a 15 arc-second interval grid (GEBCO 2023). When compared to shallow-water shipborne multi-beam bathymetry, GEBCO 2023 represents a lower resolution database, with a spatial resolution of about 1 km at the equator. Nevertheless, it is the only bathymetric/topographic dataset with global coverage and sufficient resolution to support a comprehensive comparative study such as presented here.

For each island/group of islands, we manually checked the location of the bathymetric lines of -50 m, -120 m, and -170 m, defined the limits of the polygon encircling these depths, and classified all selected islands/groups of islands as Solum, Soror, or Moliones. The -120 m bathymetric line is related to the deepest level attained during the Last Glacial Maximum (Miller et al. 2011), while the -170 m line is related to the littoral area during this period (i.e., the first submerged 50 m). In total, 31 bookmarks were employed (cf. Fig. 6), including all the islands/groups of islands that were previously selected. For each bookmark, shapefiles were used to extract bathymetric data from the global GEBCO database. A GDAL algorithm (“gdal_contour”; GDAL/OGR version 3.8.3) was employed to generate the polygon used to calculate the littoral area between the two sea levels (reference level and -50 m). In total, the GDAL algorithm was run 4,681 times (31 bookmarks × 151 ka), making it possible to determine the size (in km²) of littoral surface areas (0–50 m depth) over the last 150 ka through time steps of 1 ka, based on the sea-level curve developed by Miller et al. (2011).

For simplicity, our approach considers that no morphological changes occurred in the islands, i.e., that the present-day topography/bathymetry of an island remained the same throughout the examined time interval. Accordingly, our analyses did not account for any potential changes in topography or bathymetry resulting from erosion, sedimentation, landslides, volcanic and tectonic activity, coral reef growth, or other processes that may have occurred during the time interval under consideration. Similarly, we did not take into account any uplift or subsidence movements in our analyses, even though it is well recognised that such isostatic movements may play a role in modifying the littoral area of an island over a significant period of time (Ramalho et al. 2013; Quartau et al. 2016, 2018; Ávila et al. 2019).

We followed an iterative process to classify all ILA variation curves, as depicted in Suppl. material 1: figs S1–S6. First, all ILA curves were log-transformed (Suppl. material 1: fig. S7). We calculated the polynomial curves for each of the 104 ILA variation curves used in this study (50 Solum ILA curves, 20 Soror ILA curves, and 34 Moliones ILA curves). The curve used was a 5^th degree polynomial of the type: ILA = ax⁵ + bx⁴ + cx³ + dx² + ex + f. We also fitted this 5^th degree polynomial curve against each ILA curve and calculated the value of r². A table was then constructed for each of the 104 ILA curves, with the values of the terms (a, b, c, d, e, f) for each polynomial curve, plus the calculated r² (Suppl. material 1: table S4). Then, Principal Components Analysis (PCA) was performed on Suppl. material 1: table S4 values (Fig. 7). Finally, we double-checked the results of the PCA against each ILA curve by visual inspection and corrected any misclassifications.

We tested if the use of the GEBCO 2023 grid, which has a poor bathymetric resolution of about 1 km at the equator, could compromise the fidelity of our data, as for some islands with narrow insular shelves, the calculated ILA curves could be rough speculations. To assess this, we compared the ILA curves for Santa Maria Island in the Azores Archipelago, an example of a Solum island where high-resolution bathymetry is available (8 m resolution; Ricchi et al. 2020), with those generated from GEBCO 2023. A similar comparison was done for the ILA curves of the Pico+Faial Soror islands (Azores Archipelago), for which we have a compilation of high resolution multibeam bathymetry (50 m resolution; cf. Quartau et al. 2015), with those obtained from GEBCO 2023.

Of the 324 islands surveyed, 50 fall into the category of Solum islands (as listed in Suppl. material 1: table S1), 77 islands form 20 groups of Soror islands (Suppl. material 1: table S2), and the remaining 197 are grouped into 34 Moliones islands (Suppl. material 1: table S3). The number of islands classified as Soror within an insular system varies, ranging from two (e.g., Pico and Faial, collectively known as “Laurinsula” (Ávila et al. 2018), in the Azores Archipelago) to a maximum of 19 islands located in the Aegean Sea (Andros, Mykonos, Tinos, Syros, Rineia, Strongyli, Despotiko, Antiparos, Paros, Naxos, Donousa, Keros, Schoinoussa, Irakleia, Ios, Kato Koufonisi, Ano Koufonisi, Sikinos, and Folegandros); see also Suppl. material 1: table S2. The number of islands within Moliones systems range from two in the Caribbean Sea (e.g., Puerto Rico and Virgin Islands) to a maximum of 34 islands also in the Caribbean Sea (Grenada, Caille Island, Ronde Island, Les Tantes, Diamond Island, Large Island, Frigate Island, Saline Island, White Island, Carricou, Mabouya, Sandy Island, Petit Dominique, Fota, Petit Martinique, Petit Saint Vincent, Palm Island, Union Island, Petit Tabac, Jamesby, Petit Tameau, Baradal, Mayreau, Catholic Island, Canouan, Savan Island, Petit Mustique, Mustique, The Pillories, All Awash Island, Baliceaux, Battowia, Isle A Quatre, and Porth Elizabeth); cf. Suppl. material 1: table S3.

Four main types of ILA curves were detected using Principal Component Analysis, depicting changes in the ILA in response to sea-level variations during glacial-interglacial cycles. Islands that conformed to the SLS prediction, as their littoral area attained maximum values during interglacial periods, were classified as I-type (Table 1). In contrast, if the littoral area reached its maximum values during glacial periods, the islands were classified as G-type. “I” stands for Interglacial and “G” stands for Glacial. The PCA detected four main ILA curves: two for I-type islands (Fig. 7I–1, I-2), and two for G-type islands (Fig. 7G–1, G-2).

When the PCA results were compared with each ILA curve by visual inspection, it was clear that one more general ILA curve (I-3) was discernible. Moreover, this comparison revealed that seven out of the 104 islands/groups of islands classified according to the PCA were not well categorised and some changes were required: three ILA curves classified by PCA as I-1 were changed to I-2; one ILA curve classified by PCA as I-2 was changed to I-3; one ILA curve classified by PCA as I-3 was changed to I-2; one ILA curve classified by PCA as I-2 was changed to G-2; and finally, one ILA curve classified by PCA as G-2 was changed to G-1 (cf. Table 1).

Thus, this study identified five general types of curves (Fig. 8), which are described below:

In the following sections, we classify Solum islands, as well as Soror and Moliones groups of islands, according to the types of curves of littoral area variation for the last 150 ka based on the PCA and visual inspection of the ILA curves.

(1) Solum Islands. Suppl. material 1: figs S1–S3 display the results of the ILA variation curves in 50 Solum islands for the last 150 ka. Of these, 39 islands conformed to the SLS dynamic model predictions and were classified as I-type, whilst the remaining 11 islands were classified as G-type (Table 1). We classified 49 of the 50 Solum islands (cf. Table 1; see Suppl. material 1: table S1 for classification of each island) according to the following types of ILA curves. The ILA curve type I-1 is exemplified by 10 out of the 49 classified Solum islands (20.4% in total; Table 1): Saint Croix, Sal, Ua Huka, Tanna, Nendo, Utupua, Naunonga, Sumba, Damar, and Teun. The ILA curve type I-2 (Fig. 8I–2) is exemplified by 21 islands (42.9%; Table 1): Saba Island, Marie Galante, Dominica, Martinica, Saint Lucia, Barbados, Santo Antão, Fogo, La Palma, El Hierro, Tenerife, São Tomé, Santiago, Moho Tani, Ua Pou, Nuku Hiva, Aneityum, Erromango, Rennel Island, Malaita, and Guam. In most islands (e.g., Saba, Martinica, Barbados, La Palma, El Hierro, and São Tomé), variations in ILA also occurred during the Last Interglacial period. The ILA curve type I-3 (Fig. 8I–3) is represented by seven islands (14.3%; Table 1): Corvo, Gran Canaria, Selvagem Pequena, Selvagem Grande, Porto Santo, Príncipe, and American Samoa. The ILA curve type G-1 (Fig. 8G–1) is represented by three islands (6.1%; Table 1): Flores, Graciosa, and Saint Helena. Finally, ILA curve type G-2 (Fig. 8G–2) is represented by eight islands (16.3%; Table 1): São Nicolau, La Gomera, Isla Floreana, Etorofu (=Iturup), Ostrov Simushir, Ketoy, Izu Oshima, and Miyake.

(2) Soror Islands. Suppl. material 1: fig. S4 displays the results of the ILA variation curves of 20 groups of Soror islands for the last 150 ka. Of these, 16 groups of islands conformed to the SLS dynamic model predictions and were classified as I-type, as their ILA attained maximum values during interglacial periods. In contrast, the littoral area of four groups of islands (Gavdos+Gavdopoula islands, and the Cyclades, Lesser Sunda, and Izu archipelagos) reached its maximum values during glacial periods, and were classified as G-type (Table 1 and Suppl. material 1: table S2). We were able to classify all of the 20 groups of Soror islands (Table 1). ILA curve type I-1 is exemplified by 10 out of the 20 groups classified as Soror islands (50.0%), curve type I-2 is exemplified by four groups of Soror islands (20.0%), curve type I-3 is exemplified by two groups of Soror islands (10.0%), curve type G-1 is exemplified by one group of Soror islands (5.0%), and finally, curve type G-2 is exemplified by three groups of Soror islands (15.0%). See Suppl. material 1: table S2 for the discrimination of islands per group of Soror islands, as well as the classification of the ILA curve types.

(3) Moliones Islands. Suppl. material 1: figs S5, S6 display the results of ILA variation curves of 34 groups of Moliones islands for the last 150 ka (Suppl. material 1: table S3). We found that 27 groups conformed to the SLS dynamic model predictions, and were classified as I-type, whereas five groups of islands were classified as G-type (Table 1). We were able to classify 32 of the 34 groups of Moliones islands (cf. Table 1). Curve type I-1 is exemplified by 13 out of the 32 groups classified as Moliones islands (40.6%), curve type I-2 is exemplified by 13 groups of Moliones islands (40.6%), curve type I-3 is exemplified by one group of Moliones islands (3.1%), curve type G-1 is exemplified by one group of Moliones islands (3.1%) and finally, curve type G-2 is exemplified by four groups of Moliones islands (12.5%). See Suppl. material 1: table S3 for the discrimination of islands per group of Moliones islands as well as the classification of the littoral area curve types.

Despite its relatively coarse resolution, the ILA curves derived from GEBCO 2023 (Fig. 9A, B) are broadly consistent with those obtained from much higher resolution data (Fig. 9C, D). Although some small differences are noticeable, crucially, the overall shape and behaviour of the ILA curves remain similar (cf. Fig. 9). While future advancements in global terrain models will likely enhance the accuracy of this methodology, GEBCO 2023 currently represents the most comprehensive and up-to-date global dataset available. Accordingly, it was the most appropriate choice for the bathymetric analysis of the islands in this study.

Sea-level changes driven by glacial-interglacial cycles increased in frequency and amplitude during the last one million years of the Pleistocene (Miller et al. 2011). This period marked a transition to lower-frequency, higher-amplitude, quasi-periodic (~100-kyr) glacial variability, starting from the middle Pleistocene onwards. The five identified ILA typologies are the consequence of these shorter-lasting interglacial periods (<20 ka) compared to the longer-lasting glacial periods (100 ka). This sea level change in regime impacted and shaped coeval island shelves (Fig. 10A, C). It is well known that the width of a marine terrace is primarily dependent on the length of time that marine erosion had to carve the insular edifice (Quartau et al. 2010, 2014; Ricchi et al. 2018). Thus, the dimensions and shape of these erosional platforms are largely related to the duration of sea-level stands (Fig. 10B, D). Our ILA analysis shows that, over the last 1 Ma (cf. Fig. 10A, B), marine erosion was especially prevalent in the following depth ranges: 50–60 m, where sea level stood for 13.93% of the last 1 Ma (~139.3 ka); and 70–90 m depth where it stood for 24.88% of the time (~248.8 ka). Minor differences were detected when this analysis was restricted to the last 500 ka, with the depth-bin interval of 30 to 40 m being the third most relevant (Fig. 10D). For isostatically stable islands, that is, islands that did not experience uplift or subsidence over the last 1 Ma, it is expected that the littoral area will increase whenever the 50–60 m and the 70–90 m depth ranges are crossed by the coeval [0–50] m depth limits, which serve as interval markers for calculating the ILA. This is because sea level remained at these marks for longer periods, increasing the likelihood of larger marine terraces forming within this range. Finally, it is worth emphasising that sea level was above the present datum during the most extreme interglacials (MIS 11 and MIS 5e) for only 2.49% of the last 1 Ma (~24.9 ka). Similarly, it was below 120 m depth for the same duration (2.49%, i.e., ~24.9 ka).

Hence, in theory, the sizes of ILA are dependent on the presence or absence of: 1) submarine terraces, which are first-order topographic features imprinted on the island shelf that can enlarge areas of specific depth ranges; 2) subaerial marine terraces carved high on the shores above the present sea level, during previous interglacials, or resulting from an uplift of the island; and 3) coral reefs fringing the island that, with time and accompanying the gradual sinking of the island, will transform into barrier reefs, then into an atoll, and finally, if and when corals are not able to cope with subsidence, a guyot/submerged reef bank (Menard 1984; Ramalho et al. 2013; Montaggioni and Martin-Garin 2020). However, numerous other geological factors control the insular littoral area-depth distribution on a shelf, which, in turn, influence the diverse shapes of the ILA curves.

For oceanic islands outside the tropics, one of the most important geological factors is island age, as the longer an island is exposed to wave erosion, the larger its shelf tends to be (Quartau et al. 2010, 2014). As depicted in Fig. 2, there is a strong correlation between island age and shelf width. However, there are other significant geological processes that also control depth distribution across the shelf, such as volcanic activity, the tectonic setting of the insular edifice (whether it is subsiding, uplifting, or isostatically stable), and the depth of the insular shelf break (Ramalho et al. 2013; Quartau et al. 2015). For subsiding islands, we can add the Darwin line effect, which occurs when tropical islands reach depths at which the temperature threshold of 20 °C required for littoral corals to continue growing is crossed. All these factors are key to explain certain patterns in littoral area curves. Additionally, these geological factors may reinforce or cancel each other’s effects, making it difficult to provide general explanations.

In the case of coral reef islands, steep slopes are typically found offshore the fringing reefs or barrier reefs (Fletcher et al. 2008). Thus, a substantial increase in littoral area is expected to occur when rising sea level lead to the flooding of the large area above the edge of these reefs. Therefore, such coral reef islands conform to ILA curve type I-1 and are exemplified by the Solum island Utupua, in the Santa Cruz Islands.

For Solum islands, the General Dynamic Model of island biogeography (Whittaker et al. 2008, 2010) and its variants such as the Glacial Sensitive Model proposed by Fernández-Palacios et al. (2016) predict that land area, in conjunction with geological time, are the main drivers of terrestrial insular diversity. Therefore, for islands located within a tropical to temperate latitudinal range (50°S-50°N), the frequency of population splitting in the terrestrial realm will be higher during glacial periods (when the terrestrial area is maximal), whereas population merging and extinction frequencies will be higher during interglacial periods (when the land area is minimal) (Fig. 11B, Table 2). In contrast, Ali and Hutchinson (2014) showed that in the case of Soror islands of the Galápagos, the cycles of fusion-fission will be the main drivers of island diversity. Thus, in this case, a “species pump” mechanism (sensu Heaney 1985) will lead to a higher frequency of population merging during glacial periods as area and gene flow are maximal when Soror islands fuse. Moreover, the resulting larger population sizes will also increase dispersion rates among islands, contributing to an increasing frequency of population merging. Conversely, the frequency of population extinction and splitting will be higher during interglacial periods when Soror islands separate once again and their area becomes minimal (Fig. 11D, Table 2). For Moliones islands, Papadopoulou and Knowles (2017) showed that, due to the very shallow waters that separate them, the periods of island connectivity have been much longer than periods of isolation during the Quaternary. Thus, the duration of connectivity is the primary factor influencing the time available for speciation (which is not fully accomplished on these islands), overcoming the effect of area. A “species-vacuum” mechanism sensu Papadopoulou and Knowles (2017) is therefore proposed to explain the ephemeral, incipient genomic divergence of species during interglacials. Hence, although the frequency of splitting is higher during interglacials, speciation duration is longer than the period of island isolation. In Moliones islands, the frequency of population extinction and merging will therefore be higher during glacial periods (Fig. 8F, Table 2).

The shift from a terrestrial point of view to a marine perspective requires a complete reversal of thought, as previously argued (e.g., Ávila et al. 2018; Freitas et al. 2019; Hachich et al. 2020; De Groeve et al. 2022). As demonstrated by Ávila et al. (2019) through their SLS dynamic model of marine island biogeography, the insular littoral area is the main driver of shallow-water marine biodiversity. Here, we expand the SLS model to include Solum, Soror, and Moliones islands and, in the process, explicitly acknowledge that connectivity in Moliones islands is fundamentally different in relation to Solum and Soror islands, as the marine connectivity within Moliones island populations of benthic, shallow water organisms is much less dependent on sea level variations when compared with Solum and Soror islands. From a marine point of view, Moliones islands behave as a single unit, and this statement is supported by the geographic distribution of Conus gastropods from Cabo Verde. In Cabo Verde, all the islands are Solum, except for a single group of Moliones islands composed of two islands (São Vicente and Santa Luzia) and two islets (Branco and Raso; cf. Suppl. material 1: table S3). Ávila et al. (2019) provided information on the diversity of the Cabo Verdean Conus species (their Suppl. material 1: table S2), reporting a total of 85 species from the archipelago, of which 70 were considered as SIME (single island marine endemics). This list was the subject of recent taxonomic re-evaluations made by the WoRMS editors, and many species were synonymized. Thus, 52 is the current number of valid Conus species accepted for Cabo Verde, 48 of which are endemic to the archipelago (cf. Table 3). All 48 Cabo Verdean endemic species have a non-planktotrophic mode of larval development, and 38 are considered as SIME, each of these species with a geographic distribution restricted to a single island, sometimes even confined to a bay within an island. In contrast, there are two species (Conus pseudonivifer Monteiro, Tenorio & Poppe, 2004 and Conus venulatus Hwass in Bruguière, 1792) that, although being endemic to the archipelago, occur on four of the islands (Sal, Boavista, Maio, and Santiago). The interesting point here is the anomalously high number of six species (Conus bellulus Rolán, 1990; Conus decoratus Röckel, Rolán & Monteiro, 1980; Conus grahami Röckel, Cosel & Burnay, 1980; Conus lugubris Reeve, 1849; Conus navarroi Rolán, 1986; and Conus saragasae Rolán, 1986) that occur in the single Moliones group of islands, taking into account the notorious tendency for almost no dispersal of these endemic Conus species. In total, 11 species of Conus endemic to Cabo Verde are reported from this Moliones group of islands (cf. Table 3): the above-mentioned six species that occur simultaneously in São Vicente and in Santa Luzia; three other SIME Conus from Santa Luzia Island, and two further SIME Conus from São Vicente Island. Notably, there is also a group of five Conus species (cf. Table 3) that are shared between Boavista and Maio (both islands classified as Solum), the reason for this being the existence of a very large and shallow seamount (João Valente Bank), which highly increases the connectivity between these two Solum islands.

For Solum and Soror islands, the marine population connectivity will be lower during interglacial episodes, while isolation plays a key role in preventing gene flow among terrestrial insular populations (Papadopoulou and Knowles 2017). However, isolation is less prevalent in the marine realm (Ávila et al. 2018, 2019), primarily due to the different means of subaqueous dispersal by marine biota (larvae and eggs, as well as rafting on floating substrata by eggs/juveniles/adult organisms), and the greater connectivity of the aqueous medium in which marine organisms reside. Following the SLS dynamic model, variation in littoral area is the most prominent variable explaining the rates of marine population splitting, extinction, and merging (Ávila et al. 2019). Accordingly, across all types of islands (i.e., Solum, Soror, and Moliones), higher merging rates (due to higher gene flow) and population extinction rates (due to a minimum littoral area) are expected during glacial periods. Conversely, a higher frequency of population splitting events is expected during interglacial periods (when littoral area reaches its maximum) (Fig. 11A, C, E, Table 2). According to our data, these expectations occur in 82 of the 102 analysed islands or groups of islands; cf. Table 1). Thus, the assumptions of the SLS dynamic model are sustained in 80.4% of the 102 analysed cases, with synergetic effects of small littoral area and high gene flow wiping out or merging incipient lineages more frequently during glacial periods. However, the gradient of variation in the islands’ littoral area during glacial and interglacial periods, as identified here by the ILA curves, has a bearing on the frequency of speciation processes (i.e., population splitting, extinction, and merging): Solum islands have comparatively higher differences in littoral area and frequency of processes driving speciation events (Fig. 11A), followed by Soror islands with intermediate levels (Fig. 11C), and, finally, by Moliones islands with lower levels of speciation processes taking place (Fig. 11E). Again, the endemic Cabo Verdean Conus support our predictions, with a higher number of species registered in the Solum islands (17 endemic Conus at Boavista Island, 15 at Maio, and 10 at Sal) when compared with the Moliones islands (9 endemic Conus at Santa Luzia and 8 at São Vicente; cf. Table 3).

Although Ávila et al. (2018) found that the changes in isolation caused by glacial-interglacial cycles are negligible for most Atlantic archipelagos, the amount of time available for speciation (i.e., the period of ILA isolation) is crucial for differentiating Solum, Soror, and Moliones islands. As previously stated, the ILA has an increasing effect on the connectivity between islands, being minimal in Solum islands, intermediate in Soror islands, and reaching a maximal effect in Moliones islands. In the case of Moliones islands, no significant changes in the high values of connectivity within these islands are expected in the marine realm, with biodiversity being essentially independent of glacial-interglacial cycles (cf. Figs 3–5). The importance of island connectivity is thus ranked as medium to low in the case of Solum islands, and very relevant in the case of Soror islands.

Nevertheless, as the sea level falls during the inception of a glacial episode, the littoral area of Moliones islands suddenly decreases when the sea level falls below the shallow depths that separate these islands. When this occurs, the initial terrestrial area of the island, along with the coastal region, becomes entirely subaerial, resulting in a sudden decrease in insular littoral area. Moreover, the littoral area of Moliones islands tends to persist at low levels for extended periods compared to Soror or Solum islands. This helps explain why the frequencies of population splitting, extinction, and merging are comparatively higher in Solum islands (Fig. 11A), intermediate in Soror islands (Fig. 11C), and lower in Moliones islands (Fig. 11E).

In marine ILA habitats, the process of speciation hinges significantly on the duration of isolation (see Fig. 1). Over evolutionary timescales, while complete speciation might not occur universally, genetic divergence among populations is expected, with its extent likely influenced by the type of island. Papadopoulou and Knowles (2015) demonstrated that the high levels of genomic divergence detected in Caribbean crickets across the Virgin Islands, which are associated with repeated connectivity cycles, do not translate into species diversification. Taking all of the above into consideration, we predict a genetic population divergence-pump effect for the ILA habitats, with higher genetic population diversity in Moliones islands, and progressively, comparatively smaller genetic diversity values in Soror islands, and even smaller in Solum islands. In contrast, we predict the frequency of population-exclusive haplotypes to be higher in Solum islands than in those experiencing fusion-fission cycles (i.e., Soror- and Moliones islands). This expectation arises from the assumption that dispersal and gene flow are expected to be lower in Solum islands due to prolonged isolation. This extended period of isolation provides ample time for new genetic variants to emerge within the population and gradually increase in frequency over time, ultimately becoming exclusive to that population. At the other extreme, Moliones islands are always connected through time, continuously exchanging individuals and likely working either as a single population or as a marine metapopulation (Kritzer and Sale 2006) during interglacial periods. This connectivity reduces the likelihood of the appearance of new exclusive haplotypes. Soror islands constitute an intermediate case, where connectivity and gene flow levels follow a pulse-like trend during glacial-interglacial cycles, alternating between continuous and fragmented habitats. During interglacial periods, these are separated by unsuitable, deep-water habitats which hamper connectivity among populations, causing a certain degree of isolation and allowing for the occasional appearance of new haplotypes within each island or population. The fusion of islands during glacial periods ensures the physical continuity of shallow-water habitats and allows individuals to be exchanged, thus increasing overall genetic diversity. When the next interglacial period starts, the habitat undergoes fragmentation once more, resulting in the interruption or, at the very least, reduction of gene flow.

Bernardi et al. (2014) studied the phylogeography of the Galápagos’ reef fishes, whereas Almada et al. (2017) studied the phylogeography of the ballan wrasse Labrus bergylta Ascanius, 1767 across its distribution range, and Baptista et al. (2021) studied the population structure of the intertidal gastropod Cingula trifasciata (J. Adams, 1800) in the Azores Archipelago. Their results support our predictions. The six insular reef fish species studied by Bernardi et al. (2014) [Lythrypnus gilberti (Heller & Snodgrass, 1903), Dialommus fuscus Gilbert, 1891, Lepidonectes corallicola (Kendall & Radcliffe, 1912), Malacoctenus zonogaster Heller & Snodgrass, 1903, Stegastes beebei (Nichols, 1924), and Stegastes arcifrons (Heller & Snodgrass, 1903)] consistently show higher frequency of exclusive haplotypes in Solum islands (e.g., San Cristobal, Española, Marchena, and Wolf), when compared with the Soror islands (Fernandina+Isabela; cf. fig. 2 of Bernardi et al. 2014). In the Azores, the populations of Labrus bergylta studied by Almada et al. (2017) also show higher frequency of exclusive haplotypes in Solum islands (e.g., Santa Maria, Corvo), when compared with the group of Soror islands (Faial+Pico (although Pico was not sampled); cf. figs 2, 3 of Almada et al. 2017). In a similar manner, the populations of the intertidal gastropod studied by Baptista et al. (2021) also show higher frequency of exclusive haplotypes in four out of the five Solum islands sampled, when compared with the group of Soror islands (Pico+Faial (although Faial was not sampled); cf. fig. 3 of Baptista et al. 2021).