Виром ДНК связан с аутоиммунными заболеваниями и COVID-19.

Main

Well-regulated induction of immune responses targeting nonself-viral epitopes is beneficial for virus clearance. On the other hand, excessive or aberrant responses to virus infection can resemble those of autoimmune diseases and occasionally pose risks of self-reactive immune responses as in Guillain–Barré syndrome 1 . Similarly, epitope indistinguishability between EBNA1, encoded by Epstein–Barr virus (EBV), and host GlialCAM likely contributes to the autoreactive lymphocyte expansion seen in multiple sclerosis (MS) 2 , 3 . In the mouse model, persistent infection by Theiler’s virus causes an autoimmune response targeting myelin via epitope spreading, where self-proteins are aberrantly presented to the immune system 4 . In systemic lupus erythematosus (SLE) plasma, antibody titers against an envelope protein encoded by endogenous retrovirus correlate with interferon (IFN)-stimulated genes (ISGs) expression 5 . In addition, the immunopathogenic response recognized as a ‘cytokine storm’ is attributable to both virus infection, as seen in severe coronavirus disease 2019 (COVID-19), and autoimmune diseases 6 . Our understanding of the relationships between virus infection and autoimmune diseases is incomplete. For example, it is still unclear which autoimmune diseases are associated with virus infection and what factors license virus infection to initiate and exacerbate autoimmune disease.

DNA viruses often establish asymptomatic latent infections in healthy individuals. For example, EBV infection is highly prevalent, reaching more than 90% of adult individuals 7 . Anelloviruses are detectable in the blood of more than 8% of healthy people 8 . Beyond exogenous viruses, about 1% of humans carry endogenous HHV-6 (eHHV-6), a heritable HHV-6 genome integrated into the subtelomere region that has been proposed to be able to reactivate as a transmissible virus 9 , 10 . Because eHHV-6 is not included in reference human genomes, it is an overlooked genetic factor that may trigger antiviral responses. High-throughput sequencing has provided new methods to identify and study such viruses, characterizable as the human ‘virome’ 8 . For example, whole-genome sequencing (WGS) of human specimens and cell lines has been used to study eHHV-6 prevalence and potential reactivation 9 , 11 , 12 . Population-scale sequencing can suggest relationships between the presence of specific human viruses and diseases or immune responses. For instance, the presence of EBV nucleic acid in the spleen and blood is associated with an increased plasma cell number 13 . The presence of an ssDNA virus family discovered by virome profiling, called Redondoviridae , is associated with periodontitis 14 . Anellovirus load correlates with immunosuppression and acute respiratory disease 15 , 16 . The presence of eHHV-6 has been associated with pre-eclampsia, angina pectoris and increased levels of HHV-6-specific antibody titers 17 – 19 . These findings highlight the need for systematic virome analysis to identify new associations between viruses and human diseases. The partially revealed associations between viruses and autoimmune diseases, for example, between EBV and MS, suggest that similar relationships may be discoverable with further investigation.

In this study, we analyzed WGS from >6,000 Japanese individuals, including patients with five autoimmune diseases, COVID-19 or healthy controls. We systematically detected eHHV-6 and anellovirus nucleic acids in these datasets and conducted case–control association studies that compared eHHV-6 prevalence or anellovirus abundance in unrelated individuals with or without these diseases (Fig. 1 ). To confirm the inherited nature of the detected eHHV-6 genome, we performed genome-wide association study (GWAS) of eHHV-6B positivity and PacBio HiFi long-read sequencing of eHHV-6B-positive patients with SLE. In addition, epitope mapping by phage immunoprecipitation sequencing (PhIP-seq) from four eHHV-6B-positive and 12 eHHV-6B-negative patients with SLE and single-cell RNA sequencing (scRNA-seq) data from four eHHV-6B-positive and 11 eHHV-6B-negative patients with SLE were compared. Fig. 1 Overview of the study design. We collected 6,321 WGS datasets from patients with five autoimmune diseases (PV, RA, SLE, PAP and MS) and COVID-19 and healthy controls. The datasets were analyzed to detect eHHV-6 and anellovirus infection. In addition, we conducted a large-scale replication study of the association between SLE and eHHV-6 using All of Us WGS data. We conducted GWAS of HHV-6B positivity. We performed long-read sequencing, PhIP-seq and scRNA-seq targeting eHHV-6B-positive patients with SLE ( Methods ).

Results

Overview

We analyzed 6,321 Japanese samples. Details of the sample characteristics, sample collection and the WGS platform are described in Methods . The WGS datasets were composed of 1,357 participants diagnosed with psoriasis vulgaris (PV), 378 with rheumatoid arthritis (RA), 251 with SLE, 216 with pulmonary alveolar proteinosis (PAP), 36 with MS, 1,164 with COVID-19 and 2,919 healthy controls. The majority of the patients with COVID-19 ( n = 1,068; 91.8%) had severe disease (requiring oxygen support, artificial respiration and/or intensive care; Supplementary Table 1 ). Two bioinformatics pipelines were used to detect eHHV-6 and anellovirus infection (Fig. 1 ). Briefly, we remapped those WGS reads that did not map to the human reference genome to their respective viral reference genomes. For anellovirus detection, 3,333 human anellovirus sequences were used to accurately detect widely divergent anellovirus genomes. Logistic regression and Fisher’s exact were used to calculate the adjusted odds ratio (OR) for the risk of each disease with eHHV-6 or three anellovirus infection states defined from normalized anellovirus viral load (viral copies per cell). In each disease with a high OR, we estimated the relationship between clinical metrics in each disease and eHHV-6 or three anellovirus infection states using linear regression or logistic regression and Fisher’s test.

eHHV-6B associates with autoimmune diseases and SLE activity

We identified eHHV-6A in 16 participants and eHHV-6B in 22 participants (0.25% and 0.35% prevalence, respectively; Table 1 , Supplementary Table 2 and Fig. 2a ). Of these participants, one healthy control with eHHV-6A and four participants with eHHV-6B (one with PV, two with PAP and one healthy control) had evidence of the ‘solo-DR’ form of eHHV-6 (ref. 9 ). Table 1 Associations of eHHV-6A/eHHV-6B with autoimmune diseases and COVID-19 in Japanese individuals Disease Number of total samples eHHV-6A eHHV-6B No. (%) Logistic regression Fisher’s test No. (%) Logistic regression Fisher’s test OR 95% CI Wald test P Mid- P value OR 95% CI Wald test P Mid- P value PV 1,357 2 0.15 0.46 0.07–1.93 0.34 0.22 6 0.44 2.2 0.62–7.47 0.22 0.20 RA 378 1 0.26 0.55 0.03–3.62 0.60 0.82 1 0.26 1.5 0.07–10.9 0.71 0.76 SLE 251 0 0.0 NA (quasi-complete separation) 0.40 4 1.6 6.9 1.61–27.2 0.0061 a 0.0060 a PAP 216 0 0.0 NA (quasi-complete separation) 0.46 3 1.4 7.0 1.43–27.7 0.0075 a 0.022 MS 36 0 0.0 NA (quasi-complete separation) 0.87 1 2.8 10.9 0.55–71.6 0.034 0.085 COVID-19 1,164 2 0.17 0.61 0.09–2.50 0.54 0.31 1 0.086 0.56 0.03–3.77 0.61 0.46 Healthy controls 2,919 11 0.38 1.0 – – – 6 0.21 1.0 – – – Overall 6,321 16 0.25 – – – – 22 0.35 – – – – Wald test and Fisher’s exact test were conducted as two-sided tests. a Wald test P value < 0.05/6 and adjusted OR > 1, or Fisher’s test mid- P value < 0.05/6 and greater than healthy controls. Fig. 2 Prevalence of eHHV-6 and SLEDAI scores for eHHV-6B-positive patients with SLE. a , Scatter plot of the prevalence of eHHV-6A versus eHHV-6B in participants with five autoimmune diseases, COVID-19 and healthy controls. The marker size indicates the sample size for each disease. The red line indicates that the prevalence of eHHV-6 = 1%. b , Box plot of SLEDAI scores for patients with SLE without and with eHHV-6B. SLEDAI scores in participants with eHHV-6B were significantly higher than those without eHHV-6B (mean SLEDAI = 30.5 and 6.0 for eHHV-6B-positive and eHHV-6B-negative patients with SLE, respectively; two-sided Wald test P = 1.3 × 10 −8 ). The histogram shows the distribution of the number of patients with SLE without and with eHHV-6B (top). The horizontal histogram shows the distribution of SLEDAI scores of patients with SLE without and with eHHV-6B, with different colors (right). Boxes denote the IQR; the median is shown as horizontal bars; whiskers extend to 1.5× IQR; outliers are shown as individual points. IQR, interquartile range.

Next, we tested for case–control associations between eHHV-6 and each disease using the logistic regression model and Fisher’s test. No significant differences were observed in eHHV-6A prevalence in patients with any disease compared with healthy controls. In contrast, eHHV-6B was more prevalent in patients with SLE and PAP than in healthy controls (1.6% and 1.4% versus 0.21% prevalence; Fisher’s test mid- P = 0.0060 and 0.022, respectively). After adjustment for age and sex, we confirmed eHHV-6B associations with SLE and PAP with an adjusted OR of 6.9 in SLE (95% confidence interval (CI) = 1.61–27.2; Wald test P = 0.0061) and 7.0 in PAP (95% CI = 1.43–27.7; P = 0.0075). These observations suggest that carrying eHHV-6B contributes to the risk of these diseases with a large effect size relative to previously known genetic and environmental factors of the diseases 20 – 22 .

To identify a possible association between eHHV-6B and clinical metrics of SLE or PAP, we used the linear or logistic regression model and Fisher’s test (Supplementary Tables 1 and 3 ). eHHV-6B had no association with low serum complement levels, immunosuppressant use or prednisolone use in SLE ( P = 0.46, 0.79 and 0.76, respectively). However, eHHV-6B had a significant positive correlation with high SLE disease activity index (SLEDAI) scores 23 (mean SLEDAI = 30.5 and 6.0 for eHHV-6B-positive and eHHV-6B-negative patients with SLE, respectively; P = 1.3 × 10 −8 ). Three of four participants with SLE carrying eHHV-6B had very high SLEDAI scores (≥18), suggesting that eHHV-6B may have a substantial clinical impact on disease activity (Fig. 2b ). Of interest, three of four eHHV-6B-positive participants with SLE were positive for anti-RNP antibodies, which are typically seen in about 33% of patients with SLE 24 , while the fourth participant was positive for anti-Scl-70 antibodies, which is quite rare in SLE. Two eHHV-6B-positive patients with SLE developed complications from other autoimmune diseases. One patient who tested positive for anti-RNP antibodies had scleroderma and Sjögren’s syndrome, and the other patient positive for anti-Scl-70 antibodies had scleroderma and antiphospholipid syndrome. Given these unusual clinical manifestations, inborn eHHV-6B positivity may predispose to atypical SLE pathogenesis. eHHV-6B was not associated with the titer of antigranulocyte-macrophage colony-stimulating factor autoantibodies in PAP ( P = 0.90), which are considered central to the pathogenesis of the adult sporadic type of PAP 25 .

eHHV-6B also associates with high SLE risk in All of Us

We conducted an independent large-scale replication study of the association between SLE and eHHV-6 using WGS data from the National Institutes of Health’s (NIH) All of Us program 26 . The participants in the All of Us cohort are genetically diverse, and the risk of SLE varies considerably by ancestry, so we performed association testing in the 121,494 individuals of European ancestry and the entire 243,247-person multi-ancestry dataset using a logistic regression model (Table 2 ). Concordant with the Japanese discovery cohort, no significant associations were observed in eHHV-6A prevalence in either population ( P > 0.15). In contrast, eHHV-6B was associated with high SLE risk, with an adjusted OR of 2.7 in European ancestry (95% CI = 1.54–4.76; P = 0.0020) and 2.4 in multi-ancestry (95% CI = 1.47–3.82; P = 0.0012). These results indicate that eHHV-6B is associated with SLE in individuals of diverse ancestry. Table 2 Replication studies of eHHV-6A/eHHV-6B associations with SLE in All of Us Disease Frequency of the participants (%) eHHV-6A eHHV-6B No. OR 95% CI Score test P No. OR 95% CI Score test P European ancestry ( n = 121,494) SLE <0.6 ≤20 a 2.3 0.90–5.79 0.15 ≤20 a 2.7 1.54–4.76 0.0020 b Controls >99.4 333 – – – 937 – – – Multi-ancestry ( n = 243,247) SLE <0.7 ≤20 a 1.4 0.54–3.51 0.68 ≤20 a 2.4 1.47–3.82 0.0012 b Controls >99.3 502 – – – 1,578 – – – Score test was conducted as a two-sided test. a Counts suppressed per the All of Us Data and Statistics Dissemination Policy. b Score test P value < 0.05/2 and an adjusted OR > 1.

Risk eHHV-6B of SLE locates on chromosome (chr) 22q

