Multi-generational koala pedigree analysis reveals rapid changes in heritable provirus load associated with life history traits

Domain: chemistry
Published: 2026-01-09
Editorial outcome: minor_revision

Abstract
Retroviruses that colonize the host germline can be passed on as inherited genetic variants. The koala (Phascolarctos cinereus) is currently experiencing germline colonization by two retroviruses, the koala retrovirus (KoRV) and phaCin-β. We analyze the integration site segregation and diversity of endogenous KoRV, phaCin-β, and the related phaCin-β-like in 111 pedigreed koalas from the San Diego Zoo Wildlife Alliance and seven European Zoos. The use of multigenerational pedigrees and the inclusion of health information for each individual koala reveal elimination of retroviruses from proto-oncogenes and the generation and spread of new germline integrations. Seven-hundred-and-fourteen integrations do not persist in the living population. For the 55 triads examined, 21 unique integrations identified in individual koalas are absent in their parents. Retroviral integrations associated with leukemia, fertility, and longevity are used to estimate genetic risk scores and develop a longevity breeding index to minimize neoplasia risk in the captive koala population.

Introduction
Retroviruses are single stranded RNA viruses that integrate a reverse transcribed DNA provirus copy of their genome into host nuclear DNA during viral replication. If this occurs in the germline, the resulting provirus can be transmitted in a Mendelian fashion. Such germline proviral sequences are called endogenous retroviruses (ERVs). The koala (Phascolarctos cinereus) has experienced multiple recent and ongoing germline colonizations by retroviruses, including the Koala Retrovirus (KoRV)1 and phaCin-β2, which may still co-exist with exogenous (infectious) retroviral lineages whereas other endogenized retroviruses, such as phaCin-β-like2 and the Phascolarctos endogenous retroelement (PhER)3 do not. KoRV is a gammaretrovirus most closely related to the gibbon ape leukemia viruses (GALVs). It was first characterized from koala lymphoma tissues and in mitogen-stimulated peripheral blood mononuclear cell cultures from 163 koalas1. KoRV was subsequently determined to be a vertically transmitted ERV that had only recently colonized the koala, as populations that were KoRV free were identified4. Based on the DNA mutation rate per base pair per year of humans (2×10−9) and mice (4.5×10−9), KoRV was estimated to have begun colonizing the koala germline between 22,200 and 49,900 years ago5. The phaCin-β and phaCin-β-like ERVs were estimated to have begun colonizing the koala germline between 0.7 and 1.5 Mya, and between 2.4 and 5.5 Mya, respectively2.

The number of KoRV and phaCin-β loci per individual varies among populations, in contrast with most other vertebrate ERVs where the number remains constant. The variation in the number of integrations is consistent with both ERVs being in the early stages of endogenization. Koalas in northern Australia generally carry more endogenous KoRVs (enKoRVs), with numbers ranging from 54 to 116 per individual. In the south, koalas carry fewer enKoRVs, ranging from 8 to 23 per individual6. The absence of enKoRVs was previously reported on Kangaroo Island, south of South Australia4. A less pronounced north-south trend was reported for phaCin-β, ranging from 48 to 90 phaCin-β integrations in the north to range of 67–81 phaCin-β in the south6. The number of ERVs per individual for phaCin-β-like ranged from 46 to 64 polymorphic phaCin-β-like integrations in the north, and from 53 to 67 phaCin-β-like integrations in the south6. The more widespread stability of copy number and location of phaCin-β-like is consistent with germline colonization predating both phaCin-β and enKoRV.

There are health consequences associated with proviral integrations. For the koala, it has been reported up to 60% die from KoRV associated leukemia and lymphoma in captive populations1,7. Provirus integrations into enhancer or promoter sequences can activate adjacent gene expression. Proviruses can also influence mRNA 3’-end processing when downstream of a gene or can inactivate genes by insertional disruption8,9. Tumor tissues accumulate somatic KoRV integrations which are significantly overrepresented in proximity to cellular proto-oncogenes10.

Most studies of KoRV and other koala retroelement endogenization have been conducted at the population level without pedigree information. Using pedigrees, it is possible to determine the frequency of novel germline integrations, how and which potentially deleterious ERVs are lost, and examine whether positive selection for protective ERVs is detectable. Pedigree data allow for the determination of occurrence and dosage of ERV alleles in families. For example, new ERV integrations can be unique to offspring but absent in their parents. Additionally, metadata on life and health histories of each individual can associate information on relatedness with impacts of specific ERVs on the fitness and survival of individuals.

In this study, we characterize the diversity of enKoRV, phaCin-β, and phaCin-β-like ERVs in koalas from the San Diego Zoo Wildlife Alliance (SDZWA), utilizing the complete pedigree information to investigate the evolutionary dynamics of ERVs for up to four koala generations. Comparisons with wild koalas indicate that the SDZWA koalas are genetically more closely related to northern Australian koalas than to koalas from South Australia. In addition, there is evident genetic relatedness between SDZWA koalas and koalas from European zoos (EUZ), as expected due to translocations between zoos in North America and Europe. We identify instances of ERVs that do not persist across generations, and the appearance of new KoRV and phaCin-β proviruses in the germline. Additionally, we identify candidate markers, including ERVs and single nucleotide polymorphisms (SNPs), associated with traits, such as leukemia, fertility, and longevity, making them important markers for consideration in captive breeding programs. This analysis forms the basis of determining genetic risk scores for the captive SDZWA koala population.

Results
Koala pedigrees from the SDZWA and EUZ
Koalas were first brought to North America between 1976 and 1981, comprising 14 founders from the Lone Pine Koala Sanctuary. After 1990, occasional imports occurred, including nine koalas from Japan, Currumbin, and Coomera. SDZWA currently keeps 30 koalas, the largest captive population outside Australia. Pedigree information for all animals comprises information about parentage data, sex, date of birth, birth location, date of death, and location and cause of death (if deceased). In addition to blood samples from all living animals, cryopreserved tissues of koalas collected over 40 years were available for DNA extraction for a total of 91 koala individuals, including three individuals with different tissue replicates. This sample collection represents four generations of koalas, dating back to animals born in 1978 (Fig. 1a), including one of the original founders (Studbook number (SB) 106). It comprised 46 triads, including 19 F2 and five F3 generations (Supplementary Data 1).

Twenty koalas from seven EUZ were also used in the investigation of enKoRV only (Supplementary Data 2). These comprised nine triads, including two F2 generations. Individuals from EUZ are related to individuals from SDZWA, with two individuals having parents from SDZWA (SB 296 and 367), and individuals having grandparents from SDZWA (SB 201, 206, 244, 264, 399, 403, 436, and 537).

enKoRV and phaCin-β polymorphic integrations
A total of 2075 ERV integrations (1232 enKoRVs, 485 phaCin-β, and 358 phaCin-β-like) were characterized across 91 koala individuals from the SDZWA population. On average, each koala had 471 ERV integrations, ranging from 379 to 510. On average, there were 72 enKoRV, 48 phaCin-β, and 351 phaCin-β-like per koala. From 1978 to 2022, the number of ERV integrations per individual was stable (Fig. 1a).

Most of the phaCin-β-like, which are the oldest of these ERVs, were fixed in the population except two integrations that continue to segregate in the living population (Fig. 1b). Alternative Allele Frequency (AAF) was calculated considering the pre-integration site as reference allele and the ERV as the alternative allele. One of the segregating phaCin-β-like was intergenic and observed in 52 individuals (AAF=0.27), including the founder individual, while the other was close to fixation (AAF=0.97), located in an intron of the oncogene SLC12A7 (solute carrier family 12 member 7). SLC12A7 encodes a kinase binding protein, with overexpression associated with aggressive adrenocortical carcinoma in humans11. Both enKoRV and phaCin-β integrations displayed significantly lower AAFs than phaCin-β-like, except for a few phaCin-β, which were observed at higher frequencies that were mostly homozygous (Supplementary Figs. 1 and 2).

Germline mutation rate and retroviral germline colonization timeline
SNPs present in koala offspring (joeys) and not present in their parents in 46 triads were used to estimate the koala mutation rate. After strict filtering of new SNPs, to avoid false-positives as described in the Methods, a mutation rate per base pair per generation of 1.03×10−8 was estimated. This is slightly lower than the estimated rate for humans (1.13×10−8) and cattle (1.17×10−8). The mutation rate per base pair per year was estimated to be 1.6×10−9, lower than the annual estimated rate for humans (2.0×10−9) and mice (4.5×10−9), used previously to estimate the initial germline colonization of KoRV of not more than 50,000 years ago, phaCin-β 1.5 million years ago (Mya) and phaCin-β-like and 5.5 Mya.

