The admixed brushtail possum genome reveals invasion history in New Zealand and novel imprinted genes

Domain: earth_science
Published: 2023-10-17
Editorial outcome: major_revision

Abstract
Combining genome assembly with population and functional genomics can provide valuable insights to development and evolution, as well as tools for species management. Here, we present a chromosome-level genome assembly of the common brushtail possum (Trichosurus vulpecula), a model marsupial threatened in parts of their native range in Australia, but also a major introduced pest in New Zealand. Functional genomics reveals post-natal activation of chemosensory and metabolic genes, reflecting unique adaptations to altricial birth and delayed weaning, a hallmark of marsupial development. Nuclear and mitochondrial analyses trace New Zealand possums to distinct Australian subspecies, which have subsequently hybridised. This admixture allowed phasing of parental alleles genome-wide, ultimately revealing at least four genes with imprinted, parent-specific expression not yet detected in other species (MLH1, EPM2AIP1, UBP1 and GPX7). We find that reprogramming of possum germline imprints, and the wide epigenome, is similar to eutherian mammals except onset occurs after birth. Together, this work is useful for genetic-based control and conservation of possums, and contributes to understanding of the evolution of novel mammalian epigenetic traits.

Introduction
The common brushtail possum (Trichosurus vulpecula) is an nocturnal arboreal marsupial native to Australia. Like other marsupials, the brushtail possum gives birth to altricial young—gestation lasts only 17.5 days. Possums at birth are hairless, without overt sexual differentiation and do not even possess full respiration capacity in the lungs1. The first 100 days of life in the pouch are spent moving very little and with eyes shut, completing organogenesis and growth while continuously suckling from a teat2,3 (Fig. 1a). As such, it could be argued that the most major developmental transition for possums is at weaning and exit from the pouch around 6 months of age. Despite this, physiological changes associated with extended lactation and weaning have not been well studied, particularly in a genome-wide manner.

In addition to studies of their unique adaptations, marsupials are often used as a model for comparative evolutionary studies4–6. Of these, one of the most interesting is genomic imprinting—a form of monoallelic gene expression that is dependent on the sex of the parent from which an allele is inherited7,8. There are approximately 100–200 imprinted genes in humans and mice, characterised by paternal-specific expression of factors enhancing growth, and maternal-specific expression of those suppressing growth9. The most popular theory explaining why these genes would forgo the benefits of diploidy is known as the kinship, or parental-conflict, hypothesis10,11. The kinship hypothesis predicts genes with paternal-specific expression will maximise the use of maternal resources and increase the fitness of fathers. In contrast, maternally-derived genes act to control growth such that total reproductive output is maximised, irrespective of who the father is.

Marsupials like possum provide an interesting test for the parental-conflict hypothesis because there is limited opportunity for paternally-derived genes to influence maternal resource provision within the short-lived placenta7,12. Indeed, most of the genes that are imprinted in humans and mice appear to be not imprinted, or are without orthologues, in marsupials7,13–15; however, genome-wide searches for marsupial imprinted genes are limited16,17. Other than genes controlling imprinted X-chromosome inactivation4, genomic imprinting is yet to be studied in brushtail possums.

In addition to their value as a model for marsupial development and understanding mammalian evolution, possums are an iconic and protected species in Australia. For Aboriginal Australians, possums represent a cultural treasure, whereby cured pelts are used as a canvas upon which individual life-stories are depicted18. Possums can thrive in urban areas and are thus the native mammals most likely to come into human contact. Consequently, possums are a focus of relocation and rehabilitation efforts in areas where they are threatened19–21. These conservation efforts are hampered, at least to some extent, by the lack of genomic resources19.

In striking contrast to Australia, possums have destroyed many native forest ecosystems in Aotearoa New Zealand, following their first importation in 185822. Ecological damage occurs mainly through defoliation of trees23, but also through opportunistic invertebrate and bird predation24. Considering their additional role as a reservoir for bovine tuberculosis25, possums are often viewed as the most significant vertebrate pest species in New Zealand26. Efforts to control and locally eradicate possums have focussed on poisoning, trapping and shooting; however, more humane and effective approaches involving the suppression of possum fertility have been hampered by a lack of genetic information about the species.

To address this, we have produced a near-complete possum genome assembly as part of the Vertebrate Genomes Project, and undertaken a comprehensive functional genomics study. In addition to providing a basic understanding of marsupial pouch young development, we have identified novel communication genes, and imprinted loci not yet discovered in other species. By quantifying admixture in New Zealand possums, and tracing diversity to their source populations in Australia, we aid both pest management in New Zealand and conservation in their native range.

Results
Chromosomal assembly of the possum genome
The main brushtail possum genome assembly was generated from the forearm muscle of a 22-day-old male pouch young we called ‘Sandy’ (Fig. 1b). In addition to Sandy, tissue was collected from his sister (Sheila), a ‘backrider’ of approximately 180 days old, and mother (Puku) following trapping. This trio and many other possum samples we collected were from the Dunedin region in southern New Zealand—we concentrated our efforts here on account of it featuring multiple early possum introductions, both from the first site of establishment in nearby Southland, but also directly from Tasmania and mainland Australia22.

We generated the assembly following the Vertebrate Genomes Project pipeline v1.527, including haplotype phased contigs from PacBio long-read sequencing, polishing nucleotide errors with short-reads, scaffolding the contigs into chromosomes with 10x Genomics linked reads, Bionano Genomics optical maps, Arima Genomics Hi-C chromatin maps, and curation for structural errors and naming of chromosomes (GenBank accession GCA_011100635.1). Almost all sequence (99.88%) mapped to chromosomes, matching the karyotype28 (Fig. 1c), and possessing strong assembly rates (scaffold N50, 442,560,073 bp; contig N50, 4,314,688 bp). Gene annotation was assisted by long-read IsoSeq RNA-sequencing from muscle, brain and testis tissue (GenBank accession: GCF_011100635.1). Our assembly registered a BUSCO score of 96.8% complete genes, better than the most commonly used marsupial genomes (e.g., Monodelphis domestica, 91.6% complete)29, with the exception of the koala (98% complete)30,31 (Supplementary Table 1). Due to hemizygosity, the X-chromosome suffered increased fragmentation relative to the autosomes. Nevertheless, we assembled approximately ~1.28 Mbp of Y-chromosome sequence from 3 separate scaffolds that were specific to males (Fig. S1). At least one gene on these scaffolds, synaptonemal complex element protein 1 (SYCE1Y), is a novel Y-linked gene in mammals with differential expression compared to its X-linked homologue, SYCE1X.

RNA sequencing reveals unique physiological adaptations associated with altricial birth, extended lactation and weaning
To assist with functional annotation of the possum genome, we undertook RNA-sequencing from 13 of Sandy’s tissues. We also produced a further 103 transcriptome reactions from tissues (mostly liver) collected over a range of developmental stages and locations in the Dunedin region (southern New Zealand; Fig. 1d). Many genes upregulated in possum pouch young relative to adults appeared characteristic of fetal life in eutherian mammals. For example, alpha-fetoprotein (AFP; LOC118852573) is the major albumin protein during fetal life in humans and mice and its expression is rapidly silenced after birth32. Yet we saw high expression in possum pouch young, including juvenile pouch young >120 days post birth (Fig. S2a). Silencing of globin gene expression33 was also observed in early pouch young, thus reflecting the transfer of blood system production from the liver to bone marrow in post-natal marsupials34 (Fig. S2b; Supplementary Data 1).

Perhaps the most striking gene expression change we observed in the liver between pouch young and adults was associated with the metabolic effects of weaning (Fig. 2a; Supplementary Data 2). Genes associated with the Leloir pathway, involved with the catabolism of galactose and other carbohydrates, were high in pouch young, but low in adults. In contrast, cytochrome P450-family (CYP) genes associated with the first major stage of metabolic breakdown (i.e., oxidation) were low in pouch young but became highly expressed in adults. These gene expression changes reflect liver maturation and a metabolic transition to processing a range of dietary molecules at weaning. Interestingly, we found that the ‘back-rider’ Sheila (Sandy’s half-sister), had a low galactase/CYP ratio, indicative of weaning. This result provides molecular proof that at least some back riders stay with their mothers post-weaning, perhaps to learn browsing techniques.

To uncover molecular signals of unique marsupial traits associated with weaning, we searched for genes associated with the galactase/CYP expression ratio. The gene with the highest correlation to the galactase/CYP ratio encoded the imprinted long non-coding RNA, H19 (LOC118853838; Fig. 2b). H19 is thought to help orchestrate the paternal-specific expression of insulin-like growth factor 2 (IGF2) in mice35, an iconic and intensively studied imprinted gene. Similar regulation of IGF2/H19 has been proposed for the tammar wallaby, an Australian macropodid marsupial36,37. To explore this further in possum, we searched for single-nucleotide polymorphisms (SNPs) within H19 and IGF2 of Sandy and other individuals we collected, and tested their expression. We found exclusive maternal expression of H19, whereas IGF2 was paternally expressed in pouch young and biallelically expressed in adults where H19 expression was silenced (Fig. 2b, c; Fig. S2c).

Genes anti-correlated with the galactase/CYP weaning signature were also identified. Two of the most anti-associated were lipocalin family members (LOC118842298 and LOC118843263), known for binding and transporting small hydrophobic molecules such as steroid hormones, odorants, retinoids, and lipids38–40. Homology searches revealed that LOC118842298 encoded vulpeculin, previously known in possum only by its protein sequence41,42, with the other being a close relative (which we have named ganderin), situated amongst a mini-cluster of lipocalin genes on chromosome 3 (Supplementary Table 2). Vulpeculin is notable as being the most abundant protein in possum urine42, implying it acts in a similar manner to a large family of related lipocalins, the major urinary proteins in rodents, which are expressed in the liver but excreted in urine and heavily implicated in pheromone biology and chemical communication39. Interestingly, we found that expression of vulpeculin and ganderin is somewhat variable in females, but is consistent and much higher in males (Figs. 2d, S2d).

Early invasion history and admixture of possums in New Zealand
Brushtail possums are widely distributed within their Australian home range, with up to 6 recognised subspecies43. Of these, Trichosurus vulpecula fuliginosus from lutruwita Tasmania, and Trichosurus vulpecula vulpecula from south-east Australia, have well-documented importation into New Zealand22. The extent to which these subspecies are hybridising has not yet been fully described—elsewhere in New Zealand, certain populations maintain significant substructure, whereas others are freely interbreeding44–48. In all cases, genetic tracing to reveal the source Australian populations remains incomplete.

We hypothesised that possums from the Dunedin area may hold considerable genetic variation as this was an early and intensive site of possum release22. Indeed, mitochondrial sequencing from 3 Dunedin locations revealed 4 distinct mitochondrial haplotypes perhaps representing separate Australian source populations. One of these haplotypes clustered with possum samples originating near Sydney in New South Wales (NSW), with another closely related to the Tasmanian reference possum we sequenced (Fig. 3a). The other two haplotypes aligned to D-loop sequence from possums sourced from separate mainland populations (99.1% and 97.0% identity, respectively)48. We assessed nuclear genetic diversity from our RNA-sequencing data and reference samples49,50, finding standardised multilocus heterozygosity (sMLH) is higher on average in Dunedin possums than in the native Australian range, even when restricting comparisons to possums sharing the same mitochondrial haplotype. When compared to RNA-sequencing data from Australian mainland possums50, and low-coverage genomic sequencing of the Tasmanian individual, we found that diversity among Dunedin possums appears to be driven by admixture between diverse sources in Tasmania and mainland Australia (Fig. 3b).

