Dopa decarboxylase is a genetic hub of parental control over offspring behavior

Dopa decarboxylase (DDC) regulates the synthesis of monoaminergic neurotransmitters and is linked to psychiatric and metabolic disorders. Ddc exhibits complex genomic imprinting effects that have not been functionally studied. Here, we investigate different noncanonical imprinting effects at the cellular level with a focus on Ddc. Using allele-specific reporter mice, we found Ddc exhibits dominant expression of the maternal allele in subpopulations of cells in 14 of 52 brain regions, and dominant paternal or maternal allele expression in adrenal cell subpopulations. Maternal versus paternal Ddc allele null mutations differentially affect offspring social, foraging and exploratory behaviors. Machine learning analyses of naturalistic foraging in Ddc−/+ and +/− offspring uncovered finite behavioral sequences controlled by the maternal versus paternal Ddc alleles. Additionally, parental Ddc genotype is revealed to affect behavior independent of offspring genotype. Thus, Ddc is a hub of maternal and paternal influence on behavior that mediates diverse imprinting and parental effects. HIGHLIGHTS Dopa decarboxylase (Ddc) allelic expression resolved at the cellular level Cells differentially express maternal versus paternal Ddc alleles Maternal and paternal Ddc alleles control distinct behavioral sequences Parental Ddc genotype affects offspring independent of mutation transmission eTOC Allelic reporter mice and machine learning analyses reveal dopa decarboxylase is affected by diverse imprinting and parental effects that shape finite behavioral sequences in sons and daughters.


HIGHLIGHTS
• Dopa decarboxylase (Ddc) allelic expression resolved at the cellular level • Cells differentially express maternal versus paternal Ddc alleles • Maternal and paternal Ddc alleles control distinct behavioral sequences • Parental Ddc genotype affects offspring independent of mutation transmission eTOC Allelic reporter mice and machine learning analyses reveal dopa decarboxylase is affected by diverse imprinting and parental effects that shape finite behavioral sequences in sons and daughters.

KEYWORDS
Genomic imprinting, parental effects, foraging, dopamine, serotonin, noradrenaline, dopa decarboxylase, behavior, machine learning, CRISPR knock-in mice, hypothalamic-pituitary-adrenal axis, HPA axis SUMMARY Dopa decarboxylase (DDC) regulates the synthesis of monoaminergic neurotransmitters and is linked to psychiatric and metabolic disorders. Ddc exhibits complex genomic imprinting effects that have not been functionally studied. Here, we investigate different noncanonical imprinting effects at the cellular level with a focus on Ddc. Using allele-specific reporter mice, we found Ddc exhibits dominant expression of the maternal allele in subpopulations of cells in 14 of 52 brain regions, and dominant paternal or maternal allele expression in adrenal cell subpopulations. Maternal versus paternal Ddc allele null mutations differentially affect offspring social, foraging and exploratory behaviors. Machine learning analyses of naturalistic foraging in Ddc -/+ and +/offspring uncovered finite behavioral sequences controlled by the maternal versus paternal Ddc alleles. Additionally, parental Ddc genotype is revealed to affect behavior independent of offspring genotype. Thus, Ddc is a hub of maternal and paternal influence on behavior that mediates diverse imprinting and parental effects.
However, we do not fully understand gene regulatory mechanisms in mammalian cells or how they shape different phenotypes, including brain functions and behaviors. Recent studies suggest important gene regulatory mechanisms remain to be uncovered at the allele level in the mammalian genome (Kravitz and Gregg, 2019). Genomic imprinting is a form of allele-specific gene regulation that evolved in mammals and flowering plants (Haig, 2000a). It involves heritable epigenetic mechanisms that cause preferential expression of the maternal or paternal allele for some genes in offspring and has important roles in the brain (Gregg, 2014;Kravitz and Gregg, 2019;Perez et al., 2016). Imprinting effects impact the phenotypic effects of a heterozygous mutation according to whether it resides in the maternal or paternal allele (Peters, 2014). The function of imprinting is debated (Haig and Haig, 2004;Spencer and Clark, 2014). A deeper understanding of imprinting could improve our understanding of maternal and paternal influences on mammalian phenotypes.
Canonical imprinting involves epigenetic silencing of one parental allele in offspring. However, using RNASeq profiling, we , and others (Andergassen et al., 2017;Babak et al., 2015;Crowley et al., 2015;Perez et al., 2015), described imprinted genes that exhibit a bias to express either the maternal or paternal allele at a higher level. We refer to these cases as "noncanonical imprinting" effects . They are more frequent in the mouse genome than canonical imprinting, exist in wild-derived populations and are enriched in the brain . Others described this phenomenon as a "parental allele bias" or "nonclassical imprinting" (Andergassen et al., 2017;Babak et al., 2015;Crowley et al., 2015;Perez et al., 2015).
However, we observed that these effects appear to be more than an allele expression bias, and may be a cell-specific form of imprinting . Additionally, we found noncanonical imprinting can interact with inherited heterozygous mutations to affect behavior . The cell populations impacted by noncanonical imprinting are largely unknown. The impacts of noncanonical imprinting on protein levels and effects on offspring behavior are also unknown.
Here, we further investigate noncanonical imprinting effects at the cellular level in mice and test different roles in behavior. We focus on dopa decarboxylase (Ddc), an enzyme required for synthesis of monoaminergic neurotransmitters, including dopamine (DA), norepinephrine (NE), epinephrine (E) and serotonin (5-HT). Ddc is located in the genome next to the canonical imprinted gene, Grb10. Grb10 exhibits dominant paternal allele expression in the brain and dominant maternal allele expression in non-neural tissues (Menheniott et al., 2008a;Plasschaert and Bartolomei, 2015).
We previously found Ddc noncanonical imprinting effects that involve biased expression of the maternal allele in the mouse brain . The effects are strongest in the hypothalamus. Ddc imprinting was first uncovered in mice with uniparental duplications and found to exhibit dominant expression of the paternal allele in the embryonic heart (Menheniott et al., 2008a).
The cellular nature and function of maternal and paternal Ddc imprinting effects are not known.
The monoamine system is fundamental. It modulates brain functions, including reward, feeding, motivation, activity, fear, arousal, sensorimotor processes, social behaviors, learning & memory and others (Gershman and Uchida, 2019;Klein et al., 2019;Okaty et al., 2019;Sara, 2009;Volkow et al., 2017). Monoamines have important roles in psychiatric and neurological disorders, including addiction, depression, sleep disorders, obesity, Parkinson's disease and others (Klein et al., 2019;Marazziti, 2017;Sara, 2009). Monoamine signaling also has roles in non-neural tissues (Arreola et al., 2016;Martin et al., 2017). In the adrenal medulla, subpopulations of DA, NE and E expressing cells modulate stress responses as part of the Hypothalamic-Pituitary-Adrenal axis (Kvetnansky et al., 2009). Additionally, 5-HT and DA regulate gene expression through histone modifications involving serotonylation (Farrelly et al., 2019) and dopaminylation (Lepack et al., 2020), and influence embryonic development (Bérard et al., 2019). Previous studies investigating the functional diversity of monoaminergic cells in the brain focused on biophysical properties, transmitter co-release, connectivity, developmental lineage and molecular profiling (Berke, 2018;Farassat et al., 2019;Gershman and Uchida, 2019;Okaty et al., 2019). Currently, we have little understanding of how allele-specific gene regulatory mechanisms contribute to the functional diversification of monoaminergic cells or how Ddc imprinting effects may shape the effects of inherited genetic mutations in the Ddc locus. In humans, DDC genetic variation is significantly linked to multiple biomedically important phenotypes (Watanabe et al., 2019).
Our study further tests the hypothesis that noncanonical imprinting is a cell-type dependent form of imprinting that shapes ethological behavior. We test whether different noncanonical imprinting effects are linked to neuronal versus non-neuronal brain cell-types. We then perform a high-resolution cellular analysis of different neural and non-neural Ddc imprinting effects using targeted knock-in allele-specific reporter mice. To investigate potentially different functional roles for the maternal versus paternal Ddc alleles in shaping behavior, we use reciprocal Ddc heterozygous mice and test how the parental origin of the mutated allele affects social, foraging and exploratory behaviors. Our approach uses behavioral and machine-learning methods we previously developed to analyze naturalistic foraging in mice . Our previous study found that foraging patterns are constructed from finite, genetically-controlled behavioral sequences that we call modules. Here, we discover that multiple forms of parental effects are transmitted through the Ddc locus and affect behavior. Our results uncover distinct functional effects of the maternal versus paternal Ddc alleles, linking each allele to aspects of social behavior and to finite foraging modules in sons and daughters. Additionally, parental Ddc genotype is revealed to cause transgenerational effects on behavior that are independent of the offspring's genotype. Our findings suggest Ddc evolved as a genetic hub mediating different parental influences on offspring behavior.

