A pool of data already existing about D-loop, i.e., the Control Region (CR) haplotypes of the mitochondrial DNA (mtDNA) of brown trout, Salmo trutta Linnaeus, 1758, tentative Adriatic trout Salmo farioides Karaman, 1938, and tentative Macedonian trout, Salmo macedonicus (Karaman, 1924), and their reconstructed phylogeography makes a good starting point for resolving their evolutionary history. That includes the dating of particular events in it. The events have hitherto been dated using the method of a molecular clock. Various calibrations were applied for the mutation rate, owing to the incongruence between the time of divergence that various authors notified and general knowledge about events in geological history and the periods in which they occurred in the Mediterranean region. Since geological history events were mandatory for setting the scene for the evolutionary history of brown trout, the incongruence between them has questioned the molecular clock calibration’s validity. From results about both the phylogeography and phylogenetic relations between native haplotypes (both partial and whole CR sequences) and the population genetics that characterized particular populations, we calculated the time of divergence between haplotypes in the regions of the western part of the Balkans: Iron Gate broader area in eastern Serbia, continental Montenegro and south-eastern Serbia. The distinct status of adjacent populations was verified by frequencies of microsatellites’ alleles and the STRUCTURE analysis that examined the significance of differences between them. In particular, we examined the populations that were clearly separated either by physical barriers, such as a waterfall in eastern Serbia (e.g., the upper and lower River Rečka supplemented by nearby rivers Vratna and Zamna), or by underground drops in Montenegro (e.g., upper and lower River Zeta, and rivers Nožica and lower River Mrtvica as isolated counterparts). We used the so far most common substitution rate of 1% in a million years’ (MY) period. The divergence times we obtained were compared to the events known for the region from available geological history data. There was a fairly good congruence between the dating obtained by the molecular clock method and that by geological history where the advanced, i.e., modern haplotypes, were concerned. In contrast, the congruence was worse for dating of divergence when more ancient haplotypes were in question, being much better if the mutational rate would be decreased to lower rates. That supported results both from the Rate Correlation Test about the independence of evolutionary rates in different lineages of brown trout, and from the Molecular Clock Test, which revealed that the evolutionary rate throughout the phylogenetic tree is not equal. That implies a difference in the speed of evolution in them, which was likely slower and faster, in the ancient, pre-Pleistocene haplotypes and the advanced, Pleistocene ones, respectively. The setting of the variable, or non-linear (i.e., logarithmic) speed of evolving seems helpful, since the early cladogenesis with the dominance of mutations was most likely combined afterwards with the acting of other evolutionary mechanisms, especially of genetic drift in populations that passed through the bottleneck episodes of the abrupt decrease in population size during the unfavourable periods of their evolutionary history.

brown trout, evolutionary history, geological history, tentative Adriatic trout, tentative Macedonian trout, molecular clock

Brown trout, Salmo trutta Linnaeus, 1758, is natively dispersed across a wide geographic area in the Northern Hemisphere and its overall variability is striking. Apart from the variability inherent in its geographic distribution, brown trout’s plasticity in life history traits in many local and often isolated populations resulted in particular ecological forms, i.e., the morphae that additionally complicated their taxonomy and nomenclature. Owing to that, by following the classification approach relying on the previously widely adopted typological species concept and employing predominantly external morphological features as taxonomic characters, more than 29 tentative taxa at both species and subspecies levels were assigned in time. Although new species concepts have been introduced since then, e.g., biological (Mayr 1969; Mayr and Ashlock 1991) and phylogenetic (Cracraft 1989), the proliferation of new locally dispersed brown trout taxa has not even slowed recently (Kottelat 1997; Kottelat and Freyhof 2007), increasing the current classification ambiguity.

Crête-Lafrenière et al. (2012) situated the origin of the genus Salmo at Oligocene, 26–29 million years (MY) ago and differentiation of the Atlantic salmon, Salmo salar Linnaeus, 1758, and the brown trout, Salmo cf. trutta Linnaeus, 1758, at mid-Miocene, 10–14 million years (MY) ago. Vladimirov (1948) and Osinov and Lebedev (2004) reported about the Salmo spp. fossil remnants occurring also in the Miocene layers (15 MY ago). There is another record of a fossil remain of Salmo immigratus^† (Kramberger, 1891) at the vicinity of Samobor (Croatia) that makes an ancestry of Salmo spp. traceable in the western part of the Balkans down to 13 million years (MY), in the mid-Miocene strata (Anđelković 1989). Artamonova et al. (2021), however, considered that fossil remnant identification inaccurate, owing to its poor preservation. The intensive cladogenesis of Salmo cf. trutta that Crête-Lafrenière et al. (2012) reported after the fossil remnants found in the layers at the Caucasus dated at 2.5–5.0 MY ago (Pliocene).

Brown trout have a wide geographical dispersal area, with many isolated populations. The explanation for their wide geographical distribution and intense biogeographical differentiation could be geological processes, primarily of orogenic nature, that led either to geographical separation of aquatic realms, or the birth of a new palaeogeographical area. For example, the ascent of the Alpine chain led to a partition of the Tethyan Ocean on two different biogeographical entities, the Mediterranean and the Paratethys Seas, yet around the Eocene/Oligocene boundary (Piller and Harzhauser 2005). These basins were further subjected to internal differentiation and changes in the paleoenvironment (sea-level and salinity fluctuations, connections, isolation of water bodies, etc.).

In contrast to evolutionary dating suggested by the fossil records of ancient brown trout, there is great variability in the dating by method of molecular clock that uses different average mutational rates. Using the molecular clock method relying on the average mutational rate of 1% at each 1 MY in the mitochondrial DNA (mtDNA), Bernatchez (2001) estimated that recent brown trout speciation started 0.5–2.0 MY ago and passed the differentiation during the Pleistocene glaciations that lasted about 0.7 MY. On sharing the same common base pairs at particular places in the structure of their mtDNA’s D-loop (i.e., the Control Region, CR) as synapomorphic characters, he defined initially five main clades that represent phylogeographic lineages - the haplogroups of brown trout, assigning the Danube (DA) haplogroup ancestral for all other haplogroups: Adriatic (AD), Mediterranean (ME), Atlantic (AT) and marmoratus (MA) (Bernatchez et al. 1992). Subsequently, two more haplogroups, the Duero (DU) (Suárez et al. 2001) and Tigris (TI) (Sušnik et al. 2005; Bardakçi et al. 2006), were introduced. The ancestral character of the DA clade of brown trout was concluded from their origin, estimated to 0.15–0.3 MY, and from their wide dispersal in southwestern Europe and adjacent regions. Using the Relaxed Bayesian molecular clock model on the nuclear loci, Pustovrh et al. (2014) estimated the emergence of the common brown trout clade’s ancestor to be around 1.4 MY ago (Pleistocene), the divergence of AT and DU haplogroups from DA haplogroup to 1.2 MY ago, while the divergence of AD and ME haplogroups from the DA haplogroup was estimated to 0.82 MY ago, i.e., to the mid-Pleistocene period. They reconstructed the sister-group relation between the AT-DU clade on the one and DA-AD-ME clade of brown trout on the other hand. The molecular diversity in each of the clades is ascertained by reconstruction of relations between CR sequences of brown trout in local populations defined as haplotypes on the autapomorphies featuring them. Using exclusively the complete CR haplotypes’ sequences in her analysis, Sanz (2018) inferred that the AD lineage is ancestral for all other brown trout lineages. She considered it reasonable to suggest that a large proportion of the ancient brown trout diversity, including the differentiation of the AD lineage, originated in the Balkans, with the alternative hypothesis of a pre-Pleistocene isolation of brown trout in the Asia Minor and diversification of DA, AD and TI lineages in its central part, Anatolia (Bardakçi et al. 2006).

