A microRNA or messenger RNA point of departure estimates an apical endpoint point of departure in a rat developmental toxicity model.

Traditional developmental toxicity testing practice examines fetal apical endpoints to identify a point of departure (POD) for risk assessment. A potential new testing paradigm involves deriving a POD from a comprehensive analysis of molecular-level change. Here, the rat ketoconazole endocrine-mediated developmental toxicity model was used to test the hypothesis that maternal epigenomic (miRNA) and transcriptomic (mRNA) PODs are similar to fetal apical endpoint PODs. Sprague-Dawley rats were exposed from gestation day (GD) 6-21 to 0, 0.063, 0.2, 0.63, 2, 6.3, 20, or 40 mg/kg/day ketoconazole. Dam systemic, liver, and placenta PODs, along with GD 21 fetal resorption, body weight, and skeletal apical PODs were derived using BMDS software. GD 21 dam liver and placenta miRNA and mRNA PODs were obtained using three methods: a novel individual molecule POD accumulation method, a first mode method, and a gene set method. Dam apical POD values ranged from 2.0 to 38.6 mg/kg/day; the lowest value was for placenta histopathology. Fetal apical POD values were 10.9-20.3 mg/kg/day; the lowest value was for fetal resorption. Dam liver miRNA and mRNA POD values were 0.34-0.69 mg/kg/day, and placenta miRNA and mRNA POD values were 2.53-6.83 mg/kg/day. Epigenomic and transcriptomic POD values were similar across liver and placenta. Deriving a molecular POD from dam liver or placenta was protective of a fetal apical POD. These data support the conclusion that a molecular POD can be used to estimate, or be protective of, a developmental toxicity apical POD.

INTRODUCTION

Conventionally, human health risk assessment of crop protection molecules has relied upon animal toxicity studies to identify the lowest point of departure (POD) among all adverse effects observed. Of several toxicity studies conducted to support crop protection molecule registration, the guideline‐driven developmental toxicity study examines effects on fetal endpoints just prior to parturition following exposure during embryofetal development (OECD, 2018). The conceptus endpoints analyzed are termed apical endpoints and include observations such as organ gross morphology, skeletal ossification/morphology, fetal body weight, and post implantation loss (embryo/fetal death) (Johnson et al., 2016). While guideline‐driven developmental toxicity studies provide a detailed analysis of potential adverse effects at a high level of biological organization, this study design is time consuming, uses a large number of animals/study (>1,000), and has been criticized for lacking sensitivity and a comprehensive analysis of potential effects (Tweedale, 2017).

A modernized safety assessment study design that leverages genome‐wide molecular profiling to determine a benchmark dose‐based transcriptome POD has been suggested as an alternative to conventional apical toxicity testing paradigms (Buesen et al., 2017; LaRocca, Johnson, LeBaron, & Rasoulpour, 2017; Mezencev & Subramaniam, 2019; Schmitz‐Spanke, 2019). In this new assessment method, a molecular‐level benchmark dose (BMD) is identified by profiling the transcriptome via a comprehensive method such as RNA sequencing (RNAseq) using BMD software designed for transcriptome data (BMDExpress) (Phillips et al., 2019). Using this process, a limited number of studies have examined the concordance between transcriptome and apical POD values. Target organ transcriptome POD values typically are within an order of magnitude of apical endpoint POD values for rodent carcinogenicity or subchronic general toxicity study designs (Bianchi et al., 2021; Chepelev et al., 2018; Gwinn et al., 2020; Jackson et al., 2014; Moffat et al., 2015; Thomas et al., 2011; Thomas et al., 2013). To date, a comparison of transcriptome and apical POD values for a developmental toxicity phenotype has not been performed.

Pregnant rat oral exposure to high dose levels of some molecules within the azole class of chemistries causes a suite of adverse apical effects in the fetus. Liver toxicity is also commonly observed in adult rats following azole exposure, and liver apical effects observed include hepatocyte hypertrophy, vacuolation, and necrosis (Heise et al., 2015; Khoza, Moyo, & Ncube, 2017). Ketoconazole, a pharmaceutical compound, represents a model developmental and liver toxicant in the azole class. While human clinical exposure to ketoconazole is not associated with developmental toxicity (Kazy, Puhó, & Czeizel, 2005), adverse fetal effects of rat ketoconazole exposure include cleft palate, decreased body weight, and fetal death (i.e., post‐implantation loss) (Amaral & Nunes, 2008; Nishikawa, Hara, Miyazaki, & Ohguro, 1984; Taxvig et al., 2008). The molecular initiating event for cleft palate and fetal death is direct inhibition of cytochrome p450 (CYP) enzymes (Marotta & Tiboni, 2010). For cleft palate, the azole toxicity mode‐of‐action (MoA) includes CYP26 inhibition and increased retinoic acid signaling in embryonic neural crest cells leading to altered craniofacial skeletal development (Marotta & Tiboni, 2010; Menegola, Broccia, Di Renzo, & Giavini, 2003; Tiboni, Marotta, & Carletti, 2009). For fetal death, the azole toxicity mode‐of‐action (MoA) includes the following precursor key events (Menegola et al., 2003; Stinchcombe et al., 2013): 1) aromatase (CYP19A1) inhibition in the maternal ovary (Stinchcombe et al., 2013); 2) reduction of circulating maternal estradiol levels during late gestation (Stinchcombe et al., 2013; Taxvig et al., 2008); and 3) placenta functional deficit observed as an increase in placenta weight and histopathology (Furukawa, Hayashi, Usuda, Abe, & Ogawa, 2008; Ichikawa & Tamada, 2016).

To examine the concordance of maternal liver and placenta epigenomic (miRNA) and transcriptomic (mRNA) PODs with maternal and fetal apical endpoint PODs, the rat ketoconazole developmental toxicity model was leveraged. The maternal apical endpoints examined were feed consumption, body weight, and liver and placenta weight and histology. The fetal apical endpoints examined were body weight, fetal death, and skeletal ossification and morphology. Using the ketoconazole model, it was hypothesized that 1) maternal epigenomic (miRNA) and transcriptome (mRNA) POD values are similar and predict the maternal apical POD and 2) the maternal liver and placenta miRNA and mRNA POD estimates the fetal apical POD within 10X.

METHODS

Exposure paradigm

The study design was chosen to model the OECD 414 developmental toxicity guideline study design (OECD, 2018). Ketoconazole (CAS number 65277‐42‐1) was purchased from Alfa Aesar (Tewksbury, MA, USA) (lot number Y02B008) with a manufacturer certificate of analysis purity of 99.5% by ultra‐high performance liquid chromatography. On each morning from GD 6–20, time‐mated female rats (n = 10/exposure group) were administered ketoconazole in corn oil vehicle (Sigma‐Aldrich; St. Louis, MO, USA; catalog number C8267) via oral gavage. Ketoconazole dose levels were 0 (corn oil vehicle alone), 0.063, 0.2, 0.63, 2, 6.3, 20, or 40 mg/kg(body weight)/day (mkd) at a dose volume of 2 mL/kg body weight.

Dose levels of 20 and 40 mkd were selected to produce treatment‐related fetal apical effects in at least two ketoconazole dose levels (Amaral & Nunes, 2008; Mineshima et al., 2012) which were expected to provide adequate effect size data for BMD analysis (Davis, Gift, & Zhao, 2011; Slob, 2014). Lower dose levels and dose spacing were selected to generate data near the apical and molecular POD values and provide data in the lower dose range for robust BMD analysis (Davis et al., 2011; Slob, 2014).

Animal care and use

Animal care and use were performed at the Toxicology and Environmental Research and Consulting laboratory of The Dow Chemical Company in Midland, MI. This test facility was fully accredited by the Association for the Assessment and Accreditation of Laboratory Animal Care International (AAALAC International). The animal experimentation protocol was approved by the Institutional Animal Care and Use Committee.

