Sample collection and sample size

Nutrias were collected through the culling efforts of local authorities (Jewish National Fund–JNF–KKL) by professional hunters using shotguns at the Hula Lake Park, Israel over several campaigns from December 2019 to March 2021. We followed the hunters, collected the culled nutrias, and dissected them within 30 min. Since nutria are an invasive species in Israel and are not protected by law, permits were not required for collecting carcasses. In this study, a total of 27 pregnant females (with a gestation age of 60–131 days) and their fetuses (N = 153) were included. Carcasses of pregnant females were weighed using a spring scale (10 kg, Pesola, Switzerland) and morphometric measurements (e.g., total length, nose to the base of the tail, shoulders to the base of the tail, tail length, tail base circumference, and right hind foot) were conducted using a measuring tape to the nearest millimeter. Upon dissection, the IUP for each fetus was determined, including the uterus horn location (left or right), and its location relative to the ovary (the closest fetus to the ovary was numbered one). Moreover, two conventions were used to label fetuses according to their proximity to other fetuses (Fig. 6). Fetuses were sexed according to external morphology using AGD (see Fig. 1, [19]). Fetuses in litters with an estimated gestational age of > 60 days were weighed using an analytical balance to the nearest 0.01 mg (Precisa, Switzerland). The gestational age was estimated using Newson’s formula [55]: Estimated gestational age = \(43.69+14.27\times \sqrt[3]\)

Fig. 6

IUP labeling example using two conventions. The Contiguity Method (vom Saal et al. [76]) refers to a fetus’ position relative to males. 0 M = fetus not next to a male; 1 M = fetus neighboring one male; 2 M = fetus neighboring two males. The second convention refers to a fetus’ position relative to a fetus of the opposite sex [19]. P0 = fetus next to fetuses of their own sex, P1 = fetus neighboring a fetus of the opposite sex on either or both sides. Figure designed by Ariel Yael, drawn by Aiden Braner

Tissue isolation and storage

Whole brains were isolated from 64 fetuses (11 litters, with a gestational age of > 85 days) and grossly dissected into the hypothalamus, prefrontal cortex, and striatum (Fig. 7). Samples were placed in sterile cryogenic tubes, snap-frozen with liquid nitrogen for several hours, and stored at − 80 ºC until analysis.

Fig. 7

Fetal nutria brain dissections. A Nutria whole fetal brain, ventral view; B Prefrontal cortex was dissected with a scalpel, at an approximately 45º angle; C Hypothalamus, located on the ventral side of the brain, was removed with sharp tweezers (D once removed); Not pictured: the striatum

RNA extraction and RNA quality analysis

Frozen brain samples (tissue weight: 22.56–49.99 mg) were thawed on ice and homogenized for 30 s using TissueLyser II (Qiagen, Germany) with two sterile 2.5 mm stainless-steel beads in 500 µL of TRIzol reagent (Bio-Lab Ltd., Israel). Following homogenization, total RNA was extracted with TRIzol according to the manufacturer’s protocol. RNA yield and purity, analyzed by the OD A260/A280 and A260/A230 ratio, were determined by the NanoDrop 2000 Spectrophotometer (Thermo Scientific).

Sequence analysis for primer design