Artamonova et al. (2021) shares the same opinion that many aquatic species originated in East or Central Asia and spread to Europe through the Paratethys, which was the migration corridor and a centre of origin of new aquatic species and genera since the Miocene. These processes were facilitated by the Late Oligocene desalination of the Paratethys (Popov et al. 2009). The Paratethys (Laskarev 1924), a vast intracontinental sea, began to form in the Late Eocene, about 37 MY ago, and was finally shaped during the Oligo-Miocene (Marović et al. 2002). Within it, distinct basins that experienced a complex pattern of changing seaways and land bridges were included, as well as those connecting them to the Mediterranean Sea and the western Indo–Pacific Ocean (e.g., Rögl 1998, 1999). The drainage area of the present River Rečka, that has recently evolved into the backcountry of the Iron Gate Gorge was within the Dacian basin of the Paratethys (Fig. 1), the smallest water body that roughly comprised the recent lowlands in southern Romania and northern Bulgaria, as well as northeastern Serbia (Jipa and Olariu 2009). The gradual infilling of the Dacian basin with the sediments derived by Carpathian Basin led to its disappearance and transformation into the accumulation area and a place of fluvial transport. That occurred during the Romanian era in Pliocene, 4.7–1.7 MY ago (Jipa 2015). According to Vasiliev et al. (2005), the transition of the Dacian Basin from brackish marine to continental domain was completed during a period about 700.000 years ago. It has been along with the whole former Paratethys subjected to vertical and a significant horizontal mobility and clockwise rotation, locally up to the Pleistocene. The established clockwise rotation for the Dacian basin was 5–20° (Marović et al. 2002).

Figure 1. 

The broader area of the River Rečka (square outline) within the former Dacian Basin (roughly corresponding to the Carpathian Foredeep Basin, Flisch Belt and area south of it), with the Inner Carpathian Orogens a source of sediments for filling the Dacian Basin.