Sexually mature, virgin female Crl:CD(SD) rats weighing approximately 200–250 g were mated at Charles River Laboratories (Raleigh, NC). The morning of observing a vaginal plug was considered gestation day (GD) 0. Time pregnant rats arrived at the test facility on GD 1 or 2. Prior to study start, all animals were evaluated by a trained veterinarian who verified the appropriate health status of each animal on study. Animals were stratified by GD 0 body weight and then randomly assigned to treatment groups using a computer program designed to increase the probability of uniform group weights and standard deviations at the start of the study.

The vivarium room was maintained at approximately 22°C with approximately 50% humidity and a 12‐h light/dark photoperiod. Animals were housed one per cage in solid bottom stainless steel cages with corn cob bedding. Animals were provided water via a pressure activated lixit valve‐type system and meal‐form feed (LabDiet Certified Rodent Diet #5002; PMI Nutrition International, St. Louis, MO) ad libitum. Enrichment included a nylon rodent chew (Animal Specialties and Provisions, Quakertown, PA; catalog number WGP100) and open areas on the cage side for visualization of other rats.

Clinical observations were conducted on all animals at least once daily. Body weights were recorded on GD 0 and daily from GD 6–21. Feed consumption was recorded for all animals every 3 days from GD 3–21. On the morning of GD 21 approximately 24 h following the last ketoconazole dosing, non‐fasted dams were anesthetized with a mixture of isoflurane vapors and medical‐grade oxygen for blood collection via the orbital sinus, further anesthetized with carbon dioxide, and euthanized by cervical dislocation and exsanguination.

Necropsy

All data collection during and subsequent to necropsy was performed with the observer blinded to treatment group. Dam liver, gravid uterus, placentae (juxtaposed to viable fetuses only), and fetal weights (viable fetuses only) were recorded. The number of implantations, dead or resorbed fetuses, and fetuses with cleft palate was recorded. The uteri of females lacking visible implantations were stained with sodium sulfide to verify pregnancy status (Kopf, Lorenz, & Salewski, 1964). Decapitated fetuses were skinned, eviscerated, preserved in alcohol, and double stained with Alcian Blue and Alizarin Red S for skeletal examination (Trueman, Jackson, & Trueman, 1999). A random number generator was used to select one placentae from each uterine horn, which were cut into quarters, pooled, and placed in RNAlater (ThermoFisher Scientific, Waltham, MA, USA). All remaining placentae adjacent to a live fetus were fixed in neutral buffered 10% formalin and processed into paraffin blocks. The upper third of the left lateral lobe of the dam liver was placed in RNAlater; the remaining dam liver was fixed in neutral buffered 10% formalin and processed into paraffin blocks. Dam blood was centrifuged at 4°C and 1200 g for 20 min, and the resulting plasma stored in salinized glass vials at −20°C for estradiol (E2) quantification.

Maternal plasma estradiol measurement

Plasma samples from non‐pregnant females or dams with completely resorbed litters were not analyzed. E2 was quantified via electrochemiluminescence immunoassay (Estradiol III kit; catalog number 06656021 190; Roche Diagnostics, Indianapolis, IN, USA) according to the manufacturer's protocol using a cobas e411 Immunoanalyzer (Roche Diagnostics, Indianapolis, IN, USA). The measuring range was 5–3000 pg/mL. Within treatment group coefficients of variation ranged from 13.01–64.54%.

Histology

Six micron thick placenta and liver sections were stained with hematoxylin and eosin and examined by a board‐certified veterinary pathologist. Histopathologic findings were graded to reflect the severity of the specific lesions. Very slight and slight grades reflected lesions of minimal severity and typically with <25% parenchyma involvement. A moderate grade was used when the severity extent (up to 50% of the parenchyma) may adversely affect organ function. A severe grade reflected a finding extensive enough to cause significant organ dysfunction or failure.

RNAseq and identification of differentially expressed microRNAs and mRNAs

Samples from non‐pregnant females or dams with completely resorbed litters were excluded from the analysis. Using a random number generator, placenta and dam liver from the same five maternal‐fetal units/group were selected for RNAseq. Placenta were minced finely with a razor blade, and a random selection of pieces was used for total RNA isolation. For gene expression analysis, total RNA (small and large) was isolated from 30 mg of placenta and dam liver using a mirVana RNA isolation kit (catalog number AM1560; ThermoFisher Scientific, Waltham, MA, USA). Total RNA integrity was determined using an Agilent Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA).

For mRNA, a sequencing library was prepared from 1 μg total RNA using the Illumina TruSeq Stranded mRNA Library Prep Kit and protocol (Illumina, San Diego, CA, USA). Briefly, mRNAs were fragmented using ZnCl2 solution, reverse transcribed into cDNA using random primers, end repaired to create blunt end fragments, 3′ A‐tailed, and ligated with Illumina paired‐end library adaptors. Ligated cDNA fragments were PCR amplified using Illumina paired‐end library primers and purified PCR products were checked for quality and quantity on the Agilent TapeStation (Agilent, Santa Clara, CA, USA). Average library insert size was 300 bp, and final sequencing concentration was determined by quantitation on the Agilent TapeStation (Agilent, Santa Clara, CA, USA).

For miRNA, miRNA was purified from the same tissue homogenate as the total RNA using a mirVana kit (catalog number AM1560; ThermoFisher Scientific, Waltham, MA, USA) and quantified using a Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Sequencing libraries from each small RNA sample were prepared using 1 μg input with the TruSeq small RNA kit according to manufacturer's instructions (Illumina, Inc., San Diego, CA, USA). Briefly, Illumina TruSeq adapters are ligated onto purified smRNA ends, reverse transcribed, and the cDNA are PCR‐amplified using Illumina indexed TruSeq primers. PCR‐amplified libraries were checked for quality and quantified using the Agilent Tapestation 4200 with D1000 screentape. Equal amounts of library were pooled based on nano‐molarity concentration and then size‐selected using a PippinHT automated DNA size selector (Sage Science Inc., Beverly, MA). The size‐selected samples were ethanol precipitated to concentrate and then checked for quality and quantity using an Agilent TapeStation 4200 with D1000 screentape.

Library sequencing targeted 30 million reads/sample for mRNA and 15 million reads/sample for miRNA. Library sequencing was completed using the standard Illumina HiSeq Rapid Mode protocol. The sequencing library was hybridized to two lanes of an Illumina flow cell at a concentration of 6.0 pM per lane with a 2% Phi‐X control library spiked in. The flow cell was amplified, blocked, linearized, and primer hybridized using the Illumina Onboard Sequencing protocol. Sequencing was completed on the Illumina HiSeq 2,500 (Illumina, San Diego, CA). All mRNA samples sequenced contained between 25 and 80 million raw reads with a raw read length of 100 base pairs. All microRNA samples sequenced contained between 9 and 24 million reads with a raw read length of 50 base pairs. All RNAseq data have been submitted to the National Center for Biotechnology Information Gene Expression Omnibus (accession number GSE200184).

The raw mRNA sequencing reads for each sample were assessed for sequencing quality metrics using fastQC (v.0.10.0) and processed with the fastq‐mcf tool from the ea‐utils software package (v1.1.0) for adapter sequence removal and trimming low quality bases (parameters for minimum remaining sequence length option set to 50 and quality threshold set to 30). The trimmed reads were subsequently mapped to the rat reference genome (version 6.0.80) using the hisat2 mapping software (v2.0.3) run using default parameter settings and in stranded mode. The read alignment BAM files were processed with the HTseq software (v0.7.2) to obtain counts of reads mapped to the antisense strand of each mRNA and tabulated across samples that was used for downstream data normalization and gene expression analysis.