The sequencing depth of the HHV-6 genome detected in this study was nearly half that of the human genome in each participant, suggesting that the HHV-6 sequences are inherited along with single human haplotypes (Extended Data Fig. 1a,b ). To understand the location of eHHV-6B in the human genome, we conducted a GWAS using common variants called the T2T-CHM13 reference. Comparison of 22 eHHV-6B-positive participants and 216 eHHV-6B-negative participants revealed a highly significant association with variants on the distal end of chr22q (Fig. 3a ). This result was consistent with a previous GWAS using Fisher’s test 9 , but using T2T-CHM13 as a reference provided a more densely mapped set of associated variants. Focusing on the distal end of chr22q, we performed logistic regression including rare variants (Fig. 3b,c and Supplementary Table 4 ). The 16 variants with P < 1 × 10 −4 spanned a 500-kb region including different gene regions and different linkage disequilibrium (LD) block structures, but all were in high LD with eHHV-6B. By using 30× PacBio HiFi long-read sequencing, we confirmed that these tagged variants were on the same haplotype in each of the four eHHV-6B-positive patients with SLE (Fig. 3d and Extended Data Fig. 1c ). These GWAS results reflect the LD relationship with the HHV-6B genome, suggesting that eHHV-6B in the Japanese population is integrated to the distal end of chr22q and has been inherited together as a linked unit. All of the tagged variants are intergenic or intronic, thus not supporting their pathogenicity, whereas all four eHHV-6B-positive patients with SLE carried the full-length HHV-6B genome implicated with viral reactivation 27 (Extended Data Fig. 1b ). eHHV-6B itself rather than the tagged variants might therefore be involved in the pathogenesis of SLE. Fig. 3 eHHV-6B in Japanese integrated into chr22q. a , Manhattan plot of common variant-based GWAS of eHHV-6B (22 cases versus 216 controls). The straight line indicates the genome-wide significance level (5 × 10 −8 ). b , Regional association plot for the distal end of chr22q. The dotted line indicates the genome-wide significance level. The 16 variants with P < 1 × 10 −4 were labeled by reference SNP ID (rs ID). Color shading represents the maximum value of LD ( r 2 ) with each of the lead variants in East Asian populations of 1KG. The circle size indicates LD ( r 2 ) with the HHV-6B genome in the present WGS data. c , LD map of the lead variants in East Asian populations of 1KG. d , Heatmap representing the dominant genotype of the 16 variants in 22 eHHV-6B-positive participants. We confirmed that the variants were on the same haplotype in four patients with SLE with eHHV-6B by 30× HiFi long-read sequencing (Extended Data Fig. 1c ). Uncorrected P values from the GWAS analysis (two-sided Wald test) are shown in a and b . LRS, long-read sequencing; EAS, East Asian.

eHHV-6B results in distinctive immunity in patients with SLE

To investigate whether eHHV-6 can result in the induction of unique antiviral responses, we performed epitope mapping of anti-HHV-6 antibodies in 16 patients with SLE (four eHHV-6B-positive and 12 eHHV-6B-negative; Supplementary Table 5 ) by PhIP-seq (Fig. 4a ) 28 , 29 . As a result of epitope mapping to 108 HHV-6B proteins, two HHV-6B virion proteins, antigenic virion protein and glycoprotein O, were frequently recognized as antigens across all participants, in agreement with the established notion that most adults have at one time been infected by exogenous HHV-6 and have acquired immunity against it 30 , 31 . To test whether there are epitopes specific to eHHV-6B carriers, we tested for an association between the presence of eHHV-6B and the antigenicity of HHV-6 peptides (Fig. 4b ). Responses against one region (amino acids 476–490) in the immediate-early A (IE-A) transactivator were strongly detected in all eHHV-6B carriers, whereas responses against this region were not as strong in all eHHV-6 noncarriers ( P = 3.5 × 10 −6 ; Fig. 4c ). This demonstrates that eHHV-6B results in distinctive immune responses against HHV-6B in patients with SLE. Because eHHV-6B originated as a virus but is now encoded in the human genome—it has the nature of both virus and self—eHHV-6B may be a heritable viral antigen to which immune responses represent autoimmunity, a phenomenon we term ‘endoimmunity’. Fig. 4 HHV-6 epitope mapping of antibodies in patients with SLE. a , Detection of epitopes on HHV-6 proteins. Means of the PhIP-seq scores of the indicated participant groups are shown. Each dot shows a score for a peptide split by a 14 amino acids sliding window applied to HHV-6A and HHV-6B proteins. Specific proteins of interest are labeled. b , Detection of epitopes on the three HHV-6B proteins. Copy numbers of immunoprecipitated phage genomes encoding HHV-6B peptides are shown as a heatmap. Three HHV-6B proteins that harbor abundant epitopes are shown. The region of the IE-A transactivator associated with eHHV-6B positivity is labeled with an orange arrow (amino acids 476–490). c , Statistical association between the presence of eHHV-6 and the copy numbers of immunoprecipitated phage genomes encoding HHV-6 peptides. IE-A transactivator amino acids 476–490 are labeled in orange. Uncorrected P values from the two-sided Wald test are shown. Red horizontal lines show the threshold of P (<0.05) after Bonferroni correction. IS, immunosuppressant use; NC, no antibody control; M, male; F, female.

eHHV-6B induces antiviral response in SLE monocytes

Participants with both eHHV-6 and SLE showed evidence of distinctive humoral immunity to HHV-6 peptides. As SLE is characterized by immune cell dysfunction, we sought to further investigate the immunopathogenesis in eHHV-6B-positive participants with SLE. Peripheral blood mononuclear cells (PBMCs) were collected from the same individuals who were investigated by PhIP-seq (4 eHHV-6B-positive patients with SLE as cases and 11 eHHV-6B-negative patients with SLE as controls; Supplementary Table 5 ).

After processing using a unified single-cell analysis pipeline ( Methods ), we obtained 66,915 high-quality cells from PBMCs of 4 cases and 11 controls. There was no detectable expression of HHV-6B genes in any of the cells. Cells were annotated using Azimuth 32 (Fig. 5a ). To identify more differential expression (DE) patterns than the pseudobulk methods or the mixed models, we applied Milo 33 and identified 2,818 neighborhoods, none of which showed evidence of differential abundance (DA; spatial false discovery rate (FDR) < 0.05), and 20 neighborhood groups consisting of neighborhoods that shared at least five cells (Fig. 5b and Extended Data Fig. 2a ). After excluding several cell groups to avoid sampling bias, DE and Gene Ontology (GO) enrichment analyses were performed (Fig. 5c and Extended Data Fig. 2b–e ). GO analyses showed that antiviral responses were enriched in cells in groups 6, 7 and 18, which were annotated as monocytes and dendritic cells (DCs). To evaluate the IFN response to viruses, we calculated antiviral ISG scores of cells in each neighborhood and compared them between cases and controls. In groups 6, 7 and 18, the antiviral ISG scores of the cases were significantly higher than those of the controls (FDR < 0.05; Fig. 5d ). Within-group comparisons showed that the mean scores in the cases were significantly higher than those in the controls in groups 6, 7 and 18, and the differences were larger than those in the other groups (Fig. 5e ). Differences in the antiviral ISGs scores described above were observed even between eHHV-6B-positive and eHHV-6B-negative patients with immunosuppressant use, suggesting that these increased antiviral responses in monocytes and DCs were not attributable to immunosuppressants (Extended Data Fig. 3a,b ). These results are consistent with the distinctive immunity to HHV-6B between eHHV-6B carriers and noncarriers revealed by the epitope mapping. Fig. 5 scRNA-seq of PBMCs from eHHV-6B-positive and eHHV-6B-negative patients with SLE. a , UMAP embedding of scRNA-seq data for all 66,915 PBMCs from the participants with SLE ( n = 4 for eHHV-6B-positive and 11 for eHHV-6B-negative). Thirty cell types were defined by RNA expression of marker genes using Azimuth. b , Graphical representation of neighborhoods identified by Milo. Nodes are neighborhoods, colored by neighborhood groups, and their sizes correspond to the number of cells in each neighborhood. Graph edges depict the number of cells shared between adjacent neighborhoods. c , Enriched biological processes by GO analysis of DE genes of each neighborhood group. The total number of DE genes in each group is indicated in brackets. Dot color indicates the statistical significance of the enrichment (adjusted P values via the Benjamini–Hochberg method), and dot size represents the gene ratio annotated to each term. d , Graphical representation of neighborhoods identified by Milo as in b . Node color indicates the difference between the mean of the antiviral ISG scores of cases and the mean of the scores of controls within the neighborhood (spatial FDR < 0.05), comparing eHHV-6B-positive patients with SLE as cases ( n = 4) and eHHV-6B-negative patients with SLE as controls ( n = 11). Neighborhoods with nonsignificant differences (spatial FDR ≥ 0.05) are colored white. Neighborhoods that contain only cells of either cases or controls are colored black. e , Comparison of differences within each neighborhood group shown in d using a two-sided paired t test. Diff is the average difference between the case and control means. Boxes denote the IQR; the median is shown as horizontal bars; whiskers extend to 1.5 times the IQR. *FDR < 1.0 × 10 −10 and **FDR < 1.0 × 10 −50 . ASDC, Axl+ Siglec-6+ dendritic cell; pDC, plasmacytoid dendritic cell; HPSC, hematopoietic stem and progenitor cell; T reg , regulatory T cell; NK, natural killer; T EM , effector memory T cell; ILC, innate lymphoid cell; MAIT, mucosal-associated invariant T; GDT, γδ T cells; CTL, cytotoxic T cell; P adj , adjusted P value; Nhood, neighborhood.

Anellovirus associates with autoimmune diseases and COVID-19

Anellovirus was analyzed in 5,597 Japanese WGS samples, excluding participants whose WGS was conducted using DNA extracted from cell lines. The anellovirus load was estimated and stratified ( Methods ). Of the 5,597 participants, 689 anellovirus-positive participants (12.3%; that is, those with anellovirus load above the detection limit), 57 participants (1.0%) with anellovirus load >1.0 considered anellovirus-infected and 16 participants (0.29%) with anellovirus load >8.0 considered anellovirus-viremic were found (Table 3 , Supplementary Table 6 and Extended Data Fig. 4 ). Table 3 Associations of anellovirus infection and viremia with autoimmune diseases and COVID-19 in Japanese individuals Disease Number of total samples Anellovirus infection Anellovirus viremia No. (%) Logistic regression Fisher’s test No. (%) Logistic regression Fisher’s test OR 95% CI Wald test P Mid- P value OR 95% CI Wald test P Mid- P value PV 1,062 7 0.66 4.0 1.01–19.8 0.061 0.025 0 0.0 – – – 1.0 RA 378 14 3.7 15.4 4.89–68.2 2.5 × 10 −5a 7.4 × 10 − 10a 3 0.79 NA (quasi-complete separation) 0.0023 a SLE 251 12 4.8 58.5 15.6–293 2.6 × 10 − 8a 3.6 × 10 − 10a 6 2.4 NA (quasi-complete separation) 5.6 × 10 − 7a PAP 216 0 0.0 NA (quasi-complete separation) 0.72 0 0.0 – – – 1.0 MS 36 0 0.0 NA (quasi-complete separation) 0.94 0 0.0 – – – 1.0 COVID-19 1,164 20 1.7 6.2 2.04–26.7 0.0042 a 2.9 × 10 − 7a 7 0.60 NA (quasi-complete separation) 3.3 × 10 − 4a Healthy controls 2,490 4 0.16 1.0 – – – 0 0.0 1.0 – – – Overall 5,597 57 1.0 – – – – 16 0.29 – – – – Anellovirus infection, participants with anellovirus load >1.0; anellovirus viremia, participants with anellovirus load >8.0. a Wald test P value < 0.05/6 and an adjusted OR > 1 or Fisher’s test mid- P value < 0.05/6 and greater than healthy controls.

Correlations between the three anellovirus detection states and each disease were assessed using the logistic regression model and Fisher’s test. The proportions of anellovirus infection were higher in patients with SLE, RA and COVID-19 than in healthy controls (4.8%, 3.7% and 1.7% versus 0.16%; mid- P = 3.6 × 10 −10 , 7.4 × 10 −10 and 2.9 × 10 −7 , respectively), and the proportions of anellovirus viremia were higher in patients with SLE, RA and COVID-19 than in healthy controls (2.4%, 0.79% and 0.61% versus 0%; mid- P = 5.6 × 10 −7 , 0.0023 and 3.3 × 10 −4 , respectively). Anellovirus infection associated positively with SLE, RA and COVID-19, with an adjusted OR of 58.5 in SLE (95% CI = 15.6–293; P = 2.6 × 10 −8 ), 15.4 in RA (95% CI = 4.89–68.2; P = 2.5 × 10 −5 ) and 6.2 in COVID-19 (95% CI = 2.04–26.7; P = 0.0042). Again, the quantitative magnitudes of the observed OR were larger than those previously reported for genetic and environmental factors of the diseases.

Next, linear or logistic regression models and Fisher’s test were used to identify a possible association between anellovirus detection states and clinical metrics of SLE, RA and COVID-19 (Supplementary Tables 1 and 3 ). In SLE, anellovirus was positively associated with low serum complement levels, immunosuppressant use and prednisolone use, with an adjusted OR of 4.4 in participants with low serum complement levels (95% CI = 1.46–19.5, P = 0.020), 3.0 in participants treated with immunosuppressants (95% CI = 1.51–6.03, P = 0.0020) and 2.2 in participants treated with prednisolone (95% CI = 1.05–5.21, P = 0.047; Extended Data Fig. 5a–d ). Anellovirus detection states were not associated with SLEDAI scores. Although anellovirus detection states were insignificantly associated with rheumatoid factor (RF) or anticitrullinated protein antibody (ACPA), all participants with RA who were viremic for anellovirus ( n = 3) were RF- and ACPA-positive. Patients with severe COVID-19 ( n = 1,068) were not statistically at higher risk of anellovirus infection, but all anellovirus-infected participants with COVID-19 ( n = 20) had severe disease. The proportion of participants with anellovirus infection was statistically higher in participants who died of COVID-19 than in those who survived (mid- P = 0.033).