Earlier reconstructions of phylogeny between brown trout haplotypes were accomplished on their partial CR sequences (Bernatchez et al. 1992; Bernatchez 2001; Bardakçi et al. 2006; Marić et al. 2006). Marić et al. (2006) and Simonović et al. (2017) worked out a reconstruction of phylogeny from partial CR sequences seizing to the poly-T block in brown trout from DA and AD haplogroups in streams at the peninsular divide between the Aegean and Black Seas in SE Serbia. The basal position and ancient, i.e., plesiomorphic character of particular novel haplotypes, e.g., Da*Vr (GenBank Accession Number #DQ318125), Da*Dž (#DQ318124) and Ad*Bož (#DQ318128) they reported was strongly supported. The Da-s6 (# MW589188) (Kanjuh et al. 2021) is of the same ancient character. In addition to molecular markers (Marić et al. 2006; Simonović et al. 2017; Škraba Jurlina et al. 2020; Kanjuh et al. 2020, 2021), the external morphology (Simonović et al. 2007) also indicated the importance of the Balkans in recent brown trout lineages’ delimitation. The presence of predominantly land-locked, resident form of brown trout, but also of the other trout (Salmo spp.) taxa originally occurring in the western part of the Balkans, is probably due to the warm climate and regression of Paratethys, as well as due to drying out of the Mediterranean in the upper Tertiary–Neogene, i.e., Late Miocene–Early Pliocene (i.e., the Messinian period), 5–6 MY ago. At the end of the Miocene, 5.5 MY ago, the Mediterranean Sea closed and almost dried (‘Messinian Salinity Crisis’) and then it was refilled with freshwater from the Paratethys Sea. This event may have facilitated the dispersal of freshwater organisms around the Mediterranean (Bianco 1990). Only a few tentative trout taxa, e.g., the lake-dwelling tentative Adriatic trout, Salmo farioides Karaman, 1938, in Lake Skadar system in Montenegro (Mrdak et al. 2006; Snoj et al. 2010; Škraba Jurlina et al. 2018, 2020), and marble trout, Salmo marmoratus Cuvier, 1829, in the Hutovo Blato marshland of the River Neretva system in Herzegovina (Snoj et al. 2010; Škraba Jurlina et al. 2020) recalled the potamodromy in their life history. Recent studies have shown that the present Lake Skadar is young and formed from initial marshland not earlier than 1200 years ago (Mazzini et al. 2015). However, the Skadar basin, with its system of sublacustrine springs, is ancient and dates back at least to the Pliocene, ca. 3 MY ago (Grabowski et al. 2018). Bošković et al. (2004) acknowledged the belief of the Italian botanist Baldacci (1929) that the permanent Lake Skadar has evolved from an ancient and very large sea bay into a marshland, due to significant supply of the alluvial sediments by the rivers Morača, Drim and Bojana, together with minor tributaries. The River Bojana cut its riverbed into the lake sediments (Bošković et al. 2004), whereas the initial rift structure of the majority of Balkan rivers, including the River Morača, which is the main inflow of Lake Skadar, as well as the rivers Struma (Strymon), Vardar (Axios) and Neretva, was formed along large faults in the mid–Late Miocene (e.g., Tzankov et al. 1996; Mulić et al. 2006).

During the Pleistocene glaciations, the stream-dwelling brown trout in the River Danube’s drainage area passed through a very dynamic evolutionary history. Recent secondary connections between the River Danube’s tributaries during the Pleistocene (Razpet et al. 2007) can partially explain the recent distribution of haplotypes in the Balkan Peninsula. Simonović et al. (2017) and Kanjuh et al. (2021) hypothesized that the occurrence of unique and vicariant brown trout haplotypes in particular river systems in the River Danube drainage area at the Western Balkans reflects different periods and areas of their emergence and diversification, when they had opportunities to spread, prior to the subsequent reduction to their recent, narrower dispersal areas. On one source, the Dacian basin of the Parathetys Sea, which would include the area of the recent eastern Serbia where the River Rečka is situated, communicated with the western, Pannonian basin over three straits until the Late Miocene (Badenian). Then the Wallachian uplift of a part of Carpathians caused the Iron Gate (Porto ferre) to retain as their only connection (Stevanović 1990). The other sources claim that the connection between the western, Pannonian and eastern, Dacian–Euxinic–Caspian basins of the Parathetys Sea diminished in the Late Miocene–Early Pliocene (i.e., Pontian–Messinian) (Anđelković 1970; Almaca 1990). Matoshko et al. (2019) estimated the final establishment of the Danube at 1.9–1.8 MY ago.

The most common native haplotypes in Serbia are of DA haplogroup, with AT and AD haplogroups detected as well (Marić et al. 2006; Tošić et al. 2016; Simonović et al. 2015, 2017). The most common DA haplotype in the brown trout populations throughout the River Danube drainage area in the western part of the Balkans is Da1, with recently even newly discovered sub-haplotypes (Kanjuh et al. 2021). However, the common haplotype from the River Crni Timok system and from the broader Iron Gate region is the modern and derived Da23c (#KC630984) (Tošić et al. 2014, 2016). It is indigenous for this area, since it occurs only in this region in the Iron Gate Gorge and the River Crni Timok system, and nowhere else. It is the most frequent haplotype there, recently occurring in 67% of brown trout. In addition to it, in the upper section of the northern fork of the River Rečka upstream of the Bledaria Waterfall, we detected only the otherwise widespread, but here vicariant, very limited and completely isolated Da1a haplotype (Tošić et al. 2016). In the western fork of the River Rečka, both Da23c and Da1 haplotypes were recorded (Marić et al. 2006). Tošić et al. (2016) and Simonović et al. (2017) reconstructed the Da1 haplotype as ancestral and the Da23c as the most derived in the modern, derived group of haplotypes of brown trout in the Western Balkans’ streams. The Da2 haplotype in streams of the River Crni Timok system that we recorded is considered non-native, i.e., stocked (Tošić et al. 2016).

Observing the novel results obtained in Montenegro, it is necessary to emphasize that the whole area of Montenegro itself displays complex geological and tectonic evolution and significant differences in the distribution and abundance of water resources ranging from arid karst areas to areas rich in both surface and groundwater. The domination of carbonate rocks in geological composition of the mountains in Montenegro enabled the development of karst process and karst relief (Bešić 1969) that have had a great influence on the subsequent glaciation leading to specific glaciers’ behaviour.

Water from the territory of Montenegro drains into two basins: the Adriatic Sea and the Black Sea, whose watersheds are separated by Triassic dolomites (Drobne et al. 1979).

Škraba Jurlina et al. (2018) reported an occurrence of the Adcs11 (#AY836340), Ad-M1 (#DQ381566) and Ad+Prz (#DQ318129) in the lower River Mrtvica, the right tributary of the River Morača in the Adriatic Sea basin of southern Montenegro. The occurrence of three native AD haplotypes there also testifies to the dynamic evolutionary history and multiple colonization events in that area. However, only the Adcs11 haplotype was detected in the resident, stream-dwelling tentative Adriatic trout from the upper River Mrtvica, supporting its ancestral character in relation to the other two AD haplotypes. The distinct character of two populations from the upper and lower River Mrtvica was explained by the partial migratory character of tentative Adriatic trout in the lower River Mrtvica, where the lake-dwelling tentative Adriatic trout ascend from Lake Skadar each year to the spawning grounds in the lower River Mrtvica (Škraba Jurlina et al. 2020). Further detailed research of brown trout in the lost streams at the Adriatic Sea basin, the River Nožica and in the upper River Zeta and in its tributary, the River Gračanica upstream (i.e., above) of their underground drops to the River Mala and lower River Zeta valley, respectively also revealed the occurrence of the Da1a haplotype (Simonović et al. 2017).

A broad karst massif, consisting of thick and permeable fluvioglacial sediments of the Würm glaciation stage, is situated between the River Morača Canyon and the lower River Zeta Plain, which considers tectonic depression formed at the Miocene/Pliocene boundary (Cvijić 1926). Groundwater from it enters either into the River Mrtvica, or into the lower River Zeta valley, as the surrounding Upper Cretaceous flysch sediments represent a physical barrier to water, particularly for streams coming from the west. Similar behaviour is displayed by the Permian–Triassic clastites in the River Gračanica valley of the upper River Zeta system. Cvijić (1926) also has pointed out the importance of faults in the catchment area of the River Zeta that allowed either the uplift of the northern block or subsidence of the southern block. He attributed tectonic movements of lower intensity within the downstream region of the River Morača to stabilizing of the Skadar depression as has been inferred from the up to 88 m reduced height between the neighbouring terraces and a thick fluvioglacial deposit.

This study aims to clarify and support the evolutionary history of brown trout populations belonging to either the different haplogroups (e.g., in the Adriatic Sea basin), or different haplotypes of the same haplogroup (e.g., in the broader Iron Gate region of the River Danube drainage area) in the Western Balkans, including the results from use of microsatellites’ loci as molecular markers. The records that we have gathered about the adjacent brown trout populations which remained entirely isolated, and phylogenetic relations between them that we reconstructed, illustrate the dispersal of brown trout of particular haplogroups and haplotypes as evolutionary units in various time periods over the western part of the Balkans. That is related to the facts from the geological history known in the region and their timing, in order to check the estimation of the timing of the events dated using the molecular clock technique. Calibration of the molecular clock using calibration points from externally derived dates, such as biogeographical and geological events, can be used for interpolation of divergence times (i.e., the event to be estimated falls within the calibration points or within the calibration point and the tip of the branch), their extrapolation (i.e., the event to be estimated falls beyond the calibration points), or combination of both (Wilke at al. 2009; Ho and Duchéne 2014). As for an application of the molecular clock in estimation of divergence time in salmonids, Osinov and Bernatchez (1996) stated the need for its re-calibration and proposed the substitution (i.e., the mutation) rate 0.8% per one million of year (MY^–1). Rambaut and Bronham (1998) already stated that molecular clock estimates of dating of evolutionary history events should be interpreted with caution owing to secondary contacts and isolation episodes between brown trout lineages. In addition, Apostolidis et al. (1997) stated that CR appears to evolve at a lower rate compared to other regions in Salmo spp. and used the substitution rate of 0.5%–0.9% MY^–1 of Martin and Pallumbi (1993), estimating the divergence time between main mtDNA lineages (AT, DA, ME, AD and MA) in the time period of 2.5–6.0 MY ago (Late Miocene, i.e., Messinian–Early Pliocene). Discrepancies that occurred in estimations of divergences in various papers and possibility of the CR’s lower rate of evolving urged Sanz (2018) to calibrate the molecular clock to its lower edge of the range commonly used for this fish group, at 0.75% MY^–1, as proposed by Shedlock et al. (1992). This resulted in estimations concordant with those obtained before, setting the linages’ differentiation at the Pliocene–Pleistocene boundary at 2 MY ago. The greatest modifying of the molecular clock 0.31% MY^–1 was proposed by Crête‐Lafrenière et al. (2012). If applied, it would set the events far back in the past, to the Miocene–Pliocene boundary. Hence, we aimed here to examine few actual cases of distinctness in a very close range in the western part of the Balkans by comparing the dating of divergence between them using the phylogeography, i.e., the evolutionary history data and molecular clock with the dating of events using the geological history data.

Since a majority of haplotypes in both areas of interest was inferred as modern, i.e., derived (Simonović et al. 2017; Škraba Jurlina et al. 2018), we added to this calibration analysis the more ancient ones from southeastern Serbia belonging to both DA and AD haplogroups (Marić et al. 2006; Simonović et al. 2017), in order to provide material for calibration in the earlier period of brown trout evolutionary history in the region.

The microsatellites’ analyses served to reveal the distinct population status of brown trout populations of the native haplogroups and haplotypes in streams in the Adriatic Sea basin in Montenegro and in the broader Iron Gate region, respectively.

In this study, three populations of brown trout from streams in Serbia and Montenegro with impassable barriers were analysed (Fig. 2). Data from River Zeta (upper and lower River Zeta) and River Nožica are novel results, while data about the genetic structure of brown trout from River Mrtvica and River Rečka (Marić et al. 2006; Tošić et al. 2014; Škraba Jurlina et al. 2018, 2020) are taken from previous studies. The upper River Zeta flows southward through the high karst plateau, the Nikšić Karst Field. At the rim of the plateau, it drops several hundred meters in elevation through crevices and caves, to emerge as the lower River Zeta. It merges with the River Morača from the Pleistocene.

Figure 2. 

Regions in Serbia and Montenegro with the broader Iron Gate Gorge area (A), Adriatic Sea basin of Montenegro (B) and southwestern Serbia (C), where the populations of brown and tentative Adriatic and Macedonian trout were sampled (1, River Rečka; green quadrate denotes the Bledaria waterfall; 2, River Vratna; 3, River Zamna; 4, River Mrtvica; 5, River Nožica; 6, upper River Zeta; 7, lower River Zeta, red rectangle marks the subterranean drop; 8, River Studenačka; 9, River Vlasina; 10, River Džepska; 11, River Vrla; 12, River Božica). The dashed line denotes the peninsular divide between Adriatic, Black and Aegean Seas.

The River Nožica is situated also at the high karst plateau in southern Montenegro, in the drainage area of the River Morača, the main tributary of Lake Skadar in Montenegro. At the western rim of the plateau, it drops through the subterranean crevice, to emerge as the short River Mala that joins the River Morača. The River Rečka from eastern Serbia is situated at the Mt. Miroč in the broader Iron Gate Gorge area. It joins the River Danube under the name River Reka. In the headwaters’ section, it is formed by the confluence of the two forks, northern and western ones. The northern fork in its most upstream part flows through the high and narrow mountain valley built by impermeable Cretaceous clastic rocks that enable a surface stream. At its end, it drops down from the (about) 12 meters high Bledaria Waterfall that is impassable for fish from the downstream section, where tufa is precipitating even today. The western fork is free of barriers upstream, all along to the karstic spring in a deep forest. The drainage area is in the southernmost part of the so-called Dževrin Greda, formed along a fault of the same name, which extends from the Tekija in the North–South direction about 18 km. Evolving as the vertical one, the fault has brought Cretaceous and the Upper Jurassic rocks over the younger, Pliocene units. However, such movements came to an end and the fault started to act until recently as a strike-slip or lateral right fault (Grubić 1992). Such activity resulted in migration of the eastern block to the south and the western block to the north, allowing the displacement of riverbeds. In the broader Iron Gate area, rivers Vratna and Zamna are closest to the River Rečka, making those populations adequate for this genetic analysis, since the microsatellite data for brown trout in the western fork of the River Rečka are lacking.

Sampling of materials for phylogeographic and genetic researches was accomplished during the period from 2004 to 2015. Brown trout anal fin clips (approximately 16 mm²) were collected by electrofishing using AquaTech device IG200⁄1 (input 12 V per maximum 15A DC, output 500 V, and frequency 65 P s^–1) in Serbia and Suzuki-Bosch (220V DC, Imax = 6A) in Montenegro and stored in 96% ethanol. All sampled fish from analysed streams were released alive immediately after the sampling.

Sequences of haplotypes from southeastern Serbia were already published in Marić et al. (2006) and Simonović et al. (2017), and those from eastern Serbia (the broader Iron Gate area) and continental Montenegro were published in Škraba Jurlina et al. (2018, 2020). In total, 144 samples were analysed for CR mtDNA haplotype, 18 samples from the former region, whereas from the latter two regions 126 samples were analysed for CR mtDNA haplotype using a partial sequence until poly-T block and 100 were analysed for eight microsatellites’ loci (Table 1). Total DNA was extracted from the tissue samples using the High-Salt Extraction technique of Miller et al. (1988) and Quick-gDNATM MiniPrep extraction kit following the manufacturer’s instructions (Zymo Research Corporation, Irvine, CA). The CR of mtDNA was amplified using the forward primer Trutta_mt_F (5′-TGAATGAACCTGCCCTAGTAGC-3′, designed by M. Brkušanin), and the reverse primer HN20 (Bernatchez and Danzmann 1993), under the following conditions: initial denaturation (95°C, 5 min) followed by 30 cycles of strand denaturation (94°C, 45 s), primer annealing (52°C, 45 s) and DNA extension (72°C, 2 min) in the programmable Applied Byosistems thermal cycler. Each PCR reaction in volume of 30 μL contained 10 μM of each primer (ThermoScientific), 10 mM dNTP, 10× PCR buffer with MgCl2 (Kapa Biosystems), 1 U of Taq polymerase (Kapa Biosystems) and 100 ng (i.e., 1 μL) of genomic DNA. Amplified DNA fragments were run on a 1.5% agarose gel using Applichem SYBR Green for visualization. PCR products were purified and sequenced at Macrogen Europe. Sequencing reactions were performed in a DNA Engine Tetrad 2 Peltier Thermal Cycler (BIO-RAD) using the ABI BigDye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems), following the protocols supplied by the manufacturer by single-pass sequencing on each template using the forward (Trutta_mt_F) primer.

Haplotypes’ sequences were aligned and edited using the Mega X (Kumar et al. 2018) software, and used for assessment of difference between them in positions, using either the whole CR (aligned to 993 bp in length), or its partial sequence up to the poly-T block (561 bp). The calculation of the time of divergence (t_D, in years) between particular haplotypes analysed here was based on the branching topology of phylogenetic tree (Wilke et al. 2009) inferred using the partial sequences of 561 bp of 13 haplotypes and applying the Maximum Likelihood (ML) method in the Mega X software package (Stecher et al. 2020). The branch length information, i.e., the time of divergence was estimated using the node depth (d_N, in number of substitutions per site), together with the substitution (i.e., the molecular clock) rate (λ, in number of substitutions per site and year):

$t_{\text{D}}=\frac{d_N}{\lambda }$

The calculation of divergence times between haplotypes followed the calibration of Bernatchez (2001) for mitochondrial DNA (mtDNA) at an average mutational rate λ₁ of 1% MY^–1 and that of Crête‐Lafrenière et al. (2012) at an average mutational rate λ₂ of 0.31% MY^–1. They served as an initial referent time of divergence for matching with the geological history events (Ho and Duchéne 2014) and for comparison with the time tree inferred using the RelTime method that does not require assumptions for linear rate variations to the supplied ML phylogenetic tree, with the Atcs1 haplotype as an outgroup, using the Mega XI software (Tamura et al. 2021). The calibration of nodes and distribution mode of the molecular clock together with the confidence interval (CI, a percentage of times expected to get close to the same estimate under the same experiment conditions) and offset (λ, distribution to the oldest fossil) were determined following Drummond et al. (2006). The testing of evolutionary rates in all brown trout lineages (i.e., haplogroups) was accomplished on 13 nucleotide sequences of haplotypes, using the Test of Molecular Clock. The test was performed comparing the ML value for the given tree topology with and without the molecular clock constraints under Kimura (1980) 2-parameter model (Kimura 1980). The independence of evolutionary rates in particular lineages was analysed using the Rate Correlation Test (Tao et al. 2019) for 13 haplotypes with the phylogeny given in phylogenetic tree reconstructed using the ML method and the best-fitting nucleotide substitution model with the Bayesian information criterion (BIC) that was determined using the jModelTest 2.1.7 software package (Darriba et al. 2012), with the 1000 bootstrap replications.

Eight microsatellite loci [SsoSL438 (Slettan et al. 1995), Ssa85 (O’Reilly et al. 1996), Str73INRA (Estoup et al. 1998, SSsp2216 (Paterson et al. 2004), Ssa410Uos (Cairney et al. 2000), SsaD190, SsaD71 (King et al. 2005), OMM1064 (Rexroad et al. 2002)] were combined in four duplex reactions with forward primer labelled with fluorescent dye. The final concentrations of PCR components were: 1× PCR buffer (Invitrogen), 0.2 mM dNTPs, appropriate concentration of primers, 1.5 mM MgCl₂, 0.5 U of Taq polymerase (Invitrogen), and ~100 ng of DNA. PCR programs were as follows: initial denaturation for 3 min on 94°C, followed by 30 cycles (35 for the 4^th duplex) of denaturation for 45 s on 94°C, annealing for 1 min on 60°C (57°C for the 4^th duplex) and elongation for 30 s on 72°C, and final elongation for 1 h on 72°C. Fragment analysis was performed using GeneScan 500 LIZ Size Standard (Applied Biosystems, USA) on an ABI-3130 Genetic Analyzer (Applied Biosystems, USA). Analysis was done using GeneMapper ID v3.2.1 software (Applied Biosystems, USA).

Alleles’ frequencies, number of alleles per locus, expected (H_exp), and observed (H_obs) heterozygosities, as well as F_st values representing a measure of differentiation between populations and gene flow (N_m) calculated (following Wright 1969) as

(1 – F_st)(4F_st)^–1

were obtained using GENETIX 4.05 (Belkhir et al. 2004), as well as factorial correspondence analysis (FCA). FSTAT 2.9.3.2 software (Goudet 2002) was used to calculate allelic richness expressed as the mean number of alleles per locus (Ā_n) and values for Fisher’s F statistics per locus (F_IS, F_IT, F_ST). Allelic richness (A_r) and private allelic richness (P_ar) expressed as a number of alleles in a population and as a number of unique, i.e., private alleles that are present only in one population, respectively were calculated using programme HP Rare (Kalinowski 2005). To detect recent bottlenecks in populations we used Wilcoxon sign-rank test and mode shift test in program BOTTLENECK 1.2.02 (Cornuet and Luikart 1996; Piry et al. 1999). Wilcoxon sign-rank test was used under stepwise mutation model (SMM) and the two-phase model (TPM) with 95% for proportion of single-step mutations and variance mutation size set to 12 (Luikart et al. 1998) for significance estimation with 10 000 iterations.

Population structure was analysed using the STRUCTURE 2.3.4 program (Pritchard et al. 2000), using the admixture model with proposed number of clusters, K = 10. Length of burn in periods was set to 50 000 with the number of Markov Chain Monte Carlo (MCMC) repeats of 7 for each K depending on convergence after burning set to 500 000. Structure Harvester software (Earl and von Holdt 2012) was applied to estimate most probable K according to Evanno et al. (2005).

Several mtDNA haplotypes were detected in analysed regions: seven of DA haplogroup, five of AD and one of AT haplogroup (Table 1). The upper Zeta carried haplotypes of DA haplogroup (93.33%) with two samples (6.67%) of AT haplogroup, while in the lower Zeta all samples carried haplotypes of AD haplogroup. All brown trout sampled from River Nožica had Da1 haplotype and the majority of samples (69.2%) from the lower River Mrtvica had Adcs11 haplotype, 23.1% had Ad*Prz and 7.7% had Ad-M haplotypes. Samples from the Northern fork of the River Rečka upstream of the impassable barrier (Bledaria Waterfall) possessed Da1a haplotype, and from the western fork haplotypes were Da1a and Da23c (Table 1).

The ML tree reconstructed (Fig. 3A) revealed the two main clades: that of the DA haplogroup and the second one comprising two sister subclades of the AT and AD haplogroups. The Test of Molecular Clock rejected the null hypothesis of equal evolutionary rate throughout the tree at a 5% significance level (with the molecular clock: +G = –35114.352, +I = 13; without the molecular clock: +G = –943.416, +I = 24; P < 0.05). The Rate Correlation Test for 13 haplotypes revealed independence (the score of 3.49 * 10^–7, P > 0.05) of evolutionary rates in particular brown trout lineages. The best-fitting nucleotide substitution model that jModelTest software found by lowest BIC values and Δ = 0 was Hasegawa–Kishino–Yano with Invariable sites (HKY + I). Calibration of nodes was accomplished using the exponential clock distribution (λ = 0.31, with the confidence interval CI = 1 – 12.9 and offset = 1).

Figure 3. 

(A) Phyogenetic tree reconstructed from 13 CR mtDNA haplotypes’ sequences of brown trout using maximum likelihood (ML) method from three (DA, AD and AT) phylogenetic lineages in the western part of the Balkans (numbers at each clade denote bootstrap probabilities); (B) Timetree inferred by applying the RelTime method to the phylogenetic tree shown in a) and the Hasegawa–Kishino–Yano substitution model, with the branch lengths below each of them and time of divergence to the right of each node.

Estimation of the time of divergence between brown trout that have the ancestral haplotypes using the partial CR haplotypes’ sequences and under the λ₁ = 1% and λ₂ = 0.31% (Table 2) revealed that a divergence between AD (Ad*Bož) and DA (Da*Vr) lineages occurred at 1.972 MY and 3.450 MY ago, respectively and that a differentiation in tentative Adriatic trout within the AD haplogroup (Ad*Bož and Ad-M1) occurred at 0.535 and 1.725 MY ago, respectively, with the most derived haplotype Ad-M1 who diverged from Adcs11 at 0.357 and 1.150 MY ago, respectively.

Calculation of divergence times by TimeTree of Mega XI software using the substitution rate λ = 0.31% (Fig. 3B) revealed that the time of divergence between AD (Ad*Bož) and DA (Da*Vr) lineages occurred at 4.740 MY ago and that a differentiation in tentative Adriatic trout within the AD haplogroup (Ad*Bož and Ad-M1) occurred at 1.698 MY ago, with the most derived haplotype Ad-M1 who diverged from Adcs11 at 0.378 MY ago.

As for the analysis of their microsatellite loci, in brown trout from the Adriatic Sea basin, the highest expected heterozygosity (H_exp) was detected in the lower Mrtvica (0.63) and the lowest one in the lower section of the River Rečka (0.20). The highest observed heterozygosity (H_obs) was in the lower River Zeta (0.64), as well, and the lowest (0.24) in the River Rečka. Generally, the majority of localities had slightly lower values for observed heterozygosity except the lower Zeta, Vratna and Rečka (Table 3)

Populations from upper and lower Zeta rivers had a greater number of alleles per locus than populations from rivers Nožica and Mrtvica (Table 3). Tentative Adriatic trout from the lower River Zeta had eight private alleles on SsaD71 and one on SsaD190. Single private alleles were identified at four loci for brown trout from the upper River Zeta (Str73INRA, Ssa410Uos, SSap2216, OM1064), three loci for tentative Adriatic trout from River Mrtvica (Ssa410Uos, SSap2216, OM1064), and one locus for brown trout from River Nožica (Ssa410Uos). The high number of private alleles indicates isolation of those populations that are confirmed by relatively small estimated gene flow (N_m) and relatively high F_st values between those populations (Table 4).

Allelic richness in the population of brown trout from northern fork of the upper River Rečka above the Bledaria Waterfall was the lowest in Iron Gate region (1.75), as well as the observed heterozygosity (H_obs = 0.24), but still higher than expected heterozygosity (H_exp = 0.20). The presence of genetic bottleneck was not detected in this stream, but the deviation from L-shape allele frequency distribution under mode-shift test was present. This test does not show statistical significance, but the shape of the curve indicates deviation from Hardy–Weinberg equilibrium (HWE) in this river. In the previous studies (Tošić et al. 2016), the population from the northern fork of the River Rečka above the Bledaria Waterfall was clearly separated from other Iron Gate populations by Factorial Correspondence (FCA), cluster and STRUCTURE analyses. Small number of alleles per locus (Table 5) and one private allele (locus Ssa410UOS) indicated long-term isolation, as well.

The STRUCTURE analysis included the novel tentative Adriatic trout samples from the upper and lower River Zeta and novel brown trout sample from the River Nožica, whose microsatellites were analysed here for the first time, and the tentative Adriatic trout samples from the lower River Mrtvica and brown trout samples from the rivers Zamna and Vratna in the Đerdap (Iron Gate) Gorge area we have already reported about (Škraba Jurlina et al. 2020), and re-analysed here together with the novel samples. Results from that comprehensive sample structuring revealed the highest Delta K value of 119.6168 for K = 5 (Supplement 1), i.e., for the five clearly distinct populations (from the rivers Nožica, upper Zeta, western fork of Rečka and two genetic clusters: lower Mrtvica with lower Zeta and Zamna with Vratna) (Fig. 4). These results clearly show grouping of brown trout populations together, as well as putative Adriatic trout.

Figure 4. 

The STRUCTURE analysis reveals the distinct population status (represented by distinct colour) of particular brown trout populations (NZ, River Nožica; UZ, upper River Zeta; REC, River Rečka with one genetic cluster; VRA + ZAM, rivers Vratna and Zamna; LM + LZ, lower River Mrtvica and lower River Zeta as one genetic cluster in putative Adriatic trout) . Y-axis presents the participation of alleles from all populations (expressed as a fraction of individual genome ranging from 0–1) in the genome of each individual in all populations, while bars in the x-axis represent individuals and colour-assigned obtained genetic clusters i.e., populations.

Comparison of dating that various authors gave for particular events in the evolutionary history of genus Salmo (Vladimirov 1948; Anđelković 1989; Osinov and Lebedev 2004; Crête‐Lafrenière et al. 2012) using both molecular and paleontological data has raised incongruence several times. Crête‐Lafrenière et al. (2012) estimated the age of the family to 59.1 MY with the Confidence Interval of 63.2–58.1 MY on 63 species they analysed and the average mitochondrial rate of molecular divergence to about 0.31% MY^–1 (with the Confidence Interval of 0.27%–0.36% MY^–1). That estimation and geological dating suggest that the values for times of divergence we calculated using the mutation rate of 1% MY^–1 (Table 2) should be enlarged about threefold, to seize back to about 6 MY ago for divergence between the most ancestral haplotypes of the AD and DA lineages in the area in concern. That value roughly corresponds to those reported using the geological dating, although one should always have in mind that the rate of evolution in early cladogenesis might not necessarily correspond to that occurring later. Analyzing both more plesiomorphic (e.g., Da*Vr, Da*Dž and Da-s6) and derived haplotypes (e.g., Da1, Da2, Da22, Da23c, etc.) in the DA haplogroup of brown trout from the River Danube basin of the Western Balkans, Simonović et al. (2017) found the uniform rate of evolution in their D-loop sequences. On the contrary, this analysis revealed both that the evolutionary rate throughout the phylogenetic tree is not equal and that evolutionary rates in particular haplogroups are independent. It seems that it is not possible to assess the difference in rate of mutation between more ancient and closer periods of brown trout diversification for certain using the mtDNA molecular clock approach, though our results testify to the slower evolutionary rate in ancient, pre-Pleistocene haplotypes and faster evolutionary rate in more derived, the ones from Pleistocene. That is in line with the findings of Crête-Lafrenière et al. (2012), who stated that the evolutionary rate is dependent upon the timing of divergence, with the faster rates more applicable to recent speciation events.

Another difficulty in calibrating the molecular clock lies in the weight, sensu Mayr (1969) and Mayr and Ashlock (1991), of the character used for that, either regarding the type of the molecule, part of the sequence of the same molecule type employed, or in the length of the sequence (i.e., a 100%) from which the time of divergence is to be estimated (Osinov and Bernatchez 1996). Values for divergence between haplotypes (Table 2) reveal remarkable difference between the ones assessed using the partial sequences (up to the poly-T block) and the whole D-loop sequences (Table 2). That difference originates from the obvious difference in variability between the two D-loop parts noticeable in comparing haplogroups AD and DA and haplotypes within each of haplogroups. While using only the part of D-loop, for example, divergence times between AD and DA haplotypes are greater than if the whole D-loop sequence is used. In contrast, that time is greater for two AD halotypes (Adcs11 and Ad-M1) when the whole D-loop sequence is used.

Earlier phylogeographic studies in eastern Serbia revealed the presence of few CR mtDNA haplotypes (Tošić et al. 2014, 2016). Nuclear markers used in this study examined matching of population genetics to earlier results. Moderate-to-high values of observed heterozygosity for trout fish were also reported by Mrdak et al. (2012). Allelic richness in almost all brown trout populations in eastern Serbia indicates high genetic variability, regardless of their small population size caused by occasionally harsh environmental conditions. The lowest genetic variability was detected in the population from headwaters of the River Rečka northern fork upstream of the Bledarija waterfall. This population was distinctive in many ways compared to other analysed populations. In addition to the exclusive presence of the private Da1 haplotype (Tošić et al. 2016), the brown trout population from that fork of the River Rečka featured a low number of microsatellite alleles (Table 3), who were all private (Škraba Jurlina et al. 2020), which determined its uniqueness in the group of populations from the Iron Gate Gorge region (Fig. 3). The occurrence of private haplotype and alleles, as well as the loss of genetic variability in that brown trout population (Table 3) testifies to its long-term isolation and its native character. Its haplotype was reconstructed as ancestral in relation to brown trout that hold more derived Da23c haplotype (Tošić et al. 2014; Simonović et al. 2017). A waterfall more than 12 m in height maintained unsurmountable isolation from brown trout in the downstream section. This obstacle has prevented introgression into the ancestral population situated upstream of the waterfall. The greatest genetic differentiation revealed by Watterson’s estimator (θ) and the lowest gene flow rate (N_m) values detected between distant brown trout populations from the River Rečka and other streams in the Iron Gate area (Table 4), as well as in River Crni Timok system (Škraba Jurlina et al. 2020) supported the very high differentiation of brown trout in the upstream section from the ones in the rivers Vratna and Zamna and most likely from those below the Bledaria waterfall. The STRUCTURE analysis for upstream situated brown trout revealed a very high level of differentiation both in this research and in Škraba Jurlina et al. (2020). The long-term isolation of the River Rečka population in the northern fork has likely favoured inbreeding and acting of the genetic drift as the most prominent evolutionary mechanisms acting there. Wilcoxon test revealed no significance for the bottleneck effect in this population (Table 3), but the deviation from the L-shaped distribution of allelic frequencies (the mode-shift test that has no statistical support) indicates its recent occurrence (Škraba Jurlina et al. 2020).

Indeed, the undertaken investigations on speleothem in the cave carved into the Barremian (Lower Cretaceous) reef limestones on the opposite, Romanian side of the Danube, showed that it was precipitated between ~75 ka and ~2 ka with at least two hiatuses (Constantin, unpublished^*). The obtained dating is the age when the stream bed has already been incised. Horvatinčić et al. (2003) found that speleothem growth started several thousand years earlier than tufa and that the majority of tufa deposits in Europe formed within the last ca. 7000 years. It should be noted that the dating of karst areas in Serbia based on correlation of heights of speleological objects with the alluvial terraces, should be taken with caution and that precise data, such as those obtained by radiocarbon dating are missing. For example, alluvial terrace of sandy clays strata next to the confluence of the River Reka (the lower section of the River Rečka) to the Danube are considered equivalents of Mindel (0.337–0.374 MY ago) and Mindel–Riss (0.3–0.337 MY ago) interglacial stages (i.e., Mid-Pleistocene) (Dimitrijević et al. 1997). That dating seems roughly comparable to the dating of divergence of the Da23c haplotype from the Da1a to at least 0.101–0.178 MY ago and at most 0.325–0.575 MY ago (Table 2), using the partial and whole CR sequences and two mutation rates, respectively.

The dating of tectonic formation of the lower River Zeta valley and River Morača to Miocene–Pliocene boundary (Cvijić 1926) provides the support for immigration of brown trout of the AD molecular lineage there. That is in agreement with the hypothesizing of the scenario that happened in the high karst plateau of the Adriatic slope of the divide much later, during the Pleistocene glaciations, when brown trout of the DA haplogroup colonized the River Nožica. During the Pleistocene period, marked by alternating largescale glaciations interfered by short warm intervals, the highest mountains in Montenegro were covered by glaciers. Most of them advanced towards the River Morača incising channels into bedrock. Nevertheless, the glacier at Širokar, situated in the present headwaters’ area of the River Tara, has been split by the so-called glacial bifurcation on two parts that have taken two different directions. The southern part has moved towards Lake Bukumir (Adriatic basin) and the northeastern part towards the River Vjeruša, a fork of the River Tara, i.e., into the Black Sea basin (Cvijić 1921; Đokić et al. 1976; Petrović 2007). Glacial bifurcation occurred as the consequence of karst relief instead of a fluvial one beneath a glacial cover, and was additionally supported by the gradual descent of the ice boundary allowing, albeit a minor part of the glacier, to reach Lake Bukumir, i.e., the Adriatic basin (Petrović 2007). Currently, this Pleistocene glacier bifurcating is the best available explanation for an occurrence (i.e., a transfer) of brown trout of the otherwise widespread Da1a haplotype in the River Nožica in the Adriatic Sea Basin, which was most likely initially devoid of any trout species. The native absence of any tentative brown trout taxon in particular streams of the Adriatic Sea Basin was already inferred, e.g., the River Vrijeka in Herzegovina and River Gacka in the high Lika Karst Field in Croatia (Simonović et al. 2017). The glaciers covering the highest mountains in Montenegro were exposed to melting during interglacial stages. Glaciations generally gave rise to erosional features in the highlands and depositional features on lowlands. Hence many of the large rivers in both seas’ basins in Montenegro (e.g., rivers Tara, Mrtvica, Morača) along with some of the greatest recent canyons (e.g., the Canyon Platije of the River Morača) owe their origin to meltwater rivers at the beginning of Holocene (Đurović and Petrović 2007). The disappearance of glaciers and the release of snow and ice weight have led to postglacial uplifting of the mountainous area in Montenegro (determined recent upthrows of 2–4 and more millimetres per year), whereas parts around Lake Skadar basin and along the Adriatic Sea coastline subsided (Ivanović 1991). A divergence of the most recent haplotypes of the AD haplogroup in the region, Ad-M1 and Adcs11, was estimated to 0.178–0.575 MY and 0.302–0.975 MY ago, respectively (Table 2) as well as to 0.378 MY ago (Fig. 3B). That fits well into the Pleistocene period. It suggests also that those haplotypes occurred before the melting of glaciers in Holocene in the former River Morača system occurring there. The physical barriers of glaciers and uplifting of mountains could explain an occurrence of only ancestral Adcs11 haplotype and the absence of the other more derived Ad-M1 haplotype in the upper River Mrtvica at high altitudes (Škraba Jurlina et al. 2018). However, the ancient occurrence of the whole recent Skadar basin with its system of sublacustrine springs that dates back at least to the Pliocene, ca. 3 MY ago (Grabowski et al. 2018) suggests the possible dating of appearance of the Ad-Prz haplotype of tentative Adriatic trout, the most plesiomorphic in the AD lineage there. That dating is congruent with the one of 1.725 MY and 0.403–1.299 MY we obtained using the molecular clock method at two rates and two CR sequence lengths, respectively (Table 2), as well as with the 2.052 MY using the TimeTree option (Fig. 3B), providing the scene for their occurrence there.

On the (1) high levels of diversity at the haplogroups’ level; (2) independent evolutionary histories and prominent endemism in tentative Adriatic and Macedonian, Salmo macedonicus (Karaman, 1924), trout at the Adriatic and Aegean Seas’ basins, respectively reported from results of molecular studies, as well as (3) ancestral character reconstructed in particular Balkans populations (Marić et al. 2006; Apostolidis et al. 2011), this region was considered important for evolving of the AD lineage. The other region that Bardakci et al. (2006) and Arslan and Bardakci (2010) suggested is the Asia Minor, positioning the start of diversification there before the Pleistocene (i.e, in the Pliocene, more than 2 MY ago), while geological history events’ dating reached the Miocene/Pliocene boundary (Cvijić 1926). Apostolidis et al. (1997) stated that the high level of endemism in AD lineage could be a consequence of ancient isolation and independent evolving of allopatric populations with limited natural contact between them. For example, the vicariant Adcs1 haplotype that Sušnik et al. (2007) have considered ancestral for the AD haplogroup was questioned by finding of the Ad*Bož haplotype of the very limited dispersal in tentative Macedonian trout living in sympatry with their conspecifics with the Adcs1 haplotype in few headwaters at the slopes of the Vlasina Plateau of the SE Serbia, which join the River Struma (Strymon) (Marić et al. 2006).

The contours of the whole Balkan area were finally shaped by geodynamic events during the Neoalpine history and the interaction of three geodynamic processes: Pannonian collapse, continuous epeirogeny upward movement and the Aegean collapse (Marović et al. 2007). The former, which started in the Ottnangian–Karpatian (18.1–17.2 MY ago), less in Badenian and Sarmatian led to development of subsidential structures, such as the River Morava Graben (e.g., Marović 1990; Marović et al. 2002). The strong epeirogenic uplifting of 500 to 1500 m proper to southwest and west attained the epeirogeny uplifting over 1500 meters in the Vlasina Plateau. Opposite it, the subsidence declined successively from south to north through the Pannonian into the Pontian, but might be renewed in some areas during the Pliocene/Quaternary and even the Upper Quaternary, as the consequence of the Aegean orogen collapse.

The course of the River Struma/Strymon is mainly controlled by the NNW–SSE trending Struma Lineament, a large tectonic structure that represents the tectonic boundary of the Serbo–Macedonian and the Rhodope massifs (Zagorchev 2007). The onset of the River Struma basin is attributed to collisional processes that have taken place from Late Oligocene to Middle–Late Miocene between the Apulia and Eurasia plates (Tranos et al. 2007; Tranos 2011). At the last glacial maximum (21 500 years BP), when the sea level was 120 m lower than today, the rivers Mesta/Nestos and Struma/Strymon formed a joint delta and found their way to the sea, as well as the Vardar/Axios (Perissoratis and Conispoliatis 2003).

In the inferred phylogeny, the ancient Ad*Bož haplotype was a sister clade of all others in the AD haplogroup (Marić et al. 2006). There was only one difference between the Adcs1 and Ad*Bož, and there were five and two different base pairs between each of them, respectively and Ad*Prz that is considered closely related, and six base pairs between each of them, respectively and the Da*Vr haplotype (Table 2), the most basal in the DA haplogroup (Fig. 3A). Cortey et al. (2004) estimated the differentiation of southern trout haplogroups using the molecular clock method at the Pliocene–Pleistocene climate cooling that occurred at 1 MY ago, with the strongest diversification in the AD haplogroup occurring at 0.2–0.15 MY ago, which is much sooner than our estimations of 2.139–6.900 MY ago between the Ad-M1 and Da23c (Table 2) and 5.730 MY ago (Fig. 3B), that is at the Pleistocene to the late Pliocene) for differentiation of southern brown trout haplogroups, and of 0.357–1.150 MY ago and 0.302–0.975 MY ago (Table 2), as well as of 0.378–2.052 MY ago for diversification in the AD haplogroup (Fig. 3B) in the Pliocene–Pleistocene boundary and in Pleistocene made using molecular clock at the mutation rates of 1% MY^–1 and 0.31% MY^–1 and using partial and whole CR sequences and TimeTree inferring, respectively. The dating of geological history events gives the scene for, i.e., supporting the ancient and ancestral character of both Da*Vr and Da*Dž on one, and Ad*Bož on the other side of the peninsular divide in the DA and AD haplogroups in the Vlasina Plateau area, respectively.

In conclusion, it seems that in addition to places of its occurrence, the dating of recent Salmo spp. evolutionary history remains vague. In addition to the phylogeography of brown trout, the reconstructed phylogenetic relations determining the ancestry and descendancy of particular evolutionary units (here referring to haplotypes) should be accounted in considering the dynamics of their evolving. The matters we again have posed here, e.g., (1) the speed of evolving in more recent, i.e., advanced and older, i.e., ancestral periods and evolutionary units (here referring to haplotypes), and (2) levels and the variability of parts of the D-loop, their usefulness in calibration of the molecular clock and their validation with the dating known from geological history, still need more research to get us closer to a more reliable estimation of dating of the divergence between brown trout lineages throughout their dispersal areas. While some dating of the more recent events of evolutionary history seems to correspond well to those of geological history and at local scale, others need more consideration in reaching that aim. It is certain that the ancient scene for evolving of ancestral brown trout haplotypes has been created from geological events either at a wider, continental geographical scale, e.g., gradual succession of the Parathetis Sea, or at the more regional scale, e.g., the formation of the river system in the Skadar Plane and during the Rhodopi Mountains uplifting at the Miocene–Pliocene boundary and mid-to-Late Miocene, respectively. The evolution and spread of the more recent haplotypes, on the other hand, has occurred on that scene in the circumstances of the Late Pliocene and Pleistocene cooling and glaciations, when glacial interconnections (e.g., at the Nikšić high karst plateau and that of the Širokar glacier bifurcating at the watershed divide between the River Vjeruša, i.e., the River Tara’s fork, and River Nožica) played a significant role in dispersal of more advanced haplotypes. It might be hypothesized that those two “acts” of brown trout evolution lasted differently, i.e., unrolled at a different speed and that the distinct molecular clock (i.e., substitution-mutational) rates should be established for them, or the polled, i.e., the logarithmic rate instead of the current uniform, linear one should be set. That period of advanced brown trout and related tentative trout taxa evolving during the Pleistocene glaciations most likely included the acting of evolutionary mechanisms other than mutations (e.g., genetic drift), when they likely used to decrease in population size and pass through the bottleneck episodes. However, the evolving of their CR, their phylogeography and phylogenetic relations give us an opportunity to trace their evolutionary history in the scene that geological history has set to them.

Work was granted by the Ministry of Education, Science and Technological Development of the Republic of Serbia (Contracts No. 451-03-68/2022-14/200178, No. 451-03-68/2022-14/200007 and No. 451-03-68/2022-14/200126 for the University of Belgrade, Faculty of Biology, Institute for Biological Research “Siniša Stanković” and Faculty of the Mining and Geology, respectively).