Similarly, the miRNA raw reads for each sample were quality assessed using fastQC (v0.10.0) and were processed for adapter sequence removal using the cutadapt bioinformatics tool (v1.16, with the minimum‐length parameter (−m) set to 14 bases and minimum‐overlap (‐O) parameter set to six bases). The trimmed miRNA reads were mapped to the rat reference genome (version 6.0.80) using the bowtie2 mapping software (v2.3.2) using default parameter settings. Subsequently, the coverageBed function of the bedtools software (v2.17.0) was used to obtain counts of miRNA reads that overlapped miRNA and mRNA features in the rat reference genome annotation (using the ‐s strandedness option enabled). The resulting mapped read counts were tabulated across samples and used for downstream data normalization and miRNA expression analysis.

Any mRNAs and miRNAs that did not meet a minimal expression threshold (at least 10 read counts in at least one sample sequenced) were removed from the downstream analyses. Differential expression analyses on the filtered mRNA and miRNA expression datasets were performed in R using the DESeq2 package (Love, Huber, & Anders, 2014). The DESeq2 normalized data sets were analyzed with principal component analysis (PCA) and a correlation‐based dendrogram to provide observations on sample variations based on mRNA or miRNA expression data as well as to aid in the identification of sample outliers in the datasets. Criteria for assigning differential gene expression were a False Discovery Rate (FDR)‐corrected p‐value ≤0.05 and an absolute fold change ≥1.5.

Total RNA integrity number scores for all placenta samples were ≥ 9.0, but five liver samples had scores <7 (Table S1). Because these five liver samples generally were separated from other replicates within the same group by principal component analysis of the union set of differentially expressed mRNAs (Figure S1A), these five liver samples were deemed outliers and removed from all liver mRNA and miRNA downstream analyses. Using the same two‐step process, no placenta samples were deemed to be outliers (Table S1 and Figure S1B).

Gene Set‐Based mRNA POD determination

BMDExpress software (version 2.2; build 0148) was used to derive a gene set‐based POD (PODGeneSet) for mRNA data (Phillips et al., 2019). Prior to gene expression data model fitting, normalized mRNA RNAseq data were filtered against a Williams trend test p‐value <0.05 and a 1.5 absolute fold change. For genes passing this filter, expression data were fit to Hill, power, linear, polynomial 2, polynomial 3, exponential 2, exponential 3, exponential 4, and exponential 5 models. A best fit model for each gene was selected using the following settings/parameters: 1) maximum iterations of 250; 2) confidence level of 0.95; 3) constant variance; 3) a nested Chi‐square test with a p‐value <0.05 to identify the best polynomial model; 5) power restricted to ≥1; 5) Hill models with a k parameter <1/3 of the lowest positive dose were flagged; when flagged, the next best model with a p‐value >0.05 was used; 6) lowest Akaike Information Criterion value; and 7) a goodness‐of‐fit p‐value >0.1. Because the response level associated with an adverse change in gene expression was unknown, the benchmark response was set to a mean response equal to one standard deviation of the control mean (Davis et al., 2011). Genes with modeled BMDs > the highest dose level or genes with upper bound (95th percentile) BMD values (BMDU) and lower bound (95th percentile) BMD value (BMDL) ratios >40 were removed from further analysis. Using the Functional Classification in BMDExpress, remaining genes were mapped to Gene Ontology Biological Process (GO‐BP) terms (GO file creation date 7/19/19). Based upon optimal concordance of liver omic and systemic apical PODs within the TG‐GATES data set (Johnson, Auerbach, & Costa, 2020), GO‐BP terms with ≤2 genes having BMD values, with <2% of genes in the term having a BMD value, or a Fisher's exact two tail test p‐value <0.05 were excluded. The GO‐BP term with the smallest median gene BMD value was identified, and the final PODGeneSet value was the median gene BMDL value of that GO‐BP term. The BMDExpress data file is provided (File S1).

Gene‐Based mRNA and miRNA POD determination

Gene‐base mRNA and miRNA POD value determination followed the same method as for PODGeneSet determination up to (but not including) the Functional Classification step. mRNAs and miRNAs having BMDU/BMDL ratios >40 were culled. Next, two methods were used to determine the final POD: a First Mode (PODMode) method and an Accumulation Plot Maximum Curvature (PODAccum) method.

The PODMode method was based on the distribution of all BMD values (Pagé‐Larivière, Crump, & O'Brien, 2019). Modes and antimodes were identified as local maxima and local minima, respectively, from distribution kernel density estimates generated by the Sheather and Jones bandwidth selection method (Sheather & Jones, 1991). A minimum probability density of 5.5% was required to be considered a mode. Once modes and antimodes were identified, the first mode was deemed the PODMode.

The PODAccum method was based on a POD identification concept promulgated by Dr. Lyle Burgoon termed Good Risk Assessment Values for Environmental Exposure (GRAVEE) (NTP, 2017). This method defines the POD as the point of maximum curvature of an accumulation plot of BMDL values (Figure S2A). Note that in this method the x‐axis (BMD/L values) in the accumulation plot is displayed in a logarithm scale with base 10; for this reason, all calculations were performed with log‐transformed BMDL values. In step one, the section of the gene BMDL accumulation plot curve that displayed a positive concavity was identified. The first antimode output by the PODMode method was used to limit the portion of the curve to be further analyzed (Figure S2B). Next, the curve was smoothed so that the POD was not restricted to the original BMD/L values (Figure S2C). To maintain the accumulative (continuously increasing) nature of the curve, a shape‐constrained additive model assuming a Gaussian distribution was used to smooth the curve using version 1.2–12 of the R package “scam” (Pya & Wood, 2015). Next, the point of maximum curvature of the smoothed accumulation plot curve was identified using the Kneedle algorithm (Satopaa, Albrecht, Irwin, & Raghavan, 2011). Since the Kneedle algorithm was intended to find the “knee” of a curved line with consistent negative concavity, the smoothed BMDL accumulation plot curve was inverted prior to applying the Kneedle algorithm. When the “knee” (point of maximum curvature of the inverted curve) was identified (Figure S2D), the corresponding point in the non‐inverted BMDL accumulation plot curve was deemed the point of maximum curvature (i.e., the PODAccum) (Figures S2E and S2F).

Apical endpoint POD determination

Both the No‐Observed‐Effect‐Level (NOEL) and BMDL were used to derive the apical POD value. Categorization of apical effects as adverse or non‐adverse was not performed since biological change at the molecular level is a precursor to both apical effects; thus, a molecular POD value can correspond to adverse or non‐adverse apical POD values. The final apical POD value was the lower of the NOEL and BMDL values. Only endpoints that were deemed to have a treatment‐related response to ketoconazole exposure were modeled to identify a BMD‐based POD value.

Current best practices from the scientific community were used for BMD modeling of apical endpoint data (Davis et al., 2011; Haber et al., 2018). BMD/BMDL values were derived using BMDS software (version 3.2) developed by the United States Environmental Protection Agency (EPA). For continuous endpoints, the following models were considered: Exponential (frequentist restricted), Hill (frequentist restricted), Linear (frequentist unrestricted), Polynomial (frequentist restricted), and Power (frequentist restricted) models. The Benchmark Response (BMR) was defined as 10% relative change in central tendency with respect to the modeled control mean for dam body weight, dam body weight change, dam feed consumption, and placenta weight. For fetal body weight, a BMR factor of 5% relative change was chosen instead, as this is the standard choice for reproductive endpoints. For dam liver and placenta dichotomous data, the BMR was defined as 10% Extra Risk, and the following models were considered: Dichotomous Hill (frequentist restricted), Gamma (frequentist restricted), Logistic (frequentist unrestricted), Log‐Logistic (frequentist restricted), Log‐Probit (frequentist unrestricted), Multistage (frequentist restricted), Probit (frequentist unrestricted), and Weibull (frequentist restricted) models. For nested dichotomous data, for which the dichotomous response is measured in the offspring of exposed animals, the Nested Logistic (frequentist restricted) model was used, considering four scenarios: 1) with both litter specific covariate and intra‐litter covariate; 2) without both litter specific covariate and intra‐litter covariate; 3) with litter specific covariate but without intra‐litter covariate; and 4) without litter specific covariate but with intra‐litter covariate. The BMR was defined as 5% Extra Risk, and litter size was used as the litter‐specific covariate. Visual inspection of the plotted dose–response curves of the BMDS results was performed. The goal of this step was to provide an additional indication of how well the model fit the observed data, allowing the user to identify potential model fitting problems that might not be captured by the statistics embedded in BMDS.