Using the lower koala mutation rate, we revised KoRV germline colonization start to not more than 312,191 years ago, concomitant with the emergence of modern koalas12. Germline colonization estimates for phaCin-β and phaCin-β-like with our revised koala mutation rate resulted in colonization times of 1.9 and 6.9 Mya, respectively. Nevertheless, younger estimated dates cannot be discarded due to the possibility of gene conversion or other recombination events13.

Origin of SDZWA and EUZ ERVs
In total, 75 wild koalas from the Koala Genome Survey (Supplementary Data 3 and Supplementary Fig. 3) were compared with captive koalas from the SDZWA. Twenty-five individuals derived from each state, Queensland (QLD), New South Wales (NSW) and Victoria (VIC), carrying 2616 enKoRVs, 1446 phaCin-β, and 359 phaCin-β-like. Of these 4421 ERVs, 383 were shared between VIC and SDZWA, 582 between NSW and SDZWA, and 837 between QLD and SDZWA (Fig. 2a). The number of shared ERVs was inflated by the 356 phaCin-β-like common to all populations (Supplementary Fig. 5). Fewer shared integrations were observed for enKoRV (Fig. 2b) and phaCin-β (Supplementary Fig. 6).

Notably, the low number of enKoRVs in VIC (n=110) is consistent with other studies demonstrating low KoRV prevalence in VIC and South Australia (Fig. 2b)6. Only two enKoRVs from VIC were shared with SDZWA. However, in comparison to SDZWA, the AAFs of phaCin-β in VIC were higher (Supplementary Fig. 7).

The number of ERVs per individual differed among populations. Koalas from QLD had on average 72 enKoRVs, 70 phaCin-β, and 326 phaCin-β-like. Similarly, NSW koalas had on average 73 enKoRVs, 82 phaCin-β, and 336 phaCin-β-like. The average numbers were similar to koalas from SDZWA, except for phaCin-β, which was more frequent in wild koalas. A different pattern was observed in VIC with an average of 8 enKoRVs, 71 phaCin-β, and 337 phaCin-β-like per koala. This result corroborates previous findings on the low prevalence of enKoRV in VIC and South Australia6. Two koalas from Bawbaw and Strathbogies (VIC) had only one enKoRV each. Nevertheless, differences in sequence coverage may have influenced detection of ERVs in these individuals. A significant Pearson correlation between number of phaCin-β-like per individual and coverage was observed for NSW (r=0.72, p-value=4.4×10-5) and VIC (r=0.69, p-value=1.2×10-4), and of enKoRV per individual and coverage for VIC (r=0.59, p-value=1.9×10-3).

The 130 shared enKoRV between QLD and SDZWA reflects the founding history of the SDZWA koala population, with 14 individuals coming from the Lone Pine Koala Sanctuary in Brisbane (QLD). Imports of individuals from NSW have occurred, accounting for the shared ERVs between SDZWA and NSW koala populations.

While koalas from SDZWA cluster separately in the PCA (Fig. 2c), three individuals show close relationship with QLD and NSW koalas, including two individuals transferred to SDZWA (SB 267 and 438) and one offspring (SB 153) of a founder koala (SB 106). Individuals from QLD and NSW show overlap of ERV integrations, likely explained by gene flow between the regions. Individuals from VIC form a distinct cluster, emphasizing the isolation of koalas from this region. However, one individual from southern NSW, Narrandera, clustered with VIC koalas rather than NSW koalas, explained by the mixed ancestry of koalas in this region reintroduced from QLD and VIC in 1970s14.

Twenty koalas from different EUZ were targeted for KoRV enrichment15, resulting in the discovery of 2310 proviral integrations that we characterized as 553 enKoRVs and 1757 somatic KoRV integrations (see Methods). On average, each individual had 97 enKoRVs, ranging from 42 to 141, and 88 somatic KoRVs, ranging from 3 to 241. While most enKoRVs in koalas from EUZ were unique, 20% (244) were shared with SDZWA. Similarly, enKoRVs were shared between koalas from QLD and EUZ and between koalas from NSW and EUZ.

The enKoRV with the highest frequency in all populations (AAF=0.43) was the intronic enKoRV located in the oncogene SLC29A1 (solute carrier family 29 member 1-Augustine blood group). Another enKoRV was located upstream of SLC29A1 (AAF=0.07). Further enKoRVs high AAFs (AAF>0.10) were located inside intronic regions of the genes PCCA (propionyl-CoA carboxylase subunit alpha), encoding the alpha subunit of the heterodimeric mitochondrial enzyme Propionyl-CoA carboxylase, and PTPRT (protein tyrosine phosphatase receptor type T), encoding a tumor suppressor protein.

For phaCin-β, regions with high AFFs were located in genes, such as TMEM169 (transmembrane protein 169), KCTD21 (potassium channel tetramerization domain containing 21), KCNK2 (potassium two pore domain channel subfamily K member 2), ORC4 (origin recognition complex subunit 4), and SNRPC (small nuclear ribonucleoprotein polypeptide C), fixed in both captive and wild populations (AAF=1.00).

Non-singleton enKoRVs (1595 enKoRVs shared by at least two individuals) were analyzed using PCA including EUZ (Fig. 2d). With EUZ data, only the detection or absence of enKoRV was used, since the sequencing data for EUZ was targeted for KoRV. Captive and wild koalas formed distinct clusters, with a clear relationship between koalas from SDZWA and EUZ.

Although most ERVs identified in our study were intergenic (4479), we identified 1273 intronic ERVs. Further, 26 enKoRV and one phaCin-β were observed in coding sequence (CDS) regions. Of the 26 enKoRV integrations in CDS regions, one was classified as a start loss of a novel gene (gene ID g21247) located on contig JAOEJA010000073.1 at 15,262,207 bp. However, this was observed in only one of the SDZWA koalas (SB 127, born in 1985) and was not inherited by its offspring. The remaining 25 ERVs in CDS regions were classified as feature elongations, and only seven were observed in more than one individual, with the highest AAF=0.05. One ERV was located in the MCM9 gene, a gene encoding a mini-chromosome maintenance (MCM) protein that is indispensable for the initiation of eukaryotic genome replication. Mutations in this gene have been associated with ovarian failure, short stature, and chromosomal instability16.

ERV integrations over generations
The AAFs of deceased and living animals were different for some enKoRV and phaCin-β. In total, 714 ERVs were absent in the living population which suggests an elimination rate of 9.58×10−10 ERVs per base pair per generation for enKoRV and 2.51×10−10 ERVs per base pair per generation for phaCin-β. For example, the previously described enKoRV in the exon of the oncogene BCL2L1 (BCL2 like 1)10 was detected in four deceased animals, all heterozygous, but in none of the living animals. The four animals in three distinct triads died of leukemia. The integration was not transmitted further and disappeared from the pedigree (Fig. 3).

Most ERVs that did not persist in the current captive population were integrations restricted to single individuals that were not transmitted to offspring. However, 37 ERVs (25 enKoRVs and 12 phaCin-β) occurred with a higher frequency (AAF >0.05) but are also now absent from the current captive population (Supplementary Data 4). For example, an enKoRV located 8 kbp downstream of the ALDH1A3 (aldehyde dehydrogenase 1 family member A3) gene was present in seven individuals. Despite the fact that those koalas reached ages between 8 and 15 years, only two of them had four and five offspring, respectively, of which three were sequenced (Fig. 4). Of those three sequenced offspring, only one individual inherited the enKoRV, but did not transmit it to the next generations since it died of peritonitis, colitis, and metabolic imbalance five months after birth. Mutations in this gene are associated with fetal malformations in humans17, including anophthalmia/microphthalmia and facial dysmorphic features.

New ERV integrations
Across the 46 triads in the SDZWA, we detected 16 joey-specific ERV integrations, 15 enKoRVs and one phaCin-β, at an estimated 0.33 KoRV and 0.02 phaCin-β germline integrations per generation. The number of new ERVs varied from one (in three cases, SB 472, 560, and 748), two (in four cases, SB 551, 657, 665, and 674), to five (in one case, SB 659) per joey. Five of 16 new ERVs (20% probability) were detected in the subsequent generation. These included one new phaCin-β (SB 472), two new enKoRVs in the same joey (SB 674)—one of which was located in an intron of the DDAH1 (dimethylarginine dimethylaminohydrolase 1) gene—and two new enKoRVs detected in another joey (SB 659). The expected number of transmitted new ERVs was seven, given a 45% inheritance rate when considering a 10% probability of missing heterozygous ERVs based on a binomial model given a coverage ≥ ten sequencing reads. The observed transmission rate of five is lower than but still consistent with the accepted Mendelian inheritance rate.