Identification of novel marsupial imprinted genes
Relative to eutherian mammals, marsupial placentation affords less opportunity for manipulation of maternal physiology. As such, it is assumed that genomic imprinting is relatively reduced in marsupials, a proposition supported by biallelic expression or absence of many eutherian imprinted genes7,14–17,51. Nevertheless, the possibility that marsupials possess a novel set of imprinted genes that relate to their unique biology has not been tested outside of M. domestica16,17.

To perform an effective imprinted gene search, the ability to distinguish parental alleles over a large number of genes is essential. High levels of heterozygosity in Sandy and the other admixed New Zealand possums (Fig. 3a) provided a unique opportunity to search for novel marsupial imprinted genes. Initially, we tried searching for allele-specific expression in RNA-sequencing datasets, however, this resulted in a large number of candidate genes, many of which were false positives involving pseudogenes and repetitive regions52. Consequently, we altered our discovery pipeline to first identify allele-specific methylation (ASM) at candidate imprinting control regions53–55 using long-read nanopore sequencing (Fig. S3a). We started with an initial list of 173 sites displaying ASM (Supplementary Data 3), which included known marsupial imprint control regions (ICRs) such as intron 12 of IGF2R56 (Fig. S3b) and the PEG10 promoter57 (Fig. S3c). We then searched for monoallelic expression of SNPs expressed in nearby genes. As heterozygous SNPs were so clearly identified in Sandy’s genome due to 60-fold sequencing coverage, we started with an analysis of his transcriptome set.

In addition to known marsupial imprinted genes such as IGF2R and PEG10, we found 3 sites of ASM (Fig. 4a–c) associated with monoallelic expression of adjacent or overlapping genes in all of Sandy’s tissues (Fig. 5a–d; Fig. S4a). This included the MLH1 and EPM2AIP1 genes linked to a single intervening region of ASM, as well as single gene loci featuring UBP1 and GPX7 (Fig. 5a–d).

To further test the monoallelic nature of these genes, we then examined our RNA-sequencing datasets from other individuals. Although significant allelic variation was found in these genes, no single individual appeared to express both alleles—a violation of the Hardy-Weinberg equilibrium (Supplementary Data 4). Amplicon sequencing was then used to definitively prove which individuals were heterozygous despite showing expression of only one allele (Fig. 5e–h; Fig. S4b–e; Supplementary Data 4).

While monoallelic expression was clear for MLH1/EPM2AIP1, there appeared to be biased expression of GPX7 and UBP1. Moreover, our experiments up to this point had not definitively shown allelic dependence on parental sex. To address this, we genotyped the mothers of heterozygous offspring at our loci of interest—in instances where a mother was homozygous, this allowed us to assign the parental origin of alleles. In every case, we confirmed that GPX7 was expressed from the maternal allele (n=7), whereas MLH1 (n=5), EPM2AIP1 (n=5), and UBP1 (n=4), were all expressed from the paternal allele (Table 1; Supplementary Data 4).

Global methylation erasure in primordial germ cells coincides with the reprogramming of imprint control regions
In humans and mice, epigenetic marks acquired in one generation are generally not passed on to the next because of genome-wide epigenetic erasure in the germline58,59. By far the most distinctive methylation erasure event occurs in primordial germ cells (PGCs), in which high starting levels of methylation (70-80%) are removed through active and passive means, leaving only 5% methylation at its lowest point60,61. In contrast, while there is an alteration of methylation patterns in the germline of non-mammalian vertebrates like zebrafish, global methylation erasure in PGCs apparently does not occur62–66. This provides for the interesting possibility that acquired traits underpinned by DNA methylation can be inherited in at least some vertebrate groups58,67,68.

As a divergent mammalian group, methylation dynamics in the marsupial germline are somewhat unclear. Global loss of PGC methylation has been observed in wallaby69, however, the immunofluorescence techniques used were not quantitative. To explore genome-wide methylation changes in possum development, we isolated PGCs from the gonads of male pouch young up to 106 days of age (Fig. 6a). Low-coverage bisulfite sequencing was used to estimate global DNA methylation levels throughout post-natal development. Methylation levels were high immediately after birth (69.7%), but dropped to 41.1% by 10 days postpartum (dpp) (Fig. 6b; Supplementary Data 5). Subsequently, remethylation occurred such that at 25 dpp, 60-70% methylation was observed (t-test, p=0.001), consistent with adult sperm.

To establish if the novel imprinting control regions we identified are reprogrammed during this period of methylation erasure, we undertook deep sequencing of two libraries (13 and 16 dpp), and compared these to libraries of control somatic cells from the same individuals (i.e. non-germline gonadal cells) and adult sperm (Fig. 6c). From this analysis, it was clear that the control somatic cells were 50% methylated over the regions we found that had ASM in Sandy (Fig. 4a–c; Fig. S3b, c and Fig. S5b). In contrast, PGC methylation in the same regions was 5.06-20.7%, implying that the erasure of methylation imprints in the germline coincides with global DNA demethylation.

Interestingly, sperm methylation at the ASM sites was even lower than PGCs for IGF2R, PEG10, MLH1 and UBP1 (1.39–1.80%); however, H19 and GPX7, were fully methylated (94.24–95.79%) (Fig. 6c). This supports our findings of paternal- and maternal-specific expression for these genes, respectively (Fig. 5a–h; Fig. S4b–e; Supplementary Data 4), and implies the ASMs we discovered to act as germline imprinting control regions (ICRs). Moreover, our data for PEG10 and IGF2R is consistent with tammar wallaby56,57.

To identify other putative ICRs with methylation patterns indicative of imprinting, we assessed methylation levels of the remaining ASM sites we identified in the PGC, somatic cell and adult sperm datasets (Fig. 6d, top panel; excluding those ASM sites with missing data, n=32). We found 41 ASM sites had somatic methylation levels that fall within the range exhibited by germline ICRs (30–57%). Of these, 9 ASM sites were associated with (+/− 1 Mb distance) 13 genes displaying monoallelic expression in at least Sandy’s brain, liver or muscle (Fig. 6d, bottom panel; Supplementary Data 3). While not as robustly tested as MLH1/EPM2AIP1, GPX7 and UBP1, we nonetheless consider these ASM sites and genes to be strong imprinted gene candidates worthy of further testing.

Methods
Sample collection and metadata
Tissue samples were sourced from freshly deceased possums killed as part of a pest-control programme, and therefore not require animal ethics oversight according to guidelines issued by the NAEAC, Occasional Paper No 2, 2009 (ISBN 978-0-478-33858-4). Additional samples were collected from unused tissue belonging to captive possums euthanised from another study. Capture, husbandry and manipulation in that study were approved by the Otago Animal Ethics Committee (AUP-19-75, AUP-20-10). Further samples for mitogenome analysis were sourced from a prior study47.

Pouch young was developmentally staged using the nomogram of Lyne and Verhagen2. Pouch young >120 days of development were challenging to stage, and for this reason, have been referred to as juveniles. Upon collection, tissue was added to RNAlater (Invitrogen, AM7020) and placed on ice until storage at +4 °C for 24 h, followed by long-term storage at −80 °C. Alternatively, tissue was processed in order to collect germ cells (see below), or samples for DNA extractions (from ears, tails or limbs) were snap-frozen in liquid nitrogen and stored at −80 °C (see below). Sperm were isolated from the epididymis of an adult possum by allowing sperm to diffuse out into warm 1x PBS (37 °C) for 1 h. Plasma was removed by centrifugation at 240 × g for 2 min and washing with 1x PBS, repeated twice. The resulting pellet was resuspended in 1x PBS, passed through a 20 µm cell strainer and stored at −80 °C. Metadata for each sample used in this work is given in Supplementary Data 5.

DNA analyses
DNA extraction and library preparation for de novo genome assembly. Ultra-high molecular weight DNA (uHMW) was extracted using the agarose plug Bionano Genomics protocol from Sandy’s forearm muscle. uHMW DNA quality was assessed by a Pulsed Field Gel assay and quantified with the Qubit DNA High Sensitivity dsDNA Kit and Qubit Fluorometer (Invitrogen, Q32854 and Q33238).

PacBio continuous long reads (CLR). Five micrograms of uHMW DNA were sheared using a 26G blunt end needle (PacBio protocol, 101-181-000, Version 05). A large-insert PacBio library was prepared using the Pacific Biosciences Express Template Prep Kit v1.0 (101-357-000) following the manufacturers’ protocol. The library was then size-selected (>15 kb) using the Sage Science Blue Pippin Size-Selection System. The final library was sequenced on PacBio Sequel smrtcells (~60x coverage).

10x Genomics linked reads. The same unfragmented uHMW DNA was used to generate a linked-reads library with the 10x Genomics Chromium (Genome Library Kit & Gel Bead Kit v2, 120258; Genome Chip Kit v2, 120257; i7 Multiplex Kit, 120262). This library was then sequenced on an Illumina Novaseq S4 (150 bp paired end reads).

Bionano genomics optical maps. Unfragmented HMW DNA was also used for Bionano Genomics optical mapping. Briefly, DNA was labelled using the Bionano Prep Direct Label and Stain (DLS) Protocol (30206E) and run on one Saphyr instrument chip flow cell. 1.1 Gb of data was generated and optical maps were assembled using Bionano Access.