All the BMDS analysis results, including visual inspection remarks, are available upon request.

Non‐BMD statistics

The litter was considered the experimental unit. Data from non‐pregnant females were excluded from analysis. Percent post‐implantation loss was calculated with the following equation: [(Number of Implantations – Live Fetuses Observed / Number of Implantations] X 100. Except for fetal death, fetal skeletal, and fetal sex ratio data, endpoints were evaluated for data normality using the D'Agostino and Pearson test. If the normality test p‐value was >0.05, then the data were analyzed by a one‐way ANOVA followed by Dunnett's post‐test. If the normality p‐value was ≤0.05, then a Kruskal–Wallis test ANOVA followed by Dunn's post‐test was used. Percent post‐implantation loss data were analyzed using a censored Wilcoxon test. Fetal sex ratios were analyzed using a binomial distribution test. Fetal bone ossification, cleft palate, and litter resorption data were analyzed using a zero‐inflated binomial model (Cohen Jr., 1966). For all endpoints, a p‐value <0.05 was considered statistically significant.

Final determination of treatment‐related effects for apical endpoints considered statistical analyses, the observed effect size, the presence of a dose–response relationship, and consistency with other biological and pathological findings.

RESULTS

Maternal observations

All mated females survived to study termination (Table 1). All mated females were pregnant except for one in the 0.063 mkd and 2 mkd dose groups, three in 6.3 mkd dose group, and two in the 40 mkd dose group, which was unrelated to ketoconazole treatment (Table 1).

TABLE 1

Summary data for reproduction and fetal observations, maternal organ weight, and maternal plasma estradiol

Dose level (mg/kg/day)	0	0.063	0.2	0.63	2	6.3	20	40
Reproduction observations
Number bred	10	10	10	10	10	10	10	10
Number pregnant	10	9	10	10	9	7	10	8
Number of dam deaths	0	0	0	0	0	0	0	0
Number of moribund dams	0	0	0	0	0	0	0	0
Number of dams removed early	0	0	0	0	0	0	0	0
Pregnancies detected by Stain a	0/0	0/1	0/0	0/0	0/1	0/3	0/0	0/2
Number of totally resorbed litters	0	0	0	0	0	0	0	2
Number of litters with a viable fetus	10	9	10	10	9	7	10	6
Number of corpora Lutea/Dam b	13.7 (1.8)	12.4 (3.0)	13.9 (0.9)	13.2 (1.7)	12.1 (2.1)	11.7 (3.2)	12.6 (2.0)	11.8 (2.1)
Number of implantations/Dam b	13.7 (1.8)	12.6 (3.1)	14.1 (1.2)	13.2 (1.7)	12.3 (2.5)	10.9 (3.1)	12.6 (2.0)	11.8 (2.1)
Number of resorptions/Litter b	0.2 (0.4)	0.6 (0.5)	0.5 (0.7)	0.8 (1.0)	0.4 (0.7)	0.4 (0.5)	0.8 (1.3)	3.8 * (3.1)
Fetal observations
DO Sternebrae (fetuses) c	0/135	NE	NE	NE	2/108	1/72	5/118	25/61 +
Cleft palate (fetuses) c	0/135	0/108	0/136	0/124	0/108	0/72	3/118	36/61 +
Cleft palate (litters) d	0	0	0	0	0	0	2	6 +
Percent post‐implantation Loss b , e	1.7 (3.6)	3.9 (3.7)	3.5 (5.2)	5.9 (7.8)	3.8 (5.9)	3.3 (4.2)	7.5 (14.0)	38.7 (38.6) ^
Viable fetuses/Litter b , f	13.5 (2.1)	12.0 (2.8)	13.6 (1.3)	12.4 (1.8)	11.9 (2.4)	10.4 (3.9)	11.8 (2.9)	7.5 # (4.8)
Sex ratio (male%:Female%)	52:48	48:52	46:54	48:52	48:52	48:52	46:54	48:52
Fetal body weight (g) b , f	5.84 (0.08)	5.90 (0.09)	5.86 (0.11)	5.95 (0.11)	5.85 (0.13)	5.99 (0.20)	5.51 (0.18)	4.44 * (0.21)
Maternal observations
Plasma estradiol (pg/ml) b	29.13 (10.62)	19.46 (12.56)	26.45 (7.59)	28.35 (8.81)	25.21 (11.83)	24.71 (11.15)	19.79 (7.38)	25.72 (3.37)
Gravid uterine weight (g) b , f	107.6 (12.1)	97.2 (18.1)	110.0 (9.0)	102.0 (12.1)	98.6 (19.1)	88.5 (30.5)	106.2 (19.5)	91.6 (15.8)
Placenta weight (g) b , f	0.78 (0.06)	0.81 (0.10)	0.79 (0.07)	0.81 (0.09)	0.86 (0.10)	0.88 (0.17)	1.11 * (0.19)	1.40 * (0.13)
Relative liver weight (g) b , g	3.41 (0.28)	3.52 (0.25)	3.44 (0.19)	3.42 (0.31)	3.60 (0.34)	3.68 (0.43)	3.61 (0.26)	3.80 (0.36)
Relative kidney weight (g) b , g	0.43 (0.03)	0.46 (0.03)	0.45 (0.03)	0.43 (0.03)	0.46 (0.04)	0.46 (0.07)	0.47 (0.04)	0.48 (0.03)

Note: Observations in bold were deemed treatment related. NE: Endpoint not examined within this group.

Censored Wilcoxon test p‐value <0.05.

Kruskal–Wallis one‐way ANOVA and Dunn's post‐test p‐values <0.05.

Ordinary one‐way parametric ANOVA and Dunnett's post‐test p‐values <0.05.

Uteri with no visible fetuses were stained with sodium sulfide to confirm pregnancy status.

Mean (standard deviation).

Number of live fetuses having the observation compared to the total number of fetuses examined.

Number of litters with a fetus having the observation.

Mean Percent/Litter calculated as: [(Number of Implantations – Live Fetuses Observed / Number of Implantations] × 100.

Sexes combined.

Grams of organ weight/100 g of body weight.

Zero‐inflated binomial model p‐value <0.05.

Ketoconazole treatment‐related feed consumption, body weight, or body weight effects were observed in dams at ketoconazole dose levels ≥6.3 mkd. Dose–responsive and treatment‐related decreased feed consumption between 16.4–38.6% was observed during the GD 6–9 interval in dams exposed to ≥6.3 mkd ketoconazole (Table 2). The lower, statistically identified GD 6–9 feed consumption value in the 2 mkd ketoconazole group was deemed spurious and unrelated to treatment due to the lack of a concomitant effect on body weight or body weight gain and the lower feed consumption observed in this group prior to exposure (GD 3–6). Dams in the 20 and 40 mkd dose groups had treatment‐related decreased feed consumption during the remainder of gestation with the largest decreases observed during the GD 18–21 interval. Ketoconazole exposure also resulted in treatment‐related decreased body weight in the 20 and 40 mkd dose groups from GD 9–21 (Table 3). Dose–responsive decreased body weight gain of >34% compared to control was seen at dose levels ≥6.3 mkd during the GD 6–9 interval (Table 4). Except for the GD 9–12 and GD 15–18 intervals, treatment‐related decreased body weight gain was observed at 20 and/or 40 mkd. The statistically identified lower body weight gain value in the 6.3 mkd group during the GD 6–21 interval was deemed spurious and unrelated to treatment because this result was driven by a low value from a single animal which had only two implantations and fetuses.

TABLE 2

Feed consumption summary

Dose level (mg/kg/day)	Days of gestation
	3–6	6–9	9–12	12–15	15–18	18–21
0	20.8 (2.3)	20.7 (1.5)	20.9 (3.0)	21.3 (2.0)	22.5 (2.0)	22.6 (3.1)
0.063	19.9 (2.4)	18.7 (1.9)	20.5 (3.0)	20.7 (2.3)	22.8 (3.5)	22.1 (3.5)
0.2	20.5 (2.4)	19.6 (1.6)	20.9 (1.5)	21.7 (1.8)	23.3 (2.9)	23.5 (2.4)
0.63	19.8 (1.3)	19.2 (1.2)	20.7 (1.6)	20.6 (2.4)	23.0 (1.6)	22.5 (1.9)
2	18.9 (2.0)	18.1 (1.7) a	20.2 (1.5)	20.3 (2.1)	21.8 (3.5)	21.8 (2.9)
6.3	19.6 (2.7)	17.3 (2.5) a	20.2 (2.3)	20.1 (2.3)	21.0 (0.9)	19.7 (1.7)
20	19.9 (1.7)	15.9 (2.0) a	17.6 (1.9) a	18.1 (1.6) a	18.7 (1.3) a	16.5 (3.6) a
40	19.4 (2.4)	12.7 (1.8) a	16.8 (1.4) a	15.4 (2.5) a	15.3 (1.1) a	7.6 (2.6) a

Note: Data are the group mean value in grams/animal/day during the 3 day interval shown. Data shown are the mean (standard deviation). Observations in bold were deemed treatment related.

Ordinary one‐way parametric ANOVA and Dunnett's post‐test p‐value <0.05.

TABLE 3

Dam body weight summary

Dose level (mg/kg/day)	Day of gestation
	0	6	9	12	15	18	21	21(C) a
0	220.0 (10.8)	253.9 (15.1)	267.6 (14.0)	291.4 (17.0)	311.8 (18.0)	352.3 (20.2)	410.0 (26.4)	302.4 (22.2)
0.063	222.7 (10.5)	249.0 (15.4)	261.2 (15.1)	284.4 (21.2)	304.0 (23.4)	341.5 (29.2)	389.9 (38.4)	292.7 (25.1)
0.2	222.3 (8.1)	252.5 (14.4)	266.7 (14.2)	287.4 (15.5)	307.4 (17.1)	347.2 (21.1)	408.7 (23.4)	298.7 (19.8)
0.63	224.2 (10.5)	255.0 (10.8)	270.1 (10.7)	292.8 (10.7)	311.4 (14.0)	348.8 (13.2)	401.5 (22.1)	299.4 (18.2)
2	221.3 (7.0)	247.7 (8.9)	261.3 (9.5)	281.6 (13.3)	302.1 (18.4)	340.1 (18.4)	392.4 (25.4)	293.8 (18.6)
6.3	220.3 (7.3)	249.0 (12.0)	257.9 (13.2)	279.0 (13.3)	296.8 (14.5)	329.3 (14.0)	376.1 (22.1)	287.6 (18.6)
20	220.1 (9.7)	245.9 (10.8)	252.2 (10.4)	269.5 b (10.3)	287.2 b (11.9)	325.9 b (16.7)	366.4 b (29.9)	260.2 b (16.9)
40	219.5 (9.4)	246.7 (12.9)	246.6 b (12.7)	267.0 b (10.2)	279.7 b (13.7)	313.6 b (18.5)	318.6 b (28.5)	232.5 b (15.5)

Note: Data shown are the group mean (standard deviation) in grams. Observations in bold were deemed treatment related.

Gravid uterine weight‐corrected GD 21 body weight.

Ordinary one‐way parametric ANOVA and Dunnett's post‐test p‐value <0.05.

TABLE 4

Dam body weight gain summary

Dose level (mg/kg/day)	Days of gestation
	0–6	6–9	9–12	12–15	15–18	18–21	6–21	6‐21(C)a
0	33.8 (7.5)	13.7 (3.2)	23.9 (4.4)	20.4 (3.2)	40.4 (4.5)	57.8 (9.8)	156.1 (15.2)	48.5 (10.5)
0.063	27.4 (7.7)	10.9 (3.4)	23.2 (8.9)	19.6 (8.4)	37.5 (7.1)	48.4 (11.5)	135.2 (26.7)	40.0 (13.9)
0.2	30.2 (8.5)	14.2 (4.0)	20.7 (4.9)	20.0 (4.9)	39.8 (5.5)	61.5 (6.7)	156.2 (11.4)	46.2 (8.9)
0.63	30.8 (7.9)	15.1 (3.4)	22.7 (3.0)	18.6 (5.1)	37.3 (7.7)	52.7 (11.5)	146.4 (16.3)	44.4 (14.2)
2	26.4 (6.1)	13.6 (2.7)	20.3 (7.2)	20.4 (5.4)	38.0 (5.3)	52.3 (9.0)	144.7 (21.1)	46.1 (15.7)
6.3	28.7 (9.1)	9.0 (4.3)	21.1 (5.1)	17.8 (6.2)	32.5 (9.0)	46.8 (14.2)	127.2* (23.9)	38.7 (12.4)
20	25.8 (6.0)	6.3* (5.0)	17.3 (5.5)	17.8 (2.8)	38.7 (7.5)	40.5* (18.3)	120.5* (23.9)	14.35* (13.9)
40	27.2 (9.7)	−0.1* (5.6)	20.5 (6.2)	12.7 (8.7)	33.9 (7.9)	5.0* (15.4)	72.0* (21.9)	–13.84* (10.0)

Note: Data shown are the group mean (standard deviation) in grams. Observations in bold were deemed treatment related. aGravid uterine weight‐corrected GD 6–21 body weight gain. *Ordinary one‐way parametric ANOVA and Dunnett's post‐test p‐value <0.05.

246.6

318.6

Ketoconazole treatment at dose levels ≥20 mkd resulted in dam liver apical effects. One dam at 20 mkd and two dams at 40 mkd had treatment‐related multifocal, extramedullary hematopoiesis (Table 5). Liver histopathologic findings in dams were limited to the 40 mkd group and included the following: 1) multifocal coagulative hepatocyte necrosis or necrosis of individual hepatocytes with or without accompanying inflammation and 2) vacuolization consistent with fatty change in centrilobular and midzonal hepatocytes. Although the dam relative liver weight value was 11% higher compared to control at 40 mkd (p‐value = 0.08) (Table 1), this observation was not considered treatment‐related since the higher value was not accompanied by hepatocellular hypertrophy. There were no treatment‐related histopathologic effects in the liver of dams given 2 mkd or 6.3 mkd ketoconazole.

