Advanced paternal age diversifies individual trajectories of vocalization patterns in neonatal mice.

Infant crying is a communicative behavior impaired in neurodevelopmental disorders (NDDs). Because advanced paternal age is a risk factor for NDDs, we performed computational approaches to evaluate how paternal age affected vocal communication and body weight development in C57BL/6 mouse offspring from young and aged fathers. Analyses of ultrasonic vocalization (USV) consisting of syllables showed that advanced paternal age reduced the number and duration of syllables, altered the syllable composition, and caused lower body weight gain in pups. Pups born to young fathers had convergent vocal characteristics with a rich repertoire, whereas those born to aged fathers exhibited more divergent vocal patterns with limited repertoire. Additional analyses revealed that some pups from aged fathers displayed atypical USV trajectories. Thus, our study indicates that advanced paternal age has a significant effect on offspring's vocal development. Our computational analyses are effective in characterizing altered individual diversity.

Introduction

Human infant crying is an innate form of social communication (Acebo and Thoman, 1992, 1995) that is used to attract attention from caregivers (Soltis, 2004) and affects cognitive control in adults (Dudek et al., 2016). It has been suggested that infant crying can serve as a marker of the behavioral and cognitive development, and the altered crying features may indicate a risk for autism spectrum disorder (ASD) (Esposito and Venuti, 2010; Sheinkopf et al., 2012; Unwin et al., 2017) and other neurodevelopmental disorders (NDDs) (Kivinummi et al., 2020; LaGasse et al., 2005). Although the etiology of NDDs remains unclear, genetic, epigenetic, and environmental factors are likely to play important roles (Moreno-De-Luca and Martin, 2021; Reichard and Zimmer-Bensch, 2021; Scattolin et al., 2021; Schaefer and Mendelsohn, 2008; Tordjman et al., 2014). Recent epidemiological studies have shown that advanced paternal age is associated with the risk of NDDs and the lower body weight in offspring (Alio et al., 2012; Hubert et al., 2011; Hultman et al., 2011; Khandwala et al., 2018; Kimura et al., 2018; Krug et al., 2020; Lundstrom et al., 2010; Reichenberg et al., 2006; Yamamoto, 2021). Thus, we have investigated whether the advanced paternal age affects infant crying and body weight in mice and whether these phenotypes may help to understand NDDs by providing predictive validity in mouse genetic studies.