Defining imprinted gene expression and imprinting effects in major hypothalamic cell classes
For most noncanonical imprinted genes, we do not know the identity of the brain cell-types that express each gene or exhibit imprinting effects. To investigate the expression of maternally expressed genes (MEGs) and paternally expressed genes (PEGs) in different hypothalamic cell types we analyzed public single-cell RNASeq data from adult mouse hypothalamus . We calculated the mean expression level for each imprinted gene in previously reported hypothalamic cell types, including 11 non-neuronal and 34 neuronal cell types . Unsupervised clustering of the results grouped MEGs and PEGs with similar expression patterns across different hypothalamic cell types (Fig. S1). The data show that expression of most imprinted genes is detected in both neuronal and non-neuronal cell populations. However, it is not known whether the imprinting effects occur equally in both cell classes.
To identify brain cell types exhibiting imprinting for canonical and non-canonical imprinted genes, we initially attempted single-cell RNA-Seq (scRNA-Seq) profiling of brain cells purified from adult F1 hybrid mice generated from crossing CastEiJ x C57Bl/6J mice to resolve allele-specific expression and imprinting at the cellular level (not shown). However, we found that ~99% of the >5000 cells profiled are non-neuronal. The cellular bias and known technical noise of the data limited our ability to draw conclusions about allelic expression in subpopulations of cells. Therefore, we developed a different strategy based on studies of gene co-expression networks that found molecular subtypes of brain cells from bulk RNASeq replicates (Kelley et al., 2018;Parikshak et al., 2013;Willsey et al., 2013). These previous methods uncovered co-expression networks by identifying genes with correlated expression patterns across biological replicates. Expanding on this idea, our approach tests whether the magnitude of the imprinting effect for a given imprinted gene correlates with the expression level of genes that are known markers of brain cell types (Fig. 1A,B). Thus, for genes with cell-type specific imprinting, our approach expects the magnitude of the imprinting effect (ie. maternal -paternal allele expression difference) to be positively correlated to the expression of marker genes for the imprinted cell-types, and not correlated to markers in cells lacking the imprinting effect ( Fig. 1A,B).
To compute Imprinting ~ Cell-type Marker Expression correlation networks, we used our published bulk RNASeq replicates for the hypothalamus (arcuate nucleus region) from adult female F1cb (CastEiJ x C57BL/6J, n=9) and F1bc (C57BL/6J x CastEiJ, n=9) hybrid mice . We calculated the allele expression difference (maternal -paternal) for each imprinted gene and replicate. We then computed the expression levels of mouse-human conserved markers of neurons, oligodendrocytes, astrocytes, endothelial cells and microglia (McKenzie et al., 2018). An imprinting ~ marker expression correlation matrix (Spearman Rank) was then calculated to relate the variance in imprinting magnitude (maternal -paternal allele expression difference) to the variance in the expression of cell-type specific marker genes across RNASeq replicates. Internal control studies and methodological details are provided in the supplemental data (Data S1) and methods. We performed the analysis for all autosomal imprinted genes in the mouse  to identify MEGs and PEGs with imprinting effects significantly linked to the expression of neuronal versus non-neuronal cell marker genes (False Discovery Rate (FDR) < 5%). We found 87% of MEGs and 88% of PEGs have significant cell type dependence for imprinting ( Fig. 1B) and defined the major cell classes linked to the imprinting effects for each gene (Fig. 1C,D). Our results define major cell classes driving tissue level detected imprinting effects.
To validate predicted imprinting effects for noncanonical imprinted genes and test whether the imprinting effect manifests as an allele-bias or as allele-silencing effect, we purified neuronal and nonneuronal cells from the adult female hypothalamus for F1cb and F1bc mice (Fig. 1E, STAR methods). We used targeted pyrosequencing of cDNA prepared from the purified cells to evaluate imprinting effects in each cell type for eight noncanonical imprinted genes. Herc3 is a noncanonical MEG and our correlation analysis indicates imprinting in neurons, but not non-neuronal cells (Fig.  1C). Pyrosequencing validated this result and revealed a significant interaction between cross and cell-type, indicating cell-type dependent imprinting with a maternal allele expression bias in neurons ( Fig. 1F). Ctsh is a noncanonical PEG predicted to exhibit imprinting and preferential paternal allele expression in non-neuronal cells (Fig. 1D). Pyrosequencing confirmed this prediction and a paternal expression bias in non-neuronal cells only (Fig. 1G). Pyrosequencing validated the predicted cell-type imprinting effects for seven of the eight noncanonical imprinted genes tested (p<0.05; cross X celltype interaction). All cases showed an allele-bias rather than complete allele-silencing. Therefore, noncanonical imprinting effects continue to manifest as an allelic bias at the level of neuronal versus non-neuronal cells, leaving open the question of whether more pronounced allelic effects occur in more refined cellular subpopulations.

Noncanonical imprinting causes dominant expression of the maternal Ddc allele in discrete subpopulations of brain cells
We next tested the hypothesis that noncanonical imprinting causes dominant expression of one parental allele in small subpopulations of cells and manifests at the protein level. We focused on Ddc and generated Ddc allele-specific reporter mice in which reporter constructs are knocked-in before the Ddc stop codon, placing a c-terminal V5 epitope tag onto the DDC protein for direct detection of one allele and a stable, nuclear eGFP reporter separated from DDC by a P2A cleavage site on the other allele ( Fig. 2A). Our Ddc eGFP line produced robust nuclear eGFP expression in cells for all major monoaminergic nuclei (Fig. S2A). Similarly, the DDC protein is detected in Ddc V5 targeted knock-in mice using immunolabeling with an anti-V5 antibody (Fig. S2A). The two reporter lines were crossed in reciprocal matings to generate compound Ddc eGFP/V5 and Ddc V5/eGFP allelic reporter mice that reveal Ddc expression from the V5 tagged allele and expression of the other allele from nuclear eGFP (Fig.   S2B). Since eGFP is stable with a half-life greater than 26 hours (Corish and Tyler-Smith, 1999), we expect any V5+ and eGFP-brain cells to have stable silencing of the eGFP allele for at least days, enabling us to infer persistent allele-specific expression in affected cells.
We first measured allelic expression in the embryonic day (E)16-17 heart. We found preferential expression of the paternal allele and silencing of the maternal allele in the Ddc eGFP/V5 (Fig.   2B-D) and Ddc V5/eGFP (Fig. 2E-G) reporter mice, confirming a previous report (Menheniott et al., 2008a). However, we also found heart cells that express both Ddc alleles ( Fig. 2D and G, yellow arrows), in contrast to others exhibiting dominant paternal allele expression ( Fig. 2D and G, white arrows). Our mice therefore resolved known Ddc imprinting effects and uncovered previously unknown allelic effects that could only be resolved with a cellular level analysis. We next began a detailed study of noncanonical Ddc imprinting effects in adults.
At the RNA and tissue levels in adults, Ddc noncanonical imprinting with a significant maternal allele bias in brain and liver, and a paternal allele bias in adrenal glands (Babak et al., 2015;. We began by investigating Ddc imprinting at the cellular and protein levels in the adrenal medulla, which contains DA, NE and E expressing cells that modulate stress responses (Kvetnansky et al., 2009). Overall, we found most cells co-express both alleles, but a general paternal allelic effect is apparent at low magnification ( Fig. 2H). At higher magnification, discrete subpopulations of adrenal cells with different allelic expression effects are revealed (Fig. 2H, insets). We found subsets of DDC+ cells that exhibit dominant paternal allele expression (Fig. 2H, blue arrows) and a less frequent subpopulation with dominant maternal allele expression (Fig. 2H, pink arrows). These allelic subpopulations of adrenal cells were observed in both the Ddc eGFP/V5 and Ddc V5/eGFP reciprocal reporter lines (Fig. 2H). Other DDC+ adrenal cells exhibited biallelic expression (Fig. 2H, yellow arrows).
Next, we analyzed the brain. The strongest Ddc imprinting effects in brain are observed in the hypothalamus . In Ddc eGFP/V5 mice, optical sectioning of the hypothalamus revealed DDC+ cells exhibiting dominant maternal allele expression (Fig. 2I, pink arrows) and DDC+ cells expressing both alleles similarly (biallelic cells) (Fig. 2I, yellow arrows). To test whether these cellular subpopulations are also detectable when the parental origin of the allelic reporters is switched, we performed the same analysis in reciprocal Ddc V5/eGFP mice (Fig. 2I). We confirmed small populations of cells exhibiting dominant maternal Ddc allele expression. We did not observe cells with dominant paternal allele expression in both crosses. Therefore, noncanonical Ddc imprinting is more complex at the cellular level. Many brain regions harbor monoaminergic neurons; we do not know the extent of DDC imprinting throughout the brain and whether some regions have imprinted cells dominantly expressing the maternal allele, while others have imprinted cells dominantly expressing the paternal allele. To answer this question, we tested 52 different brain regions in adult female mice for the presence versus absence of cells with dominant maternal or paternal DDC allele expression. We imaged coronal sections from the entire rostral-caudal axis of Ddc eGFP/V5 and Ddc V5/eGFP mouse brains. Every fluorescence image was compared to an adjacent, parallel Nissl-stained section to define the anatomical location of the brain region(s) according to the Allen Brain Reference Atlas. We captured 783 images of DDC+ cell populations from Ddc eGFP/V5 and Ddc V5/eGFP mice to score whether 1) each image contained imprinted neurons (Fig. 3A, AVPV, pink arrows) or not (Fig. 3B, VTA, yellow arrows), and 2) which reporter allele was dominant (eGFP vs. V5). All scoring was performed blinded to brain region and reporter cross (Ddc eGFP/V5 or Ddc V5/eGFP ).
To test whether a brain region is significantly enriched for dominant maternal versus paternal allele expressing cells, the number of images containing eGFP allele dominant cells was compared between the two crosses (Ddc eGFP/V5 and Ddc V5/eGFP ) using a Fisher's Exact Test. The same was done for images containing V5 allele dominant cells (Fig. 3C). The results for each brain region and cross are plotted in Fig. 3C (colored by major brain divisions). Regions scored with dominant maternal allele expressing cells are plotted with a positive p-value and regions scored with dominant paternal allele expressing cells are plotted with a negative p-value. We uncovered 14 brain regions significantly enriched for cells with dominant expression of the maternal allele (Fig. 3C, upper-right quadrant). 38 brain regions were not enriched for imprinted DDC+ neurons (Fig. 3C), revealing brain region specificity for the imprinting that supports and extends our previous findings . No regions exhibited dominant paternal allele expressing cells. Of the 14 brain regions enriched for neurons with dominant maternal Ddc allele expression, nine are located in the hypothalamus (9 of 20 tested; 45%), one is in the pallidum (1 of 2 tested), one in the thalamus (1 of 2 tested) and three reside in the midbrain (3 of 12; 20%). We did not find impacted hindbrain regions (0 of 5 regions in the pons, and 0 of 9 regions in the medulla; Fig. 3C and see full atlas in Data S2 and brain region definitions in Table S1).
To further test the findings from our atlas, we performed an independent and quantitative assessment of allele dominant DDC+ neurons in the anterior ventral periventricular area (AVPV), a top region enriched for maternal allele dominant DDC+ cells (Fig. 3C). Images of the AVPV were taken from Ddc eGFP/V5 and Ddc V5/eGFP mice (n=5 mice per cross), and an investigator blind to the reporter cross (and subject) of each image scored every eGFP+ and/or V5+ cell in the image into one of five allelic categories: eGFP dominant, eGFP biased, biallelic, V5 biased or V5 dominant (Fig. 3D).
We found a significant interaction between transgenic cross and cell category (Fig. 3D, p=0.002, twoway ANOVA), indicating a parent-of-origin effect on the proportion of cells in each allelic category. This result indicates that the relative proportion of eGFP versus Ddc-V5 allele dominant cells depends on the parental origin of the allele. A relative distribution shift toward eGFP+ dominant/biased cells occurs in the Ddc eGFP/V5 reporter cross and a reciprocal distribution shift in the Ddc-V5 monoallelic/biased cells occurs in the Ddc V5/eGFP cross, indicating dominant maternal allele expression ( Fig. 3D). We did observe some eGFP+ monoallelic cells in the Ddc V5/eGFP cross, indicating dominant expression of the paternal allele. However, these cells are rare in the reciprocal Ddc eGFP/V5 cross (~ 3% of cells), indicating they only arise in one cross and likely result from persistent nuclear eGFP expression from an earlier developmental time point. Overall, Ddc imprinting in adults provides maternal and paternal influences over the monoamine system at the cellular level in the brain and adrenal medulla, respectively.