TABLE 5

Liver and placenta histopathology summary

Dose level (mg/kg/day)		0	2	6.3	20	40
Number of rats examined/dose level		10	9	7	10	6
Liver
Extramedullary hematopoiesis; multifocal
	‐very slight	0	0	0	1	2
Necrosis; coagulative; hepatocyte; centrilobular/midzonal; with accompanying inflammation; multifocal
	‐moderate	0	0	0	0	1
Necrosis; individual cell; hepatocyte; centrilobular; multifocal
	‐very slight	0	0	0	0	1
Necrosis; individual cell; hepatocyte; centrilobular/midzonal; with accompanying inflammation; multifocal
	‐slight	0	0	0	0	1
Vacuolization; consistent with fatty change; hepatocyte; centrilobular/midzonal; multifocal
	‐slight	0	0	0	0	1
Vacuolization; consistent with fatty change; hepatocyte; individual cells; multifocal
	‐very slight	1	1	1	1	3
Placenta
Necrosis; with accompanying inflammation; decidua; multifocal or diffuse
	‐slight	8	7	2	1	0
	‐moderate	2	2	4	7	1
	‐severe	0	0	1	2	5
Necrosis; with accompanying inflammation; trophospongium; multifocal or diffuse
	‐very slight	7	8	0	0	0
	‐slight	3	1	4	0	0
	‐moderate	0	0	3	4	0
	‐severe	0	0	0	6	6
