Sex chromosome diversification in the smooth snake Coronella austriaca (Reptilia, Serpentes)

The smooth snake Coronella austriaca is a widespread Palearctic colubrid species. The species has been the subject of several molecular and phylogeographic studies which highlighted the occurrence of distinct genetic lineages in different areas of the species distribution, but scarce cytogenetic data are currently available on the species. In this paper we present a molecular and karyological study performed with several banding, staining methods and NORFISH on samples of C. austriaca from different geographical areas (Italy and Greece) of the species distribution. The molecular and phylogenetic analysis unambiguously placed the studied samples in different clades with a clear geographical pattern. The karyotype of the two female samples studied was composed of 2n = 36 chromosomes with 16 macroand 20 microchromosomes and a mix of plesiomorphic and derivate chromosome features. All macrochromosomes were biarmed with the exception of pair 5 that was telocentric. NORs were detected on a microchromosome pair. In both females, the pair 4 was heteromorphic (and completely heterochromatic after C-banding in the Italian female), representing the first report of a ZZ/ZW sex chromosome system with female heterogamety in C. austriaca. In addition, the W chromosome showed a different morphology between the two female studied (submetacentric and subtelocentric), highlighting the occurrence of a chromosomal diversification among distinct geographical areas of the species distribution and further supporting that the species contains different diverging evolutionary clades.


