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.