Dominant expression of the maternal Ddc allele in the hypothalamus is linked to subpopulations of GABAergic neurons
Monoaminergic neurons are anatomically, molecularly and functionally diverse. Some are glutamatergic, while some are GABAergic (Trudeau and Mestikawy, 2018). This dictates excitatory versus inhibitory neurotransmission. To explore the nature and function of imprinted DDC+ neurons, we tested whether dominant maternal allele expression is significantly linked to glutamatergic or GABAergic neuronal identity in the mouse hypothalamus. We applied the imprinting ~ marker expression correlation approach introduced above, using marker genes previously found for 18 different GABAergic and 15 different glutamatergic neuron types in the mouse hypothalamus . We first aggregated the marker genes into two collections, all GABA neurons versus all glutamatergic neurons, and created a Ddc imprinting ~ marker expression correlation matrix. The relative expression of the maternal versus paternal Ddc alleles depends significantly on GABA versus glutamate neuron marker gene expression (P = 0.0051, chi-square test) (Fig. S3A). Furthermore, maternal Ddc allele expression is positively associated with GABAergic neuron markers, indicating preferential maternal allele expression in GABA neurons (Fig. S3A). Immunohistochemical triple labeling of GABA, eGFP and Ddc-V5 in Ddc V5/eGFP transgenic mice confirmed that hypothalamic cells exhibiting dominant expression of the maternal Ddc allele co-express GABA (Fig. S3B, white arrows and zoomed images). Expression of GABA is not unique to the imprinted DDC+ neurons, since GABA co-expression was also observed in biallelic DDC+ cells (Fig. S3B, yellow arrows).
We next tested whether expression of the maternal versus paternal Ddc allele depends significantly on the specific subtypes of GABAergic or glutamatergic neurons. We constructed Ddc imprinting ~ expression correlation matrices using the marker genes for each subtype of neuron . We found significant dependence on GABAergic and glutamatergic neuronal subtype (Fig. S3C, P=0.014, chi-square test). The majority (71%) of neuron types positively associated with maternal Ddc allele expression are GABAergic (Fig. S3C), of which GABA13 and GABA12 are GABAergic neurons are most strongly linked to maternal Ddc allele expression. We found Glu13, Glu6, Glu7, Glu15, Glu1 and GABA10 neurons are negatively associated with maternal Ddc allele expression (Fig. S3C). Thus, our results link Ddc imprinting effects to subtypes of GABAergic hypothalamic neurons.

Ddc imprinting shapes offspring social behaviors in a sex dependent manner
We do not know the functions of the Ddc imprinting effects described above. Our data suggest that Ddc imprinting may shape maternal and paternal influences on the Hypothalamic-Pituitary-Adrenal axis, and other neural systems. Thus, Ddc imprinting could affect behavior; thereby impacting the behavioral effects of inherited heterozygous mutations in Ddc depending upon the parental origin of the mutated allele. Given our finding that some brain cell populations exhibit dominant expression of the maternal Ddc allele, our primary hypothesis is that loss of maternal Ddc allele function in offspring significantly affects offspring behavior and the effects differ from loss of paternal allele function.
Alternative, haploinsufficiency and loss of the maternal Ddc allele might have no effect, or loss of the maternal versus paternal allele may cause the same phenotype in the offspring regardless of the parent-of-origin. Imprinted genes are proposed to affect social behaviors (Isles et al., 2006;Ubeda and Gardner, 2010) and we therefore began by testing the effects of Ddc imprinting on sociability and social novelty seeking behaviors.
We generated heterozygous Ddc mutant mice using a reciprocal cross strategy (Fig. 4A,B). In the first cross, a heterozygous mother and wildtype father generated offspring that are either wildtype (Ddc +/+ ) or have a mutated maternal allele (Ddc -/+ ) (Fig. 4A). In the reciprocal cross, the father is heterozygous and the mother is wildtype, generating wildtype (Ddc +/+ ) or mutated paternal allele (Ddc +/-) offspring (Fig. 4B). This design allows us to contrast the effects of maternal and paternal allele deletion. We tested sociability and social novelty seeking in these offspring using the threechamber social interaction test (Moy et al., 2004). To evaluate sociability, we compared the time spent in the chamber with the jail containing a male conspecific (Conspecific) to time in the chamber with an empty jail (Fig. 4C). Male Ddc -/+ mice with a maternal mutant allele spend proportionately less cumulative time in the chamber with the male conspecific relative to the empty chamber compared to their Ddc +/+ littermates (Fig. 4D, mat) (p=0.042; genotype X chamber interaction, mixed linear model).
In contrast, significant effects were not observed in Ddc +/males with a paternal mutant allele compared to littermate controls (Fig. 4D, pat). Other measures, including the cumulative duration spent near the conspecific jail compared to the empty jail, support this result (Fig. 4D). We did not observe significant effects in females. Thus, loss of the maternal Ddc allele significantly reduces sociability with male conspecifics in sons.
Mice prefer to investigate novel mice in their environment compared to familiar mice (Moy et al., 2004). We tested this preference by measuring the time spent in the chamber with a novel mouse (Stranger) versus a familiar mouse (Familiar) (Fig. 4E). We found that Ddc -/+ males with a mutant maternal allele, exhibit a significantly increased preference for a novel versus familiar conspecific compared to their Ddc +/+ littermates (Fig. 4F, mat, p=0.018). A significant effect was not observed for Ddc +/males with a mutant paternal allele (Fig. 4F, pat). Our finding is supported by other measures, including the cumulative duration spent near the novel male jail compared to the familiar male jail ( Fig. 4F, p=0.022). Significant effects were not observed in females. Overall, our results support our primary hypothesis that loss of the maternal Ddc allele affects offspring social behavior and the effects differ from loss of the paternal allele. In addition, we found sex dependent effects in that social behaviors in sons are the most strongly affected.
Testing ethological roles for Ddc imprinting: The modular architecture of foraging uncovered in reciprocal Ddc heterozygous mutants and littermate controls In addition to social behaviors, monoaminergic signaling modulates neural systems linked to reward, anxiety, stress, feeding, learning & memory, arousal, sensorimotor functions and other aspects of behavior. Foraging is a rich ethological behavior involving these neural systems. The developmental transition from maternal care to independent foraging changes offspring demands on maternal resources -a proposed driver of imprinting evolution (Haig, 2000a;Lee, 1996). We recently introduced a behavioral paradigm and statistical methodology to study the mechanistic basis of foraging at the level of finite behavioral sequences that we call modules ).
Here, we tested the hypothesis that loss of the maternal Ddc allele significantly affects offspring foraging and the effects differ from loss of the paternal allele. Moreover, we tested the secondary hypothesis that the maternal and paternal Ddc alleles control the expression of distinct foraging modules in offspring. Alternatively, loss of either parental allele may have a generalized effect on all types of foraging sequences or the maternal and paternal alleles affect the same foraging sequences, but the nature of the effect differs.
In our paradigm, the mouse home cage is attached to a foraging arena by a tunnel and mice are permitted to forage spontaneously in two 30-minute phases. During the Exploration phase, naïve mice express behavioral sequences related to the exploration of a novel environment and the discovery and consumption of food in a food patch (pot 2) (Fig. 5A). Four hours later, the mice are again given 30 min of access to the arena in the Foraging phase, in which the seeds are now buried in the sand in a new location in pot 4 (Fig. 5B). In this second phase, mice express behavioral sequences related to their expectation of food in the former food patch in the now-familiar environment, and also to the discovery of a new, hidden food source (pot 4). By deeply analyzing round trip excursions from the home using an unsupervised machine-learning framework we call "DeepFeats", we uncover finite, reproducible behavioral sequences, which we call modules (Fig. 5C).
Here, we investigated how Ddc imprinting effects shape foraging module expression in offspring.
We captured data for 12,631 round trip foraging excursions from the home and DeepFeats analysis found 19 measures (aka. features) that maximize the detection of candidate modules ( Fig. 5D and S4A-C). Moving forward with these measures, we partitioned the data into training and test datasets balanced by sex, genotype and cross to identify significantly reproducible clusters of similar behavioral sequences that denote modules (Fig. S4D). Our analysis found 170 significant modules from the 12,631 excursions (q<0.1, In-Group Proportion (IGP) permutation test) ( Fig. 5D and S4C,D).
Each module is numbered based on the training data clustering results. Having defined the modular architecture of the cohort's foraging patterns, we analyzed the results to test our hypotheses.