Necrosis; with accompanying inflammation; artery; multifocal or diffuse
	‐very slight	3	3	3	0	0
	‐slight	7	5	2	4	0
	‐moderate	0	0	1	5	4
	‐severe	0	0	0	0	2
Degeneration; cystic; labyrinth; focal or multifocal
	‐very slight	0	0	4	1	1
	‐slight	0	0	0	5	0
	‐moderate	0	0	1	3	2
	‐severe	0	0	0	1	3
Thickened; trophospongium; diffuse
	‐slight	1	1	2	10	6
Thickened; labyrinth; diffuse
	‐slight	0	0	0	7	5

Note: Only findings with treatment‐related changes are shown. Observations in bold were deemed treatment related.

Apical effects in the placenta were observed at ketoconazole dose levels ≥6.3 mkd. Increased placenta weights of 42.2% and 79.0% were observed in the 20 mkd and 40 mkd dose groups, respectively (Table 1). Placentae from five of seven dams given 6.3 mkd had very slight or moderate cystic degeneration of the labyrinth (Table 5 and Figure 1). Three of seven dams given 6.3 mkd ketoconazole had moderate necrosis with accompanying inflammation of the trophospongium, and five of seven dams at this dose level had moderate or severe necrosis with accompanying inflammation of the decidua. The incidence and severity of treatment‐related placenta histopathologic effects was increased at 20 mkd and 40 mkd ketoconazole (Table 5 and Figure 1). All dams at these two dose levels had diffuse thickening of the trophospongium, and most of these dams also had diffuse thickening of the labyrinth as compared to controls. Cystic degeneration of the labyrinth, ranging in severity from very slight to severe, was present in all dams given 20 mkd or 40 mkd ketoconazole. The cystic degeneration of the labyrinth was characterized by variably sized cystic spaces that were filled with blood and/or fibrin. All of the dams exposed to Ketoconazole at 40 mkd had severe, multifocal or diffuse necrosis with accompanying inflammation of the trophospongium. Eight of ten dams given 20 mkd had moderate or severe, multifocal or diffuse, necrosis with accompanying inflammation of the trophospongium. All dams given 40 mkd and nine of ten dams given 20 mkd had moderate or severe necrosis with accompanying inflammation of the decidua. Five of ten dams given 20 mkd and all dams given 40 mkd had moderate or severe necrosis with accompanying inflammation of the uteroplacental artery. There were no treatment‐related histopathologic effects in the placenta of dams given 2 mkd ketoconazole, and, therefore, histological analyses were not performed at lower dose levels.

FIGURE 1

Placenta histopathology. The following treatment‐related placenta histopathology was observed after ketoconazole exposure: diffuse thickening of the trophospongium (T) and labyrinth (L); cystic degeneration of the labyrinth (arrowheads); diffuse necrosis with accompanying inflammation of the trophospongium and decidua (dashed arrow); and necrosis with accompanying inflammation of the uteroplacental artery (arrow). See Table 5 for incidences

No treatment‐related change in GD 21 dam plasma estradiol concentration or relative kidney weight was observed at any ketoconazole dose level tested (Table 1).

Maternal apical endpoint point of departure

Apical POD values for all treatment‐related effects in the dam ranged from 2 mkd to 38.6 mkd (Table 6). In general, NOEL and BMD‐based POD values were similar for all endpoints; however, dam body weight BMD‐based POD values were up to approximately 6X higher than the NOEL due to use of a 10% BMR. Placenta had the lowest apical POD value derived from the NOEL at 2 mkd for three histopathology endpoints that could not be modeled using BMDS software. The lowest placenta POD that was modeled using BMDS was 2.2 mkd for labyrinth degeneration. The dam liver apical POD was 6.0 mkd for extramedullary hematopoiesis. Among systemic toxicity endpoints, the lowest POD of 2 mkd was for feed consumption during the GD 6–9 exposure interval and was derived from the NOEL. The lowest BMD‐based systemic toxicity POD was 2.7 mkd for GD 6–21 dam body weight gain corrected for uterine weight.

TABLE 6

Apical Endpoint Point of Departure Values

	BMDL	BMD	BMDU	NOEL	LOEL	BMDS note
Dam liver
Extramedullary hematopoiesis	6.0	20.0	NI	6.3	20	BNP
Hepatocyte coagulative necrosis	9.1	40.0	NI	20	40	BNP
Hepatocyte individual cell necrosis	9.1	40.0	NI	20	40	BNP
Multifocal hepatocyte vacuolization	9.1	40.0	NI	20	40	BNP
Hepatocyte individual cell vacuolization	14.5	22.7	NI	20	40	BNP
Placenta
Weight	4.3	5.2	6.4	6.3	20
Decidua necrosis	NI	NI	NI	2	6.3	NSM
Trophospongium necrosis	NI	NI	NI	2	6.3	NSM
Arterial necrosis	NI	NI	NI	2	6.3	NSM
Labyrinth degeneration	2.2	5.0	5.6	2	6.3
Thickened Trophospongium	2.8	5.8	NI	6.3	20	BNP
Thickened labyrinth	3.6	8.4	13.6	6.3	20
Dam feed consumption
GD 6–9	NI	NI	NI	2	6.3	PGF
GD 9–12	13.7	17.8	25.2	6.3	20
GD 12–15	10.9	13.2	16.7	6.3	20
GD 15–18	NI	NI	NI	6.3	20	VAF
GD 18–21	4.0	4.7	5.5	6.3	20
Dam body weight
GD 9	38.6	53.7	87.0	6.3	20
GD 12	35.1	47.9	74.7	6.3	20
GD 15	23.4	55.1	NI	6.3	20	BNP
GD 18	28.7	38.1	56.1	6.3	20
GD 21	17.6	21.4	31.3	6.3	20
GD 21 (corrected)	13.8	16.4	20.0	6.3	20
Dam body weight gain
GD 6–9	3.4	3.9	6.9	2	6.3
GD 12–15	6.6	11.3	31.6	20	40
GD 18–21	5.1	7.3	15.8	6.3	20
GD 6–21	7.3	9.9	16.6	6.3	20
GD 6–21 (corrected)	2.7	3.0	4.8	6.3	20
Fetal
Number of totally resorbed litters	10.9	36.3	39.3	20	40
Number of litters with a viable fetus	NI	NI	NI	20	40	PGF
Number of resorptions/litter	11.7	22.1	NI	20	40	BNP
Viable fetuses/litter	NI	NI	NI	20	40	PGF
Post‐implantation loss	14.0	27.0	NI	20	40	BNP
Fetal body weight	18.7	19.4	20.3	6.3	20
DO Sternebrae (fetuses)	20.3	26.0	NI	6.3	20	BNP
Cleft palate (fetuses)	19.4	22.6	NI	6.3	20	BNP

Note: Values are in units of mg/kg/day.

Abbreviation: BNP, BMDU Not Provided by BMDS; NDR, No Dose Response identified among dose levels; NI, Not Identified; NSM, Not Suitable for Modeling with BMDS; PGF, Poor Global Fit; VAF, Variance Assumption Failed ‐ variance was neither homogeneous nor could it be adequately modeled.

Fetal observations

Treatment‐related effects on fetal parameters were observed at ketoconazole dose levels ≥20 mkd. Treatment‐related increased conceptus resorption was observed at 40 mkd ketoconazole as determined by an increased number of totally resorbed litters and number of resorptions/litter, a decreased number of litters with a viable fetus and viable fetuses/litter, and an increase in mean percent post‐implantation loss (Table 1). Following a complete skeletal examination, treatment‐related findings included a delay in sternebrae ossification and an increase in cleft palate incidence at 20 and 40 mkd ketoconazole (Table 1). Treatment‐related decreased fetal body weights of 6% and 24% were observed at 20 mkd and 40 mkd ketoconazole, respectively (Table 1). No treatment‐related effects were observed for gravid uterine weight or fetal sex ratio at any tested dose level (Table 1).

