Competitive dietary and morphological divergence among co-occurring species are fundamental aspects of ecological communities, particularly on islands. Cabo Verde (~570 km west of continental Africa) hosts several endemic reptiles descended from common ancestors, with sympatric species exhibiting wide morphological variation and competing for limited resources. To explore the mechanisms of resource partitioning between coexisting species, DNA metabarcoding was used to compare the diets of large and small skinks, Chioninia vaillantii and Chioninia delalandii in sympatric and allopatric contexts on Fogo Island and in a more competitive context on the small and resource-poor Cima Islet. The morphological variation of all populations was also examined to test the character displacement hypothesis and to compare the effect of different competitive scenarios. Results showed significant differences in diet and linear measurements between species and populations. The two sympatric populations of C. delalandii on Fogo and Cima showed similar changes in head morphology compared to the allopatric population, supporting character displacement. The effect of higher competitive pressure on Cima was evidenced by the increased morphological and dietary variation observed. This study demonstrates how sister species develop dietary adaptations/ morphologies to maintain stable coexistence, especially in highly competitive scenarios, providing useful insights for effective conservation strategies.

Study area

Fogo Island and the Rombo’s Islets are located in the Leeward group of Cabo Verde (Figure 1B). Fogo has 476 km², reaching the highest altitude of the country at 2829 m. The landscape is characterised by the active volcano, last erupted in 2014, and large agricultural and pasture areas within its wide caldera. Its climate varies from dry – in the lower areas – to sub-humid – in higher altitude areas. In 2020, it was declared a biosphere reserve by UNESCO due to its unique ecosystem and biodiversity.

The uninhabited Rombo’s are composed of two main islets, Grande and Cima, surrounded by the smaller islets of Luiz Carneiro, Rei and Sapado (<0.25 km²). These are classified as Integral Nature Reserve and constitute Important Bird Areas. Cima, approximately 15 km west of Fogo, is a dry flat islet with about 1.5 km² and a small hill of approximately 77 m. This islet is a key biodiversity hotspot for seabird species. However, its biodiversity is severely affected by Mus musculus Linnaeus, 1758 increasing populations.

Study species

Delalande’s skink C. delalandii (Cd hereon) is medium-sized (adults between 52–92 mm snout-vent length, SVL), and easily recognized by the presence of a black dot on the axilla and yellow eyelids. It is often found in stonewalls and beneath rocks within agricultural and livestock areas. On Fogo, it is also found in urban areas, including at higher altitudes (Figure 1C;). Despite its abundance and conspicuousness, little is known about its ecology, except for anecdotal reports of omnivory. It is classified as Least Concern according to the International Union for the Conservation of Nature (IUCN) Red List of Threatened Species, and as Low Risk on the National Red List at a national and main islands level, but Data Deficient on Rombo’s.

The Vaillant's skink of Fogo C. vaillantii xanthotis (Cv hereon) is large-sized (adults between 87.5–105 mm SVL), and distributed on Fogo and Rombo’s, more precisely on Cima [23]. It can be identified by its three well-contrasted stripes along the body and the bright yellow scales around the ear opening. On Fogo, it is restricted to the north-western side (Figure 1C) and less abundant than Cd. It is typically found in stonewalls at high altitude in vegetated humid habitats. It is viviparous and its diet is thought to be composed mostly of insects and plants. It is classified as Endangered according to the IUCN and as Undetermined on a national level and on Fogo, and as Data Deficient on Rombo’s in the National Red List.

Sampling

Sampling took place on Fogo and Cima in May 2019, with additional samples of Cd collected on Fogo in September 2021 in the allopatric area. On Fogo, 20 samples of Cd were collected in sympatric areas and 23 in allopatric areas, and 17 of Cv in sympatry. On Cima, a total of 10 and 17 faecal samples of Cd and Cv were collected, respectively. Specimens were captured using either hoop fishing or bucket traps with bait (banana and mango). Faecal samples were collected by abdominal massage and immediately preserved in 96% ethanol tubes. These were stored at -20ºC until further processing. Individuals were sexed, photographed and geolocated using a GPS device. Before release, morphological linear measures of both species were taken to the nearest 0.1 mm using digital callipers (Figure 2): snout-vent length (SVL), trunk length (TL), head length (HL), width (HW), and height (HH), eye-snout (ES) and eye-ear distance (EE), eye diameter (ED), and front (FLL) and hind limb length (HLL). Animals were placed in separate opaque tissue bags while others were processed to reduce stress.

In addition, pitfall traps with water were placed close to the areas where skinks were captured to collect invertebrates to build a DNA reference database of available food items. Captured specimens were preserved in falcon tubes with 96% ethanol until examination under a magnifying lens with an assembled camera. Individuals were sorted into different high-level taxonomic groups based on morphological identification and photographed for higher-resolution identification by experts.Diet analysis

Library preparation