Across the nine EUZ koala triads, we detected five joey-specific enKoRV germline integrations. The number of new ERVs varied from one (three cases, SB 395, 429, and 475) to two (one case, SB 369) integrations per joey, representing 0.56 germline integrations per generation. Only one new enKoRV was transmitted to the F2 generation (from SB 395 to SB 475). However, only two joeys had data available for the F2 generation and it is thus not possible to determine if the transmission is consistent with Mendelian inheritance. Two of the enKoRVs were located in the flanking regions of two genes, PTPN21 (protein tyrosine phosphatase non-receptor type 21), as part of the protein tyrosine phosphatases family, known to regulate a variety of cellular processes including cell growth, differentiation, mitotic cycle, and oncogenic transformation18, and MGAT5 (alpha-1,6-mannosylglycoprotein 6-beta-N-acetylglucosaminyltransferase) a mediator of cell migration in tumor development19,20.

enKoRV can dysregulate gene expression
Three genes containing commonly shared intronic enKoRV in the captive populations were differentially expressed in comparison to koalas without these enKoRVs. RNA-seq data was available for six of the EUZ koalas and five publicly available KoRV-free wild koalas10,21. Four of the EUZ individuals containing the most common intronic enKoRV in the SLC29A1 gene showed significantly higher SLC29A1 expression than individuals without the enKoRV (Fig. 5a). A full length KoRV-A was integrated in intron 9 of SLC29A1.

Two other genes upregulated due to the presence of an intronic enKoRV were LONP2 (Lon Peptidase 2) and STAP1 (Signal Transducing Adapter Family Member 1). LONP2 encodes for a peroxisome enzyme responsible for selectively degrading oxidatively damaged proteins, and is found to be upregulated in various cancers, in particular cervical cancer22. STAP1 is a gene that encodes for a protein involved in intracellular signaling by facilitating interactions between other proteins within signal transduction pathways. STAP1 has been linked to B-cell receptor signaling, which is crucial for the host immune response23. Moreover, STAP1 expression in glioma-associated microglia was positively correlated with the degree of malignancy and poor prognosis of glioma24. In SDZWA, 18% of the koalas carried the enKoRV in the LONP2 gene, while 60% of the EUZ koalas carried the same enKoRV. In the case of the enKoRV in the STAP1 gene, a prevalence of 11% in SDZWA and 25% in EUZ was observed. The enKoRVs were absent from wild populations examined.

enKoRV and phaCin-β are associated with increased cancer rates in koalas
Based on cause of death information for each koala, cancer rates were compared among deceased animals (n=58). In total, 27 individuals died of neoplasia, including diagnosis of lymphoma (n=13) and/or leukemia (n=11), sarcoma (n=6), osteochondroma (n=4), and adenocarcinoma (n=1). Overall, 34 enKoRVs showed an increased prevalence of more than 10% in cancer cases compared to other causes of death (Supplementary Data 5). Those 34 enKoRVs were located inside or in the flanking regions (±10 kb) of 15 genes, including known proto-oncogenes, such as BCL2L1, SLC29A1, and LZTS1 (leucine zipper tumor suppressor 1). In the case of phaCin-β, 14 showed a higher prevalence in neoplasia cases flanking three genes: DOCK8 (dedicator of cytokinesis 8), FAM216B (family with sequence similarity 216 member B), and Vmn2r65 (vomeronasal 2, receptor 65).

To estimate the direction and magnitude of ERV effects on neoplasia, a logistic regression was performed in a case-control model testing leukemia cases versus other death causes (associations with p≤0.05, Table 1). Seven alleles had an odds ratio (OR) >1 indicating a higher probability of developing leukemia. Interestingly, individuals without the enKoRV located in the intron of SLC29A1 gene, which has the highest enKoRV AAF in the koala population, showed 4.9-fold higher chances of developing leukemia than individuals carrying this enKoRV. The enKoRV in the exon region of BCL2L1, which is no longer present in the current captive population, showed a 21-fold higher chance of developing leukemia. The enKoRV in the intron region of the proto-oncogene LZTS1 had an OR of 9.3. Two phaCin-β integrations were potentially protective (OR <1), one located downstream of GIN1 (gypsy retrotransposon integrase 1) gene.

ERVs and SNPs are associated with fertility and longevity
Reproductive success (presence or absence of offspring) was tested using logistic regression. For longevity estimate, age at death was tested using a linear regression. In total, 24 deceased koalas never had offspring, while 34 had at least one. No SNP was significantly associated with the reproductive success trait, but 18 ERVs suggested association (associations with p≤0.05, Table 2). Two ERVs were associated with higher reproductive success, one enKoRV in the proximity of the genes RAB13 (RAB13, member RAS oncogene family), NUP210 (nucleoporin 210), and RPS27 (ribosomal protein S27), and one intergenic phaCin-β. The remaining ERVs showed decreased likelihood of reproductive success. The eliminated enKoRV downstream of ALDH1A3 also showed an OR <1, meaning lower chance of reproductive success. Another phaCin-β in the Vmn2r65 gene, close to the phaCin-β with 2.3-fold OR for leukemia in the intron of the same gene, also decreased the chances of reproductive success. However, it was less frequent in the population (AAF=0.05) than the phaCin-β associated with leukemia (AAF=0.44). Although not significant, the intronic enKoRV in the LONP2 (Lon Peptidase 2) gene, which caused an upregulation of the LONP2 gene, showed a p-value of 0.06 for having offspring with an OR=0.21 (Supplementary Data 6). Eight SNPs, 48 enKoRVs, and 13 phaCin-β were suggestively associated with age at death (Supplementary Fig. 9). Identified SNPs were in the introns of SPATA21 (spermatogenesis associated 21), NECAP2 (NECAP endocytosis associated 2), and CROCC (ciliary rootlet coiled-coil, rootletin) genes. Though not reaching statistical significance, the ERVs were found in the flanking region of genes, such as FANCI (FA complementation group I), associated with fanconi anemia25, and FAP (fibroblast activation protein alpha) and CX3CL1 (C-X3-C motif chemokine ligand 1), candidates for polyarteritis nodosa, a disease causing inflammation of the arteries26. CX3CL1 has also been linked to inflammatory diseases and cancer27.

Genetic risk score estimation
For two binary traits analyzed, leukemia and reproduction success, a genetic risk score was calculated based on the markers detected with a suggestive association. An average prediction accuracy of 0.64 and 0.80 was detected for reproduction success and leukemia, respectively, applying 100 bootstrapping repeats. The cutoffs (Fig. 6) were determined based on the highest accuracies generated.

The genetic risk score for living animals was also calculated. Based on the score, eight animals were identified as having high risk for developing leukemia. Similarly, four animals from the breeding program had low scores for reproductive success. An overlap of one individual was observed between the two traits analyzed. During the time this study took place, one individual (SB 451) determined as high risk died of leukemia.

In the case of age at death, since it is a continuous trait rather than a case-control variable, a simple breeding index was calculated, similarly to the genetic risk score (Supplementary Fig. 12). This breeding index for longevity should result in animals that could potentially live longer and remain healthier. This should be used together with breeding strategies to maintain genetic diversity, as a tool to support developing a healthier captive population.

Methods
Sample collection
In total, 186 koala individuals were included in this study: 91 from SDZWA, 20 from EUZ, and 75 wild koalas. Blood from living koalas from SDZWA were sampled in accordance with the Institutional Animal Care and Use Committee (IACUC) under IACUC proposals 15-017, 18-024 and 21-019 “Opportunistic Sample Collection during Veterinary Procedures”. Samples from deceased animals were taken opportunistically upon death of the animal and stored in the frozen tissues collection of the SDZWA. The 91 koalas from SDZWA (Supplementary Data 1) were shotgun sequenced with Illumina 150 bp paired-end reads at the Illumina laboratory Services. All 91 SDZWA individuals were represented by one sample except three individuals (SB 153, 169, and 308) which had duplicate samples. Sample types included blood (n=56), liver (n=23), lung (n=4), lymph nodes (n=4), nasal mass (n=3), spleen (n=3), and intestine (n=1). Samples were collected between 1992 to 2023 and preserved frozen at -80°C until DNA extraction. DNA was extracted using the Qiagen QIAamp DNA mini kit blood and tissue protocol. Samples with concentrations greater than 1 mg/µL were diluted using AE buffer. DNA quantitation was performed using a Qubit 4 fluorometer and Qubit dsDNA BR assay kit (Thermo Fisher Scientific, Waltham, MA, USA) and Quant-iT PicoGreen assay (Thermo Fisher Scientific, Waltham, MA, USA). Extracted DNA was kept at −20°C. Sample type and location are listed in Supplementary Data 1. Complete pedigree information was available for all SDZWA individuals, including cause of death for 58 deceased koalas. From the 58 deceased koala individuals, 27 were neoplasia related mortalities, including lymphoma, leukemia, and lymphosarcoma (n=20), sarcoma (n=3), osteochondroma and osteoma (n=3), and adenocarcinoma (n=1).