Fetal apical endpoint point of departure

Fetal apical POD values for all treatment‐related effects were similar and ranged from 10.9 to 20.3 mkd (Table 6). The endpoint with the lowest POD value was the number of totally resorbed litters. The POD value for cleft palate was 19.4 mkd.

Dam liver and placenta differential mRNA and miRNA expression

A treatment‐related change in mRNA expression was observed in dam liver. There were 12,423 mRNAs expressed in dam liver (Table S2). A robust increase in mRNA differential expression was observed at ketoconazole dose levels ≥6.3 mkd (Table 7). Although two mRNAs passed the differential expression criteria at 0.063 mkd ketoconazole, these mRNAs were not considered differentially expressed due to a subset of samples within the group driving the fold change for ENSRNOG00000017775 and the lack of differential expression of these two mRNAs at all higher ketoconazole dose levels. Therefore, the dam liver mRNA No Observed Transcriptional Effect Level (NOTEL) was considered to be 2 mkd.

TABLE 7

Number of differentially expressed mRNAs and miRNAs

Organ	Molecule	Ketoconazole dose level (mg/kg/day)
		0.063	0.2	0.63	2	6.3	20	40
Liver	mRNA	0	0	0	81	3,680	3,490	1,212
Liver	miRNA	0	2	3	34	56	34	15
Placenta	mRNA	5	28	0	0	15	258	878
Placenta	miRNA	0	0	0	2	30	41	75

Values in bold were deemed treatment related. Values in italics were deemed the No Observed Transcriptional Effect Level (NOTEL). Differential expression criteria were a False Discovery Rate‐corrected p‐value <0.05 and an absolute fold change ≥1.5.

A treatment‐related change in miRNA expression was observed in dam liver. There were 276 miRNAs expressed in dam liver (Table S3). The number of miRNAs passing the differential expression criteria at 0.2 mkd and 0.63 mkd ketoconazole were 1 and 3, respectively. The liver miRNA NOTEL was considered to be 2 mkd given the lack of a definitive response below this dose level.

A treatment‐related change in mRNA expression was observed in placenta. There were 14,857 mRNAs expressed in placenta (Table S4). Although five mRNAs passed the differential expression criteria at 0.063 mkd ketoconazole, these mRNAs were considered not differentially expressed since the expression change lacked a dose–response and was driven by a single control sample with a higher expression value (Table S4). The mRNAs passing the differential expression criteria at 0.2 mkd ketoconazole lacked differential expression at adjacent dose levels (Table S4); therefore, differential expression of these mRNAs was considered spurious and unrelated to treatment, and the placenta mRNA NOTEL was considered to be 2 mkd.

A treatment‐related change in miRNA expression was observed in placenta. There were 385 miRNAs expressed in placenta (Table S5). Two miRNAs passed the differential expression criteria at 2 mkd ketoconazole, and no miRNAs passed these criteria at lower dose levels (Table 7). Although only a modest response was observed at 2 mkd, the placenta miRNA NOTEL was considered to be 0.63 mkd.

Dam liver and placenta mRNA and miRNA point of departure

Dam liver mRNA PODGeneSet, PODMode, and PODAccum values were 0.34, 0.69, and 0.62 mkd, respectively (Table 8). The reasoning behind generating molecular POD values using three methods is described in the discussion section. The dam liver mRNA and miRNA accumulation plot of individual molecule BMDL values (used as the starting point to derive the PODAccum) is shown in Figure 2. Dam liver miRNA PODMode and PODAccum values were 0.61 and 0.43 mkd, respectively. The two GO‐BP gene set terms driving the mRNA PODGeneSet were GO:0016441 and GO:0035194. The range of the five lowest GO‐BP gene set BMDL values was 0.34–0.77 mkd, and the range of all GO‐BP gene set BMDL values was 0.34–20.56 mkd. The range of all individual mRNA and miRNA BMDL values was 0.06–37.64 and 0.01–39.65 mkd, respectively.

TABLE 8

Liver and placenta mRNA and miRNA point of departure values

Endpoint	Liver	Placenta
mRNA PODGeneSet	0.34	3.29
mRNA PODMode	0.69	3.96
mRNA PODAccum	0.62	2.53
miRNA PODMode	0.61	6.83
miRNA PODAccum	0.43	5.97
Median BMDL/BMD/BMDU of gene set driving the PODGeneSet	0.34/1.09/3.13	3.29/4.05/10.03
Gene set terms driving the PODGeneSet	GO:0016441 and GO:0035194	GO:0035606
Range of five lowest gene set BMDL values	0.34–0.77	3.29–5.13
Range of all gene set BMDL values	0.34–20.56	3.27–25.54
mRNA first mode BMDL/BMD values	0.69/1.94	3.96/7.14
mRNA accumulation plot maximum curvature BMDL/BMD values	0.62/1.81	2.53/3.34
Range of all individual mRNA BMDL values	0.06–37.64	0.01–38.10
Range of all individual mRNA BMD values	0.31–43.95	0.03–50.88
miRNA first mode BMDL/BMD values	0.61/1.79	6.83/7.51
miRNA accumulation plot maximum curvature BMDL/BMD values	0.43/1.35	5.97/6.88
Range of all individual miRNA BMDL values	0.01–33.61	0.43–32.37
Range of all individual miRNA BMD values	0.05–39.65	2.15–43.74

Numerical values are in units of mg/kg/day.

FIGURE 2

mRNA and microRNA BMDL Accumulation Plots for Dam Liver and Placenta Overlayed with mRNA, microRNA, and Apical Endpoint POD Values. BMDL accumulation plot for all dam liver (red accumulation line) and placenta (green accumulation line) mRNAs (a) and microRNAs (b). Red and green shaded vertical lines are the range of PODGeneSet, PODMode, and/or PODAccum values for liver and placenta, respectively. Blue, purple, orange, and black vertical dashed lines are the lowest apical endpoint POD values for placenta, dam systemic toxicity, dam liver, and fetal effects, respectively

Placenta mRNA PODGeneSet, PODMode, and PODAccum values were 3.29, 3.96, and 2.53 mkd, respectively (Table 8). The placenta mRNA and miRNA accumulation plot of individual molecule BMDL values (used as the starting point to derive the PODAccum) is shown in Figure 2. Placenta miRNA PODMode and PODAccum values were 6.83 and 5.97 mkd, respectively. The two GO‐BP gene set term driving the mRNA PODGeneSet was GO:0035606. The range of the five lowest GO‐BP gene set BMDL values was 3.29–5.13 mkd, and the range of all GO‐BP gene set BMDL values was 3.29–25.54 mkd. The range of all individual mRNA and miRNA BMDL values was 0.01–38.10 and 0.43–32.37 mkd, respectively.

DISCUSSION

Since the pioneering work developing quantitative measures of transcriptome change in toxicology studies (Thomas et al., 2007; Yu et al., 2006), mRNA‐based transcriptome POD values have been compared to traditional apical endpoint POD values across a number of general toxicity studies. These comparisons have led to a developing consensus that a transcriptome POD from a short‐term to subchronic exposure of adult rodents approximates with reasonable accuracy a concurrent and/or chronic exposure apical endpoint POD (Bianchi et al., 2021; Gwinn et al., 2020; Jackson et al., 2014; Johnson et al., 2020; LaRocca, Costa, Sriram, Hannas, & Johnson, 2020; Moffat et al., 2015; Thomas et al., 2011; Thomas et al., 2013). To our knowledge, comparison of a maternal transcriptome POD to an embryofetal apical POD within a developmental toxicity study design and derivation of a POD based upon changes in miRNA have not been published. Thus, goals of this study were to determine if 1) a maternal miRNA and mRNA POD from a surrogate organ (liver) and/or a developmental toxicity key event organ (placenta) would estimate or be protective of a fetal apical endpoint POD and 2) POD values from two omic molecular data types (mRNA and miRNA) were similar.

