Natural toxic impact and thyroid signaling interplay orchestrates riverine adaptive divergence of Dolly Varden


 Background Adaptive radiation in fishes has been actively investigated over the last decades. Along with numerous well-studied cases of lotic radiation, some examples of riverine sympatric divergence have been recently discovered. In contrast to the lakes, the riverine conditions do not provide evident stability in the ecological gradients. Consequently, external factors triggering the radiation, as well as developmental mechanisms underpinning it, remain unclear. Herein, we present the comprehensive study of external and internal drivers of the riverine adaptive divergence of the salmonid fish Salvelinus malma. In the Kamchatka River, this species splits in the reproductively isolated morphs that drastically differ in ecology and morphology: the benthivorous Dolly Varden (DV) and the piscivorous Stone charr (SC). To understand why and how these morphs originated, we performed a series of field and experimental work, including common-garden rearing, comparative ontogenetic, physiological and endocrinological analyses, hormonal "engineering" of phenotypes.Results We revealed that the type of spawning ground acts as the main external factor driving the radiation of S. malma . In contrast to DV spawning in the leaf krummholz zone, SC reproduces in the zone of coniferous forest, which litter has a toxic impact on developing fishes. SC enhances resistance to the toxicants via metabolism acceleration provided by the elevated thyroid hormone content. These physiological changes lead to the multiple heterochronies resulting in a specific morphology and SC's expansion into a piscivorous niche.Conclusions S. malma represents a notable example of how the thyroid axis contributes to the generation of diverse phenotypic outcomes underlying the riverine sympatric divergence. Our findings, along with the paleoecology data concerning spruce forest distribution during the Pleistocene, provide an opportunity to reconstruct a scenario of S. malma divergence. Taken together, obtained results with the data of the role of thyroid hormones in the ontogeny and diversification of fishes contribute a resource to consider the thyroid axis as a prime director orchestrating the phenotypic plasticity promoting evolutionary diversification under the changing environmental conditions.