In addition, blood was taken from 20 live koalas in the European collections during regular veterinary care, according to each national regulations and institution code of ethic, or opportunistically upon the death of an individual koala. The collections were conducted upon the request of the coordinator (Kerstin Ternes) and the veterinary advisor (Baptiste Mulot) of the European Association of Zoos and Aquaria (EAZA) Ex Situ Program (EEP), in compliance with EAZA ethical standards. Samples were frozen immediately at -20°C, transported on dry ice and then stored at -80°C. The DNA was extracted from 100 µl of blood sample using the DNeasy Blood and Tissue Kit (Qiagen, Hilden, NRW, Germany) following the non-nucleated blood cell protocol. The DNA concentration and quality was then determined using an Agilent Tapestation (Agilent Technologies, Santa Clara, CA, USA) with Genomic Screen Tapes and reagents. Koalas from EUZ were sequenced using a PacBio based inverse PCR protocol as described previously10,36. Shotgun Illumina data for 75 koalas from wild populations in VIC (n=25), QLD (n=25), and NSW (n=25) was obtained from the open access Koala Genome Survey and scientific and ethics documentation from that project are listed within Hogg et al.37.

Paired-end library preparation and sequencing – SDZWA koalas
Paired-end (PE) libraries were generated for the 94 SDZWA samples (91 koalas, including three duplicates) using the TruSeq DNA PCR-Free Sample Prep kit (Illumina, San Diego, CA, USA). Prior to fragmentation, gDNA was purified using paramagnetic sample purification beads (Beckman Colter). DNA was then fragmented and libraries were end-repaired prior to size selection using TruSeq DNA PCR-Free beads for even coverage of areas that are traditionally difficult to sequence. The libraries yield was quantified using a real-time qPCR assay, and Illumina DNA Standards with primer master mix qPCR kit (KAPA Biosystems, Roche, Basel, Switzerland). Following library quantitation, 94 DNA libraries were denatured and normalized to 1.5 nM prior to be being loaded on a NovaSeq 6000 instrument where onboard clustering and sequencing occurs on 151 bp paired-end run, generating a minimum of 30X genome coverage per sample.

ERV Detection of PE reads – SDZWA and wild koalas
Reads were deduplicated, adapters removed, and filtered for an average quality score of 20 using fastp v0.22.038. Bases on the read tails were trimmed if their quality was below 20, retaining only reads with a minimum length of 35 bp. Subsequently, the processed reads were aligned to all KoRV strains (AB721500.1, KC779547.1, KP792564.1, KX587950.1, KU533853.1, KU533852.1, KX587961.1, KX587979.1, KX587966.1, and AB822553.1), phaCin-β, and phaCin-β-like 2 sequences using BWA-mem v2-2.2.139. Unmapped reads with a mate aligned to a viral sequence were retrieved using samtools v1.740. Moreover, soft-clipped regions (exceeding 20 bp) from reads only partially aligned to a viral sequence were extracted, preserving the ERV locations. Due to short length of PE reads, an accurate differentiation between KoRV strains was not possible. Therefore, all reads aligning to a KoRV strain were treated and named KoRV in further analyses; 98.2 and 95.6% of the reads aligning to KoRV in the koalas from SDZWA and in the wild koalas, respectively, aligned to KoRV-A, while the remaining reads aligned to KoRV-B. Clusters of reads containing both KoRV-A and KoRV-B were present in the dataset, which probably indicates incorrect strain assignment.

Both unmapped and soft-clipped reads were subsequently aligned to a retrovirus-masked koala assembly (accession number JAOEJA000000000)41 using BWA-mem v2-2.2.139, with alignments filtered to include only those with a minimum mapping quality score of 30. The masking was performed using RepeatMasker v4.1.2-p142, for the same viral sequences used in the first alignment step. Furthermore, ERV locations within the koala genome were identified by clustering at least 20 reads, requiring a maximum inter-cluster read-position difference of 500 bp for KoRV and 9 kbp for other viruses – since the JAOEJA000000000 assembly is enKoRV-free, but not free of other ERVs, which are normally around 8 kbp long. This approach suggests the proximity of an ERV based on the alignment between reads in the koala genome and their respective read-pair mates aligned to the retrovirus (Supplementary Fig. 11).

The exact ERV positions were determined based on the prevalence of overlaps among soft-clipped reads, distinguishing between novel and known integrations. Novel integrations were defined as those not identified in the JAOEJA000000000 assembly. This is the case for all KoRV integrations, as this assembly originates from an enKoRV-free southern Australian koala. The classification into novel or known integrations was based on the analysis conducted by RepeatMasker, with novel integrations pinpointing to a single exact position. In contrast, known integrations previously identified in the JAOEJA000000000 assembly were characterized by both start and end positions.

In total, before filtering, 9178 KoRV, 9040 phaCin-β, 83601 phaCin-β-like, were detected for all 94 koala samples. ERVs were considered shared between individuals when their cluster locations overlapped, based on the start and end positions. A final list of shared ERVs was used to call zygosity.

ERV zygosity
The zygosity of the integrations was estimated after initially aligning all trimmed reads to the unmasked koala JAOEJA000000000 assembly. This estimation was based on comparing the distribution of soft-clipped reads against reads spanning the ERV positions. Reads at these positions were deemed soft-clipped if they were truncated by BWA and terminated precisely at, or within 10 bp of, the ERV positions—accounting for minor mismatches due to the Target Site Duplication (TSD) resulting from the proviral integration mechanism. Conversely, reads were classified as crossing if they spanned the integration sites, with start and end points situated beyond the TSD margins (Supplementary Fig. 12).

Zygosity was assessed differently for novel insertions compared to those previously observed in the JAOEJA000000000 assembly. For novel variants, a ratio threshold of 30% of the lesser-read group count was applied to classify insertions as heterozygous. An insertion was deemed homozygous if over 70% of the pertinent reads were soft-clipped; otherwise, an ERV with less than 30% soft-clipped reads was categorized as having undetermined zygosity, possibly indicating false positives or somatic insertions. For known insertions, both start and end positions were used for classification, with homozygosity determined by over 70% crossing reads.

At the population level, any integration not detected was assumed homozygous reference (indicating the absence of an integration). For this study, the viral integration represents the alternative allele, while its absence signifies the reference allele. In addition, the Mendelian error rate was calculated for all available triads or pair parent-offspring. In total, 63 koalas had at least one parent genotyped. AAF, percentage of undetermined zygosity, and Mendelian error per ERV were used for filtering. Only ERVs with an AAF>0, percentage of undetermined zygosity ≤0.1, and Mendelian error rate per ERV ≤0.05 were kept for further analyses. An average Mendelian error rate of 0.1±0.77% was observed per animal after filtering. In addition, genotype concordance was also checked for the three replicates, ranging from 96.0 to 98.0%. For subsequent analyses, a consensus between the replicates of each of the three koalas was used, being inconsistent genotypes coded as undetermined (coded as -1).

SNP calling of PE data
Reads filtered for quality control as previously described in the ERV detection step were aligned to the koala assembly JAOEJA000000000 using dragmap v1.3.043, and sorted and merged with samtools v1.740. SNPs and indels were called using GATK v4.4.044 per sample on HaplotypeCaller tool and then recalled population wide on GenomicsDBImport and GenotypeGVCFs tools. Finally, SNPs and indels were filtered with the VariantFiltration tool of GATK using several criteria. Variants with an RMS Mapping Quality (MQ) <40.0, a Fisher Strand Bias (FS) >60.0, a Quality by Depth (QD) <2.0, a Strand Odds Ratio (SOR) >3.0, a Mapping Quality Rank Sum Test (MQRankSum) <-12.5, a Read Position Rank Sum Test (ReadPosRankSum) <-8.0, and a Depth (DP) <10 were removed from the analysis. In addition, only variants with a call rate >90%, totalizing 21,083,491 SNPs and indels were used for further analysis.

Furthermore, SNPs for enKoRV were called similarly to the host SNPs, using the mapped reads to KoRV-A reference sequence (AB721500). Filtering of reads followed the same thresholds, totaling 111 SNPs.