Like the analysis detailed by Matas et al. [47], assembled nutria genome sequences were downloaded from the NCBI site (https://www.ncbi.nlm.nih.gov/assembly/GCA_004027025.1/) and used to create a BLAST database. Using tBLASTx, Rattus norvegicus HPA (GR and MR) and HPG (GnRH, GNRHR, AR, aromatase, and ESR1) axes genes were used as a query to identify sequences of corresponding nutria orthologs (SI Table 1C). The best hit both in terms of coverage and similarity was selected. Seqret from the EMBOSS package was used to extract the corresponding genomic loci. The mRNA sequence was predicted using GENSCAN, then Splign was run using the predicted mRNA sequence and the corresponding genomic loci to ascertain the exact exon–intron splice sites. Identified coding sequences for HPA and HPG axes components are documented in SI Table 1B. The housekeeping gene glycogen synthase 1 (Gys1) was also identified (SI Table 15). Identified protein sequences of these nutria genes were analyzed by multiple sequence alignment (MSA) with other homologs using CLUSTAL Omega algorithm. Similarity and coverage were determined by NCBI blast. Identification of known domains was carried out by NCBI Search for conserved domains within a protein.

Reverse transcription and real-time PCR

Genomic DNA was removed from nutria fetal brain RNA samples by DNA-free DNA Removal kit (Invitrogen, Thermo Fisher Scientific, MA, USA), and reverse transcription of the purified RNA was performed using the qScript cDNA Synthesis Kit (Quantabio, MA, USA). The resulting cDNA was used as a template for real-time PCR (qPCR) with CFX- 96 (Bio-Rad Laboratories, Hercules, CA, USA) using SYBR Green FastMix/ROX qPCR Master Mix (Quantabio, MA, USA), according to the manufacturer’s instructions. We measured mRNA expression levels of the genes for GR, MR, GnRH, GNRHR, AR, aromatase, and ESR1 in the prefrontal cortex, hypothalamus, and striatum. We also quantified the expression of the housekeeping gene Gys1 for normalization of target genes. Forward and reverse primer sequences are reported in SI Table 16. Primers were synthesized by Syntezza Bioscience Ltd. (Jerusalem, Israel) from exon–intron junction sequences. All primers were validated for non-specific amplification and primer-dimer formation by melt curve analysis. Efficiencies of all primers range from 90.3 to 100.5%. qPCR reactions were run in triplicates using cDNA diluted according to the validations. Thermocycling conditions were: 95 ℃ for 30 s and 55 ℃ for 30 s, for 40 cycles. A final melting phase of 65–95 ℃ (by increments of 0.5 ℃) for 30 s was run to confirm the single-product specificity of each sample.

Hair testosterone and cortisol extraction and quantification

Using hair to quantify steroids provides a long-term steroid accumulation profile over the time of growth and reflects free, unbound, and presumably bioavailable hormones [37]. Moreover, hair steroid levels are not sensitive to the acute stress of culling [19, 37]. This method has been validated in this system as an indicator of maternal status during pregnancy and of the condition of the fetuses in utero [19–21]. Fetal hair follicles in the nutria appear at 85–90 days of gestation, around the beginning of the last trimester of their ~ 130 day-long gestation [21]. Hair samples were collected from 25 fetuses past the gestational age of 111 days. As detailed in Fishman et al. [21], all fetuses were washed under the tap with deionized water to remove blood and amniotic fluid, then dried for at least 24 h. Once dried, hair was shaved from the entire fetus using an electric razor (Moser Chromini 1591 0.7 mm). Hair samples were washed twice with isopropanol (Romical, Israel) on a shaker for three minutes to remove external contaminants and dried again for a minimum of 24 h in an open petri dish (10 mm). Dried hair samples were weighed, and steroids were extracted by methanol (Romical, Israel), sonicated for 30 min in a sonication bath (Elma, Germany), and then incubated overnight at 50 ℃ for 19 h with gentle shaking. Methanol was then separated from hair, centrifuged to remove small particles, and dried under a stream of nitrogen gas in a Dry-block heater at 45 ℃ (Techne, UK).

Hair cortisol and testosterone were quantified with a commercial enzyme-linked immunosorbent assay (ELISA) kit (Salimetrics, Europe, Newmarket, UK) using a protocol that was developed for wildlife and previously validated for nutria [19–21]. For testosterone, the manufacturer reported antibody cross-reactivity of 36.4% with dihydrotestosterone, 21.02% with 19-nortestosterone, 1.9% with 11-hydroxytestosterone, 1.157% with androstenedione, and < 0.49% with all other steroids. For cortisol, reported antibody cross-reactivity was 19.2% with dexamethasone and < 0.568% with all other steroids. Validations showed linearity (in the range of 5–40 mg hair for testosterone and 0.5–10 mg hair for cortisol) and parallelism between serially diluted hair extracts and kit standards (slope covariance P = 0.91 for testosterone and P = 0.36 for cortisol). Intra-assay CV for six repeats on the same plate was 7.23% for testosterone and 5.29% for cortisol. Inter-assay CV was 0.6% for testosterone across four plates and 13.1% for cortisol across five plates. Recovery was quantified by spiking hair samples with a known steroid amount and was calculated as 100.7% for testosterone and 90.9% for cortisol.

LC–MS/MS validation

The presence of cortisol and testosterone in a pooled nutria hair extract was confirmed via LC–MS/MS using a similar sample preparation method as a recent study of cortisol in muskox qiviut [15]. The hair pool was gently washed by hand rotation in 20 mL ice-cold tap water with Neutrogena TM water for 3 min and then rinsed by tap water for 0.5 min. After being quickly rinsed by 20 mL ice-cooled HPLC grade IPA, the sample was dried by paper towel and placed in a fume hood overnight. The sample was then transferred into a 13 × 100 mm culture test tube. As an additional validation, another test tube contained 200 mg of the burial matrix that surrounded the individual hair samples. To both test tubes, 100 μL of deuterium labeled internal standard solution and 9 mL cold methanol were added. The test tubes were capped and stored in a 4 ℃ fridge for 20 h. After the sample was removed from the test tube, the extract was evaporated to dryness under N2 at 40 ℃ by use of Techne Sample Concentrator and reconstituted with 150 µL of H2O/MeOH (50/50, v/v). The solution was centrifuged at 14,000 rpm (Legend micro-21R, Thermo Scientific) for 20 min and 120 µL of supernatant was submitted to LC–MS.

All samples were analyzed by using an Agilent 1200 binary liquid chromatography (LC) system connected with an AB SCIEX QTRAP® 5500 tandem mass spectrometer equipped with an atmospheric pressure chemical ionization (APCI) source. LC separation was performed on an Agilent Poroshell 120 C18 column (50 × 3 mm, 2.7 µm particle size) at 45 ℃. The mobile phase A was H2O/MeOH (75/25, v/v) and the mobile phase B was MeOH/IPA (90/10, v/v). The 8.5 min gradient was 20–40% B (0–1.0 min), 40–60% B (1.0–5.0 min), 60–100% B (5.0–5.5 min), 100% B (5.5–6.5 min), 100–20% B (6.5–7.0 min), and held at 20% B (7.0–8.5 min). The flow rate was 0.6 mL/min and the injection volume was 20 µL. The analytes were ionized under positive APCI mode and data were acquired via multiple reaction monitoring (MRM). The lowest limit of quantification for cortisol was 0.25 ng/mL and for testosterone 0.05 ng/mL, which were the lowest concentrations that gave < 20% CV and <  ± 20% error. Nutria hair pool steroid levels showed 20.26 pg cortisol/mg hair and 0.51 pg testosterone/mg hair. More details about the quantification protocol are provided in the SI (Tables 17–19).

Amniotic fluid (AF) testosterone and cortisol extraction and quantification

AF was sampled from 52 fetuses (8 litters, with gestational age > 85 days) whose amniotic sac was intact. AF was extracted from the amniotic sacs of each individual fetus using a sterile syringe and needle. AF samples were flash frozen in liquid nitrogen and stored in -80ºC until analysis. Cortisol was quantified in AF samples following dilution in assay diluent, which was provided with the commercial ELISA kit (see below). For testosterone measurement, AF samples were extracted with ethyl acetate, vortexed, and centrifuged at 13,300 g for 10 min. The supernatant was collected and evaporated under a stream of nitrogen. Samples were then reconstituted in 10% methanol and 90% assay diluent.

Testosterone and cortisol were quantitated in AF extracts using commercial ELISA, according to the manufacturer’s recommendations (Salimetrics; Ann Arbor, MI, USA, item no. 1-3002-5, for cortisol, and item no. 1-2402-5 for testosterone). The cross-reactivity of the kits is reported above, in the quantitation in hair section. Testosterone validation using serial dilutions of AF pool (N = 20) showed linearity between 10 and 200 µL AF (equivalent to 15–240 pg/ml testosterone) and parallelism with the provided kit standards (univariate analysis of variance in SPSS; P = 0.805). Intra-assay repeatability was determined using 6 replicates of the pool on the same ELISA plate. The calculated coefficient of variation was 12.6%. Inter assay variability was determined by running duplicates of the pool on 4 different plates (n = 8). The coefficient of variation was 18.32%. Efficiency of 116.44% was retrieved using exogenous testosterone. Serial dilutions of AF pool (n = 12 samples) for cortisol validation showed linearity between 2 and 22 µL AF (equivalent to 0.111–3 μg/dL cortisol). Dilutions of the pool were parallel to the kit standards (univariate analysis of variance in SPSS; P = 0.071). Intra-assay repeatability using 6 replicates of the pool on the same ELISA plate showed coefficient of variation of 6.546%. Inter assay variability was determined by running 6 duplicates of the pool (n = 12). The coefficient of variation was 7.11%. Recovery of 84.5% was demonstrated using exogenous cortisol.

Statistical analyses

All statistical analyses were performed using the JMP software (v.16). We used general linear mixed models (GLMM) with standard least squares to test whether sex (M/F), IUP (P0/P1), and their interaction (sex*IUP) affect HPA and HPG axes components separately for each brain region. We tested normality with Shapiro-Wilks tests, and when necessary, we used a log transformation (hypothalamus: GnRH, AR, ESR1, GR, and MR; cortex: GnRH, GNRHR, AR, aromatase, GR, and MR; striatum: GNRHR, AR, aromatase, ESR1, and MR). When the log-transformed values were still not normally distributed, the SHASH transformation (hypothalamus: GNRHR; cortex: ESR1; striatum: GnRH, MR) or the Johnson Sb transformation (hypothalamus: aromatase) were used.

Maternal identity and year were included as random factors in all models to account for variability between individual mothers and sampling over the 2 years (2020 and 2021). Litter size and season were included in all models. Estimated gestational age and fetal weight are interchangeable variables since the calculation of estimated gestational age depends on the average fetal weight for each litter [55]. For HPG axis components, we tested models including either estimated gestational age, fetal weight, or a residual of estimated gestational age and fetal weight, along with litter size and season, using model selection (AICc value; Tables SI 20–22). For HPA axis genes, we tested models including either estimated gestational age, fetal weight, a residual of estimated gestational age and fetal weight, or shoulders to the base of the tail length (SBL), or a residual of SBL and estimated gestational age, along with litter size and season, using model selection (AICc value; Tables SI 23–25). Fishman et al. [20] observed that nutria fetuses neighboring an opposite-sex fetus in utero were longer from shoulder to base of the tail, regardless of sex, implicating better lung development. Therefore, for analyses of HPA axis components (i.e., GR and MR mRNA expression level), we included SBL in the model criteria test.

In a subset of fetuses in the last trimester of gestation (i.e., with hair samples), we tested whether HPA axis components expression is affected by sex, IUP, sex*IUP, and fetal hair cortisol levels and whether HPG axis components expression is affected by sex, IUP, sex*IUP, and fetal hair T levels. To test the bidirectional interaction between the HPA and HPG axes, we analyzed the effect of sex, IUP, sex*IUP, and fetal hair T levels on HPA axis components and sex, IUP, sex*IUP, and fetal hair cortisol levels on HPG axis components.

To examine the relationship between gene expression across brain regions, we analyzed pairwise correlations using Pearson’s correlation coefficient (R). Correlation tables for all genes in all brain regions were split by age according to fetal weight and SBL averages (Group 1: 85–110 days; Group 2: 111–130 days). The split was verified using Tukey–Kramer analysis that showed significant differences between developmental groups.