Vomeronasal Receptors Associated with Circulating Estrogen Processing Chemosensory Cues in Semi-Aquatic Mammals

In numerous animals, one essential chemosensory organ that detects chemical signals is the vomeronasal organ (VNO), which is involved in species-specific behaviors, including social and sexual behaviors. The purpose of this study is to investigate the mechanism underlying the processing of chemosensory cues in semi-aquatic mammals using muskrats as the animal model. Muskrat (Ondatra zibethicus) has a sensitive VNO system that activates seasonal breeding behaviors through receiving specific substances, including pheromones and hormones. Vomeronasal organ receptor type 1 (V1R) and type 2 (V2R) and estrogen receptor α and β (ERα and ERβ) were found in sensory epithelial cells, non-sensory epithelial cells and lamina propria cells of the female muskrats’ VNO. V2R and ERα mRNA levels in the VNO during the breeding period declined sharply, in comparison to those during the non-breeding period, while V1R and ERβ mRNA levels were detected reversely. Additionally, transcriptomic study in the VNO identified that differently expressed genes might be related to estrogen signal and metabolic pathways. These findings suggested that the seasonal structural and functional changes in the VNO of female muskrats with different reproductive status and estrogen was regulated through binding to ERα and ERβ in the female muskrats’ VNO.


Introduction
Information exchange can transmit information between different animals, which is essential for animal social behavior and physiology involving sexual recognition, individual survival and mating [1][2][3]. Chemical communication is an essential way of information exchange between individuals, which is usually less influenced by external environmental factors [4,5]. In mammalians, the olfactory and vomeronasal systems can receive a chemosensory signal that mediates their particular pattern of behaviors [6]. The main olfactory system is known to be closely related to the limbic system and has a crucial role in emotions and behavior. The vomeronasal system (VNS), considered an accessory olfactory system, which is associated with specific behaviors and hormonal responses to pheromone stimuli, is related to reproductive behavior, for instance, sexual attraction, maternal behavior or aggression [7][8][9]. The occurrence and development of the vomeronasal system varies greatly among mammals, depending on the environment in which they live [10]. Most rodents and land mammals have a VNS, but marine mammals and birds' VNS disappeared [11]. It would appear that the vomeronasal system is of less use for airborne vertebrates and marine mammals [12]. In most mammals, the vomeronasal organ (VNO) comprises three parts: ducts, tissue and cartilage [13,14]. The vomeronasal duct (VND) is a bilateral mucus-filled superior duct at the bottom of the nasal septum and is typically crescent shaped with a thick sensory epithelium (SE) in the middle and ciliated non-sensory epithelium (NSE) on the sides [15]. The VND is surrounded by cartilage and lating estrogen could regulate female muskrat response to different chemical signals and individual behaviors in different reproductive statuses, which might play a significant part in the physiological foundation of semi-aquatic mammals' response to chemosensory cues.

Morphological Observations and Histological Features of the VNO of Female Muskrats
Morphological observations and histological features of the female muskrats' VNO are presented in Figure 1. The VNO of female muskrats was observed as a paired tubular organ, which is situated at the basal part of the nasal septum (NS) ( Figure 1A,B). Therefore, odor molecules are received through the nasal or oral cavity by the VNO ( Figure 1B). Histologically, the VNO of the female muskrat is crescent-shaped and surrounded by the vomeronasal bone (Vo) on the periphery and separated by the NS in the middle ( Figure 1C). The VNO consists of a vomeronasal duct (VND) and the lamina propria, including a vomeronasal gland (VG), blood vessel (BV), smooth muscles and connective tissue ( Figure 1C). In both periods, the vomeronasal epithelium is separated into the sensory epithelium (SE) at the protruding part of the VND and the non-sensory epithelium (NSE) at the concave part of the VND ( Figure 1D,E). The epithelium of the VNO consists of three parts: long columnar ciliated supporting cells (SC), sensory cells (SC) and basal cells (BC) near the lamina propria (LP) ( Figure 1D,E). The ratio of surface epithelium to total NSE was calculated and is shown in Figure 1F. Based on the boxplot, the proportional contribution of the surface epithelial thickness was significantly higher in May than in December. remains unclear. Our study investigates the expression patterns of VRs and estrog ceptors (ERs) in the female muskrats' VNO to reveal that the morphological chan the VNO associated with circulating estrogen could regulate female muskrat respo different chemical signals and individual behaviors in different reproductive st which might play a significant part in the physiological foundation of semi-aquatic mals' response to chemosensory cues.

