Deer sample collection and study area

Between 1 November and 31 December 2021, adult and yearling free-ranging white-tailed deer were sampled as part of the Ontario Ministry of Natural Resources and Forestry’s (MNRF) annual Chronic Wasting Disease (CWD) surveillance programme. Samples were collected from hunter-harvested deer in Southwestern and Eastern Ontario and included nasal swabs and RPLNs. Samples were collected by staff wearing masks and disposable gloves. Nasal swabs were stored in individual 2 ml tubes with ~1 ml of universal transport medium (Sunnybrook Research Institute (SRI)), and RPLN tissues were stored dry in 2 ml tubes. After collection, samples were immediately chilled on ice packs then transferred to a −20 °C freezer where they were held for up to 1 week. Samples were then transferred to a −80 °C freezer where they were held until analysis. Location, date of harvest, and demographic data (age/sex) were recorded for each animal when available.

PCR screening and detection

RNA extractions and PCR testing of samples collected from deer were performed at the SRI in Toronto, Ontario. RNA extractions were conducted with 140 µl of nasal swab sample spiked with Armored RNA enterovirus (Asuragen; https://www.asuragen.com) via the Nuclisens EasyMag using Generic Protocol 2.0.1 (bioMérieux Canada) according to the manufacturer’s instructions; RNA was eluted in 50 µl. Tissue samples were thawed, weighed, minced with a scalpel and homogenized in 600 µl of lysis buffer using the Next Advance Bullet Blender (Next Advance) and a 5 mm stainless steel bead at 5 m s−1 for 3 min. RNA from 30 mg tissue samples was extracted using Specific Protocol B 2.0.1 via Nuclisens EasyMag; RNA was eluted in 50 µl. RT–PCR was performed using the Luna Universal Probe One-Step RT–qPCR kit (New England BioLabs, NEB; https://www.neb.ca). A SARS-CoV-2 5′ UTR and E specific multiplex RT–PCR were used for RNA detection47. Quantstudio 3 software (Thermo Fisher Scientific; https://www.thermofisher.com) was used to determine the cycle threshold (Ct). All samples were run in duplicate, and samples with Ct <40 for both RT–PCR targets and Armored RNA enterovirus in at least one replicate were considered presumed positive. For tissue samples, the presence of inhibitors was assessed by a 1:5 dilution of one of the replicates. Samples were considered inconclusive if no Armored enterovirus was detected or if only one RT–PCR target was detected and re-extracted for additional analysis. Samples were considered indeterminate if inconclusive after re-extraction or if no original material was left. Presumed positive samples were further analysed for human RNAse P to rule out potential human contamination9. Original material from presumed positive samples detected at SRI were sent to the Canadian Food Inspection Agency (CFIA) for confirmatory PCR testing. The MagMax CORE Nucleic Acid Purification Kit (Thermo Fisher Scientific) and the automated KingFisher Duo Prime magnetic extraction system was used to extract total RNA spiked with Armored RNA enterovirus. The enteroviral armored RNA was used as an exogenous extraction control. A SARS-CoV-2 E and nucleocapsid (N) specific multiplex RT–PCR was used for confirmatory RNA detection7. Master mix for qRT–PCR was prepared using TaqMan Fast Virus 1-step Master Mix (Thermo Fisher Scientific) according to the manufacturer’s instructions. Reaction conditions were 50 °C for 5 min, 95 °C for 20 s and 40 cycles of 95 °C for 3 s then 60 °C for 30 s. Runs were performed by using a 7500 Fast Real-Time PCR System (Thermofisher, ABI). Samples with Ct <36 for both RT–PCR targets were considered positive.

WGS

WGS was performed at both SRI and CFIA using independent extractions and sequencing methods. At SRI, DNA was synthesized from extracted RNA using 4 μl LunaScript RT SuperMix 5× (NEB) and 8 μl nuclease free water, and was added to 8 μl extracted RNA. Complementary DNA synthesis was performed under the following conditions: 25 °C for 2 min, 55 °C for 20 min, 95 °C for 1 min and holding at 4 °C.

The ARTIC V4 primer pool (https://github.com/artic-network/artic-ncov2019) was used to generate amplicons from the cDNA. Specifically, two multiplex PCR tiling reactions were prepared by combining 2.5 μl cDNA with 12.5 μl Q5 High-Fidelity 2× Master Mix (NEB), 6 μl nuclease-free water and 4 μl of respective 10 μM ARTIC V4 primer pool (Integrated DNA Technologies). PCR cycling was then performed in the following manner: 98 °C for 30 s followed by 35 cycles of 98 °C for 15 s and 63 °C for 5 min.

PCR reactions were then both combined and cleaned using 1× ratio sample purification beads (Illumina) at a 1:1 bead to sample ratio. The amplicons were quantified using the Qubit 4.0 fluorometer using the 1× dsDNA High Sensitivity (HS) Assay Kit (Thermo Fisher Scientific) and sequencing libraries prepared using the Nextera DNA Flex Prep kit (Illumina) as per the manufacturer’s instructions. Paired-end (2 × 150 bp) sequencing was performed on a MiniSeq with a 300-cycle reagent kit (Illumina) with a negative-control library with no input SARS-CoV-2 RNA extract.

WGS performed at CFIA used extracted nucleic acid quantified using the Qubit RNA HS Assay Kit on a Qubit Flex Fluorometer (Thermo Fisher Scientific). Eleven microlitres or 200 ng of total RNA was subject to DNase treatment using the ezDNase enzyme (Thermo Fisher Scientific) according to the manufacturer’s instructions. DNase-treated RNA was then used for library preparation and target sequence capture according to the ONETest Coronaviruses Plus Assay protocol (Fusion Genomics48). The enriched libraries were quantified using the Qubit 1× dsDNA HS Assay Kit on a Qubit Flex Fluorometer (Thermo Fisher Scientific) and subsequently pooled in equimolar amounts before fragment analysis on 4200 TapeStation System using the D5000 ScreenTape Assay (Agilent). The final pooled library was sequenced on an Illumina MiSeq using a V3 flowcell and 600 cycle kit (Illumina).

Human specimens are received at Public Health Ontario Laboratory for routine SARS-CoV-2 diagnostic testing (RT–PCR) from multiple healthcare settings, including hospitals, clinics and coronavirus disease 2019 (COVID-19) assessment centres. The human sample (ON-PHL-21-44225) was sequenced at Public Health Ontario Laboratory using an Illumina-based ARTIC V4 protocol (https://doi.org/10.17504/protocols.io.b5ftq3nn), similar to the deer sequencing methods. Briefly, cDNA was synthesized using LunaScript reverse transcriptase (NEB). Amplicons were generated with premixed ARTIC V4 primer pools (Integrated DNA Technologies). Amplicons from the two pools were combined, purified with AMPure XP beads (Beckman Coulter) and quantified. Genomic libraries were prepared using the Nextera XT DNA Library Preparation Kit (Illumina), and genomes were sequenced as paired-end (2 × 150 bp) reads on an Illumina MiSeq instrument.

Genomic analysis

Paired-end illumina reads from ARTIC V4 and Fusion Genomics sequencing were initially analysed separately with the nf-core/viralrecon Nextflow workflow (v2.3) (refs. 49,50,51) that ran: FASTQC (v0.11.9) (ref. 52) read-level quality control, fastp (v0.20.1) (ref. 53) quality filtering and adapter trimming, Bowtie2 (v2.4.2) (ref. 54) read mapping to Wuhan-Hu-1 (MN908947.3) (ref. 55) SARS-CoV-2 reference, Mosdepth (v0.3.1) (ref. 56)/Samtools (v.1.12) (ref. 57) read mapping statistics calculation, iVar (v1.3.1) (ref. 58) ARTIC V4 primer trimming, variant calling and consensus generation; SnpEff (v5.0) (ref. 59)/SnpSift (v4.3t) (ref. 60) for variant effect prediction and annotation; and Pangolin (v3.1.20) (ref. 61) with PangoLEARN (2022-01-05), Scorpio (v0.3.16) (ref. 62), and Constellations (v.0.1.1) was used for PANGO lineage63 assignment. iVar primer trimmed soft-clipped read alignments were converted to hard-clipped alignments with fgbio ClipBam (http://fulcrumgenomics.github.io/fgbio/). Reads from hard-clipped BAM files were output to FASTQ files with ‘samtools fastq’. nf-core/viralrecon was re-run in amplicon mode without iVar primer trimming on the combined Fusion Genomics and ARTIC V4 primer trimmed FASTQ reads to generate the variant calling, variant effect and consensus sequence results used in downstream analyses. Additional quality-control steps to check for negative-control validity, drop-out, sample cross-contamination and excess ambiguity were performed using ncov-tools v1.8.0 (ref. 64). The mutations identified by the Nextclade (v1.10.2) (ref. 65) with 2022-01-05 database and xlavir (v0.6.1) report were manually searched in outbreak.info’s ‘Lineage | Mutation Tracker’ (on 2 February 2022) (ref. 66) to get information on the prevalence of observed mutations globally and within Canada. Mutations were also investigated for presence in specific lineages including VOCs, Michigan mink samples and other animal samples. Finally, mutations were searched in GISAID (on 2 February 2022) to tally the number of non-human hosts each mutation had been observed in.

Some limitations in genome quality and coverage existed that may have resulted in failure to detect additional mutations. All B.1.641 samples had missing terminal domains and contained internal regions with no or low coverage when sequenced using the ARTIC v4 amplicon scheme. This is a widespread issue that may explain the rarity of the 3′ proximal ORF10:L37F in GISAID. Significantly in our samples this meant there was no or <10× coverage in all five deer-derived sequences from ~27000 to 27177 (drop-out of ARTICv4 amplicons 90-91), which includes regions of the M gene. However, by combining the ARTIC v4 sequencing with additional sequencing using probe-based enrichment we were able to compensate for this drop-out and generate high coverage and completeness (<100 positions with no coverage in all deer and <100 positions with <10× coverage in 3/5 deer genomes; Supplementary Table 7).

Phylogenetics

To evaluate possible sampling biases due to the poorly defined and diverse B.1 and B.1.311 lineages and select closely related publicly available sequences for further phylogenetic analysis, a phylogenetic placement analysis based on UShER (https://genome.ucsc.edu/cgi-bin/hgPhyloPlace)34 was performed using the 7,217,299 sample tree (derived from UShER placement of GISAID, GenBank, COG-UK and CNCB onto 13-11-20 sarscov2phylo ML tree) via the SHUShER web-portal (shusher.gi.ucsc.edu). Phylogenetic analyses were performed using CFIA-NCFAD/scovtree Nextflow workflow (v1.6.0) (https://github.com/CFIA-NCFAD/scovtree/) with the consensus sequences contextualized with closely related sequences identified by UShER and randomly sampled representative sequences from major WHO SARS-CoV-2 clades from GISAID12 (downloaded 10 February 2022). This workflow generated a multiple sequence alignment using Nextalign CLI (v1.10.1) (ref. 65) and inferred an ML phylogeny using IQ-TREE (v2.2.0_beta)67 using the general time-reversible (GTR) model for visualisation with Phylocanvas68 via shiptv (v0.4.1) (https://github.com/CFIA-NCFAD/shiptv) and ggtree69. Divergence times for the inferred global ML topology were estimated using BEAST v1.10.4 (ref. 70). This analysis used a coalescent model with constant population size, a Hasegawa–Kishono–Yano substitution model with four discrete gamma categories, and a log-normal distributed strict molecular clock rate of 9.5 × 10−4 substitutions per site per year (based on a tip-to-root regression performed using TempEst71). Internal node heights and root node height were sampled by BEAST over 50 million MCMC generations (recorded every 1,000 iterations) before collation using BEAST’s LogCombiner. The final maximum clade credibility tree was generated using BEAST’s TreeAnnotator with node heights set to median values before final visualization in FigTree and Inkscape.

A subset of 157 taxa from an ancestral clade of B.1.641 were selected from the global phylogenetic tree shown in Fig. 2 to generate the phylogenetic tree shown in Fig. 3. Multiple sequence alignment of this subset of sequences was performed with MAFFT (v7.490) (ref. 72). An ML phylogenetic tree was inferred with IQ-TREE (v2.2.0_beta) using the GTR model and 1,000 UFB replicates73. Nextclade (v1.10.2) analysis was used to determine amino acid mutations and missing or low/no coverage regions from the sample genome sequences. Amino acid mutation profiles were determined relative to the B.1.641 samples, discarding mutations that were not present in any of the Ontario samples. Taxa with duplicated amino acid mutation profiles were pruned from the tree, keeping only the first observed taxa with a duplicated profile.

Recombination analyses were performed using 3Seq (v1.7) (ref. 30) and Bolotie (e039c01) (ref. 31). Specifically, 3Seq was executed with B.1.641 sequences and the most recent example of each lineage found in Canada and closest samples in GISAID in subtree (n = 595). Bolotie was executed using the B.1.641 sequences and two datasets, the provided pre-computed 87,695 probability matrix and a subsample of the earliest and latest example of each lineage in GISAID with all animal-derived samples and closest UShER samples (n = 4,688). Additionally, HyPhy’s32 (v2.5.31) genetic algorithm recombination detection method33 was applied to local alignment and ML phylogeny (Fig. 3) for all possible sites. Phylogenies were inferred using IQTree for segments either side of the identified putative breakpoint, and the B.1.641 clade and local topology was unchanged. Sequence statistics such as C > T rate were directly calculated from nextclade results (v1.10.2 with 2022-01-05 database). Additional figures were generated and annotated using BioRender74 and Inkscape75.

A phylogenetic approach with HyPhy32 (v2.5.31) was used to investigate signatures of selection within B.1.641 relative to the wider B.1 background lineage. To ensure necessary genomic completeness for codon alignment, all B.1 sequences in GISAID (as of 8 March 2022) were filtered to those with <0.1% N with full date information. Genomes with 0% N were removed to avoid biases from consensus workflows that replace undetermined sequence with reference genome. From this, all animal-derived (49 mink and 1 cat) sequences and 100 randomly sampled human B.1 sequences were extracted (n = 150). Finally, the WH0-1 reference genome was added to this alignment along with the five complete Ontario deer-derived genomes, associated human sequence, and the two most closely related Michigan human samples (MI-MDHSS-SC23517 and M-MDHSS-SC22669). Virulign76 (v1.0.1) was then used to generate codon alignments for E, M, N, S, ORF1ab, ORF3a, ORF6, ORF7a, ORF7b, ORF8 and ORF10 genes relative to the Wuhan-Hu-1 (MN908947.3) reference. ML phylogenies were inferred for these alignments using raxml-ng (v1.0.2) (ref. 77) with the GTR model and three parsimony-based starting trees. These phylogenies were manually inspected and rooted on Hu-1, and the B.1.641 branches were labelled using phylowidget37,78. Genes for which the phylogeny did not have a resolved B.1.641 clade (ORF7a and ORF7b) or a viable codon alignment without any internal stop codons (ORF8 and ORF10) were excluded from further analyses. For each gene, signatures of positive selection were evaluated using HyPhy’s adaptive branch-site random effects likelihood (aBSREL) method36 and signatures of gene-wide episodic diversification were evaluated using the branch-site unrestricted statistical test for episodic diversification (BUSTED) method37 with ten starting points. Finally, evidence of intensification or relaxation of selection was investigated using the RELAX method79 with ten starting points and synonymous rate variation. Additional code for divergence dating, recombination and selection analyses can be found under https://doi.org/10.5281/zenodo.7086599.

Analysis of mutational spectrum

The mutational spectra were created using a subset of 3,645 sequences used to create the high-quality global phylogeny. Included in this dataset are the seven unique deer samples from Ontario (samples 4538, 4534, 4662, 4649, 4581, 4645 and 4658; the five high-quality genomes plus two genomes with lower coverage), three white-tailed deer samples from Quebec (samples 4055, 4022 and 4249) and one human sample from Ontario (ON-PHL-21-44225), and any remaining human, mink and deer sequences from North America. The counts for each type of nucleotide change, with respect to the reference strain, were compiled and used to create a 12-dimensional vector. This subset was then filtered to remove sequences with fewer than 15 nucleotide changes. The counts were converted into the mutation spectrum by simply dividing each count by the sum of the counts in each sample10,41. As the mutation spectrum summarizes how host factors act upon the genome of SARS-CoV-2, it can potentially be used as an independent source of evidence supporting the evolution of the virus in a different host10,41. To investigate this, we conducted experiments using a dbWMANOVA80. If a significant difference between hosts was detected, a pairwise distance-based Welch t-test was used to identify which pair of hosts differed81. This approach was used because it is more robust on unbalanced and heteroscedastic data80,81. In the first experiment, samples from all lineages were used. The second experiment used only the subset of samples belonging to Nextstrain clade 20C (which contains the B.1.641 sequences). An analysis using the entire Pangolin B.1 lineage was not appropriate since this lineage also includes Nextstrain clades 20A and 20B and the inclusion of these samples could potentially obfuscate important patterns since the results would have to be interpreted in a context wider than necessary.

Codon usage analysis

Consensus sequences of SARS-CoV-2 samples from this and previous studies and additional sequences gathered from public databases were used. The sequences include the reference SARS-CoV-2 Wuhan-Hu-1 (NCBI NC_045512), SARS-CoV-2 mink/Canada/D/2020 (GISAID EPI_ISL_717717), SARS-CoV-2 mink/USA/MI-20-028629-004/2020 (GISAID EPI_ISL_2834697), Cervid atadenovirus A 20-5608 (NCBI OM470968)82, EHDV serotype 2/strain Alberta (NCBI AM744997 - AM745006), epizootic haemorrhagic disease virus, EHDV serotype 1/New Jersey (NCBI NC_013396 - NC_013405), EHDV 6 isolate OV208 (NCBI MG886400 - MG886409) and elk circovirus Banff/2019 (NCBI MN585201) (ref. 83) were imported into Geneious (v.9.1.8) (ref. 84). Annotations for the coding sequences of SARS-CoV-2 samples were transferred from the reference sequence SARS-CoV-2 Wuhan-Hu-1 (NC_ 045512) using the Annotate from Database tool. The coding sequences were extracted using the Extract Annotations tool for all viral sequences. An annotated file of the coding sequences for the Odocoileus virginianus texanus isolate animal Pink-7 (GCF_002102435.1) genome was downloaded from NCBI (https://ftp.ncbi.nlm.nih.gov/genomes/all/annotation_releases/9880/100/GCF_002102435.1_Ovir.te_1.0/). Coding sequences were input into CodonW (http://codonw.sourceforge.net/) with settings set to concatenate genes and output to file set to codon usage. Codon usage indices were set to the effective number of codons (ENc), GC content of gene (G + C), GC of silent third codon position (GC3s), silent base composition, number of synonymous codons (L_sym), total number of amino acids (L_aa), hydrophobicity of protein (Hydro) and aromaticity of protein (Aromo).

Virus isolation

For virus isolation, T25 flasks were seeded to confluency 1 day before infection with cathepsin L knock-out Vero E6 cells overexpressing TMPRSS2. The following day, swab samples were vortexed and spun down and 200 μl of the swab samples medium was combined with 16 μg ml−1 working concentration of TPCK-treated trypsin (NEB), 2× A/A/p/s antifungal/antibiotic solution (Wisent) and a 0.1% working concentration of BSA (Thermo Fisher Scientific) and added to the cell monolayer after removal of the medium. Samples were subjected to a 45 min adsorption with rocking every 5 min, after which the inoculum was removed and discarded and the monolayer was either washed once with 2 ml of D1 to remove blood cells present in the samples (4581, 4649 and 4676) or not washed (4645, 4658 and 4662) and 5 ml of DMEM with 1% FBS and antibiotics was added to the flask and incubated at 37 °C with 5% CO2. At 4 days post-infection, samples with visible cytopathic effect (partial, 50% or less rounded or detached cells) were harvested followed by collection and centrifugation at 4,000g for 10 min at 20 °C. The harvested supernatants were aliquoted and stored at −80 °C or inactivated and removed from the CL3 laboratory and RNA was extracted with the QIAamp Viral RNA Mini Kit (Qiagen), and stored at −20 °C until downstream analyses were carried out. All infectious work was performed under biosafety level 3 conditions.

Codon-optimized spike constructs, cells, sera and antibodies

Expression constructs of S mutants corresponding to samples 4581/4645 (S:H49Y, S:T95I, S:Δ143-145InsD, S:F486L, S:N501T, S:D614G), 4658 (S:T22I, S:H49Y, S:T95I, S:Δ143-145InsD, S:S247G, S:F486L, S:N501T, S:D614G) and ON-PHL-21-44225 (S:H49Y, S:T95I, S:Δ143-145InsD, S:F486L, S:N501T, S:Q613H, S:D614G) were generated by overlapping PCR as described previously85. S:D614G and S Omicron (BA.1) constructs were described elsewhere86. All constructs were cloned in pCAGGS and verified by Sanger sequencing.

HEK293T cells (ATCC) were cultured in DMEM supplemented with 10% FBS (Sigma), 100 U ml−1 penicillin, 100 µg ml−1 streptomycin and 0.3 mg ml−1l-glutamine (Invitrogen) and maintained at 37 °C, 5% CO2 and 100% relative humidity.

Serum samples were obtained from consenting participants in several cohort studies with sample collection and sharing for this analysis approved by the Sinai Health System Research Ethics Board (#22-0030-E). Plasma of SARS-CoV-2 naïve, naïve-vaccinated (28–40 days after two or three doses of BNT162b2) and unvaccinated SARS-CoV-2 Delta previously infected donors was collected (Supplementary Table 8), heat inactivated for 1 h at 56 °C, aliquoted and stored at −80 °C until use. The conformation-independent monoclonal anti-S2 CV3-25 from a convalescent individual was described and produced as described previously87,88. The goat anti-human IgG conjugated with Alexa Fluor-647 was purchased from Invitrogen (A21445).

Spike binding assays

Hek293T cells seeded in a 10 cm Petri dish at 70% confluency were transfected with 10 µg of SARS-CoV-2 spike protein plasmid, 1 µg of lentiviral vector bearing green fluorescent protein (GFP) (PLV-eGFP) (gift from Pantelis Tsoulfas, Addgene plasmid number 36083) (ref. 89) using Jetprime transfection reagent (Polyplus, catalogue number 101000046) according to the manufacturer’s instructions. At 16 h post-transfection, the cells were stained with sera samples (1:250 dilution) for 45 min at 37 °C. Alexa Fluor-647-conjugated goat anti-human IgG (H + L) was used to detect plasma binding of the treated cells following 1 h incubation at room temperature. Samples were washed once with PBS, fixed in 1% paraformaldehyde and acquired using BD LSR Fortessa Flow cytometer (BD Biosciences). The seropositivity threshold was defined on the basis of mean fluorescence intensity (MFI) for naïve samples plus three standard deviations. The data were normalized by surface expression on the basis of the MFI of the monoclonal antibody CV3-25 (5 μg ml−1). The data analysis was performed using FlowJo 10.8.1 (Extended Data Fig. 6). For each set of sera, binding was compared across samples using Welch’s (unequal variance) one-way ANOVA procedure and a post-hoc Tukey’s honestly significant difference test (using a family-wise error rate of 0.05) via the statsmodel library (v0.14.0) (ref. 90).

Pseudotype production and neutralization assays

HEK293T seeded in 10 cm dishes were co-transfected with lentiviral packaging plasmid psPAX2 (gift from Didier Trono, Addgene number 12260), lentiviral vector pLentipuro3/TO/V5-GW/EGFP-Firefly Luciferase (gift from Ethan Abela, Addgene number 119816) and plasmid encoding the indicated S construct at a 5:5:1 ratio using jetPRIME transfection reagent according to the manufacturer’s protocol. Twenty-four hours post-transfection, media were changed, and supernatants containing lentiviral pseudotypes were harvested 48 h post-transfection, filtered with a 0.45 µM filter and stored at −80 °C until use.

HEK293T stably expressing human ACE2 (293T-ACE2, kind gift of Hyeryun Choe, Scripps Research91) were seeded in poly-d-lysine-coated 96-well plates. The next day, supernatants containing lentiviral pseudotypes were incubated with sera (serially diluted by five-fold, from 1:50 to 156,250) for 1 h at 37 °C and then added to cells in the presence of 5 µg ml−1 polybrene. Seventy-two hours later, media were removed, and cells were rinsed in phosphate-buffered saline and lysed by the addition of 40 µl passive lysis buffer (Promega) followed by one freeze–thaw cycle. A Synergy Neo2 Multi-Mode plate reader (BioTek) was used to measure the luciferase activity of each well after the addition of 50–100 µl of reconstituted luciferase assay buffer (Promega) as per the manufacturer’s protocol. Neutralization ID50 was calculated using Graphpad Prism and represents the plasma dilution that inhibits 50% of pseudotype transduction in 293T-ACE2. For each set of sera, neutralization was compared across samples using Welch’s (unequal variance) one-way ANOVA procedure and a post-hoc Tukey’s honestly significant difference test (using a family-wise error rate of 0.05) via the statsmodel library (v0.14.0) (ref. 90).

Reporting summary

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