Our findings, along with the paleoecology data concerning spruce forest distribution during the Pleistocene, provide an opportunity to reconstruct a scenario of S. malma divergence. Taken together, obtained results with the data of the role of thyroid hormones in the ontogeny and diversification of fishes contribute a resource to consider the thyroid axis as a prime director orchestrating the phenotypic plasticity promoting evolutionary diversification under the changing environmental conditions. Background Sympatric speciation commonly arising from competing for resources is widely distributed throughout various fish taxa (Robinson & Wilson, 1994;Bolnick & Fitzpatrick, 2007;Via, 2009;Nosil, 2012  In order to ascertain external and internal factors driving riverine adaptive radiation, we studied the 'flock' of charrs (Salmonidae) dwelling in the Kamchatka River. This basin is one of the most ancient and the largest water net of the North Pacific, which was not glaciated during the late Pleistocene (Chereshnev, 1998; Elias & Brigham-Grette, 2013; Barr & Solomina, 2015). In the absence of glaciers, the Kamchatka River acted as a refugium providing for a diversity burst in the lineage of the charr Salvelinus malma (Brunner et al., 2001;Oleinik et al., 2015). Along with the partially anadromous Dolly Varden (DV) there are several scarcely investigated sympatric non-anadromous morphs: i) the benthivorous morph inhabiting the upper tributaries; ii) the light-colored piscivorous morph in the lower course; and iii) the marble-colored piscivorous morph in the middle course of the river (Savvaitova, 1970 andSavvaitova & Maksimov, 1970;Pavlov & Savvaitova, 1991;Bugaev et al., 2007;Salmenkova et al., 2009).
The morphological and ecological peculiarities of the latter, the so-called 'stone charr' (SC), were described in detail (Savvaitova & Maksimov, 1970;Pavlov & Savvaitova, 1991; Melnik et al., 2020). As it was established, SC is defined by the traits typical of a piscivorous fish: a big head, elongated jaws, the presence of basibranchial teeth and barbs on gill rakers, caudally shifted fins, and reduced pyloric caeca number. It also has a specific pigment pattern and demonstrates an accelerated growth rate as contrasted to DV. The studies of microsatellite polymorphism, mtDNA and nuclear DNA marker variability demonstrated that SC is an isolated morph, which diverged from DV in the late Pleistocene -the early Holocene (Melnik et al., 2020).
The spawning sites of DV and SC are spatially segregated. The former migrates to spawn upstream to the leaf krummholz zone, while the latter breeds in the lower sections located in the coniferous forest. These sites do not display substantial discrepancies in the annual temperature regimes and water velocity, but drastically differ in the composition of bottom substrate and water chemistry. DV spawns on the pure pebble substrate, whereas SC reproduces on the substrate covered with a layer of conifer needles (Melnik et al., 2020). The conifer litter displays lower decomposition rate than the leaf one (Graca, 2001;Laitung et al., 2002;Ormerod et al., 2004). It also deoxygenates, acidifies the water, and continually supplies the environment by potentially toxic water-soluble and suspended compounds (Tremolieres, 1988 Multiple phenotypic novelties arising in rather a short evolutionary period coupled with a low level of genetic changes indicate the impact of an internal factor influencing numerous critical decisions over the course of SC life history. Among others, thyroid hormone (TH) level seems to be the most likely candidate due to their pleiotropic effect on the development and phenotype of fishes. TH influence the metabolic rate, the rate and timing of various ontogenetic events including skeletal and pigment patterning, affect the morphology, behavior and reproduction (Janz, 2000 Following them, we propose thyroid signaling as an internal factor driving the occurrence of SC phenotypic novelties in response to the toxic environment produced by coniferous litter. To test this hypothesis we performed a series of field and laboratory works: (1) assessed the resistance of both morphs to the conifer litter decomposition products ; (2) compared the metabolic and developmental rates, somatic growth, temporal characteristics of skull ossification, and ontogenetic dynamics of TH level in DV and SC; and (3) examined the contribution of TH in the occurrence of phenotypic novelties and resistance to the conifer litter decomposition products.

Normal development
The diameter of SC and DV eggs was similar: 0.58 ± 0.13 cm and 0.50 ± 0.15 cm, respectively (t-test df = 98 p = 0.889). In SC group, the hatching occurred at 335370 (50% achievement ~ 350) day*degree post fertilization, whereas in DV group it occurred later, at 360400 (~ 380) day*degree post fertilization. The size and weight characteristics of the newly hatched SC and DV were commensurable, FL = 1.61.9 cm and W = 6580 mg. The transition to larva stage proceeded in SC at 105165 (50% ~150) day*degree post hatching (= ddph), in DV at 130190 (~ 170) ddph. The onset of the fry stage occurred in SC at 450530 (~ 510) ddph, and 490570 (~ 540) ddph in DV (Fig. 1, upper inset). For the description of the postnatal developmental stages see Methods section (also Supplementary Fig. S1). The significant differences in FL and W between the experimental SC and DV were found at 500 ddph ( Fig. 1a and b). These discrepancies increased 215 dd later at the end of the experiment, and FL / W of SC fry were larger about 15% / 40%, respectively, than those characteristics of DV fry. In nature the differences in FL and W were found as well. The wild SC parrs (~ 850 ddph) were significantly larger than the same-aged DV parrs (Table   S2).
The total whole-body triiodothyronine (T 3 ) content was similar in the newly hatched SC and DV. Then, prior to the onset of external feeding, both groups displayed a synchronous drop of the hormonal level. In SC, the decrease was more pronounced (H-W test H 5;131 = 21.27 p = 0.0160 for SC and H 5;186 = 18.32 p = 0.0242 for DV). As the result, at 160 ddph SC and DV prelarvae significantly differed in T 3 content (Fig. 2). In further development, until the transition to the fry stage, SC and DV demonstrated similar positive dynamics of T 3 . At the fry stage, the hormonal level within the groups did not display any serious fluctuations, but differences between the groups increased. At the end of the experiment, SC fry had significantly higher level of T 3 . In nature SC parrs demonstrated the significantly higher values of T 3 content relative to DV parrs. The T 3 content values in both SC and DV wild parrs were 1.31.4 times higher than those detected for the experimental fry (U test df = 39 p = 0.0052 for SC and df = 59 p = 0.0094 for DV) (Fig. 2).
The ontogenetic dynamics of the whole-body thyroxin (T 4 ) content was similar in the experimental SC and DV. Both morphs displayed a gradual increase of T 4 level from the hatching to the fry stage (H 3;17 = 6.58 p = 0.0310 for SC and H 4;20 = 8.11 p = 0.0397 for DV). Nevertheless, we found a significantly lower hormone level in SC larvae and fry compared with the same-aged DV (Fig. S2).
We did not reveal any significant fluctuations of the biochemical parameters within the groups of experimental fry (H tests p > 0.05 for all parameters). However, the stable differences between the artificially reared SC and DV fry, as well as between the wild SC and DV parrs were noted (Fig. 3). In the wild and experimental DV, the blood glucose level was significantly higher than in the wild and experimental SC. The phospholipid content was significantly higher in the wild DV than in the wild SC. The same but insignificant differences in the phospholipid content were detected for the experimental DV and SC. In contrast, the tissue antioxidant activity in the experimental DV was significantly lower than in SC.
We failed to find distinguishable discrepancies in the head and body shape between the experimental SC and DV at the early life stages: in prelarvae Procrustes ANOVA The plot of principal component scores depicted the ontogenetic channels of SC and DV, and the timing once morphological differences appeared (Fig. 4a). At the early stages (prelarva, larva and late larva) DV and SC had partially overlapping channels. During the transition to the fry stage, the channels of both morphs changed the direction and, at the fry stage (> 650 ddph), became separate. According to the landmark loadings on PC2 (Table S1), the most expressed differences between SC and DV fry occurred in the eye diameter, maxilla length, head height, and ventral fins position. These morphological discrepancies were reflected in the 'consensus' body shapes stretching along the canonical root ( Fig. 4b and Table S1), which clearly demonstrated the head enlargement and the caudal shift of the fins in SC fry relative to DV fry. Similar morphometric differences were revealed for the wild parrs of SC and DV (Fig. 4c). The additional data concerning the main differences in the head proportions and fin position in experimental fry and wild parrs are present in Table S2.
SC and DV differed in the ossification rate of skull bones. The MANCOVA showed that, despite the strong effect of age (p ≤ 0.0321), a 'group' predictor was significant in determining the ossification rate of supraethmoid, frontal bone, maxilla, vomer and preopercle (Table). Two 'modules' of cranial ossifications -neurocranium and teeth-armed bones -displayed significant differences in the developmental rates between SC and DV.
Thus, the initial advancement in ossification of the neurocranium in SC prelarvae was replaced by the retardation during the subsequent developmental stages (Fig. 5a). As the result, DV demonstrated a significantly more ossified neurocranium starting from the larva stage. In contrast, the teeth-armed bones displayed an accelerated ossification in SC from the early stages (Fig. 5b). The ossification rate of preopercle resembled the ossification pattern of the neurocranium bones being advanced in DV. The intraoral tongue-bite apparatus did not demonstrate a significant difference in ossification rate. We discovered a minor advancement in the ossification of these bones in the SC late fry (Fig. 5c)  The goitrogen treatment and administration of T 3 during the larva-fry transition (from 160 to 515 ddph) changed the thyroid status of experimental fish ( Fig. 2 and S2, in boxes). In In both SC and DV reared under the 0.2 g l − 1 and 0.5 g l − 1 solutions of thiourea, the mortality rate was 10%, and 20%, respectively. The mortality rate in DV reared under T 3 treatment was higher: 25% in the 1.0 µg l − 1 solution, and 55% in the 5.0 µg l − 1 solution.
In the control groups of DV and SC, the mortality rate was drastically lower, 1% of individuals only.
The goitrogen-caused decrease in the biochemical parameters was statistically insignificant in both groups (K-W test p > 0.50), but, in fact, more pronounced in SC than in DV: H 2;58 = 2.85 / H 2;81 = 0.37 for the blood glucose, = 2.62 / 0.33 for the tissue antioxidant activity, = 0.19 / 0.05 for the phospholipid content, respectively. T 3 administration provoked two times stronger relative to the thiourea treating but still insignificant physiological response in DV (Fig. 3, in boxes).
The goitrogen negatively affected the fish growth rate ( Fig. 1c and d). The significant differences in size were observed between normal and thiourea-treated SC (Table S3; ANOVA for length F 2;110 = 19.04 p ≤ 0.0250). In DV, the reaction to the thiourea treatment was less pronounced; the significant size differences were revealed between the control fish and the fish reared under the 0.5 g l − 1 solution only (F 2;143 = 4.86 p ≤ 0.0011). In contrast to thyroxin, T 3 administration did not affect growth rate ( Fig. 1c and d; Table S3).
The alterations of thyroid status resulted in the developmental abnormalities. In the majority of hypothyroid SC late larvae and early fry (90%), the body walls did not close on the belly around the remains of the yolk sac. However, we failed to find similar defects in The supraethmoid ossification in the hyperthyroid DV did not exceed that in the control group.
The geometric morphometrics revealed that the hypo-and hyperthyroid DV displayed different vectors of morphological changes and were aligned along the first canonical root in CV morphospace in the following order: 0.5 g l − 1 of thiourea -0.2 g l − 1 of thioureacontrol -1.0 µg l − 1 of T 3 -5.0 µg l − 1 of T 3 (scoliotic individuals were excluded from the analysis). The vector of morphological changes in the thiourea-treated SC coincided with that in the thiourea-treated DV (Fig. 6a). The appropriate loading on landmarks (Table S1) and stretching of the 'consensus' body shape along the first root (Fig. 6b) Table S4.
SC demonstrated additional morphometric shift along the second canonical root (Fig. 6a).
The appropriate loading on the landmarks (Table S1) and the stretching of the 'consensus' body shape along the second root (Fig. 6c)  The goitrogen and T 3 treatments altered DV response to the impact of 1.5% coniferous infusion. In the euthyroid fry, the mean mortality rate was 23% in the infusion and 0% in the pure water. In the survived euthyroid fish, an increase of T 3 content (+ 1.91 ng g − 1 ) was detected, but intra-sample variance remained stable. In the hyperthyroid fry, the mean mortality rate was 10% only. The survived individuals demonstrated a two-times increase of T 3 content accompanied by a splash in the intra-sample variance. In the hypothyroid fry, the mean mortality rate was 54%. They also manifested elevated T 3 content (+ 0.34 ng g − 1 ) and 1.5-fold drop in the variance. The alterations of T 4 content were also detected (Fig. S2).

Discussion
The obtained results support the proposed scenario where spawning grounds' contamination and thyroid signaling act as drivers of the Kamchatka River S. malma diversification. We have managed to reveal the toxic influence of coniferous litter on the developing charrs and the difference in the resistance to this impact between DV and SC.
The former displays a conspicuous physiological response resulting in a significant increase of the mortality. The SC reaction is less pronounced. This finding seems to address the question of different spawning sites of DV and SC. We suggest that toxic environment in the coniferous forest zone progressively built the major reproductive  . Therefore, intensification of the thyroid signaling activity seems an essential adaptation at sites within the coniferous forest zone.
The accelerated thyroid activity observed in nature and in the experiment inherent to SC in comparison with DV supported well the last proposition. Moreover, the experimental data (rearing SC and DV in the pure water) clearly indicate that T 3 content difference between morphs is not a consequence of phenotypic plasticity in response to toxic impact only, but genetically determined hormonal phenotypes: ancestral DV and hyperthyroid SC.
As it was mentioned above, TH regulate charrs' metabolic rate. The comparison of the biochemical parameters revealed that SC is characterized by a higher metabolic rate, both under natural and experimental conditions. Thus, unlike DV, SC has a physiological adaptation to the water containing excessive amount of coniferous decomposition products.
However The evident morphological differences between the morphs emerge during the larva-fry transition, i.e. simultaneously with the hormonal discrepancies. This finding indicates that TH initiates the morphological divergence between SC and DV. It also backs up the hypothesis that the main mechanism driving sympatric microevolution in fish is an early ontogeny tuning promoting the phenotype divergence (Carroll, 2008 Table S5).
The fertilized eggs of SC and DV were obtained via the dry method. The crossing of each group was obtained by mixing the sexual products collected from three females and three males to obtain no less than 1000 viable larvae. The rate of fertilization was close to 90% in both SC and DV. The spawners from both groups were selected to be the same size FL = 4146 cm. In two days (48 h), the eggs of each group were delivered to the laboratory.

Temporal characteristics of development
To compare the temporal developmental characteristics of the SC and DV ontogeny we used a traditional day*degree scale (= accumulated daily temperature). For the convenience, we divided the ontogeny in two periods: pre hatching and post hatching.
Therefore, we used day*degree post fertilization to compare the timing of hatching and day*degree post hatching (ddph) to compare the timing of further development. Six postnatal developmental stages typical of salmonids, i.e. free embryo, prelarva, larva, late larva, fry and parr were identified following Balon (1985) (see Fig. S1).
The experimental rearing lasting for 36 weeks (715 ddph) was finished at the stage of the late fry. In accordance with the data of temperature loggers installed for a year at the SC spawning site (see Melnik et al. (2020)), the age of the wild parrs was about 850 ddph.
The annual dynamic of temperatures in nature is presented in Fig. S4.

Experimental design
The fertilized eggs were incubated in 250-l tanks filled with UV-treated soft water To study the normal development, the role of TH in ontogeny, and acute oxidative effects of coniferous substances we used three experimental sets. First, a common-garden experiment focused on the development of control groups of SC and DV from the free embryo to late fry stage. The fish were reared in cages in a single tank and sampled for the analysis six times (Table S6). As a criterion of sampling timing, we used the transition of 50% of fish from one stage to another. The specimens were collected randomly. The experimental fry were compared with the parrs collected from the nature.
The second experimental set was aimed to assess the TH involvement in SC and DV development. For this purpose, we reared the fish since 50% switching to external feeding for ten weeks (350 dd) under the hypo-and hyperthyroid conditions. In the first case we Data processing

Morphometrics
The left side of every individual was photographed in the orthogonal projection with the camera DTX90 (Levenhuk). To avoid any optical distortions the image width was always at least twice the fish length, and the same zoom was applied to all images. Fork length (FL) was measured using tpsUtil and tpsDig2 v2.16 software (Rohlf, 2010). The total dry weight (W) was assessed with an accuracy of ± 0.001 g. The size-weight characteristics were compared using a t-test and one-way ANOVA (Tukey test) in Statsoft v.10.0; the significance was considered at the p < 0.05 level.
Thirteen landmarks obtained from digital images were used for the description of the head and body shape in the experimental larvae, fry and wild parrs (Fig. S5). The geometric morphometrics technique was applied for the analysis in MorphoJ v1.06d. We used the In this version of the analysis, the first component reflects size variation, while the second one describes the main part of body shape variability, so the space of the first two components displays the "channels" representing allometric relationships among fish of different size (Tissot, 1988;Mina, 2001). calculations were based on the molar extinction coefficient ε = 2220 l mol − 1 cm − 1 (Goth, 1991

Consent for publication
The manuscript does not contain any individual person's data in any form

Availability of data and materials
All data generated or analyzed during this study are included in this published article [and its supplementary information files]

Competing interests
We declare the absence of competing interests

Funding
The transportation costs, as well as expenses for essential expendables, aquaria chemistry and equipment, chemicals for ELISA and biochemical analysis, food for larvae and juveniles spent during the field and laboratory works, were covered by the Russian Additional file 3. Figure S3. -preopercle (0, 1 -fused scute, 2 -superior plate with seismosensory apertures, 3formed walls of seismosensory channel, 4 -channel is closed, 5 -inferior process protrude over the opercle edge, crests are formed).
Position of the bones is shown additionally on the skull Additional file 7. Tables S1-S6. Table S1. Loadings of landmarks on the principal component axes and canonical roots in the analysis of fish body shape variability. Table S2. Linear morphometric measurements for experimental fry and natural parrs (> 650 ddph). Table S3. Significance level of differences (ANOVA) for the body length (on the top) and weight (below) among the groups.    Developmental dynamics of morphs body shape in common-garden experiment.
Principal component (PC) scores for body shape of SC (red dots) and DV (brown triangles) plotted together with growth (a). Deformation of the grid represents the change in body shape along the canonical root for older fry in the experiment (b) and for parrs in the wild (c).  Deformation of the grid represents the change in body shape along the first (b) and the second (c) canonical root.

Supplementary Files
This is a list of supplementary files associated with the primary manuscript. Click to