Infant gut anellovirus composition associates with early life factors and childhood atopic disease

Domain: neuroscience
Published: 2025-12-12
Editorial outcome: major_revision

Abstract
Anelloviruses are single-stranded circular DNA viruses that infect eukaryotic cells. They are acquired in infancy and can be found in most body compartments. The Anelloviridae family is highly prevalent and genetically diverse, but factors shaping the composition and the relationship with disease remain elusive. Here, we characterise faecal viromes of 647 one-year-old children from the deeply phenotyped Copenhagen Prospective Studies on Asthma in Childhood (COPSAC 2010) mother-child cohort. In this study, 111 unique anellovirus species are identified and reveal immense prevalence and individual variability. The anellovirus and bacteriophage compositions are slightly correlated, while no correlation is observed with the bacterial composition. Early-life exposures such as attending daycare and having older siblings are associated with a reduced anellovirus abundance and diversity, and influence the anellovirus composition. Additionally, genetic secretors of the FUT2 enzyme have a higher anellovirus abundance compared to non-secretors. Minor associations are observed between anelloviruses and later allergic rhinitis, but no associations are found for concurrent acute infections or the development of other atopic traits. Our findings indicate that anelloviruses are a commensal component of all infants’ faecal virome, and further investigations into how immune-training factors influence their composition longitudinally are pertinent.

Introduction
The gut microbiome is composed of millions of microorganisms, including bacteria and viruses. It is established in infancy and is influenced by both intrinsic and extrinsic factors, such as delivery mode, nutrition, host genetics and environmental exposures. An increasing body of research is finding links between the gut microbiome and acute and chronic diseases. Current research is largely focused on the bacterial component of the gut microbiome and, to some extent bacteriophages. However, the members of the highly prevalent Anelloviridae family infect human cells and showcase a considerably high abundance in some individuals. Despite this, it is still largely unknown how these viruses colonise the host, which factors influence this colonisation, and how they may influence health and disease in childhood.

Anelloviruses are single-stranded circular DNA viruses that infect eukaryotic cells and exhibit an immense genetic diversity. They are highly prevalent, believed to be present in nearly all humans across the earth, and have been found in a multitude of sample types ranging from skin, saliva and blood, to faeces. Coinfections with anelloviruses of multiple different genotypes in one individual also commonly appear, especially in early life.

Anelloviruses are acquired in infancy, potentially already through transplacental transmission, considering that anelloviruses have been identified in umbilical cord blood, amniotic fluid and newborns, though the findings are not entirely concordant. Mothers and their infants shared some, but not all, anelloviruses in the blood. Similarly, adult twins share less faecal anelloviruses than child twins. Several potential sources of transmission have been suggested, such as breast milk or contact with skin, hair or saliva from infected individuals, though there is conflicting evidence. Moreover, birth mode may influence the anellovirus composition, as one study found that children born by spontaneous vaginal delivery exhibited larger gut anelloviral diversity than children born by emergency Caesarean section. Anelloviral richness in blood and faeces has been shown to be low during the first few weeks and months of life, possibly because unknown factors in mother’s milk protect against infection. Subsequently, with weaning, both richness and abundance increase between 6 to 15 months of age. Thereafter, it appears that the faecal anelloviral load decreases, demonstrated by a lower faecal anelloviral richness in mothers compared to infants, which could reflect infant immunity development that is partly dependent on maternal antibodies and the host’s own immune maturation.

Due to their high prevalence, anelloviruses are not believed to be causal in acute or chronic disease development, but this remains an intense topic of investigation. Furthermore, it has been speculated that direct interactions with the developing immune system in childhood could potentially affect the risk of developing immune disorders such as asthma and allergy. On one hand, anelloviruses from several sample types have reportedly been associated with a range of diseases such as IBD (faeces), fever (plasma), cancer (cancer tissue) and asthma (airways, plasma). On the other hand, another study reported no association of intestinal anelloviruses with type 1 diabetes development or seroconversion. In summary, there is a lack of concrete evidence that anelloviruses are causal in disease, and it remains unknown whether certain anelloviral species have more pathogenic potential than others.

In this study, we characterise the infant gut anellovirus composition in 647 children from the deeply phenotyped Copenhagen Prospective Studies on Asthma in Childhood (COPSAC 2010) cohort at age 1 year. The faecal viromes were previously sequenced and meticulously catalogued and described, identifying hundreds of novel viral families. The concurrent bacterial and bacteriophage compositions have previously been associated with asthma development. With this study, we aim to describe the composition and intrinsic and extrinsic shaping factors of the anelloviruses in early life and its association with infection susceptibility and development of the atopic diseases asthma, eczema and allergy.

Results
Characterisation of the early-life virome
We analysed the sequencing data from DNA viruses extracted from faecal samples of 647 1-year-old children in the COPSAC 2010 cohort (Table 1). The early-life faecal virome consisted mainly of bacteriophages, but also contained two families of eukaryotic DNA viruses, namely the adenoviruses and anelloviruses. We have previously presented a deep characterisation of the bacteriophageome and its relation to asthma, as well as a report on the associations of environmental exposures on the overall gut virome in the same cohort. Even though the bacteriophages accounted for the vast majority of relative abundance in most samples (median = 97%, range = 7−100%), the anelloviruses were the most abundant eukaryotic viral family (median = 2%, range = 0−93%), while adenoviruses dominated in just a few individuals (median = 0%, range = 0−76%) (Fig. 1A). More than 9000 unique species of bacteriophages were identified, whereas the anelloviruses consisted of 111 unique species (970 viral OTUs or vOTUs) and the adenoviruses only comprised five unique species (Fig. 1B).

Infant gut anelloviruses are taxonomically diverse
vOTUs were delineated by virome assembly and deduplication at 95% nucleotide identity, where after anelloviral vOTUs were defined based on encoded protein similarity. After aligning the resulting 970 anelloviral vOTUs against anellovirus type species from the ICTV (International Committee for Taxonomy of Viruses website) and constructing an average nucleotide identity (ANI)-based tree (Fig. 2), we found that most (711, 73%) of the vOTUs belonged to the genus betatorquevirus, while alpha- and gammatorquevirus comprised a smaller fraction of the vOTUs. Furthermore, we found that 22 vOTUs formed two genus-level clades of their own that were related to, but separate from alpha- and gammatorquevirus respectively. Here, we propose naming these novel anelloviral genera phitorquevirus and tautorquevirus. They were present in 39% and 49% of the children, respectively (Supplementary Fig. 1A). The 970 anelloviral vOTUs were further grouped into species, using the recent ICTV-ratified anellovirus species cutoff at >35% average nucleotide divergence (100-ANI = AND). This yielded 111 distinct species. In order to bridge the large gap between the more than 100 species and the few genera, we defined anellovirus subgenera by cutting the tree at a depth corresponding to 50%, resulting in 36 anelloviral subgenera, which were named according to the genus to which they belonged. In further analyses, we calculated alpha and beta diversity at the species level; however, focused taxa-specific analyses at the subgenus and genus level for improved interpretability and to reduce the multiple testing burden.

Anelloviruses do not associate with common childhood infections
We tested the associations of acute infections (gastroenteritis, cold, and fever) on the total relative abundance of the bacteriophages, adenoviruses and anelloviruses in the COPSAC 2010 cohort. Participants with reported gastroenteritis symptoms within 7 days of faecal sampling had a significantly higher total adenovirus relative abundance (log2 fold change = 2.4, p < 0.047, Wilcoxon rank-sum test, Fig. 1C). The total relative abundance of anelloviruses and bacteriophages were not related to any acute symptoms (p > 0.05, Wilcoxon rank-sum test, Fig. 1C). Since the role of adenoviruses on infections is already well established we excluded them from further analysis, and focused on the unexplored anelloviruses.

Anelloviruses are uniquely yet highly prevalent in the 1-year-old infant gut
As we demonstrated previously for bacteriophages, the anelloviruses were highly unique for each individual at the vOTU level, with 80% of the anellovirus vOTUs characterised as rare (present in 1−20% of participants), and only four vOTUs (0.4%) were part of the core anelloviruses (present in >50% of participants) (Fig. 1D). However, at the species, subgenus and genus level, the anelloviruses were increasingly shared among participants, and all samples contained at least one anellovirus. The anellovirus subgenera ranged in prevalence from 4% of participants (subgenus gamma 3 virus), to the subgenus betatorquevirus 1 (beta 1) being present in 100% of participants (Fig. 3A). The prevalence, median relative abundance and species richness of the anelloviral subgenera and genera were strongly associated (Fig. 3A and Supplementary Fig. 1A). However there were exceptions, such as the subgenus tautorquevirus 1 (tau 1) with a median relative abundance of 0, a prevalence of 49% and consisting of only 1 observed species (Fig. 3A). Interestingly, there did not appear to be a competitive (inverse) relationship among most anelloviruses, and there were some subgenera that were highly positively correlated (Fig. 3B and Supplementary Data 1). The alpha, beta, and gamma genera were all present in 100% of the samples, while the two previously undescribed genera we identified in this study, tautorquevirus and phitorquevirus, were present in 49 and 39% of the children, respectively (Supplementary Fig. 1A). Similar to the subgenera, the anelloviral genera showed largely positive correlations, in particular, the beta and alpha genera were highly correlated with each other, but also gamma and beta were significantly correlated (Supplementary Fig. 1B and Supplementary Data 1). Of the 111 anellovirus species detected in the cohort, each child carried on average 43 unique species (range = 9−86). The most prevalent anellovirus species was the alphatorquevirus 3, which was present in 631 (98%) of the children, whereas the least common species was the phitorquevirus 57, which was only present in 11 (2%) of the children in our cohort. Furthermore, we found that the anellovirus composition was significantly associated with the bacteriophageome (Procrustes, Supplementary Fig. 2A), but it was not associated with the bacterial composition as determined by 16S rRNA gene amplicon sequencing (Supplementary Fig. 2B).

Extrinsic and intrinsic factors associated with the anellovirus alpha diversity and abundance
How anelloviruses colonise the host in early life and what factors influence the anellovirus composition is still not fully known. To shed light on this, we analysed the total anellovirus relative abundance and Faith’s phylogenetic diversity (PD) in relation to several social, environmental, genetic, pre-, peri- and postnatal factors. Using a generalised linear model adjusted for technical covariates, we found several significant links (Fig. 4A and Supplementary Data 2). The pregnant mothers in the COPSAC 2010 cohort took part in a randomised placebo-controlled fish oil intervention, which has been described previously. In the current study, we found that children born to mothers who received the fish oil intervention during pregnancy had a higher total relative anelloviral abundance (model estimate = 1.18 [1.03, 1.35], p = 0.05) and an increased Faith’s PD at 1 year of age (1.16 [1.00, 1.35], p = 0.02). There were no significant associations between high-dose Vitamin D intervention, delivery mode, or breastfeeding duration and anellovirus diversity or abundance. However, samples collected in the summer had a lower total anellovirus abundance and Faith’s PD compared to samples collected in other seasons (abundance: model estimate = 0.78, CI = [0.64, 0.96], p = 0.02, Faith’s PD: 0.80 [0.67, 0.97], p = 0.02). Similarly, children who had started daycare before faecal sampling had a lower anellovirus abundance and alpha diversity (abundance: 0.75 [0.62, 0.90], p = 0.002, Faith’s PD: 0.61 [0.52, 0.72], p = 1e-08), and the number of older siblings associated with a lower anellovirus abundance (abundance: 0.92 [0.86, 1.00], p = 0.04).

The FUT2 gene encodes the fucosyltransferase 2 enzyme, which determines a person’s histo-blood group antigen secretor status. Noroviruses and rotaviruses have been shown to bind to secreted histo-blood group antigens in the gut mucosa, and non-secretors demonstrate a decreased susceptibility to infection, thus we hypothesised that there may be a similar relationship for anelloviruses. Analysis of secretor status revealed that children who were secretors (GA or GG allele, N = 475) based on the SNP rs601338 had a significantly higher total relative anelloviral abundance compared to non-secretors (AA allele, N = 101) (1.28 per SD increase in abundance, [1.04, 1.59], p = 0.02), though there was no effect on Faith’s PD.

Extrinsic and intrinsic factors associated with the anellovirus beta diversity
When comparing all the previously mentioned extrinsic and intrinsic parameters against the overall anellovirus composition using several beta diversity metrics, we again identified significant associations (Fig. 4B and Supplementary Data 3). Some perinatal factors, such as the mother’s antibiotic use at birth (Aitchison distance R2 = 0.2, p = 0.02) and whether the vaginal birth was induced (Aitchison distance R2 = 0.3, p = 0.002, Jaccard distance R2 = 0.4, p = 0.006) were significantly associated with the overall anellovirus composition. Concordant with the results on alpha diversity and overall relative anelloviral abundance, postnatal factors that may be involved in immune system development, such as whether children had started daycare before the faecal sampling, is associated with a significantly different anellovirus composition compared to children who had not yet started daycare (Bray-Curtis dissimilarity and Jaccard and Aitchison distance R2 = 0.3, 0.5 and 0.2, respectively, all p < 0.05).

Extrinsic and intrinsic factors associated with the anellovirus taxa abundances
Finally, we associated the relative abundance of anelloviral subgenera and genera with all previously mentioned extrinsic and intrinsic parameters to investigate taxa-specific associations. Here we again found several strong signals in the same direction as for the analyses on overall anelloviral abundance (Fig. 4A). Particularly, being in daycare and the number of older siblings were negatively associated with anellovirus abundance across many subgenera, as were several birth-related factors such as sample date, season and induction (Fig. 5 and Supplementary Data 4, 5).

The anellovirus composition is linked to the gut bacterial maturation
Due to the numerous findings relating to the decreased anellovirus abundance and diversity with daycare start, we employed additional analyses to elucidate this relationship. We found that children who had entered daycare had a higher incidence of cold and fever (Chi-squared test, OR 1.63 [1.09−2.47], p = 0.02 and 2.02 [1.27−3.31], p = 0.003), though not gastroenteritis (1.25 [0.67−2.53], p = 0.48). Thus, we considered whether the lower abundance of anellovirus in children in daycare was a compositional decrease driven by an increased abundance of other eukaryotic viruses, such as adenovirus. However, we did not find an association between daycare attendance and the overall adenovirus abundance. We have previously calculated microbiota-by-age z-scores that capture the maturation of the bacterial component of the gut microbiome of the children in the COPSAC 2010 cohort. There was a slight negative correlation between the gut microbial maturation score and the overall relative anelloviral abundance (Spearman’s correlation rho = −0.08, p = 0.04).