Unlike other published GD21 rat data on azole fungicides including ketoconazole (Stinchcombe et al., 2013; Taxvig et al., 2008), a decrease in dam plasma estradiol at GD21 was not observed in the current study. The main difference between the current study and published studies was the length of time between administration of the final dose and plasma estradiol measurement. The interval was 24 h in the current study, whereas prior published work administered the azole on the day of estradiol quantification. The lack of an observed decrease in dam plasma estradiol in the current study is likely due to a combination of the short half‐lives of ketoconazole (1.5 h) and estradiol (2–8 h) in the rat and the 24 h interval between the final ketoconazole dosing and plasma estradiol quantification (Jagger, Chow, & Chambers, 1996; Sjöberg, Ekman, & Lundqvist, 1988).

To derive a mRNA and miRNA POD from whole transcriptome data, three methods were used. Three methods were examined since there is not a current scientific consensus on the best practice, and the most common practice of identifying a gene set‐based molecular POD is not possible with miRNA data. The more common approach is one which uses biological knowledge in the form of pathways; individual mRNA POD values are mapped to GO‐BP gene sets, and the final mRNA POD is a gene set‐based POD (NTP, 2018; Phillips et al., 2019). The other two methods (PODMode and PODAccum) have the advantage of not requiring mapping of individual molecules to gene sets, which makes it feasible to generate a POD value for any molecular data type including miRNA. In addition, molecular POD values using the PODMode and PODAccum methods can be derived in species with poor biological annotations. Even in a species like the rat with better gene set annotation, the PODMode and PODAccum methods remove the need to map genes to gene sets which might add unnecessary complexity to the POD derivation method. The PODMode is hypothesized to represent a molecular‐level MoA POD (Pagé‐Larivière et al., 2019). The PODAccum method was based upon the concept that a biological system is an integrated network of molecular and higher order components (Vidal, Cusick, & Barabási, 2011). Such a system transitions between different states via switch‐like changes involving a concerted change in multiple molecular components (Atay, Doncic, & Skotheim, 2016). It was hypothesized that the molecular POD would correspond to initial the point along the dose response continuum where a switch‐like increase in the rate of molecular change was observed. This point (the PODAccum) was defined as the point of maximum curvature along the miRNA or mRNA BMDL accumulation plot line.

The three methods used to derive a mRNA or miRNA POD all generated similar values. For mRNA and miRNA data, POD values using the three methods were all within 2X of each other (Table 8 and Figure 2). This was observed for both liver and placenta data. As concluded by others (Farmahin et al., 2017), this suggests that a transcriptome POD is robust to the method used to derive it. Because the PODAccum method was restricted to the first mode of POD gene level values, the method produced a slightly smaller POD value than the PODMode method. The similarity of the POD value using biological pathway information (PODGeneSet) to POD values not informed by biological pathways (PODMode and PODAccum) suggests that there may be no advantage to using pathway‐based gene sets to obtain a POD value. It is hypothesized that all three methods identify concerted molecular change that is required to initiate perturbations at higher levels of biological organization such as a cell, organ, or organism.

The data reported here support the hypothesis that high content data from different types of molecular endpoints (e.g., transcriptomic and epigenomic) that are collapsed to a single value can be used interchangeably to derive a molecular POD for estimation of an apical endpoint POD. For mRNA‐based POD derivation there are numerous examples of the transcriptome POD to apical POD concordance for general toxicity study designs in the literature (Chepelev et al., 2018; Gwinn et al., 2020; Jackson et al., 2014; Johnson et al., 2020; LaRocca et al., 2020; Moffat et al., 2015; Thomas et al., 2011; Thomas et al., 2013). For epigenetic data such as miRNA, the scientific literature is much sparser. In a study exposing mice for seven days to di(2‐ethylhexyul)phthalate, BMD‐derived values for individual liver miRNAs were higher than the di(2‐ethylhexyul)phthalate mouse liver tumor POD; however, some miRNA and tumor POD values appeared to be within 10X of each other (Chorley, Carswell, Nelson, Bhat, & Wood, 2020). A larger data set comparing miRNA‐based POD values and apical endpoint POD values will be needed to critically examine concordance between these data types.

Critically, the data from this study supports the hypothesis that a rat maternal mRNA‐ or miRNA‐based POD value is protective of a developmental toxicity apical POD. Additional studies will be required to determine if this is a general conclusion or specific to ketoconazole‐induced developmental toxicity. Exposure of rat dams to some azole chemistries (including ketoconazole) causes two types of developmental toxicity via different modes of action: 1) fetal death (i.e., post‐implantation loss or resorptions) via the key events of ovarian aromatase inhibition leading to decreased circulating dam estradiol and subsequent placenta histopathology and dysfunction (Furukawa et al., 2008; Ichikawa & Tamada, 2016; Stinchcombe et al., 2013; Taxvig et al., 2008) and 2) cleft palate via inhibition of cytochrome P450 family 26 (CYP26) protein in embryonic neural crest cells causing altered craniofacial development and cleft palate (Marotta & Tiboni, 2010; Menegola et al., 2003; Tiboni et al., 2009). In placenta, the range of mRNA or miRNA POD values was 2.53–6.83 mkd, and the placenta apical endpoint POD value of 2 mkd was within 2 – 4X of these miRNA or mRNA POD values (Figure 2). Since placenta pathology is a key event in ketoconazole‐induced fetal death, it is biologically consistent that a placenta molecular POD value is protective of a fetal death POD. For the cleft palate fetal endpoint, however, a molecular POD was not determined in a key event biological compartment; nonetheless, the range of dam mRNA and miRNA POD values from both liver and placenta were protective of the cleft palate apical POD (19.4 mkd). Although a maternal molecular POD was protective of the cleft palate POD in this study, it remains to be determined if a maternal molecular POD will be protective of all potential developmental toxicity POD values for developmental toxicants directly targeting embryofetal processes. It may be that determination of a molecular POD in fetal (or neonatal) organs would be required to perform a health‐protective risk assessment for all potential developmental toxicants.

The concordance of dam molecular and apical POD values reported here align with published data showing adult rat general toxicity molecular POD values are typically within an order of magnitude of general toxicity apical POD values (Chepelev et al., 2018; Gwinn et al., 2020; Jackson et al., 2014; Johnson et al., 2020; LaRocca et al., 2020; Moffat et al., 2015; Thomas et al., 2011; Thomas et al., 2013). Interestingly, the dam liver molecular POD values were approximately 30X lower than the dam liver apical POD values (Figure 2). The reason for this 30X lower value is unknown, but it is speculated that the liver molecular POD value may reflect a chronic liver toxicity apical POD. It has been suggested that the liver might serve as a “sentinel” organ for an organism‐wide apical POD which could be used to derive a molecular POD protective of numerous apical POD values (Gwinn et al., 2020; Johnson et al., 2020), and the data reported here support this hypothesis. However, additional data are needed across a wider chemical, organ, and study design space before a definitive conclusion can be made.

Under the conditions of this study, the following conclusions were drawn. Similar to data from rodent models of general toxicity, a rat maternal molecular POD estimates (or is protective of) a developmental toxicity apical endpoint. Molecular POD values derived from different molecular classes (mRNA and miRNA) are similar. These data support the conclusion that a molecular POD can be used to estimate (or be protective of) an apical POD across multiple toxicity study designs and life stages. These data support continued focus within the scientific community on ways to utilize an in vivo molecular POD in industrial chemical and agrochemical human health safety assessments.

FUNDING INFORMATION

This work was funded by the Dow Chemical Company and Corteva Agriscience™.

CONFLICT OF INTEREST

The authors declare that they have no conflicts of interest.

Figure S1

Figure S1B

Figure S1C

Figure S1D

Click here for additional data file.

Figure S2

Click here for additional data file.

Table S1

Click here for additional data file.

Table S2

Click here for additional data file.

Table S3

Click here for additional data file.

Table S4

Click here for additional data file.

Table S5

Click here for additional data file.

Appendix S1

Click here for additional data file.