Reconsidering reproductive patterns in a model dissociated species, the red-sided garter snake: Sex-specific and seasonal changes in gonadal steroidogenic gene expression

Sex steroid hormones are powerful regulators of reproductive behavior and physiology in vertebrates, and steroidogenesis has distinct sex- and season-specific patterns ultimately dictated by the expression of key enzymes. Most comparative endocrinology studies, however, focus only on circulating levels of sex steroids to determine their temporal association with life-history events in what are termed associated reproductive patterns. The red-sided garter snake (Thamnophis sirtalis parietalis) is a notable exception; this species exhibits maximal sex behavior decoupled from maximal sex steroid production and gametogenesis in what is termed a dissociated reproductive pattern. And while this is true for male red-sided garter snakes and their production of testosterone, females have maximal estradiol production during peak breeding (spring) but only immediately after mating. Here, we demonstrate that expression of ovarian aromatase (conversion of androgens to estrogens) matches the established seasonal hormone pattern in females. Additionally, steroidogenic gene expression in the ovary is broadly reduced if not suppressed compared to the testis throughout the active year. Bizarrely, male red-sided garter snakes demonstrate an unexplained pattern of steroidogenic gene expression in the testis. StAR (import of cholesterol to steroidogenesis) is maximally expressed in spring, yet Hsd17b3 expression (conversion of androstenedione to testosterone) is highest in summer, with the latter matching the established summer peak in male testosterone. The function of elevated StAR in spring is unknown, but our results suggest a decoupling between maximal StAR expression and testosterone biosynthesis (Hsd17b3 expression). We also purport that the reproductive pattern binary should be reassessed given its lack of fit for many vertebrate species that demonstrate seasonal, mixed patterns of (a)synchrony between circulating sex hormones and reproductive behavior.

Sex steroid hormones are powerful regulators of reproductive behavior and physiology in vertebrates, and steroidogenesis has distinct sex-and seasonspecific patterns ultimately dictated by the expression of key enzymes. Most comparative endocrinology studies, however, focus only on circulating levels of sex steroids to determine their temporal association with life-history events in what are termed associated reproductive patterns. The red-sided garter snake (Thamnophis sirtalis parietalis) is a notable exception; this species exhibits maximal sex behavior decoupled from maximal sex steroid production and gametogenesis in what is termed a dissociated reproductive pattern. And while this is true for male red-sided garter snakes and their production of testosterone, females have maximal estradiol production during peak breeding (spring) but only immediately after mating. Here, we demonstrate that expression of ovarian aromatase (conversion of androgens to estrogens) matches the established seasonal hormone pattern in females. Additionally, steroidogenic gene expression in the ovary is broadly reduced if not suppressed compared to the testis throughout the active year. Bizarrely, male red-sided garter snakes demonstrate an unexplained pattern of steroidogenic gene expression in the testis. StAR (import of cholesterol to steroidogenesis) is maximally expressed in spring, yet Hsd17b3 expression (conversion of androstenedione to testosterone) is highest in summer, with the latter matching the established summer peak in male testosterone. The function of elevated StAR in spring is unknown, but our results suggest a decoupling between maximal StAR expression and testosterone biosynthesis (Hsd17b3 expression). We also purport that the reproductive pattern binary should be reassessed given its lack of fit for many vertebrate species that demonstrate seasonal, mixed patterns of (a)synchrony between circulating sex hormones and reproductive behavior.