DNA from faecal samples was extracted using the Stool DNA Isolation Kit, following the manufacturer's instructions. Each batch of 23 samples included an extraction negative control. Extracted DNA was amplified using previously validated markers [47] for three diet groups: invertebrates with the IN16STK-1F-mod/IN16STK-1R-mod primers, targeting the mitochondrial 16S gene; plants with g/h primers targeting the short P6-loop of chloroplast trnL (UAA) gene and vertebrates with the 12sv5F/12Ssv5R primer pair, targeting V5-loop within the mitochondrial 12S gene. To prevent the amplification, with the 12S primers, of host DNA, a blocking primer was designed (5'-TCCTCTAGGTCGGTATGGGGCACCGCCA(C3 spacer)-3'). For this, 12S sequences of both Chioninia species were obtained from GenBank and aligned with sequences from birds, fish and reptiles known to occur in Cabo Verde. Its effectiveness for Cd and Cv was validated by comparing PCR amplification with and without the blocking primer, using previously extracted DNA from tail tissues. PCR products (see details in Supplementary material) were cleaned using Agencourt AMPure XP beads, quantified using Nanodrop, diluted to 15nM, and pooled for each marker. The three libraries were then quantified using qPCR and pooled to obtain a final library at 4nM, that was sequenced using the MiSeq Reagent Kit V2 for an expected average of 20,000 paired-end reads per PCR product.

Invertebrate DNA for the reference library was extracted from legs or whole individuals (when very small), using saline extraction methods following and amplified using two primer sets. The IN16STK-mod primers, which follow the previously described methods and allow matching to dietary items, and the LCO1490/HC02198 primers, which target the animal barcode region of cytochrome oxidase I (COI). PCR conditions used for the latter are described in previous studies. Barcode sequencing was undertaken to assist with taxonomic assignments.

Bioinformatics filtering

Sequences were bioinformatically processed using the OBITtools software, which in summary performs the alignment of the sequences obtained and the filtering of PCR/sequencing errors to obtain molecular operational taxonomic units (MOTUs; see details in Supplementary material). The LULU R package was used to combine MOTUs with >84% similarity, removing remaining errors. Based on PCR blank counts, PCR products with less than 2000 reads were removed and within the remaining samples, MOTUs representing less than 1% of the total reads of that PCR product were also excluded. The remaining MOTUs were then compared against the NCBI Nucleotide Database and our reference collection, using the BLAST+ software, and classified to the lowest possible taxonomic rank. If a sequence matched more than one species with equal similarity values all belonging to the same genus, a genus-level assignment was considered. When a sequence matched more than one species or genus with similar probabilities, only species or genera known to occur in the archipelago were considered, otherwise a higher ranking would be attributed (e.g., family). If different MOTUs corresponded to the same taxon, a number would be attributed to each one of them (e.g., Araneae_1; Araneae_2). MOTUs detected in extraction blanks, considered as contaminations or of sampling baits (banana and mango) were removed from the corresponding batch of samples.

Data analysis

The frequencies of occurrence (FO) of each of the three diet item groups (plant, invertebrate, vertebrate) were estimated for the sympatric populations from Fogo (Cd.F.s and Cv.F.s from hereon) and Cima Islet (Cd.c.s and Cv.c.s from hereon), and the allopatric population of Cd from Fogo (Cd.F.a from hereon). To test for differences in proportions of each diet item group between populations, FOs were compared using two-proportion z-tests with the function prop.test in R software version 4.3.2.

Family-level identification was used for the analysis due to the inability to achieve species-level taxonomic resolution for several MOTUs. To visualize the dietary networks at the family level for the five populations, bipartite networks were generated using the function plotWeb from bipartite package v2.19. Adjacency plots were calculated at the order level using the function metaComputeModules and visualized using plotModuleWeb. Rarefaction curves were used to compare resource richness at the family level among the five populations. Richness values were estimated for the double of species with the lower sample size and differences were considered significant if the 95% confidence intervals did not overlap. The functions iNEXT and ggiNEXT from the INEXT package were used to perform and visualize this analysis, respectively.

Dietary niche overlap within the two sympatric populations (Cd.F.s vs. Cv.F.s and Cd.c.s vs. Cv.c.s) was calculated using Pianka’s index at the family level. This index ranges from 0 to 1, where 0 indicates no common resources and 1 indicates complete resource overlap. To test whether the dietary niche overlap differed significantly from what would be expected by chance, null models were generated based on the RA3 randomization algorithm, generating 10,000 null matrices that were compared to observed data, using the EcoSimR package. Non-metric multidimensional scaling (NMDS) was used to evaluate prey-use differences among populations. Samples were ordinated in a two-dimensional space according to diet dissimilarity, using the function metaMDS from the vegan package and based on a Jaccard distance matrix.

Finally, generalised linear models (GLMs) were performed to test for significant differences in diet composition at the MOTU, family and order levels in the following contexts: within sympatric species (Cd.F.s vs. Cv.F.s and Cd.c.s vs. Cv.c.s); between species from different islands (Cd.F.s vs. Cd.c.s and Cv.F.s vs. Cv.c.s), and between Cd allopatric populations (Cd.F.s vs. Cd.F.a). GLMs for multivariate presence/absence data were fitted by implementing the manyglm function of the mvabund package. The complementary log-log distribution provided the best fit and was used for subsequent tests. The significance of the GLMs was tested using the anova.manyglm function that was then implemented with the argument p.uni= 'adjusted' to perform univariate tests for identifying prey items responsible for differences in diet composition between species and among populations.