Methods
A total of 103 adult lizards from five different species were sampled in September 2020: Podarcis bocagei (males = 22; females = 9),Podarcis lusitanicus (males = 6; females = 2), Podarcis siculus (males = 13; females = 6), Podarcis virescens (males = 16; females = 6) and Teira dugesii (males = 7; females = 13).
All these lacertid species are small-sized, from about 4 cm to 6 cm snout-vent length in P. virescens , P. bocagei and P. lusitanicus , and up to 8 cm and 9 cm in the introduced species P. siculus and Teira dugesii (Arnold & Burton, 2002).Podarcis bocagei and P. lusitanicus were collected from a semi-natural habitat in Moledo, northern Portugal (Fig. 1a) (41°50’19.2”N 8°52’24.5”W), where they live in syntopy (i.e., occurrence of two species in the same habitat at the same time). This location has limited human disturbance and has extensive vegetation with natural and artificial shelters (e.g., walls of agricultural properties) that can be used by lizards. Ecological adaptation is considered a major factor favoring the isolation between these two species; P. lusitanicuslives more on rocks, while P. bocagei is ground-dwelling (Carretero et al. 2015). The diet of these two species is mainly composed by prey belonging to Hemiptera, Coleoptera, Diptera, Hymenoptera and Araneae, with minimal differences between species or sexes (Kaliontzopoulou et al. 2011). Podarcis siculus andP. virescens were collected in Lisbon, at Parque das Nações (Fig. 1b, c) (38°76’22.4”N, 9°09’44.3 W), where both live in sympatry (sharing habitat type). This is a highly urbanized area near the Tejo River, characterized by large residential and commercial areas, with considerable daily human disturbance. While P. virescens is native to this location, P. siculus is an invasive species introduced about two decades ago (González de la Veja et al. 2001). Its plasticity in spatial use of habitat, morphology, behaviour, and diet explains its successful colonization of multiple locations outside its native range (Vervust et al. 2010; Carretero & Silva-Rocha, 2015; Damas-Moreira et al. 2019; Damas-Moreira et al. 2020). This invasive species can present a more versatile diet, as it can also consume fruits and nectar (Mačát et al. 2015; Vervust et al. 2010), while P. virescens is known to be insectivorous and to feed mainly on individuals of the class Arachnida and the orders Hymenoptera, Hemiptera, Coleoptera and Diptera (Juan, 1997). Finally, we collected Teira dugesii in a nearby area in Lisbon, in the Alcantara docks, close to the city port area (38°70’33.8“N, 9°16’54.1“W). Similar to the other Podarcis spp. captured in Lisbon, T. dugesii occupies an anthropogenic area, although less busy, close to railway tracks with limited vegetation cover (Fig. 1d). This species is thought to have been accidentally introduced via transport ships from Madeira Island three decades ago, in 1992 (Sá-Sousa, 1995). Teira dugesii feeds preferentially on insects but also on small fruits (Sadek, 1981).
All individuals were captured using nooses. Lizards were carefully immobilized, avoiding any human contact with the cloaca. We quickly inserted a sterile cotton swab into the entrance of the cloaca to obtain individual microbial samples. The tips of the swabs were cut into individual tubes and stored in ice boxes in the field, and then stored at -80°C upon arrival in the laboratory. After microbial sampling, each lizard was sexed, and the snout-vent length was measured (SVL; from head to cloaca) using a digital caliper (± 0.01mm error).
In the laboratory, DNA was extracted from the swabs using the DNeasy® PowerSoil® Kit (QIAGEN, Hilden, Germany) according to the manufacturer’s instructions. DNA concentration and quality were measured with the Epoch™ Microplate Spectrophotometer (BioTek Instruments, Inc.; United States of America). DNA was shipped in dry ice to the Centre for Microbial Systems at the University of Michigan Medical School (USA) where the V4 region of the 16S rRNA gene (~ 250 bp) of the bacterial communities was amplified for each sample, along with the extraction blanks, PCR controls and a mock community (D6306 ZymoBIOMICS Microbial Community DNA Standard, Zymo Research, USA) using the primers 515F (5′-GTGCCAGCMGCCGCGGTAA-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT −3′) and following the protocol of Kozich et al. 2013. The V4 region of this gene is widely used to characterize bacterial communities in various taxa, including reptiles (e.g. Colston and Jackson, 2016; Chiarello et al. 2018). Amplicons were sequenced in a single Illumina MiSeq run using a MiSeq Reagent Kit V2 500 cycles.
All analyses were performed using the R Software v.4.1.1 (Team, 2020). Raw FASTQ files were denoised using the DADA2 pipeline (Callahan et al. 2016). After an assessment of read quality plots, the parameters for trimming and filtering were set as: trimLeft = 20, truncLen = c(220, 200), maxN = 0, maxEE = c(2, 2), truncQ = 2. The SILVA 138 database (Pruesse et al. 2007; Quast et al. 2013) was chosen for taxonomic assignment. After quality control and taxonomic assignment, sequences identified as Archaea, Eukaryota, Mitochondria, Chloroplast, as well as sequences unassigned to bacteria were removed from the dataset. An Amplicon Sequence Variant (ASV) frequency table was constructed using the R package phyloseq(McMurdie & Holmes, 2013). Normalized read counts were obtained using the negative binomial distribution implemented in DESeq2 (Love et al. 2014; McMurdie and Holmes, 2014). ASVs not only with a count of less than 0.001% of the total number of reads (3586752 [total number of reads] x 0.001% = 36) but which were also present in a single sample were removed (Appendix 1). The composition and abundance of taxa in the mock community were similar to those described by the manufacturer.
Bacterial taxonomic diversity (alpha-diversity, calculated intra-sample) and structure (beta-diversity, calculated as the dissimilarity or distance between pairs of samples) were estimated using thephyloseq and the picante packages (McMurdie & Holmes, 2013; Kembel et al, 2010). Alpha-diversity was estimated using the number of observed ASVs, the Shannon index and Faith’s Phylogenetic Diversity (PD). Beta-diversity was measured using the Bray-Curtis index and the Unifrac phylogenetic weighted and unweighted distances. Principal Coordinate Analysis (PCoA) were used to visually assess dissimilarity among groups.
First, statistical differences in alpha-diversity between localities were assessed using species as a random factor using a linear mixed effects model (lmer(alpha-diversity ~ locality + (1 | species)). Given the significant effect of locality on alpha-diversity (see results section), differences in alpha-diversity among species and between sexes were further assessed using another linear mixed effects model with locality as a random factor (lmer(alpha-diversity ~ species + species:sex + (1|locality)) using the package lme4 (Bates et al. 2014). The effects of locality and species on microbial beta-diversity were assessed using a permutational analysis of variance (PERMANOVA) with 9999 permutations, with the adonis2 function of the Rvegan package (Oksanen et al. 2013), using the formula (adonis2(beta-diversity ~ locality + species)). Since both locality and species significantly affected beta-diversity, the pairwise effects of species and sex were tested for each locality separately using the pairwise.adonis2 function (Arbizu, 2020) using the model (pairwise.adonis2(beta-diversity ~ species + species:sex)). P-values for multiple comparisons were adjusted with the Bonferroni correction. Differences in the proportions of the most abundant taxa at the phyla and genera levels (represented by ≥ 3% on average of all sequences) were assessed between species and sex for each locality separately using a linear model (lm(bacterial taxa ~ species + species:sex)). Correlations between individual size and bacterial alpha-diversity were also tested using the Pearson correlation test for each species, using the ggpubrpackage (Kassambara & Kassambara, 2020).
To further understand the levels of similarity between sympatric and syntopic species, bi-directional bacterial transmission between each pair of species from Moledo and Parque das Nações was estimated using the FEAST software (Shenhav et al. 2019), by testing the contribution of each species (source) to the microbial diversity to its sympatric congener (sink). To this end, the non-normalized ASV frequency table was used and, due to differences in the number of samples between P. bocagei and P. lusitanicus, only a fraction of the individuals of P. bocageiwere included (the ones with the most similar sex and SVL ratios to theP. lusitanicus samples as possible), following the FEAST developers’ recommendations to avoid overestimation of transmission.