Ultrasonic vocalizations (USVs) in rodents, particularly in mice, have mostly been investigated in relation to the neurobiology of vocal communication (Mooney, 2020; von Merten et al., 2021). The separation of a pup from its mother and littermates induces the emission of USVs consisting of various sound elements (i.e., syllables), and the USVs can trigger maternal approach and retrieval behavior (D'Amato et al., 2005; Hahn and Lavooy, 2005). Neonatal USVs in mice could model some aspects of human infant cries (Esposito et al., 2017) and, interestingly, exhibit gradual postnatal developmental changes in their acoustic features and compositions (Grimsley et al., 2011). Mouse models for NDDs have been shown to exhibit various differences in USV parameters, including reduced number, higher or lower syllable frequency, and shorter durations (Dougherty et al., 2013; Scattoni et al., 2008; Shu et al., 2005; Wohr, 2014). Previously, we reported that an advanced paternal age alters offspring’s behavioral phenotypes relevant to NDDs in mouse models (Yoshizaki et al., 2016, 2021). It leads to deficiencies in the number, composition, and properties of USVs at postnatal day 6 (P6) (Yoshizaki et al., 2021) in C57BL/6J mice, the most common inbred strain for behavioral and genetic analyses in preclinical studies.

Here, we comprehensively evaluated the developmental trajectory of early vocal communication in postnatal C57BL/6J mice derived from young and aged fathers, using semi-automatic, supervised, and unsupervised computational methods. Our data suggest that an advanced paternal age diversifies the developmental trajectories of vocalization, resulting in an increased proportion of atypical individuals with lower body weight gain.

Results

Offspring were obtained by mating young (3-month-old) female C57BL/6J mice with either young (3-month-old) or aged (20-month-old) male C57BL/6J mice (YFO and AFO, respectively, Figure 1A). To record the USVs, each pup was separated one by one from its mother and littermates. The USV sonograms were transferred from the recorded acoustic waveforms and analyzed by three tools: a semi-automatic procedure (USVSEG) (Tachibana et al., 2020), a supervised machine learning approach (VocalMat) (Fonseca et al., 2021), and an unsupervised modeling algorithm (variational autoencoders, VAE) (Goffinet et al., 2021) (Figure 1B). The USVSEG requests manual syllable classification and noise inspection. To confirm and quantify the USV alterations without a bias from a human perspective, we also performed the VocalMat and VAE.

Figure 1

The experimental design and methods of syllable analyses

(A) Workflow chart of the experimental design. Offspring were obtained by mating young C57BL/6J female mice with either young (3-month-old) or aged (20-month-old) C57BL/6J male mice (YFO and AFO groups, respectively). At P3, P6, P9, and P12, each pup was separated from its mother and littermates for a 5-min ultrasonic vocalization (USV) recording.

(B) The maternal-separation-induced USVs were processed and analyzed using USVSEG, VocalMat, and variational autoencoders (VAE).

Advanced paternal age altered syllable properties and composition during development

The USVSEG revealed that the AFO group emitted fewer USVs (fathers’ age, F (1, 411) = 45.879, p < 0.001; postnatal day, F (3, 411) = 6.508, p < 0.001; body weight, F (1, 411) = 2.065, p = 0.151, fathers’ age × postnatal day, F (3, 411) = 1.811, p = 0.145, mixed model) with shorter durations (fathers’ age, F (1, 411) = 34.787, p < 0.001; postnatal day, F (3, 411) = 21.567, p < 0.001; body weight, F (1, 411) = 16.836, p < 0.001, fathers’ age × postnatal day, F (3, 411) = 1.754, p = 0.155, mixed model) (Figures 2A-2B). No significant difference was found in the maximum frequency or amplitude between the groups (Figures 2C-2D). To further characterize call types, we classified syllables into 12 types (Figure 2E) based on their spectrotemporal patterns (Hori et al., 2020). Overall, the AFO group presented with a reduced number of syllables in three types (mixed model with post hoc comparison Tukey-Kramer test, see Table S1 for results of mixed model). Also, this group of pups had shorter durations in one jump syllables and higher maximum amplitudes in two types of syllables (chevron and wave), but no significant difference was found in the maximum frequency (Figure S1 and Table S1).

Figure 2

Advanced paternal age altered syllable properties and composition during development (data were obtained from USVSEG)

(A) The number of calls in pups sired by young (3-month-old) or aged (20-month-old) male mice (YFO and AFO groups, respectively). ∗∗∗p < 0.001 indicates a significant fixed effect of fathers’ age (mixed model). Data are shown as the mean ± SEM for each group.

(B) Syllable duration in the YFO and AFO groups. ∗∗∗p < 0.001 indicates a significant fixed effect of fathers’ age (mixed model). Data are shown as the mean ± SEM for each group.

(C) The syllable maximum frequencies in the YFO and AFO groups. Data are shown as the mean ± SEM for each group.

(D) The syllable maximum amplitudes in the YFO and AFO groups. Data are shown as the mean ± SEM for each group.

(E) Twelve distinct types of syllables (Hori et al., 2020) were obtained. ∗p < 0.05 and ∗∗∗p < 0.001 indicate a significant difference between the two groups (Tukey-Kramer test, see Table S1 for the results of a mixed model). Data are shown as the mean ± SEM for each group. Scale bar: 50 kHz, 20 ms.

(F) Pie graphs illustrating the percentage composition of the 12 types of syllables from P3 to P12 in the YFO and AFO groups. Colors indicate the syllable types shown in (E). ∗∗p < 0.01 and ∗∗∗p < 0.001 indicate a significant difference between the two groups (multivariate ANOVA following Benjamini-Hochberg correction).

(G) The syllable ratio was calculated as the difference between the number of pink spectrum syllables minus blue spectrum syllables and the number of pink spectrum syllables plus blue spectrum syllables, as follows: . The blue dot indicates that the ratio reached zero. ∗∗p < 0.01 indicates a significant fixed effect of fathers’ age (mixed model). Data are shown as the mean ± SEM for each group.

(H) Developmental trajectories of syllable types in the YFO and AFO groups. ∗∗∗p < 0.001 indicates a significant fixed effect of fathers’ age (mixed model). Data are shown as mean ± SEM for each group.

(I) Entropy scores for syllables in the YFO and AFO groups. ∗∗∗p < 0.001 indicates a significant fixed effect of fathers’ age (mixed model). Data are shown as the mean ± SEM for each group.

(J) Trajectories of the one-note ratio in the YFO and AFO groups. ∗∗p < 0.01 indicates a significant fixed effect of fathers’ age (mixed model). Data are shown as mean ± SEM for each group.

The syllable analysis showed that in the YFO group all types of syllables changed the percentage from P3 to P12: four types of syllables (downward, flat, one jump, and harmonics) decreased the percentages, whereas other types of syllables increased the percentages (Figure S2). We next analyzed the developmental transition of increasing and decreasing syllables; pink and blue spectrum colors from dark to light indicate the syllable percentages increase and decrease from high to low in YFO at P12, respectively (Figure 2F). The syllable types of blue spectrum colors were gradually replaced with those of pink ones from P3 to P12 in both YFO and AFO groups (Figure 2F). Although no significant difference was detected at P3, the syllable compositions significantly differed between the AFO and YFO groups from P6 to P12 (fathers’ age, P6: F (11, 93) = 2.687, p = 0.006; P9: F (11, 93) = 6.266, p < 0.001; P12: F (11, 93) = 4.150, p < 0.001, multivariate ANOVA [MANOVA]) (Figure 2F). To investigate the dynamics of syllable composition during development, we used a syllable ratio to evaluate how quickly pups’ syllables transited from the initial period that the syllable proportion of blue spectrum color was enriched to the later period that the syllable proportions of blue and pink spectrum colors were equal (Figure 2G). At P3, the syllable ratio was almost identical between the YFO and AFO groups. Although the YFO group showed the syllable ratio reaching to zero before P12, the AFO group did not (fathers’ age, F (1, 411) = 10.957, p = 0.001; postnatal day, F (3, 411) = 11.739, p < 0.001; body weight, F (1, 411) = 0.476, p = 0.490, fathers’ age × postnatal day, F (3, 411) = 2.448, p = 0.063, mixed model).

Finally, we addressed the syllable diversity. Although the YFO group emitted USVs that contained an increasing number of syllable types with age, the AFO group produced USVs with a significantly smaller number of syllable types from P3 to P12 (fathers’ age, F (1, 411) = 75.150, p < 0.001; postnatal day, F (3, 411) = 4.104, p = 0.007; body weight, F (1, 411) = 14.985, p < 0.001; fathers’ age × postnatal day, F (3, 411) = 0.845, p = 0.470, mixed model) (Figure 2H). To quantify the syllable properties, entropy score, which ranged from 0 to 1, was calculated to indicate production uniformity (Grimsley et al., 2011) (Figure 2I). The score approached 1 when a pup uniformly (or diversely) produced all the syllable types and approached 0 when a pup preferred to produce only one syllable type (limited repertoire). The entropy score increased continuously in both groups during development. However, the AFO group consistently showed lower entropy scores than the YFO group across all postnatal stages (fathers’ age, F (1, 390.1) = 60.200, p < 0.001; postnatal day, F (3, 384.0) = 2.748, p = 0.043; body weight, F (1, 372.4) = 12.899, p < 0.001; fathers’ age × postnatal day, F (3, 393.8) = 1.841, p = 0.139, mixed model). For reflecting the complexity of syllable usage from a different angle, we calculated the one-note ratio according to a previous study (Klenova et al., 2021). A higher one-note ratio indicated the less complex syllable usage. As we expected, the AFO group constantly showed a higher one-note ratio from P3 to P12 (fathers’ age, F (1, 411) = 10.200, p = 0.002; postnatal day, F (3, 411) = 7.939, p < 0.001; body weight, F (1, 411) = 2.713, p = 0.100; fathers’ age × postnatal day, F (3, 411) = 3.015, p = 0.030, mixed model) (Figure 2J).

The aforementioned data demonstrate that advanced paternal age reduced the number and duration of specific call types, delayed development of syllable compositions, less complex usage, and limited syllable repertoire in offspring. Thus, the AFO group had overall a limited vocalization repertoire.

Advanced paternal age increased postnatal USV variation in offspring

To compare the developmental trajectories among individual mice, we applied principal component analysis (PCA) to extract the syllable features of the USVSEG data (Figure 3A) using five syllable parameters. The loading plot showed correlations between the syllable parameters (i.e., the number of USV, the number of syllable types, duration, maximum frequency, and maximum amplitude) and the two principal components (i.e., PC1 and PC2) (Figure 3B). The number of USV, the number of syllable types, duration, and maximum amplitude positively contributed to the PC1, whereas maximum frequency contributed negatively. The number of USV, the number of syllable types, and maximum frequency positively contributed to the PC2, whereas duration and maximum amplitude contributed negatively. Pups of 90 percentile in each group were grouped using ovals created by JMP software. At P3, two ovals were wide and overlapped in both groups, indicating that individual P3 pups were diverse in vocal repertoire regardless of paternal age. From P6 to P12, the YFO group showed a developmentally convergent pattern with gradually smaller oval areas in YFO (gray ovals in Figure 3A). In contrast, a relatively wider individual difference remained in the AFO group, even at P12 (pink ovals in Figure 3A). This longitudinal analysis clearly indicated the different developmental trajectories of the YFO and AFO groups (fathers’ age, F (1, 418) = 106.864, p < 0.001, postnatal day, F (3, 421) = 20.132, p < 0.001; fathers’ age × postnatal day, F (3, 412) = 2.412, p = 0.066, MANOVA).

Figure 3

Advanced paternal age influenced individual trajectories of vocal development and resulted in an increased number of atypical individuals (data were obtained from USVSEG)

(A) Principal component analysis (PCA) summarized the patterns of individual pups sired by young (3-month-old) or aged (20-month-old) male mice (YFO and AFO groups, respectively) from P3 to P12. Each point represents a syllable pattern of a single pup. Shorter distances between points indicate greater similarity in the syllable patterns. The gray and pink ovals contain 90% of the population of the YFO and AFO groups, respectively. Ovals are based on the two-variate normal distribution generated by JMP software.

(B) The loading plot of the PCA shows correlations between the original ultrasonic vocalization (USV) parameters and two principal components.

(C) Clustering analysis separated individual offspring into five clusters according to PC1 and PC2. Thick black and red dotted lines represent the average of YFO and AFO, respectively.

(D) Based on the cluster analysis according to PC1 and PC2, the pups of Clusters 1–4 shown in graded purple colors were classified as typical individuals, and the pups of Clusters 5 in orange were classified as atypical individuals. Typical and atypical pups are marked on the PCA spatial map.

To better understand the individuality of USVs during development, clustering analyses were performed using Gaussian mixed models (GMMs) along with the Akaike information criterion (AIC) (Table S2). The number of clusters was objectively determined using the minimum AIC. PC1 and PC2 were loaded together to apply clustering analyses and were clustered into five different clusters that shared a common pattern (Figure 3C). The different cluster patterns (i.e., the proportion of individuals in each cluster) were apparent between the YFO and AFO groups; the YFO group was separated into four clusters (see black lines, the thick dotted line represented the average) and enriched in Cluster 1 (55.9%), whereas the AFO group was separated into five clusters (see red lines, the thick dotted line represented the average) and was enriched in Cluster 5 (26.1%, χ2 = 22.993, p < 0.001, chi-squared independence test). A similar phenotype was observed in the additional cluster analyses using other USV parameters. The number of calls was positively correlated with syllable duration (YFO: Pearson’s correlation coefficient r = 0.512, n = 59, p < 0.001; AFO: r = 0.547, n = 46, p < 0.001), and these two factors were clustered and separated into five clusters (Figure S3A). The developmental trajectories of the number of calls and duration observed in the YFO group were distributed among four clusters and enriched in Cluster 1, whereas those in the AFO group were distributed among five clusters and enriched in Cluster 3 to 5. The cluster patterns significantly differed between the YFO and AFO groups (χ2 = 14.810, p = 0.005, chi-squared independence test).

Next, we categorized the developmental trajectories of the number of syllable types into five clusters (Figure S3B), in which the YFO group was again enriched in Cluster 1, whereas the AFO group was enriched in Clusters 3 to 5. There was a significant difference in the cluster patterns between the YFO and AFO groups (χ2 = 29.642, p < 0.001, chi-squared independence test). Finally, to examine the individual developmental patterns of syllable properties, we clustered the entropy scores. Again, the cluster distribution significantly differed between the YFO and AFO groups (χ2 = 31.304, p < 0.001, chi-squared independence test). The YFO group occupied three clusters and was dominant in Cluster 1, whereas the AFO group occupied five clusters and was dominant in Cluster 3 to 5 (Figure S3C). Therefore, the longitudinal developmental patterns significantly differed between the YFO and AFO groups, with some patterns unique only to individuals in the AFO group.

A unique trajectory pattern (i.e., Cluster 5) was observed in the clustering analysis (Figure 3C). Because this cluster only comprised with AFO, we evaluated the pups in Cluster 5 as “atypical” individuals. Conversely, pups in Clusters 1 to 4 were judged to be “typical” ones. These typical and atypical pups were labeled with graded purple and orange colors in the PCA data (Figures 3D and S4). From P3 to P12, the distance among the typical and atypical pups gradually decreased, which showed converged but different patterns (Figures 3D and S5). Although the pups were mixed at P3, the typical and atypical groups gradually diversified as they developed. The convergent phenotype was more apparent in the YFO group because this group contained only typical individuals (Figures S4 and S5). On the other hand, atypical individuals with a unique pattern were originated from 26.1% of the AFO group.

Advanced paternal age altered the postnatal development of syllable composition in pups

To confirm and quantify the alteration of syllable composition without a bias from a human perspective, we next applied a supervised machine learning method, VocalMat, for the automated detection and classification of syllables (Fonseca et al., 2021). After classifying the syllables into 11 categories, VocalMat visualized the three-dimensional probability distribution of the syllables using diffusion maps, which showed a dynamic change from P3 to P12 (Figure S6). Differences in syllable composition between the YFO and AFO groups were also visually verified (Figure 4A). To quantitatively evaluate these differences, the pairwise distance between the centroids of syllable types within each manifold structure was computed as a heatmap (Figure 4B). For a summary and direct comparison of the manifolds between the groups across the four postnatal days, matrix correlations were calculated within and between groups (Figures 4C and 4D). The three-dimensional syllable compositions underwent gradual changes within both groups and the highest similarity was observed between P6 YFO and P9 AFO when comparing the two groups (Figure 4D). These results confirm that the AFO group emitted USVs with different and developmentally delayed syllable compositions.

Figure 4

Advanced paternal age led to a different three-dimensional distribution of syllable composition (data were obtained from VocalMat)

(A) Syllables from pups sired by young (3-month-old) or aged (20-month-old) male mice (YFO and AFO groups, respectively) were mapped in a diffusion Map after classification from P3 to P12.

(B) The pairwise distance between the centroids of syllable types within each manifold structure.

(C) The matrix of Pearson’s correlation within groups. The number indicates the correlation coefficient.

(D) The matrix of Pearson’s correlation between groups. The number indicates the correlation coefficient.

VAE may support the convergent pattern of syllable development

To objectively identify important USV variations, we further applied an unsupervised modeling algorithm, VAE (Goffinet et al., 2021), to characterize subtle changes without a human bias in syllable classification. Each point representing one syllable was mapped to the inferred latent space and visualized (Figure 5A). A unique syllable distribution at each developmental stage in both the YFO and AFO groups seems to show that vocalization patterns developed through a dynamic process, which is consistent with our results from USVSEG and VocalMat. However, it is difficult to quantify and compare these syllable patterns directly between the YFO and AFO groups. In the AFO map, some areas with lower densities were noted. Therefore, we analyzed distributions of syllables for each pup and visualized using maximum mean discrepancy (MMD) (Figure 5B). Lighter colors indicated more similar syllable repertoires. The MMD matrices showed that the syllable repertories became similar between individuals from P3 to P12 in the YFO group, whereas the syllable repertories of pups of AFO remained different from each other. We further performed t-distributed stochastic neighbor embedding (t-SNE) to the VAE data to reduce the dimensions (Figure 5C). Again, the vocal repertoires of the YFO group demonstrated developmental convergence despite different trajectories, whereas those of the AFO group were relatively wide scattered during development.

Figure 5

Variational autoencoders confirmed that the advanced paternal age affected the convergent pattern of syllable development (data were obtained from variational autoencoders)

(A) Variational autoencoders (VAE) mapped syllables from pups sired by young (3-month-old) or aged (20-month-old) male mice (YFO and AFO groups, respectively). Each point represents a UMAP projection of one syllable. The distance between points indicates neighbor similarity, with closer points being more similar.

(B) Similarity matrix between syllable repertoires for each pup from P3 to P12 based on MMD. Lighter values correspond to more similar syllable repertoires (lower MMD).

(C) t-SNE summarized the patterns of individual pups from the YFO and AFO groups from P3 to P12. Each point represents the syllable pattern of one pup. Shorter distances between points indicate greater similarity in syllable patterns. Gray and pink clusters contain 90% of the population of the YFO and AFO groups, respectively.

The consistent results of the USVSEG and VAE findings were as follows: (1) the inferred spatial positions of the YFO and AFO groups were clearly similar at P3 but became significantly different in later developmental stages; and (2) a relatively divergent pattern of syllable development was observed in the AFO group, in contrast to the convergent pattern noted in the YFO group.

Advanced paternal age affected the body weight gain

To explore the implications of vocal communication in the physical development of neonatal mice, we measured the body weight after each USV recording. Lower body weight gain was consistently observed in the AFO group than the body weight gain in the YFO group (fathers’ age, F (1, 412) = 35.888, p < 0.001; postnatal day, F (3, 412) = 1842.759, p < 0.001; fathers’ age × postnatal day, F (3, 412) = 2.196, p = 0.088, mixed model) (Figure 6A). The aforementioned USVSEG data showed the body weight affected the syllable duration, the number of syllable types, and entropy score in a mixed model. To determine whether syllable development was directly associated with body weight gain in each stage, Pearson correlation was applied and visualized as a heatmap (Figure 6B). Unexpectedly, a greater number of syllable parameters was significantly correlated with body weight in the AFO group than the one in the YFO group. Further, in the YFO group, most syllable parameters showed significant correlations with body weight on P6; but in the AFO group, P9 was the equivalent developmental stage.

Figure 6

Advanced paternal age affected the body weight gain

(A) Body weight of pups sired by young (3-month-old) or aged (20-month-old) male mice (YFO and AFO groups, respectively). ∗∗∗p < 0.001 indicates a significant fixed effect of fathers’ age (mixed model). Data are shown as the mean ± SEM for each group.

(B) Pearson’s correlation between body weight and ultrasonic vocalization (USV) parameters that were obtained from USVSEG was calculated and presented from P3 to P12 as a heatmap. ∗p < 0.05 and ∗∗p < 0.01 indicate a significant correlation between the two variables.

(C) Body weight of pups of five clusters according to Figure 3D. ∗∗∗p < 0.001 indicates a significant fixed effect of cluster (mixed model). Different letters (a and b) next to the legends indicate a significant difference (p < 0.05) tested by the Tukey-Kramer test. Data are shown as the mean ± SEM for each group.

(D) Body weight of typical and atypical pups. ∗∗∗p < 0.001 indicates a significant fixed effect of pups’ developmental pattern (mixed model). Data are shown as the mean ± SEM for each group.

Finally, we examined whether the body weight difference was also presented between typical and atypical pups. Different body weight gain was seen among five clusters based on the above clustering analyses of PC1 and PC2 (cluster, F (4, 400) = 7.856, p < 0.001; postnatal day, F (3, 400) = 1400.183, p < 0.001; cluster × postnatal day, F (12, 400) = 0.271, p = 0.993, mixed model) (Figure 6C). The body weight gain of Cluster 5 was significantly different with the body weight gain of Clusters 1, 3, and 4 (Tukey-Kramer test, see Table S3 for the results). A significant lower body weight gain was also demonstrated in the atypical pups compared with the typical pups (pups’ developmental pattern, F (1, 412) = 13.288, p < 0.001; postnatal day, F (3, 412) = 695.246, p < 0.001; pups’ developmental pattern × postnatal day, F (3, 412) = 0.548, p = 0.649, mixed model) (Figure 6D).

Thus, the atypical pups in the AFO group displayed lower body weight gain. The AFO group also showed a delayed postnatal time point, although a specific correlation between body weight and USVs was difficult to evaluate.

Discussion

Using the three sophisticated computational analyses, we report that an advanced paternal age diversifies the developmental trajectories of USVs in mice at early postnatal stages. One of our main findings was that the advanced paternal age causes alterations in early vocal behavior and increases the number of offspring with atypical developmental patterns. This phenotype recapitulates clinical evidence that advanced paternal age is a risk factor for the atypical development observed in children with NDDs (Hultman et al., 2011; Lundstrom et al., 2010; Reichenberg et al., 2006) and suggests that the effect of advanced paternal age could be detected in early infancy.

Previous studies reported that USVs from mouse pups might be analogous to crying in human infants and are thus one of the few methods used to investigate behavioral development during the early postnatal period (Branchi et al., 2001; Scattoni et al., 2009). In our study, the YFO group emitted various types of syllables during the postnatal period, and this finding is consistent with that of a previous study demonstrating normal syllable development in CBA/CaJ mouse pups (Grimsley et al., 2011). Compared with the YFO group, the AFO group emitted a narrower spectrum of syllable types with a reduced number of syllables, which is reminiscent of the poorer vocal repertoire in some ASD patients (Patten et al., 2014). A reduction in syllable repertoire has also been described in genetic NDD models, such as Cd157 KO mice (Lopatina et al., 2017) and Tbx1 heterozygous mice (Hiramoto et al., 2011; Takahashi et al., 2016). In addition, we observed an altered composition of syllables in the AFO group, which has previously been reported in other models, such as the Reelin mutant (Romano et al., 2013), fmr1 knockout (Lai et al., 2014), ScSn-Dmd/J mutant (Miranda et al., 2015), and Tbx1 heterozygous (Hiramoto et al., 2011; Takahashi et al., 2016) mice. Thus, our model of advanced paternal age showing impairments in syllable properties and composition also fits the NDD-like phenotype.

There is a consensus that early childhood development typically follows a series of developmental milestones within a similar age range. However, 10%–17.8% of children demonstrate mild to severe developmental delays, including NDDs, where delayed and/or atypical development has been considered an early warning sign (Boyle et al., 1994; Rice et al., 2014; Rosenberg et al., 2008; Zablotsky et al., 2019). Our mouse model revealed that individual pups in the AFO group showed deviated patterns becoming most significantly different at P12, whereas the YFO group exhibited typical development of vocal repertoires gradually converging with age. In addition, 26.1% of the AFO group displayed a unique developmental pattern and was identified as atypical individuals. This phenotype suggests that the impact of advanced paternal age is detectable in offspring during an early neonatal period in mice. We notice that the unique pattern of AFO was not only displayed in the clustering analysis based on PCA data (Figure 3C) but also was shown in other cluster analyses based on syllable parameters (Figure S3). These results may elaborate the heterogeneity of NDDs and mirror children with NDDs showing a variety of atypical developmental behaviors (American Psychiatric Association, 2013; Grzadzinski et al., 2013; Wozniak et al., 2017), including nonuniformly impairment that becomes more distinct with age (Ozonoff et al., 2010; Piven et al., 2017; Zwaigenbaum et al., 2005). It would be interesting to elucidate the underlying mechanisms and determine whether those atypical pups show behavioral deficiencies in the adulthood stage in future studies. A limitation here is that no standard criteria can be used to determine the typical and atypical individuals in neonatal USVs. A previous study showed that the one-note ratio could reflect the complexity of syllable usage, which may provide a potential to identify the atypical USV pattern (Klenova et al., 2021). As we expected, the atypical pups indeed demonstrated a higher one-note ratio (i.e., less complex syllable usage) (Figure S7). Therefore, more evidence is needed to establish a uniform standard. At this moment, we do not rule out the possibility that atypical individuals also exist in the YFO group.

Because sensitivity to behavioral signs is low during early development and the phenotypic characteristics of NDDs are complex and heterogeneous in nature, early diagnosis of NDDs is challenging at the clinical level. Recently, dimensional and/or clustering approaches have been applied to identify developmental changes and to screen individuals in general or at-risk populations (Astle et al., 2021). Moreover, a mouse model of 16p11.2 copy number variation, a risk for NDDs, has shown that dimensional features of peripubertal social behaviors can be computationally predicted by developmental trajectories of neonatal USVs (Nakamura et al., 2021). In the present study, PCA and clustering analysis were integrated to extract the complicated features of USVs and identify individual pups with atypical USV patterns. These approaches may provide a translational perspective to capture the underlying heterogeneity of NDDs in the population at large and offer valuable suggestions for early diagnosis and intervention in children with NDDs.

Epidemiological studies have consistently demonstrated that the advanced paternal age could be a cause of smaller fetuses in gestational age, as well as lower body weight (Alio et al., 2012; Khandwala et al., 2018). In humans, the association between physical growth, including body weight gain, and neurodevelopmental delay has been examined in general populations of term-born infants; it was found that infants with poor growth in early infancy are at an increased risk of neurodevelopmental delay (Sanefuji et al., 2021; Skuse et al., 1994). We noticed that during the first two postnatal weeks, the AFO group exhibited lower body weight gain, which was also observed in a rat model of advanced paternal age (Krug et al., 2020). We observed a similar tendency in our previous study, in which we only recorded USV and body weight at P6 (Yoshizaki et al., 2021). The lower body weight gain was consistently shown in the atypical pups. In other ASD genetic models such as Clstn2-KO pups, abnormal physical growth (i.e., lower body weight and the increased head-to-body ratio) was also reported (Klenova et al., 2021). Therefore, our model exhibited a similar physical alteration to other mouse models of advanced paternal age and ASD.

There are two possible explanations for the lower gain of body weight in the AFO group. First, the advanced age of a father could directly cause a developmental delay in the physical conditions of the offspring. However, we observed a significant developmental delay in USV parameters (Figures 2 and 3), and three USV parameters (i.e., duration, the number of syllable types, and entropy score) were affected by body weight but robust deviation in individual difference or syllable repertoire complexity (Figures 4 and 5) in the AFO group. The group differences in physical developmental delays, therefore, could not simply reflect the impairments of USV development. The transgenerational influence of advanced paternal age is complex. Our previous study found that hypo-methylation of sperm DNA can be a key molecular feature modulating neurodevelopmental programs in offspring (Yoshizaki et al., 2021). The genetic and/or epigenetic factors associated with the advanced paternal age might trigger lower body weight gain affecting the other factors irreversibly during postnatal development. Another possibility is that body weight gain is influenced by altered maternal behaviors due to differences in the communication of the offspring, including USVs. It is well known that pup USVs have communicative significance (Ehret and Bernecker, 1986; Hashimoto et al., 2001; Kikusui and Hiroi, 2017; Liu et al., 2003, 2006), which is crucial for maternal care behaviors, preserves the social bonds between mother and infant, and is essential for the healthy development of offspring (Mogi et al., 2011, 2017). Previous studies have shown that USVs of Tbx1 heterozygous pups do not induce maternal approach behaviors (Kato et al., 2021; Takahashi et al., 2016). More work is needed to further determine whether USVs emitted by AFO elicit less maternal approach behavior and thus lower body weight. Although underlying mechanisms are not fully understood, an advanced paternal age affects the physical condition and neurodevelopment of pups, possibly via a combination of genetic, epigenetic, and environmental factors.

Although approaches to analyzing USV have advanced dramatically during the past decade, it is undeniable that any method has its limitations. The semi-automatic approach (USVSEG) provides a reliable segmentation performance for USV syllables (Tachibana et al., 2020), whereas the manual syllable classification and noise inspection were time-consuming and affected by human bias. VocalMat, a supervised machine learning method, is able to automatically detect and classify USVs (Fonseca et al., 2021) but still has room to improve its accuracy. On the other hand, an unsupervised machine learning algorithm (VAE) can capture data variability but must be trained on each dataset, with the limitation that mouse USV syllables cannot be categorized into distinct types (Goffinet et al., 2021).

In the current study, we applied a composite of these three approaches, understanding the strengths and limitations of each, and collected reliable data from different angles. USVSEG demonstrated the detailed alterations of acoustical features and composition that were obtained after a half year’s work (Figures 2 and 3). VocalMat contributed to the high-throughput detecting difference and similarity of syllable repertoire within days (Figures 4 and S6), although it may not cover all the syllable categories. VAE explored the developmental transition of syllable variability and individual difference with minimizing human bias, although the bulk of syllable difference was relatively ambiguous between groups (Figure 5). Even though these different approaches do not share the same computer algorithms, we were able to reach a unanimous conclusion from these comprehensive USV analyses; USVs in the AFO group showed altered composition and different developmental trajectories. Technology advances in extracting information from USVs can solve shortcomings of existing methods and, most importantly, unify the criteria to identify translational relevant traits.

Our findings provide evidence for the hypothesis that the advanced paternal age expands distinct individuals reflecting “neurodiversity” in offspring. Moreover, as in modern societies in which the age when individuals give birth is increasing, advanced paternal age may represent a risk factor for neurodevelopmental disorders (Centers for Disease Control and Prevention, 2021). In these cases, epigenetic changes are the most likely mechanisms to underly this phenomenon (Osumi and Tatehana, 2021). Another possibility is that advanced paternal age causes the accumulation of spontaneous de novo mutations (Kong et al., 2012; O'Roak et al., 2012), possibly leading to neurodevelopmental disorders. Last but not least, increased longevity occurred during human evolution. The investigation of sophisticated developmental traits through epigenetic inheritance is necessary to account for the process of civilization.

Limitations of the study

Our study has two limitations: first, additional experiments are required to determine whether atypical pups exhibit behavioral deficits in communication and social interaction in adulthood. Second, to thoroughly analyze physical growth, it would be advantageous to include additional biological developmental data in addition to body weight.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact Prof. Noriko Osumi (osumi@med.tohoku.ac.jp).

Materials availability

This study did not generate new unique reagents.

Experimental model and subject details

20-months-old male C57BL/6J mice were used as aged fathers because 18–24 months of age in mice correlates to 56–69 years of age in humans (Flurkey et al., 2007), which meets the definition of “aged”. Eight young (3-month-old) and seven aged male mice were crossed with young (3-month-old) virgin female C57BL/6J mice. After mating, each female mouse was separated from the male mouse and isolated to minimize possible confounding factors related to offspring behavior. In this study, 59 offspring from young fathers and 46 offspring from aged fathers were used. Our exploratory preliminary analysis did not detect a statistically significant effect for the sex. This result is consistent with previous studies (Rieger and Dougherty, 2016; Scattoni et al., 2008). Therefore, the sex effect was excluded from our analysis, and both male and female pups were used. Offspring that died during the experimental period were excluded from this study (mortality rate, YFO group: 4.8% [3/62] vs. AFO group: 8% [4/50]). At P3, each offspring was tattooed using the Aramis Animal Microtattoo System (Natsume Co., Ltd., Tokyo, Japan) for individual recognition after the USV test (described below). The average litter size and the number of pups did not significantly differ between the YFO (7.75 ± 1.16) and AFO (7.14 ± 1.86) groups (t = 0.744, p = 0.474, t-test).

All animals were housed in standard cages in a temperature- and humidity-controlled room with a 12-h light/dark cycle (lights on at 08:00) and had free access to standard laboratory chow and tap water. All experimental procedures were approved by the Ethics Committee for Animal Experiments of Tohoku University Graduate School of Medicine (#2014-112), and the animals were treated according to the National Institutes of Health Guidance for the Care and Use of Laboratory Animals.

Method details

USV recordings

According to previously described protocols (Shu et al., 2005; Yoshizaki et al., 2016, 2017, 2021), each pup was separated from its mother and littermates, placed on a transparent plastic dish with woodchip bedding, and transferred to a sound-attenuating chamber for USV recordings at P3, P6, P9, and P12. An ultrasound microphone (Avisoft-Bioacoustics CM16/CMPA) was inserted through a hole in the middle of the cover of the chamber, approximately 10 cm above the offspring, to record vocalizations. The recorded vocalizations were transferred to an UltraSound Gate 416H detector set (Avisoft Bioacoustics, Germany) at 20–125 kHz. After a 5-min recording session, the body weight of the pups was measured, and they were returned to their home cage. This procedure was repeated until all pups had been recorded. The room temperature was maintained at 22°C.

Semi-automatic analysis of syllable properties

Acoustic waveforms were processed using USVSEG, a GUI-based MATLAB (MathWorks Inc., MA, USA) script, originally developed for segmenting USVs emitted by rodents (Hori et al., 2020; Tachibana et al., 2020). Briefly, the script computed the spectrograms from each waveform (60 s/block), applied a threshold to remove noise from the signal and detected syllables within a frequency range of 20–120 kHz. The segmentation criteria for identifying syllables were as follows: a minimum gap of 10 ms should separate two syllables, and each syllable should have a minimum duration of 2 ms. This script segmented each syllable and exported them as individual jpeg files. The duration, peak frequency (frequency at maximum amplitude), and maximum amplitude of each syllable were calculated automatically. Based on previously published criteria (Hori et al., 2020), by visual inspecting jpeg files, segmented syllables were manually classified into 12 types or excluded as noises (false positive). The pink/blue syllable ratio was calculated as the difference between the number of pink spectrum syllables minus blue spectrum syllables and the number of pink spectrum syllables plus blue spectrum syllables, as follows:

Based on a previous study (Klenova et al., 2021), seven types of syllables (i.e., downward, flat, chevron, upward, complex, wave, and short) were classified into the one-note syllable. The one-note ratio was calculated as follows:

The distance between the typical and atypical groups in Figure S5 was calculated using PC1 and PC2 data as follows:

Calculated the centroid (x, y) of each group using arithmetic mean as follows:

The radius of the circles in Figure S5 was calculated from the average Euclidean distance of each pup to the centroid.

Supervised automatic syllable classification

The supervised learning software VocalMat was used for the automatic detection and classification of syllables. Following the fixed threshold and parameters of the originally published method for syllable segmentation, syllables were classified into 12 categories (11 syllable types + noise) (Fonseca et al., 2021). VocalMat then visualized the syllable distribution using diffusion maps to reduce the 11 dimensions to three dimensions. The pairwise distance between the centroids of syllable types within each group at each stage was calculated and output as matrices using MATLAB.

Unsupervised modeling and characterization of vocalizations

The VAE, which is an unsupervised learning method for extracting essential features without any instructional input, was used to characterize the recorded vocalizations. We used the VAE to reduce the dimensions of the raw data and quantify subtle changes in behavioral variability to preserve as much information as possible, according to an established procedure (Goffinet et al., 2021). The VAE algorithm, which was implemented in PyTorch (v1.1.0), was trained to maximize the evidence lower bound as the standard setting. Because the audio segmentation method implemented in the original program detected large amounts of noise and was unable to detect some subtle syllables, we decided to use USVSEG for audio segmentation prior to using VAE. The syllables were detected and segmented by USVSEG using the same parameters as those mentioned above. After training, a latent map of all syllables was created using a uniform manifold approximation projection. To visualize the differences between individuals, t-SNE was applied to reduce the dimensions for each pup based on the maximum mean discrepancy (MMD) (Goffinet et al., 2021). We excluded pups that emitted two or fewer calls (two in the AFO group at P3, three in the AFO group at P6, and three in the AFO group at P9) from the VAE analysis because the MMD cannot be applied in such datasets.

Quantification and statistical analysis

A mixed model with the fixed effects (father’s age, postnatal day, and body weight), random effects (litter), and repeated measures was applied to examine the statistical significance of overall syllable parameters, syllable ratio, number of syllable types, and entropy. To detect the statistical significance of the number, duration, frequency, and amplitude of 12 distinct syllables, the mixed model with the fixed effects (father’s age, postnatal day, and syllable type) was performed. Tukey-Kramer test post hoc test was applied if the interaction was significant. The significance of body weight was determined by a mixed model with the fixed effects (father’s age or cluster, and postnatal day), random effects (litter), and repeated measures. MANOVA was performed to detect differences in syllable composition on each postnatal day with the effect of fathers’ age. The differences in principal components and t-SNEs between the YFO and AFO groups were detected using MANOVA with the father’s age and postnatal day as two independent effects. The post-hoc comparison was performed using the F test if the interaction was significant. The Pearson’s correlation coefficient was used to determine the correlation between body weight and USV parameters.

The entropy score was calculated using the information theory toolbox in MATLAB. We excluded some offspring (one from the YFO group and one from the AFO group at P3, three from the AFO group at P6, and two from the AFO group at P9) whose total number of syllables and number of syllable types produced in 5 min was insufficient (three or fewer syllables and only one syllable type) from the entropy calculation. To understand the development at the individual level, clustering analysis with GMMs was applied, where the data dimension of eight corresponded to the number and duration of syllables at the four time points. We fitted GMMs with diagonal covariance Gaussian components using the MATLAB function fitgmdist. The number of clusters was determined by minimizing the AIC (Konishi and Kitagawa, 2008; McLachlan and Peel, 2000). Based on the fitted GMM, each individual mouse pup was categorized into a cluster with the maximum posterior probability. The chi-square independence test was subsequently applied to determine whether the cluster distribution was significantly different between the two groups. A PCA was performed to objectively characterize the typical syllable patterns of individual offspring. In the present study, syllable data, including the syllable number, number of syllable types, duration, maximum frequency, and maximum amplitude, were inputted to the PCA to generate principal components. For all comparisons, the significance level was set at α = 0.05. JMP13 Pro software (SAS Institute, Cary, NC, USA) was used for statistical analyses. Data are presented as mean ± standard error of the mean for each group.

All statistical details can be found in the Results section. For comparisons, three asterisks (∗∗∗) represent p values of less than 0.001; two asterisks (∗∗) represent p values of less than 0.01; one asterisks (∗) represent p values of less than 0.05.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Deposited data
USV raw data and curated data obtained from USVSEG	OSF	https://doi.org/10.17605/OSF.IO/PFVTS
Experimental models: Organisms/strains
C57BL/6J mice	Charles River Laboratories Japan (The Jackson Laboratory Japan)	RRID:IMSR_JAX:000664; MGI: 3028467
Software and algorithms
Avisoft RECORDER	Avisoft Bioacoustics	http://www.avisoft.com/recorder/
MATLAB 2021a, 2021b	MathWorks	https://jp.mathworks.com/products/matlab.html
Python 3.7.6	Python Software Foundation	https://www.python.org/downloads/release/python-376/
JMP Pro 15	JMP Statistical Discovery LLC	https://www.jmp.com/ja_jp/home.html