Introduction
In vertebrates, both reproductive behavior and physiology are predominantly controlled by sex steroid hormones (1). Sex steroid hormones are primarily synthesized in the gonads from cholesterol then refined by key cytochrome P450 enzymes into androgens and/or estrogens to exert their effects in target tissues (2). Sex hormone signaling has manifold effects across reproductive functions in vertebrates, including gametogenesis, vitellogenesis, growth and development, and mating behavior (e.g., [3][4][5][6][7]. In circulation, sex steroid levels fluctuate seasonally and often coinciding with major reproductive events in species with distinct breeding seasons (e.g., 8). The relationship between sexual behavior, circulating sex steroid hormones, and gametogenesis is often described as one of three types of reproductive patterns. Most vertebrate species display an associated reproductive pattern, wherein mating behavior, maximal gametogenesis, and peak circulating sex steroid concentrations coincide temporally (9). Contrastingly, species with dissociated reproductive patterns exhibit mating behavior independent of or out of sync with maximal gametogenesis and sex steroid secretion (9) (Figure 1). Many vertebrates, however, do not fit strictly into either of these most common categories, especially species with more than one breeding season where both associated and dissociated patterns are observable in a given year (10). Third is the opportunistic pattern where gametogenesis and mating are either constitutive or paused until a mating opportunity is present, a pattern described in many vertebrates living in unpredictable environments (11)(12)(13). To determine the strength of the temporal association between reproductive behavior, gametogenesis, and sex steroid production, the majority of studies focus only on either circulating levels of sex steroids or gametogenesis and compare that to documented instances of mating in a given species, and, worse, most studies typically assume an associated reproductive pattern.
The red-sided garter snake (Thamnophis sirtalis parietalis) is a species that provides a unique opportunity to examine reproductive physiology as it is one of the most extensively studied non-model, seasonally breeding vertebrates with a dissociated reproductive pattern (14)(15)(16). After an 8-month winter dormancy in the northern extent of their range, red-sided garter snakes emerge in the tens of thousands from limestone hibernacula in spring, and males immediately begin displaying intense courtship behavior toward singly-emerging females across the 4-week breeding season (17). Though male courtship behavior is maximal in spring, the testes are in their most regressed state, circulating testosterone is either at or declining toward baseline, and neither castration nor supplemental testosterone influence courtship at this time (14,18). Courtship behavior is extinguished by the end of spring mating, and several months later the testes eventually recrudesce and circulating testosterone peaks in late summer (15,19). Similar to males, females are maximally receptive upon emergence in spring, yet the ovaries are typically regressed (20-22). However, unlike males, female garter snakes have significant variation in both receptivity and ovarian state in spring (23). Immediately post-mating, females produce a significant, short surge in circulating estradiol (E 2 ) from the ovaries and eventually undergo ovarian recrudescence, summer gestation, and parturition (23). Mating is not a requirement for FIGURE 1 Dissociated (A) vs. associated reproductive patterns (B) based on the (a) synchrony of behavior with gonadal activity. Redrawn from (9). Lincoln et al. 10.3389/fendo.2023.1135535 Frontiers in Endocrinology frontiersin.org ovarian recrudescence and parturition, though: females who do not mate in a given spring can still become gravid and give birth that year (22, 24). Further, E 2 is crucial for regulating female reproduction in this species, such as vitellogenesis, receptivity, and attractivity/ pheromone production (21,(24)(25)(26). During the late summer, prior to dispersal back to the den, gametogenesis occurs in both sexes in preparation for the next spring (9, 15, 17). Therefore, the red-sided garter snake is a textbook example of the dissociated reproductive strategy, making this species valuable for understanding cyclical shifts in gene expression underlying broad physiological phenomena (27). Dissociated reproduction in T. s. parietalis is considered an adaptative mechanism that resulted from the environmental pressures of living at extreme latitudes (9, 28). The chief limitation in such environments is temporal: the amount of time these ectotherms can be safely active during the year is approximately six months at most (May to October). Worse, weather events during the active period of the year can be random and severe, such as a mass mortality event one spring caused by a sudden snowstorm and freeze (29). The frequency and intensity of extreme climatic events are likely to increase due to climate change, increasing the need to understand the range and limitations of phenotypic plasticity (30,31). In the redsided garter snake, selective pressures shaping reproduction under such stochastic conditions with narrow mating windows likely favored seasonal modulation of gene expression in steroidogenesis (32). The enzymes responsible for regulating sex steroid production are sensitive to many external modulators that regulate their expression and activity and, thus, downstream hormonal products, especially in seasonally breeding vertebrates (e.g., 33). Indeed, seasonal variation in expression or activity of steroidogenic enzymes has been observed across mammals, birds, fish, and reptiles (e.g., [33][34][35][36][37], including some studies highlighting mismatches between enzyme expression and circulating sex steroids (e.g., 38,39). Because of the relationship that seasonal stimuli have with gonadal activity, environmental modulation of endocrine systems has been posited as a major generator of adaptive phenotypic plasticity, especially in species living in extreme conditions (e.g., 40). Therefore, measurements of seasonal changes in gene expression, not just the circulating output of said genes (e.g., sex steroids), can inform predictions about the physiological substrates shaped by environmental selection.
In this study, we aimed to determine whether mRNA expression of four key steroidogenic enzymes varies seasonally in accordance with established behavioral and hormonal patterns in both sexes of the red-sided garter snake. The four principal enzymes of interest are StAR, Cyp17a1, Hsd17b3, and Cyp19a1 (aromatase) as they comprehensively span the length of the step-wise steroidogenic pathway and cover four major steroid products therein (Figure 2), and these enzymes have also been examined in another excellent FIGURE 2 Simplified vertebrate steroidogenesis pathway focused on principle sex steroids (modified from (41)). Steroid hormones and precursors are in gray, enzyme gene names are in italics, and enzymes in bold indicate target genes for which primers could be designed. Genes are defined below pathway.
reptile model for comparative reproductive endocrinology, the green anole (Anolis carolinensis) (35). mRNA expression was quantified via qPCR using cDNA synthesized from both male and female gonadal tissue in the red-sided garter snake across three sequential phases during the active period of the same year. To date, the only steroidogenic enzyme studied in this species is aromatase, and even then, only the sexual and seasonal variation in protein within specific nuclei of the brain related to sex behavior (42). Therefore, this study offers insight at the molecular level on the seasonal steroidogenic capacity of the gonads in a model reptile with dynamic reproduction and life history.

Animal husbandry and seasonal conditions
A total of n=60 red-sided garter snakes (n=30 females, n=30 males) were collected from the wild at a known hibernaculum in Inwood, Manitoba, Canada (N 50°30' 25.5416", W 97°29' 49.7003", 277 m). The snakes were collected in May 2018 during the peak of the spring breeding season. All females were newly emerged from the hibernaculum (body temperatures~2-4°C), were extremely attractive to males (mated females become unattractive to males) and had no mating plug or remnants of a plug indicating they had not mated that spring. Females were further determined to be attractive based on whether they received repeated, daily courtship from sexually active, wild males in the den. All males were sexually active and vigorously courting females in the den at the time of capture.
Snakes (n=10 females, n=10 males) were euthanized for tissue collection at three different times: spring (May), summer (July), and fall (October). Spring snakes were euthanized at the field station after being separated by sex and held in outdoor arenas (e.g., 19) at ambient conditions for four days and provided water ad libidum. The remaining snakes were brought to the laboratory at James Madison University, Harrisonburg, Virginia, and housed in 37.8 L glass aquariums (n=3 per tank; sexes separated) equipped with basking lamps, shedding rocks, water bowls, hide objects, and shredded cardboard bedding (e.g., [43][44][45]. All were maintained within an environmental chamber programmed to simulate the summer conditions of the Interlake region in Manitoba (16 h:8 h L:D; 25:15°C ). The chamber was kept at these summer conditions from June through August, after which the program was transitioned to simulate fall field conditions (11 h:13 h L:D; 15:8°C) from September until October. The housing, temperature, and photoperiod conditions were established in previous studies on red-sided garter snakes (43,44), and these seasonal conditions are sufficient to cause shifts in gene expression in T. s. parietalis (45) and induce natural steroid hormone cycles that are established regulators of gametogenesis, sex behavior, and reproduction in this species (e.g., 14,15,44). Snakes were fed a mixed diet of earthworms and brook trout once per week throughout the summer, and water was provided ad libitum. Housing densities (n=3 for males; n=2 for females) enabled sufficient animal care to ensure all individuals were eating each week in compliance with animal care requirements, and snakes were also weighed and measured once per month to track condition. As part of the transition from summer to fall, the snakes were not fed the month prior to fall tissue collection (45); in the wild, red-sided garter snakes become aphagic by fall prior to winter dormancy (19). All procedures involving the use of vertebrate animals were approved by Manitoba Conservation and the IACUC of James Madison University.

Tissue collection
After euthanization, gonads (along with many other tissues) were collected at the three timepoints (spring, summer, fall) following protocols in Ashton et al. (45). Summer samples were collected one month after snakes had acclimated to simulated summer field conditions, and fall samples were collected one month after snakes had acclimated to fall conditions. All snakes were euthanized with an overdose of brevital sodium (10 mg/kg) by intraperitoneal injection followed by decapitation. Tissues were then excised, placed in cryotubes, snap frozen in liquid nitrogen, and stored at -80°C until use in RNA extraction.

Primer design, RNA extraction, and cDNA synthesis
Methods for all molecular work followed Ashton et al. (45), with modification to the RNA extraction and cDNA synthesis protocols. Primers for all control and target genes were designed in NCBI BLASTR. The target genes were: steroidogenic acute regulatory protein (StAR), P450 17a-hydroxylase/C17-20lyase (Cyp17a1), 17 beta-hydroxysteroid dehydrogenase type 3 (Hsd17b3), and cytochrome P450 19A1 (aromatase/Cyp19a1). Control genes were glyceraldehyde-3-phosphate dehydrogenase (Gapdh) and TATA box binding protein (TBP). See Table 1 for further details about each primer set. Prior to designing primers, protein sequences for each target gene were accessed across ten species, including the red-sided garter snake, and then aligned in Clustal Omega to determine regions of homology where designed primers in Thamnophis sirtalis had highest probability of authentically amplifying the target gene. Visualization of alignments also helped predict intron-exon boundaries for quantitative PCR primer design. All primer sets were validated using traditional PCR and gel electrophoresis before performing quantitative PCR (Figure 3).
To prepare tissues for RNA extraction, single, frozen, whole gonads were weighed then homogenized in liquid nitrogen using clean mortar and pestles (Cryo-Cup). Homogenized tissues were kept frozen (-80°C) until use in RNA extraction. All gonadal tissues were weighed prior to RNA extraction, however, only whole, anterior testis wet weight could be accurately reported; the entire string of ovarian follicles was not collected per female. Wet mass of the anterior testis (mean ± s.d.) for each season was: 65.2 ± 22.6 mg (spring), 96.6 ± 37.3 mg (summer), and 104.2 ± 43.5 mg (fall). Further, sperm were not visible in the ductus deferens in summer males but were visible in the fall males, indicating that spermiation had occurred by fall tissue collection (R. Parker, pers. obs.). An attempt was made to thaw and individually isolate follicles to determine mass, but follicles ruptured upon thaw and the smallest could not be individually isolated (~1mm diameter). Follicles were widely variable in size within and among individuals, but fall follicles were flaccid, visibly regressed, and filled with transparent fluid as is typical for wild females in fall in this population (e.g., 46). Spring and summer follicles were turgid and contained opaque fluid indicative of vitellogenin and other yolk-type proteins but were not fully vitellogenic (e.g., follicle diameter > 1 cm; 22).
RNA was extracted using Quick-RNA Miniprep Plus kits (Zymo) and followed the manufacturer protocol. RNA quality was assessed A B C FIGURE 3 All primers were validated in traditional PCR using tissue-specific cDNAs from red-sided garter snakes.  using a Gen5 Microplate Reader; absorbance ratio (260nm/280nm) and concentration (ng/mL) were recorded, and concentration was used to standardize input RNA concentrations to ensure parallel cDNA synthesis. Absorbance ratios averaged 2.035 (0.056 s.d.). Prior to cDNA synthesis, in-tube DNase I treatment (Invitrogen) was used to digest genomic DNA following the manufacturer protocol. cDNA was synthesized using the RevertAid First Strand kit (ThermoScientific) following the manufacturer protocol. A standard amount of sample RNA (1mg) was used in all cDNA syntheses and stored at -20°C until further use. cDNAs were validated using traditional PCR with primers for Gapdh.

Quantitative PCR
Primers were tested to confirm primer efficiencies using control tissue (pooled summer ovarian cDNA) before progressing to testing of experimental tissues. Experimental primer efficiencies were 110.3% for StAR, 96.7% for Cyp17a1, 110.6% for Hsd17b3, and 108.1% for Cyp19a1 ( Figure S1). All primer sets resulted in single products from melt curves. Methods followed those of Ashton et al. (2018). cDNAs were diluted to a 1:10 ratio for all qPCR reactions, and reactions were conducted following the SYBR Green protocol (ThermoFisher). Each sample was run in triplicate on 96-well plates in a CFX384 Touch Real-Time PCR System (Bio-Rad). The thermal cycles were as follows: activation at 95°C for 30 s, 45 cycles at 95°C (5 s each), and 60°C for 30 s followed by 31 s at 65°C. Melt curve analysis followed the completion of each run (65°C to 95°C in 0.5°C increments for 5 s). Samples across all three seasons (spring, summer, and fall) and each sex were included on every 96-well plate to control for interassay variation. Two target genes (StAR and aromatase; Cyp17a1 and Hsd17b3) and both control genes (Gapdh and TBP) were run per plate. All plates also included an internal control (pooled [n=6] summer ovarian cDNA).

Data analysis and visualization
For determining the effect of season on gene expression, one-way ANOVAs were conducted for each gene followed by post-hoc comparisons (Tukey tests). For determining the effect of sex and season on gene expression, two-way ANOVAs were conducted for each gene to determine any interactions prior to sex-within-season pairwise comparisons (Tukey tests). Correlations in gene expression were first targeted from a Pearson correlation matrix then analyzed using simple linear regression. Analysis and visualization were carried out in SigmaPlot 14 (Systat Software Inc.). For all statistical analyses, alpha was set at 0.05 and both significant (p<0.05) and marginal (0.05<p<0.1) differences were reported.
In ovary, there was significant seasonal variation in Cyp17a1 and aromatase expression but not StAR or Hsd17b3 ( Figure 5). Cyp17a1 expression was seasonally variable (F 2,28 = 4.71, p=0.018), significantly elevated in summer compared to spring (q=4.28, p=0.015), and marginally higher in fall compared to spring (q=2.68, p=0.069) ( Figure 5B). Summer and fall were not different (q=1.66, p=0.24). Aromatase expression in the ovary was also seasonally variable (F 2,28 = 3.38, p=0.049), but there was a different seasonal pattern compared to Cyp17a1 ( Figure 5D). Aromatase expression was significantly elevated in spring compared to fall (q=3.02, p=0.042) and marginally higher in summer compared to fall (q=3.30, p=0.068). Summer and spring did not differ (q=0.36, p=0.79). No other comparisons were statistically significant.

Discussion
Contrary to expectations based on the timing of spermatogenesis and the known peak in circulating testosterone for male red-sided garter snakes, StAR was maximally expressed in the testes during the spring mating season. In vertebrates, StAR is widely purported as the rate-limiting enzyme in the steroidogenic pathway (47)(48)(49)(50). The next enzyme analyzed in the steroidogenic pathway, Cyp17a1, also showed highest expression in spring testis, though this was not statistically significant. These results match what would be predicted for a species with an associated reproductive pattern: our observations suggest active investment in steroidogenesis when mating behavior is maximal in male Thamnophis sirtalis parietalis. However, it is paradoxical in the context of a species exhibiting a dissociated reproductive strategy. Male red-sided garter snakes do not undergo spermatogenesis during the spring breeding season, nor do they have significantly elevated testosterone relative to the rest of the active season. Moreover, male courtship behavior in this species is not induced by elevated testosterone (18); rather, it is winter dormancy that alters neuronal steroidogenesis to upregulate sexually dimorphic regions of the hypothalamus and activate sex behavior (42). Peak spermatogenesis and increased circulating testosterone will not occur until late summer, thus the enzymes in this pathway should be upregulated in sync with or tightly coupled to these events (15,51). Our data only partially support this: Hsd17b3 was elevated in summer and remained high into fall. Though StAR expression was highest in spring, expression of StAR in summer was strongly, positively correlated with the peak of Hsd17b3 expression, indicating that steroidogenesis in summer is maximizing testosterone production in sync with gametogenesis as predicted by previous research (9, 15, 52). Lastly, we also observed greatest individual variation in Hsd17b3 expression compared to the other steroidogenic genes, which suggests this gene may be a source of phenotypic variation in male red-sided garter snakes.
At emergence from low-temperature dormancy in spring, the testes of male Thamnophis sirtalis parietalis are regressed, and thus no major steroidogenic activity was predicted. However, our results strongly suggest that steroidogenesis is occurring-and much earlier than expected. Perhaps males are front-loading precursor androgens in anticipation of summer gametogenesis. Circulating testosterone levels can be elevated in male garter snakes upon emergence from hibernation, but as soon as they actively begin courting females at the den these residual androgens are cleared from circulation (15, 16, 53). The testis may be able to store weaker androgens (e.g., DHEA, androstenedione) as substrates for immediate conversion to more potent androgens (e.g., testosterone) once downstream enzymes are upregulated/active. Thus, steps in androgenesis may occur seasonally in male red-sided garter snakes. Tissue-specific measurement of sex steroids would illuminate this idea, as has been demonstrated recently in brown anoles (Anolis sagrei) (39). Our results also beg the question of how environmental stimuli may act differentially on specific enzymes within the same tissue, analogous to temperaturedependent allozymes in cold-adapted fish species (54). Suites of steroidogenic genes may be tuned to specific temperature optima in red-sided garter snakes to result in distinctly seasonal, step-wise regulation of steroidogenic enzymes. Future research in our laboratory will address this possibility. Female red-sided garter snakes appear to be limited in the synthesis of estrogens based on our results, with maximal expression of the terminal enzyme, aromatase, occurring during the predicted peak in estrogen synthesis (spring). The rate-limiting enzyme in steroidogenesis, StAR, was maximally expressed in spring and was positively correlated with aromatase expression, suggesting estrogen biosynthesis is most active during the mating season. Females produce a mating-induced E 2 surge in spring that is aromatase-dependent (23). Further, the E 2 surge can be induced in ex vivo preparations of ovarian fragments from female garter snakes in spring via prostaglandin treatment (55,56). Seasonal, anticipatory upregulation of ovarian aromatase would be beneficial in maximizing E 2 synthesis given the important roles of estrogen signaling in female reproductive physiology. However, the broadscale reduction in aromatase expression we observed is unusual. The females in our study were unmated; how and to what extent mating in spring may shift gene expression in the ovary is an outstanding question. Perhaps we have observed facultative quiescence in steroidogenic gene expression, especially as there are many examples of social behavior inducing changes in gene expression at each level of the vertebrate hypothalamic-pituitary-gonadal axis (reviewed in 56). By extension, female garter snakes may thus have significant phenotypic plasticity in steroidogenesis that is dependent upon spring mating. As with the testes, ovarian Cyp17a1 was upregulated in summer compared to spring. This did not match our expectation that StAR would be the major enzyme with seasonal variation. Previous work in this species indicated that steroidogenesis was mating dependent: E 2 peaks immediately but briefly after mating (approximately 6 hours postcopulation; returns to baseline by 24 hours) and then remains relatively low for the remainder of the year (23,55). Through the examination of upstream mechanisms (steroidogenic enzymes), we have partially explained one limitation for female production of estrogens in red-sided garter snakes.
The seasonal pattern of steroidogenesis in females is likely tied to the unresolved mechanisms regulating vitellogenesis in garter snakes. Vitellogenesis is crucial for gestation and parturition of offspring because vitellogenins serve as lipid transport proteins and egg-yolk precursors that facilitate lipid circulation for follicular development by the ovaries (24, 57). The exact mechanisms that control and induce vitellogenesis in garter snakes are unclear, as unmated females in a given year can become vitellogenic without experiencing the matinginduced E 2 surge (23). Instead, it was hypothesized that these females mated the previous fall and experienced a post-mating E 2 surge, initiated then halted vitellogenesis prior to winter dormancy, and eventually resumed vitellogenesis the next spring (20, 23). Further, in the ovarian ex vivo experiment mentioned above, E 2 synthesis was significantly variable: the ovaries of some females were indicated as "unresponsive" due to no change in E 2 secretion following A B D C FIGURE 5 Normalized steroidogenic gene expression in ovarian tissue from female red-sided garter snakes across three seasons (spring, summer, and fall). Genes: prostaglandin treatment (56). Other studies have shown that body condition at emergence may also predict vitellogenesis; however, this was not consistent as some females in poor condition still became vitellogenic (58). These separate studies help explain the general variation in steroidogenic gene expression we observed between individual females for both Cyp17a1 and aromatase, and this variation was not as pronounced for StAR or Hsd17b3.
Female pheromone production in garter snakes is directly associated with estrogen availability. Across the spring breeding season, females become lose their attractivity because of a reduction in pheromone quality and quantity associated with declining circulating estradiol (59). If ovariectomized in spring, females have reduced attractivity the next spring (21). Systemic aromatase inhibition across the active year (spring to fall) reduces female pheromone production and attractivity the next spring (26). Further, male garter snakes given estradiol implants produce abundant female sex pheromone and are attractive to wild males, an effect lost once the implants are removed which indicates an activational role for estradiol (25). However, the seasonal relationship between ovarian status, circulating estradiol, and female pheromone production is still unclear. Female garter snakes ovariectomized in fall or during hibernation remain attractive the next spring (60), mating renders females unattractive despite the reflexive surge in estradiol (20, 23), and low temperature dormancy increases pheromone quality and quantity (61). Female garter snakes may thus have constitutive production of sex pheromone from a skin that is organized at some point during early life to be relatively independent of circulating estradiol (25).
Our study is confined by our focus on mRNA expression, which does not always correspond with protein production or activity. For example, at least one previous study indicated a mismatch between mRNA expression and gene activity for steroidogenic genes (38). Our analyses only utilized cDNA, but protein quantification and functional assays should also be pursued to determine the outcomes of the observed seasonal and sexual patterns of mRNA expression. Measurements of steroidogenic conversion of specific hormones in gonadal and extra-gonadal tissues would be particularly informative. We observed substantial variation in gene expression among individuals, especially females; as stated above, that variation is biological and informative. The gonadotropins (LH, FSH) strongly influence the rate of gonadal steroidogenesis (e.g., 62), but molecular identification and purification of gonadotropins in squamates (snakes, lizards) has been challenging and, to our knowledge, remains elusive in snakes (e.g., 63, 64). The gene for FSHb is annotated for Thamnophis sirtalis (NCBI reference sequence NW_013658033.1), and recent work has demonstrated a key role for GnRH in neural plasticity of garter snakes (65); thus, additional questions could be pursued across the hypothalamic-pituitary-gonadal axis to understand the relative influence of tropic hormones on steroidogenesis. Lastly, we did not focus on quantifying steroid hormone levels in our study because the patterns of testosterone and estradiol production are well-established in Thamnophis sirtalis parietalis; however, a survey of additional sex steroids, especially DHEA, is warranted to support and extend our conclusions.
Our data suggest that the reproductive pattern of red-sided garter snakes, as is the case for many other snake species, is a mix of dissociation and association between steroidogenesis and mating activity. Many pitviper species show mixed reproductive patterns, with spermatogenesis occurring in late summer but breeding seasons varying per species and even population (10,66). Eastern copperheads (Agkistrodon contortrix) and pygmy rattlesnakes (Sistrurus miliarius) have bimodal breeding seasons but only unimodal gametogenesis; thus, these species exhibit both associated and dissociated reproductive patterns in a single annual cycle (67, 68). Contrastingly, European asps (Vipera aspis) and western diamondback rattlesnakes (Crotalus atrox) have unimodal increases in circulating sex steroids uncoupled from their unimodal breeding seasons (69, 70). Lastly, the type of reproductive pattern can even be sex-specific, as with the cottonmouth (Agkistrodon piscivorus) (66). Our results further call into question the associated-dissociated binary when describing hormone-behavior patterns in vertebrate reproduction, as others have suggested (71). Such labels are restrictive, and it may be more pertinent to treat reproductive patterns as a continuum rather than a binary.
Our study shows that steroidogenic enzyme expression in the gonads is seasonally and sexually variable in red-sided garter snakes. A notable question centers on the importance of winter dormancy in the regulation of the steroidogenic pathway. The growing body of evidence suggests that beyond the central control of male courtship behavior, prolonged low-temperature dormancy is a critical period of significant physiological shifts in the annual cycle of red-sided garter snakes (e.g., 44, 61,65,72). How temperature is influencing gene expression in the gonads also holds implications for climate changenamely how reduced periods of dormancy and subsequent, diffuse spring breeding seasons may relax selection on steroidogenesis in this temperature-dependent system.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material further inquiries can be directed to the corresponding author/s.

Ethics statement
The animal study was reviewed and approved by JMU Institutional Animal Care and Use Committee.
Author contributions JL, MB, and MP conceived of and designed the study. JL and MB conducted all of the laboratory work and collected and organized all data. HR maintained the animals and collected tissues. MP performed statistical analyses. MB and JL wrote the first draft of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.

Funding
The authors would like to thank the Honors College, Hillcrest Foundation, and the College of Science and Mathematics at James Madison University for funding. This research was also funded in part by 4-VA, a collaborative partnership for advancing the Commonwealth of Virginia.