The gut anellovirus composition does not associate with the development of childhood atopic disease
While the anellovirus composition was not associated with infections in the COPSAC 2010 cohort, previous studies have suggested potential links to chronic diseases and immune system development. Thus, we investigated the relationship between the anelloviruses and the development of asthma, allergic rhinitis or atopic dermatitis by age 6 years. The total relative anellovirus abundance and individual anellovirus diversity (alpha diversity, as measured by multiple metrics) was not associated with atopic dermatitis or asthma development (Figs. 6A, B, 7A, B and Supplementary Data 6). However, the anellovirus Shannon index and evenness were lower in children who developed allergic rhinitis by 6 years of age (Shannon index: estimate 0.6 [0.41, 0.98], p = 0.04 and evenness: 0.1 [0.01, 0.67], p = 0.01, Fig. 8A, B and Supplementary Data 6). We delved further into allergic sensitisation to food and inhaled allergens at ages 6 months, 18 months or 6 years. Here we found suggestive results of aeroallergen sensitisation at 18 months associated with evenness (estimate 0.04 [0.00, 1.20], p = 0.067) and Shannon index (0.54 [0.25, 1.09], p = 0.089), and food sensitisation at 6 years associated with inverse Simpson index (1.04 [1.00, 1.08], p = 0.058) (Supplementary Data 6). All other tests on allergic sensitisations were not significant. Likewise, the overall anellovirus composition (beta diversity, as measured by multiple dissimilarity metrics) was not associated with atopic dermatitis or asthma development (PERMANOVA adjusted for technical covariates, Figs. 6C, 7C and Supplementary Fig. 3A, B, 3C, D and Supplementary Data 7). Nevertheless, Bray-Curtis dissimilarity was significantly associated with allergic rhinitis development (F = 1.62, R2(%) = 0.26, p = 0.01, PERMANOVA adjusted for covariates, Fig. 8C and Supplementary Fig. 3E, F, and Supplementary Data 7). When investigating the beta diversity in relation to allergic sensitisation, we found several metrics to be significantly different between cases and controls of any sensitisation at 6 years and aeroallergen sensitisation at 6 years and a suggestive association with aeroallergen sensitisation at 6 months (Supplementary Data 7). All other tests on sensitisation were not significant.

To investigate the relationship between the gut anellovirus composition and atopic disease at a more detailed level, we tested the differential relative abundance of anelloviral genera, subgenera, and species for all three atopic diseases. None were significantly different (Supplementary Figs. 4A, D, G, 5A, D, G for subgenera and genera, respectively, Supplementary Data 8). In addition to this, we employed an indicator species analysis, aiming to find genera, subgenera, and species that were indicative of cases or controls for all three atopic diseases. None of the results were significant after adjusting for multiple testing (Supplementary Figs. 4B, C, E, F, H, I, 5B, C, E, F, H, I for subgenera and genera, respectively, Supplementary Data 9). Moreover, we employed two separate multivariate methods, namely sparse partial least squares regression (sPLS) and random forest, for investigating the relationship between the faecal anellovirus relative abundance at the species or subgenus level at age 1 year and the development of these atopic diseases. These methods did not yield sufficiently discriminatory results, with a median area under the curve (AUC) for the test set below 0.562 (Supplementary Data 10).

Methods
Study population
The study is composed of samples collected in the Danish population-based Copenhagen Prospective Studies on Asthma in Childhood (COPSAC 2010) mother-child cohort. Samples were collected from 700 children and their mothers, starting at 24 weeks of pregnancy up to 6 years of age.

Inclusion and ethics
The study was conducted in accordance with the guiding principles of the Declaration of Helsinki and was approved by the National Committee on Health Research Ethics (H-B-2008-093) and the Danish Data Protection Agency (2015-41-3696). Both parents gave written informed consent before enrolment.