Discussion

Here 6,321 Japanese WGS datasets were used to test for associations between six immune-mediated diseases and two constituents of the human blood virome, eHHV-6 and anellovirus. We present robust evidence of positive associations between these viruses and several autoimmune diseases and COVID-19. The observed risks of the diseases were relatively larger than those of the previously known genetic and environmental factors of the diseases, although their prevalence remains rare, suggesting that these viruses substantially impact disease biology and are relevant to clinical medicine.

We used WGS data, a robust method for detecting eHHV-6 (ref. 11 ), to estimate the higher prevalence of eHHV-6B in patients with SLE and PAP. eHHV-6 prevalence among healthy controls was 0.59%, similar to those in previous studies of healthy Japanese participants 34 , 35 . Our previous survey of eHHV-6 in a disease-oriented Japanese biobank found a higher prevalence of HHV-6B than HHV-6A, the opposite of that we observed here using healthy controls 9 . This is consistent with the present finding that eHHV-6B, but not eHHV-6A, is associated with an increased risk of several diseases and should motivate additional biobank-scale studies. It is also consistent with the independent replication data from the All of Us cohort.

Multiple facets of our analysis support a meaningful association between SLE and eHHV-6B. Patients with SLE with eHHV-6B had higher SLEDAI scores, suggesting a role for HHV-6B in SLE pathogenesis and clinical activity, and half of them had other autoimmune complications. We identified the presence of distinctive immune responses against HHV-6B in eHHV-6B carriers with SLE. These findings are consistent with previous studies, reporting that people with eHHV-6 have higher antibody levels against HHV-6-encoded protein IE1 (ref. 19 ) and show different adaptive immune responses to other herpesviruses and immune gene expression patterns 36 . It is possible that the antibodies against the IE-A transactivator can act as autoantibodies; however, we cannot exclude the possibility that these are a consequence of broader immune activation against exogenous HHV-6 aided or triggered by eHHV-6B. Herpesviruses often encode immunomodulatory viral chemokine receptors 37 , and HHV-6 encodes the chemokine receptors U12 and U51 (ref. 38 ). The expression of these viral proteins may also contribute to pathological immune activation by eHHV-6B. A previous study found an association between SLE and circulating HHV-6 infection 39 . Single-cell analysis of PBMCs showed increased antiviral responses in monocytes and DCs. scRNA-seq data showed no eHHV-6 expression in PBMCs; however, this is far from mutually exclusive with the activation of humoral immunity against eHHV-6B we observed. eHHV-6 genes are weakly expressed in multiple human tissues other than blood cells, such as the brain and esophagus in the genotype-tissue expression, which recruited individuals without respect to specific diseases 19 . It may be possible that eHHV-6B carriers of SLE acquired antibodies to the HHV-6B protein after the onset of SLE. Although we cannot identify the timing of B cell activation against IE-A transactivator from the present dataset, it may be possible that eHHV-6B is one factor contributing to severity, rather than the initiating step of SLE. Further investigation into the ‘endoimmunity’ hypothesis as well as T cell activation characteristic of eHHV-6B carriers will shed additional light on the immunopathogenicity of eHHV-6B in SLE.

The variants in high LD with eHHV-6B found by the comprehensive analysis using T2T-CHM13 were unlikely to be pathogenic for SLE. Based on the present associations alone, especially in a single population, we cannot rule out the possibility that a variant that is in high LD with eHHV-6B is the mutation causative of increased SLE risk. However, eHHV-6B found in the Japanese and European populations is inserted in different chromosomes 9 . This implies that the variants in high LD with eHHV-6B in the Japanese population are different from those in the European population. Considering that eHHV-6B is associated with SLE in both populations, it is unlikely that a variant in high LD with eHHV-6B has a causal role in SLE. Furthermore, eHHV-6B is a polymorphism that has been overlooked in previous GWASs because it is not in the reference genome, complementarily enhancing our knowledge of human disease risk variants.

Anelloviruses are characterized by enormous genetic diversity 40 – 42 , which must be considered when detecting and estimating their viral load from WGS data. Here we did so by comprehensively using more than 3,000 anellovirus reference sequences. However, as with any mapping-based approach, we are limited by the degree to which these reference genomes represent the viral diversity in this population. Unlike germline integration of HHV-6, anellovirus status in the patient is not expected to match that of a cell line derived from the patient. Thus, 5,597 Japanese blood WGS samples were analyzed, and the proportions of three continuous anellovirus detection states were assessed. Elevated anellovirus load was frequently found in the blood of participants with SLE, RA and COVID-19. Interestingly, patients with COVID-19 were anellovirus-positive in a similar proportion as healthy controls but included significantly higher anellovirus-infected and anellovirus-viremic than healthy controls. A possible scenario is that the atypical immune states characteristic of these diseases (or treatments) allow anelloviruses to reach higher levels in the blood 43 , 44 . This would also be supported by the observed association between anellovirus positivity and immunosuppressant and/or prednisolone use in SLE. An association between anellovirus positivity and low complement levels in SLE suggests that, similar to other viral infections, complement deficiency predisposes patients with SLE to anellovirus infection 45 , 46 . It is also possible that higher anellovirus levels contribute to the immune responses and inflammation leading to these diagnoses. Previous studies have shown that anellovirus DNA can induce proinflammatory cytokines 47 , 48 , suggesting that the higher load of anellovirus is a driver of severe inflammation during acute respiratory virus infection and, potentially, autoimmune response. While high-quality evidence of a causal relationship between anellovirus and human disease is limited, clonally expanded CD4 + T cells recognizing an anellovirus, torque teno virus, have been isolated from the cerebrospinal fluid of a patient with MS 49 . Notably, our conclusion is similar to that of ref. 50 that the prevalence of detectable (at any level) torque teno virus DNA is not higher in patients with SLE or RA than in healthy controls. We did observe, however, that participants with SLE and RA are more likely to have elevated levels of anellovirus DNA, supporting the hypothesis that anellovirus replication may be involved in autoimmune diseases. Longitudinal rather than cross-sectional studies are needed to adequately assess this. Notably, in our results, the increased number of anellovirus-infected participants among those with fatal COVID-19 is consistent with the hypothesis that anellovirus load may be higher in patients with more severe COVID-19.

Although only a small number of individuals carry eHHV-6 or a high load of anellovirus, positivity appears to have a striking impact on disease risk and clinical phenotype, especially compared with previously identified genetic and environmental factors of these diseases 20 – 22 . This suggests that the presence of these viruses in an individual’s virome is a potentially clinically useful biomarker for these diseases, enabling a personalized medicine approach to prevent disease. In addition to disease risk, eHHV-6 may be an important factor in clinical practice, such as cell therapies 51 . With the expanding availability of clinical WGS, incorporating bioinformatic pipelines that assay the virome may be a cost-effective way to screen for these individuals. The link between endogenous viruses or viral infection and immune-related diseases is a longstanding mystery in basic and clinical science. This study provides part of the missing pieces and outlines a path through which the study of the human virome can be integrated into the genotype–phenotype association analysis for which human WGS are often collected. Additional clinical and translational studies are needed to better assess the pathologically relevant mechanisms through which eHHV-6 and anellovirus modulate immunity.

Methods

Ethics

We have complied with all relevant ethical regulations. This study was approved by the ethical committees of Osaka University Graduate School of Medicine and affiliated institutes. All the participants provided written informed consent, as approved by the ethical committees of the affiliated institutes (Osaka University, Kyoto University and Keio University). No compensation was provided to the participants.

Study populations

We enrolled 2,401 Japanese individuals comprising 1,062 PV, 378 RA, 251 SLE, 216 PAP, 36 MS, 94 COVID-19 patients and 364 healthy controls at Osaka University and related institutes. We also used CRAM files from WGS of 2,126 healthy controls generated by the Nagahama Study 52 and related studies at Kyoto University and BAM files from WGS of 1,070 patients with COVID-19 generated by the Japan COVID-19 Task Force 53 . For all patients, blood was drawn when the disease was active. Especially for patients with COVID-19, blood samples were drawn immediately after diagnosis. Also, we used the genomic DNA of 295 patients with PV and 429 healthy controls obtained from EBV-transformed B-lymphoblast cell lines of unrelated Japanese individuals established by the Japan Biological Informatics Consortium (JBIC) 54 .

We conducted principal component analysis (PCA) against HapMap3 data 55 using SNP data of the same individuals to confirm the East Asian genetic background.

WGS data preprocessing

Genomic DNA was extracted from the whole blood of individuals at Osaka University and related institutes, with all samples processed at Osaka University. The DNA samples from Osaka University and those from the JBIC cell lines were sequenced at Macrogen Japan Corporation. Sequencing libraries were prepared using reagents from the TruSeq DNA PCR-free Library Prep kit (Illumina) following the manufacturer’s protocols. Also, 2× 150 bp paired-end WGS was conducted on HiSeqX or NovaSeq (Illumina). The BWA-MEM (v0.7.13 or v0.7.8) 56 was used to align the sequenced reads to the human reference build hs37d5 (ref. 57 ). The reads in the raw alignments were flagged for duplication using Picard MarkDuplicates (v2.10.10 or v2.17.2). The number of average aligned reads was 386 million per sample, corresponding to an average genomic coverage of 18×. WGS CRAM files from the Nagahama Study and related studies contained 370 million 2× 150 bp paired-end reads per sample, aligned to hs37d5 with BWA-MEM (v0.7.13), and followed by Picard MarkDuplicates (v2.10.10), with a mean coverage of 17×. WGS BAM files from the Japan COVID-19 Task Force contained 751 million 2× 150 bp paired-end reads per sample, aligned to hs38DH 58 with BWA-MEM (v0.7.15), and followed by GATK4 MarkDuplicates 59 (v4.1.0.0), with a mean coverage of 26×. In total, 6,321 BAM or CRAM files were prepared for downstream analyses.

Detection of eHHV-6A/eHHV-6B from WGS

The 6,321 WGS datasets were analyzed to detect eHHV-6 using a previously developed bioinformatics pipeline for eHHV-6 detection and reconstruction 11 . Briefly, we extracted reads that did not map to the human reference genome and mapped these reads to the reference HHV-6A (National Center for Biotechnology Information (NCBI) accession NC_001664.4 ) and HHV-6B ( NC_000898.1 ), followed by flagging for duplication to remove PCR duplicates. Then, the mapping depth and coverage of each virus were calculated. When 5% of the HHV-6 genome was covered at more than 2× read depth, we considered the sample to reflect the carriage of eHHV-6. Details are described in a previous paper 11 . The workflow described here was compiled as a Python pipeline and is available from the following GitHub repository: https://github.com/shohei-kojima/integrated_HHV6_recon (10.5281/zenodo.12528938) 60 . The ‘main.py’ in the repository wraps all the workflow described here.

Detection of ‘anellovirus reads’ from WGS

The 5,597 WGS datasets were analyzed to detect reads derived from anelloviruses, excluding 724 participants from JBIC cell lines because cell lines are unlikely to reflect anellovirus infection status in the participants. Due to the extensive sequence divergence among the anellovirus genomes 41 , 42 , we modified our previously reported pipeline to accurately detect and quantify anelloviruses. We downloaded human anellovirus FASTA files from NCBI Virus on 21 April 2022. At this time, 27,780 anellovirus sequences were deposited in NCBI. To find human anelloviruses with complete genomes, we referred to the metadata of the deposited sequences and extracted 3,333 complete human anellovirus sequences that were annotated as follows: ‘host’, ‘ Homo sapiens ’ and ‘Nuc_completeness’, ‘complete’. We found that some anellovirus genomes contain simple repeat sequences and regions with high similarity to human transposons. To avoid potential mismapping of human sequences to those regions, we masked simple repeats and human transposons by RepeatMasker 61 (v4.1.1) with ‘-species 9606 -s -no_is’ options. As a result, 1.09% of the anellovirus genomes were masked. The majority of the masked regions (103,757 bp) were simple repeats, but 69 and 1,724 bp were masked as long interspersed nuclear element-1 and long terminal repeat transposons, respectively, potentially reflecting mis-assemblies of viral genomes from metagenome sequencing data. We used the masked 3,333 genomes as anellovirus reference genomes. The human genome-unmapped reads were aligned to these 3,333 human anellovirus genomes using hisat2 (ref. 62 ) with ‘—mp 2,1 –no-spliced-alignment –all --secondary’ options, followed by flagging for duplication using Picard MarkDuplicates. Then, the read pairs mapped to at least one anellovirus sequence were counted as ‘Anellovirus reads’. Read pairs that were flagged as duplicates were excluded from ‘Anellovirus reads’. We considered a read pair to be an ‘Anellovirus read’ when at least one of the read pairs was mapping to at least one anellovirus genome. The workflow described here was compiled as a Python pipeline and is available from the following GitHub repository: https://github.com/shohei-kojima/human_anellovirus_detection (10.5281/zenodo.12528938) 60 . The ‘main.py’ in the repository wraps the workflow described here.

