Musk secretion in muskrats ( Ondatra zibethicus L .): association with lipid and cholesterol metabolism-related pathways

: Male muskrats ( Ondatra zibethicus L.) secrete musk from their scent glands during musk secretion season. Musk plays an important role as a communication pheromone during the breeding season. In this study, gas chromatography – mass spectrometry (GC – MS) was used to analyze the main components of musk. The GC – MS results after methyl esteri ﬁ cation showed that 71.55% of the musk is composed of fatty acids. The other components of muskrat musk include cholesterol (9.31%) and other organics. Transcriptome comparison between musk secretion and non-secretion seasons showed signi ﬁ cant changes in the scent glands for 53 genes involved in fatty acid and cholesterol synthesis and metabolism regulatory pathways, which include fatty acid biosynthesis, elongation, and metabolism; steroid biosynthesis; steroid hormone biosynthesis pathways. A reverse transcription – polymerase chain reaction analysis con ﬁ rmed these detected changes. Overall, our results indicated that lipid synthesis and metabolism play important roles in musk compound synthesis by providing energy for musk production, and the produced musk provides a mechanism for male muskrats to communicate with females during the breeding season.


Introduction
The muskrat (Ondatra zibethicus L.) is a rodent that lives in a semi-aquatic environment. Adult male muskrats have a pair of abdominal glands that secrete musk during the musk secretion season (April-May), which is also the breeding season for muskrats. Musk secreted by scented glands plays an important role in communication with females (van Dorp et al., 1973). Male muskrats attract females during the breeding season by secreting this musk. These scent glands experience seasonal cyclical changes (Shereshkov et al., 2006). During the musk secretion season, the glands swell, and secretion fills the glands. This secretion is a milky, viscous fluid that has a relatively strong scent . During the non-secretion season (October-February), the glands decrease in size (Chen et al., 1996a) and do not secrete musk.
Previous chemical analyses showed that muskrat musk contains various components, such as fatty acids, esters, cholesterol, cyclic ketones, and alcohols (Chen et al., 1996b(Chen et al., , 1998. Different studies obtained varying results on the levels of lipids in muskrat musk. In 1994, Bao-Tang et al. (1994 employed gas chromatography-mass spectrometry (GC-MS) techniques to examine muskrat musk after methyl esterification, and found that it contained 12 different types of fatty acids, of which unsaturated fatty acids accounted for more than 60% of the fatty acid content. In 1998, Chen et al. (1998) employed GC-MS to analyze muskrat musk after methyl esterification, and found that it contained five types of fatty acids, which accounted for 50.98% of the musk composition, whereas cholesterol accounted for 7.21% of the musk composition. In 1998, Chen et al. (1998) employed GC-MS techniques and found six types of fatty acids, which accounted for 22.45% of the total content. In 2009, Zhao et al. (2009) employed GC-MS techniques and found that muskrat musk contained five types of fatty acids, which accounted for 29.5% of the musk composition. These findings indicate that fatty acids are the main components of muskrat musk and that cholesterol may be one of the components. However, as there were differences in the sample processing methods and test conditions of these studies, muskrat musk fatty acid quantitation results showed variability. Alternatively, there are very few studies on the cholesterol content of muskrat musk. However, the regulatory mechanism of the synthesis of these components has remained unclear. Therefore, we studied the role of the lipid and steroid synthesis pathways in muskrat secretion using RNA sequencing (RNA-Seq) technology. These findings will provide a basis for understanding how musk is synthesized by muskrats during the musk secretion period.
The aforementioned studies showed that, during the musk secretion season, musk secreted by muskrats produce large amounts of fatty acids, cholesterol, and other substances. However, the synthesis mechanisms and functions of these components are still unclear. In this study, we applied GC-MS to quantify muskrat musk components and analyze lipid components. We combined RNA-Seq techniques to study gene expression changes in the lipid synthesis-and metabolism-related pathways. This study provides evidence that can help understand the synthesis mechanism of muskrat musk components.

Research animals
In this study, all animals were treated in accordance with the National Animal Welfare Legislation, and all experimental procedures were carried out in accordance with the guidelines established by Beijing Forestry University. Six healthy adult male muskrats were purchased from a muskrat breeding center in Xinji City, Hebei Province. All muskrats were euthanized by aeroembolism (inject air into the vein to kill the animals) after anesthesia.

Sample preparations
Sampling of the scent glandular tissues was conducted on 7 November 2016 (during the non-secretion season), and sampling of both scent glandular tissues and secreted musk was conducted on 7 May 2017 (musk secretion season). The collected tissues were photographed and stored in liquid nitrogen for subsequent experimental analyses. The scent gland samples were dehydrated in an ethanol series and embedded in paraffin wax. Serial sections (4-6 μm) were mounted on slides coated with 3-aminopropyltriethoxysilane (APES; ZSGB-BIO, Beijing, China). Some sections were stained with hematite hematoxylin (Solarbio, Beijing, China) for general histological observation.

GC-MS analysis
For the process of methyl esterification, 0.5 g of musk was weighed and dissolved in 5 mL of hexane. The solution was mixed for 10 min to make sure that it is completely dissolved, and 1 mL of 0.5 mol/L NaOH-CH 3 OH solution was added. Heat the mixture to 50°C after swaying for about 30 min to complete the methyl esterification process. To detect the components other than fatty acids, 50 mg of musk was dissolved in benzene and mixed with a whirlpool mixer at 4000 rpm/min for 1 min. The supernatant was placed in 5 mL flasks, from which 2-mL samples were extracted into vials, which were then placed in the GC-MS system to conduct measurements (Shimadzu, Kyoto, Japan). GC-MS system parameters were as follows: DB-1 fused silica capillary column (28.0 m × 250.0 μm; 0.25 μm); precolumn pressure, 50 kPa (splitless injection); inlet temperature, 270°C; interface temperature, 280°C; temperature program: starting point at 100°C for 2 min, followed by increments of 6°C/min up to 280°C, and holding at 280°C for 10 min; carrier gas, He; injection volume, 1 μL; ionizing voltage, electronic ionization; ion source temperature, 250°C; gas flow rate, 1.0 mL/min; electron bombardment energy, 70 eV; electron multiplier voltage, 1.7 kV; and mass scan range, 20.0-400.0 amu. The results of the fragment ion analyses for the mass spectrum of each component were compared with the possible structure obtained from the National Institute of Standards and Technology library. The results were then combined with the retention time and data from previous studies to determine the musk components (Fan et al., 2018).
RNA extraction, library construction, and RNA-seq Total RNA of each scent gland tissue was isolated using Trizol reagent (Qiagen, Valencia, CA, USA), characterized on 1.8% agarose gel, and examined with a NanoDrop ND1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA). The total RNA of each sample was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer's instructions and treated with DNase I (Qiagen, Mississauga, Ontario, Canada). The quality and quantity of total RNA were analyzed using an Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA). The concentrations of total RNA were supposed to be more than 200 ng/μL, and the quantity should be more than 10 μg. RNA integrity was checked by agarose gel electrophoresis (1% agarose gel). Suitable RNA samples were selected for cDNA synthesis.
The RNA insert sizes of these samples were assessed using a Qubit 2.0 (Thermo Fisher Scientific Inc., Carlsbad, CA, USA) and Agilent 2100 Bioanalyzer (Santa Clara, CA, USA). The construction of the libraries and RNA-Seq were performed by the Biomarker Biotechnology Corporation (Beijing, China). mRNA was enriched and purified with Oligo(dT)-rich magnetic beads and then broken into short fragments with fragmentation buffer. Using these cleaved mRNA fragments as templates, first-strand cDNA was synthesized with random hexamers; then, buffer, deoxyribonucleotide triphosphates, RNase H, and DNA polymerase I, were added to synthesize second-strand cDNA. The resulting cDNAs were then purified by AMPure XP beads (Thermo Fisher Scientific Inc., Carlsbad, CA, USA). The purified cDNAs were then subjected to endrepair, and an overhang 'A' base was inserted at the 3' ends of the cDNA fragments. AMPure XP beads were used to select fragment sizes. Next, PCR amplification was performed to enrich the purified cDNA template. Finally, the six libraries were sequenced using an Illumina HiSeq™ 2500, which generated 100-nt paired-end reads.

De novo transcriptome assembly and annotation
After the clean reads were filtered from the raw reads, the reads were then assembled de novo using the Trinity platform (http:// trinityrnaseq.sourceforge.net/). For each library, short reads were first assembled into longer contigs; then, components were generated using the overlap between the contigs. Finally, the components were identified using a De Bruijn graph and read information to obtain 22,680 unigenes. Unigenes were aligned to a series of protein databases using BLASTx (E-value ≤ 10−5), including the NCBI Nr, KEGG (http:// www.genome.jp/kegg/kegg2.html), and GO (http://wego. genomics.org.cn/cgi-bin/wego/index.pl) databases. The deduced amino acid sequences of uni-transcripts were then compared with the Pfam (protein family) database using HMMER (E-value ≤ 10−10) to obtain unigene annotations.

Expression annotation
All usable reads were compared with each unigene using Bowtie (http://bowtie-bio.sourceforge.net/manual.shtml), then the expression level was estimated by RSEM (http:// deweylab.github.io/RSEM/) according to the comparisons. Fragments Per Kilobase of transcript per Million mapped reads (FPKM) values were used to determine the expression abundance of each unigene.

Statistical analysis
Results are presented as means ± standard error of the mean. To investigate any differences between the secretion and nonsecretion seasons, we performed t-tests in Sigma Plot v. 12.5 ; 0.01 < p < 0.05 was considered significantly different; 0.001 < p < 0.01 was considered highly significant. Unigene abundance differences between the samples were calculated using the Benjamini-Hochberg method based on the ratio of the FPKM values, and the false discovery rate (FDR) control method was used to identify the threshold of the p-value. Here, only unigenes with an absolute value of fold change ≥ 2 and FDR < 0.01 were used for subsequent analyses.

Histology
During the musk secretion season, the scent glands are dilated to an average length of 38.3 ± 3.97 mm. In contrast, during the non-secretion season, the glands show an average length of 15.3 ± 6.66 mm. The gland volume is significantly greater during the musk secretion season than during the nonsecretion season (p < 0.01; Fig. 1a). The diameter of musksecreting cells is significantly (p < 0.01) larger during the secretion season (13.93 ± 1.64 μm) than during the nonsecretion season (5.45 ± 0.93 μm; Figs. 1b and 1c).
GC-MS results showed that fatty acids and cholesterol were both major components of muskrat musk, with methyl esterified fatty acids accounting for 71.55% (Tab. 1) and non-methyl esterified cholesterol accounting for 9.31% (Tab. 2) of muskrat musk.

Transcriptome analysis
The Illumina HiSeq 2500 high-throughput sequencing system was adopted to sequence the transcriptome of six samples of muskrat scent glands, which generated 27.84 Gb of clean data. Each sample produced 3.57 Gb of FIGURE 1. (a) An image of the muskrat scent glands during the secretion and non-secretion seasons. The volume of the scent gland is significantly higher in the secretion season than in the non-secretion season (p < 0.01). Left, secretion season; Right, non-secretion season. (b, c) Morphological sections in muskrat scent glands during the secretion and non-secretion seasons. The diameter of the secreting cells is significantly higher in the (b) secretion season than in the (c) non-secretion season (p < 0.01). Abbreviations: GC, glandular cells; IC, interstitial cells; GCa, glandular cavity. Scale bars represent 50 μm. clean data, of which Q30 bases accounted for at least 88.98%. Clean reads were assembled into a total of 14,917,239 contigs. As the length increased, the number of contigs decreased accordingly.
Database comparison provided functional annotation for all unigenes. We compared all the unigenes to the NCBI non-redundant (Nr) database using BLAST with an E-value cutoff ≤ 10−5 and a HMMER parameter E-value ≤ 10−10. As a result, a total of 22680 unigenes with annotated information were obtained. A total of 6013 unigenes were aligned to the COG database, which accounted for 26.51% of all unigenes; 13147 unigenes were aligned to the Gene Ontology (GO) database (57.97% of all unigenes); 9699 unigenes were aligned to the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (42.76% of all unigenes); 13568 unigenes were aligned to the KOG database (59.82% of all unigenes); 16551 unigenes were aligned to the Swiss-Prot database (72.98% of all unigenes), and 22217 unigenes were aligned to the Nr database (97.96% of all unigenes).
The expression changes of genes that regulate the fatty acid metabolic signaling pathways are shown in Tab. 3, the expression changes of genes that regulate the synthesis of cholesterol and related hormone signaling pathways are shown in Tab. 4, and the number of genes that showed significant expression changes in related pathways are shown in Tab. 5. The significant changes in gene expression levels (p < 0.05) during the musk secretion and nonsecretion seasons are shown in Figs. 4 and 5. The changes (p < 0.05) in some of the gene expression levels of the scent glands during the musk secretion and non-secretion seasons are shown in Figs. 6 and 7.

Quantitative real-time PCR results
The RT-PCR results showed that genes with significant changes in expression in the ko00140, ko00900, ko00120, ko00510, ko00020, ko00061, ko00062, ko00071, and ko00561 pathways were detected in the muskrat scent glands. These changes were consistent with the gene expression changes detected in the transcriptome analysis (Figs. S2 and S3).
RNA-Seq results showed that the expression levels of seven genes in pathways related to fatty acid biosynthesis and fatty acid elongation (Fig. S1) underwent significant changes. Of these seven genes, six showed significant increases (Tabs. 3b and 3c; Figs. 4b, 4c, 5b and 5c). Expression levels of five genes related to cholesterol synthesis pathways (Fig. 3) significantly increased (Tab. 4a; Figs. 6a and 7a). The expression levels of many genes in the terpenoid backbone biosynthesis pathway, which is the precursor for cholesterol synthesis (Fig. S1), also significantly increased (Tab. 4e; Figs. 6e and 7e).
qRT-PCR demonstrated that changes in lipid biosynthesis-and metabolism-related gene expression are consistent with the transcriptome analysis results (Figs. S2 and S3), which verified the reliability of the RNA-Seq results.

Discussion
The morphology results showed that, during the musk secretion season, muskrat scent gland size and glandular cell quantity were significantly greater than during the nonsecretion season (Fig. 1), which is consistent with the findings of previous studies (Pushkala and Gupta, 2001;Zhang et al., 2017). Different forms indicate differences in function, and the associated changes in morphology and function require further research. The GC-MS results are consistent with those of Chen et al. (1998). RNA-Seq results showed that the gene expression levels in the pathways related to fatty acid biosynthesis and fatty acid elongation (Fig. S1), cholesterol synthesis pathways (Fig. 3), and the terpenoid backbone biosynthesis pathway (Fig. S1) underwent significant changes. From this, we suggest that these genes and metabolic pathways play important roles in the synthesis of major components of secreted muskrat musk. In addition, other studies have found that the musk fluid or musk secreted by Chinese forest musk deer, civet cats, and other musk-secreting animals also contains large amounts of fatty acids and cholesterol (Mookherjee et al., 2004;Li et al., 2016).
In the pathways that we hypothesized were relevant (Tab. 5), genes in the citrate cycle (TCA cycle), fatty acid biosynthesis, metabolism-related, steroid biosynthesis, and terpenoid backbone biosynthesis pathways significantly changed. These pathways showed a similar trend that actively expressing during the musk secretion season. Fatty acid biosynthesis can be divided into three steps: the de novo synthesis of saturated fatty acids, fatty acid carbon chain elongation, and unsaturated bond formation. During this process, fatty acid carbon chain elongation mainly occurs in the ER and mitochondria. In the ER, carbon-chain elongations are completed with fatty acyl-CoA as the starting point in a four-step cycle of condensation, reduction, dehydration, and re-reduction. Fatty acids are usually stored as triglycerides in organisms (Tvrdik et al., 2000;Leonard et al., 2004;Jakobsson et al., 2006). RNA-Seq results showed that, in the glycerolipid (Fig. S1) metabolic pathway, fatty acid degradation (Fig. 2), and citrate cycle (Fig. S1), 23 genes showed significant changes in gene expression, of which 19 were significantly increased (Tabs. 3a, 3d and 3e; Figs. 4a, 4d, 4e, 5a, 5d, and 5e). Triglycerides are important energy storage material, and fatty acid degradation can release large amounts of energy (Cheung and Wong, 2008). Triglycerides can undergo hydrolysis to generate fatty acids and glycerol. Fatty acid oxidation is a catabolic reaction that produces CO 2 and H 2 O and can provide the organism with a large amount of energy. In animals, fatty acid catabolism occurs in mitochondria and peroxisomes (Zhang et al., 2002). Fatty acids are converted to fatty acyl-CoA for β-oxidation. There are four β-oxidation steps: oxidation, hydration, FIGURE 3. The KEGG pathway of steroid biosynthesis. The red marks show genes that are up-regulated. This KEGG pathway map is adapted here from http://www.kegg.jp/kegg/kegg1.html. The KEGG database has been described previously (Kanehisa and Goto, 2000;Kanehisa et al., 2017Kanehisa et al., , 2019 re-oxidation, and thiolysis. These four cyclic steps continue until all fatty acyl-CoA have become acetyl-CoA, which goes through the citrate cycle to produce CO 2 and a large amount of energy (Cheung and Wong, 2008). Musk secretion is a physiological activity that consumes large amounts of energy. During the musk secretion season, the provision of abundant fodder to muskrats can significantly increase the amount of musk secreted by muskrats (Chen, 1988;Huang, 1998). Studies have found that adiponectin plays an important regulatory role during musk secretion in muskrats  and causes a significant increase in the activity of intracellular mitochondria in muskrat scent gland cells, during the musk secretion season (Chen et al., 1996b;Chen, 2007). Therefore, we propose that fatty acid oxidation and related metabolic pathways produce large amounts of energy that are required for muskrat musk secretion.
Cholesterol is a small lipid molecule that serves an important physiological role in organisms. It is a precursor for the synthesis of many bioactive molecules, such as bile acid and steroid hormones (Hanukoglu, 1992;Liu and Song, 2013). Cholesterol synthesis is mostly based on the acetyl-CoA in the cytosol, which, under enzyme catalysis, sequentially becomes terpenoids, such as mevalonate, isoprene units, squalene, and lanosterol. After the removal of three carbons, cholesterol is formed (Zhang et al., 2002). Steroid hormones include hormones from the adrenal cortex and sex hormones, which are secreted from the adrenal cortex and gonads (Pushkala and Gupta, 2001). Different enzymes can then stimulate cholesterol to produce different bioactive molecules. Bile acids are the major end-product of cholesterol metabolism in the body and can assist in lipid absorption and transport, and promote fatty acid oxidation (Davis, 2007;Russell, 2009). As noted above, fatty acid metabolism provides large amounts of energy for musk secretion in muskrats. As an important regulatory factor in fatty acid metabolism, primary bile acids play an important role in musk secretion. Transcriptome results showed that the expression levels of five genes in the primary bile acid biosynthesis pathway (Fig. S1) showed significant increases (Tab. 4b; Figs. 6b and 7b). These findings showed that this pathway is enhanced during the musk  The amount of genes that showed significant changes in expression between the secretion and non-secretion seasons in related pathways secretion season compared with during the non-secretion season, and the increase in primary bile acid synthesis supported this hypothesis. A study showed that steroid hormones are an important factor for scent gland development and secreting musk (Li et al., 2011). Previous experimental results have revealed that muskrat scent gland cells have the ability to self-synthesize male hormones, and the autocrine and paracrine secretions of sex hormones jointly regulate musk secretion (Han et al., 2017). The transcriptomics results showed that the expression levels of six genes in the steroid hormone biosynthesis pathway (Fig. S1) showed significant changes, of which four were significantly increased (Tab. 4c; Figs. 6 and 7c). These results are consistent with the hypothesis that muskrat scent gland cells can carry out the autocrine secretion of sex hormones, and the expression levels of many crucial genes involved in steroid hormone synthesis are significantly increased during musk secretion season. Our transcriptomic results provide an important basis for studying the autocrine steroid hormones secreted by muskrat scent glands that regulate musk secretion.

Name of pathways
In animals, terpenoids can be used for N-glycan synthesis in addition to cholesterol synthesis (Cheung and Wong, 2008). However, there are seven genes in the N-glycan biosynthesis pathway (Fig. S1) that showed significant changes, of which three were significantly increased (Tab. 4d; Figs. 6d and 7d). These results showed that the N-glycan synthesis pathway in muskrats, compared with the cholesterol synthesis pathway, may not play a substantial role during the musk secretion season.
Recent studies have found that the levels of cytokines such as adiponectin, which regulate lipid metabolism, are significantly increased during the musk secretion season, and jointly participate in regulating musk secretion. This research applies GC-MS and RNA-Seq techniques to study the functions of lipid biosynthesis and metabolism regulation in muskrat musk secretion. Gene expression in lipid and steroid biosynthesis pathways provide the main components of musk during musk secretion season. Musk is an important pheromone that males use to attract and communicate with females. During the mating season, musk production and secretion requires a lot of energy. Gene expression in the lipid metabolic FIGURE 6. Transcript accumulation measurements of cholesterol biosynthesis-related genes involved in the muskrat musk secretion process. The ID of each enzyme name and expression pattern is indicated above. The expression pattern of each uni-transcript is shown on two grids, with the left one representing the FPKM values in the secretion seasons, and the right one representing the FPKM value in the non-secretion season. The grids with sixteen different grayscale levels show the absolute expression magnitude of secretion seasons, with the RPKM values 0-10, 10-20, 20-30, 30-40, 40-50, 50-60, 60-70, 70-80, 80-90, 90-100, 100-300, 300-500, 500-700, 700-900, 900-1100, and over 1100 represented by grayscale levels 1-16, respectively (FDR < 0.01, FC ≥ 2). (a) Steroid biosynthesis, (b) primary bile acid biosynthesis (c), steroid hormone biosynthesis, (d) N-glycan biosynthesis, (e) terpenoid backbone biosynthesis.
pathway provides the energy to support this important physiological process. Musk components and secretion regulation are complex issues, and there is still a substantial amount that is unknown about this unique physiological phenomenon.

Acknowledgement:
The experimental facilities were provided by the School of Ecology and Nature Conservation Innovation Laboratory, Beijing Forestry University. We thank Mallory Eckstut, Ph.D., from Liwen Bianji, Edanz Editing China (www.liwenbianji.cn/ac), for editing the English text of a draft of this manuscript.

Data Availability Statement:
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest:
The authors declare that they have no conflicts of interest to report regarding the present study. FIGURE S1. KEGG reference mappings for citrate cycle (TCA cycle), fatty acid biosynthesis, fatty acid elongation, glycerolipid metabolism, primary bile acid biosynthesis, steroid hormone biosynthesis, N-glycan biosynthesis and terpenoid backbone biosynthesis pathways. The red marks shows genes that are up regulated and the green ones shows down regulated. The blue one means some of the genes are up regulated and some are down regulated. FIGURE S2. Grey scale measurements of fatty acid metabolism-related genes and GAPDH in the scent gland by RT-RCR. Data are presented as the mean ± SEM from at least three independent experiments, and error bars indicate SEM; Ã , ÃÃ and ÃÃÃ Denote statistically significant differences between secretion and non-secretion seasons ( Ã , 0.01 < p < 0.05; ÃÃ , 0.001 < p < 0.01, ÃÃÃ , p < 0.001).