The INMA project is a Spanish population-based mother–child multicentre cohort set up in 2003 in seven areas: Ribera d’Ebre, Menorca, Granada, Asturias, Valencia, Sabadell and Gipuzkoa [15]. This study uses data from the INMA Valencia, Sabadell, Asturias, and Gipuzkoa cohorts. The aim of the INMA project is to study exposure to the most important environmental pollutants in air, water and diet during pregnancy and early life and their effects on child growth and development. Over time, multiple lines of research have been developed, including the effect of pets on child development.

Recruitment procedures and inclusion criteria were described elsewhere [16]. Briefly, mothers were recruited during their first prenatal visit to their reference hospital before week 13 of gestation. The inclusion criteria were at least 16 years of age, 10–13 weeks of gestation, singleton pregnancy, intention to undergo follow-up and delivery at the corresponding reference hospital and no communication impediments.

Baseline participants were collected between November 2003–June 2005 in Valencia (n = 855), June 2004–September 2006 in Sabadell (n = 657), April 2004–June 2007 in Asturias (n = 494), and May 2006–February 2008 in Gipuzkoa (n = 638).

Follow-up visits and sample evolution are described in Fig. 1 for the joint cohorts and in Supplementary Fig. 1 for each cohort separately. Briefly, information was collected in follow-ups from pregnancy, age 1 year, age 4–5 years and age 7–8 years. The follow-up was approved by local institutional ethical review boards (Dirección General de Salud Pública; Centro Superior de Investigación en Salud Pública, Parc de Salut Mar; Regional Clinical Research Ethics Committee, Comité de Investigación del Principado de Asturias; and the Euskadi Clinical Research Ethics Committee, respectively), and the participants provided their written and informed consent to participate. This study was carried out in accordance with the principles of the Declaration of Helsinki.

Follow-up visits and sample evolution

Pet companionship was collected by an interviewer at child ages 1 and 4–5 years for five types of animals: any kind of pet, dog, cat, bird or other animal, including small mammals/rodents (rabbits, rats, hamsters, squirrels, guinea pigs), fish and reptiles/amphibia (turtles, tortoises, geckos, salamanders). For each category, presence in the previous year (no/yes) was reported. This yielded five exposure variables (any kind of pet, dog, cat, bird or other animal). A combined variable for each type of animal was created considering the two follow-ups involved, with the following categories: never/always/only at age 1/only at age 4–5 years.

Mental health in the previous six months was assessed via the Strengths and Difficulties Questionnaire (SDQ) [17]. This is a brief behavioural screening questionnaire with 25 items employed in children aged 2–17 years. It was completed by parents when the children were 7–8 years of age. The parental form presents good psychometric properties, with Chronbach’s alpha of 0.73 in the original [18] version and 0.76 in the Spanish version [19]. The SDQ comprises five subscales: (1) emotional symptoms (worry, fear, nervousness and feeling sad or unhappy); (2) conduct problems (disobedience, tantrums, fighting, lying and stealing); (3) hyperactivity/inattention (restlessness, fidgeting, distractibility, impulsivity and attention span); (4) peer problems (bullying, being picked on, having few friends or being solitary); and (5) prosocial behaviour (kind, helpful, empathic, generous) [20]. Each subscale includes five questions rated on a 3-point Likert scale [not true (0), somewhat true (1), and certainly true (2)], and each subscale can therefore be scored between 0 and 10. Higher scores indicate more behavioural problems. The prosocial subscale reflects strengths rather than difficulties, so we reversed its scoring to align the interpretation direction with the other subscales (i.e., higher scores consistently indicate greater difficulties). In this work, we employed the broadband scales of internalizing and externalizing scales as our main analysis. This can be calculated by summing subscales I + iv and ii + iii, respectively. In addition, we explored further subscales in supplementary analyses.

Family, parental, perinatal and child characteristics were collected through medical records and structured questionnaires at different follow-up visits (pregnancy, birth and ages 1, 4–5, and 7–8 years).

Parental sociodemographic variables were collected separately for each parent during pregnancy and consisted of occupational social class (lower/middle/upper), defined via a Spanish adaptation of the British social class classification [21] and education level (up to primary/secondary/university) defined by the International Standard Classification of Education 1997 [22]. Age, country of birth and parity were collected during pregnancy. Biological sex, preterm birth, small for gestational age and Apgar score, were obtained at birth. Type of breastfeeding and duration of exclusive breastfeeding were recorded at child’s age 1 and main care provider and nursery attendance at the child’s age 2.

Data on parental employment status were collected at the child’s age of 4–5 years. Parental intelligence was assessed with the similarity subtest of the Wechsler Adult Intelligence Scale-III [23], parental use of toxicants (tobacco and alcohol), family structure, number of siblings, type of dwelling and type of area were collected at the child’s age 4–5 years. Clinical conditions such as rhinitis or lifestyle factors such as the duration of physical activity were considered at the age of 7–8 years.

For descriptive analyses, frequencies and percentages were used for categorical variables, while medians and interquartile ranges were used for continuous variables. Chi-square P values were calculated for differences among cohorts. For bivariate analyses, we applied Mann‒Whitney U and Kruskal‒Wallis tests to assess possible relationships between categorical covariates and pet ownership with SDQ internalizing and externalizing scores. For continuous variables, the possible relationships between SDQ scores and covariates were assessed via Spearman’s correlations.

The relationships between pet ownership and behavioural and emotional symptoms were assessed via generalized negative binomial regression models. Cohort and child age and sex were included in all models regardless of their statistical significance. The final models were constructed following a three-step procedure. First, univariate models were implemented with the covariates and SDQ scores. Covariates significantly related to SDQ scores at P values < 0.20 in the likelihood ratio test were retained and included in the multivariate models. Second, multivariate models were implemented and a variable selection method was used. Thus, variables with a P value < 0.10 were selected, resulting in the core models. Third, the exposure variable (pet ownership) was included in each core model, yielding a total of five multivariate models for each outcome. Tukey’s multiple comparisons were applied in each model to determine differences between categories and post hoc analyses with false discovery rates were also developed to correct for multiple testing.

Six sensitivity analyses were performed: the first was to control for sample attrition using the inverse probability weighting (IPW) method [24] (Supplementary Text); the second, included farm animals in the “other animals” exposure variable; the third, excluded preterm children; the fourth analysis, disaggregated by subscale, was to check the robustness of our findings in the broadband scales (internalizing and externalizing problems). In the fifth sensitivity analysis, we adjusted each model for the remaining exposures to control for the residual effect of the ownership of other pets. The sixth analysis consisted of a comparison between pet ownership and pet cohabitation in a subsample (n = 1528) when these data were available (age 4–5 years).

Statistical analyses were performed using the IBM SPSS Statistics package version 26, R and RStudio (versions 4.1.3 and 2022.02.3 + 492, respectively) with the MASS, haven, foreign, ggplot2, mgcv, multcomp and sjPlot packages. Descriptive plots were developed using datawrapper [25]. Flowcharts were designed using draw.io by JGraph Ltd. [26].