Foraging reveals that Ddc mediates multiple forms of parental influence on offspring behavior
To determine whether Ddc imprinting affects the overall expression of foraging modules, we counted the number of modules expressed by each mouse, which is a measure of the number of times they enter the arena, and analyzed the data according to sex and genotype. Unexpectedly, we found that Ddc -/+ and Ddc +/+ male offspring generated from the cross between Ddc heterozygous mothers and C57BL/6J (B6) fathers express significantly fewer foraging modules than Ddc +/and Ddc +/+ males generated from the reciprocal cross (p=0.003, likelihood ratio test for cross effect, Gamma distribution  (Fig. 5E). Indeed, a significant interaction effect between sex and cross was observed after absorbing variance due to Ddc genotype (het vs wt) in the offspring (p=0.03, likelihood ratio test). Therefore, Ddc heterozygosity in the parents causes a significant phenotypic effect in sons independent of the son's genotype (heterozygous or wildtype). We refer to this form of parental effect as a "cross effect" and it could be caused by the in utero environment or behavior of Ddc heterozygous mother or father.
The presence of significant cross effects on offspring foraging requires that we absorb variance due to the cross effect prior to testing for imprinting effects. We used this nested modeling strategy to test our primary hypothesis that loss of the maternal Ddc allele causes behavioral effects that differ from loss of the paternal allele in sons and daughters. We found that the main effect of the paternal allele genotype (wildtype (+) versus mutant (-) allele) on module expression is not significant (sons: p=0.9, daughters: p=0.2; likelihood ratio test), indicating that loss of the paternal allele does not significantly affect overall module expression levels (Fig. 5E). In contrast, after absorbing variance due to cross and paternal allele effects with a nested model, we found a significant main effect of the maternal allele genotype (+ versus -) in sons (p=0.04, Fig. 5E). Significant effects were not observed in daughters (p=0.3, Fig. 5E), though this does not preclude the possibility of effects on specific modules that are not revealed from overall module expression. Overall, the results support our primary hypothesis that loss of the maternal Ddc allele in offspring significantly affects foraging behavior and has phenotypic effects that differ from the paternal allele. Additionally, we found Ddc imprinting effects on behavior depend upon the sex of the offspring and uncovered Ddc-mediated cross effects on behavior (Fig. 5F).

Ddc imprinting and cross effects control the expression of groups of related foraging modules
Foraging modules are discrete behavioral sequences that differ in terms of their form, function, expression timing, context of expression and regulation . However, some modules have similarities. For example, multiple different modules can be linked to food patch exploitation, each involving a different approach to a common goal -food consumption. Our data reveal Ddc mediates multiple different parental effects on offspring, including maternal and paternal imprinting and cross effects. Here, we determine whether the 170 modules uncovered in these mice include groups of related modules and whether Ddc-mediated imprinting and/or cross effects significantly impact entire groups of related module. Alternatively, different Ddc-mediated paternal effects might only act at the level individual modules, rather than affecting module groups. This analysis begins to test our secondary hypothesis that Ddc imprinting effects impact specific modules and the maternal and paternal Ddc alleles control distinct subsets of foraging modules. Additionally, we test whether cross effects impact specific modules.
We began by defining groups of related modules. For each module, we computed the centroid, defined as the average values of the 19 measures that cluster behavioral sequences together in a module. Centroids summarize how the different behavioral measures delineate different module types and we display them in a heatmap (Fig. 6A). Unsupervised hierarchical clustering of the centroid data revealed 23 module groups (Fig. 6A) and seven partitions that are singleton modules that did not group with others ( Fig. 6A, eg. module-74). For each module group, tested whether the aggregated expression of the modules in the group is significantly affected by cross effects, the genotype of the paternal allele and/or the maternal allele. We found that 27% of groups exhibit a significant main effect of cross in males, and no cross effects in females ( Fig. 6B; p<0.05 in green columns; summary in Fig. 6C). These results reveal that cross effects impact the expression of specific groups of modules in sons. For example, the expression of module Group-1 is decreased in Ddc +/+ and -/+ sons born from Ddc heterozygous mothers compared to sons from wildtype mothers (Fig. 7A). In contrast, there are no Group-1 cross effects in daughters (Fig. 7A) and some other module groups do not have cross effects in either sex, such as Group-5 ( Fig. 7B and see 6B, p>0.05, green columns).
After absorbing variance due to cross effects, we uncovered module groups significantly affected by Ddc imprinting. From males and females, we found 7 groups significantly affected by the paternal Ddc allele genotype (-versus +) ( Fig. 6B; p<0.05, blue columns, main effect of paternal allele) and 5 different groups are affected by the maternal Ddc allele genotype ( Fig. 6B; p<0.05, red column, main effect of maternal allele) (totals in Fig. 6C). Thus, the maternal and paternal Ddc alleles significantly impact different groups of modules. Moreover, we found that the affected module groups differ between males and females (Fig. 6B). For example, the expression of Group-13 is affected by the genotype of the paternal Ddc allele, but not the maternal allele, and the effect is sexually dimorphic (Fig. 7C). Group-13 expression is decreased in Ddc +/sons compared to littermate controls, but significantly increased in Ddc +/daughters (Fig. 7C). On the other hand, Group-30 is affected by the maternal allele genotype, but not the paternal allele, and the parental effect occurs in sons, but not daughters (Fig. 7D).
Interestingly, Group-13 modules are distinguished by strong interactions with the Foraging phase food patch (pot 4) (Fig. 6A, increased "pot4_visits"; and increased "pot2_mean duration"). However, we did not find significant paternal or maternal allele effects on these groups of Exploration phase feeding modules in sons or daughters ( Fig. 6B; p>0.05, blue and red columns in males and females). These results suggest the paternal Ddc allele affects feeding-related behaviors in the Foraging phase, but not the Exploration phase, revealing the importance of the context for revealing phenotypic effects. Moreover, the maternal Ddc allele genotype did not significantly affect any of these feeding-linked module groups, and cross effects significantly impacted Group-2 modules in sons. Our data show that Ddc-mediated maternal and paternal imprinting effects and cross effects differentially affect groups of related modules in sons and daughters, leaving open the question of how these parental effects extend to the level of individual modules.

Identification of individual foraging modules linked to specific Ddc-mediated parental effects
Individual modules are discrete behavioral sequences that we expect are linked to different neural and molecular mechanisms and serve different functions in different contexts. Parents could shape offspring behavior patterns by modulating the expression of individual modules. Alternatively, parental effects might only affect behavior broadly, by affecting module groups or overall behavioral drives, such as hunger. A simple change to hunger, for example, would presumably manifest as a net increase or decrease to all feeding-related modules, rather than changes to specific module groups or individual modules. We tested how Ddc imprinting and cross effects manifest at the level of the individual modules. To begin, we investigated the possibility that parental effects cause novel modules to form in offspring. However, we found that all 170 modules are expressed by Ddc +/+ males and females derived from the cross with a Ddc +/+ wildtype mother, which shows that none of the modules are only expressed in response to a particular parental effect. Thus, parental effects do not appear to cause the formation of new modules in offspring. Instead, the parental effects might shape the expression of a baseline set of modules.
We first assessed whether individual modules within module groups are differentially expressed and therefore behave like distinct units of behavior. In support, our generalized linear model found that 78% of module groups have a significant main effect of individual "module type" in males and/or females, indicating that individual modules within these groups differ significantly in their expression frequency (Fig. S5B,C, p<0.05, main effect of module, dark green Module column). This result suggests the individual modules are distinct units of behavior. To determine how individual modules are impacted by Ddc-mediated parental effects, we examined interactions between module type and different parental effects on module expression, including cross effects, paternal allele genotype and maternal allele genotype (Fig. S5B,C). Post-tests revealed the identities of the affected modules. We found that XXX% of module groups have a significant interaction between module type and one or more parental effects (Fig. S5C), revealing differential effects at the level of individual modules. In one example, Group-9 modules have significantly increased aggregated expression in Ddc -/+ in sons and daughters (Fig. 8A,B). However, at the level of individual modules, we found a significant interaction between module type and maternal allele genotype in sons, revealing module specific effects (Fig. 8C). Post-tests showed that modules -16, -157 and -158 in Group-9 exhibit increased expression in Ddc -/+ sons compared to littermate controls (Fig. 8C). We also observed that module-94 is uniquely impacted by cross effects (Fig. 8C). These individual modules differ in terms of their structure and expression context (ie. phase), suggesting different functions (Fig. 8D).
Overall, Ddc cross, paternal allele and maternal allele effects extend to the level of individual modules and shape offspring behavior patterns, in part, by affecting these discrete units of behavior.
We tallied the total numbers of modules significantly impacted by each of the different Ddc-mediated parental effects (Fig. 8E). Finally, in an extended behavior analysis presented in a Supplemental Data section (Data S3 and S4), we uncovered additional Ddc cross, paternal allele and maternal allele effects on specific aspects of social behavior (Data S3) and in standard lab tests of exploratory behaviors in males and females (Data S4A-D). We also found that none of the Ddc-mediated parental effects significantly alter offspring body weights, indicating no obvious growth or obesity effects (Data S4E). These supplemental data further support our major finding that Ddc mediates multiple different parental effects on offspring behavior.

DISCUSSION
Genomic imprinting is an epigenetic mechanism through which mothers and fathers differentially affect gene expression in offspring. However, the cellular nature and functional effects of noncanonical imprinting (aka. parental allele biases), and the complex imprinting effects impacting Ddc, are not well understood. Here, we defined noncanonical imprinted genes with imprinting effects linked to neuronal versus non-neuronal brain cells and found an allele bias persists at this level of cellular analysis. We then used allelic reporter mice to get finer cellular resolution of noncanonical imprinting for Ddc and uncovered complexities. In adults, dominant expression of the maternal Ddc allele was revealed in subpopulations of cells in 14 brain regions, most of which are in the hypothalamus and we found an association with GABAergic neurons. Dominant expression of the paternal and maternal Ddc alleles was found in different subpopulations of adrenal cells. Paternal allele expression was also observed in subpopulations of embryonic cardiac cells. We found DDC+ cells that co-express both parental alleles in all examined tissues. Functional studies determined that Ddc imprinting effects interact with inherited heterozygous Ddc mutations and shape offspring social, foraging and exploratory phenotypes, but not body weight. Machine learning dissections of foraging in Ddc -/+ , Ddc +/and Ddc +/+ offspring uncovered 170 discrete modules of economic behavior and parental effects at different levels of analysis. The maternal and paternal Ddc alleles are revealed to affect different module groups and individual modules. Additionally, we uncovered Ddc cross effects in which specific modules are affected in sons according to the Ddc genotype of the parents independent of the offspring genotype. Sexually dimorphic cross and imprinting effects were uncovered in social, foraging and exploratory behavior tests. Our results are foundational in that they reveal noncanonical imprinting creates functional cellular diversity, Ddc mediates multiple different parental effects on behavior and how these effects shape ethological behavior in sons and daughters at different levels.

Regulation
Noncanonical imprinted genes that exhibit a bias to express one allele higher than the other have been described by us, and others (Andergassen et al., 2017;Crowley et al., 2015;Perez et al., 2015). It is debated whether these effects are functional or an epiphenomenon.
Our current study supports a functional role for these effects by revealing that (1) the Ddc maternal allele bias in the adult mouse brain involves dominant maternal allele expression in discrete subpopulations of brain cells and manifests at the protein level, (2) the effects are linked to GABAergic neurons, (3 we uncovered subpopulations of adrenal cells with maternal and paternal Ddc imprinting, and (4) loss of the maternal Ddc allele affects behavior differently from loss of the paternal allele. In some cells, we observed an allele bias and it is unclear whether these cellular allele biases are functional and contribute to behavioral changes or whether the behavioral effects are related to cells with allele silencing effects. Cell-type dependent imprinting is a known phenomenon (Gregg, 2014;Perez et al., 2016) and, as our understanding of the cellular regulation and function of noncanonical imprinting grows, the value of differentiating these effects from canonical imprinting may change. Determining the connectivity patterns and physiological properties of monoaminergic cell populations distinguished by different imprinting effects will deepen our understanding of the functional cell types and circuits. Finally, studies profiling imprinting effects in mice differ in terms of their power to detect different effects and the total numbers of imprinted genes are debated (Andergassen et al., 2017;Babak et al., 2015;Crowley et al., 2015;Perez et al., 2015). With large sample sizes and a sensitive statistical method, we, and others, found noncanonical imprinting is more prevalent in the mouse genome than canonical imprinting Crowley et al., 2015). However, the cellular nature and functions of most cases have yet to be studied.

Effects within the Monoamine System
Previous studies uncovered various mechanisms of parental influence on offspring phenotypes, including the effects of maternal and paternal behavior (Alter et al., 2009;Curley and Champagne, 2016), the in utero environment (Lindsay et al., 2019), maternally imprinted genes and paternally imprinted genes (Haig, 2000a;Keverne, 2001;Perez et al., 2016). Our data indicate that Ddc is a single genetic locus that contributes to many, and perhaps all of these different parental effects. We found Ddc exhibits imprinting of the paternal and/or maternal allele in subpopulations of DDC+ adult brain cells, adrenal cells, and developing cardiac cells. Grb10, which is located next to Ddc in the genome, exhibits a different imprinting pattern, involving dominant expression of the paternal allele in the brain and the maternal allele in non-neural tissues (Arnaud et al., 2003;Blagitko et al., 2000;Hitchins et al., 2001;Miyoshi et al., 1998). Most imprinted genes exhibit imprinting of the same parental allele across all affected tissues and developmental stages (Andergassen et al., 2017;Crowley et al., 2015;Perez et al., 2015). Therefore, Ddc and Grb10 imprinting effects are atypical, suggesting they are subject to a unique crucible of maternal and paternal evolutionary pressures. In addition to these unusual imprinting effects, we found Ddc-mediated cross effects on behavior. Collectively, the data suggest that Ddc has evolved to be a hub mediating multiple different forms of parental influence on offspring behavior in mice. imprinting effects to specific behavioral phenotypes and modules. Conditional genetic approaches will also help to further dissect Ddc imprinting effects from cross effects.
Ddc-mediated cross effects on offspring behavior could be caused by changes to maternal physiology, the in utero environment and/or parental behaviors. A previous study found that Tph1 -/mothers lacking the enzyme synthesizing peripheral 5-HT have smaller embryos independent of the genotype of the offspring (Côté et al., 2007). The authors did not report an effect in Tph1 heterozygous mothers, which were used as controls. Our study found that Ddc heterozygosity in parents significantly affects offspring behavior without significant effects on body weight. Future studies are needed to determine whether Ddc cross effects are mediated by 5-HT, DA, NE and/or E signaling. Roles for maternally-derived 5-HT in embryonic development are the best understood (Brummelte et al., 2017). Finally, imprinted genes can significantly affect parental behaviors, such as Peg1/Mest (Lefebvre et al., 1998) and Peg3 Li et al., 1999) (but see (Denizot et al., 2016)). Therefore, Ddc-mediated cross effects might also involve changes to maternal or paternal behaviors that in turn shape offspring behavior. Future studies are needed to determine the mechanistic basis of Ddc cross effects.
Ddc cross effects most strongly impacted the behaviors of male mice. DDC modulation of in utero 5-HT levels might influence the development of neural circuits that are sensitive to the effects of sex hormones; thereby shaping behavior in a sexually dimorphic manner. Our data raise the possibility that genetic variants impacting other enzymes required for monoamine biosynthesis could also cause transgenerational effects on offspring behavior. Deep analyses of rich foraging patterns in mouse models could help to test the transgenerational behavioral effects of some known genetic variants associated with human phenotypes (Watanabe et al., 2019). There are potential biomedical implications for our observation because 5-HT levels are proposed to influence autism risks, which is three-fold more prevalent in males than females (Muller et al., 2016;Veenstra-VanderWeele et al., 2012).
Of potential importance is that Ddc has been reported to show imprinting in mouse extraembryonic tissues involving paternal allele expression in the visceral yolk sac epithelium and yolk sac (Andergassen et al., 2017). If loss of Ddc expression in extraembryonic tissues causes transgenerational effects on offspring behavior, this mechanism could link Ddc cross effects and imprinting effects into a unified transgenerational genetic axis of parental influence. Thus, further studies are needed to determine whether the parental origin of a Ddc heterozygous mutation in the parents significantly influences Ddc-mediated cross effects on offspring behavior and whether environmental factors, like diet, impact cross effects. There may be adaptive advantages to shaping offspring foraging through parental effects on monoamine biosynthesis pathways (Brummelte et al., 2017). Such effects might shape offspring survival and reproductive success by changing the expression of discrete foraging modules to improve offspring foraging success in the environment.
Indeed, foraging patterns are under strong selective pressures that shape life histories and must be well adapted to environmental conditions (Lee, 1996;Nislow and King, 2006;Stephens et al., 2007).
The functions of parental effects on offspring social behaviors have been considered by others (Haig, 2000b;Isles et al., 2006;Ubeda and Gardner, 2010). Further studies are needed to uncover the potential adaptive advantages of different Ddc-mediated parental effects.
Finally, DDC imprinting appears to occur in humans, though follow up studies are needed to confirm and better characterize the effects (Babak et al., 2015;Baran et al., 2015). Large-scale genome-wide association studies have linked genetic variation at the Ddc locus to risks for multiple major diseases, including bipolar disorder, schizophrenia, depression, addiction, anorexia nervosa, type 2 diabetes, obesity and other biomedically-important behavioral, metabolic, cardiovascular, endocrine and immunological phenotypes (Watanabe et al., 2019). These disorders are phenotypically variable and the mechanisms involved are unclear, but could involve DDC imprinting and cross effects.