INTRODUCTION
The smooth snake Coronella austriaca Laurenti, 1768, is a small sized ovoviviparous colubrid with a widespread distribution across the western Palearctic (Strijbosch, 1997). The species occurs in western, central and southern Europe reaching as far as the Ural Mountains and the Caspian Sea, up to northern Iran and western Kazakhstan. Coronella austriaca is absent from European and Mediterranean islands, except for Southern England, Sicily, the Island of Elba and the Island of Krk (Engelmann, 1993;Strijbosch, 1997;. In its large distribution range, the species have a relatively uniform morphology and only two subspecies are currently recognised: C. a. acutirostris Malkmus, 1995, from the north-western Iberian Peninsula and the nominal subspecies C. a. austriaca which occupies the rest of the geographical range of the species. A third subspecies, C. a. fitzingeri (Bonaparte, 1840) was formerly described from southern Italy and Sicily, but it has been recently synonymised with the nominal subspecies (Speybroeck et al., 2016).
The species has been the object of several molecular and phylogeographic studies (Santos et al., 2008;Llorente et al., 2012;Galarza et al., 2015;Jablonsky, 2019) and rep-resents an ideal candidate to assess patterns of intraspecific genetic diversity among different populations as well as the occurrence of different refugial areas across its distribution range. These studies evidenced that C. austriaca comprises several distinct clades showing a complex haplotype structure and deep genetic divergence. These evidences do not reflect the current taxonomy of the species and a revision is probably required to better describe its molecular and geographic intraspecific diversity.
On the other hand, the available chromosome data on C. austriaca are dated and refer to the studies by Matthey (1931) and Kobel (1967) that described the karyotype of two males from Switzerland, but no information is currently available from other, geographically and genetically distinct clades of the species. The karyotype described by Matthey (1931) and Kobel (1967) was composed of 2n = 36, with 16 macrochromosomes and 20 microchromosomes. All the macrochromosomes were biarmed, excluding the elements of 5 th pair that were telocentric. Furthermore, these studies were performed only with standard staining methods. No information is currently available on the sex chromosome system, location of nucleolar organizer regions (NORs) or other chromosome markers which would be useful for karyotype comparison among different population of the species or with other phylogenetically related colubrid species. In fact, cytogenetic inferences, especially when linked to molecular data in a phylogenetic perspective, are useful tools to detect plesio-and apomorphic states and to reconstruct evolutionary trends in the studied species (Odierna et al., 1987;Mezzasalma et al., 2013Mezzasalma et al., , 2017Mezzasalma et al., a, 2017bFuller et al., 2018).
In this work we present the results of a karyological study performed on smooth snake samples belonging to different populations and geographical areas of the species distribution (Italy and Greece) using several standard staining and banding methods and molecular cytogenetics. The aim of this study was to evidence the possible presence of a chromosome sex determination system, the occurrence and location of heterochromatic regions and of NORs. We also performed a molecular and a phylogenetic analysis on the studied samples, adding our data to those available from the literature, in order to evaluate their genetic diversity and place our chromosome data in a phylogenetic context.

Sampling
We studied preserved tissue and cell suspensions of two samples of C. austriaca obtained from existing collections and no aminal was collected during the realization of this study. Studied samples include one female from Italy (Picentini Mountains, Avellino), and one female from Greece (Peloponnese) both hosted at the Dipartimento di Biologia Evolutiva e Comparata, Università degli Studi di Napoli Federico II, since 1972 (sample numbers CA0201-CA0202). These two samples were both used in a preliminary molecular analysis and in the karyolgical study as described below. The samples used in this study were already used in previous analyses (Mezzasalma et al., 2014(Mezzasalma et al., , 2016a.

Molecular and phylogenetic analysis
In order to properly identify the studied samples, assess their genetic diversity and establish the taxonomic affinities with other available sequences on the species, a molecular analysis was performed using a fragment of the mitochondrial Cytochrome b (Cytb) gene. This mtDNA gene was chosen considering the number of available sequences for several populations of C. austriaca (Llorente et al., 2012;Santos et al., 2008Santos et al., , 2012Galarza et al., 2015;Jablonsky et al., 2019).
Genomic DNA was extracted from both samples using the standard phenol-chloroform method by Sambrook et al. (1989) and then used in the PCR amplification of a fragment of about 600 bp of the Cytb. The primers used were FORWARD: 5' AACTTCGGATCCATACTACTAA 3' and REVERSE: 3' TAAA-GATGTTAGGGGTGAATGA 5' , and the PCR parameters those reported by Mezzasalama et al. (2015a).
PCR products were sequenced on an automated sequencer ABI 377 (Applied Biosystems, Foster City, CA, USA) using BIGDYE TERMINATOR v3.1 (ABI). Sequences were blasted in GenBank and chromatograms were checked and edited using CHROMAS LITE 2.1.1 and BIOEDIT 7.2.6.1 (Hall, 1999). The best-fitting substitution model (GTR+I+G) was chosen using JMODELTEST 2.1.7 (Darriba et al., 2012), under the corrected Akaike information criterion (AICc). In the phylogenetic analysis we used the two newly determined sequences (Accession Numbers: MW861682-MW861683) and available homologous sequences from Jablonski et al. (2019), downloaded from Gen-Bank, choosing when possible, the longest available sequences for each major clade (see Fig. 1). As outgroup we used an available homologous sequence of C. girondica (AN: AF471088). Phylogenetic analysis with Bayesian inference (BI) was performed using MRBAYES 3.2.7 (Ronquist et al., 2012), with two parallel runs of 8 000 000 generations and four incrementally heated Markov chains (using default heating values), with a burn-in of 25% and sampling the chains every 500 generations. A 50% majority-rule consensus tree was retrieved from the post-burn-in samples and chain convergence was checked with convergence diagnostic values (average standard deviation of split frequencies >> 0.01, potential scale reduction factor).

Chromosome analysis
Chromosomes were obtained from tissue samples and cell suspensions using the air-drying method, as described in Mez-zasalma et al. (2019). The chromosome analysis was performed with traditional staining (5% Giemsa solution at pH 7 for 10 min) and several chromosome staining and banding techniques: Quinacrine staining, Chromomycin A 3 -methyl green staining (CMA/MG) (Mezzasalma et al., 2015b), C-banding (Sumner, 1972), and sequential C-banding + CMA + DAPI (Sidhom et al. 2020a, b). Karyotype reconstruction was performed after scoring at least five metaphase plate from each sample studied. Nucleolus organizing regions (NORs) were located following the Ag-NOR banding protocol reported by Howell and Black (1980) and fluorescence in situ hybridization (NOR-FISH) as described in Sidhom et al. (2020a), using as probe the PCRamplified and biotinylated 18S rRNA gene of the gekkonid Tarentola mauritanica. The detection of FISH signals was carried out using ExtrAvidin FITC (Sigma Aldrich), counterstained with propidium iodide (PI).
Metaphase plates were detected and recorded using an epifluorescent microscope (Axioscope Zeiss) equipped with an image analysis system.

Molecular and phylogenetic analysis
A a fragment of about 500 bp of the Cytb was successfully amplified and sequenced from the two original samples, respectively from Italy and Greece (see Sampling), used in the cytogenetic analysis. No interruptions of the reading frame were detected in either sequences. Nucleotide identity samples was about 99% between the newly sequenced sample from Italy (Monti Picentini) and an available sequence from Italy (AN: MH382919), and 100% between the newly determined sequence from Greece and an available sequence from Greece (AN: EU022647).
In the phylogenetic analysis, the final Cytb alignment contained 39 sequences and 1031 nucleotide positions. The resulting tree (Fig. 1) retrieved all the main clades and a topology similar to that reported by Jablonski et al. (2019). The relative evolutionary relationships among the main clades mostly reflect an East-West scenario of geographical diversification by distance. Our tree retrieved, four major clades, comprising several smaller subclades: Iranian 1, Sicilian + Iberian 1-3 + Western 1-2, Anatolian 1-2 +Balkan + Central European, and Iranian 2 + Transcaucasian + Eastern, respectively (Fig. 1). In our tree, the Easternmost major clade (excluding the only sample of the Iranian 1 clade), comprising the Iranian 2 + Transcaucasian + Eastern subclades, is the outgroup of all the remaining clades. The two original samples from Italy and Greece clustered unambiguously within the Western 2 and Balkan clades, respectively (Fig. 1). Statistical support values were generally high at terminal nodes while varied from high to low at deeper nodes (Fig. 1).

Chromosomal analysis
The Italian and Greek samples of C. austriaca showed a very similar karyotype composed of 2n = 36 chromosomes, with 16 macrochromosomes and 20 microchromosomes (Fig. 2). Among macrochromosomes, pairs 1, 3, 6 and 8 were metacentric, pairs 2 and 7 were submetacentric, and pair 5 was telocentric. The pair 4 was heteromorphic, carrying the female sex chromosomes, ZW. In particular, the chromosome Z was always metacentric while the W resulted subtelocentric in the female from Picentini Mountains (Italy) ( Fig. 2A) and submetacentric in the female from Peloponnese (Greece) (Fig. 2B).
Other chromosome stainings and banding techniques were performed only on the Italian sample of C. austriaca, as quality and quantity of metaphase plates of the Greek individual were adequate just for its karyotype description.
Quinacrine stained evenly both macro-and microchromosomes (Fig. 3A), while CMA/MG evidenced two microchromosomes and the telomeric regions of all macrochromosomes and the telomeric regions almost all microchromosomes (Fig. 3B). A microchromosome pair resulted specifically marked by both Ag-NOR staining and NOR-FISH, (Fig. 3C and D).
After C-banding, heterochromatic bands were hardly visible, if not absent, on most macro-and microchromosomes, with the exclusion of the W chromosome, which was completely heterochromatic (Fig. 4A), but not evidenced with CMA or DAPI (Fig. 4B, C).

DISCUSSION
Our phylogenetic inference produced similar results to those obtained by Jablonski et al. (2019), retrieving all the 14 clades previously described. The main differences concern the position of one of the major clades (comprising Iranian 2 + Transcaucasian + Eastern sub-  clades) which is not involved in a basal polytomy with the other major clades but appear to be the outgroup of all the remaining clades. This results further supports the hypothesis of an East-West geographical diversification process in C. austriaca. Similar East-West differentiation processes are also shared by different European and Palearctic reptile taxa, including distinct major snake lineages of Asiatic origin (see, e.g., Utiger et al., 2002;Nagy et al., 2004;Mezzasalma et al. 2015aMezzasalma et al. , 2018. However, the phylogenetic position and the evolutionary relationships of some clades (such as the Iranian 1 clade) remain to be better determined and this hypothesis should be further tested. In fact, as already highlighted by Jablonski et al. (2019), the uncertain phylogenetic position of some subclades as well as the low statistical support values at some nodes is probably due to a number of missing haplotypes and a more inclusive sampling with the addition of more molecular markers is probably required to better assess the evolutionary relationship of some subclades.
The karyological formula of the two females of C. austriaca here studied is consistent with that previously described by Mattey (1931) and Kobel (1967) for two males from Switzerland. However, the comparison among the pair 4 of the different karyotypes allowed us to detect for the first time the ZZ/ZW sex chromosome system in C. austriaca. In particular, the metacentric element of the 4 th pair was identified as the Z sex chromosome and its heteromorphic counterpart as the W chromosome. In addition, the W chromosome showed a different morphology among the two studied females, suggesting the occurrence of a chromosomal diversification among distinct molecular clades. Overall, the karyotype of C. austriaca shows a mixture of plesiomorphic and derivate chromosomal characters. In particular, the chromosomal characters that can be considered plesiomorphic in colubrids include the karyological formula, the microchromosomal localization of the NORs, and the pair 4 representing the ZZ/ZW sex chromosomes. A chromosome complement composed of 2n = 36 chromosomes with 16 macro-and 20 microchromosomes is displayed by both primitive (Henophidia) and advanced (Caenophidia) snake lineages and is supposed to represent the ancestral snake karyotype (Gorman and Gress, 1970;Singh, 1972;, 1986, , Oguiura et al., 2009Mezzasalma et al., 2014Mezzasalma et al., , 2016bMezzasalma et al., , 2019. Concerning the localization of NORs, the methods here used show their presence on two microchromosomes. Ag-NOR banding suggests that both loci of NORs are active. In fact, Ag binds to proteins essential to nucleolar structure and therefore to the transcriptional activity of ribosomal cistrons during the previous interphase (Howell, 1977;Jiménez et al., 1988). The occurrence of NORs on a microchromosome pair is not unusual among snakes, being exhibited in representatives of various families, including Colubridae (Olmo and Signorino, 2006), and it is considered a primitive condition in Squamata (Porter et al., 1991;Aprea et al., 2006). Derivate chromosome characters in the karyotype of C. austriaca can be considered the morphology of the pair 5 (telocentric in all the studied samples), and the occurrence of an heteromorphic, completely heterochromatic W sex chromosome, which is also morphologically differentiated betwenn the two studied samples from different geographic regions. In fact, in the putative ancestral snake karyotype of 2n = 36 all the macrochromosome pairs are biarmed (meta-or submetacen- tric), though a telocentric morphology of the pair 5 has been observed in different colubrid species belonging to independent evolutionary lineages, such as different species of the genera Elaphe and Hierophis (see also Singh, 1972;Kobel 1967;Mezzasalma et al., 2015b). According to Singh (1972), because the fifth pair is biarmed in most of the other colubrids, a simple pericentric inversions can be assumed to explain the different morphology of this pair. The ZZ/ZW sex determination system was supposed to be a plesiomorphic state in snakes (Mengden, 1981;Matsubara et al., 2016), however recent evidences suggest that different sex chromosome systems evolved multiple times, independently, in different snake lineages, including species with either female (ZZ/ZW) or male (XX/XY) heterogamety along with a discrete number of species with undifferentiated sex chromosomes (Gamble at al., 2017;Mezzasalma et al., 2019). This makes the suborder Serpentes, and more in general the whole order Squamata, which also includes various taxa with temperature-dependent sex determination (Gamble, 2010;Gamble et al., 2017;Pallotta et al. 2017;Alam et al., 2018), a unique study system to analyze the evolution and diversification of different mechanisms of sex determination. Nevertheless, in the family Colubridae, the fourth macrochromosome pair is usually composed of the ZW elements: the Z is metacentric and conserved in most species, while the W is often heteromorphic compared to the Z and largely heterochromatic (see e.g., Mengden, 1981;Mezzasalma et al., 2015a;Rovatsos et al., 2015;Matsubara et al., 2016). In C. austriaca, the studied Greek and Italian samples display a different morphology of the W chromosome, resulting submetacentric and subtelocentric, respectively. In the studied Italian female, the W chromosome is completely heterochromatic but not evidenced with fluorochromes (DAPI and CMA 3 ). The lack of data on the chromatin composition and distribution of the W chromosome from the Greek female here studied, does not allow us to establish if the differences between the W chromosomes of the Italian and the Greek sample also concern the heterochromatin pattern. Nevertheless, the different W morphology here found among Greek and Italian samples of C. austriaca highlight the occurrence of a karyological diversification among different clades of the species (Balkan and Western 2, see Fig. 1 and Jablonski et al. 2019), further supporting that the species contains different diverging evolutionary lineages. From a biogeographic point of view, as already documented for other Palearctic reptiles, the Quaternary climatic oscillations had an important role in shaping the current diversity of extant species, mainly through the contraction and re-expansion of their distribution ranges and the isola-tion of populations in different "refugia within refugia" (see, e.g., Ursenbacher et al., 2008;Gvoždík et al., 2010;Kindler et al., 2013;Mezzasalma et al., 2018). Unstable climatic conditions, fragmentation of the distribution range, small population size and isolation in distinct glacial refugia are also particularly favorable conditions for the fixation of chromosome mutations (see also Mezzasalma et al., 2015aMezzasalma et al., , 2017a for similar examples in different amphibian and reptile species), which may happen in different populations independently.
In conclusion, this paper provides the first record of a ZZ/ZW sex chromosome system in C. austriaca, with the occurrence of different morphologies of the W chromosome in different clades (Western and Balkan). More inclusive molecular and cytogenetic data from other areas of the wide distribution of C. austriaca would be useful to characterize the chromosome variability of different molecular clades the European smooth snake, helping to better assess their taxonomy.