To evaluate the mismapping of human sequences, we mapped ‘Anellovirus reads’ in 2,715 participants (cohorts excluding COVID-19) to CHM13. We found no mismapping of human reads to anellovirus genomes in a majority of the samples ( n = 2,627). In 88 samples (3.2% of participants), we observed mismapping of human reads to anelloviruses, but the number of mismappings was at most one read. This indicates that mismapping of human sequences (potentially nonreference sequences) is extremely rare, but a small amount is still present.

Calculation of anellovirus load

We calculated anellovirus load (anellovirus copies per cell) to normalize anellovirus reads by WGS coverage as below: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{array}{l}{\rm{Anellovirus}}\,{\rm{load}}=\left(\frac{{\rm{anellovirus}}\,{\rm{reads}}}{{\rm{anellovirus}}\,{\rm{genome}}\,{\rm{length}}}\right)\\\qquad\qquad\qquad\qquad\left/\left(\frac{{\rm{human}}\,{\rm{reads}}}{{\rm{diploid}}\,{\rm{human}}\,{\rm{genome}}\,{\rm{length}}}\right)\right.\end{array}.$$\end{document} Anellovirus load = anellovirus reads anellovirus genome length human reads diploid human genome length . ‘Anellovirus genome length’ is the mean sequence length of 2,907.3 bp, calculated from 3,333 anellovirus FASTA files. ‘Human reads’ is the number of properly paired read pairs obtained using the SAMtools flagstat 63 . ‘Diploid human genome length’ is 5,618,659,254 bp in hs37d5 and 5,777,320,084 bp in hs38DH, calculated as the sum of 22 autosomal pairs and sex chromosome lengths of the human reference genome excluding Ns and gaps. The sex of the participants was not considered in this calculation.

Thresholds for anellovirus infection states: positive, infection and viremia

Considering the possibilities of human-derived reads mappable on anellovirus sequences and index hopping, we did not simply define anellovirus-positive in the situation when a single (or small number of) ‘anellovirus read’ was detected. Therefore, we robustly defined ‘anellovirus reads’ >5 (anellovirus load > 0.08) as anellovirus-positive (Extended Data Fig. 6 ). To be sure, we tentatively loosened the threshold a little and evaluated ‘anellovirus reads’ >3 (anellovirus load > 0.05) as anellovirus-positive, but the conclusion remained the same (Supplementary Table 7 ). When anellovirus load = 1, there is one dsDNA genome per cell on average; therefore, participants with anellovirus load >1 were defined as anellovirus-infected. The maximum value of anellovirus load in the 2,490 healthy controls was 7.4. Thus, we defined participants with anellovirus load >8 as those demonstrating an atypical state of anellovirus viremia.

Statistical analysis

The logistic regression model was used to calculate the adjusted case–control OR with 95% CI for the risk of each disease with eHHV-6 or each anellovirus infection state. We also used a two-sided Fisher’s exact test to evaluate the hypothesis of no association between the prevalence of eHHV-6A/eHHV-6B and the proportions of each anellovirus infection state in each disease and healthy controls. We applied the Bonferroni correction to adjust the P value for multiple test corrections. Then, the multivariable logistic regression model was performed with age and sex as independent variables. Anellovirus viremia was undetected in patients with PV, PAP or MS; therefore, no regression model was used. Wald test P value < 0.05 and Fisher’s test mid- P value < 0.05 were considered significant. All analyses were conducted using R (v4.0.5) with the R package MASS (v7.3.58) 64 and exact2x2 (v1.6.6) 65 . In each disease with a high OR, the linear or logistic regression model was used to estimate the relationship between the clinical metrics in each disease and eHHV-6 or each anellovirus infection state, with age and sex as independent variables. When the clinical metrics were binary, we also used Fisher’s test. Under regression analysis, listwise deletion of cases with missing values was applied.

Replication in All of Us

To find potential eHHV-6 carriers from biobank-scale WGS datasets, we developed a fast-screening method. This screening uses unmapped reads and detects HHV-6 reads by aligning with the reference genomes of HHV-6. It first extracts unmapped reads stored at the end of a BAM or CRAM file using the ‘samtools view’ command and converts reads to the fastq format. The number of unmapped reads varies between individuals and is influenced by mapping algorithms. To minimize the computational burden, our screening approach uses no more than 1 million unmapped reads by default. To avoid mismapping of telomere reads to the telomere-like repeats in the reference genomes of HHV-6, it filters out reads containing at least one unit of the telomere repeat (TAACCC). After filtering potential telomere reads, it aligns reads with the HHV-6A and HHV-6B reference genomes, NCBI accessions NC_001664.4 and NC_000898.1 , respectively, using hisat2 with the relaxed mismatch penalties ‘--mp 2,1 --no-spliced-alignment’. To reduce the computational burden, it treats paired-end reads as single-end reads when mapping to HHV-6. Then, it counts reads successfully mapped to HHV-6. To avoid counting non-HHV-6 reads partially aligned with HHV-6, it counts reads without clipping and insertion–deletion (InDels). Finally, it outputs the numbers of reads screened and those mapping to HHV-6. This method is available as a script ‘quick_check.py’ deposited in GitHub— https://github.com/shohei-kojima/integrated_HHV6_recon (10.5281/zenodo.12528938) 60 . Because WGS reads were mapped to the human reference genome by DRAGEN, most of the CRAM files contain more than 1 million unmapped reads (usually around 15 million unmapped reads). To reduce the computational burden, we used 1 million reads stored at the end of CRAM files for screening. In total, we screened 245,472 datasets in the All of Us Controlled Tier Dataset v7. For all datasets with six or more HHV-6 reads using the screening method, we applied the HHV-6 reconstruction method, ‘main.py’ in the integrated_HHV6_recon package to detect individuals carrying eHHV-6.

To define participants with SLE, we used diagnosis by the International Classification of Diseases, Ninth Revision (ICD9) (710.0) and ICD10 (M32) codes. To improve phenotype credibility, participants having two or more diagnoses on unique dates that were at least 1 year apart were defined as cases. As for PAP, another autoimmune disease with significant eHHV-6B risk in the Japanese population, we did not conduct a replication study given its ultra-rare disease prevalence (<10 cases per million population worldwide) 66 . We used all available participants for the association test, except for ones carrying ‘solo-DR’ and presumably two or more copies of full-length eHHV-6. We used the logistic regression model below. To avoid high type I errors, we performed a score test using saddle point approximation implemented in the R package SPAtest ( https://cran.r-project.org/web/packages/SPAtest/index.html ; v3.1.2). \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\rm{SLE}} \sim {\rm{eHHV{\hbox{-}}6}}+{\rm{age}}+{{\rm{age}}}^{2}+{\rm{sex}}+{\rm{top}}\,10\,{\rm{genetic}}\,{\rm{PCs}}.$$\end{document} SLE ~ eHHV- 6 + age + age 2 + sex + top 10 genetic PCs .

GWAS of eHHV-6B positivity

From each disease, we randomly selected ten times as many eHHV-6B-negative participants as eHHV-6B-positive cases (22 eHHV-6B-positive versus 220 eHHV-6B-negative). The FASTQ reads were aligned to T2T-CHM13v2.0 (ref. 67 ) with BWA-MEM (v0.7.27), followed by GATK4 MarkDuplicates and Base Quality Score Recalibration (v4.2.6.1) according to the GATK Best Practice. Then, we performed per-sample SNP and InDel calling using GATK4 HaplotypeCaller and joint genotyping using GATK4 GenomicsDBImport and GenotypeGVCF. We conducted LD-based genotype refinement for low-confidence genotypes and missing sites in WGS data using BEAGLE (v5.4) 68 with default settings. Variant quality control (QC) was accomplished by keeping autosomal SNPs and InDels with minor allele frequency (MAF) ≥ 0.05, while removing multi-allelic sites, spanning deletions and variants deviating from the Hardy–Weinberg equilibrium ( P < 1 × 10 −10 in cases and P < 1 × 10 −6 in controls). Individuals were excluded if they showed conflicting sex assignments between genetically inferred sex by variants and WGS coverage, deviating heterozygosity rate (±3 s.d.) or cryptic relatedness (pi-hat > 0.2). Finally, 6,464,509 variants and 238 participants (22 eHHV-6-B-positive versus 216 eHHV-6B-negative) passed QC.

For GWAS, logistic regression analysis was conducted using PLINK (v2.0) 69 to test the differences in the dominant genotype between eHHV-6B and controls. We used the dominant genotype because eHHV-6B was designated as binary without distinguishing the number of alleles. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\rm{eHHV{\hbox{-}}6B}} \sim {\rm{variant}}+{\rm{sex}}+{\rm{PC}}1+{\rm{PC}}2.$$\end{document} eHHV- 6 B ~ variant + sex + PC 1 + PC 2 . Focusing on the distal end of chr22q, we performed logistic regression including rare variants (MAF < 0.05). Variant QC other than MAF filtering, sample QC and PCs followed the GWAS described above.

PacBio HiFi long-read WGS