Resolving and Regulating the Modular Architecture of Complex Foraging Patterns
Vertebrate economic choices that shape reward, effort and risk evolved under selective pressures for foraging and function to help balance caloric intake with the costs of energetic demands and predation risks (Stephens et al., 2007). Understanding the architecture and regulation of these complex economic behavior patterns is important and fundamental. The field of neuroeconomics proposes that many chronic health problems have roots in reward, effort and risk decisions (DeStasio et al., 2019). In addition, mental illnesses, like bipolar disorder, schizophrenia and depression, have been framed as disorders of neuroeconomic processes (Sharp et al., 2012). However, the mechanisms shaping complex economic behavior patterns are not well understood. Genomic imprinting evolved in mammals and the nursing to foraging transition at weaning uniquely shaped mammalian life histories and brain development (Lee, 1996), suggesting that different parental and imprinting effects may strongly shape mammalian foraging and economic behaviors. Previously, we found support for this idea and discovered foraging patterns in mice are constructed from finite, genetically controlled modules defined from round trip excursions from the home . Here, we have begun to dissect the parental, molecular and cellular mechanisms controlling foraging modules. So far, we found Ddc-mediated imprinting and cross effects do not cause new modules to form in offspring, but affect the expression frequency of specific module groups and individual modules in a sex dependent manner. Future studies are needed to test potential effects on module expression timing and sequential order to better understand how Ddcmediated effects shape higher level economic patterns. We also do not know whether complex social behaviors can be dissected into finite modules controlled by Ddc or what foraging modules might be revealed in richer and more naturalistic foraging environments. Our understanding of how different parental, molecular and cellular mechanisms shape complex behavior, health and disease phenotypes will improve with deeper machine-learning dissections of behavior in the future.