Faecal sample collection, virome extraction and computational processing
Faecal samples from 660 children at 1 year of age were collected either by the parents following comprehensive instructions, or at the clinic. Samples were sent to the laboratory by post and mixed with 1 ml glycerol broth (10% vol/vol glycerol) on arrival and stored at −80°C. Detailed protocols on extraction and sequencing have been described previously. In short, to obtain viromes the faecal samples were filtered to enrich for viral particles and the DNA was extracted using the QIAmp viral RNA mini kit and underwent multiple displacement amplification for 30 min. The virome libraries were prepared following the manufacturer’s directions for the Illumina Nextera XT kit (FC-131-1096) and sequenced on the Illumina HiSeq X platform (2 × 150 bp paired-end reads), resulting in an average depth of 3 Gb per sample. The sequencing was successful for 647 of the 660 samples. The reads were filtered, trimmed using Fastq Quality Trimmer/Filter v0.0.14 (options -Q33 -t 13 -l 32 -p 90 -q 13), and adaptors were removed using cutadapt (v2.0). The processed reads were then de-replicated using VSEARCH (v2.4.3) “derep_prefix” and assembled with Spades (v3.10.1) using the meta flag with hamming disabled. Resulting contigs were deduplicated across all samples at the 95% average nucleotide identity (ANI)-level to yield vOTUs. Virus family-level clusters of vOTUs were first defined based on encoded protein similarity in order to identify all vOTUs belonging to the anellovirus family. This protein-level identification step was to ensure sensitivity towards novel anelloviruses distant from those already found in public databases. The workflow is described in detail on GitHub (https://github.com/shiraz-shah/VFCs). In short, pairwise alignments (from FASTA v. 36) between all vOTU-encoded open reading frames were used to (1) construct an all-against-all vOTU distance matrix for a neighbour-joining tree that was cut deeply, at the virus family-level, and (2) for Markov Clustering in order to define viral gene families, that had their functions resolved by profile-profile alignment against the Pfam database. Finally, the family-level vOTU cluster that carried the anelloviral ORF1 gene family was selected as the anellovirus family in our data. Next, whole-genome Smith−Waterman alignments were made at the DNA-level, between all identified anelloviral vOTUs in addition to public genome representatives for all ICTV-ratified genera, all the way from alpha- to zayintorquevirus. The whole-genome nucleotide alignments were used to define pairwise distances for a neighbour-joining tree. The tree was rooted using iotatorquevirus (known to infect non-human mammals such as pigs) as an outgroup. Subsequently, cut distances from the root were measured that would reproduce all the ICTV genera. This revealed two additional and undescribed genera that were named based on the vacant Greek letters tau and phi. A similar species-level root cut distance was derived based on the reproduction of select ICTV species belonging to alpha, beta and gammatorquevirus. The representative ICTV species were chosen because they had been revised recently, thus likely comprising borderline cases with a rigorous and up-to-date application of the most recent ICTV species-level cutoff. Cutting the tree at a distance reproducing those species enabled comprehensive species-level resolution for the entire dataset. While inspecting the tree distances, an additional local minimum was found between the genus and species-level root distance, and this was used as the cut distance for defining anellovirus subgenera.

Bacterial extraction and computational processing
Details on the processing of faecal samples for 16S rRNA sequencing and processing have been published previously. In short, the bacterial DNA was extracted from the faecal samples using MoBio PowerSoil kits targeting the V4 region of the 16S rRNA gene and sequenced using the v2 kit on the Illumina MiSeq platform to generate paired-end 250 bp reads. Samples with more than 2000 reads were kept before primers were removed, the reads were trimmed and ASVs were inferred.

Clinical phenotypes
Asthma. The asthma diagnosis was assessed prospectively based on a quantitative symptom algorithm published previously. A case in this study defined as having ever had an asthma diagnosis by the age of 6 years.

Atopic dermatitis. The atopic dermatitis diagnosis was based on the criteria of Hanifin and Rajka and with clinical assessment throughout childhood. To be labelled as a case, at least three of four major criteria and at least 3 of 23 minor criteria were needed, and a case was considered as having ever had an atopic dermatitis diagnosis by the age of 6 years.

Allergic rhinitis. Allergic rhinitis at age 6 years was diagnosed based on allergic sensitisation and clinical interviews of the parents on history of significant nasal congestion, sneezing and/or runny nose outside periods with common cold and in the relevant period of exposure to the sensitised aeroallergen. If the participant did not self-report rhinitis, the participant was labelled as a control. If the participant self-reported rhinitis and at least one of the allergens was positive, the participant was labelled as a case. If the participant self-reported rhinitis and all allergens were negative, the participant was labelled as a control. If the participant self-reported rhinitis and the above-mentioned criteria for being a case was not fulfilled, and at least one of the allergens or antibodies were blank, the participant case status was set as missing. If the participant self-reported rhinitis but the IgE and skin prick test were inconclusive, the participant's case status was set as missing.

Infections. Symptoms of infections included colds, fever, and gastroenteritis within (+-) 7 days of the faecal sample were collected by daily diaries.

Covariates. Information on physician-diagnosed asthma in the mother and father, maternal age at birth, smoking during pregnancy, antibiotics in pregnancy, preeclampsia, gestational age, birth weight, C-section or vaginal birth, induction of birth, parents atopy, antibiotics administered to the child, pet ownership, older children in the home, introduction of foods (cow’s milk, eggs, fish and hypoallergenic formula), daycare start and duration of exclusive and total breastfeeding period was obtained during the scheduled visits to the research clinic. C-section was either encoded as 'C-section': any C-section (acute or planned) compared to any vaginal delivery (induced or spontaneous, or as 'Emergency C-section': acute C-section compared to spontaneous vaginal delivery. Information on intrapartum antibiotics was validated against birth records from hospitals, and information on antibiotic use during pregnancy and childhood was validated against national registries. The living environment was based on home address at birth using the CORINE satellite-based land cover database. Fish oil and high-dose vitamin D intervention groups in pregnancy were included.

Statistical analyses
All analyses were conducted in the programming language R version 4.2.2 (Vienna).

Preparation and normalisation of relative abundance data. We calculated the relative abundance per sample, including all viruses detected in our samples, because we wanted to investigate the role of anelloviruses relative to the entire virome community, partly due to the potential biological relationship between these viral compartments, as demonstrated by the procrustes analysis described in the results section 'Anelloviruses are uniquely yet highly prevalent in the 1-year-old infant gut'.

Technical covariates. Several analyses were adjusted for covariates as noted; the selected covariates were the proportion of viral reads, age at faecal sampling, sequencing lane and sequencing depth. These were chosen based on their relationship with the anellovirus composition as measured by several alpha and beta diversity measures.

Alpha and beta diversity calculations. Alpha diversity was calculated at the species level without filtering, with Faith’s phylogenetic diversity (PD) (package picante version 1.8.2), Shannon index, inverse Simpson’s index (phyloseq package version 1.44.0, estimate_richness function) and evenness (Shannon index divided by the number of observed species). Beta diversity was calculated at the vOTU level without filtering using Aitchison (pseudocount of half of the lowest non-zero value in the data), Bray-Curtis, Jaccard, Canberra (function vegdist from vegan package version 2.6-4) and Spearman distance (function separman.dist from bioDist package version 1.70.0), as well as weighted UniFrac and unweighted UniFrac distance (unifrac function from phyloseq package version 1.44.0). The difference in beta diversity between cases and controls was tested using a permutational multivariate analysis of variance using distance matrices (PERMANOVA, adonis2 function from the vegan package version 2.6-4), with 9999 permutations, sequencing lane set as strata, and covariates as described above.

Indicator species analysis. The relative abundance data was filtered to remove species with a prevalence < 50%, all subgenera and genera were included for analyses. Indicator species were identified with the multipatt function (999 permutations) from the indicspecies package version 1.7.14). This allows for the identification of particular taxa that are indicators of certain sample groups, such as cases or controls of clinical outcomes. Indicator species are defined as either being specific (positive predictive value, the probability that the sample belongs to a group due to the species being present) or being sensitive (the probability of finding the species in samples belonging to a group). Species displaying high specificity to a group occur only in samples from that group (though not necessarily present in all samples) and are missing from the samples of the other group. Species displaying high sensitivity to a group occur in all samples belonging to that group and are largely missing from samples of the other group.

Differential abundance analysis. Differential abundance of anellovirus taxa was tested with a generalised linear model (glm, family = 'binomial’) adjusted for the above-mentioned covariates, with the relative abundance zero mean unit variance scaled.

Secretor status analysis. Based on the genome-wide genotyping of the blood cell fraction of plasma using the Illumina OmniExpress imputed with the Hapmap and 1000 genomes reference populations, the rs601338 SNP was encoded to either non-secretors (AA allele, N = 101) or secretors (GA or GG allele, N = 475). The statistical test was performed using a generalised linear model (glm, family = 'binomial’) adjusted for the above-mentioned covariates.

Multivariate analyses using sPLS. Sparse partial least squares (sPLS) from the mixOmics package was used to test the multivariate association between the anelloviruses and clinical outcomes. This methodology maximises the covariance between latent variables of the data, rather than the variance explained on the outcome side. This results in a robust model that is resistant to overfitting and can handle high-dimensional data as both dependent and independent variables, even with small sample sizes. This model identifies linear combinations of the variables from the datasets as a measure to reduce the dimensionality of the data. The microbial data agglomerated at the species level was filtered to remove species with a prevalence < 30%, while at the subgenus level, all taxa were kept, before log transforming with a pseudocount of half of the minimum non-zero value. The data was then split into a training (80%) and testing (20%) set ten times (termed stability runs) using the function createDataPartition from the caret package version 6.0-94 to achieve stable estimates of the model accuracy, thus avoiding spurious results from certain sample combinations. The model was trained on the training set using the train function from the caret package with scale and centre as pre-processing, using fivefold repeated cross-validation with three repeats using the trainControl function. The hyperparameters used were keepX (number of taxa to keep) which started at 5 to the maximum number of taxa, with a step size of 4 for subgenera, and a step size of 8 for species, and ncomp (number of components) which was set to 1:4. The model hyperparameter set with the best median area under the curve (AUC) calculated across all repeats was selected and evaluated on the testing set. The median test AUC of the ten stability runs was used to evaluate the model performance.

Multivariate analyses using random forest. A random forest model was additionally used to test the multivariate association between the anelloviruses and clinical outcomes. This type of methodology can identify nonlinear and multivariate relationships between a high number of independent variables to a single dependent variable. Given a sufficient sample size, and with appropriate training and testing splits, these models are robust and can be operated to avoid overfitting. This model was applied to all taxa agglomerated at the subgenus level, or at the species level after removing species with a prevalence < 30%. This was performed on the relative taxa abundance, and 80% of the data were used for training and 20% for testing, over ten separate stability runs to obtain stable estimates of the model accuracy and avoid spurious results from certain sample combinations. The hyperparameters were tested using the train function (method = 'ranger’, importance = 'impurity’) from the caret package version 6.0-94, where the control was set to fivefold repeated cross-validation with three repeats and random search using the trainControl function. The hyperparameter grid used was: ntrees = (500 and 2000), nodesizes = (20, 50 and 80), splitrule = (variance, extratrees and beta), and mtry started at 5 to the maximum number of taxa, with a step size of 4 for subgenera, and a step size of 8 for species. The model parameters with the best median AUC calculated across all repeats was used to evaluate the model on the test set, using the function ranger (importance = impurity) from the package ranger version 0.15.1. The median test AUC of the ten stability runs was used to evaluate the model performance.

Procrustes analysis. Procrustes is a statistical method that compares two ordinations or matrices from, e.g. environmental and microbial data, by scaling (stretching) and rotating them to minimise the sum of squared differences. The protest function thereafter performs permutation of the Procrustes analysis repeatedly to estimate the significance of the Procrustes statistic and returns a correlation-like statistic from the symmetric Procrustes sum of squares. This was used to test associations between the bacteriophage composition and the anellovirus composition (N = 647), and the association between the anellovirus composition and the 16S rRNA gene amplicon-sequencing-based bacterial composition (N = 611). The data were first re-normalised separately, e.g. the anelloviruses and bacteriophages were re-normalised to 100%. The ordination was created using the distance-based redundancy analysis (dbrda) function, conditioned on lane and sequencing depth, and the association between the two ordinations was assessed with the procrustes function (symmetric = TRUE) and tested with the protest function (9999 permutations), all from the vegan package version 2.6-4.

SpiecEasi analysis. SPIEC-EASI (SParse InversE Covariance Estimation for Ecological Association Inference) is a method for identifying ecological networks among microbes (Kurtz et al. 2015). The SpiecEasi package version 0.1.4 was applied to all taxa agglomerated at the subgenus level, with the following settings: method = ‘mb’, lambda.min.ratio = 1e-1, nlambda = 50.

P value threshold for significance. All p values were considered significant at the nominal p < 0.05 level, or at the multiple testing adjusted (p.adjust function with method fdr) p < 0.05 level as noted. All statistical tests were two-tailed.

No statistical method was used to predetermine sample size. No data were excluded from the analyses.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This paper by Boulund, et al, describes fecal anellovirus populations in 1 year old infants from a very large study of 647 one year old children from the Copenhagen Prospective Studies on Asthma in Childhood study. This is an incredible resource in terms of samples and clinical data for analysis.\n\nThe general process used for analyzing the samples is as follows: Sequence data was generated from filtered stool samples in order to analyze viruses present in the stool. Sequences were assembled and analyzed as viral OTUs. Prevalence and relative abundance of anelloviruses were associated with various clinical correlates using a variety of statistical tests.\nWhile these steps are generally reasonable, there is an absence of critical experimental details that are needed to interpret the results.\n\nLine 282: “…previous studies have suggested potential links to chronic diseases and immune system development.”\nCitation needed.\n\nCitations are given for the extraction and sequencing are given, but more detail is needed to describe what exactly is done in this paper.\n- What amount of stool was used for virome extraction?\n- How was DNA extracted? Kit?\n- Was anything done to determine how effective the viral enrichment was? What proportion of the sequence reads generated were viral? Eukaryotic viruses?\n- What proportion of sequences were not included in assemblies/analysis of vOTUs?\n- Average depth, 3Gbp. Range? Was read depth used to scale results at any point?\n- Should add into discussion biases that would be introduced by using MDA. You would be more likely to get small circular viruses like anelloviruses from MDA amplified samples than other DNA viruses that might be present. It seems surprising to only have anelloviruses and adenoviruses representing eukaryotic viruses in your data set. Is there other support for this observation or should it be explained by technical choices? Are other eukaryotic viruses detected in the individual sequencing reads but not assembled or identified in vOTUs?\n- How were reads de-replicated?\n- How were reads assembled?\n- How were condigs deduplicated?\n- If vOTUs are based on gene clusters, is there an overestimate of different viral species? For example, are vOTU from gene 1 and vOTU from gene 2 counted as two different taxa since we do not know if they are from the same genome or not? How are OTUs from different genes incorporated into the neighbor joining trees?\n- Methods for tree needed\n- How is relative abundance calculated? By contigs or underlying reads? Based on eukaryotic viruses? All viruses? Viruses and bacteria?\n- How are diversity metrics calculated? Contigs #’s? All viruses? What taxonomic level?\n\nFigures are not legible, so it is not possible to fully evaluate the results. Font is too small and blurry.\n\nSome connections between points in results and discussion are a little unclear. For example:\nDiscussion: “By contrast, the anellovirus and the bacteriome compositions were not associated.”\nResults: “There was a slight negative correlation between the gut microbial maturation score and the overall relative anelloviral abundance.”\nBoth things may be true, but a little clarification of exactly what is meant in each statement is needed.\n\nThe use of relative abundance sometimes makes interpretation difficult. For example, in the discussion: “The most consistent and strong effect on the anellovirus composition was whether the child had started day care before the sample collection.” Is there any way this could be explained by changes in abundance or diversity in bacteria or phage composition? It seems likely that it would have a big effect since daycare exposure is likely to increase diversity in the gut, as has been previously shown.\n\nWith some additional work and clarification, this work will be of interest to the field of human virome research.", "reviewer_2": "Early life factors shaping infant gut anellovirus composition and associations with childhood atopic disease\nGeneral Comments\nThe strength of the study is the sample size and that this is part of a well-characterised cohort. I find that the most interesting output is that the Anelloviridae composition was associated with whether the child had started day care before the sample collection, which was also significantly associated with a reduced Anelloviruses alpha diversity and abundance. This is because there is a well characterised link between day care, infections, and the shift of the eukaryotic virome.\n\nI find the inclusion of bacteriophages and Adenoviruses not efficiently justified in this study and I suggest that they are removed. The discussion around the relative abundance of Anelloviruses within the total virome is not supported scientifically since the data have not been appropriately normalised to allow for such analysis. For example, genome size can significantly affect relative abundances since fragments of larger genomes have higher probability of detection by shotgun sequencing. Also, the authors have not analysed the differential effect that sample processing, DNA extraction and MDA has in different virus groups. E.g. Differences in the GC, AT, and nucleotide sequence repeats of Anelloviruses and bacteriophages. I do not think the inclusion of 16S data is relevant either since a link between Anelloviruses and the bacteriome is not easily justified.\n\nThe authors have not adequately described the analytical pipelines and they often reference previous studies which I find problematic. My main scepticism is for the fact that the authors have not included negative controls and whole process negative controls in their strategy nor have they screened the reagents for Anelloviruses. In addition, they have not verified the presence of the Anelloviruses species using a second approach to provide confidence on the efficiency of their methods. There is an available Anelloviruses nester PCR available that could use in combination with sequencing to verify results at least in a fraction of samples. Or, have they been able to isolate any Anelloviruses species from the samples?\n\nThe metanalysis could also be significantly improved by unsupervised clustering of the samples based on Anelloviridae composition and then look for statistical association with outcomes. Dirichlet clustering could work nicely because it takes into account the lost biosphere effect of metagenomic datasets. Finally, different arbitrary thresholds are used throughout the metanalysis that confuses the reader and it is difficult to evaluate their meaning.\n\nThe authors need to be very careful when they discuss about Anelloviruses infection versus presence and clarify every time these terms are used.\n\nOverall, I believe that the manuscript cannot be published as it is and will need major restructuring and additional experiments.\n\nSpecific comments:\nIntroduction\nL51: Expand, what kind of diseases?\nL52: Not true, the gut virome is increasingly being recognised as an important component of the gut microbiome in health and disease\nL54: Infect: has this been demonstrated in vivo?\nL54: Which individuals? In what diseases? In which tissue/organ?\nL55: Colonisation: very different biology from infection mentioned above\nL64-68: Gut microbiota are transmitted vertically. There are many good studies about this\nL70: Define infected\nL75: Define infection\nL77: Richness is not the same as load. Load refers to copy numbers, richness in how many different species are present\nL78-79: Not necessarily, especially for Anelloviruses that do not elicit a robust immune response?\nL81-82: Define causal, and conflict with above discussion regarding infection by Anelloviruses.\nL87: Good references that looked into Anelloviruses and their immunology in the airways: PMID: 37221274, PMID: 36595290\nMethods:\nL493-494: Upload code in github\nL496: Describe inclusion and exclusion criteria.\nL505: Virome extraction: Describe the process of sampling processing and DNA isolation including what measures you have taken to reduce reagent contamination. Are whole process (sampling to sequencing), and DNA isolation negative controls included? What is the effect of the reagent “kitome” in your data? Has this been evaluated and considered during sequencing data analysis and metanalysis? Are Anelloviruses present in the reagents used?\nL510: What is the effect of MDA on Anelloviruses detection? What is the efficiency of the amplification across samples?\nL513: Provide sequencing read data and descriptive statistics. Provide descriptive of sequencing read assignment to taxonomic levels. How many of your sequencing reads are assigned to Order, family, genus, subgenus, species?\nL514: Why these samples failed?\nL514,515: Are the sequencing reads decontaminated from host sequences?\nL517: What is the reason for analysing amino acid sequences and not genomic?\nL518: Provide statistics on counts, contig size, coverage, and gene information. Can you infer function? Is there any normalisation based on contig size and/or coverage performed?\nL524: Define the genome sequence similarity thresholds for assigning contigs to strain, species, genera etc.\nL534: I do not understand why the analysis of the bacteriome, bacteriophages and adenoviruses were included in the study. There are no data to support interaction with Anelloviruses. Data produced in other studies and are being used in this study should be efficiently described in the current manuscript.\nL557: What are these periods exactly? Especially the last 5 years they have changed from what they were before 20 years. I do not understand the criteria for assigning positive, negatives etc. Justify this decision.\nL567: Expand\nL571: How are covariates included in the statistical analysis? Sex, seasonality, and geographical location should be included in the covariate list and accounted for in all statistical testing\nL589: How is the virome identified? Provide details on the computational pipelines. Anelloviruses have very small genomes and circular genomes. How does this affect sequencing detection and relative abundance calculation? Using an outcome as a reason to justify a decision to perform a specific analysis is not scientific. Where there any indications from other studies that Anelloviruses can interact with bacteriophages and Adenoviruses that could inform your decision? The only reasonable way to link Anelloviruses with other viruses in the gut is through Enteroviruses present especially in kids that go to nursery. This could be potentially be linked to Anelloviruses both directly and indirectly.\nL596: Sex, seasonality, and geographical location should be included in the covariate list and accounted for in all statistical testing\nL603: Define terminology for ecological indexes.\nL611: Dispersion should also be included in the analysis. There might be dispersion effects rather than shifts.\nL616: This is wrong. Why are you filtering out species observed in less than half of your population? I do not understand the justification. Biologically how is this justified? Provide a graph demonstrating the relationship between Species frequency of detection (occurrence), number of reads and number of contigs. What fraction of the total Anelloviruses dataset species with <50% occurrence represent? This is around 80%? I find Indicator species analysis not informative and with many limitations. A much better way to do this analysis is to use probabilistic clustering that take into account the power law distribution of metagenomic data to assign metacommunities, then group samples and simply investigate within each group the ratio of health versus disease or any other clinical outcome.\nL641, L663, and 683: Here you use different thresholds. There needs to be consistency when using arbitrary assigned thresholds across the different tests. As it is now it looks as fishing for a significant association. Describe the filtered subsets? How your results change from non-filtered to 100% with increments of 10%? What will mean if you see differences?\nL691: There is no biological justification in generating a virome component of Anelloviruses and bacteriophages. Because of this I do not agree in generating a total relative abundance based on them alone. Especially since you have not provided data of how DNA isolation, MDA, and sequencing differentially affect bacteriophages and Anelloviruses detection. What normalisation technique justifies the co-analysis of bacteriophages and Anelloviruses?\nTable 1:\nAre these correct? Explain the table contents. Controls have asthma?\nAsthma ever by 6 years (%)\nCase 140 (21.6)\nControl 487 (75.3)\nMissing 20 (3.1)\nAtopic dermatitis ever by 6 years (%)\nCase 201 (31.1)\nControl 432 (66.8)\nMissing 14 ( 2.2)\nRhinitis by 6 years (%)\nCase 41 (6.3)\nControl 576 (89.0)\nMissing 30 (4.6)\nResults:\nResults should focus on the description of Anelloviruses prevalence in the cohort, their ecology and genetics.\nFigures: Font sizes are too small to read. All figures need reformatting.\nFigure 1D: thresholds are arbitrary\nFigure 1E: Why are you analysing subgenera and not species?\nFigure 1F: Very weak associations. What does this mean?\nFigure 2: should be Figure 1\nL102: Remove bacteriophages, and any other component of the virome and focus on Anelloviruses. These data have been published before and they do not add value in the current report.\nL116: This section described the presence of Anelloviruses. It should be expanded to include contig information, genome coverage and what genes are encoded. Are all encoded proteins known? What is the role of the genes in the life cycle of the virus> Do we know if they interact with human signalling pathways? Present the species and contig distribution data in your population. Compare species accumulation rare fraction curves in the population and in the different clinical groups. How much of the estimated true diversity have you captured with your approach?\nL174: Expected since the variables are linked!\nL178: Spearman association do not reflect ecological competition or collaboration. There are appropriate packages for this type of analysis, e.g. https://github.com/zdk123/SpiecEasi\nDiscussion:\nL387-L398: You have not demonstrated Anellovirus infection or any other infection in your cohort? So how superinfection exclusion is relevant?", "reviewer_3": "The manuscript titled “Early Life Factors Shaping Infant Gut Anellovirus Composition and Associations with Childhood Atopic Disease” by Boulund et al. focuses on factors that may influence the anelloviral community in the gastrointestinal tract of young children and explores potential links with atopic disease. The study identifies several associations between social, environmental, and genetic factors and the anellovirus community; however, none were found with atopic disease, except for a weak association with rhinitis and allergic sensitization. The key takeaway from a clinical perspective is that the child cohort under investigation is large enough to support the credibility of these negative findings.\n\nMajor Comment:\nFrom the perspective of the readership of Nature Communications, it is questionable whether the presented data is comprehensive enough for publication in this forum. In its current form, the manuscript primarily reports microbial data without correlating it to clinical or immune functional data, which should be available in the cohort.\n\nWhile this is a highly interesting and important field with many unresolved questions about the influence of anellovirus on physiology and immune competence, it would be of utmost interest to learn more about the potential role that anellovirus might play in the maturation of the microbiome, its possible association with host immune responses, and whether anellovirus could serve as a biomarker of immune function.\n\nMinor Comments:\nIn lines 121 and 128, please clarify for non-bioinformatician readers what the functon(\"AND\" >35%) means.\nDo the anellovirus profiles in the children resemble those of their mothers?\nLine 294: The results in the preceding lines were also not significant. Therefore, it is not accurate to state in this sentence that “no other results were significant.”"}, "author_response": "Copenhagen Prospective Studies on Asthma in Childhood\nJakob Stokholm\nstokholm@copsac.com\nOctober 22, 2024\nNature Communications\nPoint-by-point response to reviewer comments\nReviewer #1 (Remarks to the Author):\nThis paper by Boulund, et al, describes fecal anellovirus populations in 1 year old infants from a\nvery large study of 647 one year old children from the Copenhagen Prospective Studies on\nAsthma in Childhood study. This is an incredible resource in terms of samples and clinical data\nfor analysis.\nThe general process used for analyzing the samples is as follows: Sequence data was\ngenerated from filtered stool samples in order to analyze viruses present in the stool.\nSequences were assembled and analyzed as viral OTUs. Prevalence and relative abundance of\nanelloviruses were associated with various clinical correlates using a variety of statistical tests.\nWhile these steps are generally reasonable, there is an absence of critical experimental details that\nare needed to interpret the results.\nAnswer: We appreciate the comments from the reviewer and have updated the manuscript to\ninclude more details regarding the experimental design and process.\n1. Line 282: “…previous studies have suggested potential links to chronic diseases\nand immune system development.” Citation needed.\nRegarding chronic diseases, we referenced in the introduction for example IBD (DOI:\n10.1016/j.cell.2015.01.002) and asthma (DOI: 10.3389/fmicb.2018.00686, 10.1186/s12866-024-\n03187-7 and 10.3390/v14020278), which we have now added to this particular sentence.\nFurthermore, there are hypotheses regarding anellovirus influence on immune system\nmodulation and development in early life (DOI: 10.2174/0929867324666171005112921,\n10.1093/femsre/fuaa007), we have now added these citations to the manuscript as well, please\nsee lines 256-257.\nCitations are given for the extraction and sequencing are given, but more detail is needed to\ndescribe what exactly is done in this paper.\nAnswer: We agree that the original version of this manuscript was too sparse in regards to the\ndetails of the methods, and as such we have added extensively to this section. Please see the\nspecific answers below.\n2. - What amount of stool was used for virome extraction?\nThis has now been detailed in the methods, please see lines 500-501. The original stool\nsamples were collected via small spoons from the collection tube, and 2 spoonfuls were\npreserved in a 10% glycerol and beef broth stock solution (1mL), before freezing in -80C. For\nthe virome extractions, 300 ul of this stock was thawed per sample, corresponding to\napproximately 150 mg faecal material.\n3. - How was DNA extracted? Kit?\nThis was mentioned in the Methods section “Faecal sample collection, virome extraction and\ncomputational processing” on lines 511-514; The virome libraries were prepared following the\nmanufacturer’s directions for the Illumina Nextera XT kit (FC-131-1096) and sequenced on the\nIllumina HiSeq X platform (2 x 150 bp paired-end reads).\n4. - Was anything done to determine how effective the viral enrichment was? What\nproportion of the sequence reads generated were viral? Eukaryotic viruses?\nWe have now added a reference to our previous description of this dataset\n(https://doi.org/10.1038/s41564-023-01345-7), and a summary of this in the methods section of\nthis paper, please see lines 508-511. In total, 51% of the reads mapped to non-viral sequence\nclusters, while 49% of the reads mapped to vOTUs. It is estimated that approximately 1-2% of\nthe DNA in a faecal sample is of viral origin. With this in mind, the enrichment factor was\nestimated to be between 25 and 50 times enrichment of viral DNA over cellular DNA when\ncompared to shotgun metagenomes for the same samples.\n5. - What proportion of sequences were not included in assemblies/analysis of\nvOTUs?\nThis has been clarified in the methods section, please see lines 516-519. While MDA\namplification allows for the detection and quantification of ssDNA viruses, it does introduce\nredundancy in the sequenced reads such that multiple reads with the same sequence occur. In\naddition, virome extractions often have short inserts, leading to sequencing adapters. For these\nreasons, on top of the standard read quality thresholds (minimum PHRED quality and length),\nwe removed reads that were shorter than 32 bp after adapter removal, and we collapsed reads\nthat were identical in sequence. On average this filtering removed 78% (+/- 13%) of the raw\nreads per sample. The filtered subset of reads, comprising on average 22% of the original\nnumber of sequenced reads, was used for both assembly and mapping.\n6. - Average depth, 3Gbp. Range? Was read depth used to scale results at any\npoint?\nThe range of sequencing depth per sample was 0.01 to 24.63Gbp (mean 3Gbps +/- 1.6 Gbps).\nWe have now added this to the methods, please see lines 511-514. The total number of reads\nper sample was included as a technical covariate in statistical analyses when possible as\ndescribed in the particular methods sections, and was also noted in the methods section\n“Technical covariates”, please see lines 617-621.\n7. - Should add into discussion biases that would be introduced by using MDA. You\nwould be more likely to get small circular viruses like anelloviruses from MDA\namplified samples than other DNA viruses that might be present. It seems\nsurprising to only have anelloviruses and adenoviruses representing eukaryotic\n2\nviruses in your data set. Is there other support for this observation or should it be\nexplained by technical choices? Are other eukaryotic viruses detected in the\nindividual sequencing reads but not assembled or identified in vOTUs?\nWe have added some discussion of this and the potential limitations in the manuscript, please\nsee lines 426-431.\nAs the reviewer correctly notes, MDA amplification does introduce biases, with circular ssDNA\nbeing preferentially amplified. We curated the assemblies manually and did not detect other\neukaryotic viruses in our data besides dsDNA adenoviruses and ssDNA anelloviruses. Most\nhuman gut viruses are RNA viruses that would not have been amplified by our protocol.\nThere are some eukaryotic DNA viruses that have been identified in infant faeces previously,\nthat we did not identify in our data, e.g. Circoviridae, Geminiviridae, Nanoviridae, Parvoviridae,\nPolyomaviridae, and Genomoviridae (DOI: 10.1038/nm.3950, 10.1073/pnas.1514285112,\n10.1128/JCM.01589-12, 10.1371/journal.pone.0240958, 10.1016/j.chom.2023.01.003,\n10.1186/s40168-019-0766-7, 10.1016/j.chom.2019.09.009). Nanoviridae have shorter genomes\nand would be removed based on our contig length cutoff of 1000 bp. There may be biological\nexplanations for the variations in viral identification compared to our data, e.g. variations in\nbreastfeeding rate and duration, geographical area, and birth mode. Furthermore, in the majority\nof the mentioned studies anelloviruses are by and large the most prevalent and abundant virus\nidentified. In several other studies (DOI: 10.1038/s41598-022-26707-9,\n10.3389/fmicb.2018.01162) anelloviruses are sometimes the only DNA virus identified in infant\nfaeces. Similarly, it is generally believed that there is a very low abundance of eukaryotic\nviruses in the gut of (adult) healthy individuals, besides anelloviridae (DOI: 10.1146/annurev-\nmicro-102215-095431). To conclude, the identification of adenoviruses and anelloviruses as the\nonly eukaryotic viruses in this dataset does not appear to be out of the ordinary in a biological\naspect.\nHowever, as the reviewer astutely has noted, there could also be technical variations that can\nexplain this. Some of the previously mentioned studies utilised PCR rather than sequencing, or\nmapped reads to viral databases. While those approaches are more sensitive to known viruses,\nthey ignore any viral dark matter, which still represents the majority of the viral diversity in\nvirome samples to date. Our method was based on assembly and conservative curation. Thus\nwe might have missed some minor groups of viruses, but any major groups were certainly\npicked up, even if they showed no similarity to any known viruses in databases. We published\nthis manually curated infant gut virome data set in the baseline paper (DOI: 10.1038/s41564-\n023-01345-7), and made all the identified viral sequences freely available as a database which\nis now “frozen” to ensure reproducibility across any studies that wish to utilise this data set. This\nstudy is based on that data set.\nAfter the publication of our baseline paper, an automated virus-prediction bioinformatics tool,\ngeNomad, was published (DOI: 10.1038/s41587-023-01982-7) that performs on par with our\nmanual curation (previous tools lacked sensitivity). Motivated by the reviewer’s question, we ran\nthis pipeline on our data and found that it does pick up additional eukaryotic viruses, but only\n3\ndetects 697 out of the 970 anelloviruses that we found manually. The additional eukaryotic\nviruses that geNomad finds were 83 CRESS viruses (including 7 circoviruses) and 42\nGenomoviridae. So some additional eukaryotic viral groups are present, albeit minor, and we\nwant to acknowledge that as a weakness of our approach. Virome annotation is still a new field\nwith ongoing development and new tools becoming available regularly. There is still no ground\ntruth or generally accepted method/pipeline that allows for an ideal balance between false\npositives and false negatives, and there is currently ongoing research into how well these newly\ndeveloped pipelines perform on different types of data (DOI: 10.1186/s13059-024-03236-4).\nThus, for consistency, we aim to keep our high-quality dataset “locked” for the purpose of the\npresent study, where we would like to maintain that our curated database was the best available\noption.\n8. - How were reads de-replicated?\nThe trimmed reads were de-replicated with VSEARCH version 2.4.3. We have added this to the\nmethods, please see lines 515-518.\n9. - How were reads assembled?\nThe reads were assembled using Spades version 3.10.1. We have added this to the methods,\nplease see lines 520-523.\n10. - How were condigs deduplicated?\nThe contigs (1.5 million) were de-duplicated if they were above 1 kb in size, to 362,000 95% ANI\n(average nucleotide identity) representatives. We have added this to the methods, please see\nlines 520-523.\n11. - If vOTUs are based on gene clusters, is there an overestimate of different viral\nspecies? For example, are vOTU from gene 1 and vOTU from gene 2 counted as\ntwo different taxa since we do not know if they are from the same genome or not?\nHow are OTUs from different genes incorporated into the neighbor joining trees?\nAfter assembling contigs from each sample, contigs across samples were merged into vOTUs if\nthey showed more than 95% average nucleotide identity (ANI), i.e. at the DNA level. From these\nvOTUs, which included bacteriophages etc., the group of vOTUs that corresponded to\nAnelloviruses was identified by their clustering at the protein-level. Although anelloviruses are\nextremely diverse in sequence, they all appear to encode the same proteins, albeit in highly\ndivergent forms. This protein-level similarity is what we used to isolate anelloviruses as a group,\nnot to infer taxonomy. The neighbour joining tree that we constructed to sort out anelloviral\ntaxonomy was based on whole-genome DNA alignments, so does not contain gene-information\nexplicitly. We hope this answers the reviewers question and are happy to elaborate here or in\nthe manuscript if needed.\n12. - Methods for tree needed\nThe methods for the tree were described in the Methods section “Faecal sample collection,\nvirome extraction and computational processing”, please see lines 526-547.\n4\n13. - How is relative abundance calculated? By contigs or underlying reads?\nBased on eukaryotic viruses? All viruses? Viruses and bacteria?\nWe have expanded this in the methods, please see lines 527-530. Reads from each sample\nwere mapped against the global vOTU set (described in question 11 above), and individual\nvOTU abundances were estimated by normalising the total read count per vOTU against the\nlength of that vOTU, as well as the total number of quality controlled reads in the sample. Thus,\nthe abundances are estimated based on read mapping to the assembled vOTUs. This was done\nfor all viruses, including bacteriophages, but not bacteria. The total relative abundance per\nsample is thus 100% for the entire virome. Hereafter, depending on the analyses, the dataset\nwas subset to e.g. adenoviruses or anelloviruses, without re-normalisation.\n14. - How are diversity metrics calculated? Contigs #’s? All viruses? What\ntaxonomic level?\nAll diversity metrics were calculated solely based on anelloviruses (or e.g. adenoviruses or\nbacteriophages in Figure 1C), however the relative abundance included all viruses, as\ndescribed above. The alpha diversity metrics were calculated on the species level, whereas the\nbeta diversity metrics were calculated on the vOTU level, please see methods section “Alpha\nand beta diversity calculations”, lines 623-644.\n15. Figures are not legible, so it is not possible to fully evaluate the results. Font is\ntoo small and blurry.\nThank you for pointing this out, the figures lost quality when submitted directly in the word\ndocument. We have now included all the figures as high-quality PDFs which are present at the\nend of the submitted document. Each panel can be zoomed in individually with high resolution.\n16. Some connections between points in results and discussion are a little unclear.\nFor example: Discussion: “By contrast, the anellovirus and the bacteriome\ncompositions were not associated.” Results: “There was a slight negative\ncorrelation between the gut microbial maturation score and the overall relative\nanelloviral abundance.” Both things may be true, but a little clarification of exactly\nwhat is meant in each statement is needed.\nWe understand that this is not entirely clear and have now added further explanation of this to\nthe manuscript, please see lines 248-250 and 307-308. The gut microbial maturation score is a\nscore calculated and presented in a previous study from the same 1-year faecal samples (DOI:\n10.1038/s41467-017-02573-2). This is a normalised score that is related to the overall\nmaturation of the gut bacterial microbiome in early life. We hypothesised that the microbial\nmaturation of the gut microbiome may be related to immune development and thus the overall\nabundance of anellovirus, which we also believe is related to immune development. However,\nthis does not necessarily mean that the actual microbial composition is related to the anellovirus\ncomposition.\n17. The use of relative abundance sometimes makes interpretation difficult. For\nexample, in the discussion: “The most consistent and strong effect on the\nanellovirus composition was whether the child had started day care before the\n5\nsample collection.” Is there any way this could be explained by changes in\nabundance or diversity in bacteria or phage composition? It seems likely that it\nwould have a big effect since daycare exposure is likely to increase diversity in\nthe gut, as has been previously shown.\nThis is an excellent point raised. As we processed and sequenced the anelloviruses,\nadenoviruses and bacteriophages together, it is impossible with the current data to fully\ndisentangle compositionally dependent effects. As you mentioned, the diversity and abundance\nof bacteria and bacteriophages is known to increase in early life and some of this may be\ninfluenced by the start of daycare (as suggested by some: e.g. DOI: 10.1038/s41522-021-\n00265-w). However, we did not report significant findings between daycare and the microbial\nbeta diversity at the 1 year fecal samples in a previous study from the same cohort (Fig 3. DOI:\nhttps://doi.org/10.1186/s41468-022-01305-z), similarly, we did not find a significant association\nin the microbial alpha diversity (Faith’s PD and richness) and daycare (p > 0.05). This suggests\nthat there is no effect of daycare on the alpha or beta diversity of the bacteriome of our samples,\nand by extension there should be little effect on the bacteriophageome as this reflects the\nbacteriome.\nThus, our results regarding daycare, gestational age, older siblings, and microbiota maturity,\nsuggest that these are important factors in the overall anellovirus abundance and diversity, and\nare likely not a remnant of a compositional effect of these factors on the bacteriophageome.\nThis is in line with findings and hypotheses from previous studies that demonstrate a\nrelationship between the immune system and anellovirome. However, as we state in the\ndiscussion (lines 342-344), this must be validated in longitudinal studies that investigate the\nanellovirome immediately before and after daycare initiation. These studies could make use of\ne.g. PCR or other quantitative methods to avoid the complications of compositionality.\nWith some additional work and clarification, this work will be of interest to the field of human\nvirome research.\nAnswer: We are very appreciative of the reviewer for the relevant questions and suggestions\non this work, we believe this has greatly improved the manuscript.\nReviewer #2 (Remarks to the Author):\nEarly life factors shaping infant gut anellovirus composition and associations with childhood\natopic disease\nGeneral Comments\nThe strength of the study is the sample size and that this is part of a well-characterised cohort. I\nfind that the most interesting output is that the Anelloviridae composition was associated with\nwhether the child had started day care before the sample collection, which was also significantly\nassociated with a reduced Anelloviruses alpha diversity and abundance. This is because there\nis a well characterised link between day care, infections, and the shift of the eukaryotic virome.\nAnswer: We greatly appreciate the reviewers notice of some of the strengths of this study, and\na shared interest in some of the strongest findings.\nI find the inclusion of bacteriophages and Adenoviruses not efficiently justified in this study and I\nsuggest that they are removed. The discussion around the relative abundance of Anelloviruses\nwithin the total virome is not supported scientifically since the data have not been appropriately\nnormalised to allow for such analysis. For example, genome size can significantly affect relative\nabundances since fragments of larger genomes have higher probability of detection by\nshotgun sequencing. Also, the authors have not analysed the differential effect that sample processing,\nDNA extraction and MDA has in different virus groups. E.g. Differences in the GC,\nAT, and nucleotide sequence repeats of Anelloviruses and bacteriophages. I do not think the\ninclusion of 16S data is relevant either since a link between Anelloviruses and the\nbacteriome is not easily justified.\nAnswer: We regret that the rationale for our analysis of the entire virome at the outset of this\nstudy was not clear. We have clarified this throughout the specific questions raised by the\nreviewer, as well as in the manuscript itself.\nThe authors have not adequately described the analytical pipelines and they often reference\nprevious studies which I find problematic. My main scepticism is for the fact that the authors\nhave not included negative controls and whole process negative controls in their strategy nor\nhave they screened the reagents for Anelloviruses. In addition, they have not verified the\npresence of the Anelloviruses species using a second approach to provide confidence on the\nefficiency of their methods. There is an available Anelloviruses nester PCR available that could\nuse in combination with sequencing to verify results at least in a fraction of samples. Or,\nhave they been able to isolate any Anelloviruses species from the samples?\nAnswer: We have now added a considerable amount of detail to describe the analytical\npipelines used in this manuscript. Additionally, we included information on the whole process\nnegative controls that were included in the analysis, including the fact that the negative controls\ndid not contain any anelloviruses. Unfortunately we have limited faecal material left over from\nthese 14 year old samples, and this prohibits us from performing PCR on the faecal samples.\nHowever, we have conducted an in silico PCR based on the method suggested (DOI:\n10.1128/JCM.01703-07). We identified 970 vOTUs to be anelloviruses based on sequencing\nand assembly, out of these 970, 869 (90%) contained one or more exact (100%) matches to the\nprimers described in the article, while 275 (28%) vOTUs contained exact 100% matches to the\nsense/antisense pairs of the universal primers. There were no children that did not have any of\nthe 275 in silico PCR-identified anelloviruses (mean 40, range 3-162 per child). The average\nprevalence among samples of these 275 in silico PCR-identified anelloviruses was 84 (median\n78, range 6-461), out of our 647 samples. It is not expected that all anelloviruses identified with\nour sequencing and computational processing will match 100% with the primers, especially\n7\nconsidering the very high genetic diversity among anelloviruses. However, based on the\nphylogenetic tree that we created with several known anellovirus species, all anelloviruses\nidentified in this study are placed to be more similar to each other compared to the outgroup:\niotatorquevirus, an anellovirus species that infects non-human mammals. With this, we believe\nthe DNA extraction, sequencing and subsequent computational processing that lead to the\nidentification of anelloviruses in our samples is reliable and robust.\nThe metanalysis could also be significantly improved by unsupervised clustering of the samples\nbased on Anelloviridae composition and then look for statistical association with outcomes.\nDirichlet clustering could work nicely because it takes into account the lost biosphere effect of\nmetagenomic datasets. Finally, different arbitrary thresholds are used throughout the\nmetanalysis that confuses the reader and it is difficult to evaluate their meaning.\nThe authors need to be very careful when they discuss about Anelloviruses infection versus\npresence and clarify every time these terms are used.\nOverall, I believe that the manuscript cannot be published as it is and will need major\nrestructuring and additional experiments.\nSpecific comments:\nIntroduction\n1. L51: Expand, what kind of diseases?\nWe have now updated this sentence with references. We recently published a report\nassociating the infant gut virome with the development of the chronic disease childhood asthma\n(DOI: 10.1038/s41591-023-02685-x). There is also an established relationship between\nadenovirus and acute gastroenteritis (DOI: 10.1136/bmj.296.6617.229-a).\n2. L52: Not true, the gut virome is increasingly being recognised as an important\ncomponent of the gut microbiome in health and disease\nWe thank the reviewer for pointing this out. We have now updated the sentence.\n3. L54: Infect: has this been demonstrated in vivo?\nIt is generally accepted and referred to that anelloviruses infect human cells, and are detectable\nby e.g. PCR (DOI: 10.3390/v14112482, 10.3390/v14061277, 10.1007/978-3-540-70972-5_1).\nWe quote from Okamoto et al. 2009 (DOI: 10.1007/978-3-540-70972-5_1) “Infection with TTV is\ncharacterised by persistent lifelong viremia in humans, with circulating levels of up to 106\ncopies/ml in the general population (Hu et al. 2005; Pistello et al. 2001). As an indication of its\ngreat replicative capacity in vivo, a study on the kinetics of clearance of TTV viremia suggested\nthat a daily production rate of >1010 virions is required to maintain the observed levels of viremia\n(Maggi et al. 2001b).” Likewise, we quote from Spezia et al. 2023 (DOI:\n10.1016/j.amolm.2023.100006): “TTV is known to infect various cell types and can be found\nthroughout the body.”.\nprevious anellovirus literature. The term “infection” has been used in studies reporting either\nPCR or sequencing data (DOI: 10.1016/j.amolm.2023.100006,\n10.3390/microorganisms9071357, 10.1007/s00705-015-2363-9). Anelloviral presence has also\nbeen referred to as colonisation in many previous studies (DOI: 10.1111/ajt.13031,\n10.1038/s41598-023-34730-7, 10.3390/v14050865, 10.1016/j.coviro.2022.101248).\nAs the precise mechanism of anellovirus infection in human cells is still unclear, and as the\nscope of this particular manuscript is not to resolve those mechanisms, we believe the use of\nthe word infection, in line with previous literature, is acceptable. We have added a clarification of\nthis in the introduction, please see lines 55-56. However, we are happy to accommodate if the\nreviewer or the editor have specific requests regarding the precise wording.\n4. L54: Which individuals? In what diseases? In which tissue/organ?\nThis is presented more in-depth throughout the introduction, please see lines 59-61, 64-66, 88-\n94.\n5. L55: Colonisation: very different biology from infection mentioned above\nPlease see our reply to your question number 3.\n6. L64-68: Gut microbiota are transmitted vertically. There are many good studies\nabout this\nIndeed, there are an increasing number of studies demonstrating both vertical and horizontal\ntransmission of gut microorganisms. In this section we present studies suggesting several\ndifferent modes of transmission for anelloviruses, e.g. vertically before birth, vertically during or\nshortly after birth, as well as horizontally in infancy.\n7. L70: Define infected\nPlease see our answer to your question number 3.\n8. L75: Define infection\nPlease see our answer to your question number 3..\n9. L77: Richness is not the same as load. Load refers to copy numbers, richness in\nhow many different species are present\nWe completely agree with the reviewers point, and we have changed the wording of this\nsentence.\n10. L78-79: Not necessarily, especially for Anelloviruses that do not elicit a robust\nimmune response?\nIndeed, it is unclear whether changes in anellovirus abundance are due to exposure to sources\nof the virus, or immune system mechanisms (from the mother or the child itself), or a\ncombination of both, as such we have updated the manuscript to reflect this, please see lines\n79-82.\n9\n11. L81-82: Define causal, and conflict with above discussion regarding infection by\nAnelloviruses.\nWe have updated the wording of this sentence, please see our reply to your questions number 3\nand 5.\n12. L87: Good references that looked into Anelloviruses and their immunology in the\nairways: PMID: 37221274, PMID: 36595290\nThank you for these suggestions, we have now updated the references used in this section.\nMethods:\n13. L493-494: Upload code in github\nThis is an excellent suggestion, as such we have now uploaded the code in a GitLab repository,\nplease see https://gitlab.com/copsac_public/infant_gut_anellovirus. We have mentioned this in\nthe methods section, please see lines 479-80.\n14. L496: Describe inclusion and exclusion criteria.\nWe agree that this is essential information that was lacking in the first version of this manuscript.\nWe have now described the inclusion and exclusion criteria for the COPSAC cohort, please\n2010\nsee lines 485-490.\n15. L505: Virome extraction: Describe the process of sampling processing and DNA\nisolation including what measures you have taken to reduce reagent\ncontamination. Are whole process (sampling to sequencing), and DNA isolation\nnegative controls included? What is the effect of the reagent “kitome” in your\ndata? Has this been evaluated and considered during sequencing data analysis\nand metanalysis? Are Anelloviruses present in the reagents used?\nThank you for pointing out this highly relevant information that was missing from the first version\nof this manuscript. We have now added information about the negative controls lab processing\nand computational processing to the methods, please see lines 501-505, and 526-527.\nThe detailed protocol for virome extraction has been published previously (DOI:\n10.3390/v11070667, also referenced in the methods). Beef broth with 10% glycerol (SSI\nDiagnostica, Hillerød, Denmark) and SM buffer (Sodium chloride/Magnesium sulphate buffer)\nwere used as negative controls. DNA from the negative controls (with replicates) were extracted\nand sequenced with the same procedures as the samples, to represent the “kitome”. Reads\ncorresponding to these negative controls are available online alongside the reads for the infant\ngut samples themselves, via the ENA accession number presented in the methods. The\nnegative controls underwent the same computational processing as for the real samples. We\ndid not observe any pattern of vOTUs conserved between the controls, or between replicates of\nthe same reagent, thus we did not correct sample vOTU abundances for the composition of the\ncontrols. Rather, the vOTUs detected in the negative controls appeared to be random faecal\nbacteriophage vOTUs also present in some of the positive samples. We suspect that such\nvOTUs could stem from aerosolic cross-contamination between neighbouring samples on a\n10\nsequencing plate, and they became amplified disproportionately owing to the little native DNA\nwithin the negative controls. For information, the controls had much lower mean richness\nestimates than the true samples, as expected (mean 25 +/- 7 vOTUs in controls, vs 147 +/- 70\nin samples). Of note, none of the negative controls contained any anellovirus vOTUs.\n16. L510: What is the effect of MDA on Anelloviruses detection? What is the efficiency\nof the amplification across samples?\nThis is an excellent point of consideration, as such we have added a paragraph regarding MDA\nin the discussion, please see lines 426-431, and we have added information about the\nenrichment in the methods, please see lines 508-511 and 523-525.\nWhat we know is that MDA introduces biases, preferentially amplifying circular ssDNA,\npotentially leading to a higher detection of anelloviruses. We also identified dsDNA\nadenoviruses in this data. Most human gut viruses are RNA viruses and would not have been\namplified by our protocol. It is estimated that approximately 1-2% of the DNA in a faecal sample\nis of viral (including bacteriophage) origin and that 5 - 10% of the sequenced viral abundance\n(after normalising for genome length) was anelloviral. The overall enrichment factor was\nestimated to be between 25 and 50 times for viral DNA over cellular DNA when compared to our\nshotgun metagenomes for the same samples. This enrichment, however, is most likely a\nreflection of the viral particle enrichment step rather than amplification. As for the amplification\nstep in particular, a recent study has found that MDA overestimates the abundance of small\nssDNA viruses by an order of magnitude compared to larger dsDNA viruses\n(https://doi.org/10.1101/2022.12.12.520144). Of note, we have previously demonstrated the\nnumber of reads that were viral in our previous manuscript on this dataset (Extended Data\nFigure 3 in Shah et al. 2022, DOI: https://doi.org/10.1038/s41564-023-01345-7). In total, 51% of\nthe reads mapped to non-viral sequence clusters, and 49% of reads mapped to vOTUs. The\n17. L513: Provide sequencing read data and descriptive statistics. Provide descriptive\nof sequencing read assignment to taxonomic levels. How many of your\nsequencing reads are assigned to Order, family, genus, subgenus, species?\nThe sequencing read data for viral sequencing and 16S rRNA are available, which can be found\nin the methods section “Data availability”, please see lines 462-466. As we described in the\nmethods section “Faecal sample collection, virome extraction and computational processing”,\nthe average sequencing depth was 3 Gbps +/- 1.6 Gbps (but ranging from 0.01 to 24.63 Gbps)\nper sample. Additionally, the viral reads were mapped to vOTUs as generated by assembly, and\nthereafter the anellovirus species, subgenus, genus and family were defined, thus the reads are\nonly mapped to vOTUs, and the proportion of reads mapping to any of the taxonomic levels are\ndeduced by aggregating the vOTU-level matches by taxonomic affiliation. Across all samples,\nthe anelloviral vOTUs recruited an average of 11000 reads, ranging from 22 to 957.000 reads,\nout of an average 4.5 million QC-passed reads per sample, ranging from 50,000 to 22 million.\n18. L514: Why these samples failed?\nWe have now updated the phrasing of this sentence, please see line 514-515. The 13 failing\nsamples were not able to be sequenced as they had insufficient DNA concentrations after\n11\nlibrary preparation. This may be due to technical issues, or the low amount of starting material\npresent in infant virome samples.\n19. L514,515: Are the sequencing reads decontaminated from host sequences?\nWe decontaminated the vOTUs by manual curation, thus after assembly. Due to the\nmethodological process employed (assembly, vOTU curation, and strict read-mapping criteria\n(95% identity) to the vOTUs), we do not expect any host reads to have mapped to our curated\nvOTUs.\n20. L517: What is the reason for analysing amino acid sequences and not genomic?\nAmino-acid sequences encoded in genes on the vOTUs were used for initial clustering of\nvOTUs into viral families because protein sequences are more conserved and allow for deeper\nclustering than DNA identity. Once all the vOTUs falling within the Anelloviridae viral family were\nidentified, all remaining resolution of anelloviral taxonomy was carried out at the DNA level.\n21. L518: Provide statistics on counts, contig size, coverage, and gene information.\nCan you infer function? Is there any normalisation based on contig size and/or\ncoverage performed?\nOn average, anelloviruses made up 9% of the total relative abundance, but this number ranged\nfrom 0.01% to 93% within individual samples. The mean anelloviral genome size was 3000 bp,\nranging from 1800 to 5600 bp. We employed a contig length cutoff of 1000 bp. As for coverage,\nacross all samples, the anelloviral vOTUs recruited an average of 11000 reads, ranging from 22\nto 957.000 reads. Relative abundances for each vOTU was calculated by taking the number of\nmapping reads and dividing by the genome size, and then again, dividing by the post-QC\nsequencing depth of the sample. Ambiguously mapping reads were distributed evenly between\nmatching vOTUs in the ratio suggested by the non-ambiguously mapping reads. This has been\ndetailed in the methods, please see lines 527-531.\nThe anellovirus genome contains several overlapping open reading frames that encode six\nproteins, and the function of these are not entirely clear (DOI: 10.1093/femsre/fuaa007).\nHowever, to our knowledge, these open reading frames are conserved across most if not all\nanelloviruses. Owing to their small genomes, anelloviruses do not have room for housing any\naccessory genes, so from a gene-functionality perspective, there should be no variation\nbetween anelloviruses in spite of their high genetic variability. Unfortunately, we have not been\nable to systematically survey the functional potential of the 970 anelloviral genomes, as\ndelineating the overlapping ORFs is not straightforward. However, for the three ORFs that were\nreliably identified on the vOTUs, we saw no functional variability.\n22. L524: Define the genome sequence similarity thresholds for assigning contigs to\nstrain, species, genera etc.\nWe have added additional details regarding this in the methods. First, reads were assembled\ninto contigs for each sample separately. Then contigs from all samples were pooled and subject\nto an all-against-all DNA sequence alignment. Sequences which exhibited an ANI over 95%\nwere clustered into vOTUs, corresponding to the species-level cutoff for bacteriophages. For\n12\nanellovirus, this would be strain-level. Protein-coding ORFs were identified on all vOTUs and\ntheir encoded protein sequences were subject to another all-against-all protein alignment to\ninfer deep vOTU clusters corresponding to viral families. The 970 vOTUs belonging to the\nanelloviridae family were extracted, pooled with 166 taxonomically resolved anelloviral\nsequences from GenBank, and subject to another round of all-against-all DNA alignment, this\ntime with a sensitive Smith-Waterman algorithm using Fasta36 package. The alignment scores\nbetween the 970 vOTUs plus 166 reference sequences were used to construct a neighbour-\njoining tree. This tree was cut at depths that reproduced the ICTV-ratified anelloviral genera\n(45% ANI) and species (63% ANI), and we made an additional cut to create subgenera (57%\nANI). These cutoffs closely match the ICTVs own criteria at 54% AND for genera and 65% AND\nfor species, where AND is the average nucleotide divergence, being the opposite of ANI (AND =\n100% - ANI). This has now been clarified in the methods.\n23. L534: I do not understand why the analysis of the bacteriome, bacteriophages and\nadenoviruses were included in the study. There are no data to support interaction\nwith Anelloviruses. Data produced in other studies and are being used in this\nstudy should be efficiently described in the current manuscript.\nWe regret that this was not clear, despite our attempts to explain the rationale throughout the\nmanuscript. We are happy to receive suggestions on how to improve the wording in the\nmanuscript to make this more clear.\nIn short, we have previously characterised the dark matter of the virome in the COPSAC\n2010\ncohort in our first study on this data, whereafter we identified a relationship between\nbacteriophages and asthma (lines 110-112). However, the eukaryotic compartment of this data\nhas never been studied before, providing the basis for this study. Here, we identified\nadenoviruses and anelloviruses, and validated previous findings on the relationship between the\nwell-known adenoviruses and gastroenteritis (lines 149-151), which emphasises the robustness\nof our data. Yet, anelloviruses remain elusive in the previous literature, though several studies\nsuggest links to chronic immune disease or immune development (lines 84-88). Thus,\nconsidering the extensive information we have on immune development and the well-validated\nchronic atopic disease diagnoses in our data, we focused the remainder of our manuscript on\nthis.\nWe do not hypothesise that there is a biological relationship between bacteriophages,\nadenoviruses and anelloviruses. However, the data was created in unison, thus there are\nseveral technical reasons for the links between the viruses in this particular dataset. As we\ndemonstrate in Supplementary Figure 2 and discuss in the manuscript (lines 308-311), we\ncannot reject a (technical) relationship between the bacteriophages and anelloviruses.\nAdditionally, we considered that the total abundance of anelloviruses, relative to the remainder\nof the virome, might be of biological interest. An anellovirus-only feature space would have\nremoved this information. However, we have conducted additional analyses in an anellovirus-\nonly feature space, please see our answer to your question number 34 below.\nand considerations in the analysis choices made. We appreciate any suggestions as to how we\ncan rearrange things to convey this better.\n24. L557: What are these periods exactly? Especially the last 5 years they have\nchanged from what they were before 20 years. I do not understand the criteria for\nassigning positive, negatives etc. Justify this decision.\nWe acknowledge that our description of allergic rhinitis criteria can be difficult to comprehend\nand have expanded and clarified now, please see lines 569-581. Periods refers to either\nepisodes of colds or peak season of allergen exposure. In short, diagnosis requires\naeroallergen sensitization and symptoms of rhinitis, which is not just representing an episode of\na common cold in the child, and these symptoms need to be present during the peak of allergen\nexposure (of the allergen which the child is sensitised towards).\n25. L567: Expand\nWe have updated the methods with additional details, please see lines 583-590. The parents of\nthe participating children filled in daily diaries from the birth of the child until 3 years of age, on\nquestions regarding the presence of symptoms of common cold (non-chronic nasal drip and\nsneezing), fever (>38C), or gastroenteritis (diarrhoea and/or vomiting). If these symptoms were\npresent and registered in the daily diaries within (before or after) 7 days of the time of the faecal\nsample collection, the participant was labelled a case (e.g. having a cold, or having\ngastroenteritis), otherwise the participant was labelled a control. This was described more in a\nprevious manuscript (DOI: 10.1016/j.jaci.2023.11.019).\n26. L571: How are covariates included in the statistical analysis? Sex, seasonality,\nand geographical location should be included in the covariate list and accounted\nfor in all statistical testing\nAs described in the section “Technical covariates” lines 617-621, we adjusted for several\ncovariates, namely the proportion of viral reads per sample, age at faecal sampling, sequencing\nlane, and sequencing depth. The method of covariate adjustment is described separately in\neach section for each statistical method used. Other biological covariates, such as sex and\nsample season were not included as covariates, as these were instead considered relevant\noutcomes and were tested for in the statistical analyses. Regarding geographical location we\nincluded “living environment” as an outcome variable in our studies, which, as described,\nencodes for urban or rural living environments. Please see our response to your question\nnumber 28 below.\n27. L589: How is the virome identified? Provide details on the computational\npipelines. Anelloviruses have very small genomes and circular genomes. How\ndoes this affect sequencing detection and relative abundance calculation? Using\nan outcome as a reason to justify a decision to perform a specific analysis is not\nscientific. Where there any indications from other studies that Anelloviruses can\ninteract with bacteriophages and Adenoviruses that could inform your decision?\nthrough Enteroviruses present especially in kids that go to nursery. This could be\npotentially be linked to Anelloviruses both directly and indirectly.\nWe have described the process for viral sequencing data processing, including vOTU\nidentification and specifically anellovirus species identification and tree-building, in our Methods\nsection “Faecal sample collection, virome extraction and computational processing”, lines 497-\n547. Please see our reply to your questions number 23 and 33.\n28. L596: Sex, seasonality, and geographical location should be included in the\ncovariate list and accounted for in all statistical testing\nThere is, to our knowledge, no research indicating a sex-specific or season-specific relationship\nfor anellovirus prevalence, diversity, richness or abundance. Furthermore, the only evidence on\ngeographical location and anellovirus prevalence or composition relate to large geographical\ndistances across countries. Our cohort is based in a small geographical area (Zealand, an\nisland in the eastern part of Denmark). Additionally, we were interested to investigate whether\nsex, sample season, or living environment (urban or rural) were related to the anellovirus\ncomposition. Therefore, we included these variables as outcomes rather than covariates. As we\ndemonstrated in our results, sex was not related to anellovirus composition. Similarly, urban\nliving environment demonstrated very little relationship to anellovirus subgenus abundance.\nFinally, sample season, which in our study is conflated with birth season (as samples were\ncollected around 1 year of age), was associated with anellovirome composition. Due to our\napproach of analysing these variables as outcomes rather than adjusting for them as covariates,\nthese findings aid future studies such that sex and urban or rural living environment are not\nrelevant covariates to include, whereas sample or birth season may be relevant covariates for\nfuture researchers to consider. We believe this approach avoids unfounded hypotheses\nbiassing the results, and provides greater biological insights which can be of use to future\nresearch.\n29. L603: Define terminology for ecological indexes.\nWe have now included a description of the ecological indices presented in the manuscript,\nplease see lines 623-44.\nFaith’s Phylogenetic Diversity (PD): Measures the total evolutionary branch lengths represented\nby taxa in a community, reflecting their phylogenetic relationships.\nShannon’s Diversity Index: Combines taxa richness and evenness to measure the uncertainty in\npredicting the taxa of a randomly selected individual, where higher values indicate greater\ndiversity.\nAitchison Distance: Measures distance based on log-ratio transformations.\nBray-Curtis Dissimilarity: Quantifies the dissimilarity based on taxa abundances, emphasising\nthe proportion of shared versus unique taxa.\nJaccard Distance: Measures the distance by comparing the number of shared species to the\ntotal number of species, ignoring abundances.\n15\n30. L611: Dispersion should also be included in the analysis. There might be\ndispersion effects rather than shifts.\nThis is a great suggestion, we have now included analyses on dispersion and added this to the\nmanuscript. Please see lines 279-282 and 641-644. None of the distances were significant for\nasthma or atopic dermatitis at six years (p > 0.05), however the Bray-Curtis distance did exhibit\na significantly different dispersion in allergic rhinitis at six years (p = 0.0017, F = 9.35). The beta\ndispersion was tested using the function betadisper from the vegan package version 2.6-1, with\n9999 permutations.\n31. L616: This is wrong. Why are you filtering out species observed in less than half\nof your population? I do not understand the justification. Biologically how is this\njustified? Provide a graph demonstrating the relationship between Species\nfrequency of detection (occurrence), number of reads and number of contigs.\nWhat fraction of the total Anelloviruses dataset species with <50% occurrence\nrepresent? This is around 80%?\nWe understand that it may seem extreme to only analyse species that were present in >50% of\nthe samples, however this filtering threshold proved a suitable balance between multiple testing\nburden and statistical power, while still analysing the vast majority of the anellovirus composition\n(80% of the average anellovirus abundance per sample). We have described this in more detail\nin the methods, please see lines 647-649. Please read below for a thorough description of the\nrationale behind this choice, including two figures to illustrate this.\nIt is standard practice to filter out taxa in microbiome analysis, with the thresholds chosen being\ndata, methodology, and hypothesis dependent. We find it important to highlight that the low-\nprevalent and low-abundant species were considered in the alpha and beta diversity analyses\n(including the procrustes analysis), as no filtering was conducted there. Regarding the\ndifferential abundance and indicator species analysis, we decided to filter out species to reduce\nthe number of variables to test, i.e. to reduce the multiple testing burden, as these two analyses\nare statistically independent tests per individual taxa. Furthermore, the power of the statistical\ntest is reduced when many samples do not have a taxon present, further reinforcing the need to\nremove species that are only present in few samples.\nThis is also important when considering the biological interpretability of the results, identifying a\npotential significant differential abundance in a species that was only present in very few\nsamples or that had a very low relative abundance can be considered of low clinical relevance,\nas it is not likely to be generalisable to other populations. However, there are of course many\nconsiderations and trade-offs between statistical power, multiple testing burden, and potential\nclinical relevance and interpretability of the findings. Removing species that were present in less\nthan 50% of the population resulted in a dataset of 38 species out of 111, and these species\ncorresponded to on average 81% of the relative anellovirus abundance per sample. We believe\nthis is an appropriate balance between reducing the feature search space to reduce multiple\ntesting burden and maintain statistical power, while still evaluating an accurate and\nrepresentative anellovirome composition in the samples. We created two figures that illustrate\n16\nthe relationship between species prevalence, mean relative abundance, and the number of\nreads that map to these species, attached here.\nFigure 1: Anellovirus species prevalence and mean relative abundance, the dashed line\nrepresents a prevalence cut-off of 50%. This demonstrates that species present in <50% of\nsamples exhibit a very small relative abundance.\nFigure 2: Number of reads mapping to anellovirus species in sample prevalence groups,\nranging from 0-100% prevalence. The number above the bars represent the number of species\nthat belong in each prevalence group. A filtering of 50% prevalence keeps 38 of 111 species,\nand represents a large majority, 81%, of the anellovirus reads.\n17\n32. I find Indicator species analysis not informative and with many limitations. A\nmuch better way to do this analysis is to use probabilistic clustering that take into\naccount the power law distribution of metagenomic data to assign\nmetacommunities, then group samples and simply investigate within each group\nthe ratio of health versus disease or any other clinical outcome.\nWe believe employing both differential abundance and indicator species analysis provide a\nmore comprehensive understanding of the biological patterns in the data. Indicator species\nanalysis not only considers the abundance but also the consistency and exclusivity of species in\nthe samples, allowing us to find species that are both faithful (consistently present in a specific\ngroup) and exclusive (primarily found in that group). This methodology identifies species that\nmay be more ecologically relevant, despite not being the most dominant. We did try to explain\nthis in the methods, please see lines 651-659, but we are happy to expand this if needed.\nFurthermore, we believe the most comprehensive understanding of an ecosystem is achieved\nby considering multiple methodologies that are complementary to each other. However, we do\nnot mind removing the findings on indicator species analysis, if the reviewer and editor believe\nthat is necessary.\nWhile we do agree that clustering can be useful in some situations, we believe that with the\ncurrent knowledge of anelloviruses there are no suggestions that they exhibit any ecological,\nfunctional, or clinically relevant groups. Furthermore, there is no evidence that anelloviruses\nexist in a symbiotic ecology with each other, further refuting the potential usefulness of\nclustering. Considering the results of the SpiecEasi analysis we added to the manuscript,\nper your suggestion (question number 44 below), the hypothesis that anelloviruses do not exist in a\nsymbiotic ecology with each other is further solidified. We are not per se averse to conducting\nclustering, however, based on the current literature we do not believe that we have an\nappropriate ecological understanding of anelloviruses to make informed choices for this type of\nanalysis at the moment. If the editor or reviewer has concrete suggestions or references that we\ncan base this analysis on, we are open to conducting this analysis.\nClustering is a non-robust methodology that is completely dependent on the samples included,\ni.e. including or excluding a single sample may have vast impacts on the results, which limits the\ngeneralisability of such methods and findings to future research and clinical interpretation. There\nare many choices to be made when performing clustering, which we do not have answers to\nfrom the previous literature. Most clustering methods are dependent on a predefined number of\nclusters, but to this end, there is no previous data that can be used to estimate the appropriate\nnumber of clusters to investigate. We could potentially generate anywhere between 2-70\nclusters, but we do not have any previous studies providing us with a direction in this choice.\nWe do not know how many anellovirus species we can expect to find in a single cluster based\non a biological relationship. Should the clustering be performed on species or subgenus\ntaxonomy level?\nTo summarise, clustering is a simplification of an ecology, and we believe that the analysis we\nhave presented in this manuscript that are conducted on alpha and beta diversity, as well as\n18\nspecies, subgenus, genus and family level of anelloviruses are more comprehensive and\ninformative to the field than clustering which is completely data and cohort dependent and would\nbe difficult to recapitulate or replicate in future studies.\n33. L641, L663, and 683: Here you use different thresholds. There needs to be\nconsistency when using arbitrary assigned thresholds across the different tests.\nAs it is now it looks as fishing for a significant association. Describe the filtered\nsubsets? How your results change from non-filtered to 100% with increments of\n10%? What will mean if you see differences?\nFirst of all, a distinction must be made when discussing subsetting of features vs. subsetting of\nsamples. Subsetting of samples was done for cross-validation, which is a common and robust,\nwidely agreed upon practice for minimising spurious findings. The 80/20 ratio is the most widely\nused across all fields that employ statistical modelling to ensure robustness and accuracy of the\nfindings to avoid spurious findings and thus cannot be considered arbitrary.\nIn the methods of L641 (sPLS) and L663 (random forest, RF), the filtering thresholds are\nentirely consistent; all subgenera are included, and species are included if their prevalence is\n>30%. For the methods of L683 (Procrustes), no filtering is applied as this is a different\nmethodology that is based on ordinations of the dataset, not on individual taxa abundances. The\napproach to this method is similar to the other methods we applied on beta diversity analysis,\nwhich are also based on the entire composition and are not suitable to apply to filtered data.\nFor the analysis on differential abundance and indicator species analysis (DA/IS), a filtering\nthreshold was applied to the species level keeping only species that were present in >50% of\nthe population. We described the motivation for this above in our response to question number\n31, where you can also see two graphs describing the relationship between mean relative\nabundance and prevalence, and the number of reads mapping to different prevalence subsets.\nThe motivation for applying a different filter to the sPLS and RF models is based on the fact that\nthese models are inherently different from DA/IS. As we mentioned, DA/IS are independent\nstatistical tests, univariate tests, where the multiple testing burden increases drastically the\nmore taxa included for testing. Furthermore, the statistical power for low-prevalent species is\nalso decreased. sPLS and RF on the other hand are multivariable models, they do not test each\nanellovirus species separately, and thus the multiple testing burden is not applicable to these\ntests. However, these methods are highly complex and may be prone to overfitting if too many\nirrelevant features are included. Thus, we decided to apply a more lenient filtering threshold to\nthese models as a means to reduce model complexity and risk of overfitting, while still\nmaintaining interpretability and utilising the power of these models, i.e. the fact that we can\nanalyse many more species in one test. The 30% prevalence filtering applied to sPLS/RF\nresults in 64 anellovirus species being included, which corresponds to on average 93% of the\nrelative abundance per sample. This results in a multivariate analysis that encompasses a great\ncoverage of the anellovirome composition per sample, while still removing low-prevalent\nfeatures that may increase overfitting and thus reduce model interpretability and clinical\n19\nrelevance of the findings. This threshold was decided a priori based on a technical\nunderstanding of the methodologies and considering the biological and clinical interpretability of\nthe findings.\nWe believe that applying the same methodology across the full range of species prevalence\n(from 0% to 100%) may not be the most appropriate approach, especially given the distinct\nstrengths and weaknesses of each model, as we outlined earlier. Including all species in the\nanalysis would likely lead to overfitting, producing spurious results with limited clinical relevance.\nTherefore, any differences observed based on the prevalence filter are more likely to reflect\nvariations in model overfitting or underfitting, rather than biologically or clinically meaningful\nfindings. As we demonstrated in the answer to your question 31, analyses on species present in\n<30% of samples would likely lead to spurious results as they cover a very small proportion of\nthe abundance.\nTo further support this point, we have reanalyzed the sPLS model at the species level for atopic\ndermatitis at 6 years of age, using different filtering thresholds. This is summarised in the table\nbelow, demonstrating that the filtering threshold applied has minimal impact on the findings\npresented in this manuscript, and does not alter the conclusions of our study at all. This is likely\ndue to the fact that the vast majority of the anellovirus abundance remains well-represented,\neven with relatively high prevalence filtering thresholds. Thus, the anelloviruses with lower\nprevalence exhibit such minimal abundance that their contribution to the model’s predictions is\nnegligible.\nFiltering prevalence Average cover of Number of species Median AUC of\nthreshold anellovirus included in model stability runs on test\nabundance per (total 111 species) set\nsample after filtering\n50% 80% 37 0.536\n40% 88% 50 0.536\n30% 93% 64 0.538\n20% 97% 80 0.536\n10% 98% 90 0.531\n5% 99.6% 101 0.523\n0% 100% 111 0.536\n34. L691: There is no biological justification in generating a virome component of\nAnelloviruses and bacteriophages. Because of this I do not agree in generating a\ntotal relative abundance based on them alone. Especially since you have not\n20\nprovided data of how DNA isolation, MDA, and sequencing differentially affect\nbacteriophages and Anelloviruses detection. What normalisation technique\njustifies the co-analysis of bacteriophages and Anelloviruses?\nPlease see our reply to your questions number 16 and 23.\nWe have re-analysed a majority of our data with the anellovirus abundance re-normalised\nseparately, so the total relative abundance of anelloviruses per sample sums to 100%, and we\npresent the comparison of these findings in the table below. We have only presented the results\nthat were significant in either the original “complete virome” normalisation, or the “anellovirus\nonly” normalisation. The results are highlighted with green if they aligned across the methods,\nand highlighted with blue if they were not aligned across the methods. The results that could not\nbe tested with a certain normalisation were highlighted with orange.\nThe comparison presented here demonstrates several things:\n1. We can test more hypotheses with the “complete virome” i.e. normalisation with all viruses\ncompared to the “anellovirus only” normalisation\n2. The results from both normalisation approaches remain to a large extent the same\n3. Any differences would not influence the conclusions of our study\nAnalysis Result “complete virome” - Results “anellovirus only” -\nnormalised with all viruses normalised with only\nanelloviruses\nAlpha diversity OR total Fish oil intervention: Fish oil intervention:\nanellovirus abundance (TAA) - Faith’s PD: Significant in - Faith’s PD: Significant in\nand environmental variables positive direction positive direction\n(Figure 3A) - TAA: Significant in - TAA: Cannot be tested\npositive direction\nGestational age: Gestational age:\n- TAA: Significant in - TAA: Cannot be tested\nnegative direction\nDaycare: Daycare:\n- Faith’s PD: Significant in - Faith’s PD: Significant in\nnegative direction negative direction\n- TAA: Significant in - TAA: Cannot be tested\nnegative direction\nNumber of older siblings: Number of older siblings:\n- TAA: Significant in - TAA: Cannot be tested\nnegative direction\nSample date spring: Sample date spring:\n- Faith’s PD: Significant in - Faith’s PD: Significant in\npositive direction positive direction\nSample date summer: Sample date summer:\n- Faith’s PD: Significant in - Faith’s PD: Significant in\nnegative direction negative direction\n- TAA: Significant in - TAA: Cannot be tested\nnegative direction\nFUT2 secretor: FUT2 secretor:\n- TAA: Significant in - TAA: Cannot be tested\npositive direction\n- Faith’s PD: Not - Faith’s PD: Significant in\nsignificant positive direction\nBeta diversity and environmental Abx birth mother: Abx birth mother:\nvariables (Figure 3B) - Aitchison: Significant - Aitchison: Significant\nInduction: Induction:\n- Aitchison: Significant - Aitchison: Significant\n- Jaccard: Significant - Jaccard: Significant\nDaycare: Daycare:\n- Aitchison: Significant - Aitchison: Significant\n- Jaccard: Significant - Jaccard: Significant\n- Bray-Curtis: Significant - Bray-Curtis: Not\nsignificant\nSample date spring: Sample date spring:\n- Aitchison: Significant - Aitchison: Significant\n- Jaccard: Significant - Jaccard: Significant\nSample date summer: Sample date summer:\n- Jaccard: Significant - Jaccard: Significant\nBreastfeeding duration: Breastfeeding duration:\n- Jaccard: Significant - Jaccard: Significant\nAlpha diversity or TAA and Asthma: Asthma:\nclinical outcomes (Figure 4A-C) - Faith’s PD: Not - Faith’s PD: Not\nsignificant significant\n- TAA: Not significant - TAA: Cannot be tested\nRhinitis: Rhinitis:\n- Faith’s PD: Not - Faith’s PD: Not\nsignificant significant\n- TAA: Not significant - TAA: Cannot be tested\nAtopic dermatitis: Atopic dermatitis:\n- Faith’s PD: Not - Faith’s PD: Not\nsignificant significant\n- TAA: Not significant - TAA: Cannot be tested\nBeta diversity and clinical Asthma: Asthma:\noutcomes (Figure 4A-C and - Aitchison: Not significant - Aitchison: Not significant\nSupplementary Figure 3A-C) - Jaccard: Not significant - Jaccard: Not significant\n- Bray-Curtis: Not - Bray-Curtis: Not\nsignificant significant\nRhinitis: Rhinitis:\n- Aitchison: Not significant - Aitchison: Not significant\n- Jaccard: Not significant - Jaccard: Not significant\n- Bray-Curtis: Significant - Bray-Curtis: Significant\nAtopic dermatitis: Atopic dermatitis:\n- Aitchison: Not significant - Aitchison: Not significant\n- Jaccard: Not significant - Jaccard: Not significant\n- Bray-Curtis: Not - Bray-Curtis: Not\nsignificant significant\nOur decision to conduct the analyses as a “complete virome” normalisation was not unfounded.\nAlthough DNA isolation, MDA, and sequencing may differentially affect bacteriophages,\nadenoviruses and anelloviruses, that is not a justification for analysing them separately. It can\n22\nbe compared to bacteriome analyses, where it is well known that different bacteria are\ndifferentially influenced by DNA extraction and PCR amplification, however, all bacteria are\ntypically normalised to a single relative abundance per sample of 100%, as such a single\nbacterium, e.g. bacterium A, is typically not re-normalised to 100% for a single sample.\nSimilarly, we analysed the virome of these samples, where bacteriophages, anelloviruses and\nadenoviruses were identified and (compositionally) quantified together. It is generally\nconsidered of interest what the total abundance of bacterium A may be in a sample relative to\nthe abundance of all the other bacteria (e.g. having a very high/low abundance of bacterium A\ncompared to the other bacteria present in the sample). Along these lines, we hypothesised that\nthe total relative anellovirus abundance could be related to clinical or environmental outcomes,\ni.e., having more or less anelloviruses in relation to the other viruses could be indicative of\ndisease or other factors (which we demonstrated to be the case in Figures 3A and 3C, and\nsimilarly we demonstrated this to be the case for adenovirus in relation to gastroenteritis in\nFigure 1C). This is a hypothesis that could not be tested if we normalised anelloviruses\nseparately from the other viruses.\nIt can be argued that there is no reason to believe that anelloviruses and bacteriophages\ninfluence each other, but this is also true for the majority of bacteria or bacteriophages in a\nsingle sample. Not every bacteria/bacteriophage will have a direct ecological or biological\nrelationship with all other bacteria/bacteriophages. However, they are still typically analysed\ncombined as a bacteriome or bacteriophageome, similar to our approach to analysing the\nvirome.\nWe are happy to make this choice more clear in the manuscript if the reviewer or editor believes\nthat is necessary.\nTable 1:\n35. Are these correct? Explain the table contents. Controls have asthma?\nAsthma ever by 6 years (%)\nCase 140 (21.6)\nControl 487 (75.3)\nMissing 20 (3.1)\nAtopic dermatitis ever by 6 years (%)\nCase 201 (31.1)\nControl 432 (66.8)\nMissing 14 ( 2.2)\nRhinitis by 6 years (%)\nCase 41 (6.3)\nControl 576 (89.0)\nMissing 30 (4.6)\nAnswer: Thank you for pointing this out, we have now updated Table 1 to be more clear. The\ncontrols for either of the diseases do not have that particular disease, but they may have either\nor both of the other diseases.\nResults should focus on the description of Anelloviruses prevalence in the cohort, their ecology\nand genetics.\n36. Figures: Font sizes are too small to read. All figures need reformatting.\nThank you for pointing this out, the figures lost quality when submitted directly in the word\ndocument. We have now included all the figures as high-quality PDFs which are present at the\nend of the submitted document. Each panel can be zoomed in individually with high resolution.\n37. Figure 1D: thresholds are arbitrary\nIndeed, as it was presented the thresholds may appear arbitrary, however they were based on a\nprevious study that identified three classes of bacteriophages, core, common and rare (DOI:\n10.1073/pnas.1601060113), which we extended with an extra class, very rare (present in <1%\nof individuals). We have now added a reference to this manuscript, please see lines 157-160.\n38. Figure 1E: Why are you analysing subgenera and not species?\nWe did conduct analyses on all taxonomic levels (total anellovirus family abundance, genus,\nsubgenus, and species), including several alpha and beta diversity metrics in this manuscript.\nWe mainly demonstrate our findings on genus or subgenus level in figures, as these taxonomic\nlevels have an appropriate number of unique taxa for interpretability via a figure (as we show in\nFigure 1B there are 5 genera and 36 subgenera). Figure 1E is presented on subgenus, and\nSupplementary Figure 1A is presented on genus level. We do not demonstrate findings on the\nspecies level via figures, as there are 111 anelloviral species, which we deemed would not be\ninterpretable via a figure, instead, the analyses on species level are reported in the text and in\nsupplementary tables 7-9. If the reviewer wishes, we can add the information presented in\nFigure 1E and Supplementary Figure 1A on the species level as an additional supplementary\ntable.\n39. Figure 1F: Very weak associations. What does this mean?\nWe have expanded this section, please see lines 168-172. There appears to be almost no\ncompetitive relationships between the anelloviruses, but also not convincingly strong\ncorrelations, suggesting that the anelloviruses are largely independent of each other. This is\nlargely confirmed with the additional SpiecEasi analysis that we have now employed upon your\nsuggestion, which identifies very few relationships between the subgenera, and the only\nrelationships that are identified are positive. Please see our answer to your question number 43.\n40. Figure 2: should be Figure 1\nThank you for this suggestion, however, the author guidelines for this journal dictate that the\norder of the figures follow the order of mention in the text.\n41. L102: Remove bacteriophages, and any other component of the virome and focus\non Anelloviruses. These data have been published before and they do not add\nvalue in the current report.\nPlease see our reply to your questions number 23 and 33.\n24\n42. L116: This section described the presence of Anelloviruses. It should be\nexpanded to include contig information, genome coverage and what genes are\nencoded. Are all encoded proteins known? What is the role of the genes in the life\ncycle of the virus> Do we know if they interact with human signalling pathways?\nPresent the species and contig distribution data in your population. Compare\nspecies accumulation rare fraction curves in the population and in the different\nclinical groups. How much of the estimated true diversity have you captured with\nyour approach?\nThe anellovirus genome consists of several open reading frames (ORFs) that generate six\nproteins, the ORF1, ORF2, ORF1/1, ORF1/2, ORF2/2, and ORF2/3 proteins. The ORF1 protein\nis believed to have a structural role such as in ssDNA packaging. The ORF2 protein is believed\nto be both a regulatory protein that can disrupt the host immune responses, and to have a\nstructural role. The ORF2/2 protein is likely involved in regulating genome replication and gene\nexpression (DOI: 10.1093/femsre/fuaa007). Anelloviruses do interact with PAMPs, TLRs and\ninflammasomes, as well as NFkB (DOI: 10.3390/microorganisms9071357,\n10.3389/fmicb.2018.00686). Please see our answer to your question number 21 for more\ninformation about the anellovirus genes.\nAs for the anellovirus vOTU genome coverage, when an anelloviral vOTU was present in the\nsample, a mean 19% (range 6%-42%) of the anellovirus genome was covered by mapping\nreads. From vOTU to vOTU within a sample, this number ranged from 2% to 100%. These\nnumbers are as expected, as the vOTUs are assembled using reads from all samples,\nwhereafter the reads are mapped from each individual sample onto all vOTUs. We have now\nadded this to the manuscript, please see lines 124-127.\nWe plotted the density of the observed vOTUs and species in our dataset, illustrated in the\nfollowing figures.\ndeep enough to garner sufficient reads for assembly. There is no reason to believe that the\nnumber of sequencing reads / sequencing depth or the number of mapped reads would be\ndifferent between the clinical groups, which we verified (all p > 0.05). However, we did adjust all\nanalyses for sequencing depth so that any variation in read depth would be addressed if it were\nthe case.\nWe have created nonpareil/rarefaction curves of all samples in the dataset and compared the\ndiversity between the clinical groups (DOI: 10.1038/ismej.2014.76,\n10.1093/bioinformatics/btt584). This demonstrates the sequencing depth and estimated true\ndiversity of all samples, please note that this is k-mer based on the quality controlled reads, thus\nthis is the estimated true diversity of all sequenced organisms in the sample, not only for\nanelloviruses. The average coverage of the entire dataset was 0.87 (range 0.62-0.97) for all\nsamples. The coverage refers to the “coverage of the estimated diversity” in the dataset, thus a\nfraction of how much of the estimated read-based k-mer diversity is covered by the sequencing\nperformed. Thus, we capture on average 87% of the estimated diversity in the samples, across\nall organisms present. We further tested if there was any significant difference in the coverage\nbetween different disease outcomes, however none were significantly different (p > 0.05).\n26\n43. L174: Expected since the variables are linked!\nIndeed, this was expected, which is why we highlighted this including the exceptions such as\nthe subgenus tautorquevirus 1 which demonstrated a very low median relative abundance but\nstill had a high prevalence despite consisting of only one species.\n44. L178: Spearman association do not reflect ecological competition or\ncollaboration. There are appropriate packages for this type of analysis, e.g.\nhttps://github.com/zdk123/SpiecEasi\nWe completely agree that SpiecEasi could add additional value to the Spearman correlations as\na representation of anellovirus relationships, thus we have now added a SpiecEasi analysis on\nthe subgenus level which was largely concordant with our previous findings, please see lines\n170-172 and Supplementary Figure 1C. The findings from the spearman heatmap were not\nindicative of compositional effects, as largely all genera and subgenera were positively\ncorrelated. Furthermore, as opposed to Spearman which identifies non-parametric correlations,\nrepresentation of anelloviral ecology. Additionally, our anellovirus data is not particularly sparse,\nwith less than 50% of the cells of the entire dataset being zero, and less than 50% of the genera\nor subgenera being considered sparse. Therefore we have opted to keep the spearman\nheatmap in the main figure while providing the SpiecEasi correlation figure as Supplementary\nFigure 1C.\nDiscussion:\n45. L387-L398: You have not demonstrated Anellovirus infection or any other\ninfection in your cohort? So how superinfection exclusion is relevant?\nWe do not mind removing this section from the manuscript, if the reviewer and editor deems that\nis more suitable. We discuss superinfection exclusion as we detected several different\nanellovirus species present in the same samples, and there is no consensus in the anellovirus\nliterature on what is deemed infection (see our answer to your questions number 3 and 5), thus\nbelieve this was a relevant topic for discussion.\nReviewer #3 (Remarks to the Author):\nThe manuscript titled “Early Life Factors Shaping Infant Gut Anellovirus Composition and\nAssociations with Childhood Atopic Disease” by Boulund et al. focuses on factors that may\ninfluence the anelloviral community in the gastrointestinal tract of young children and explores\npotential links with atopic disease. The study identifies several associations between social,\nenvironmental, and genetic factors and the anellovirus community; however, none were found\nwith atopic disease, except for a weak association with rhinitis and allergic sensitization. The\nkey takeaway from a clinical perspective is that the child cohort under investigation is large\nenough to support the credibility of these negative findings.\nMajor Comment:\nFrom the perspective of the readership of Nature Communications, it is questionable whether\nthe presented data is comprehensive enough for publication in this forum. In its current form, the\nmanuscript primarily reports microbial data without correlating it to clinical or immune functional\ndata, which should be available in the cohort.\nWhile this is a highly interesting and important field with many unresolved questions about the\ninfluence of anellovirus on physiology and immune competence, it would be of utmost interest to\nlearn more about the potential role that anellovirus might play in the maturation of the\nmicrobiome, its possible association with host immune responses, and whether anellovirus\ncould serve as a biomarker of immune function.\nAnswer: We believe that this manuscript covers a lot more than reports on microbial data, as\nwe report several findings on clinical data, both acute and chronic outcomes, as seen in Figure\n28\n1C and Figure 4, and supplementary Figures 3-5, including reports on microbial maturation\nscore in section “The anellovirus composition is linked to the gut bacterial maturation”. The\nfocus of this manuscript was characterisation of the gut anellovirus composition at one year of\nage, as well as the relationship between anelloviruses (composition, diversity, abundance) and\nseveral host-related factors (genetic, social, environmental, early-life factors and pre-, peri-, and\npostnatal factors), in addition to the analyses on acute and chronic disease. Unfortunately, we\ndo not have immune functional data collected around the time of the faecal sample collection (1\nyear of age), thus we are not able to conduct analyses directly relating to immune function.\nHowever, we agree with the reviewer that it would have been an incredibly interesting analysis\nthat we would like to have included. To conclude, we believe that this manuscript provides a\ngreat contribution to the field of anellovirus research as it characterises the anellovirus\ncomposition in possibly the largest number of samples published to date, alongside extensive\nand robust analyses into the potential links between anelloviruses and immune development,\nmicrobiome maturation, as well as chronic and acute disease in the host.\nMinor Comments:\n1. In lines 121 and 128, please clarify for non-bioinformatician readers what the\nfuncton(\"AND\" >35%) means.\nAND is an abbreviation for average nucleotide divergence, as was written in line 121. This is a\nmeasurement of the average number of nucleotide differences between two DNA sequences,\ni.e. it is the opposite of ANI (average nucleotide identity, AND = 100-ANI), which may be a more\nknown concept. We have now changed the cutoffs to ANI instead of AND and explained it\nbetter, please see lines 135-137.\n2. Do the anellovirus profiles in the children resemble those of their mothers?\nUnfortunately, we do not have virome sequencing data (or faecal samples) of the mothers in this\ncohort, so we are unable to test this otherwise very interesting and relevant hypothesis.\n3. Line 294: The results in the preceding lines were also not significant. Therefore, it\nis not accurate to state in this sentence that “no other results were significant.”\nYou are correct, we have now removed the reports of the suggestive results on sensitization\nand simply report that none of the results on sensitization were significant, please see lines\n265-270.\ncorresponding responses below, in blue:\nReviewer #1 (Remarks to the Author):\nThe authors have made a reasonable attempt to address the reviewer's comments.\nWe appreciate that the reviewer acknowledges our effort.\nThere is still some lack of detail or clear explanation of some methodological elements. In\nsome cases methods are in the results, which complicates things.\nThanks for these suggestions. We have addressed these concerns to improve the clarity.\nExample issues with methods:\n-sample collection and dna extraction for the virome analysis is not clear. How are samples\nhandled before glycerol and -80 storage? How is enrichment carried out with the sample in\nglycerol? What kit is used to extract dna from the enriched samples? Kits HAVE been given\nfor library construction and extraction of DNA from bacteria for community analysis; not for\nthe virus prep, though.\nFecal samples were collected 1 year after birth, either at the research clinic or by the parents\nat home using detailed instructions. Thereafter the samples were sent to the laboratory by\npost. Each sample was mixed on arrival with 1 mL of 10% vol/vol glycerol broth (Statens\nSerum Institut, Copenhagen, Denmark) and frozen at −80 °C. The faecal sample was\nsuspended in a buffer containing 10% glycerol, so this did not interfere with enrichment. The\nQIAmp viral RNA mini kit (Qiagen, Hilden, Germany) was used for viral DNA/RNA extraction.\nWe have now added these details to the corresponding methods section. Please see lines\n467 to 472 within the revised manuscript version.\n-the CONTIG deduplication is explained, but not the READ deduplication method\nThank you for pointing this out. The workflow has now been expanded, please see lines 475\nto 480 in the revised version.\n-the move to protein sequence for anellovirus clustering is not clear. Anelloviruses are not\nthe only virus group with sequence diversity but gene conservation (surely this is also seen\nin the phage population). However, even if protein clustering is justifiably only appleid to\nanelloviruses, the methods by which the proteins sequences are clustered are not clear -\nmethods and cutoffs? When is nucleotide sequence used vs. protein?\nIndeed, the reviewer raises an excellent point that has not been straightforward to\ncommunicate. Our already-published bacteriophage work on the same samples utilised DNA\nsimilarity to identify viral species, and protein similarity to identify viral families and determine\nviral gene functions. Here we leveraged that existing protein data to select the anelloviruses\nfor this study, based on protein similarity to anelloviral ORF1. The advantage of this\napproach is its sensitivity to any anelloviral “dark matter”, not represented in databases. A\nDNA-similarity search would only have identified anelloviruses similar to those already found\nin databases. Next, for delineating anelloviral taxonomy, we had to revert back to DNA\nsimilarity, because the ICTV cutoffs for defining anelloviral subspecies, species and genera\nare given in DNA similarity terms.\nWe acknowledge that this back and forth seems convoluted, but we do believe it was the\nmost effective way to ensure reference-independent sensitivity and taxonomic resolution. In\nan attempt to improve clarity, we have now modified the methods section at lines 480 to 490,\nand hope that this better conveys the approach and rationale.\n- The figures are clearer but the font is still very tiny.\nThank you for this suggestion. We have tried our best to improve legibility now.\n- Figure 2 - what is the difference between dark and light blue? Why are tau and alpha\ngroups clustered together?\nThis is because we chose to color the terminal tips in alternating shades to denote subgenus\nboundaries within each genus. This explanation is now given in the figure legend. Thanks for\npointing this out. As for Tau and Alpha, we can see how those two colors (“maroon” and red)\nwere too similar. We’ve now made Tau purple instead, which makes it stand out better.\nHope this helps!"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have made a reasonable attempt to address the reviewer's comments. There is still some lack of detail or clear explanation of some methodological elements. In some cases methods are in the results, which complicates things.\n\nExample issues with methods:\n-sample collection and dna extraction for the virome analysis is not clear. How are samples handled before glycerol and -80 storage? How is enrichment carried out with the sample in glycerol? What kit is used to extract dna from the enriched samples? Kits HAVE been given for library construction and extraction of DNA from bacteria for community analysis; not for the virus prep, though.\n-the CONTIG deduplication is explained, but not the READ deduplication method\n-the move to protein sequence for anellovirus clustering is not clear. Anelloviruses are not the only virus group with sequence diversity but gene conservation (surely this is also seen in the phage population). However, even if protein clustering is justifiably only appleid to anelloviruses, the methods by which the proteins sequences are clustered are not clear - methods and cutoffs? When is nucleotide sequence used vs. protein?\n- The figures are clearer but the font is still very tiny.\n- Figure 2 - what is the difference between dark and light blue? Why are tau and alpha groups clustered together?", "reviewer_3": "No further comments."}, "author_response": "We thank the responding reviewer for their final comments. Please find their queries and our corresponding responses below, in blue:\n\nReviewer #1 (Remarks to the Author):\nThe authors have made a reasonable attempt to address the reviewer's comments.\nWe appreciate that the reviewer acknowledges our effort.\n\nThere is still some lack of detail or clear explanation of some methodological elements. In\nsome cases methods are in the results, which complicates things.\nThanks for these suggestions. We have addressed these concerns to improve the clarity.\n\nExample issues with methods:\n-sample collection and dna extraction for the virome analysis is not clear. How are samples\nhandled before glycerol and -80 storage? How is enrichment carried out with the sample in\nglycerol? What kit is used to extract dna from the enriched samples? Kits HAVE been given\nfor library construction and extraction of DNA from bacteria for community analysis; not for\nthe virus prep, though.\nFecal samples were collected 1 year after birth, either at the research clinic or by the parents\nat home using detailed instructions. Thereafter the samples were sent to the laboratory by\npost. Each sample was mixed on arrival with 1 mL of 10% vol/vol glycerol broth (Statens\nSerum Institut, Copenhagen, Denmark) and frozen at −80 °C. The faecal sample was\nsuspended in a buffer containing 10% glycerol, so this did not interfere with enrichment. The\nQIAmp viral RNA mini kit (Qiagen, Hilden, Germany) was used for viral DNA/RNA extraction.\nWe have now added these details to the corresponding methods section. Please see lines\n467 to 472 within the revised manuscript version.\n\n-the CONTIG deduplication is explained, but not the READ deduplication method\nThank you for pointing this out. The workflow has now been expanded, please see lines 475\nto 480 in the revised version.\n\n-the move to protein sequence for anellovirus clustering is not clear. Anelloviruses are not\nthe only virus group with sequence diversity but gene conservation (surely this is also seen\nin the phage population). However, even if protein clustering is justifiably only appleid to\nanelloviruses, the methods by which the proteins sequences are clustered are not clear -\nmethods and cutoffs? When is nucleotide sequence used vs. protein?\nIndeed, the reviewer raises an excellent point that has not been straightforward to\ncommunicate. Our already-published bacteriophage work on the same samples utilised DNA\nsimilarity to identify viral species, and protein similarity to identify viral families and determine\nviral gene functions. Here we leveraged that existing protein data to select the anelloviruses\nfor this study, based on protein similarity to anelloviral ORF1. The advantage of this\napproach is its sensitivity to any anelloviral “dark matter”, not represented in databases. A\nDNA-similarity search would only have identified anelloviruses similar to those already found\nin databases. Next, for delineating anelloviral taxonomy, we had to revert back to DNA\nsimilarity, because the ICTV cutoffs for defining anelloviral subspecies, species and genera\nare given in DNA similarity terms.\nWe acknowledge that this back and forth seems convoluted, but we do believe it was the\nmost effective way to ensure reference-independent sensitivity and taxonomic resolution. In\nan attempt to improve clarity, we have now modified the methods section at lines 480 to 490,\nand hope that this better conveys the approach and rationale.\n\n- The figures are clearer but the font is still very tiny.\nThank you for this suggestion. We have tried our best to improve legibility now.\n\n- Figure 2 - what is the difference between dark and light blue? Why are tau and alpha\ngroups clustered together?\nThis is because we chose to color the terminal tips in alternating shades to denote subgenus\nboundaries within each genus. This explanation is now given in the figure legend. Thanks for\npointing this out. As for Tau and Alpha, we can see how those two colors (“maroon” and red)\nwere too similar. We’ve now made Tau purple instead, which makes it stand out better.\nHope this helps!"}]

Decision letter