Arima Genomics Hi-C proximally ligated DNA. Hi-C libraries were generated by Arima Genomics (https://arimagenomics.com/) using the Arima-HiC kit v1 (A510008) and sequenced on a HiSeq X (Illumina) at ~147x coverage following the manufacturer’s protocols.

DNA extraction for all other analyses. DNA was extracted from tissue (~3 mm tissue punches—liver, ear or tail) using a Bio-On-Magnetic-Beads (BOMB) protocol that utilises solid-phase reversible immobilisation (SPRI) carboxyl-coated Sera-Mag Magnetic SpeedBeads (Cytvia, 45152105050250)84. Tissue was lysed in 300 µL TNES buffer (100 mM Tris pH 8.0, 25 mM NaCl, 10 mM EDTA, 1% w/v SDS) with 5 µL Proteinase-K (20 mg/mL) at 55 °C overnight. Lysates were then mixed with 6 M Guanidinium Thiocyanate (GITC), SPRI beads suspended in 1x TE buffer (10 mM Tris pH 8.0, 1 mM EDTA), and isopropanol in a ratio of 1:2:2:4, respectively. For DNA extraction from sperm, the frozen pellet (described above) was resuspended in GITC supplemented with 10% v/v ß-mercaptoethanol as a reducing agent85. For all tissues, after a 5 min incubation to allow DNA to bind to the magnetic beads, tubes were placed on a neodymium magnetic rack for ~5 min until the solution was clarified. The supernatant was removed and the DNA bound to the magnetic beads was washed once with isopropanol and twice with 70% ethanol, and then left to air dry on the magnetic rack. DNA was eluted from the magnetic beads in nuclease-free water or 1x TE buffer and concentrations were measured using the Qubit DNA High Sensitivity dsDNA Kit and Qubit Fluorometer.

PCR amplification. Amplification of DNA for Y-chromosome analyses, mitochondrial DNA library preparations and amplicon sequencing for imprinted SNP analysis (see below) was performed using Phusion High-Fidelity DNA Polymerase (New England Biolabs, M0530L) according to the manufacturers’ protocol. Approximately 10–25 ng of DNA was used in the following PCR mix: 1x HF Buffer with 1.5 mM MgCl2, 0.2 mM of each dNTP, 0.5 µM of each primer, and 0.02 U/µL units of Phusion. Primer sequences and PCR cycling conditions used for each amplification are specified in Supplementary Table 3.

Preparation and sequencing of mitochondrial DNA libraries. The mitochondrial genomes not assembled from RNA-sequencing data were amplified using five overlapping primer pairs (Supplementary Table 3). For each sample, the five PCR products were pooled in equivolume amounts and cleaned up and size-selected using 0.8X ratio of SPRI beads diluted in standard PEG buffer (18% w/v polyethylene glycol 8000 (PEG), 1 M NaCl, 10 mM Tris (pH 8.0), 1 mM EDTA, 0.05% v/v Tween-20)84. The concentration of each resulting PCR pool was determined using the Qubit DNA High Sensitivity dsDNA Kit and Qubit Fluorometer. The NEBNext Ultra II FS DNA Library Prep Kit for Illumina (New England Biolabs, E7805L) was used with 50 ng of each PCR pool to make a DNA library with the following modifications: the reaction volumes were reduced to a ¼ of the standard reaction volume, a fragmentation/end prep time of 12 min, 1.5 µM Y-adapter (see below) was used for adapter ligation. The resulting libraries were amplified with 0.4 µM indexed Truseq-type oligos for Illumina sequencing with the following conditions: 98 °C for 30 s, followed by 5 cycles of 98 °C for 10 s, 65 °C for 75 s, and a final elongation of 65 °C for 5 min. The unique barcode for each sample is given in Supplementary Data 5. The final libraries were pooled in equimolar amounts and sequenced on the iSeq100 (Illumina) to generate 150 bp paired-end reads.

The Y-adapter includes the PE_Yadapter_1 (/5Phos/GATCGGAAGAGCACACGTCTGAACTCCAGTC) and PE_Yadapter_2 (ACACTCTTTCCCTACACGACGCTCTTCCGATC*T) which were pooled 1:1 in a 100 µL volume to make a 50 µM stock. Added to this pool was 1 µL of 5 M NaCl (50 mM final) for adapter annealing with the following conditions: 95 °C for 2 min, followed by ramping 95 °C to 25 °C over 45 min, at 0.1 °/s (or 6 °/min). The annealed Y-adapter was then diluted to 1.5 µM.

Preparation and sequencing of nanopore libraries. Following DNA purification from Sandy’s forearm tissue (BOMB extraction protocol, above), nuclease-free water was added to a final volume of 400 μL in order to reduce contaminants. The sample was treated with 2 μL of 10 mg/mL RNase A (to remove any remaining RNA) and placed at 37 °C for 5 min. To remove the RNase A and any remaining proteins, the sample was treated with 2 μL of 20 mg/mL Proteinase-K at 50 °C for 5 min. Chloroform:isoamyl alcohol (IAA) (24:1) was then added (400 μL) to the DNA sample, which was mixed by gentle rotation at room temperature for 5 min. The sample was then spun at 15,320 × g for 5 min and the top phase (~350 μL) was transferred to a new tube containing 350 μL Chloroform:IAA (24:1). The sample was mixed again by gentle rotation at room temperature for 5 min, followed by centrifugation at 15,320 × g for 5 min. The top phase (~300 μL) was transferred to a new tube and 0.1X volume (30 μL) of 3 M NaOAc pH 5.2 followed by 330 μL of 100% ice-cold ethanol was added and the tube was gently inverted. The sample was then spun at 15,320 × g for 5 min and the supernatant was removed. The resulting pellet was washed with 1 mL ice-cold 70% ethanol and spun at 5200 × g for 5 min. The ethanol wash was removed and the pellet was air-dried before resuspending the DNA in 100 μL 1x TE buffer. Sandy’s DNA was sequenced on five R9.4.1 MinION flow cells (Oxford Nanopore Technologies) prepared with the Rapid Sequencing Kit (SQK-RAD004; run 1, with 400 ng of DNA) or the Ligation Sequencing Kit (SQK-LSK110; run 2-5, with 1 μg of DNA). For run 2-5 with DNA prepared with the Ligation Sequencing Kit, ~15 fmol of DNA was loaded onto the MinION flow cell.

Amplicon sequencing for imprinted SNP analysis. A dual-indexing, four-primer PCR-based assay was used for amplicon sequencing86. Primers contain a ‘handle’ sequence such that during amplification the handle is incorporated into the amplicon (Supplementary Table 3). Following the first round of PCR amplification (as described above), ~¼ of the PCR reaction was cleaned up and size-selected using a 0.9x SPRI beads diluted in standard PEG buffer84. The DNA was eluted in 10 µL nuclease-free water and 2–5 µL was used as a template in a second round of PCR amplification (as described above) with 0.2 µM of each second step primer (indexed Truseq-type oligos for Illumina sequencing; Supplementary Table 3). PCR cycling parameters for the second step PCR were: 98 °C for 3 min, 5 cycles of 98 °C for 10 s, 62 °C for 20 s, and 72 °C for 15 s. A final elongation step was performed for 2 min at 72 °C. The final libraries were pooled in equimolar amounts, cleaned-up and size-selected using a 0.9x SPRI beads diluted in standard PEG buffer and sequenced on the iSeq100 (Illumina) to generate 150 bp paired-end reads.

Post bisulfite adapter tagging (PBAT) of sperm DNA. Bisulfite-converted genomic libraries with sperm DNA were prepared using a modified PBAT method87. Approximately 50 ng of DNA was bisulfite-converted using the EZ Methylation-Direct MagPrep kit (Zymo, D5044), but reducing the reaction volume to one ¼ of that stated in the protocol (to allow for all reaction and clean-up steps to be performed in PCR tubes). Bisulfite-converted DNA was eluted from the magnetic beads in 20 µL elution buffer and transferred to a new tube for first-strand synthesis with a 5′-biotinylated adapter primer bearing seven random nucleotides at its 3′ ends (BioP5N7, biotin-ACACTCTTTCCCTACACGACGCTCTTCCGATCTNNNNNNN). The first strand synthesis product was then purified using Dynabeads™ M-280 Streptavidin beads (Invitrogen, 11205D) and alkaline denaturation. Second-strand DNA synthesis involved the immobilised first-strand DNA and a second adapter primer also bearing seven random nucleotides at its 3′ ends (P7N7, GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCTNNNNNNN). The dsDNA product was then PCR amplified (13 cycles) using 1X HiFi HotStart Uracil+Mix (Roche, KK2802) with 0.4 µM primer (indexed Truseq-type oligos for Illumina sequencing), so that sample-specific barcodes and sequences required for Illumina flow-cell binding were added to the library. The unique barcode for this sample is given in Supplementary Data 5. The integrity of the library was assessed by agarose gel electrophoresis before size-selecting using standard PEG-diluted SPRI beads (0.8X). The pooled library was sequenced on the NextSeq 2000 (Illumina) at the Otago Genomics Facility to generate 50 bp single-end reads.

Gonad collection and germ cell labelling. Gonads were collected from male pouch young (n=17; Supplementary Data 5) and germ cells were isolated from the gonads as previously described88. Briefly, gonads were torn using needles in a drop let of TrypLE Express (Gibco, 12605010) and incubated at 37 °C for 20 min. Tissue clumps were triturated manually by gentle strokes with a P1000 pipette tip and then with a P200 pipette tip. TrypLE Express was diluted with 5% FBS (in PBS), filtered through a 40 µm cell strainer, and washed twice in PBS (centrifugation 240 × g for 2 min). The resulting pellet was resuspended in Zombie Near-Infra Red Viability Dye (BioLegend, 423105) and incubated at room temperature for 20 min. The lysate was washed in 5% FBS (in PBS), then resuspended in SSEA1 monoclonal antibody (Invitrogen, MC-480) conjugated to DyLight 488 (1:300 in 5% FBS (in PBS); Invitrogen, MA1-022-D488) on ice for 30 min. The lysate was then washed, resuspended in sort buffer (5% FBS, 0.5 mM EDTA in PBS), and kept on ice prior to cell sorting. Validation of the SSEA1 monoclonal antibody for germ cell specificity in brushtail possum was previously undertaken74.

Germ cell isolation using fluorescence-assisted cell sorting (FACS). Germ cells were isolated from the cell lysate using a BD FACSAria Fusion which detects GFP using a 488 argon laser (BD Biosciences). A gonadal cell profile was produced using forward scatter (FSC) and side scatter (SSC) and was gated to exclude doublets and dead cells. Germ cells were identified as a small sub-population that was very bright (GFP > 104) with higher SSC (Fig. 6a and Supplementary Data 6; approx. 0.1% of total cells). This population was gated conservatively, sorted into 0.2% Triton-X and stored at −80 °C. The number of GFP-positive sorted cells varied between 3 and 516 per sample (Supplementary Data 5). GFP-negative cells (somatic cells) were also sorted from each individual as negative controls.

Low cell input bisulfite sequencing. Bisulfite conversion and library preparation followed the single-cell bisulfite sequencing protocol89, substituting modified PEG-diluted SPRI beads for AMPure XP beads. This protocol is similar to PBAT (above) but carries SPRI beads through most steps with the addition of modified PEG buffer (18% w/v PEG 8000, 2.5 M NaCl, 10 mM Tris (pH 8.0), 1 mM EDTA, 0.05% v/v Tween-20) for size-selection steps. GFP-positive germ cells and GFP-negative somatic cells from each individual were incubated with 0.5 µL of 20 mg/mL Proteinase-K on ice for 15 min followed by the addition of modified PEG-diluted SPRI beads for a further incubation of 20 min at room temperature. Tubes were placed on a neodymium magnetic rack for ~5 min until the solution was clarified, the supernatant removed and beads washed twice in 90% ethanol. The DNA was eluted into 10 µL of 10 mM Tris (pH 8.0), but kept with the beads to avoid loss of the DNA. Bisulfite conversion was then performed as for PBAT (see above), carefully adding the conversion reagent to the DNA/bead suspension to prevent disturbing the beads. The resulting bisulfite-converted DNA underwent five rounds of first-strand synthesis using the 5′-biotinylated adapter primer (BioP5N7) with intermediate heat denaturation of the DNA fragments89. Exonuclease I (New England Biolabs, M0293S) was added to the first strand synthesis product to remove any remaining adapter primer, and the product was cleaned up and size-selected with modified PEG-diluted SPRI beads (0.8X). Second-strand synthesis was performed in the presence of SPRI beads and was otherwise the same as for PBAT. The second strand product (50 µL) was stimulated to bind the SPRI beads by the addition of 50 µL of 10 mM Tris-HCl and 80 µL of modified PEG buffer (0.8x), followed by on-bead PCR amplification as for PBAT (14 cycles for 100 cells, 16 cycles for <100 cells). The unique barcode for each sample is given in Supplementary Data 5. PCR products were size selected by addition of modified PEG buffer (0.8x) and integrity of each library was assessed by agarose gel electrophoresis. Libraries were then pooled in equimolar amounts and sequenced on the iSeq100 to generate 150 bp paired-end reads. Samples of GFP-positive cells (and their GFP-negative somatic cells) (n=2; 13 and 16 dpp) were further sequenced (high coverage) on the NextSeq 2000 as for the sperm sample (above).

RNA analyses
RNA extraction and library preparation for de novo assembly. In order to assist with annotation of the genome, long-read RNA-sequencing was performed on muscle, brain and testis tissue using the PacBio IsoSeq method as per VGP pipelines27.

RNA extraction and preparation for all other analyses. Total RNA was extracted from tissue samples (~50 mg) that were stored in RNAlater using TRIzol Reagent (Invitrogen, 15596026) following the manufacturer’s guidelines. RNA samples were treated to remove any contaminating DNA using the Turbo DNA-free kit (Invitrogen, AM1907) according to the protocol. The concentration of DNase-treated RNA was determined using the Qubit RNA High Sensitivity Kit (Invitrogen, Q32852) with the Qubit Fluorometer.

Library preparation of short-read RNA-sequencing libraries. Poly(A)-enriched, strand-specific RNA-sequencing libraries were prepared using either the NEXTFLEX Rapid Directional RNA-Seq Kit 2.0 and NEXTFLEX Poly(A) Beads (2.0) (Perkin Elmer, NOVA-5198 and NOVA-512991) or the NEBNext Ultra II Directional RNA Library Prep Kit for Illumina and NEBNext Poly(A) mRNA Magnetic Isolation Module (New England Biolabs, E7760L and E7490L) according to the manufacturers’ guidelines with the following modifications. For the NEXTFLEX kit, 400 ng of DNase-treated RNA was used as input, fragmentation time was 12 min, 0.48 µM of NEXTFLEX RNA-Seq 2.0 Unique Dual Index Barcodes (Perkin Elmer, NOVA-512920) was used for adapter ligation and 12 cycles for PCR amplification. The Unique Dual Index for each sample is given in Supplementary Data 5. The final libraries were pooled in equimolar amounts and run on 5x lanes of the HiSeq 2000 (Illumina) at the Otago Genomics Facility to generate 100 bp single-end reads. For the NEB kit, the reaction volumes were reduced to a ¼ of the standard reaction volume with 250 ng of DNase-treated RNA used as input, fragmentation time was 12 min, 3.125 µM of NEXTFLEX RNA-Seq 2.0 Unique Dual Index Barcodes (Perkin Elmer) was used for adapter ligation and the following modified PCR protocol because a non-NEB adapter was used: 98 °C for 30 s, followed 11 cycles of 98 °C for 10 s, 60 °C for 30 s, 65 °C for 45 s, and a final elongation of 65 °C for 5 min. The Unique Dual Index for each sample is given in Supplementary Data 5. The final libraries were pooled in equimolar amounts and run on the NextSeq 2000 at the Otago Genomics Facility to generate 50 bp single-end reads.

cDNA synthesis and reverse transcription-PCR (RT-PCR) amplification. For cDNA synthesis, 500 ng of DNase-treated RNA was used with the RevertAid First Strand cDNA Synthesis Kit (Thermo Scientific, K1632) according to the manufacturers’ guidelines for random hexamer primed synthesis and GC-rich RNA templates. For all RNA samples processed, a minus reverse transcriptase (-RT) control was included to determine if any resulting amplification could be from genomic DNA. No amplification was seen in any of the -RT control samples tested.

Amplification of cDNA for RT-PCR analyses was performed using Phusion High-Fidelity DNA Polymerase (New England Biolabs) according to the manufacturers’ protocol. Each cDNA sample was diluted 2-fold with nuclease-free water and 2 µL was used in a 20 µL volume in the following PCR mix: 1X HF Buffer with 1.5 mM MgCl2, 0.2 mM of each dNTP, 0.5 µM of each primer, and 0.02 U/µL units of Phusion. Primer sequences and PCR cycling conditions used for each amplification are specified in Supplementary Table 3. In order to sequence the resulting PCR products and determine allele-specific expression, the amplicons were cleaned-up size-selected, and processed as for amplicon sequencing for imprinted SNP analysis (see above).

Bioinformatic analyses
Reference samples (RNA-sequencing and whole-genome sequencing (WGS)) were included in various bioinformatic analyses and are outlined in Supplementary Data 5. These included samples from the Australian locations Karta Pintingga (Kangaroo Island; KAN), lutruwita (Tasmanian; TAS), and New South Wales (NSW)50, as well as Aotearoa New Zealand (NZ)49.

Whole-genome sequencing analysis. Paired FASTQ files for Sandy’s WGS data were obtained from GenomeArk (https://genomeark.s3.amazonaws.com/index.html?prefix=species/Trichosurus_vulpecula/mTriVul1/genomic_data/10x/). For the TAS WGS sample, an additional initial step removing 10X adaptors was performed using an in-house script. Paired reads were trimmed using Trim Galore! (v0.6.7) in a two-step process: first, to remove adapters, and second, to remove low-quality base calls (Phred score < 20) (https://github.com/FelixKrueger/TrimGalore). The reads were mapped to the possum reference genome GCA_011100635.1 with the BWA-MEM algorithm (v0.7.17)90. The SAM files obtained were sorted and converted to BAM files using SAMtools (v1.15.1)91. To determine the sequencing depth of Sandy’s WGS data, non-overlapping 1 kilobase windows were created for all chromosomes and contigs using the R package GenomicRanges (v1.38.0) and the number of reads mapping to each window was determined with the coverage command from BEDtools (v2.30.0)92.

RNA-sequencing expression and coverage analysis. Raw reads were trimmed using Trim Galore! (v0.6.7) to remove adapters and low-quality base calls (Phred score < 20) (https://github.com/FelixKrueger/TrimGalore). Trimmed reads were mapped to the possum reference genome GCA_011100635.1 using HISAT2 (v2.2.1)93.

For expression analysis, the HISAT2 generated BAM files were imported into SeqMonk (v1.48.0) (https://www.bioinformatics.babraham.ac.uk/projects/seqmonk/) for visualisation. The SeqMonk RNA-seq quantitation pipeline was employed with uniquely mapped reads to generate a set of probes covering every gene in the genome, which were quantitated based on the number of reads falling within this annotation. The resulting log-transformed counts were corrected for the total number of sequences in each dataset to generate an RPKM (log2) value for each gene. Two samples (020721_PY01 and 210521_PY01_liver) were removed from further expression analyses due to poor mapping efficiencies and RNA QC plot analysis assessed via SeqMonk (v1.48.0)).

SNP calling using short reads. Short-read SNP calling of RNA-sequencing and WGS-mapped data (see mapping details above) was performed following the GATK’s Best Practices (v4.2.6.1)94. Duplicate reads (from the BWA-MEM mapping) were marked using the MarkDuplicates GATK module and read group tags were added using the addreplacerg function from SAMtools (v1.15.1)91. The SNPs were called into gVCF files for each sample using the HaplotypeCaller GATK module. Individual gVCFs were combined using the CombineGVCFs GATK module and genotypes of the samples were obtained with the GenotypeGVCFs GATK module. High-quality SNPs were filtered with VCFtools (v0.1.15) using the following parameters: --minQ 30 --min-meanDP 10 --max-missing 0.9595. After running this pipeline, a genome-wide dataset comprising 92 samples and 2310 SNPs was obtained.

Mitochondrial genome assembly. Reads from mitochondrial DNA libraries were trimmed using Trim Galore! (v0.6.7) to remove adapters, low-quality base calls and short reads (Phred score < 20, length < 35) (https://github.com/FelixKrueger/TrimGalore). Trimmed reads from mitochondrial DNA libraries and RNA-sequencing libraries (see above) were then mapped to the possum reference mitochondrial genome (NC_003039.1) using BWA-MEM algorithm (v0.7.17)90. The SAM files obtained were sorted and converted to BAM files using SAMtools (v1.15.1)91. For RNA-sequencing derived assemblies a FASTQ consensus sequence was generated from each of these files using the mpileup command from BCFtools (v1.15.1) and seqtk (v1.3) was used to mask (replace with ‘N’) all base calls with Phred score values less than 20. Consensus sequences of assembled mitochondrial genomes are given in Supplementary Data 7.

Mitochondrial genome phylogenetic and spatial analyses. Assembled mitogenomes derived from RNA-sequencing and mitochondrial DNA libraries were aligned together with KAN, TAS and NSW reference samples using the default Geneious v2022.2.2 (https://www.geneious.com) algorithm. Because samples derived from RNA-sequencing had significantly higher levels of missing data for the control region/origin of the light strand, all tRNAs, and both rRNAs, phylogenetic analyses were restricted to protein-coding regions (although some genes had significantly more missing data in RNA-sequencing samples, the maximum missing data for any sample for any protein-coding gene was 1.24%). Protein-coding genes were scrutinised to ensure no premature stop codons/frame-shift mutations were present, indicative of inadvertent incorporation of nuclear-embedded mitochondrial DNA sequences (NUMTs), and maternally-related individuals were screened to ensure they had identical haplotypes96.

Phylogenetic trees were built utilising the codon-partitioned mitochondrial DNA protein-coding alignments. Maximum likelihood (ML) trees were constructed via RAxML v8.2.1297 with 100 bootstrap replicates98. Two separate runs were used to confirm convergence, utilising the GTRCAT model of substitution. After confirming topological convergence based on the bipartitions file, we utilised Run 1 for analyses as it had the highest likelihood. We verified our ML results using BEAST v2.6.699, again implementing two runs to check the convergence of 250,000,000 states, sampling every 10,000 states. We implemented a gamma site model (4 categories, initial shape 1.0, estimated), and TN93 with Kappa 1 and Kappa 2 estimated (initially set at 2.0). We also estimated site frequencies and a strict clock as we did not expect marked rate variation within species. We used a Coalescent Bayesian Skyline100 as our tree model due to uncertainty about the population demography and linked the tree model across all partitions while allowing site and clock to vary. Following the removal of the initial 10% of states as a burn-in, we checked for convergence of parameters in Tracer v 1.7.2101 and of topologies of the consensus trees constructed with TreeAnnotator (maximum clade credibility tree, common ancestor heights). After confirming convergence, the two runs were combined, checked for adequate ESS values, and a combined consensus tree was created.

The reference samples were used to define the origin of mitochondrial haplotypes within our Dunedin samples as either NSW, TAS, KAN, vs. ‘other well-supported clades’ not represented in our Australian full-length reference samples (note, although these latter clades did not match full-length mitochondrial reference sequence we had, they could be identified as being most similar to D-loop sequence originating from south-east mainland Australia48). Only one discrepancy between runs/methods was detected: Peninsula2018_022 fell within the NSW clade in both RAxML runs, but differed between the BEAST runs (NSW and TAS), with the combined BEAST run finding this sample falling within the TAS clade. However, the split between NSW and TAS clades was not strongly supported in our analyses (<61% bootstrap support). Given the overall concordance between methods, we based downstream analyses (i.e. haplotyping) on the RAxML Run 1 tree. We used R v4.2.1 and the tidyverse v1.1.2, ggmap v3.0.0, and ggrepel v0.9.1 packages to examine the spatial distribution of mitochondrial haplotypes in the greater Dunedin region98.

Nuclear population genomics. Standardised multilocus heterozygosity (sMLH) was calculated for each of the samples using variant calls derived from the RNA-sequencing samples with R utilising the inbreedR, vcfR, reshape, and future.apply packages. Comparisons of sMLH between possums with different mitochondrial haplotypes were carried out in R98 with the tidyverse and rstatix v0.7.0 packages.

To calculate levels of admixture among Dunedin possums we used the high-quality SNP data obtained from the RNA-sequencing and WGS datasets. First, we removed closely related individuals using three consecutive rounds of sample filtering, setting a threshold of PI_HAT scores > 0.25 obtained by PLINK v1.90b6.21102. This filtering was conducted separately for the NSW and Dunedin datasets. After sample removal, we obtained a dataset comprising 60 individuals and 2248 autosomal SNPs. Principal component analysis (PCA) was run on the filtered dataset using the pca command in PLINK. The number of genetic clusters best fitting the data was calculated by performing a 10-fold cross-validation in ADMIXTURE v1.3.0 (from k=2 through k=10), resulting in the highest support for two genetic clusters (k=2)103. The results were drawn using the ggplot v3.3.6 package in the R language.

Identification of allele-specific methylation by nanopore sequencing data. The raw FAST5 files obtained for each sequencing run were base-called using Guppy v6.2.1 with the following parameters: guppy_basecaller -i FAST5 -c dna_r9.4.1_450bps_sup.cfg -x “cuda:0” -s FASTQ. Sequencing reads were trimmed to remove adapters and low-quality reads using Porechop (v0.2.4) and Filtlong (v0.2.0), respectively (https://github.com/rrwick/Porechop; https://github.com/rrwick/Filtlong). Reads with internal adapters, average Phred score < 9 or length < 500 were discarded. Trimmed reads were mapped to the possum genome reference GCA_011100635.1 with Minimap2 (v2.24) using the “map-ont” parameter104 resulting in 14.05X coverage (Supplementary Table 4). The bedcov tool from SAMtools (v1.15.1)91 was used to calculate the proportion of the genome that was covered at ≥ 5x and 10x depth, which was estimated to be 98% and 80.1%, respectively. CpG methylation calls were extracted using the call-methylation module from Nanopolish (v0.13.3)105. To perform methylation phasing we first identified SNPs using Clair (v2.1.1)106. Clair parameters for SNP detection were set at a threshold of 0.2 and the 122HD34 was selected as the precomputed model. SNVoter (v.1.0) was used to improve SNP detection from low-coverage regions54. Subsequently, SNPs were phased using WhatsHap (v1.1) with default parameters and the --ignore-read-groups option107. CpG methylation calls, variant calls and phasing information were used to obtain phased CpG methylation and mock bis-sequencing BAM files using NanoMethPhase (v1.0)54. Mock bis-sequencing BAM files were visualised with IGV (v2.12.3) using the bisulfite colour alignment mode108. Differentially methylated CpG sites between haplotypes were identified using the DMLtest function of the DSS package (v2.34.0). Finally, differentially methylated CpG sites, called here allele-specific methylation (ASM) sites, were filtered at a threshold of ≥ 20 CpG and annotated with the closest gene/mRNA (± 20 kb) using SeqMonk (v1.48.1) (https://www.bioinformatics.babraham.ac.uk/projects/seqmonk/).

Monoallelic analysis of known and candidate genes. To investigate the monoallelic expression of known and novel imprinted genes, SNPs within these genes were identified following SNP analysis of the nanopore datasets (see above; Supplementary Table 5) and confirmed by PCR and amplicon sequencing (see above). The SAMtools (v1.9)91 mpileup utility (max depth 0) was then used to provide a summary of the depth of mapped reads at the SNP for each RNA-sequencing library (Supplementary Data 4).

Identification of genome-wide monoallelic SNP expression. SNP calling to identify heterozygous sites in the reference genome was performed using a similar approach to RNA-sequencing SNP calling, with the additional step of using the high-quality SNPs obtained from RNA-sequencing data for variant recalibration and filtering. In total, we obtained 11,481,987 variants for Sandy’s reference genome (0.35%). Next, we used the CollectAllelicCounts function from the GATK toolkit (v4.2.6.1) to quantify the number of reads with the reference and alternative alleles in the RNA-sequencing data94. Allelic count matrices were obtained for brain, liver and muscle. To identify monoallelic expressed SNPs, we compared the observed counts with the expected counts assuming biallelic expression (1:1). The thresholds were a minimum depth of 10x in RNA-sequencing data and a Fisher’s exact P-value < 0.05 after Bonferroni correction for the observed/expected values. SNPs were annotated using SeqMonk (v1.48.0) (https://www.bioinformatics.babraham.ac.uk/projects/seqmonk/) and the most recent annotation reference (GCF_011100635.1). Using this approach we identified 2871, 2562 and 1826 potential monoallelic expressed SNPs in brain, liver and muscle, respectively. These SNPs were then overlapped with the ASM sites, using a +/− 1 Mb cut-off, which resulted in 477, 407 and 347 monoallelic expressed SNPs in brain, liver and muscle, respectively (corresponding to 91, 61 and 62 candidate genes, respectively). IGV (v2.12.3) was then used to visualise the SNPs in RNA-sequencing libraries of Sandy’s tissues and candidate genes were removed if biallelic expression of additional SNPs was observed108. This resulted in 13 genes (17 including known imprinted genes and genes characterised in this study), associated with 9 different ASM sites, containing SNPs displaying monoallelic expression that we classify as candidate imprinted genes (Supplementary Data 3).

DNA methylation analysis. For low-coverage bisulfite sequencing, raw reads were trimmed using Trim Galore! (v0.4.2; https://github.com/FelixKrueger/TrimGalore) by removing 10 bp from the 5’ end of all reads and removing low-quality base calls (Phred score < 20). Bismark (v0.19.0)109 with the -pbat option was used to map reads to the possum genome reference GCA_011100635.1. From the Bismark output report, global methylation in the CG dinucleotide context was calculated as the proportion of total methylated cytosines over total cytosines. The bisulfite non-conversion rate was estimated by calculating the proportion of non-CG methylation; by this measure, all libraries must have had a bisulfite conversion efficiency of at least 96.02% (Supplementary Data 5). PGC libraries producing < 400 CG calls or < 35% mapping efficiency were excluded from the analysis. Based on CpG methylation levels, PGC samples were categorised into three groups: ‘pre’ methylation loss (< 10 dpp; first blue dot, Fig. 6b; n=1), ‘demethylated’ (13–16 dpp; orange dots, Fig. 6b; n=3) and ‘remethylated’ (> 16 dpp; subsequent blue dots, Fig. 6b; n=12). A two-sample t-test with unequal variance was performed on CpG methylation values for n=15 individual samples (i.e., demethylated v remethylated).

For high-coverage whole-genome bisulfite sequencing, raw reads were trimmed using Trim Galore! (v0.6.7; https://github.com/FelixKrueger/TrimGalore) by removing adapters, 10 bp from the 5’ end of all reads and removing low-quality base calls (Phred score < 20). Trimmed reads were mapped to the possum genome reference GCA_011100635.1, deduplicated, and methylation calls in the CG dinucleotide context were extracted using Bismark v0.22.3109 with the -pbat option. The resulting BAM files were imported into SeqMonk (v1.48.0; https://www.bioinformatics.babraham.ac.uk/projects/seqmonk/), 13 dpp and 16 dpp germ cell samples were grouped into a replicate set, and probes were made over features for ASM sites (Supplementary Data 3) and H19 (Fig. S5b) with removal of duplicates. The feature probes were then quantitated using the Difference Quantitation method with minimum read count=1, and % methylation data was extracted from the annotated probe report.

Amplicon sequencing analysis. Sequencing reads were imported into Geneious v2022.2.2 (https://www.geneious.com) and trimmed using the BBDuk plugin. Trimmed reads were then mapped to the ‘in silico’ PCR product, generated by extracting the corresponding sequence from the possum genome reference GCA_011100635.1; see Supplementary Table 5 for genomic coordinates for amplicons). The ‘Find Variations/SNPs’ feature within Geneious was then used to determine the nucleotide sequence(s) at the SNP site of interest (Supplementary Table 5), with a minimum variant frequency of 0.2. The coverage (total number of reads), reference raw frequency and variant raw frequency were extracted for each SNP site per amplicon. This information was used to calculate the variant frequency for each sample (Supplementary Data 4).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Bond et al. have produced a near complete brushtail possum genome for the first time. They have\nused RNAseq data from over 100 samples to study the transition from pouch young to\nindependence, highlighting novel communication genes. They have also used their dataset to\ninvestigate genomic imprinting in brushtail possum for the first time, identifying four genes that\nhave not been reported in any other species.\nTheir RNAseq data has also enabled them to perform a lineage tracing analysis on possums in\nsouthern New Zealand, showing them to be an admix of distinct populations from Tasmania and\nmainland Australia.\nThis is an interesting study that sheds further light on the evolution of genomic imprinting. It\nhighlights the need for de novo analysis to identify parent-of-origin specific gene expression in\nmore distantly related species, rather than just looking at genes known to be imprinted in\neutherians.\nComments to the author.\nThis is an interesting paper, but it lacks cohesion. The story jumps around with no flow. The\nsection on lineage tracing whilst interesting feels out of place (especially in the middle of the\npaper) and should be removed and form a starting point for a population study paper.\nThe section on transition to life out of the pouch is compelling . However, it was not immediately\nclear that the data plotted on Figure 2a were average RPKM for the different paralogues. This\nassumes that all the paralogues are regulated in the same way. The authors should instead plot\neach gene separately in supplementary material to show each behaves the same and make it clear\nthat this figure shows the paralogue average. They should also plot the galactase/CYP expression\nratio as a graph if they are going to correlate other genes' expression with it.\n\"H19 is thought to orchestrate paternal-specific expression of\ninsulin-like growth factor 2 (IGF2) in mice\" the reference given for this statement only looks at\nhuman data. Leighton PA, et al. 1995. Nature 375:34–39 should be used instead.\nThe data on the changes of H19, Vulpeculin and ganderin expression between PY and adult is\nnoteworthy but I would like to have seen the data for other genes that follow the same signature.\nIf the authors perform a differential analysis between PY and adults they would be able to identify\nthe genes involved in this critical transition period. They could then cluster the genes to give a\nmore global view of what is going on.\nThe imprinting analysis is the most robust part of this study. Using differential methylation as a\nstarting point is sensible but I think they should also run an allele specific expression analysis in\ntandem. This would allow them to identify other candidate imprinted genes. Imprinting control\nregions in eutherians can influence gene expression over megabases and the cutoff of 20kb they\nuses would miss many of the candidates.\nThey should also comment on why they did not identify the H19-DMD in this analysis.\nThey have discovered two putative novel maternally methylated DMRs that are located over\npromoters of three paternally expressed genes. The majority of eutherian ICRs (23) share this\npattern with methylation on the maternal promoter of a paternally expressed gene. What is more\nnoteworthy is the discovery of a putative paternally methylated intronic DMR. Only three\npaternally methylated ICR are know in eutherians and each of those are intergenic. The authors\nshould comment on this in the discussion.", "reviewer_2": "In the manuscript entitled \"The admixed genome of brushtail possum reveals invasion history in\nNew Zealand and novel marsupial imprinted genes\", the first near-complete sequence of a\nbrushtail possum genome was produced using nanopore sequencing technology. Combined with\nmultiple transcriptome data, the authors revealed following three main findings. (i) Some\nmetabolic and chemo-sensory genes such as GAL, CYP, Vulpeculin and Ganderin showed striking\ngene expression change between pre- and post-weaning stages. (ii) Southern New Zealand\npossums showed high levels of genetic diversity driven by admixture between diverse sources in\nTasmania, NSW and unidentified source. (iii) Four genes, GPX7, EPM2AIP1, MLH1 and UBP1 were\nappeared to be imprinted in brushtail possums and three potential imprinted control regions (one\npaternal and two maternal differentially methylated regions, DMRs) were found from these gene\nloci.\nIdentification of novel marsupial imprinted genes and DMRs by unbiased method is important to\nadvance understanding how differently imprinting has been evolved in another mammalian lineage\nthan eutherians. It is also notable that the parental-conflict hypothesis seems to be applied for\nGPX7, EPM2AIP1 and UBP1. These findings would be a significant impact on the research field of\ngenomic imprinting because most works on marsupial genomic imprinting have been restricted to\northologues of eutherian imprinted genes for a long time.\nThe text in the manuscript is clear and easy to read and the data are high quality. There is no\ncritical flaw for the data interpretation and conclusions. Below are my suggestions that could help\nstrengthen the work.\n#1 In figure legend for Fig.1, (b) is appeared twice. The second (b) should be (c) and c) should be\n(d).\n#2 In Fig.1d, there is a description for bisulfite-sequencing of 13 individuals, 2 tissues, but it is\nunclear for me where these data were used in this research.\n#3 For Fig.S2a, there is a description in the text “we saw high expression in possum pouch young,\nbut not juvenile or adults”. I cannot agree with this description because I see that AFP expression\nlevels of most juveniles are still high in Fig.S2a.\n#4 For Fig.2a, the text says that GAL were high in pouch young but low in adults and juveniles.\nBut most dots in juvenile stage show similar expression levels as pouch young stage in Fig.2a.\nAlso, the text says that CYP became highly expressed in backriders and adults. In the graph, the\ndata from adult certainly show higher expression of CYP, but most data in juvenile stage show\nsimilar low expression levels as seen in pouch young stage.\n#5 Also for Fig.2c, the text says that IGF2 was biallelically expressed in adults and juveniles.\nIn the graph, I see that dots for adult samples are certainly close to the 50% line, but the data for\njuveniles are similar as that of pouch young stage.\n#6 Fig.2c shows allelic ratio of IGF2 expression and it becomes biallelically expressed in adult\nstage. How about the change of IGF2 expression level? Is expression level increased upon the gain\nof biallelic expression?\n#7 For the novel 4 marsupial imprinted genes and 3 DMRs, imprinted expression and differential\nmethylation are marsupial specific? I suggest confirming it because these data can be easily\nobtained from public transcriptome and DNA methylation data in mouse and human. For\nmonotremes, it might be difficult to determine their allelic expression, but checking methylation\nstatus of the orthologous CpG islands for the 3 DMRs will not be difficult if genomic DNA of a\nmonotreme species is available. There is a hypothesis that emergence of DMRs is associated with\nTE insertions. Is there any marsupial-specific TE insertion event nearby the 3 DMRs which might\nbe associated with the acquisition of differential methylation?\n#8 There are allelic expression data of brain in Fig.5b, c, d, but only not in Fig.5a. Is this just a\nmistake? If there is any reason, please provide an explanation.\n#9 In figure legend for Fig.6, “c)” should be “(c)”.\n#10 In Fig.6c, these 6 DMRs are hypomethylated in PGCs, nearly 50% methylated in somatic cells\nand hypo- or fully methylated in adult sperm. Is this methylation pattern characteristic only for\nDMRs associated with imprinting? I suggest checking methylation levels of the 173 ASM sites\n(remaining 167 sites?) in PGCs, somatic cells and adult sperm as the authors could possibly detect\nmore DMRs associated with imprinting from the remaining ASM sites.\n#11 In the third paragraph in the section entitled “Global methylation erasure...”, there is a\nsentence “however, H19 and GPX7, were fully methylated (94.24-95.79%) (Fig. 5c)”. As there is\nno methylation data in Fig.5c, (Fig. 5c) must be (Fig. 6c). In addition, there is a description\n“Supplementary Fig. 4a-f;”, but I cannot find Fig.S4f.\n#12 In discussion, the authors state that MLH1 does not have an obvious link to the parental\nconflict hypothesis. For the reason of this, it can be thought that imprinting of MLH1 is by-product\nof EPM2AIP1 imprinting. Because these two genes share promoter region, selective pressure for\nthe acquisition of EPM2AIP1 imprinting might affect MLH1 expression coincidentally.\n#13 In discussion, there is a sentence “Our analysis is almost certainly conservative – we have\nonly tested 4 of the 173 sites with ASM...”. It is confusing because I believe that the authors\ntested 173 sites and found 3 novel DMRs.\n#14 Previous studies have found out the tendency that paternal DMRs are located at non-\npromoter regions and have relatively lower CpG density while maternal DMRs are located at\npromoter regions and have higher CpG density. It is interesting that this tendency is completely\napplied to the three DMRs discovered in this study. How about mentioning this somewhere in\ndiscussion?", "reviewer_3": "In this manuscript, Bond et al. sequenced and assembled the brushtail possum (Trichosurus\nvulpecula) genome. The authors combined PacBio long-read data, 10x Genomics linked-read data,\nas well as Hi-C and bionano data to achieve super high quality and continuity. The authors\nanalyzed transcriptome data across developmental stages and identified gene expression signature\nof weaning in brushtail possum. Through SNPs called from RNAseq data, they are able to show the\ngenetic admixture of Dunedin possums. The authors also searched for novel imprinted genes in\nheterozygous individuals using allele-specific methylation profiling, and allele-specific gene\nexpression in the RNA-seq data.\nMajor comments:\nPage 2, the authors mentioned “however, unbiased searches for marsupial imprinted genes are\nlimited.” It might be true at the time of submission, but a recent survey of marsupial imprinted\ngenes in Monodelphis domestica was published https://doi.org/10.1093/molbev/msad022. The\nauthors may want to compare their results with Monodelphis, or at least cite the paper.\nThe evidence for monoallelic expression for EPM2AIP1 and MLH1 is strong. However, variability\nwas observed for GPX7 and UBP1, with some individuals showing various degree of biallelic\nexpression (120320_MR_SDR1, 030620_PY01, 280720_PY01, 200421_MB_SDR1,\n250521_MG_PF1). The authors need to discuss this in the manuscript.\nMonoallelic expression can be due to genomic imprinting, random monoallelic expression, or\nstrong cis-eQTL effects (https://doi.org/10.1038/hdy.2014.18). To exclude the possibility of\nrandom monoallelic expression and confirm imprinting, ideally, the parental transmission direction\ncan be validated in reference-allele-expressing individuals and alternative-allele-expressing\nindividuals. For IGR2R, GPX7, and EPM2AIP1, the authors were able to confirm the parental\ntransmission direction for samples expressing the reference and the alternative alleles, providing\nstrong evidence of parent-of-origin effect. However, for UBP1 and MLH1, only one direction was\nconfirmed and the a large proportion (>90%) of the informative individuvals lack parent-of-origin\ninformation. I would call these two genes (UBP1 and MLH1) as candidate imprinted genes because\nof the possibility of random monoallelic expression.\nThe authors need to be cautious to claim that the study is an unbiased genome-wide search of\nimprinted genes. “Our discovery pipeline relied on first identifying allele specific methylation”\nimplies a hypothesis that all imprinted genes are associated a DMR (differentially methylated\nregion). However, in eutherian mammals, only a third of the imprinted genes have a DMR.\nTherefore, discovery based on ASM will miss the imprinted genes regulated by other epigenetic\nmechanisms, such as histone modifications.\nSpecific comments:\n1) I would appreciate if the authors can add line numbers in the manuscript.\n2) Page 5, “chromosome assembly”, this assembly has very high continuity and completeness. The\nauthors should include the genome assessment stats (contig N50, scaffold N50 etc) in the results\nto highlight the strength of the manuscript.\n3) ASM discovery: what proportion of the genome was covered by >10x depth of Nanopore reads?\n4) Figure 5, please add the total RNA-seq read depth to each SNP/tissue.\n5) Page 12, global methylation dynamics: bisulfite sequencing may subject to low library\ncomplexity tissue, especially when the coverage is low. Please state the number of replicates\nclearly in the methods or results. Also, if the authors would like to claim a dramatic demethylation\nand remethylation, some sort of formal statistical analysis is needed. Simply comparing the overall\nmethylation percentages is not sufficient."}, "author_response": "RESPONSE TO REVIEWERS’ COMMENTS\nReviewer #1 (Remarks to the Author):\nBond et al. have produced a near complete brushtail possum genome for the first time. They have\nused RNAseq data from over 100 samples to study the transition from pouch young to\nindependence, highlighting novel communication genes. They have also used their dataset to\ninvestigate genomic imprinting in brushtail possum for the first time, identifying four genes that\nhave not been reported in any other species.\nTheir RNAseq data has also enabled them to perform a lineage tracing analysis on possums in\nsouthern New Zealand, showing them to be an admix of distinct populations from Tasmania and\nmainland Australia.\nThis is an interesting study that sheds further light on the evolution of genomic imprinting. It\nhighlights the need for de novo analysis to identify parent-of-origin specific gene expression in more\ndistantly related species, rather than just looking at genes known to be imprinted in eutherians.\nWe thank reviewer 1 for favourably identifying the novel aspects of our work, and its range of\nimpact from population genetics to developmental biology and epigenetics. We also agree that our\nde novo search for imprinted genes is both timely and interesting, as well as providing much needed\ninsight into the evolution of imprinting.\nComments to the author.\nThis is an interesting paper, but it lacks cohesion. The story jumps around with no flow. The section\non lineage tracing whilst interesting feels out of place (especially in the middle of the paper) and\nshould be removed and form a starting point for a population study paper.\nFollowing the direction from the editor, we have retained the lineage tracing section in this\nmanuscript. While lineage tracing is an important part of the population genetics story, in our\nopinion, the admixture aspect feeds nicely into genomic imprinting studies, because admixture\nprovides greater imprinted gene search resolution.\nThe section on transition to life out of the pouch is compelling . However, it was not immediately\nclear that the data plotted on Figure 2a were average RPKM for the different paralogues. This\nassumes that all the paralogues are regulated in the same way. The authors should instead plot each\ngene separately in supplementary material to show each behaves the same and make it clear that\nthis figure shows the paralogue average. They should also plot the galactase/CYP expression ratio as\na graph if they are going to correlate other genes' expression with it.\nTo address this point we have provided the GAL/CYP ratio as a line plot in the main text (Fig. 2a), and\nmade our axis labels and figure legend clearer regarding the nature of this ratio (i.e. the fact it is an\naverage of multiple paralogues rather than representing individual genes).\nIn addition, the expression levels (RPKM - log2 values) for each GAL and CYP paralogue used for this\nanalysis is given in Supplementary Table 3. We have also added plots with the expression profiles for\nthese individual paralogues to this file.\n\"H19 is thought to orchestrate paternal-specific expression of insulin-like growth factor 2 (IGF2) in\nmice\" the reference given for this statement only looks at human data. Leighton PA, et al. 1995.\nNature 375:34–39 should be used instead.\nWe have now corrected this error.\nThe data on the changes of H19, Vulpeculin and ganderin expression between PY and adult is\nnoteworthy but I would like to have seen the data for other genes that follow the same signature. If\nthe authors perform a differential analysis between PY and adults they would be able to identify the\ngenes involved in this critical transition period. They could then cluster the genes to give a more\nglobal view of what is going on.\nTo be clear about our analysis, we initially found the GAL/CYP ratio by comparing PY and adults, and\nthen we searched for other genes that showed co-regulation with this change in liver metabolism at\nweaning. In doing so, we uncovered H19, vulpeculin and ganderin as these show the greatest\ncorrelation (or anti-correlation) to the GAL/CYP ratio - i.e., the metabolic weaning ‘signature’ we\nuncovered. We have since updated the manuscript to be clearer about how we found the GAL/CYP\nratio (i.e. differential analysis between PY and adults):\n“Perhaps the most striking gene expression change we observed in liver between pouch young and\nadults was associated with the metabolic effects of weaning (Fig. 2a).”\nNevertheless, we now realise that readers could be interested in weaning associated genes in\naddition to those we have discussed. For this reason we have supplied the expression levels (RPKM -\nlog2 values) for all genes, and sorted them based on their correlation to the GAL/GYP ratio. This data\nis supplied in Supplementary Table 3.\nThe imprinting analysis is the most robust part of this study. Using differential methylation as a\nstarting point is sensible but I think they should also run an allele specific expression analysis in\ntandem. This would allow them to identify other candidate imprinted genes. Imprinting control\nregions in eutherians can influence gene expression over megabases and the cutoff of 20kb they\nuses would miss many of the candidates.\nWhile this was not detailed in the original manuscript, previously we have searched for imprinted\ngenes in possums using only allele-specific expression analysis (Reese, 2021, Thesis). We were not\nsatisfied with this approach as most candidate genes were false positives, often caused by non-\nexpressed pseudogene DNA sequences providing the false impression of an imprinted expression in\nits paralogue. Our strike-rate identifying imprinted genes improved greatly once we used allele\nspecific methylation identified via nanopore sequencing. This was the primary mechanism to identify\ncadidate imprinted regions as we searched for mono-allelelic expression nearby.\nIn the initial submission, we did this in a rather ad hoc manual manner (focusing on genes next to\nASM sites), however, during revision we performed a search for genome-wide mono-allelic\nexpression of SNPs in Sandy’s brain, liver and muscle, and then overlapped this with our ASMs using\na +/- 1-Mbp cut off - the average size of the largest imprinted domains in mammals (MacDonald and\nMann, 2020, PMID: 32760061; Hubert and Demars, 2022, PMID: 35368671). This resulted in a list of\n91, 61 and 62 genes with potential mono-allelic expression in RNA-sequencing data. The list was\nthen filtered based on evidence of biallelic expression of additional SNPs in the genes of interest\nand/or removal of false positives. As a final ‘sanity-check’ for the candidate imprinted regions\npassing these filters, the associated ASMs were examined in germ cells undergoing epigenetic\nreprogramming. Of these, we found 8 ASMs showing the expected methylation pattern for\nimprinting (unmethylated in germ cells, around ~50% in somatic cells, either ~100% or ~0%\nmethylation in adult sperm) (Fig 6d, lower panel).\nThese high confidence candidates will ultimately need comprehensive verification in subsequent\nstudies, but provides a useful starting point for further analysis and the opportunity to compare with\nother studies (e.g. opossum - see below).\nThey should also comment on why they did not identify the H19-DMD in this analysis.\nThe animal that formed the basis of our genome sequence (Sandy) was homozygous at the H19\nlocus, hence it was not detected in the ASM analysis (which relies on hetorzygous SNPs to generate\nthe methylation haplotypes, or allele specific methylation patterns). We have now added an\nexplanation of this in the discussion of our manuscript.\nThey have discovered two putative novel maternally methylated DMRs that are located over\npromoters of three paternally expressed genes. The majority of eutherian ICRs (23) share this\npattern with methylation on the maternal promoter of a paternally expressed gene. What is more\nnoteworthy is the discovery of a putative paternally methylated intronic DMR. Only three paternally\nmethylated ICR are know in eutherians and each of those are intergenic. The authors should\ncomment on this in the discussion.\nReviewer 1 makes an interesting observation here that we had overlooked. This has now been\nadded to the discussion following their suggestion:\n“Associated with these genes was parent-specific methylation - two instances where the maternal\nchromosome was methylated (like most imprinted loci), and one example of paternal methylation\n(GPX7). Interestingly, this site of differential methylation is located in an intron, similar to other known\npaternally methylated imprint control regions, and unlike maternal imprint control regions which are\ngenerally located at transcription start sites70.”\nReviewer #2 (Remarks to the Author):\nIn the manuscript entitled \"The admixed genome of brushtail possum reveals invasion history in New\nZealand and novel marsupial imprinted genes\", the first near-complete sequence of a brushtail\npossum genome was produced using nanopore sequencing technology. Combined with multiple\ntranscriptome data, the authors revealed following three main findings. (i) Some metabolic and\nchemo-sensory genes such as GAL, CYP, Vulpeculin and Ganderin showed striking gene expression\nchange between pre- and post-weaning stages. (ii) Southern New Zealand possums showed high\nlevels of genetic diversity driven by admixture between diverse sources in Tasmania, NSW and\nunidentified source. (iii) Four genes, GPX7, EPM2AIP1, MLH1 and UBP1 were appeared to be\nimprinted in brushtail possums and three potential imprinted control regions (one paternal and two\nmaternal differentially methylated regions, DMRs) were found from these gene loci.\nIdentification of novel marsupial imprinted genes and DMRs by unbiased method is important to\nadvance understanding how differently imprinting has been evolved in another mammalian lineage\nthan eutherians. It is also notable that the parental-conflict hypothesis seems to be applied for\nGPX7, EPM2AIP1 and UBP1. These findings would be a significant impact on the research field of\ngenomic imprinting because most works on marsupial genomic imprinting have been restricted to\northologues of eutherian imprinted genes for a long time.\nThe text in the manuscript is clear and easy to read and the data are high quality. There is no critical\nflaw for the data interpretation and conclusions. Below are my suggestions that could help\nstrengthen the work.\n#1 In figure legend for Fig.1, (b) is appeared twice. The second (b) should be (c) and c) should be (d).\nWe have now corrected this error.\n#2 In Fig.1d, there is a description for bisulfite-sequencing of 13 individuals, 2 tissues, but it is\nunclear for me where these data were used in this research.\nThe bisulfite sequencing data was used in the epigenetic reprogramming studies. We have now\nclarified this in the Fig. 1 legend as follows:\n“A range of functional genomics data supports the assembly, including: RNA-sequencing (mostly\nliver) from a wide population of possums and 13 tissues from Sandy. Also performed was DNA\nmethylation profiling, featuring nanopore and bisulfite sequencing for genomic imprinting and\ngermline reprogramming analysis, respectively, and a sequence assembly of a Tasmanian possum.”\n#3 For Fig.S2a, there is a description in the text “we saw high expression in possum pouch young, but\nnot juvenile or adults”. I cannot agree with this description because I see that AFP expression levels\nof most juveniles are still high in Fig.S2a.\nThis was an oversight on our behalf. We were not clear that juveniles are actually pouch young >120\ndays old (for which we were not able to accurately date). As such we have corrected this sentence to\nthe following:\n“Yet we saw high expression in possum pouch young, including juvenile pouch young >120 d post birth\n(Fig. S2a).”\nWe have also updated our methods to be clearer about how developmental staging was undertaken,\nand what the terms ‘juvenile’ and ‘backrider’ actually mean:\n“Pouch young >120 days of development were challenging to stage, and for this reason have been\nreferred to as ‘juveniles’ (i.e. pouch-young greater than >120 d, but younger than ~180 days, upon\nwhich time they leave the pouch). The sole ‘back-rider’ was Sheila, who was captured with her mother\n(Puku) soon after having left the pouch.”\n#4 For Fig.2a, the text says that GAL were high in pouch young but low in adults and juveniles. But\nmost dots in juvenile stage show similar expression levels as pouch young stage in Fig.2a. Also, the\ntext says that CYP became highly expressed in backriders and adults. In the graph, the data from\nadult certainly show higher expression of CYP, but most data in juvenile stage show similar low\nexpression levels as seen in pouch young stage.\nWe have corrected this sentence to be congruent with the pouch-young, juvenile and adult\ndefinitions as for the comment above:\n“Genes associated with the Leloir pathway, involved with the catabolism of galactose and other\ncarbohydrates, were high in pouch young, but low in adults (Supplementary Table 3).”\n#5 Also for Fig.2c, the text says that IGF2 was biallelically expressed in adults and juveniles. In the\ngraph, I see that dots for adult samples are certainly close to the 50% line, but the data for juveniles\nare similar as that of pouch young stage.\nWe have corrected this sentence to be congruent with the pouch-young, juvenile and adult\ndefinitions as for the comment above:\n“We found exclusive maternal expression of H19, whereas IGF2 was paternally expressed in pouch\nyoung and biallelically expressed in adults where H19 expression was silenced (Fig. 2b,c;\nSupplementary Fig. S2d; Supplementary Table 4).”\n#6 Fig.2c shows allelic ratio of IGF2 expression and it becomes biallelically expressed in adult stage.\nHow about the change of IGF2 expression level? Is expression level increased upon the gain of\nbiallelic expression?\nWe found that IGF2 expression decreased during development/aging (see figure below), a result\nconsistent with our understanding of the insulin like growth factors in human - IGF2 is\npredominantly a fetal growth factor that gets silenced during early postnatal life in humans, with its\nexpression replaced by IGF1.\nIGF2 in possum\nIGF1 in possum\nWe feel this result does not justify a new figure in the manuscript, however, we could provide one if\nthe reviewer/editor believes we should.\n#7 For the novel 4 marsupial imprinted genes and 3 DMRs, imprinted expression and differential\nmethylation are marsupial specific? I suggest confirming it because these data can be easily\nobtained from public transcriptome and DNA methylation data in mouse and human. For\nmonotremes, it might be difficult to determine their allelic expression, but checking methylation\nstatus of the orthologous CpG islands for the 3 DMRs will not be difficult if genomic DNA of a\nmonotreme species is available. There is a hypothesis that emergence of DMRs is associated with\nTE insertions. Is there any marsupial-specific TE insertion event nearby the 3 DMRs which might\nbe associated with the acquisition of differential methylation?\nThis is an interesting question which we are actively trying to address. MLH1/EPM2AIP1 expression\nand methylation has been extensively studied in humans (Lynch et al., 2015), and is not imprinted.\nNevertheless, there are indications from cattle that parent-of-origin effects might exist at this locus\n(Kenny et al., 2022; PMID: 35910233), implying that at least imprinting of MLH1/EPM2AIP1 is not\nmarsupial specific. Experiments are underway to definitively test MLH1/EPM2AIP1 imprinting (and\nthat of UBP1 and GPX7) in other eutherian, marsupial and monotreme species using tissue samples\ncollected for this purpose, as well as publicly available datasets. While this work is exciting, it is\nongoing, and we believe beyond the scope of this manuscript.\n#8 There are allelic expression data of brain in Fig.5b, c, d, but only not in Fig.5a. Is this just a\nmistake? If there is any reason, please provide an explanation.\nFigure 5a shows mono-allelic expression of GPX7 in ‘Sandy’ the individual sequenced in the genome.\nThe number of RNA-seq reads for brain is <5, so it is not presented (as stated in the figure legend).\n#9 In figure legend for Fig.6, “c)” should be “(c)”.\nThis has been corrected.\n#10 In Fig.6c, these 6 DMRs are hypomethylated in PGCs, nearly 50% methylated in somatic cells and\nhypo- or fully methylated in adult sperm. Is this methylation pattern characteristic only for DMRs\nassociated with imprinting? I suggest checking methylation levels of the 173 ASM sites\n(remaining 167 sites?) in PGCs, somatic cells and adult sperm as the authors could possibly detect\nmore DMRs associated with imprinting from the remaining ASM sites.\nWe have looked at the 173 ASM sites as suggested. We have added a heatmap to Fig. 6 (i.e. Fig. 6d -\nshown above. The top panel shows methylation in PGCs, sperm and somatic cells for all ASM sites\n(excluding those with missing data, n=32) sorted by somatic methylation. We identified 41 sites that\nhave somatic methylation between 30 and 57% (i.e. the approximate range of somatic methylation\nexhibited by sites with confirmed mono-allelic expression in Fig. 6c). Of these 41 sites, 32 (including\nthe 6 ASMs presented in Fig. 6c) display methylation for PGCs and sperm within acceptable ranges to\nsuggest imprinting (<15% for PGCs, <5% or >94% for sperm; see Supplementary Table 6). The bottom\npanel of Fig. 6D shows the ASM sites that are associated with genes that have confirmed imprinted,\nmono-allelic expression or are candidate imprinted genes (n = 13).\n#11 In the third paragraph in the section entitled “Global methylation erasure...”, there is a\nsentence “however, H19 and GPX7, were fully methylated (94.24-95.79%) (Fig. 5c)”. As there is\nno methylation data in Fig.5c, (Fig. 5c) must be (Fig. 6c). In addition, there is a description\n“Supplementary Fig. 4a-f;”, but I cannot find Fig.S4f.\nThese were typographical oversights which we regret, and have now been corrected.\n#12 In discussion, the authors state that MLH1 does not have an obvious link to the parental conflict\nhypothesis. For the reason of this, it can be thought that imprinting of MLH1 is by-product of\nEPM2AIP1 imprinting. Because these two genes share promoter region, selective pressure for the\nacquisition of EPM2AIP1 imprinting might affect MLH1 expression coincidentally.\nThis possibility that MLH1 is a bystander of EPM2AIP1 imprinting is a pertinent point, which we have\nnow made clearer in the discussion:\n“The MLH1 gene does not have an obvious link to the parental conflict hypothesis, and therefore may\nbe a ‘bystander’ caught in the imprinted regulation of closely linked EPM2AIP1.”\n#13 In discussion, there is a sentence “Our analysis is almost certainly conservative – we have\nonly tested 4 of the 173 sites with ASM...”. It is confusing because I believe that the authors\ntested 173 sites and found 3 novel DMRs.\nWe apologise for not being clearer here - of the 173 sites identified as having allele-specific\nmethylation (i.e. ‘ASMs’), we have only comprehensively tested 3 sites for parent-specific expression\nin neighbouring genes. This is detailed in the revised manuscript as follows:\n“It is likely more imprinted genes in possum will be uncovered - of the 173 ASM sites identified by\nnanopore sequencing, we have performed experiments validating parent-specific expression on genes\nassociated with just 3 sites.”\n#14 Previous studies have found out the tendency that paternal DMRs are located at non-promoter\nregions and have relatively lower CpG density while maternal DMRs are located at promoter regions\nand have higher CpG density. It is interesting that this tendency is completely applied to the three\nDMRs discovered in this study. How about mentioning this somewhere in discussion?\nAs for Reviewer 1, Reviewer 2 makes an interesting observation here that we had overlooked. This\nhas now been added to the discussion following their suggestion:\n“Associated with these genes was parent-specific methylation - two instances where the maternal\nchromosome was methylated (like most imprinted loci), and one example of paternal methylation\n(GPX7). Interestingly, this site of differential methylation is located in an intron, similar to other known\npaternally methylated imprint control regions, and unlike maternal imprint control regions which are\ngenerally located at transcription start sites70.”\nReviewer #3 (Remarks to the Author):\nIn this manuscript, Bond et al. sequenced and assembled the brushtail possum (Trichosurus\nvulpecula) genome. The authors combined PacBio long-read data, 10x Genomics linked-read data, as\nwell as Hi-C and bionano data to achieve super high quality and continuity. The authors analyzed\ntranscriptome data across developmental stages and identified gene expression signature of\nweaning in brushtail possum. Through SNPs called from RNAseq data, they are able to show the\ngenetic admixture of Dunedin possums. The authors also searched for novel imprinted genes in\nheterozygous individuals using allele-specific methylation profiling, and allele-specific gene\nexpression in the RNA-seq data.\nMajor comments:\nPage 2, the authors mentioned “however, unbiased searches for marsupial imprinted genes are\nlimited.” It might be true at the time of submission, but a recent survey of marsupial imprinted\ngenes in Monodelphis domestica was published https://doi.org/10.1093/molbev/msad022. The\nauthors may want to compare their results with Monodelphis, or at least cite the paper.\nThe reviewer is correct - the mentioned Monodelphis domestica paper came out during the\nmanuscript review process. We have now cited this work in the introduction and results section (Fig.\n6d), and also described some of their results in our discussion:\n“In addition, a recently-published analysis of RNA-sequencing in Monodelphis domestica estimates\nthere may be up to 60 imprinted genes in the Amercian opossum17.”\nThe evidence for monoallelic expression for EPM2AIP1 and MLH1 is strong. However, variability was\nobserved for GPX7 and UBP1, with some individuals showing various degree of biallelic expression\n(120320_MR_SDR1, 030620_PY01, 280720_PY01, 200421_MB_SDR1, 250521_MG_PF1). The\nauthors need to discuss this in the manuscript.\nWe have changed our description of UBP1 and GPX7 expression to include ‘bias’ rather than strict\nmono-allelic expression. For example, the results now reads:\nWhile monoallelic expression was clear for MLH1/EPM2AIP1, there appeared to be biased expression\nof GPX7 and UBP1.\nAnd the discussion now reads:\nUsing long-read methylation sequencing as a starting point, we discovered 4 genes with parent-specific\nimprinted expression across a wide range of tissues, albeit with two genes (UBP1 and GPX7)\nshowing some evidence for leaky monoallelic expression, as previously reported for some marsupial\nimprinted genes72.\nMonoallelic expression can be due to genomic imprinting, random monoallelic expression, or strong\ncis-eQTL effects (https://doi.org/10.1038/hdy.2014.18). To exclude the possibility of random\nmonoallelic expression and confirm imprinting, ideally, the parental transmission direction can be\nvalidated in reference-allele-expressing individuals and alternative-allele-expressing individuals. For\nIGR2R, GPX7, and EPM2AIP1, the authors were able to confirm the parental transmission direction\nfor samples expressing the reference and the alternative alleles, providing strong evidence of\nparent-of-origin effect. However, for UBP1 and MLH1, only one direction was confirmed and the a\nlarge proportion (>90%) of the informative individuvals lack parent-of-origin information. I would call\nthese two genes (UBP1 and MLH1) as candidate imprinted genes because of the possibility of\nrandom monoallelic expression.\nOur population analysis from heterozygous individuals clearly shows expression is monoallelic, but\nthat both versions of the SNP can be expressed in the population. Thus, cis-eQTL can not be the\nexplanation for monoallelic expression for our 4 novel imprinted geens (MLH1, EPM2AIP1, GPX7 and\nUBP1).\nRandom monoallelic expression is common at certain gene clusters, like olfactory and immune\ngenes, or genes subject to transcriptional ‘bursting’ (Eckersley-Maslin, 2014, PMID: 24780084).\nImportantly, random monoallelic expression features the expression of alternative alleles within the\nsame individual and tissue (hence the term, ‘random’, and resulting in mosaic expression patterns,\nlike the calico cat). For it to be confused with imprinted gene expression like we saw (with\nmonoallelic expression in all tissues), the ‘random’ inactivating event would need to have occured in\nthe zygote or somehow affect the same allele in all of the progenitor cells giving rise to the embryo.\nThere is to our knowledge, no example of such a non-imprinted, random monoallelically expressed\ngene that affects all tissues of an individual. Hence we consider this possibility very unlikely.\nNevertheless, to bolster our data, we did identify an additional duo for MLH1 and UBP1, and saw the\nsame directionality of imprinting in both cases. This data has been added to Supplementary Table 4.\nThe authors need to be cautious to claim that the study is an unbiased genome-wide search of\nimprinted genes. “Our discovery pipeline relied on first identifying allele specific methylation”\nimplies a hypothesis that all imprinted genes are associated a DMR (differentially methylated\nregion). However, in eutherian mammals, only a third of the imprinted genes have a DMR.\nTherefore, discovery based on ASM will miss the imprinted genes regulated by other epigenetic\nmechanisms, such as histone modifications.\nWe respect the reviewer’s position on this point, and have now removed all references to an\n‘unbiased’ search for imprinted genes.\nSpecific comments:\n1) I would appreciate if the authors can add line numbers in the manuscript.\nThis has now been updated.\n2) Page 5, “chromosome assembly”, this assembly has very high continuity and completeness. The\nauthors should include the genome assessment stats (contig N50, scaffold N50 etc) in the results\nto highlight the strength of the manuscript.\nThis has now been updated in the manuscript as follows:\nThis resulting assembly had 99.88% of the sequence mapped to chromosomes, matching the\nkaryotype (Fig. 1c), and possessing strong assembly rates (scaffold N50, 442,560,073 bp; contig N50,\n4,314,688 bp).\n3) ASM discovery: what proportion of the genome was covered by >10x depth of Nanopore reads?\nWe used the tool bedcov from the package SAMtools (v1.16.1) to calculate the coverage at 5X and\n10X for the nanopore reads. Using this methodology we estimated a coverage of 98% and and\n80.1%, respectively. We included this information in the methods of the revised manuscript.\n4) Figure 5, please add the total RNA-seq read depth to each SNP/tissue.\nAs in the original manuscript, RNA-sequencing read depth for all SNPs analysed is given in\nSupplementary Table 4. To further help readers asssess the strength of data, we have added the\nRNA-sequencing read depth for the heterozyous individuals presented in Fig. 5, Fig. S2d and Fig. S4.\nWe have also explained that these numbers represent the read depth values (≥5) in the legends for\nthese figures.\n5) Page 12, global methylation dynamics: bisulfite sequencing may subject to low library complexity\ntissue, especially when the coverage is low. Please state the number of replicates clearly in the\nmethods or results. Also, if the authors would like to claim a dramatic demethylation and\nremethylation, some sort of formal statistical analysis is needed. Simply comparing the overall\nmethylation percentages is not sufficient.\nIn response to this comment, we have stated the replicates used more clearly in the methods\nsection (“DNA methylation analysis”). We agree that a statistical comparison is appropriate,\nhowever, since we have only one sample in our dataset representing PGCs before methylation\nreduction (see first blue dot, Fig. 6b), we could not include this time point in a statistical test and\nhave removed the word ‘dramatically’ associated with demethylation in the results. We instead\nperformed a t-test between the average CpG methylation in the ‘demethylated’ samples (orange\ndots, Fig. 6b; n=3) vs ‘remethylated’ samples (subsequent blue dots, Fig. 6b; n= 12) and found that\nmethylation differed significantly (p = 0.00002906), as expected. We have included details of this\nstatistical comparison in the methods and results sections.\nREVIEWERS' COMMENTS\nReviewer #1 (Remarks to the Author):\nThis is my second time in reviewing this article. The authors have addressed most of my\ncomments after the first review.\nI am pleased that the authors are now able to include a comparison of their findings with those\nfrom Cao et al. Of particular interest is their finding that 4 of their novel ASMs are located near\ngenes\nCao et al identified to have imprinted expression in opossum. This suggests that marsupial may\nhave their own set of imprinted domains. This will be of interest to others in the marsupial\nimprinting field.\nReviewer #2 (Remarks to the Author):\nThe manuscript has been revised well and the responses to my comments and suggestions are\nacceptable.\nRESPONSE TO REVIEWERS’ COMMENTS\nReviewer #1 (Remarks to the Author):\nThis is my second fime in reviewing this arficle. The authors have addressed most of my comments\nafter the first review.\nI am pleased that the authors are now able to include a comparison of their findings with those from\nCao et al. Of parficular interest is their finding that 4 of their novel ASMs are located near genes\nCao et al idenfified to have imprinted expression in opossum. This suggests that marsupial may have\ntheir own set of imprinted domains. This will be of interest to others in the marsupial imprinfing\nfield.\nWe are also pleased with the outcome of our comparison with Cao et al. and agree that it will be of\ninterest to others in the marsupial imprinfing field and those interested in the evolufion of this\nimportant mechanism. Thank you for your comments and taking the fime to review our manuscript.\nReviewer #2 (Remarks to the Author):\nThe manuscript has been revised well and the responses to my comments and suggesfions are\nacceptable.\nThank you for reviewing our manuscript. We are pleased that our responses to your comments are\nsuggesfions are acceptable."}]

Decision letter