SUPPLEMENTAL INFORMATION
The supplemental information includes five supplemental figures, one supplemental table and a supplemental data results section.

ACKNOWLEDGEMENTS
The statistical methodology in the study was supervised and developed by statisticians in the University of Utah Biostatistics Core Facility, including Drs. Greg Stoddard and Alun Thomas. We

DECLARATION OF INTERESTS
A patent has been filed on the DeepFeats algorithm. C.G. is a co-founder of Storyline Health Inc., which is building artificial intelligence technologies for behavior analysis.       (E-F) The data show phenotypic effects uncovered in the three-chambered test for social novelty (E).
Significant phenotypic effects were uncovered in Ddc -/+ maternal allele mutant heterozygous males compared to their Ddc +/+ littermates (mat) (F), but not in Ddc +/paternal allele mutant heterozygous males (F, pat). Ddc -/+ males spend relatively more cumulative time in the chamber (F, top graph, mat) and in the jail area (F, bottom graph, mat) with the novel stranger compared to the familiar conspecific. Statistical analyses involved a generalized linear mixed effects model testing for an interaction between chamber/jail and genotype. *P<0.05, N.S., not significant. (F) Schematic depicts the cross between Ddc heterozygous mothers X wildtype fathers to generate Ddc -/+ and Ddc +/+ littermate offspring (top), and the reciprocal cross between Ddc heterozygous fathers X wildtype mothers to generate Ddc +/and Ddc +/+ offspring (bottom). This paradigm enables us to uncover the behavioral effects of Ddc imprinting and link loss of maternal versus paternal allele function to specific behavioral phenotypes (red box), and reveal and account for cross effects due to heterozygosity of the parents independent of offspring genotype (green boxes).     (E) The pie charts show the total numbers of foraging modules significantly impacted by parental effects, including the maternal Ddc allele genotype (red), paternal Ddc allele genotype (blue) and cross effects (orange) in males and females. The counts include all modules in module groups with a significant parental effect and the specific modules impacted for module groups that have a significant interaction effect between module type and parental effect (see Figures 6 and S5). Tables   Table S1.

