Methods
A total of 103 adult lizards from five different species were sampled in
September 2020: Podarcis bocagei (n = 33), Podarcis
lusitanicus (n = 8), Podarcis siculus (n = 20), Podarcis
virescens (n = 22) and Teira dugesii (n = 20). All these
lacertid species are small-sized, diurnal, mostly insectivorous, and
exhibit sexual dimorphism, with males usually being larger than females.
Podarcis bocagei and P. lusitanicus were collected from a
semi-natural habitat in Moledo, northern Portugal (Fig. 1d)
(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 lots of 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
[23]. 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 [24].Podarcis siculus and P. virescens were collected in
Lisbon, at Parque das Nações (Fig. 1a, b) (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. WhileP. virescens is native to this location, P. siculus is an
invasive species introduced about two decades ago [25]. Its
plasticity in spatial use of habitat, morphology, behaviour, and
versatile diet explains its successful colonization of multiple
locations outside its native range [26-29]. This invasive species
can present a more versatile diet, as it can also consume fruits and
nectar [30] and have a more herbivorous diet [e.g. 26], whileP. virescens is known to be insectivorous and to feed mainly on
individuals of the class Arachnida and the orders Hymenoptera,
Hemiptera, Coleoptera and Diptera [31]. Finally, we collectedTeira 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. dugesiioccupies an anthropogenic area, although less busy, close to railway
tracks with limited vegetation cover (Fig. 1c). This species is thought
to have been accidentally introduced via transport ships from Madeira
Island three decades ago, in 1992 [32]. Teira dugesii feeds
preferentially on insects but also on small fruits [33].
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.
Sampling of cloaca exudate was preferred to using fecal samples as it is
known to provide a closer representation of the microbial communities
residing in the lower gut of lizards [22].
After the 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). These protocols and research were approved by the
Portuguese Institute for Conservation of Nature and Forests (ICNF)
(License 703/2021/CAPT).
In the laboratory, we extracted the DNA 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 and PCR controls, following the protocol of Kozich et
al. [34]. Amplicons were sequenced in a single Illumina MiSeq run.
All analysis were performed using the R Software v.4.1.1. Raw FASTQ
files were denoised using the DADA2 pipeline [35] in R with the
parameters for filtering and trimming: trimLeft = 20,
truncLen = c(220,200), maxN = 0, maxEE = c(2,2), truncQ = 2; and the
SILVA 138 database [36, 37] 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 [38] and normalized read counts were
obtained using the negative binomial distribution implemented in DESeq2
[39]. ASVs with a count of less than 0.001% of the total number of
reads (3586752 [total number of reads]x 0.001% = 36) and that were
present in a single sample were also removed.
Bacterial taxonomic alpha-diversity (intra-sample) and beta-diversity
(inter-sample) were estimated using the phyloseq package.
Alpha-diversity was estimated using the number of observed ASVs, and the
Shannon, 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.
Statistical differences in alpha-diversity between locality/habitat,
species and sex were analyzed using a linear model (lm(alpha-diversity
~ locality + species + sex)). Given the significant
effect of locality (which also corresponded to semi-natural and
urbanized habitats) on alpha-diversity (see results section),
differences in the proportions of the most abundant taxa at the phyla
and genera level (represented by ≥ 1% on average of all sequences) were
assessed between species and sex for each locality/habitat separately
using a linear model (lm (bacterial phyla/genus ~
species * sex)). The effects of locality, species and sex on microbial
beta-diversity were assessed using permutational analysis of variance
(PERMANOVA) with 10000 permutations, implemented using the adonisfunction of the R vegan package [40] (adonis(beta-diversity
~ local + species + sex)). Correlations between
individual size and bacterial alpha-diversity were also tested using the
Pearson correlation test for each species, using the ggpubrpackage [41].
Bacterial transmission between each pair of species from sympatric
populations living in Moledo and Parque das Nações was estimated using
the FEAST software [42], 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 andP. lusitanicus, only a fraction of the individuals of P.
bocagei was included (with the most similar sex and SVL ratios to theP. lusitanicus samples as possible), following the FEAST
developers’ recommendations to avoid overestimation of transmission.