DNA samples from four eHHV-6B-positive patients with SLE were sheared targeting the size of 20 kb using Megaruptor 3 (Diagenode). HiFi sequencing libraries were prepared using SMARTbell Expres Template Prep Kit 2.0 according to the manufacturer’s protocols. The libraries were further size-selected electrophoretically using either the SageELF or Blue Pippin Systems (Sage Science). The libraries were sequenced on the Sequel II (Pacific Bioscience) system using the Sequel II Binding Kit 2.0 and Sequel II Sequencing Kit 2.0. Based on the sequenced subreads, PacBio HiFi reads were generated using the circular consensus sequence tool in SMART Link v9.0.0 (Pacific Bioscience). The HiFi reads were aligned to T2T-CHM13v2.0 using pbmm2 (v1.13.1; https://github.com/PacificBiosciences/pbmm2 ). DeepVariant (v1.6.0) 70 was used for SNP and InDel calling, and WhatsHap (v2.2) 71 was used for haplotype phasing. Phased genotypes were visually inspected using Integrative Genomics Viewer (v2.8.9) 72 .

HHV-6 peptide phage library design

Open reading frame (ORF) annotations for HHV-6A strain U1102 ( NC_001664.4 and KC465951.1 ) and HHV-6B strain Z29 ( NC_000898.1 , AF157706.1 and MW536483.1 ) were downloaded from NCBI. In total, 88 nonredundant ORFs for HHV-6A and 108 for HHV-6B were annotated. These ORFs were translated into amino acid sequences and then split into 64 amino acid peptide fragments with 50 amino acid overlaps. Fragments shorter than 64 amino acids at the C-terminus of a protein were extended with poly-alanine to reach the 64 amino acid length. This process yielded a total of 5,745 peptide fragments. For a positive control for the PhIP-seq library, we used a total of 248 human microbiota peptide antigens commonly detected in healthy individuals, as obtained from a previous study 73 . Using ICOR 74 , the peptide sequences were reverse-translated into DNA sequences optimized for Escherichia coli codon usage. HindIII and XhoI sites within the codon-optimized sequences were replaced with codons frequently used in E. coli (frequency >5% for each amino acid) in such a way that the protein sequence would be preserved.

DNA oligo library synthesis, cloning and phage library synthesis

DNA oligo pools encoding the peptide sequences of HHV-6 were synthesized by Twist Bioscience. Sequences of the DNA oligos are listed in Supplementary Data 1 . DNA oligos were PCR-amplified, digested with HindIII and XhoI and subsequently cloned into the HindIII/XhoI site of the T7 Select 10-3b phage vector. A ligation reaction was performed using T4 DNA Ligase. The resultant phage genomes were then packaged into T7 bacteriophage using the T7 Select Packaging Kit and amplified for use in PhIP-seq according to the manufacturer’s protocol.

PhIP-seq

We adapted previously established PhIP-seq protocols, with modifications 28 , 29 . IgG concentrations in serum or plasma samples were quantified using the Human IgG ELISA Kit. To minimize nonspecific binding, 1.5 ml tubes and Dynabeads protein A/G were blocked overnight at 4 °C on a rotator with 3% BSA in TBS-T. Sera or plasma containing 5 μg IgG and the HHV-6/human microbiota phage library diluted to 2 × 10 5 -fold representation (2 × 10 10 plaque-forming units for a library representing 10 5 peptide antigen repertoires) were added to preblocked tubes. Three technical replicates were performed for each sample, alongside 12 technical replicates for a negative control lacking sera or plasma. Phage-antibody mixtures were incubated overnight at 4 °C on a rotator, and 20 µl of Dynabeads protein A/G mix was added to each reaction. Following immunoprecipitation at 4 °C for 1 h on a rotator, beads were washed three times with 1,500 μl of PhIP-seq wash buffer (50 mM Tris–HCl (pH 7.5), 150 mM NaCl and 0.1% Tween 20) and resuspended in lysogeny broth with antibiotics. E. coli culture (BLT5403 strain) was then added to resuspended beads and incubated at 37 °C with shaking until lysis. Phage lysates made 0.5 M in NaCl were centrifuged, and the supernatants diluted to 2 × 10 5 -fold representation were used for the second round of immunoprecipitation. Immunoprecipitated beads were subsequently washed three times, resuspended and lysed at 95 °C for 10 min. Finally, the genomes of immunoprecipitated phages were amplified by PCR for Illumina sequencing. A total of 2,000 copies of spike-in DNA (enhanced green fluorescent protein sequence) for normalization were also included in each PCR reaction. Sequences of the spike-in DNA and primer sequences are listed in Supplementary Data 2 and 3 . Following the PCR, technical replicates from the same sample were pooled at equimolar concentration and size-selected using SPRIselect. DNA concentrations were determined by quantitative PCR using the KAPA Library Quantification Kit.

PhIP-seq data analysis

Phage genomes were amplified and sequenced on the Illumina NextSeq platform for 2× 150 bp paired-end sequencing, yielding at least 0.24 M spots per sample. Amplicons should contain unique oligo sequences flanking phage genome sequences. First, we counted reads containing the flanking phage sequences. We allowed one base mismatch between next-generation sequencing (NGS) reads and the flanking phage genome sequences. Across all samples, at least 84% of reads contained phage genome sequences flanked by unique oligo sequences, validating that most PCR amplicons are derived from phages. Next, we detected the unique oligo sequences from the reads carrying the flanking phage sequences. We allowed one base mismatch between NGS reads and target oligos. Reads 1 and 2 containing the same oligo sequence were assigned to that oligo. This process enabled the assignment of at least 74% of the total reads to the oligo library across all samples.

To check whether most of the phages displaying HHV-6 peptides are included in the phage library, we amplified and sequenced them by the Illumina NextSeq platform. The 4.4 M spots were obtained, and, on average, 713 reads were detected for each oligo. In total, 313 oligos were not abundantly detected; the read counts were less than 5% of the average (<36 reads). We excluded those oligos from the data analysis.

The immunoprecipitated phage counts were calculated based on the number of reads detecting the spike-in oligo. Across all samples, the number of reads detecting the spike-in oligo was enough (>14,000 reads) compared to the actual spike-in oligo copy number (that is, 2,000), ensuring the accuracy of estimation of immunoprecipitated phage copy numbers. The immunoprecipitated phage counts were then normalized by the abundance of each oligo in the phage library. For this normalization, we used the normalization factor given below: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm{facto{r}}_{i}=\frac{{\sum }_{n}C}{n\times{C}_{i}},$$\end{document} factor i = ∑ n C n × C i , where n is the number of oligos (that is, 5,432), C i is the read count of the i th oligo detected in the sequencing of the phage library and factor i is the normalization factor for the i th oligo.

The PhIP-seq score was calculated for each 14 amino acid region by averaging over the technical replicates and over the log-transformed normalized immunoprecipitated phage counts encoding the 14 amino acid region. The phage oligo detection from the NGS reads, normalization and the calculation of the PhIP-seq score were done by a custom Python script.

The statistical association between the presence of eHHV-6 and the PhIP-seq score was assessed by linear regression. We fit a model below, and the P for the β estimate of eHHV-6 was calculated by the Wald test (H0: β = 0). The model was fitted by the Python statsmodels module. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{array}{l}\mathrm{log{\hbox{-}}transformed}\,\mathrm{PhIP{\hbox{-}}seq}\,\mathrm{score} \sim \mathrm{eHHV{\hbox{-}}6}+\mathrm{age}\\+\mathrm{sex}+\mathrm{assay}\,\mathrm{date}+\mathrm{hospital}+\mathrm{immunosuppressant}\,\mathrm{use}\end{array}$$\end{document} log- t r a n s f o r m e d PhIP- s e q score ~ eHHV- 6 + age + sex + assay date + hospital + immunosuppressant use

scRNA-seq

Four eHHV-6B-positive patients with SLE and 11 eHHV-6B-negative patients with SLE who were close to them in age, sex and treatment history were selected. For both, blood was collected onto heparin tubes, and PBMCs were isolated using Leucosep (Greiner Bio-One) density gradient centrifugation according to the manufacturer’s instructions.

Single-cell suspensions were processed through the 10X Genomics Chromium Controller following the protocol outlined in the Chromium Single Cell V(D)J Reagent Kits (v1.1 Chemistry) User Guide. Chromium Next GEM Single Cell 5′ Library and Gel Bead Kit v1.1 (PN-1000167), Chromium Next GEM Chip G Single Cell Kit (PN-1000127) and Single Index Kit T Set A (PN-1000213) were applied during the process. Oil droplets of encapsulated single cells and barcoded beads (gel bead-in-emulsions (GEMs)) were subsequently reverse-transcribed in a Veriti Thermal Cycler (Thermo Fisher Scientific), resulting in cDNA tagged with a cell barcode and unique molecular index (UMI). Next, cDNA was then amplified to generate single-cell libraries according to the manufacturer’s protocol. Quantification was made with an Agilent Bioanalyzer High-Sensitivity DNA assay (High-Sensitivity DNA Kit; Agilent, 5067-4626). Subsequently amplified cDNA was enzymatically fragmented, end-repaired and polyA-tagged. Cleanup/size selection was performed on amplified cDNA using SPRIselect magnetic beads (Beckman-Coulter; SPRIselect, B23317 ). Next, Illumina sequencing adapters were ligated to the size-selected fragments and cleaned up using SPRIselect magnetic beads. Finally, sample indices were selected and amplified, followed by a double-sided size selection using SPRIselect magnetic beads. Final library quality was assessed using an Agilent Bioanalyzer High-Sensitivity DNA assay. Samples were then sequenced on an Illumina NovaSeq 6000 in paired-end mode.

scRNA-seq data processing

Raw sequencing data were first filtered and then mapped to a custom reference (GRCh38 + HHV-6A/HHV-6B) using Cell Ranger (v.7.0.1; 10x Genomics). Filtered expression matrices were used. Then, cell doublets were removed by scds (v1.10.0) 75 for each sample. QC was performed to remove low-quality cells with <200 genes expressed, >6,000 genes expressed, >60,000 UMIs, < 5% of reads from ribosomal genes or >10% of reads from hemoglobin genes. Thresholds for filtering out low-quality cells with a high proportion of mitochondrial reads were determined using miQC (v1.2.0) 76 .

The R package Seurat (v4.3.0) was used for data scaling, transformation, clustering, dimensionality reduction and visualization 32 . Data were scaled and transformed using the SCTransform() function. Then, the first 2,000 highly variable genes were selected for PCA on gene expression using the SelectIntegrationFeatures() function. The batch effect among different patients was removed by the Harmony (v1.0) 77 . Uniform Manifold Approximation and Projection (UMAP) dimension reduction was generated based on the first 30 harmony-adjusted principal components. The nearest neighbor graph was calculated using the FindNeighbors() function, followed by clustering using the FindClusters() function. Cell types were annotated using the Azimuth human PBMC reference 32 , 78 .

DE testing and GO enrichment

We used Milo (v1.2.0) 33 to test for the DE genes in cells within defined neighborhoods between the following two conditions: eHHV-6B-positive with SLE (cases) and eHHV-6B-negative with SLE (controls). We constructed a k -nearest neighbor graph and defined the representative neighborhoods on the k -nearest neighbor graph with default parameters. Then, we grouped overlapping neighborhoods that shared at least five cells and ignored their differences in abundance, although there was almost no difference. The neighborhood groups in which the number of cells in one case accounted for more than 60% of the total number of cells in the four cases were excluded from further analysis to avoid sampling bias. We performed DE analysis within these aggregated groups using the testDiffExp() function. We included age and sex in the model as covariates. Statistically significant DE genes were defined with FDR < 0.01 and |fold change | >2. To find the function of DE genes, GO biological process analysis with a P value < 0.01 and an FDR < 0.05 was performed using the enrichGO() and compareCluster() implemented in clusterProfiler (v4.7.1.3) 79 .

Antiviral ISGs scores

To evaluate IFN response to the virus, antiviral ISG scores were calculated using the AddModuleScore() function implemented in Seurat with default parameters. Thirty-three antiviral ISGs were selected based on previous reports 80 – 82 (Supplementary Table 8 ). The ISGs scores of cells in each neighborhood were used to compare case–control, and the differences between their respective mean of the ISG scores in cells from two conditions were calculated. If there were two or more case samples and two or more control samples with cells in the neighborhood, Welch’s t test was performed. To control for multiple testing, we adapted the spatial FDR implemented in Milo and used 5% spatial FDR as a threshold for significance. The spatial FDR and the differences in the ISG scores between two conditions in each neighborhood were used for visualization. Two-sided paired t test was performed to compare the mean of the ISG scores between two conditions in the neighborhoods within each neighborhood group.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Online content

Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at 10.1038/s41588-024-02022-z.

Supplementary information

Supplementary Information Contents of Supplementary Tables. Reporting Summary Supplementary Tables Supplementary Tables 1–8. Supplementary Data 1 Sequences of the oligo pool used for HHV-6 epitope mapping by PhIP-seq. Supplementary Data 2 Spike-in DNA sequence used in PhIP-seq. Supplementary Data 3 Primer sequences used for PhIP-seq library preparation.

Extended data

Extended Data Fig. 1 HHV-6 genome coverage and haplotypes on the distal end of chr22q in four eHHV-6B-positive patients with SLE. a , b , The mean of HHV-6A/B genome coverages normalized by half of the human genome coverage. The gray shade area represents the 95% confidence interval. Typically, eHHV-6-positive subjects carry the whole HHV-6 genome, but some carry only one of the terminal direct repeats (DRs)—solo-DR—and not the unique (U) region of the HHV-6 genome 9 . All four eHHV-6B-positive patients with SLE carried the whole HHV-6B genome, and two out of three eHHV-6B-positive patients with PAP carried only solo-DR. The mean normalized coverage of full-length HHV-6 was approximately 1, which supports the HHV-6 genomes we found are heritable endogenous viruses. Since the aligner could not distinguish between DR-L and DR-R sequences, solo-DR-derived reads were mapped to both. The mean normalized coverage of solo-DR was approximately 0.5, which also supports the endogeneity of the HHV-6 genome. c , Phased genotypes of the 16 variants with uncorrected P < 1 × 10 −4 in GWAS, derived from long-read sequencing of 4 eHHV-6B-positive patients with SLE. Dark red and black boxes indicate haploid visually confirmed with IGV. Light red and blue indicate genotypes from WGS joint genotyping. IGV, Integrative Genomics Viewer.

Extended Data Fig. 2 DA analysis, DE analysis and GO enrichment analysis of scRNA-seq data of PBMCs from eHHV-6B-positive and eHHV-6B-negative patients with SLE. a , Beeswarm plot showing the distribution of adjusted log 2 (fold change) in abundance between cases and controls in neighborhoods, with eHHV-6B-positive patients with SLE as case ( n = 4) and eHHV-6B-negative patients with SLE as control ( n = 11). b , The DE analysis between cases and controls in neighborhood groups as in a . DE genes (FDR < 0.01 and |fold change| > 2) are colored in red and labeled by gene symbols. c – e , Left: gene-concept network (cneplot) in neighborhood groups 6, 7 and 18. It shows the linkage of the DE genes and the concepts (GO biological process). The color of the genes represents the fold change values, and the size of the GO terms represents the number of associated genes. Right: directed acyclic graph for GO enrichment analysis in neighborhood groups 6, 7 and 18. Circles represent the GO function, colored by their FDR. Arrows indicate the relationship between the upstream and downstream GO functions.

Extended Data Fig. 3 Comparison of the antiviral ISG scores in patients with SLE on immunosuppressants. a , Graph representation of neighborhoods identified by Milo as in Fig. 5b . Node color indicates the difference between the mean of the antiviral ISGs scores of eHHV-6B-positive patients (case; n = 3) and eHHV-6B-negative patients (control; n = 6) with SLE on immunosuppressants (spatial FDR < 0.05). Neighborhoods with non-significant differences (spatial FDR ≥ 0.05) are colored white. Neighborhoods that contain only cells of either cases or controls are colored black. b , Comparison of differences within each neighborhood group shown in a using a two-sided paired t -test. Diff is the average difference between the case and control means. Boxes denote the IQR; the median is shown as horizontal bars; whiskers extend to 1.5 times the IQR. *FDR < 1.0 × 10 −10 , **FDR < 1.0 × 10 − 50 .

Extended Data Fig. 4 Reverse cumulative distributions of anellovirus load. The reverse cumulative distribution plot of anellovirus load in subjects with five autoimmune diseases and COVID-19 and healthy controls. Y-axis indicates the proportions of subjects with anellovirus load > values (x-axis). The proportions of anellovirus infection with anellovirus load >1 were significantly higher in patients with SLE, RA and COVID-19 than in healthy controls. The proportions of anellovirus viremia with anellovirus load >8.0 were significantly higher in patients with SLE, RA and COVID-19 than in healthy controls (dashed lines and asterisks).

Extended Data Fig. 5 Distribution of anellovirus load among various clinical metrics. Two histograms show the distribution of each axis. a , Strip plot of anellovirus load for patients with SLE without and with low serum complement levels. Anellovirus-positive subjects associated significantly with low serum complement levels (OR 4.4, Wald test uncorrected P = 0.020). NA, not available. b , Strip plot of anellovirus load for patients with SLE treated without and with immunosuppressants. Anellovirus-positive and anellovirus-infected subjects had a significantly positive correlation with immunosuppressant use (OR 3.3 and 7.3, Wald test uncorrected P = 0.0050 and 0.012, respectively). c , Strip plot of anellovirus load for patients with RA stratified by RF and ACPA. All anellovirus-viremic subjects were RF-positive and ACPA-positive. d , Strip plot of anellovirus load for patients with COVID-19 treated without and with oxygen therapy. All subjects with a high anellovirus load were treated with oxygen therapy. Horizontal dashed lines in a – d represent thresholds for anellovirus positivity (0.08), anellovirus infection (1.0) and anellovirus viremia (8.0).

Extended Data Fig. 6 The detection limit of the anellovirus load. Histogram of anellovirus load below 0.2. The colors of the bars indicate the number of ‘anellovirus reads’ in each sample. We set the detection limit of the anellovirus load to 0.08 (red dashed line). Subjects with ‘anellovirus reads’ >5 (anellovirus load >0.08) were defined as anellovirus-positive. Subjects with anellovirus load between 0 and 0.08 might be false-positive.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

A list of authors and their affiliations appears at the end of the paper.

These authors contributed equally: Noah Sasa, Shohei Kojima, Rie Koide.

Extended data

is available for this paper at 10.1038/s41588-024-02022-z.

Supplementary information

The online version contains supplementary material available at 10.1038/s41588-024-02022-z.

Acknowledgements

We thank all the participants involved in this study. We gratefully acknowledge All of Us participants for their contributions, without whom this research would not have been possible. We also thank the NIH’s All of Us Research Program for making available the participant data examined in this study. This research was supported by the Japan Society for the Promotion of Science KAKENHI (22H00476 and 24K02697 to Y.O., 20H05682 and 21H02972 to N.F.P., 22K15385 to S. Kojima and 24K18476 to R. Koide), Japan Agency for Medical Research and Development (AMED) (JP22ek0410075, JP23km0405211, JP23km0405217, JP23ek0109594, JP23ek0410113, JP223fa627002, JP223fa627010, JP233fa627011, JP23zf0127008, JP24wm0625504 and JP24gm1810011 to Y.O. and JP24gm7010004 to S. Kojima), JST Moonshot R&D (JPMJMS2021 and JPMJMS2024 to Y.O.), Takeda Science Foundation, Bioinformatics Initiative of Osaka University Graduate School of Medicine, Institute for Open and Transdisciplinary Research Initiatives, Center for Infectious Disease Education and Research, Center for Advanced Modality and Drug Delivery System (DDS), Osaka University, Cluster for Pioneering Research under the Hakubi fellowship program, the discretionary budget of the Director of the RIKEN Center for Integrative Medical Sciences and the RIKEN Special Postdoctoral Researcher Program. The funders had no role in study design, data collection and analysis, decision to publish or preparation of the paper.

Author contributions

N.F.P. and Y.O. designed and supervised the study. N.S., S. Kojima and Y.O. conducted the data analysis. S. Kojima, N.S. and N.F.P. developed methodology. R. Koide performed PhIP-seq. N.S., S. Kojima, R. Koide, A.G.B., N.F.P. and Y.O. wrote the paper. All the other co-authors collected and managed the samples and data construction.

Peer review

Peer review information

Nature Genetics thanks Mike Nalls and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

Data availability

Individual-level eHHV-6 and anellovirus data of the cases are publicly deposited at the Japanese Genotype–phenotype Archive (JGA) with accession codes JGAS000741 / JGAD000876 and JGAS000739 / JGAD000874 under controlled access. These individual-level eHHV-6 and anellovirus data can also be accessed through an application at the National Bioscience Database Center (NBDC) with the accession codes hum0197 and hum0343 . The GWAS summary statistics are publicly deposited at the NBDC with the accession code hum0197 . The WGS data of the controls are available from the Human Genetic Variation Database (HGVD) ( https://www.hgvd.genome.med.kyoto-u.ac.jp/ ). The All of Us data are available as a featured workspace to registered researchers of the All of Us Researcher Workbench (see details at https://www.researchallofus.org/ ). Short/long-read WGS of eHHV-6B-positive participants are publicly deposited at the European Genome–Phenome Archive (EGA) with the accession code EGAS00001007886 under controlled access. PhIP-seq data of eHHV-6B-positive and eHHV-6B-negative participants with SLE are publicly deposited at EGA with the accession code EGAS00001007955 under controlled access. scRNA-seq data of eHHV-6B-positive and eHHV-6B-negative participants with SLE are publicly deposited at the JGA with accession codes JGAS000741 / JGAD000876 under controlled access. This scRNA-seq data can also be accessed through the application at the NBDC with the accession code hum0197 ( https://humandbs.dbcls.jp/en/hum0197-latest ). The 3,333 human anellovirus complete genomes were downloaded as FASTA files from the NCBI nucleotide database ( https://www.ncbi.nlm.nih.gov ) on 21 April 2022.

Code availability

The codes to detect virome, integrated_HHV6_recon and human_anellovirus_detection, and to perform WGS data preprocessing, GWAS, PhIP-seq analysis and scRNA-seq analysis are deposited in Zenodo (10.5281/zenodo.12528938 (ref. 60 ) and 10.5281/zenodo.13953637 (ref. 83 )).

Competing interests

The authors declare no competing interests.

References

1. Hardy TA Blum S McCombe PA Reddel SW Guillain–Barré syndrome: modern theories of etiology Curr. Allergy Asthma Rep. 2011 11 197 204 21451970 10.1007/s11882-011-0190-y Hardy, T. A., Blum, S., McCombe, P. A. & Reddel, S. W. Guillain–Barré syndrome: modern theories of etiology. Curr. Allergy Asthma Rep. 11 , 197–204 (2011). 21451970 10.1007/s11882-011-0190-y

2. Bjornevik K Longitudinal analysis reveals high prevalence of Epstein–Barr virus associated with multiple sclerosis Science 2022 375 296 301 35025605 10.1126/science.abj8222 Bjornevik, K. et al. Longitudinal analysis reveals high prevalence of Epstein–Barr virus associated with multiple sclerosis. Science 375 , 296–301 (2022). 35025605 10.1126/science.abj8222

3. Lanz TV Clonally expanded B cells in multiple sclerosis bind EBV EBNA1 and GlialCAM Nature 2022 603 321 327 35073561 10.1038/s41586-022-04432-7 PMC9382663 Lanz, T. V. et al. Clonally expanded B cells in multiple sclerosis bind EBV EBNA1 and GlialCAM. Nature 603 , 321–327 (2022). 35073561 10.1038/s41586-022-04432-7 PMC9382663

4. Miller SD Persistent infection with Theiler’s virus leads to CNS autoimmunity via epitope spreading Nat. Med. 1997 3 1133 1136 9334726 10.1038/nm1097-1133 Miller, S. D. et al. Persistent infection with Theiler’s virus leads to CNS autoimmunity via epitope spreading. Nat. Med. 3 , 1133–1136 (1997). 9334726 10.1038/nm1097-1133

5. Tokuyama M Antibodies against human endogenous retrovirus K102 envelope activate neutrophils in systemic lupus erythematosus J. Exp. Med. 2021 218 e20191766 34019642 10.1084/jem.20191766 PMC8144942 Tokuyama, M. et al. Antibodies against human endogenous retrovirus K102 envelope activate neutrophils in systemic lupus erythematosus. J. Exp. Med. 218 , e20191766 (2021). 34019642 10.1084/jem.20191766 PMC8144942

6. Fajgenbaum DC June CH Cytokine storm N. Engl. J. Med. 2020 383 2255 2273 33264547 10.1056/NEJMra2026131 PMC7727315 Fajgenbaum, D. C. & June, C. H. Cytokine storm. N. Engl. J. Med. 383 , 2255–2273 (2020). 33264547 10.1056/NEJMra2026131 PMC7727315

7. Ikuta K Satoh Y Hoshikawa Y Sairenji T Detection of Epstein–Barr virus in salivas and throat washings in healthy adults and children Microbes Infect. 2000 2 115 120 10742683 10.1016/s1286-4579(00)00277-x Ikuta, K., Satoh, Y., Hoshikawa, Y. & Sairenji, T. Detection of Epstein–Barr virus in salivas and throat washings in healthy adults and children. Microbes Infect. 2 , 115–120 (2000). 10742683 10.1016/s1286-4579(00)00277-x

8. Moustafa A The blood DNA virome in 8,000 humans PLoS Pathog. 2017 13 e1006292 28328962 10.1371/journal.ppat.1006292 PMC5378407 Moustafa, A. et al. The blood DNA virome in 8,000 humans. PLoS Pathog. 13 , e1006292 (2017). 28328962 10.1371/journal.ppat.1006292 PMC5378407

9. Liu X Endogenization and excision of human herpesvirus 6 in human genomes PLoS Genet. 2020 16 e1008915 32776928 10.1371/journal.pgen.1008915 PMC7444522 Liu, X. et al. Endogenization and excision of human herpesvirus 6 in human genomes. PLoS Genet. 16 , e1008915 (2020). 32776928 10.1371/journal.pgen.1008915 PMC7444522

10. Wood ML Variation in human herpesvirus 6B telomeric integration, excision, and transmission between tissues and individuals eLife 2021 10 e70452 34545807 10.7554/eLife.70452 PMC8492063 Wood, M. L. et al. Variation in human herpesvirus 6B telomeric integration, excision, and transmission between tissues and individuals. eLife 10 , e70452 (2021). 34545807 10.7554/eLife.70452 PMC8492063

11. Kojima S Kamada AJ Parrish NF Virus-derived variation in diverse human genomes PLoS Genet. 2021 17 e1009324 33901175 10.1371/journal.pgen.1009324 PMC8101998 Kojima, S., Kamada, A. J. & Parrish, N. F. Virus-derived variation in diverse human genomes. PLoS Genet. 17 , e1009324 (2021). 33901175 10.1371/journal.pgen.1009324 PMC8101998

12. Aswad A Evolutionary history of endogenous human Herpesvirus 6 reflects human migration out of Africa Mol. Biol. Evol. 2021 38 96 107 32722766 10.1093/molbev/msaa190 PMC7782865 Aswad, A. et al. Evolutionary history of endogenous human Herpesvirus 6 reflects human migration out of Africa. Mol. Biol. Evol. 38 , 96–107 (2021). 32722766 10.1093/molbev/msaa190 PMC7782865

13. Kumata R Ito J Takahashi K Suzuki T Sato K A tissue level atlas of the healthy human virome BMC Biol. 2020 18 55 32493363 10.1186/s12915-020-00785-5 PMC7269688 Kumata, R., Ito, J., Takahashi, K., Suzuki, T. & Sato, K. A tissue level atlas of the healthy human virome. BMC Biol. 18 , 55 (2020). 32493363 10.1186/s12915-020-00785-5 PMC7269688

14. Abbas AA Redondoviridae, a family of small, circular DNA viruses of the human oro-respiratory tract associated with periodontitis and critical illness Cell Host Microbe 2019 25 719 729.e4 31071295 10.1016/j.chom.2019.04.001 PMC6510254 Abbas, A. A. et al. Redondoviridae, a family of small, circular DNA viruses of the human oro-respiratory tract associated with periodontitis and critical illness. Cell Host Microbe 25 , 719–729.e4 (2019). 31071295 10.1016/j.chom.2019.04.001 PMC6510254

15. Abbas AA The perioperative lung transplant virome: torque teno viruses are elevated in donor lungs and show divergent dynamics in primary graft dysfunction Am. J. Transpl. 2017 17 1313 1324 10.1111/ajt.14076 PMC5389935 27731934 Abbas, A. A. et al. The perioperative lung transplant virome: torque teno viruses are elevated in donor lungs and show divergent dynamics in primary graft dysfunction. Am. J. Transpl. 17 , 1313–1324 (2017). 10.1111/ajt.14076 PMC5389935 27731934

16. Maggi F TT virus loads and lymphocyte subpopulations in children with acute respiratory diseases J. Virol. 2003 77 9081 9083 12885924 10.1128/JVI.77.16.9081-9083.2003 PMC167221 Maggi, F. et al. TT virus loads and lymphocyte subpopulations in children with acute respiratory diseases. J. Virol. 77 , 9081–9083 (2003). 12885924 10.1128/JVI.77.16.9081-9083.2003 PMC167221

17. Gravel A Inherited chromosomally integrated human herpesvirus 6 as a predisposing risk factor for the development of angina pectoris Proc. Natl Acad. Sci. USA 2015 112 8058 8063 26080419 10.1073/pnas.1502741112 PMC4491735 Gravel, A. et al. Inherited chromosomally integrated human herpesvirus 6 as a predisposing risk factor for the development of angina pectoris. Proc. Natl Acad. Sci. USA 112 , 8058–8063 (2015). 26080419 10.1073/pnas.1502741112 PMC4491735

18. Gaccioli F Fetal inheritance of chromosomally integrated human herpesvirus 6 predisposes the mother to pre-eclampsia Nat. Microbiol. 2020 5 901 908 32367053 10.1038/s41564-020-0711-3 PMC7610361 Gaccioli, F. et al. Fetal inheritance of chromosomally integrated human herpesvirus 6 predisposes the mother to pre-eclampsia. Nat. Microbiol. 5 , 901–908 (2020). 32367053 10.1038/s41564-020-0711-3 PMC7610361

19. Peddu V Inherited chromosomally integrated human herpesvirus 6 demonstrates tissue-specific RNA expression in vivo that correlates with an increased antibody immune response J. Virol. 2019 94 e01418-19 31597766 10.1128/JVI.01418-19 PMC6912112 Peddu, V. et al. Inherited chromosomally integrated human herpesvirus 6 demonstrates tissue-specific RNA expression in vivo that correlates with an increased antibody immune response. J. Virol. 94 , e01418-19 (2019). 31597766 10.1128/JVI.01418-19 PMC6912112

20. Sakaue S Genetic determinants of risk in autoimmune pulmonary alveolar proteinosis Nat. Commun. 2021 12 1032 33589587 10.1038/s41467-021-21011-y PMC7884840 Sakaue, S. et al. Genetic determinants of risk in autoimmune pulmonary alveolar proteinosis. Nat. Commun. 12 , 1032 (2021). 33589587 10.1038/s41467-021-21011-y PMC7884840

21. Okada Y A genome-wide association study identified AFF1 as a susceptibility locus for systemic lupus eyrthematosus in Japanese PLoS Genet. 2012 8 e1002455 22291604 10.1371/journal.pgen.1002455 PMC3266877 Okada, Y. et al. A genome-wide association study identified AFF1 as a susceptibility locus for systemic lupus eyrthematosus in Japanese. PLoS Genet. 8 , e1002455 (2012). 22291604 10.1371/journal.pgen.1002455 PMC3266877

22. Ogawa K Next-generation sequencing identifies contribution of both class I and II HLA genes on susceptibility of multiple sclerosis in Japanese J. Neuroinflammation 2019 16 162 31382992 10.1186/s12974-019-1551-z PMC6683481 Ogawa, K. et al. Next-generation sequencing identifies contribution of both class I and II HLA genes on susceptibility of multiple sclerosis in Japanese. J. Neuroinflammation 16 , 162 (2019). 31382992 10.1186/s12974-019-1551-z PMC6683481

23. Bombardier C Gladman DD Urowitz MB Caron D Chang CH Derivation of the SLEDAI. A disease activity index for lupus patients. The Committee on Prognosis Studies in SLE Arthritis Rheum. 1992 35 630 640 1599520 10.1002/art.1780350606 Bombardier, C., Gladman, D. D., Urowitz, M. B., Caron, D. & Chang, C. H. Derivation of the SLEDAI. A disease activity index for lupus patients. The Committee on Prognosis Studies in SLE. Arthritis Rheum. 35 , 630–640 (1992). 1599520 10.1002/art.1780350606

24. Arbuckle MR Development of autoantibodies before the clinical onset of systemic lupus erythematosus N. Engl. J. Med. 2003 349 1526 1533 14561795 10.1056/NEJMoa021933 Arbuckle, M. R. et al. Development of autoantibodies before the clinical onset of systemic lupus erythematosus. N. Engl. J. Med. 349 , 1526–1533 (2003). 14561795 10.1056/NEJMoa021933

25. Latzin P Anti-GM-CSF antibodies in paediatric pulmonary alveolar proteinosis Thorax 2005 60 39 44 15618581 10.1136/thx.2004.021329 PMC1747161 Latzin, P. et al. Anti-GM-CSF antibodies in paediatric pulmonary alveolar proteinosis. Thorax 60 , 39–44 (2005). 15618581 10.1136/thx.2004.021329 PMC1747161

26. All of Us Research Program Investigators et al. The ‘All of Us’ research program. N. Engl. J. Med . 381 ,668–676 (2019). 10.1056/NEJMsr1809937 PMC8291101 31412182

27. Huang Y Human telomeres that carry an integrated copy of human herpesvirus 6 are often short and unstable, facilitating release of the viral genome from the chromosome Nucleic Acids Res. 2014 42 315 327 24057213 10.1093/nar/gkt840 PMC3874159 Huang, Y. et al. Human telomeres that carry an integrated copy of human herpesvirus 6 are often short and unstable, facilitating release of the viral genome from the chromosome. Nucleic Acids Res. 42 , 315–327 (2014). 24057213 10.1093/nar/gkt840 PMC3874159

28. Mohan D PhIP-seq characterization of serum antibodies using oligonucleotide-encoded peptidomes Nat. Protoc. 2018 13 1958 1978 30190553 10.1038/s41596-018-0025-6 PMC6568263 Mohan, D. et al. PhIP-seq characterization of serum antibodies using oligonucleotide-encoded peptidomes. Nat. Protoc. 13 , 1958–1978 (2018). 30190553 10.1038/s41596-018-0025-6 PMC6568263

29. Shrock EL Shrock CL Elledge SJ VirScan: high-throughput profiling of antiviral antibody epitopes Bio Protoc. 2022 12 e4464 35937932 10.21769/BioProtoc.4464 PMC9303818 Shrock, E. L., Shrock, C. L. & Elledge, S. J. VirScan: high-throughput profiling of antiviral antibody epitopes. Bio Protoc. 12 , e4464 (2022). 35937932 10.21769/BioProtoc.4464 PMC9303818

30. De Bolle L Naesens L De Clercq E Update on human herpesvirus 6 biology, clinical features, and therapy Clin. Microbiol. Rev. 2005 18 217 245 15653828 10.1128/CMR.18.1.217-245.2005 PMC544175 De Bolle, L., Naesens, L. & De Clercq, E. Update on human herpesvirus 6 biology, clinical features, and therapy. Clin. Microbiol. Rev. 18 , 217–245 (2005). 15653828 10.1128/CMR.18.1.217-245.2005 PMC544175

31. Takemoto K Time-series analysis comparing the prevalence of antibodies against nine viral species found in umbilical cord blood in Japan Jpn. J. Infect. Dis. 2016 69 314 318 26567842 10.7883/yoken.JJID.2015.366 Takemoto, K. et al. Time-series analysis comparing the prevalence of antibodies against nine viral species found in umbilical cord blood in Japan. Jpn. J. Infect. Dis. 69 , 314–318 (2016). 26567842 10.7883/yoken.JJID.2015.366

32. Hao Y Integrated analysis of multimodal single-cell data Cell 2021 184 3573 3587 34062119 10.1016/j.cell.2021.04.048 PMC8238499 Hao, Y. et al. Integrated analysis of multimodal single-cell data. Cell 184 , 3573–3587 (2021). 34062119 10.1016/j.cell.2021.04.048 PMC8238499

33. Dann E Henderson NC Teichmann SA Morgan MD Marioni JC Differential abundance testing on single-cell data using k-nearest neighbor graphs Nat. Biotechnol. 2022 40 245 253 34594043 10.1038/s41587-021-01033-z PMC7617075 Dann, E., Henderson, N. C., Teichmann, S. A., Morgan, M. D. & Marioni, J. C. Differential abundance testing on single-cell data using k-nearest neighbor graphs. Nat. Biotechnol. 40 , 245–253 (2022). 34594043 10.1038/s41587-021-01033-z PMC7617075

34. Tanaka-Taya K Human herpesvirus 6 (HHV-6) is transmitted from parent to child in an integrated form and characterization of cases with chromosomally integrated HHV-6 DNA J. Med. Virol. 2004 73 465 473 15170644 10.1002/jmv.20113 Tanaka-Taya, K. et al. Human herpesvirus 6 (HHV-6) is transmitted from parent to child in an integrated form and characterization of cases with chromosomally integrated HHV-6 DNA. J. Med. Virol. 73 , 465–473 (2004). 15170644 10.1002/jmv.20113

35. Miura H Chromosomally integrated human herpesvirus 6 in the Japanese population J. Med. Virol. 2018 90 1636 1642 29905966 10.1002/jmv.25244 Miura, H. et al. Chromosomally integrated human herpesvirus 6 in the Japanese population. J. Med. Virol. 90 , 1636–1642 (2018). 29905966 10.1002/jmv.25244

36. Kumata R Ito J Sato K Inherited chromosomally integrated HHV-6 possibly modulates human gene expression Virus Genes 2020 56 386 389 32067146 10.1007/s11262-020-01745-5 Kumata, R., Ito, J. & Sato, K. Inherited chromosomally integrated HHV-6 possibly modulates human gene expression. Virus Genes 56 , 386–389 (2020). 32067146 10.1007/s11262-020-01745-5

37. De Munnik SM Smit MJ Leurs R Vischer HF Modulation of cellular signaling by herpesvirus-encoded G protein-coupled receptors Front. Pharmacol. 2015 6 40 25805993 10.3389/fphar.2015.00040 PMC4353375 De Munnik, S. M., Smit, M. J., Leurs, R. & Vischer, H. F. Modulation of cellular signaling by herpesvirus-encoded G protein-coupled receptors. Front. Pharmacol. 6 , 40 (2015). 25805993 10.3389/fphar.2015.00040 PMC4353375

38. Isegawa Y Ping Z Nakano K Sugimoto N Yamanishi K Human herpesvirus 6 open reading frame U12 encodes a functional β-chemokine receptor J. Virol. 1998 72 6104 6112 9621074 10.1128/jvi.72.7.6104-6112.1998 PMC110416 Isegawa, Y., Ping, Z., Nakano, K., Sugimoto, N. & Yamanishi, K. Human herpesvirus 6 open reading frame U12 encodes a functional β-chemokine receptor. J. Virol. 72 , 6104–6112 (1998). 9621074 10.1128/jvi.72.7.6104-6112.1998 PMC110416

39. Broccolo F Selective reactivation of human herpesvirus 6 in patients with autoimmune connective tissue diseases J. Med. Virol. 2013 85 1925 1934 23983182 10.1002/jmv.23670 Broccolo, F. et al. Selective reactivation of human herpesvirus 6 in patients with autoimmune connective tissue diseases. J. Med. Virol. 85 , 1925–1934 (2013). 23983182 10.1002/jmv.23670

40. Worobey M Extensive homologous recombination among widely divergent TT viruses J. Virol. 2000 74 7666 7670 10906223 10.1128/jvi.74.16.7666-7670.2000 PMC112290 Worobey, M. Extensive homologous recombination among widely divergent TT viruses. J. Virol. 74 , 7666–7670 (2000). 10906223 10.1128/jvi.74.16.7666-7670.2000 PMC112290

41. Cebriá-Mendoza M Deep viral blood metagenomics reveals extensive anellovirus diversity in healthy humans Sci. Rep. 2021 11 6921 33767340 10.1038/s41598-021-86427-4 PMC7994813 Cebriá-Mendoza, M. et al. Deep viral blood metagenomics reveals extensive anellovirus diversity in healthy humans. Sci. Rep. 11 , 6921 (2021). 33767340 10.1038/s41598-021-86427-4 PMC7994813

42. Arze CA Global genome analysis reveals a vast and dynamic anellovirus landscape within the human virome Cell Host Microbe 2021 29 1305 1315.e6 34320399 10.1016/j.chom.2021.07.001 Arze, C. A. et al. Global genome analysis reveals a vast and dynamic anellovirus landscape within the human virome. Cell Host Microbe 29 , 1305–1315.e6 (2021). 34320399 10.1016/j.chom.2021.07.001

43. De Vlaminck I Temporal response of the human virome to immunosuppression and antiviral therapy Cell 2013 155 1178 1187 24267896 10.1016/j.cell.2013.10.034 PMC4098717 De Vlaminck, I. et al. Temporal response of the human virome to immunosuppression and antiviral therapy. Cell 155 , 1178–1187 (2013). 24267896 10.1016/j.cell.2013.10.034 PMC4098717

44. Shibayama T Inverse relationship between the titre of TT virus DNA and the CD4 cell count in patients infected with HIV AIDS 2001 15 563 570 11316992 10.1097/00002030-200103300-00004 Shibayama, T. et al. Inverse relationship between the titre of TT virus DNA and the CD4 cell count in patients infected with HIV. AIDS 15 , 563–570 (2001). 11316992 10.1097/00002030-200103300-00004

45. Stoermer KA Morrison TE Complement and viral pathogenesis Virology 2011 411 362 373 21292294 10.1016/j.virol.2010.12.045 PMC3073741 Stoermer, K. A. & Morrison, T. E. Complement and viral pathogenesis. Virology 411 , 362–373 (2011). 21292294 10.1016/j.virol.2010.12.045 PMC3073741

46. Jung J-Y Suh C-H Infection in systemic lupus erythematosus, similarities, and differences with lupus flare Korean J. Intern. Med. 2017 32 429 438 28490724 10.3904/kjim.2016.234 PMC5432804 Jung, J.-Y. & Suh, C.-H. Infection in systemic lupus erythematosus, similarities, and differences with lupus flare. Korean J. Intern. Med. 32 , 429–438 (2017). 28490724 10.3904/kjim.2016.234 PMC5432804

47. Zheng H Torque teno virus (SANBAN isolate) ORF2 protein suppresses NF-κB pathways via interaction with IκB kinases J. Virol. 2007 81 11917 11924 17686849 10.1128/JVI.01101-07 PMC2168763 Zheng, H. et al. Torque teno virus (SANBAN isolate) ORF2 protein suppresses NF-κB pathways via interaction with IκB kinases. J. Virol. 81 , 11917–11924 (2007). 17686849 10.1128/JVI.01101-07 PMC2168763

48. Rocchi J Torquetenovirus DNA drives proinflammatory cytokines production and secretion by immune cells via toll-like receptor 9 Virology 2009 394 235 242 19765789 10.1016/j.virol.2009.08.036 Rocchi, J. et al. Torquetenovirus DNA drives proinflammatory cytokines production and secretion by immune cells via toll-like receptor 9. Virology 394 , 235–242 (2009). 19765789 10.1016/j.virol.2009.08.036

49. Sospedra M Recognition of conserved amino acid motifs of common viruses and its role in autoimmunity PLoS Pathog. 2005 1 e41 16362076 10.1371/journal.ppat.0010041 PMC1315278 Sospedra, M. et al. Recognition of conserved amino acid motifs of common viruses and its role in autoimmunity. PLoS Pathog. 1 , e41 (2005). 16362076 10.1371/journal.ppat.0010041 PMC1315278

50. Seemayer CA Prevalence of TTV DNA and GBV-C RNA in patients with systemic sclerosis, rheumatoid arthritis, and osteoarthritis does not differ from that in healthy blood donors Ann. Rheum. Dis. 2001 60 806 809 11454648 10.1136/ard.60.8.806 PMC1753793 Seemayer, C. A. et al. Prevalence of TTV DNA and GBV-C RNA in patients with systemic sclerosis, rheumatoid arthritis, and osteoarthritis does not differ from that in healthy blood donors. Ann. Rheum. Dis. 60 , 806–809 (2001). 11454648 10.1136/ard.60.8.806 PMC1753793

51. Lareau CA Latent human herpesvirus 6 is reactivated in CAR T cells Nature 2023 623 608 615 37938768 10.1038/s41586-023-06704-2 PMC10999258 Lareau, C. A. et al. Latent human herpesvirus 6 is reactivated in CAR T cells. Nature 623 , 608–615 (2023). 37938768 10.1038/s41586-023-06704-2 PMC10999258

52. Yano, M., Matsuda, F., Sakuntabhai, A., & Hirota, S. (eds). Socio-Life Science and the COVID-19 Outbreak. Economics, Law, and Institutions in Asia Pacific (Springer, 2022).

53. Tanaka H Clinical characteristics of patients with coronavirus disease (COVID-19): preliminary baseline report of Japan COVID-19 task force, a nationwide consortium to investigate host genetics of COVID-19 Int. J. Infect. Dis. 2021 113 74 81 34601141 10.1016/j.ijid.2021.09.070 PMC8482546 Tanaka, H. et al. Clinical characteristics of patients with coronavirus disease (COVID-19): preliminary baseline report of Japan COVID-19 task force, a nationwide consortium to investigate host genetics of COVID-19. Int. J. Infect. Dis. 113 , 74–81 (2021). 34601141 10.1016/j.ijid.2021.09.070 PMC8482546

54. Hirata J Variants at HLA-A, HLA-C, and HLA-DQB1 confer risk of psoriasis vulgaris in Japanese J. Invest. Dermatol. 2018 138 542 548 29031612 10.1016/j.jid.2017.10.001 Hirata, J. et al. Variants at HLA-A, HLA-C, and HLA-DQB1 confer risk of psoriasis vulgaris in Japanese. J. Invest. Dermatol. 138 , 542–548 (2018). 29031612 10.1016/j.jid.2017.10.001

55. Altshuler DM Integrating common and rare genetic variation in diverse human populations Nature 2010 467 52 58 20811451 10.1038/nature09298 PMC3173859 Altshuler, D. M. et al. Integrating common and rare genetic variation in diverse human populations. Nature 467 , 52–58 (2010). 20811451 10.1038/nature09298 PMC3173859

56. Li, H. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. Preprint at 10.48550/arXiv.1303.3997 (2013).

57. 1000 Genomes Project Consortium et al. A global reference for human genetic variation. Nature 526 ,68–74 (2015). 10.1038/nature15393 PMC4750478 26432245

58. Church DM Extending reference assembly models Genome Biol. 2015 16 13 25651527 10.1186/s13059-015-0587-3 PMC4305238 Church, D. M. et al. Extending reference assembly models. Genome Biol. 16 , 13 (2015). 25651527 10.1186/s13059-015-0587-3 PMC4305238

59. McKenna A The genome analysis toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data Genome Res. 2010 20 1297 1303 20644199 10.1101/gr.107524.110 PMC2928508 McKenna, A. et al. The genome analysis toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 20 , 1297–1303 (2010). 20644199 10.1101/gr.107524.110 PMC2928508

60. Kojima, S. Software of ‘Blood DNA virome associates with autoimmune diseases and COVID-19’. Zenodo 10.5281/zenodo.12528938 (2024). 10.1038/s41588-024-02022-z PMC11735405 39753770

61. Smit, A., Hubley, R. & Green, P. RepeatMasker Open-4.0. www.repeatmasker.org (2013).

62. Kim D Paggi JM Park C Bennett C Salzberg SL Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype Nat. Biotechnol. 2019 37 907 915 31375807 10.1038/s41587-019-0201-4 PMC7605509 Kim, D., Paggi, J. M., Park, C., Bennett, C. & Salzberg, S. L. Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nat. Biotechnol. 37 , 907–915 (2019). 31375807 10.1038/s41587-019-0201-4 PMC7605509

63. Danecek P Twelve years of SAMtools and BCFtools Gigascience 2021 10 giab008 33590861 10.1093/gigascience/giab008 PMC7931819 Danecek, P. et al. Twelve years of SAMtools and BCFtools. Gigascience 10 , giab008 (2021). 33590861 10.1093/gigascience/giab008 PMC7931819

64. Venables, W. N. & Ripley, B. D. Modern Applied Statistics with S. Fourth Edition Vol. 53 (Springer, 2002).

65. Fay MP Confidence intervals that match Fisher’s exact or Blaker’s exact tests Biostatistics 2010 11 373 374 19948745 10.1093/biostatistics/kxp050 PMC2852239 Fay, M. P. Confidence intervals that match Fisher’s exact or Blaker’s exact tests. Biostatistics 11 , 373–374 (2010). 19948745 10.1093/biostatistics/kxp050 PMC2852239

66. Inoue Y Characteristics of a large cohort of patients with autoimmune pulmonary alveolar proteinosis in Japan Am. J. Respir. Crit. Care Med. 2008 177 752 762 18202348 10.1164/rccm.200708-1271OC PMC2720118 Inoue, Y. et al. Characteristics of a large cohort of patients with autoimmune pulmonary alveolar proteinosis in Japan. Am. J. Respir. Crit. Care Med. 177 , 752–762 (2008). 18202348 10.1164/rccm.200708-1271OC PMC2720118

67. Rhie A The complete sequence of a human Y chromosome Nature 2023 621 344 354 37612512 10.1038/s41586-023-06457-y PMC10752217 Rhie, A. et al. The complete sequence of a human Y chromosome. Nature 621 , 344–354 (2023). 37612512 10.1038/s41586-023-06457-y PMC10752217

68. Browning BL Zhou Y Browning SR A one-penny imputed genome from next-generation reference panels Am. J. Hum. Genet. 2018 103 338 348 30100085 10.1016/j.ajhg.2018.07.015 PMC6128308 Browning, B. L., Zhou, Y. & Browning, S. R. A one-penny imputed genome from next-generation reference panels. Am. J. Hum. Genet. 103 , 338–348 (2018). 30100085 10.1016/j.ajhg.2018.07.015 PMC6128308

69. Hill A Stepwise distributed open innovation contests for software development: acceleration of genome-wide association analysis GigaScience 2017 6 1 10 28327993 10.1093/gigascience/gix009 PMC5467032 Hill, A. et al. Stepwise distributed open innovation contests for software development: acceleration of genome-wide association analysis. GigaScience 6 , 1–10 (2017). 28327993 10.1093/gigascience/gix009 PMC5467032

70. Poplin R A universal SNP and small-indel variant caller using deep neural networks Nat. Biotechnol. 2018 36 983 987 30247488 10.1038/nbt.4235 Poplin, R. et al. A universal SNP and small-indel variant caller using deep neural networks. Nat. Biotechnol. 36 , 983–987 (2018). 30247488 10.1038/nbt.4235

71. Patterson M WhatsHap: weighted haplotype assembly for future-generation sequencing reads J. Comput. Biol. 2015 22 498 509 25658651 10.1089/cmb.2014.0157 Patterson, M. et al. WhatsHap: weighted haplotype assembly for future-generation sequencing reads. J. Comput. Biol. 22 , 498–509 (2015). 25658651 10.1089/cmb.2014.0157

72. Robinson JT Integrative genomics viewer Nat. Biotechnol. 2011 29 24 26 21221095 10.1038/nbt.1754 PMC3346182 Robinson, J. T. et al. Integrative genomics viewer. Nat. Biotechnol. 29 , 24–26 (2011). 21221095 10.1038/nbt.1754 PMC3346182

73. Vogl T Population-wide diversity and stability of serum antibody epitope repertoires against human microbiota Nat. Med. 2021 27 1442 1450 34282338 10.1038/s41591-021-01409-3 Vogl, T. et al. Population-wide diversity and stability of serum antibody epitope repertoires against human microbiota. Nat. Med. 27 , 1442–1450 (2021). 34282338 10.1038/s41591-021-01409-3

74. Jain R Jain A Mauro E LeShane K Densmore D ICOR: improving codon optimization with recurrent neural networks BMC Bioinformatics 2023 24 132 37016283 10.1186/s12859-023-05246-8 PMC10074884 Jain, R., Jain, A., Mauro, E., LeShane, K. & Densmore, D. ICOR: improving codon optimization with recurrent neural networks. BMC Bioinformatics 24 , 132 (2023). 37016283 10.1186/s12859-023-05246-8 PMC10074884

75. Bais AS Kostka D scds: computational annotation of doublets in single-cell RNA sequencing data Bioinformatics 2020 36 1150 1158 31501871 10.1093/bioinformatics/btz698 PMC7703774 Bais, A. S. & Kostka, D. scds: computational annotation of doublets in single-cell RNA sequencing data. Bioinformatics 36 , 1150–1158 (2020). 31501871 10.1093/bioinformatics/btz698 PMC7703774

76. Hippen AA miQC: an adaptive probabilistic framework for quality control of single-cell RNA-sequencing data PLoS Comput. Biol. 2021 17 e1009290 34428202 10.1371/journal.pcbi.1009290 PMC8415599 Hippen, A. A. et al. miQC: an adaptive probabilistic framework for quality control of single-cell RNA-sequencing data. PLoS Comput. Biol. 17 , e1009290 (2021). 34428202 10.1371/journal.pcbi.1009290 PMC8415599

77. Korsunsky I Fast, sensitive and accurate integration of single-cell data with Harmony Nat. Methods 2019 16 1289 1296 31740819 10.1038/s41592-019-0619-0 PMC6884693 Korsunsky, I. et al. Fast, sensitive and accurate integration of single-cell data with Harmony. Nat. Methods 16 , 1289–1296 (2019). 31740819 10.1038/s41592-019-0619-0 PMC6884693

78. Stuart T Comprehensive integration of single-cell data Cell 2019 177 1888 1902 31178118 10.1016/j.cell.2019.05.031 PMC6687398 Stuart, T. et al. Comprehensive integration of single-cell data. Cell 177 , 1888–1902 (2019). 31178118 10.1016/j.cell.2019.05.031 PMC6687398

79. Yu G Wang L-G Han Y He Q-Y clusterProfiler: an R package for comparing biological themes among gene clusters OMICS 2012 16 284 287 22455463 10.1089/omi.2011.0118 PMC3339379 Yu, G., Wang, L.-G., Han, Y. & He, Q.-Y. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 16 , 284–287 (2012). 22455463 10.1089/omi.2011.0118 PMC3339379

80. Schoggins JW Rice CM Interferon-stimulated genes and their antiviral effector functions Curr. Opin. Virol. 2011 1 519 525 22328912 10.1016/j.coviro.2011.10.008 PMC3274382 Schoggins, J. W. & Rice, C. M. Interferon-stimulated genes and their antiviral effector functions. Curr. Opin. Virol. 1 , 519–525 (2011). 22328912 10.1016/j.coviro.2011.10.008 PMC3274382

81. Schoggins JW Interferon-stimulated genes: what do they all do? Annu. Rev. Virol. 2019 6 567 584 31283436 10.1146/annurev-virology-092818-015756 Schoggins, J. W. Interferon-stimulated genes: what do they all do? Annu. Rev. Virol. 6 , 567–584 (2019). 31283436 10.1146/annurev-virology-092818-015756

82. Chemudupati M From APOBEC to ZAP: Diverse mechanisms used by cellular restriction factors to inhibit virus infections Biochim. Biophys. Acta Mol. Cell Res. 2019 1866 382 394 30290238 10.1016/j.bbamcr.2018.09.012 PMC6334645 Chemudupati, M. et al. From APOBEC to ZAP: Diverse mechanisms used by cellular restriction factors to inhibit virus infections. Biochim. Biophys. Acta Mol. Cell Res. 1866 , 382–394 (2019). 30290238 10.1016/j.bbamcr.2018.09.012 PMC6334645

83. Sasa, N. & Kojima, S. Scripts for ‘Blood DNA virome associates with autoimmune diseases and COVID-19’. Zenodo 10.5281/zenodo.13953637 (2024). 10.1038/s41588-024-02022-z PMC11735405 39753770