Key to atlas of adult brain regions containing DDC+ neurons with dominant maternal allele expression. Key to Figures 3 and Data S2
The table provides definitions for points plotted in Figs. 3C and S2. Abbreviated Region -A brain region(s) is labelled according to the online Allen Mouse Brain Reference Atlas version 1: coronal (https://mouse.brain-map.org/experiment/thumbnails/100142143?image_type=atlas).
Images of DDC+ cell clusters captured by the 20X coronal field of view and grouped together for statistical analysis often contained more than one anatomically defined brain region in the Reference Atlas and are therefore named together with hyphenations. Brain Region -Definitions for the brain region (

Lead Contact
• Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Christopher Gregg (chris.gregg@neuro.utah.edu).

Materials Availability
• Plasmids generated in this study are available from the Gregg lab and will be deposited into AddGene in the future with peer-reviewed publication of the paper.
• Mouse lines generated in this study are available from the Gregg lab and will be deposited into JAX labs in the future with peer-reviewed publication of the paper.

Data and Code Availability
• All code is available from the Gregg lab • All data is available from the Gregg lab and will be deposited into public repositories in the future with peer-reviewed publication of the paper.

Housing and husbandry
All experiments were conducted in compliance with protocols approved by the University of breeders (6weeks to 1year of age) were paired continuously, and pups were weaned at postnatal day 21 (P21) and cohoused with up to five same-sex littermates or similar aged samesex mice of the same line; mice were never singly housed. As needed, ear punches were taken at P7 (Ddc Hets) or P17-P21 for both genotyping biopsy samples and mouse identification purposes. Before dissections of brain and embryonic heart tissue, mice were put to sleep with isoflurane gas anesthesia and decapitated.

METHOD DETAILS), and the University of Nebraska Medical Center Mouse Genome
Engineering Core Facility (Omaha, NE) used these constructs to perform CRISPR mediated homology directed repair for targeted insertion of the reporters immediately before the stop codon of the gene dopa decarboxylase (GRCm38/mm10; chr11:11815230) into C57BL/6J mice (see Easi-CRISPR Targeted Mutagenesis in METHOD DETAILS). Targeted insertion into the genome of F0 founder mice was confirmed by PCR amplifications using primers that flanked the 5' and 3' genome integration sites, and sanger sequencing confirmed that the reporter was inframe with the Ddc coding sequence, contained no indels, and no non-synonymous amino-acid substitutions.
F0 founder mice were shipped to the University of Utah where they were backcrossed into the C57Bl/6J background strain for five generations to reduce propagation of any potential unknown off-target CRISPR mutations before producing homozygous reporter lines. We found that the mRuby reporter was not detectable and directly labeled for the V5 protein tag to detect expression in the line referred to as Ddc V5 . Reciprocal crosses of Ddc eGFP dams X Ddc V5 sires and Ddc V5 dams X Ddc eGFP sires produced Ddc eGFP/V5 and Ddc V5/eGFP offspring, respectively, for microscopy studies. All brains and adrenal glands from reporter mice were collected from adult females (Atlas, P65; AVPV P79-197), and expression of the reporters was restricted to brain regions with known monoaminergic cell populations. Embryonic heart was collected from both sexes between E16-18.

Genotyping
Ear punches taken at P7-P21 were lysed in 75µL of 25mM NaOH + 0.25mM EDTA with a 1 hour incubation in a thermalcycler at 98˚C. Lysates were then pH neutralized with an equal volume of 40mM Tris.HCl, pH5.5. Two µL of lysates were then added to make 20µL PCR reactions with DreamTaq Green Master Mix (ThermoFisher, K1081) and 0.5µM primers (Table: Genotyping Primers).

Pyrosequencing validation of cell-type specific imprinting
Adult female F1bc and F1cb reciprocal hybrid mice (see above) were euthanized, brains were Pyrosequencing Allelic Quantification (AQ) analysis measured the ratio of maternal to paternal allele expression of imprinted genes from neuron and non-neuron cell fractions isolated from F1cb and F1bc mice using the Pyromark Q24 System (Qiagen, 9001514) according to manufacturer protocols and our published methods . PyroMark Assay Design Software designed amplification and sequencing primers (Table: Pyrosequencing AQ Assay Primers) for AQ assays that contain strain distinguishing SNPs (Cast vs B6) within the genes. Four to six samples per strain (F1cb and F1bc) and cell-type fraction (neuron vs. nonneuron), four groups in total, were PCR amplified to measure allelic expression for each gene.
An interaction effect between strain (F1cb and F1bc) and cell type (neuron and non-neuron) in two-way ANOVAs validated cell-type specific imprinting effects in the brain.

Tissue preparation
Under deep isoflurane anesthesia, adult reporter mice were transcardially perfused with PBS, to clear blood, followed by ~25-50ml of 4% paraformaldehyde to fix tissues. After perfusion, brains were dissected from the skull and adrenals from the abdomen, post-fixed in 5ml of 4% paraformaldehyde for ~48hrs at 4˚C, cryoprotected by immersing serially in 15% and then 30% sucrose at 4˚C until the tissues sank, embedded and frozen into OCT in cryomolds using a methanol dry-ice bath, and stored at -80˚C until sectioning. Coronal brain sections from regions of interest were then cut on a cryostat at 20µm thickness to directly mount onto Superfrost Plus slides (Fisher, #12-550-15), air-dried, and stored in a slide box at -80˚C until Nissl stained or immunolabelled. Brain and adrenal tissue labelled as floating sections were processed as above except they were cut at 30µm thickness, collected into vials of antifreeze solution (300g sucrose, 10g polyvinylpyrrolidone 300ml ethylene glycol, and 0.1M phosphate buffer to 1L) and stored at -20˚C until immunolabeling.
For embryonic heart tissue, pregnant mothers were euthanized under isoflurane anesthesia before dissection out the uteri containing embryos into ice-cold PBS in a petri dish.
Embryos were then dissected from the uterus and extraembryonic tissue and transferred to a clean dish with fresh PBS. Embryonic hearts were micro-dissected under a stereo microscope, immersed in 4% paraformaldehyde for ~48hrs at 4˚C, exchanged with 15% and 30% sucrose, frozen in OCT molds made from aluminum foil, and stored at -80˚C. 10µm embryonic heart sections were cut on a cryostat, mounted onto Superfrost Plus slides, and stored in a slide box at -80˚C until staining.

Nissl Stain
For the brain Atlas of Ddc maternal dominant expressing cells, adjacent parallel sections were Nissl stained to identify neuroanatomical landmarks according to the Allen Mouse Brain Reference Atlas. Slides with 20µm sections were immediately submerged in 1:1 ethanol/chloroform overnight, then rehydrated through 100%, 95%, 70%, 50% ethanol to distilled water. Sections were then stained for 10 min. in prewarmed 0.1% cresyl violet made fresh with 0.3%v/v glacial acetic acid in a 45˚C oven. Slides were then quickly rinsed in distilled water and differentiated for 2-30 min. in 95% ethanol while checking microscopically for best result. Sections were then dehydrated in 100% ethanol 2 X 5min., cleared in xylenes 2 x 5min, and cover-slipped with EcoMount (Biocare Medical, EM897L) mounting medium.

Immunolabelling
Slides with 20µm brain sections (Figs. 2I, S2, 3) were removed from -80˚C storage, air-dried and hydrophobic barrier pen was drawn around sections. By Immersing slides in Copland jars on an orbital shaker with gentle-agitation, sections were hydrated in PBS, permeabilized for 10min. in PBS + 0.2% triton-X 100 and washed 3 X 5 min. with PBS + 0.025% triton-X 100. In a humidified chamber sections were covered and blocked with 10% normal donkey serum (Lampire Biological Products #7332100) and 1% BSA in PBS for 2hr. at room temp., and then incubated overnight at 4˚C with goat polyclonal anti-V5 primary antibody (abcam ab9137) diluted 1:1000 in primary buffer (1% BSA in PBS). The next day, sections were washed 3 X 5 min. with PBS + 0.025% Triton-X 100 in Copeland jars at room temp. with gentle agitation.

DDC Brain Atlas
The entire fixed brains of Ddc GFP/V5 and Ddc V5/GFP adult female mice were sectioned at 20µm thickness and mounted onto positive charged slides in five parallel series containing representative coronal sections along the entire rostral caudal axis. Every section from one series for both crosses was fluorescently immunolabled for V5 as described above to detect cells expressing the Ddc V5 allele, while expression of the Ddc GFP allele was visualized by native GFP fluorescence concentrated in the nucleus. A second series with parallel sections to the immunolabelled series was Nissl stained. Epifluorescence images of DDC+ cell clusters were captured on a Zeiss AxioImager 2 as optical slices using ApoTome2 structured illumination, a 20X objective, fluorescence filters (DAPI, GFP/488, TexasRed/568, 647) and Zen2 software.
Nissl stained sections were captured with brightfield illumination, a 5X objective, and multiple images capturing entire coronal sections were stitched together using a motorized stage and Zen2 tiling. The fluorescence images captured in the camera's 20X field of view were compared to an adjacent, parallel Nissl stained section and the Allen Brain Reference Atlas to identify and label the anatomical location(s) of the DDC+ cells (Fig. 3A

Foraging
Foraging behavior testing was performed as detailed previously .
The detailed protocol published in the supplementary material for this previous study was followed.

Open field
Two 40x40cm open field arenas enclosed by 40cm high walls, made from white melamine laminated (impervious to liquids) pressboard, were placed adjacent to each other on the testing room stage. Arenas were illuminated to 170lux with two track-lighting, bendable gooseneck LED fixtures positioned to minimize both shadows and glare from light-source reflections.
Ethovision Experimental Settings used two cameras at a resolution of 800 x 600 with a merged image of both arenas and center-point detection of subject. Detection Settings used "automatic setup" with a C57Bl/6J mouse exploring the arena, settings were adjusted so that missed samples were minimized (i.e to zero), and different test days used duplicated detection settings with an updated image capture of the arena background (arenas without a subject). Using a background image Arena Settings were adjust to defined arena areas, zones, and calibrated distances to include a 12cm wide wall-zone perimeter surrounding a 16x16cm 2 center-zone.
New trials were defined for every two mice tested simultaneously.

Elevated Zero Maze
A single elevated zero-maze (Stoelting, 60190) for mice was placed on the test-room stage, and was recorded in the dark with infrared illumination and a single camera mounted above. The zero maze consisted of a circular platform with a 50cm inner-diameter and outer-diameter, that is held above the stage with leg supports. The maze is divided into equal quadrant sectors: two closed sectors on opposite sides of the maze enclosed by high inner and outer walls, and inbetween closed sectors were two equally sized open sectors without walls. The inner walls were made with infrared transparent plastic to track the mouse with the camera through the walls.
Ethovision Experiment Settings were set to 800 X 600 resolution with one arena and center-

Light Dark Box
A single light-dark box (Stoelting 63101) was set on the testing-room stage underneath a camera, LED white lights, and infrared lights mounted above. The box was separated into light and dark chambers by a central divider with a door in the middle for the mouse to pass freely between the two sides. The dark chamber outer walls, lid, and central divider were made of plastic that blocked white light but transmitted infrared light for tracking mice in the dark with the infrared sensitive camera. The light chamber outer walls were made from transparent plastic and had no lid so that it could be illuminated from above by aversive white light (~100-200 lux).
Ethovision Experiment Settings were set to 800 x 600 resolution with one arena and centerpoint detection. Arena settings defined the arena area, the light and dark chamber zones, and calibrated distance. Trial Control was set to start tracking when the subject was detected in the maze for two seconds and stop tracking after 20 min. Detection Settings were adjusted with "automatic setup" to minimize missed and subject not found samples while a C57Bl/6J mouse explored the maze; different test days used duplicated detection settings with an updated background image. Smoothing was set to "high". A single subject was tested in each trial of the Trial List, which defined the subjects and their independent variables. After capturing the background image and starting the trial under Acquisition, the subject mouse was lifted by the tail from its' home-cage and gently placed into the center of the light side facing away from the researcher with the door to the dark side on the right side of the mouse. The curtains were then closed around the stage and the experimenter left the room for the 20 min. trial. After testing all subjects, the Ethovision Analysis Profile was set to export the following measures: total distance moved (cm), maximum velocity (cm/s), frequency of entries into the light side (n), cumulative duration in the light side (s), latency to enter the dark side (s), number of transitions from the dark side to the light side (n), and mean, minimum and maximum acceleration (cm/s 2 ).

Social Preference and Social Novelty
A standard three-chambered box (Stoelting 60450), with center, left and right chambers separated by doors, was used to test social preference and social novelty behavior. In the week prior to testing, social stimulus conspecific adult C57Bl/6 males were habituated to being confined to the round jail-cells used in the three chambered box for 30 min. for 3-5 days. Tests were recorded in the dark with infrared illumination, and two mice were tested in separate mazes at a time. In the social preference task one chamber (left or right) was occupied with a conspecific confined to a jail-cell, while the opposite chamber (right or left) contained an empty jail-cell. To begin the social preference task, the test mouse was habituated to the center chamber for five minutes with the doors closed to the left and right chambers. The experimenter then started Ethovision recording and removed the doors between the chambers to allow the test mouse to investigate the chamber and jail-cell containing the conspecific, and the empty chamber and jail-cell, for 10 min. Immediately following the 10 min. social preference test, the test mouse was returned to the center chamber with the doors closed to begin the social novelty task. For the social novelty task, a novel "stranger" conspecific was placed into the previously unoccupied jail cell while the now "familiar" conspecific used in the preceding social preference task remained in its' jail cell. The experimenter then removed the doors while Ethovision recorded the test mouse investigating the chambers and jail cells with either the "familiar" or "stranger" conspecific for an additional 10min. trail. Conspecific stimulus mice were not reused in consecutive tests, and use of the left or right chamber for the first/familiar and empty/stranger conspecifics was counterbalanced by sex and genotype. Ethovision Analysis Profile exported the following measures for the social preference task: total distance moved (cm); maximum velocity (cm/s); cumulative duration (s), frequency of entries (n), mean duration per entry (s), and latency to enter (s) the conspecific chamber; cumulative duration (s), frequency of entries (n), mean duration per entry (s), and latency to enter (s) the empty chamber; cumulative duration (s), frequency of entries (n), and mean duration per entry (s) in center chamber; cumulative duration (s), frequency (n), mean duration (s), and latency(s) for nose to investigate conspecific jail (s); cumulative duration (s), frequency (n), mean duration (s), and latency (s), for nose to investigate empty jail; and fecal boli (n). The same measures were exported for the social novelty task, except the chambers and jails contained either a "familiar" or "stranger" conspecific.

Imprinting ~ Expression Correlation Network Analysis
Published RNASeq maternal and paternal allele expression data and gene expression data for the arcuate nucleus region of the adult female hypothalamus for F1cb (n=9) and F1bc (n=9) mice was obtained . For each biological replicate, the difference in the expression of the maternal and paternal allele was computed for each gene (maternal -paternal). The correlation of the maternal and paternal allele expression difference to the expression of published conserved human-mouse cell-type marker genes in the data was computed. The top 1000 most specific marker genes for mouse-human conserved markers of neurons, oligodendrocytes, astrocytes, endothelial cells and microglia were used (McKenzie et al., 2018). We then tallied the number of positively and negatively correlated marker genes for each category and a Chi-square test of the dependence of the allelic expression difference on the cell type category was computed. A mosaic plot was computed using the vcd package in R, which revealed the pearson residuals and associations between cell-types and parental allele expression.

Single Cell RNASeq Data Analysis
Published adult mouse hypothalamus single cell RNASeq data generated by  were downloaded from GEO (gene expression omnibus). We extracted the gene expression data for each cell assigned to each cell-type, according to definitions in the published study. The mean expression level for each gene was computed across all cells for a given type. We then extracted the data for all cell types and imprinted genes, based on imprinted genes identified in our previous study .
Unsupervised hierarchical clustering of imprinted gene cellular expression profiles was performed in R using Euclidean distance and the Ward.D2 method.

Pyrosequencing
For each AQ assay, the percent expression from C57Bl/6J and Castaneous alleles were converted to log ratios; log(B6/Cast). Data was then analysed by two-way ANOVA in R with Cross (levels: F1cb, F1bc) and Cell-type (levels: neuron, non-neuron), and a significant interaction effect Cross X Cell-type validated cell-type specific imprinting in the brain. Cross + Chamber + Cross:Chamber + (1|Mouse)) with data collapsed for levels of Genotype and measured separately for each Sex (Fig. S6A).

Behavioral Analyses with Cross Effects
Cross effects in the data were absorbed in a deep analysis of maternal and paternal allele effects on aspects of social behavior using generalized linear models; and also on analyses of Open Field, Elevated Zero Maze, and Light-Dark Box exploratory and anxiety related behaviors.
Generalized linear models were fit using the glm function in R with factors of Cross (levels: maternal het, paternal het), Paternal Allele (levels: wt, het), and Maternal Allele (levels: wt, het); i.e. glm(Variable ~ Cross + PatAllele +MatAllele). The glm models used the following distributions depending on the type of data: gamma (link="inverse) for latency data, quasipoisson (link="log") for over-dispersed count data, and gaussian (link="identity") for all other measures. The anova function in R was then used to generate Analysis of Deviances Tables with p-values to determine which terms significantly reduced deviance in the models.

Foraging Behavior
Excursion Data Capture Our DeepFeats approach for analyzing modularity in foraging was performed as previously described  with some modifications and advances. In our study, mice were tracked with Noldus Ethovision software. Noldus settings were used to define regions of interest in the foraging arena and indicated when the mouse was in each area. To ensure the tracking is equivalent across different mice, a Procrustes transformation of the XY coordinates was performed to put every tracking file in the same coordinate space. The track coordinates were zero'd to the center of the tunnel to the home cage. We then generated custom code in R to parse the raw Noldus tracking files into discrete, round trip home base excursions from the home cage tunnel. Each excursion is assigned a unique ID key that we call the Concise Idiosyncratic Module Alignment Report (CIMAR) string key. It stores the coordinates of the excursion in the data and the CIMAR string includes metadata regarding the mouse number, excursion number, sex, age, genotype and phase. Next, custom code compares the CIMAR coordinates to the raw Noldus data files and constructs a new dataset that extracts 57 measures from the Noldus output, which we use to initially statistically describe each excursion.
The 57 measures are presented in  and are designed to capture a relatively comprehensive array of different behavioral and locomotor parameters, as well as describe interactions with food and non-food containing patches and exposed regions in the environment. These measures consist of shape, frequency, order and location statistics of an animal's X and Y movements, numbers of visits and time spent at different features in the arena, including food patches (Pots#2 and 4), non-food containing patches (Pots#1 and 3), the tunnel zone, wall zone and center zone of the arena and data describing locomotor patterns, including velocity, gait and distance traveled. The 57 measures for each excursion are scaled (normalized and centered across excursions) because they are in different units.

Identification of a Set of Behavioral Measures to Resolve Modularity in Excursions
This section details the methods to define the set of behavioral measures that best resolve candidate modules. A data matrix was constructed in which the rows are excursions performed by the mice, labeled by CIMAR keys, and the columns are the 57 behavioral measures. A correlation matrix was constructed from the data using the Pearson correlation statistic. The measures were then systematically filtered from the data as detailed in the main text based on different correlation thresholds using the "findCorrelation" function in the caret package in R.
With this approach, the absolute values of pairwise correlations are considered. If two variables have a high correlation, the function looks at the mean absolute correlation of each variable and removes the variable with the largest mean absolute correlation. We systematically threshold the data in r = 0.5+ value increments as shown in the main text to identify the best set of measures to resolve clusters of excursions. At each threshold, the retained measures are used in an unsupervised clustering analysis to define clusters of excursions. We used the Ward.D2 minimum variance method implemented using the "hclust" function in R to perform the clustering and define compact, spherical clusters. We then statistically define discrete excursion clusters from the results using the Dynamic Tree Cut algorithm (Langfelder et al., 2008). This is a powerful approach because it is adaptive to the shape of the dendrogram compared to typical constant height cutoff methods and offers the following advantages: (1) identification of nested clusters; (2) suitable for automation; and (3) can combine the advantages of hierarchical clustering and partitioning around medoids, giving better detection of outliers. We detect clusters using the "hybrid" method and use the DeepSplit parameter set to 4 and the minimum cluster size set to 20. The total number of clusters detected is quantified at each correlation threshold. Conceptually, more relaxed correlation threshold cutoffs could reduce cluster detection by retaining redundant measures that mask important effects from other measures.
On the other hand, thresholds that are too stringent could reduce cluster detection by pruning informative measures. Our objective is to identify the threshold that uncovers the most informative and sensitive set of measures for resolving different clusters of excursions, setting the stage for the discovery of potential modules.

Measures
In our study, Dynamic Tree Cut will deeply cut branches in a dendrogram generating large numbers of small clusters if there are few bona fide relationships in the data. Thus, to test whether bona fide clusters of excursions exist in the data we implemented a random sampling procedure in R in which we randomly sample from the matrix of the retained behavioral measure data to break the relationships between the excursions and the measures. The sampled null data matrix is then subjected to the same clustering and quantification procedure to determine the number of clusters found by Dynamic Tree Cut. A null distribution is created from 10,000 iterations and compared to the observed number of clusters, which is expected to be significantly less than the null due to bona fide biological relationships between the excursions and set of retained measures. A lower tailed p value was computed to test this outcome. In a modification compared our previous study,

IGP Permutation Test to Identify Significant Modules of Behavior
To test whether reproducible modules of behavior exist in the data for the foraging excursions, we use the in-group proportion (IGP) statistical method for testing for reproducible clusters between two datasets (Kapp and Tibshirani, 2007). We built a modified version of this function for parallelized computing to speed the analysis for large number numbers of permutations. The excursion data for the mice is separated into a training data and test data partition for reproducibility testing. A balanced partition was generated according to genotype, sex and phase factors using the "createDataPartition" function in the caret package in R. Unsupervised hierarchical clustering was performed on the Training data partition excursions and clusters were defined using Dynamic Tree Cut. Next, the centroids for each training data cluster were computed as the mean values of the behavioral data for the excursions in the cluster. The training data centroids were then used to compute the IGP statistic for each training data cluster based on the test partition data, thereby evaluating the reproducibility of each cluster.
In a modification compared to our previous study , we created a custom IGP permutation test that is based on a distance calculation, rather than the correlation implementation in the clusterRepro R package. The distance IGP testing framework was written in C++ and speeds the permutation test by many fold and is a more accurate replication of the clustering parameters used in the test data. We used this approach to compute p values for each cluster to determine whether the IGP value is greater than chance. False positives due to multiple testing errors were controlled using the q-value method (Storey, 2002).
Modules are thus defined as significantly reproducible training partition excursion clusters (q < 0.1). Each module detected was assigned an ID number and individual excursions in the data were annotated based on the module they match to. This approach facilitated quantifications of module expression frequency by the mice.

Frequency
To statistically evaluate the factors that significantly affect module expression frequency, we used generalized linear modeling functions implemented in R. The full and nested models tested are detailed in the text for each analysis. The generalized linear modeling is performed using a Poisson or quasipoisson distribution, as indicated in the text. We tested the goodnessof-fit of the full model with a chi-square test of the residual deviance and degrees of freedom in R. If the test result was not significant (p > 0.05), we concluded the model fit the data and proceeded based on Poisson-distributed errors. If the Poisson model was not a good fit, we used a quasipoisson model. The likelihood ratio test comparing the full and nested models was performed using the anova() function in R with the additional option test = "Chisq." Our statistical framework is based on a primary and secondary hypothesis testing strategy to manage multiple testing (Bender and Lange, 2001), as indicated in the main text. Likelihood ratio testing for full versus nested models was performed to define significant main and interaction effects in the data. Individual mice were treated as random replicates of the parental and genetic factors. The following full model was used to test aggregated module expression data for males and females (Fig. 5): glm(expression ∼ cross + paternal allele genotype + maternal allele genotype). The following full model was used to test expression data for grouped modules with similar centroid patterns for males and females ( Fig. 6 and S5): glm(expression ∼ module + cross + paternal allele genotype + maternal allele genotype + module:cross + module:paternal allele genotype + module:maternal allele genotype)

DATA AND CODE AVAILABILITY
All raw data and code are available from the Gregg lab. The raw data will be deposited in a public database following publication in a peer-reviewed journal.