Morphological Observations and Histological Features of the VNO of Female Muskrat
Morphological observations and histological features of the female muskrats are presented in Figure 1. The VNO of female muskrats was observed as a paired t organ, which is situated at the basal part of the nasal septum (NS) ( Figure 1A,B). fore, odor molecules are received through the nasal or oral cavity by the VNO (Figu Histologically, the VNO of the female muskrat is crescent-shaped and surrounded vomeronasal bone (Vo) on the periphery and separated by the NS in the middle ( 1C). The VNO consists of a vomeronasal duct (VND) and the lamina propria, inclu vomeronasal gland (VG), blood vessel (BV), smooth muscles and connective tissue ( 1C). In both periods, the vomeronasal epithelium is separated into the sensory epith (SE) at the protruding part of the VND and the non-sensory epithelium (NSE) at th cave part of the VND ( Figure 1D,E). The epithelium of the VNO consists of three long columnar ciliated supporting cells (SC), sensory cells (SC) and basal cells (BC the lamina propria (LP) ( Figure 1D,E). The ratio of surface epithelium to total NS calculated and is shown in Figure 1F. Based on the boxplot, the proportional contri of the surface epithelial thickness was significantly higher in May than in Decembe

Immunohistochemical Locations of V1R, V2R, ERα and ERβ in the VNO of Muskrats
The immunohistochemical localizations of V1R, V2R, ERα and ERβ in the VNO of female muskrats are depicted in Figure 2. The immunoreactivity of V1R and V2R was localized mainly in the cytoplasm of the SE and NSE in the VNO in May and December (Figure 2A-D). The semi-quantitative positive staining of V1R was remarkedly higher in December than that of May ( Figure 3A, Table 1), while the semi-quantitative positive staining of V2R was remarkedly higher in May than that in December ( Figure 2C,D and Figure 3A, Table 1). The immunoreactivity of ERα and ERβ were localized mainly in the nucleus of the SE and NSE of the VNO in May and December ( Figure 2I-L). The positive staining of ERα in May was significantly higher than that in December ( Figure 2I,J and Figure 3A, Table 1). Meanwhile, relatively high positive staining for ERβ was observed in December ( Figure 2K,L and Figure 3A, Table 1). In the negative controls, there was no signal visible ( Figure 2E-H,M-P). (D,E). The error bars represent means ± SEM (n = 6, each period). * Statistically significant values (** p < 0.01).

Immunohistochemical Locations of V1R, V2R, ERα and ERβ in the VNO of Muskrats
The immunohistochemical localizations of V1R, V2R, ERα and ERβ in the VNO of female muskrats are depicted in Figure 2. The immunoreactivity of V1R and V2R was localized mainly in the cytoplasm of the SE and NSE in the VNO in May and December (Figure 2A-D). The semi-quantitative positive staining of V1R was remarkedly higher in December than that of May ( Figure 3A, Table 1), while the semi-quantitative positive staining of V2R was remarkedly higher in May than that in December (Figures 2C,D and 3A, Table 1). The immunoreactivity of ERα and ERβ were localized mainly in the nucleus of the SE and NSE of the VNO in May and December ( Figure 2I-L). The positive staining of ERα in May was significantly higher than that in December (Figures 2I,J and 3A, Table 1). Meanwhile, relatively high positive staining for ERβ was observed in December (Figures 2K,L and 3A, Table 1). In the negative controls, there was no signal visible ( Figure 2E-H,M-P).    The immunohistochemical staining was determined as negative (−), positive (+), strongly (++) and very strongly positive (+++). Staining that was weak but higher than that of the con set as positive (+). The highest intensity staining was set as very strongly positive (+++). A intensity between + and +++ was set as strongly positive (++). No signal was set as negativ means sensory epithelium; NSE means non-sensory epithelium.

The mRNA Expressions of v1r, v2r, erα and erβ in the VNO of the Muskrats
The relative mRNA expression levels of v1r, v2r, erα and erβ were examined time quantitative (qRT-PCR) ( Figure 3B). The levels of v1r and erβ transcripts in th of female muskrats declined remarkedly from December to May (p < 0.01) (Figu However, the mRNA expression levels of v2r and erα in the VNO of female musk creased significantly from May to December (p < 0.01) ( Figure 3B). The ratios of V ERs were calculated according to the mRNA expression levels in the VNO of female rats ( Figure 3C,D). The ratio of V1R to V2R was notably higher in December comp May (p < 0.01) ( Figure 3C). However, the ratio of ERα to ERβ was higher in May th in December (p < 0.01) ( Figure 3D).  The immunohistochemical staining was determined as negative (−), positive (+), strongly positive (++) and very strongly positive (+++). Staining that was weak but higher than that of the control was set as positive (+). The highest intensity staining was set as very strongly positive (+++). A staining intensity between + and +++ was set as strongly positive (++). No signal was set as negative (−). SE means sensory epithelium; NSE means non-sensory epithelium.

The mRNA Expressions of v1r, v2r, erα and erβ in the VNO of the Muskrats
The relative mRNA expression levels of v1r, v2r, erα and erβ were examined by realtime quantitative (qRT-PCR) ( Figure 3B). The levels of v1r and erβ transcripts in the VNO of female muskrats declined remarkedly from December to May (p < 0.01) ( Figure 3B). However, the mRNA expression levels of v2r and erα in the VNO of female muskrats increased significantly from May to December (p < 0.01) ( Figure 3B). The ratios of VRs and ERs were calculated according to the mRNA expression levels in the VNO of female muskrats ( Figure 3C,D). The ratio of V1R to V2R was notably higher in December compared to May (p < 0.01) ( Figure 3C). However, the ratio of ERα to ERβ was higher in May than that in December (p < 0.01) ( Figure 3D).

Concentration Profiles of 17β-Estradiol in the Female Muskrats
Female muskrats' circulating levels of 17β-estradiol during both seasons were measured using ELISA ( Figure 3E). The circulating 17β-estradiol concentration level in female muskrats reached the maximum in May and dropped remarkedly in December (p < 0.01).

Transcriptome Data Characterization
All transcriptome experiments were performed on both periods according to the data using sequencing with Illumina. Figure 4 displays the findings of the analysis of the transcriptome data in female muskrats' VNO. A total of 18,339 gene fragments were detected to be co-expressed during both seasons. Respectively, 689 and 967 gene fragments were observed in May and December ( Figure 4A). The volcano plot depicted that 883 differentially expressed genes (DEGs) were detected in the VNO, of which 554 genes were down-regulated and 329 genes were up-regulated ( Figure 4B). The results of the DEG enrichment analysis are demonstrated in Figure 4C-F. The significantly differentially expressed pathways consisted of epithelial tube morphogenesis, mesenchyme development, regulation of epithelial proliferation and gland development, which includes some epidermal growth factors, such as the Sox family, fgfr2, epha2 and snal2; these were identified via GO enrichment analysis and are shown in the chord diagram ( Figure 4C). Bubble plots for the top 10 GO-enriched pathways are shown in Figure 4D. Consistently, GO enrichment analysis of all DEGs revealed that they are composed of the epidermal morphology pathway and gland development pathway ( Figure 4D). There were six pathways related to odorant receptors: olfactory receptor activity, sensory perception of smell, odorant binding, pheromone receptor activity, response to stimulus and detection of chemical stimulus involved in sensory perception of smell, identified via cluster analysis of DEGs. The cluster of up-regulated genes comprised aphrodisin, v1r, v2r and obp2b, and down-regulated genes included or10h2, lcn3 and vn1r17p ( Figure 4E). Through performing KEGG enrichment analysis of DEGs, enriched pathways included the metabolic pathway, steroid hormone biosynthesis and the estrogen signaling pathway. Olfactory transduction is demonstrated in Figure 4F, implying that the steroid hormone and metabolic pathways may be involved in the regulation of seasonal differential expression of vomeronasal organs in female muskrats.

Molecular Docking Stimulation
The top five compounds detected via liquid and gas chromatography-mass s trometry (LC-MS and GC-MS) in male muskrats' scent glands are listed in Table 2 as viously reported [34][35][36][37]. In brief, the top five compounds with the highest concentra identified via the GC-MS method were selected. In the LC-MS method, the top five mone-related compounds with the highest concentration were selected. The bindin finities and interactions between the selected molecular candidates and the V1R and main protease are shown in Figure 5. The compounds' binding affinities were scored ing AutoDock Vina (1.1.2 for Windows, San Diego, CA, USA), and the value with hig negative is shown in Table 3. The highest negative values were selected after calcula the binding affinities. Among them, all binding ability was less than −4.0 KJ/mol; the b ing ability of V2R and Androsterone sulfate was the strongest, and the binding affi

Molecular Docking Stimulation
The top five compounds detected via liquid and gas chromatography-mass spectrometry (LC-MS and GC-MS) in male muskrats' scent glands are listed in Table 2 as previously reported [34][35][36][37]. In brief, the top five compounds with the highest concentration identified via the GC-MS method were selected. In the LC-MS method, the top five hormone-related compounds with the highest concentration were selected. The binding affinities and interactions between the selected molecular candidates and the V1R and V2R main protease are shown in Figure 5. The compounds' binding affinities were scored using AutoDock Vina (1.1.2 for Windows, San Diego, CA, USA), and the value with highest negative is shown in Table 3. The highest negative values were selected after calculating the binding affinities. Among them, all binding ability was less than −4.0 KJ/mol; the binding ability of V2R and Androsterone sulfate was the strongest, and the binding affinity was −8.1 KJ/mol within interacting residues SER-84, THR-83 and GLU-14 (Table 3). Meanwhile, the binding ability of V2R and Muscone was −7.3 KJ/mol within interacting residues TYP-57. The binding ability of V1R, Muscone and Androsterone sulfate were −5.5 KJ/mol and −6.6 KJ/mol within interacting residues GLU-394 and GLN-217, respectively (Table 3). was −8.1 KJ/mol within interacting residues SER-84, THR-83 and GLU-14 (Table 3). Meanwhile, the binding ability of V2R and Muscone was −7.3 KJ/mol within interacting residues TYP-57. The binding ability of V1R, Muscone and Androsterone sulfate were −5.5 KJ/mol and −6.6 KJ/mol within interacting residues GLU-394 and GLN-217, respectively (Table 3).

Discussion
In mammals, due to multiple behavioral and reproductive strategies, the morphological characteristics of the VNO as well as the reaction to pheromones in the VNS are different in different species [2,38,39]. The VNO of hedgehogs consists of hyaline cartilage, a capsule and epithelial lumen soft tissue with an incisive duct directly into the nasal and oral epithelium, suggesting that odorants are mainly received through the nasal and oral cavity in the VNO of hedgehogs [14]. The brown bear's VNO comprises the cartilage vomeronasal sensory epithelia (VNE), NSE, VG and soft tissue including the opening of the vomeronasal lumen into the mouth, raising the possibility that bears employ the VNO to detect pheromones in the environment [13]. In giraffes, the SE and NSE were described in the lateral and medial sections of the VNO lumen, which is covered by veins and numerous thin-walled vessels, indicating that the VNO has a potential role in regulating female giraffes' function [40]. On the other hand, it is reported that individual status may affect the morphological characteristics of the vomeronasal epithelium [38,41]. For example, in rats, sex and age influence the morphological characteristics of the VNE [42]. In black bears, the mean thickness of VNE shows an age dependence, suggesting that VNO morphology may be affected by age and gender [13]. Our morphological characteristics and histological results indicated that the female muskrat's VNO was well-developed, indicating that this organ was structurally able to receive external pheromones through the nasal or oral cavity, like most mammals. Consistently, the proportion of vomeronasal epithelium in female muskrats' VNO in May was higher than that in December. Meanwhile, transcriptomic GO analysis depicted that DEGs were enriched into the epithelial morphological differential pathway. Therefore, the current results provided more evidence that seasonal variation of VNO morphology and histology in the female muskrat might closely be mediated by external chemical signals.
It is known that the VNO mainly binds different molecular chemicals through VRs to efficiently participate in response to the sensory environment as well as mediate different behavioral responses, for example, aggression and conspecific attraction [12,43]. Previous studies showed that odor directly induced changes in the expression of neurons and activated corresponding neuronal activation in the VNO epithelium, suggesting that the mammalian VNS can activate ligands of different VRs to detect odors secreted by other individuals [2,44]. Our results identified that the immunohistochemical expressions as well as gene expressions of V1R and V2R were shown in the female muskrats' VNO, indicating that V1R and V2R might regulate the functions of the VNO in female muskrats. The immunoreactivity and transcriptional levels of V2R were notably higher in May compared to December. However, the expression levels of V1R were higher in December. Consistently, different chemical molecules had different binding abilities with V1R and V2R through the molecular docking analysis. The liquid and gas phase molecules androsterone sulfate and muscone had strong binding ability with V2R and V1R, respectively, indicating that V1R and V2R bind different ligands to respond to different environments. Differences in the structures and functions of V1R and V2R had also been found in other animals [12,45]. According to previous studies, differences in the expression of V1R and V2R in vertebrates may be related to the environment in which the organism lives [46,47]. In semiaquatic beavers, the V1R and V2R genes were described according to gene repertory analysis, suggesting that volatile and non-volatile substances can be received by beavers using the vomeronasal system [44]. In the rat and mouse, V1R and V2R counts were dynamic during different developmental periods, indicating that VRs might fulfill different behavior-related functions within the developmental process in rodents [48]. Several V2R genes have been found in both frogs and red-legged salamanders, while V1R repertoires were reduced, demonstrating that amphibians are well suited to terrestrial and aquatic habitats and could use water-soluble and volatile substances as chemical cues via various VRs [49,50]. Our results that the increase of V2R was expressed in the female muskrats' VNO during the breeding period indicate that V2R is suspected of acting as the domain subtype VRs in the female muskrats' VNO and playing a crucial role in chemical communication, including the recognition of water-soluble molecules and chemicals floating on the water surface. Thus, differential expression patterns of V1R and V2R during both periods in the female muskrats' VNO might be consistent with the amphibious living environment of the muskrats via binding different molecular ligands responsible for essential functions.
In mammals, the VNO is regarded as a sensory organ and has been reported to be highly sensitive to steroid pheromones, which would be involved in mating and sexual behavior [8,51]. Recently, several studies found through immunohistochemical staining that mRNA levels of steroid hormone transport proteins and steroid hormone receptors are expressed in the SE of rodents' VNO, indicating that steroid hormones projected the limbic system through the VNO [25,52]. In terrestrial salamanders, circulating steroid hormone levels were found to stimulate the morphological characteristic of chemosensory cues in the VNO during the mating season [53]. In this study, ERs were colocalized with VRs in the NSE and SE cells of the VNO, indicating possible cross talk between the estrogen signaling pathway and chemosensory cue processing. Circulating peptides or steroid hormones could modulate odorant responses to alter the behavior of animals [25,54]. In particular, 17β-estradiol may be synthesized locally in the VNO to readily regulate the VSN-mediated odor response [25]. Although the expression pattern of factors related to the estrogen signaling pathway was not completely consistent with VRs, it cannot be excluded that estrogen may be involved in regulating VR expression through modulating the expression of different types of ERs. Furthermore, molecular docking results show that V1R and V2R have different affinities for different ligands. It has been reported that estradiol may mediate the action of the VNO through reducing the response of calcium ions in urine to steroid sulfate. In the most recent research, nuclear and membrane estrogen receptors can synergistically regulate estrogen's numerous functions [55]. In addition to its regulatory role through nuclear estrogen receptors, estrogen can also rapidly regulate estrogenic signaling through membrane estrogen receptors, particularly in the nervous system. Among them, membrane-associated estrogens α and β, as well as GPER1 and ER-X, can act through the second messengers, thus acutely regulating estrogenic signaling throughout the organism [55]. Rapid responses to external stimuli via membrane estrogen receptors can initiate the transcription of diverse estrogen-regulated genes. The role of estrogen receptors in the vomeronasal organ of female muskrats may also be through the reception of membrane estrogen receptors, resulting in rapid signaling and reception in the vomeronasal organ. Therefore, using the semi-aquatic muskrat as a model, our results elucidate that circulating estrogen may cooperatively regulate the differential expression of VRs through the locally expressed ERs, thus affecting chemosensory processing in the VNO. Together, the results indicate that the female muskrats' VNO was a main target of estrogen, and estrogen might have a vital role in pheromone-dependent behavioral alterations through regulating the VRs and ERs, including more activity in receiving chemical signals to prepare for reproduction.

Animals
The Beijing Forestry University Experimental Animal Ethics Committee gave its approval to all animal experiment procedures. All female muskrats were obtained in May (the breeding, B, n = 12) and December (the non-breeding season, NB, n = 12) from Xinji Muskrats Farm in Hebei Province, China. The heads were quickly dissected from female muskrats after being paralyzed using diethyl ether. The samples were promptly frozen at −80 • C for the transcriptome analysis and the molecular experiments; other tissues were fixed in paraformaldehyde for 24 h and then switched to 70% alcohol for the following procedures. For the subsequent hormonal study, serum was isolated from blood by means of centrifuging at 3000× g for 20 min. The animal study protocol was approved by the Experimental Animal Ethics Committee of Beijing Forestry University. Approval Code EAWC_BJFU_2021005. Approval Date: 11 March 2021.

Macroscopic Anatomy and Histological Process
One of the female muskrats' heads was dissected freshly for macroscopic anatomical evaluation of the VNO prior to the histological process outlined below. In accordance with the accepted procedure, the VNO was taken out of the nasal cavity and put into paraffin. In brief, the samples were calcified and dehydrated in succession of alcohol and xylene series and embedded in paraffin wax. The samples were subsequently sliced into transverse sections that were 5 µm thick for the general histological and immunohistochemical observations. Hematoxylin-eosin (H&E) staining was used to perform histological analysis in some sections.

Immunohistochemical Process
To lessen background staining brought on by the second antibody, 10% normal goat serum was incubated with the serial slices of the VNO tissues. Then, the primary antibody (1:200 dilution) against VMN1R41 (ARP96665_P050, Aviva Systems Biology, San Diego, CA, USA), VN2R1P (OASG07552, Aviva Systems Biology, San Diego, CA, USA), AR (sc-7305, Santa Cruz Biotechnology, Santa Cruz, CA, USA), ERα (sc-542, Santa Cruz Biotechnology, Santa Cruz, CA, USA), ERβ (sc-390243, Santa Cruz Biotechnology, Santa Cruz, CA, USA) for 12 h at 4 • C. Instead of using the main antibodies, normal rabbit IgG was applied at a ratio of 1:10000 to the control sections. All sections were incubated with a second antibody, goat anti-mouse or goat anti-rabbit IgG, conjugated using a Kit as we described before, after which the coloring was observed using DAB as the chromogen. The slides incubated with anti-VMN1R41 and anti-VN2R1P were counterstained with hematoxylin solution. Analysis was performed using a relative quantitative positivity score (ranging from negative − to very strong positive +++) [31]. The positive optical density of all cells stained with target antibodies in the muskrat VNOs was assessed and analyzed with Fiji software (version 2.13.1) and GraphPad Prism 8.0 (GraphPad Software Inc., San Diego, CA, USA).

RNA Isolation and Quantitative Real-Time Polymerase Chain Reaction (qPCR)
Following the manufacturer's instructions, total RNA was extracted from each sample using the TRIzol reagent (Invitrogen Co., Carlsbad, CA, USA). A reverse transcription kit (Promega Corporation, Madison, WI, USA) was used to synthesize the first-strand cDNA from total RNA. The amplification primer sequences for qPCR were summarized in Table 4. The qPCR reaction conditions were as follows: preheating for 10 min at 95 • C, followed by 40 cycles of 30 s at 95 • C, 30 s at 60 • C, 30 s at 72 • C and ending with melting curves. Each sample were conducted in triplicate in individual experiments. The target and reference genes' PCR efficiency was similar, and the intra-assay variation was less than 10%. Using the 2 −∆∆Cq method, the relative expression level of each target mRNA in relation to β-actin was determined.

Hormone Assay
The female muskrats' serum was centrifuged for 15 min at 1000× g at 4 • C. Using ELISA Kit (CSB-E05109m for 17β-estradiol, Cusabio Biotech Co., Ltd., Wuhan, China), the concentrations of 17β-estradiol were determined. Finally, a microplate reader (PT 3502G, Beijing Potenov Technology Co., Ltd., Beijing, China) was used to color the assay plate and read it at 450 nm within 10 min. Each sample was measured in duplicate. The parallelism between the standard curve and the series diluted sample curve was examined in order to establish the validity. The assay's sensitivity was 40 pg/mL. The co-efficients of variation for the intra-and inter-assays were 1.25% and 3.87%, respectively.

Library Preparation for Transcriptome Sequencing
The RNA sample preparations employed a total of 1.5 g RNA per samples (n = 3, each time). The manufacturer's proceeding states that the NEBNext ® UltraTM RNA Library Prep Kit for Illumina ® (NEB, Ipswich, MA, USA) was used to create the RNA-sequencing libraries. Using poly-T Oligo-attached (dT) magnetic beads, the mRNA was extracted and enriched from total RNA. M-MuLV Reverse Transcriptase (RNase II), DNA Polymerase I and random hexamer primer were used to create the double-stranded cDNA from the RNA. The AMPure XP beads were used to purify the double-stranded cDNA, and then PCR amplification was used to create a cDNA library for later library detection and sequencing via synthesis.

Transcriptome Analysis
Adapters and low-quality reads were removed from raw reads to get clean data. Using the previously described method, the measurement of gene expression level was obtained [56]. Using the DESeq R package, the differentially expressed genes (DEGs) of two periods were examined and identified. The genes that were differentially expressed were identified using an adjusted p-value of 0.05. Further implementation of enrichment analyses of DEGs using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) was performed with KOBAS software (version 3.0, accessed on 1 January 2023) [57].

Protein-Ligand Interaction Analysis
Autodock Vina (1.1.2 for Windows, San Diego, CA, USA) was employed to analyze the binding affinities and modes of interaction between the chemicals and vomeronasal receptors [58]. The molecular structures of volatile chemicals and water-soluble chemicals (androsterone sulfate compound CID 159663 and muscone compound CID 10947) were required from PubChem Compound (https://pubchem.ncbi.nlm.nih.gov/, accessed on 1 January 2023). Using Swiss-model software (https://swissmodel.expasy.org/, accessed on 1 January 2023), 7dd6.1.A and 7l0p.1.A were used as models to predict the tertiary structure of V1R and V2R in the female muskrats' VNO, and the sequence with high similarity was selected for the next analysis. After removing all water molecules and adding polar hydrogen atoms, the protein and ligand files were transformed into the pdbqt format. To support free molecular motion and cover the structure of each protein, a cubic grid box with dimensions of 25 × 25 × 34 Å was constructed. Protein-ligand interaction analysis was carried out using AutoDock Vina (1.1.2 for Windows, San Diego, CA, USA), and model visualization and analysis were performed using PyMol (version 2.5.0) [59].

Morphometry and Statistical Analysis
Using the arbitrary distance approach, the thickness of the NSE was measured in micrometers (µm) from the basement membrane to the surface. All measurements were observed under a 10× objective lens with an H&E-stained section, and each measurement was taken 10 times to calculate the relative mean value of NSE thickness. GraphPad Prism 8.0 was used for the statistical analyses after the trial produced quantitative data. To assess group differences, Student's t-test was employed. Statistical values of p < 0.05 were regarded as significant.

Conclusions
To date, this study is the first report to elucidate the morphological and functional alterations of the VNO associated with various reproductive statuses in semi-aquatic mammals. This study highlighted that the female muskrats' VNO was the target organ for estrogen, and estrogen might regulate the female muskrats' VNO response to chemical signals through ERs in coordination with VRs in May and December. Although the reason for the interaction mechanism of ERs and VRs in the VNO of female muskrats remains unclear, the present results suggest the VNO might be the primary target of estrogen in the response to volatile and water-soluble molecules via ERs and VRs to regulate their physiological role in different reproductive statuses. We will further identify the neuronal circuits regulated by natural pheromones coordinating with sex steroids in the brain, which can facilitate our understanding of chemical communication in mammals.