iPCR and PacBio sequencing—EUZ koalas
DNA was sheared and circularized10, for the 20 samples from the EUZ. Sonication inverse PCR (SIP) reactions were then performed on all samples (fragmented and circularized) using the following LTR primers: forward 5’-ATTTGCATCCGGAGTTGGT-3’ and reverse 5’-AGGGGCACCCTAGAAACTGT-3’. Libraries were then created for each sample using the PacBio (Pacific Biosciences, Menlo Park, CA, USA) 5 kb template preparation protocol and the SMRTbell™ Template Prep Kit 1.0 following the manufacturer’s guidelines. Sequencing on the PacBio Sequel II platform was performed using the MagBead Standard protocol, C4 chemistry and P6 polymerase on a single v3 Single-Molecule Real-Time (SMRT) cell with 1×180 min movie for each sample. Every sample was run across two Single-Molecule Real-Time (SMRT) cells to increase the coverage. The reads from the insert sequence were processed within the SMRT® Portal browser (minimum full pass = 1; and a Minimum Predicted Accuracy of 90).

enKoRV detection of long-reads – EUZ koalas
KoRV integrations were detected as previously described10. Briefly, parts of the PacBio reads mapping to any KoRV strains were identified using BLAST v2.12.045 and removed so that only the koala genomic sequence (flanking regions of the provirus) remained. The trimmed reads were aligned to the South Australian koala genome (GCA_030178435.1) using Minimap v2.146 without secondary alignments. Matches longer than 50 bp, with ≥90% of matching bases, and with a mapping quality ≥ 30 were retained. KoRV integrations were detected by finding peaks of alignment coverage >10 reads per breakpoint with an overlap within 10 bp (as of the TSD) between the breakpoints, except for cases of masked sequences, where the maximum masked sequenced for recKoRV of 1,411 bp was allowed between breakpoints. ERV and somatic KoRV were defined based on coverage, with a cutoff of 500 reads established as the minimum threshold to define enKoRV. An exception was made for one koala (SB 369), where the threshold was set at 1000 reads (see Supplementary Fig. 13).

RNA-Seq and visualization of ERVs
RNA-Seq data from five publicly available KoRV-free koala individuals was used (PRJEB21505) in addition to RNA-Seq from six of the EUZ koalas sequenced with PacBio. RNA-Seq data was processed following the open-source snakemake pipeline47, trimming reads with fastqc v0.12.148, aligning and sorting reads with hisat2 v2.2.149 and samtools v1.1740, respectively, and counting reads with htseq-count v2.0.850. Gene count was normalized using the DESeq2 v4.1.3 package in R v4.1.1. For each of the genes found with an ERV in our study and present in the RNA-Seq two groups were formed, one with the ERV and the other without it. The two groups were compared with a two-side Mann-Whitney test using the SciPy v1.11.4 package51 in Python v. 3.12. ERVs were visualized using either the reference genome sequence, or alignment of assembled reads in the Geneious Prime v2025.0.2 software (www.geneious.com).

PCA
Principal component analysis (PCA) of both SNP and ERV data were performed using the Scikit-learn v. 1.552 package in Python v. 3.12. In the case of enKoRV data alone, data was first transformed using standard scale. In addition, only non-missing markers with an AAF≥0.05 were used.

Detection of joey specific ERVs
New ERV integrations were assessed in 55 triads, 46 triads from the SDZWA population and nine from the EUZ population. In the case of SDZWA triads, new ERV integrations were represented by heterozygous high-coverage ERVs unique to the joey and absent in both parents (coverage ≥40 reads, corresponding to a 0.6% probability of missing heterozygous ERVs based on a binomial model for a 30–70% heterozygous ratio). New ERVs with lower coverage were very likely somatic integrations detected. Since 19 F2 were available for those triads, new ERVs could be validated in the next generation. For these, the presence of one copy of the ERV in the offspring with a coverage ≥ 10 reads was used as validation. In the case of EUZ, only KoRV integrations were available, and at much higher coverage, given the inverse PCR methodology used15. The same approach used for the detection of new ERVs in the SDZWA was used for EUZ koalas, except for the zygosity estimation and the definition of somatic integrations. In order to remove somatic integrations enriched, a coverage threshold at ≥500 reads for both breakpoints of an integration was established, keeping only integrations with a coverage above this threshold. The coverage for the joey SB 369 was set to 1000 reads due to higher coverage rate (manually curated based on the coverage distribution of unique versus shared ERVs – Supplementary Fig. 13).

Mutation rate and initial germline colonization estimates
The germline mutation rate in koala was estimated based on the 46 triads from the SDZWA. High-quality unique SNPs not detected in either parent were considered a de novo mutation, using GATK v4.444 (FindMendelianViolations). A more stringent quality control for de novo SNPs was applied based on recommendations of the results from the Mutationathon benchmark53. Those include a minimum depth of 20 reads, minimum genotyping quality of 20, and minimum allele balance of 0.4 at sites that are heterozygous in the offspring. Mutation rate per base pair per generation (mutations_site_generation) was calculated as the total number of de novo SNPs divided by the number of triads multiplied by the diploid assembly length (6,469,964,576 bp). Mutation rate per base pair per year (mutations_site_year) was calculated in the same way, except that the total number of SNPs was estimated per year, summing up the number of detected SNPs per year per joey divided by the arithmetic mean age of the parents.

The initial germline colonization times of KoRV, phaCin-β and phaCin-β-like were estimated based on the 5’-LTR and 3’-LTR divergence within the same provirus. Since the 5’-LTR and 3’ LTR are identical at the time of integration, any variations accumulated over time work as a molecular clock that can be used to estimate earliest time of integration. The enKoRV located in the intron of the SLC29A1 gene was used to determine the LTR divergence, since it is the enKoRV with the highest frequency in the koala population and very likely one of the most ancient, which shows divergence (p-distance) of 0.001 (one insertion/deletion observed in 1000 bases-Supplementary Fig. 4). In the same way, germline colonization times of phaCin-β and phaCin-β-like were updated based on the median LTRs divergence previously reported2, following the formula13:

Colonization time = LTR divergence / (mutations site year * 2)

Genome-wide association
Using a logistic regression in PLINK v2, SNPs and ERVs were analyzed in a case-control model for leukemia, leukemia strict control (without other cases of neoplasia as controls), neoplasia, neoplasia without leukemia, and neoplasia without leukemia strict control, and fertility (in terms of total offspring: 24 koalas with no offspring versus 34 with offspring ranging from 1 to 14). In addition, ERVs and SNPs were tested for longevity determined as the age at death using a linear regression in PLINK v2. Markers with Minor Allele Frequency (MAF) ≥0.01 were tested for the additive effects. AIC (Akaike information criterion) was used to test which available cofactors including age at death, total offspring, sex, death location, birth location, cause of death, and PC1–PC3 were contributing to each regression model (Delta AIC <-5). The cofactor total offspring was used in the longevity model, while the cofactor age was used in the fertility model. For the other investigated phenotypes, no cofactor contributed to the models.

For ERVs and SNPs, thresholds for significant (p<0.05) associations were estimated using Bonferroni multiple testing correction considering 953,591 independent SNPs derived from linkage disequilibrium based pruning in PLINK v2 with the parameters window size=500, step size=100, and r2=0.5. A log base 10 threshold of 1.3 (p<0.05 without multiple testing correction) was used to detect suggestive ERVs and of 6 for suggestive SNPs. Genes were annotated for associated SNPs and ERVs, and their flanking region of ±10 kb, using previously predicted gene annotation41. Scaffolds were assigned to chromosomes as defined in Supplementary Note 1. ERV and SNP consequences were predicted using Ensembl Variant Effect Predictor release 10554. Association to viral SNPs were also tested as described in Supplementary Note 2, Supplementary Fig. 15, and Supplementary Data 9.

Genetic risk score
Genetic risk score (GRS) and the longevity breeding index were estimated to assess genetic contributions to traits. The GRS was computed using odds ratio (OR) values obtained from suggestive markers associated with leukemia and fertility, where each individual’s score was derived by summing the products of these OR coefficients and their respective genotype counts. The longevity breeding index, since it is a continuous trait, similarly utilized beta values.

To estimate prediction accuracy of the GRS estimates, we employed a bootstrapping approach with 100 iterations, randomly splitting the dataset into training (75%) and testing (25%) sets. A logistic regression model was trained on the GRS to predict outcomes. We calculated accuracy as the proportion of correctly classified cases and averaged the results from all iterations, providing a robust estimate of prediction accuracy. Different cutoffs of the GRS were tested, and the cutoff which generated the highest accuracy was selected. The GRS calculated for living animals was used then to detect animals falling into the cutoff estimates, as individuals one would recommend to be avoided in SDZWA koala breeding program.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This paper reports a detailed analysis of the endogenous retroviruses found in koalas. It is novel and interesting but very\nhard going because of the vast amount of data reported. I can’t help wondering whether it might have been better to focus on\nthe enKoRVs, leaving a comparison with the other elements for a follow up paper. Certainly, it would be easier on the reader\nand I doubt whether the main conclusions would change. In any event it is an impressive body of work.\nSome specific comments\nL57. Insert ‘per base pair per year’ after ‘of humans’\nL93-5. Relatedness to what?\nL116. (ST1) Why are entries to the Studbook not listed in numerical order?\nL127. The number of phaCin-β-like elements is much higher than previous estimates (L71-6. Some comment would appear\nindicated.\nL162. If the phaCin-β and phaCin-β-like colonization dates increase by c25% why does that of the enKoRV elements\nincrease several fold? Is there overreliance on the SLC29A1 provirus?\nL180. Are there really numerous segregating phaCin-β-like elements?\nL217. Isn’t it surprising that neither parent of SB153 shows this relationship?\nL222. There seems to have been some confusion in quoting ref 14 (“Given the historical population bottleneck of Narrandera\nkoalas in the 1890s and their subsequent reintroduction from Queensland and Victoria in the 1970s, we designated….”)\nL269. Is there no information on the cause of death of SB21, 24, 111, 196?\nL278. Or seven? Does SB268 not count?\nL273-284. I wonder what conclusion on is supposed to draw from this section given that the provirus clearly has been\ntransmitted to offspring? Is there any evidence for fetal malformations or other harm in koalas?\nL484. Given that these animals have a decent lifespan how would such negative selection work? Reduced fertility?\nL531. One might argue that the first priority should be to try to identify individuals that show no evidence of reinsertions.\nL541. How similar was the proviral content of duplicate samples from the same individual?", "reviewer_2": "In this manuscript, analysis of the endogenous retroviruses (ERVs) in the Koala was carried out. Virtually all higher\norganisms carry multiple copies of ERVs that result from germline infection (and insertion of viral DNA), followed by\ntransmission of the ERVs to subsequent generations. The majority of ERVs result from ancient infections in progenitors to\nthe modern day species, and the great majority of ERVs have become fixed in the genomes. The retroviruses that gave rise\nto the ERVs are generally no longer present as infectious agents for most species. The koala is a unique case in that it\ncarries two infectious retroviruses, KoRV and phaCin-beta that are currently generating new ERVs. This is evident from the\nfacts that many ERVs of these two classes are not fixed in the koala genome, and that new ERVs (particularly endogenous\nKoRVs) can be detected in offspring. Study of the ERVs in koalas provides a window into early phases of ERV\nendogenization which cannot be studied in other species (including humans).\nThis report continues this group’s studies of koala retroviruses and ERVs, focusing on koala colonies from European zoos,\nand from the San Diego Zoo Alliance. They combine a broad range of techniques, including genomic sequencing and\nRNAseq, with family pedigrees (over as many as 4 generations) and health records of 91 animals. They have identified by\nDNA sequencing more than 2000 ERVs in the population. The great majority of the endogous KoRV and phaCin-beta ERVs\nwere not fixed in the genome, while an older phaCin-beta-like ERV was largely fixed. The investigators applied population\ngenetic analyses which provided insight into the impacts of new koala ERVs on risk for\nneoplasia and reproduction.\nThis study is a powerful combination of genetic and molecular studies with genealogy and health records in two well-studied\npopulations that provides important new insights into the effects of ERVs as they invade a species. Implications are clear for\nour general understanding of ERVs of other species such as humans where ERVs make up ~8% of the genome. The first\nstudies in the paper recalibrating when KoRV and phaCin-beta entered the germlines of koalas or their ancestors, as well\ncharacterization of the distribution of ERVs in the zoo populations compared to wild koalas in northern and southern\nAustralia are convincing.\n1. In the section on ERVs through the generations and the examples in Figs 3 and 4, the implication is that the ERVs are\nresponsible for the failure of the progeny to survive or breed. One possibllity could be that these ERVs are producing\ninfectious KoRV which could lead to the neoplasia or loss of fertility. Do animals with these inherited ERVs have higher\nKoRV viral loads? Perhaps future RNAseq or qRT-PCR assays could shed light on this. For the pedigree in Fig 3, if tumor\ntissues are available, were there somatic KoRV integrations near other proto-oncogenes that might collaborate with the\nBCL2L1 insertion? Also, since the ERVs in both Fig 3 and Fig 4 were originally present in more than one animal, a block to\ngenetic transmission was obviously not absolute.\n2. One of the most interesting results from the study are the quite high rates of new endogenizations observed over just one\ngeneration (i.e. the triads analyzed). However, the results for the SDZWA and EUZ analyzed should be reported similarly.\nFor EUZ, 4 of the 8 triads showed new ERVs; for the SDZWA animals the results are reported in terms of the number of new\nERVs (16) in the 46 triads, but it appears that these new ERVs were confined to 8 of the 46 triads. In any event, the high\nfrequency (and frequent instances of more than one new ERV in an F1) bears comment in the Discussion. Is there evidence\nof high KoRV or phaCin-Beta expression or infection in tissues of the female or male reproductive tract?\n3. Another strength of this work is that the effects of koala ERVs (in aggregate as well as individually) on cancer rates, fertility\nand longevity were determined. As mentioned above, it will be important to distinguish between high levels of infectious\nKoRV (or phaCin-beta) vs. effects of the proviral insertion on the host gene. The effects of a given ERV may reflect either or\nboth of these processes. The ERVs associated with altered cancer incidence in Table 1 and those associated with\nreproductive success in Table 2 are interesting. The written narrative (351-413) is impeded by description of the functions of\nmany of the host genes at the insertion sites. It would be clearer to focus on a few ERVs where there is a strong case for\nbiological importance, e.g. the insertions in SLC29A1 and BCL2L1. The other host sites could be grouped into proto-\noncoges, tumor suppressor genes, etc. and simply listed by name. For instance the description of the function of LZTS1 on\nlines 358-360 could be eliminated.\n4. Finally, it would be helpful for the authors to set their studies into a larger biological/evolutionary framework in the\nDiscussion. On the one hand their results demonstrate significant effects of ERVs (particularly KoRVs) on cancer and\nreproductive capacity on koalas, as well as high ongoing endogenization of new ERVs. On the other hand, they now\ncalculate that KoRV (the most recent koala ERV) may have begun to infect at least some koalas as far back as 300,000\nyears ago, concomitant with emergence of modern koalas. And KoRV as an infectious agent may be spreading from north to\nsouth in Australia. What are the Implications for the wild koala population?", "reviewer_3": "Neumann and team have presented a manuscript that significantly advances our understanding of key retroviruses in the\nkoala. Currently the koala is unique in that it is experiencing germline invasion by two viruses, KoRV and phaCin-B. The\nstudy of the process of invasion and inheritance by these viruses provides key understandings to not only the koala but\nmuch more broadly. The authors have significantly advanced previously published work by using a powerful set of 111\nkoala samples from family trees of koalas with up to three known generations. This enabled them to follow inheritance of key\ninsertion sites in the koala genome. The use of these multi-generational pedigrees combined with key health data on both\nliving and deceased koalas enabled the authors to link key insertional sites with health predictors. Finally, the key retroviral\ninsertions could be associated with health signals such as leukemia, fertility and longevity, enabling the authors to develop a\ngenetic risk score, which will be very useful for breeding of koalas in captivity in particular.\nThe authors have presented a significant amount of raw data, both in the main manuscript and also in Supplementary files,\nand this data strongly supports their study and their conclusions. The science is high quality and the methodology is strong\nand very well structured.\nThis study is the first comprehensive analysis of such a large and well structured pedigree to further characterise koala\nretrovirus integrations and their positive and negative selections.\nOne minor error on Line 445. The Lone Pine Koala Sanctuary is located in Brisbane, not on the Gold Coast of Queensland\nAustralia.\nI strongly support the publication of this manuscript."}, "author_response": "RESPONSE TO REVIEWERS’ COMMENTS\nReviewer #1 (Remarks to the Author):\nThis paper reports a detailed analysis of the endogenous retroviruses found in koalas. It\nis novel and interesting but very hard going because of the vast amount of data reported.\nI can’t help wondering whether it might have been better to focus on the enKoRVs,\nleaving a comparison with the other elements for a follow up paper. Certainly, it would\nbe easier on the reader and I doubt whether the main conclusions would change. In any\nevent it is an impressive body of work.\nResponse: We thank the reviewer for their thorough and positive assessment of the\nwork. We do understand the concern about the amount of data, but the authors decided\nto report it all together as it contributes to the same story line and strongly support the\nresults. In particular, the genetic risk scores can be used to implement changes in\ncaptive koala (and perhaps at some point managed wild koala) breeding programs.\nTherefore, it is important they have all available data to support their decisions.\nSome specific comments\nL57. Insert ‘per base pair per year’ after ‘of humans’-\nResponse: Done. Placed before humans, since it is the case for both humans and mice\nL93-5. Relatedness to what? –\nResponse: we changed to “Comparisons with wild koalas indicated that the SDZWA\nkoalas are genetically more closely related to northern Australian koalas. In addition,\nthere is evident genetic relatedness between SDZWA koalas and koalas from European\nzoos as expected due to translocations between zoos in North America and Europe.”\nL116. (ST1) Why are entries to the Studbook not listed in numerical order? –\nResponse: We thank the reviewer for noticing this error in presentation. We changed to\nnumerical order.\nL127. The number of phaCin-β-like elements is much higher than previous estimates\n(L71-6. Some comment would appear indicated.-\nResponse: The numbers reported previously were only for polymorphic phaCin-β-like,\nwhile in our study we report all integrations, including fixed phaCin-β-like. This has to do\nwith the difference in methodology. Since phaCin-β-like are fixed and present in the\nkoala reference genome, Lillie et al. only detected deletions, meaning polymorphic\nreference ERV loci. So that means that our numbers are much smaller. Only six phaCin-\nβ-like were segregating in SDZWA and 14 when we also consider the 75 wild koalas,\nout of 361 phaCin-β-like in total. If we would use all 430 koalas, as previously, the numbers\nwould probably increase, out of 469 phaCin-β-like previously detected in total (Lillie et\nal, 2024). We added to the discussion now the following sentence: L465 “The difference\nin the number of phaCin-β-like insertions per individual between the current study and\nprevious report (Lillie et al., 2024) is due to methodological differences. Previously, only\npolymorphic reference ERVs were examined. Since phaCin-β-like insertions are\ngenerally fixed and present in the reference koala genome, most phaCin-β-like\ninsertions were not reported by Lillie et al.”.\nL162. If the phaCin-β and phaCin-β-like colonization dates increase by c25% why does\nthat of the enKoRV elements increase several fold? Is there overreliance on the\nSLC29A1 provirus?\nResponse: The difference in the increase for enKoRV is due to a methodological\ndifference in the calculation. We opted to keep the same approach applied in the study\nof phaCin-β and phaCin-β-like (Lillie et al., 2024), as referred to in the Methods.\nL180. Are there really numerous segregating phaCin-β-like elements? -\nResponse: We thank the reviewer for the observation. Numerous was not an\nappropriate term. We decided to remove this statement, keeping only the information\nabout the AAFs of phaCin-β: “However, in comparison to SDZWA, the AAFs of phaCin-\nβ in VIC were higher (Supplementary Figure 7).”.\nL217. Isn’t it surprising that neither parent of SB153 shows this relationship? –\nResponse: It is indeed a bit surprising. This was likely caused by the lower sequence\ncoverage of the liver sample. In this case, we had duplicates: one sample from lung\ntissue and one from liver tissue. Since we decided to use the consensus genotype\ninstead of the genotype from the sample with the highest coverage, the number of ERVs\nfor this individual was reduced to 417, compared to an average of 470 ERVs per koala in\nthe SDZWA population. Nevertheless, the low Mendelian error rates of 0.12% for the\nlung and 0.04% for the liver validate the relationship within the triad. Additionally, based\non the called SNPs, this individual is correctly placed within the SDZWA clade, in direct\nassociation with the parents (not reported in the manuscript).\nL222. There seems to have been some confusion in quoting ref 14 (“Given the historical\npopulation bottleneck of Narrandera koalas in the 1890s and their subsequent\nreintroduction from Queensland and Victoria in the 1970s, we designated….”)\nResponse: We thank the reviewer for catching this mistake. We changed the time range\nfrom the 1890s to 1970s.\nL269. Is there no information on the cause of death of SB21, 24, 111, 196?\nResponse: Unfortunately not for all of them. SB21 and SB24 died in Australia. SB111\nwas apparently healthy, and had anesthetic death (no histology report available). SB196\nhad a spinal fracture. We color only the death by leukemia. We change the figure\ncaption to “…Death or euthanasia due to leukemia (where know) is indicated by color\npattern”.\nL278. Or seven? Does SB268 not count?\nResponse: True, it is indeed seven. We did not count SB268, as it was the only observed\noffspring that inherited the ERV, but considering the way sentence is written, we should\nalso count it.\nL273-284. I wonder what conclusion on is supposed to draw from this section given that\nthe provirus clearly has been transmitted to offspring?\nResponse: There is no clear evidence for fetal malformations, although, it is known in\nkoalas. We believe this ERV could be associated with decreased fertility, in other words,\nan increased predisposition to abortion, for instance. That is not necessarily a\nmendelian trait, rather a complex one. Based on the graph, we could infer that SB108,\nSB210, and SB240 probably carried the same ERV. None of them seemed to have\nstruggled to have offspring, having 4, 5, and 6, respectively. Nevertheless, we believe, it\nis worth the observation, since, other people working with koalas could eventually\ndetect the same ERV.\nL484. Given that these animals have a decent lifespan how would such negative\nselection work? Reduced fertility?\nResponse: The idea is that the carriers of this ERV die younger, reducing the potential\noffspring number, and reducing the chances of transmitting the ERV. The average age at\ndeath of our cohort was 9 years old. From the four carries, two females died at ages 5\nand 7 years old, with total offspring 3 and 5 joeys each. While, the males lived both 10\nyears, and had an offspring of 3 and 9 joeys each.\nL531. One might argue that the first priority should be to try to identify individuals that\nshow no evidence of reinsertions.\nResponse: The problem is that the reinsertions only show up during the developing of\nneoplasia, as previously observed comparing tumor tissues versus health tissues\n(Gayle et al, 2021). When the tumor is observed, it is already too late and the animal is\neuthanized.\nL541. How similar was the proviral content of duplicate samples from the same\nindividual?\nResponse: Since no enrichment was performed, most proviruses detected were ERVs.\nSome somatic viral insertions detected were probably discarded due to low coverage\nduring zygosity calling. The consensus between the duplicates varied from 93% to 96%\nfor raw data, and 96% to 98% after filtering (see Methods, L673).\nReviewer #2 (Remarks to the Author):\nIn this manuscript, analysis of the endogenous retroviruses (ERVs) in the Koala was\ncarried out. Virtually all higher organisms carry multiple copies of ERVs that result from\ngermline infection (and insertion of viral DNA), followed by transmission of the ERVs to\nsubsequent generations. The majority of ERVs result from ancient infections in\nprogenitors to the modern day species, and the great majority of ERVs have become\nfixed in the genomes. The retroviruses that gave rise to the ERVs are generally no longer\npresent as infectious agents for most species. The koala is a unique case in that it\ncarries two infectious retroviruses, KoRV and phaCin-beta that are currently generating\nnew ERVs. This is evident from the facts that many ERVs of these two classes are not\nfixed in the koala genome, and that new ERVs (particularly endogenous KoRVs) can be\ndetected in offspring. Study of the ERVs in koalas provides a window into early phases of\nERV endogenization which cannot be studied in other species (including humans).\nThis report continues this group’s studies of koala retroviruses and ERVs, focusing on\nkoala colonies from European zoos, and from the San Diego Zoo Alliance. They combine\na broad range of techniques, including genomic sequencing and RNAseq, with family\npedigrees (over as many as 4 generations) and health records of 91 animals. They have\nidentified by DNA sequencing more than 2000 ERVs in the population. The great majority\nof the endogous KoRV and phaCin-beta ERVs were not fixed in the genome, while an\nolder phaCin-beta-like ERV was largely fixed. The investigators applied population\ngenetic analyses which provided insight into the impacts of new koala ERVs on risk for\nneoplasia and reproduction.\nThis study is a powerful combination of genetic and molecular studies with genealogy\nand health records in two well-studied populations that provides important new insights\ninto the effects of ERVs as they invade a species. Implications are clear for our general\nunderstanding of ERVs of other species such as humans where ERVs make up ~8% of\nthe genome. The first studies in the paper recalibrating when KoRV and phaCin-beta\nentered the germlines of koalas or their ancestors, as well characterization of the\ndistribution of ERVs in the zoo populations compared to wild koalas in northern and\nsouthern Australia are convincing.\nWe thank the reviewer for their positive assessment and for providing valuable insights\nregarding its potential outcomes.\n1. In the section on ERVs through the generations and the examples in Figs 3 and 4, the\nimplication is that the ERVs are responsible for the failure of the progeny to survive\nor breed. One possibllity could be that these ERVs are producing infectious KoRV\nwhich could lead to the neoplasia or loss of fertility. Do animals with these inherited\nERVs have higher KoRV viral loads? Perhaps future RNAseq or qRT-PCR assays could\nshed light on this. For the pedigree in Fig 3, if tumor tissues are available, were there\nsomatic KoRV integrations near other proto-oncogenes that might collaborate with\nthe BCL2L1 insertion? Also, since the ERVs in both Fig 3 and Fig 4 were originally\npresent in more than one animal, a block to genetic transmission was obviously not\nabsolute.\nResponse: This is a very interesting idea. Unfortunately, we do not have any data on viral\nloads for this population. However, we agree that this information could be very\ninformative, and we plan to conduct more RNA sequencing in the future, as tumor\ntissues from these animals were not available. Nevertheless, there is clear support from\nour previous work (Gayle et al., 2021) showing reintegrations near other proto-\noncogenes in tumor tissues. Checking whether the same regions are enriched in tumor\ntissues of koalas from SDZWA would be very interesting.\n2. One of the most interesting results from the study are the quite high rates of new\nendogenizations observed over just one generation (i.e. the triads analyzed).\nHowever, the results for the SDZWA and EUZ analyzed should be reported similarly.\nFor EUZ, 4 of the 8 triads showed new ERVs; for the SDZWA animals the results are\nreported in terms of the number of new ERVs (16) in the 46 triads, but it appears that\nthese new ERVs were confined to 8 of the 46 triads. In any event, the high frequency\n(and frequent instances of more than one new ERV in an F1) bears comment in the\nDiscussion. Is there evidence of high KoRV or phaCin-Beta expression or infection in\ntissues of the female or male reproductive tract?\nResponse: We thank the reviewer for the comment. We changed the results to report\nboth results similarly. Indeed, the detection of new ERVs was unequal in both\npopulations. We speculate this might be due to differences in the methods, as data on\nenKoRV for EUZ was enriched, causing difficult in the differentiation between somatic\nand endogenous insertions. Regarding the high rates, we now added to Discussion “The\nhigh rates of new ERV integrations are consistent with previous observations of an\nactive ERV element in cattle (Tang et al., 2024). Although the transmission rate was not\nreported, the ERVK[2-1-LTR] element in cattle exhibited an average mobilization rate of\n1 de novo ERV per approximately 150 sperm cells, with more than a 10-fold difference in\nrates between individual animals. The de novo ERVK[2-1-LTR] elements tended to\npreferentially insert near telomeric ends, in GC-rich regions, and within genes.”.\nIn the future, one could attempt to analyze koala sperm to quantify the rates of new\ninsertional ERV polymorphisms in the male germline.\n3. Another strength of this work is that the effects of koala ERVs (in aggregate as well as\nindividually) on cancer rates, fertility and longevity were determined. As mentioned\nabove, it will be important to distinguish between high levels of infectious KoRV (or\nphaCin-beta) vs. effects of the proviral insertion on the host gene. The effects of a given\nERV may reflect either or both of these processes. The ERVs associated with altered\ncancer incidence in Table 1 and those associated with reproductive success in Table 2\nare interesting. The written narrative (351-413) is impeded by description of the\nfunctions of many of the host genes at the insertion sites. It would be clearer to focus on\na few ERVs where there is a strong case for biological importance, e.g. the insertions in\nSLC29A1 and BCL2L1. The other host sites could be grouped into proto-oncoges, tumor\nsuppressor genes, etc. and simply listed by name. For instance the description of the\nfunction of LZTS1 on lines 358-360 could be eliminated.\nResponse: We thank the reviewer for this suggestion. We changed results accordingly.\n4. Finally, it would be helpful for the authors to set their studies into a larger\nbiological/evolutionary framework in the Discussion. On the one hand their\nresults demonstrate significant effects of ERVs (particularly KoRVs) on cancer\nand reproductive capacity on koalas, as well as high ongoing endogenization of\nnew ERVs. On the other hand, they now calculate that KoRV (the most recent\nkoala ERV) may have begun to infect at least some koalas as far back as 300,000\nyears ago, concomitant with emergence of modern koalas. And KoRV as an\ninfectious agent may be spreading from north to south in Australia. What are the\nImplications for the wild koala population?\nResponse: We changed discussion accordingly adding a paragraph on implication for\nwild population: “Considering that the process of endogenization has been occurring in\nkoalas for at least the last 300,000 years, and that this is a natural event responsible for\nshaping 5 to 10% of mammalian genomes, it does not immediately warrant intervention\nin the wild. However, when other threats—such as climate change and habitat loss\ncaused by human activity—are added to the equation, the situation worsens. Therefore,\nmonitoring is recommended to allow genetic drift and natural selection to proceed.\nFollowing this, further evaluation of whether intervention is possible or practical should\nbe conducted, especially for critically endangered populations. Nonetheless, for\ncaptive populations that serve as gene reserves for the species, it becomes crucial to\nprevent an increase in the frequency of deleterious ERVs.”\nReviewer #3 (Remarks to the Author):\nNeumann and team have presented a manuscript that significantly advances our\nunderstanding of key retroviruses in the koala. Currently the koala is unique in that it is\nexperiencing germline invasion by two viruses, KoRV and phaCin-B. The study of the\nprocess of invasion and inheritance by these viruses provides key understandings to not\nonly the koala but much more broadly. The authors have significantly advanced\npreviously published work by using a powerful set of 111 koala samples from family\ntrees of koalas with up to three known generations. This enabled them to follow\ninheritance of key insertion sites in the koala genome. The use of these multi-\ngenerational pedigrees combined with key health data on both living and deceased\nkoalas enabled the authors to link key insertional sites with health predictors. Finally,\nthe key retroviral insertions could be associated with health signals such as leukemia,\nfertility and longevity, enabling the authors to develop a genetic risk score, which will be\nvery useful for breeding of koalas in captivity in particular.\nThe authors have presented a significant amount of raw data, both in the main\nmanuscript and also in Supplementary files, and this data strongly supports their study\nand their conclusions. The science is high quality and the methodology is strong\nand very well structured.\nThis study is the first comprehensive analysis of such a large and well structured\npedigree to further characterise koala retrovirus integrations and their positive and\nnegative selections.\nOne minor error on Line 445. The Lone Pine Koala Sanctuary is located in Brisbane, not\non the Gold Coast of Queensland Australia.\nResponse: Thank you for the comment. We corrected that to Brisbane.\nI strongly support the publication of this manuscript.\nWe thank the reviewer for taking the time to review our manuscript and for their support\nfor publication."}, {"round": 2, "reviewer_comments": {"reviewer_1": "I am happy to support publication of this manuscript in revised form although I still feel that a degree\nof simplification by focussing on the KoRVs would be advantageous. However, I do fully understand\nthe authors’ wish to include all the available data in one place.\nI have a few minor comments/corrections:\nL84. Insert ‘potentially’ before deleterious\nL94. More than?\nL306. ‘Expected’ or ‘accepted’ rather than excepted\nL529. Where does the figure of 300,000 come from? The phaCin-(like) elements are clearly older\nwhile the KoRV data seems rather soft (“younger estimated dates cannot be discarded”-L166)", "reviewer_2": "In their revised manuscript, the authors have addressed the questions raised in my previous review\nand modified the presentations where appropriate. I recommend publication of this work. It is\ncomprehensive, novel and appropriate for this journal."}, "author_response": "RESPONSE TO REVIEWERS' COMMENTS\nReviewer #1 (Remarks to the Author):\nI am happy to support publication of this manuscript in revised form although I still feel that a degree\nof simplification by focussing on the KoRVs would be advantageous. However, I do fully understand\nthe authors’ wish to include all the available data in one place.\nI have a few minor comments/corrections:\nL84. Insert ‘potentially’ before deleterious\nL94. More than?\nL306. ‘Expected’ or ‘accepted’ rather than excepted\nL529. Where does the figure of 300,000 come from? The phaCin-(like) elements are clearly older\nwhile the KoRV data seems rather soft (“younger estimated dates cannot be discarded”-L166)\nWe thank the reviewer for their kind revision and apreciate their understanding in our decision to\nkeep all the results together. All suggested modifications were applied to the final version of the\nmanuscritp.\nThe figure of 300,000 years comes from our current estimation of 312,191 years for KoRV. That is\ntrue, our conclusion was miswritten when considering that endogenization of phaCin-(like)\nelements is much older. We changed it now to 6.9 Mya. And we do agree that the estimation to\nKoRV colonization is still soft and is open for discussion. But most importantly is that it is a recent\nevent in comparison to the other retroviruses.\nReviewer #2 (Remarks to the Author):\nIn their revised manuscript, the authors have addressed the questions raised in my previous review\nand modified the presentations where appropriate. I recommend publication of this work. It is\ncomprehensive, novel and appropriate for this journal.\nThank you."}]

Decision letter