A chromosome-level genome of Antechinus flavipes provides a reference for an Australian marsupial genus with male death after mating.

The 15 species of small carnivorous marsupials that comprise the genus Antechinus exhibit semelparity, a rare life-history strategy in mammals where synchronized death occurs after one breeding season. Antechinus males, but not females, age rapidly (demonstrate organismal senescence) during the breeding season and show promise as new animal models of ageing. Some antechinus species are also threatened or endangered. Here, we report a chromosome-level genome of a male yellow-footed antechinus Antechinus flavipes. The genome assembly has a total length of 3.2 Gb with a contig N50 of 51.8 Mb and a scaffold N50 of 636.7 Mb. We anchored and oriented 99.7% of the assembly on seven pseudochromosomes and found that repetitive DNA sequences occupy 51.8% of the genome. Draft genome assemblies of three related species in the subfamily Phascogalinae, two additional antechinus species (Antechinus argentus and A. arktos) and the iteroparous sister species Murexia melanurus, were also generated. Preliminary demographic analysis supports the hypothesis that climate change during the Pleistocene isolated species in Phascogalinae and shaped their population size. A transcriptomic profile across the A. flavipes breeding season allowed us to identify genes associated with aspects of the male die-off. The chromosome-level A. flavipes genome provides a steppingstone to understanding an enigmatic life-history strategy and a resource to assist the conservation of antechinuses.

INTRODUCTION

For every animal species, the most critical activity is mating and the raising of offspring. Reproduction is a costly process, and animals invest energy into strategies that maximize the number and survival of offspring (Collett et al., 2018). Broadly, these strategies span the so‐called iteroparity–semelparity continuum: from a single (semel) to repeated (itero) reproductive episode(s) before death (reviewed in Hughes, 2017; Krajewski et al., 2000). In mammals, semelparity is observed in some Australian dasyurids (males only) and South American didelphids (both sexes) (Fisher et al., 2013; Leiner et al., 2008). Obligate male semelparity is the most extreme form of semelparity, occurring in the marsupial subfamilies Phascogalinae (15 Antechinus species and three Phascogale species; Bradley, 2003; Naylor et al., 2008) and Dasyurinae (the kaluta, Dasykaluta rosamondae; Hayes et al., 2019). In these dasyurid species, all males die soon after a 1‐ to 3‐week mating period (“die‐off” at 11.5 months of age), and a proportion of females survive to breed in the subsequent year (Baker & Dickman, 2018).

Breeding is a stressful life‐history event, but the stress response enables most animals to cope and to mobilize energy (Romero, 2004). Hormones play a vital role in this process. While in antechinus, both sexes have elevated stress hormone levels (chiefly the corticosteroid cortisol) during breeding, males also exhibit surging plasma testosterone levels (Bradley, 2003; Braithwaite & Lee, 1979; Croft, 2003; McAllan, 2006; Naylor et al., 2008). This interferes with negative feedback mechanisms and reduces levels of corticosteroid binding globulin (CBG), increasing the levels of free (i.e., unbound), active corticosteroids in the blood (Bradley, 2003; Naylor et al., 2008). As the breeding season progresses, corticosteroid levels continue to rise; this ultimately contributes to systemic collapse and death of all males. Corticosteroids are anti‐inflammatory and, at high concentrations, immunosuppressive, triggering a decline of body condition (Manoli et al., 2007; Perogamvros et al., 2012). Stress, combined with malnutrition, results in a significant negative energy balance and is associated with numerous comorbidities (Argiles et al., 2016). Cortisol is also elevated in females throughout the short breeding season, but in the absence of elevated testosterone, free cortisol levels remain relatively constant.

There is mounting evidence that testosterone and cortisol contribute to differences in ageing and related diseases (Austad & Fischer, 2016; Gems, 2014; Lemaitre et al., 2020; Martocchia et al., 2016). Moreover, ageing and obligate semelparity manifest as a decline in immune function—immunosenescence (Childs et al., 2017; Froy et al., 2019). In contrast to current mammal models of ageing, antechinuses are “natural” models of physiological senescence (see Diamond, 1982; McAllan, 2006). For example, the brain of the ageing wild‐type mouse and rat does not accumulate amyloid‐beta (Aβ) (Jucker, 2010). Aβ plaques, which are associated with chronic stress and elevated cortisol levels (e.g., see Lesuis et al., 2018), accumulate at the end of the antechinus lifespan in antechinuses—mirroring Alzheimer's disease neuropathology and potentially making them a suitable disease model (McAllan, 2006; Naylor et al., 2008). A high‐quality genome assembly is critical for future genetic experiments (e.g., CRISPR). Moreover, considering the strong effects of cortisol and testosterone on male obligate semelparity, identifying the expression of genes associated with these hormones and its consequences may help to understand the evolution of this extreme reproductive strategy.

Recent work has uncovered several new antechinus species, bringing the total number of species to 15 (Baker et al., 2013; Baker et al., 2014; Baker et al., 2015; Baker et al., 2012; Mutton et al., 2019). In 2018, Antechinus arktos and A. argentus were listed as Endangered on the Australian Government's threatened species list and were included in the 20 Australian mammals most likely to go extinct in the next two decades (Geyle et al., 2018). Moreover, areas of their habitat were recently ravaged by extensive wildfires and drought. It is estimated that nearly 3 billion animals died in these events in 2020 (Chris Dickman, 2020, University of Sydney, pers. comm.). Genomic resources are urgently required to support research on the population structure, speciation and recovery of antechinus species.

To investigate the phylogeny, demographic history, and genomic diversity of Antechinus flavipes and to pave the way to understanding the evolution of semelparity, high‐quality genomic data are essential. A genome assembly of the brown antechinus (Antechinus stuartii) was recently reported (Brandies et al., 2020a). However, a chromosome‐level genome is not yet available. PacBio long‐read sequencing combined with the chromosome conformation capture method Hi‐C resolves complex regions (such as repeat regions) and can greatly improve downstream analyses (Schloissnig et al., 2021). Here, we employed third‐generation sequencing technologies to generate the first chromosome‐level genome of an antechinus, that of the yellow‐footed antechinus (A. flavipes), as well as draft genomes of three related dasyurids (Figure 1a). We highlight genome features, phylogeny and historical demography, and a gene expression profile of A. flavipes male die‐off.

FIGURE 1

Overview of sequenced species and reference genome assemblies generated in this study. (a) Photographs of the five marsupial species sequenced. A yellow‐footed antechinus (Antechinus flavipes), silver‐headed antechinus (Antechinus argentus), and black‐tailed dusky antechinus (Antechinus arktos) from Australia; and a black‐tailed dasyure (Murexia melanurus) from Papua New Guinea. Antechinus spp. and M. melanurus photographs are courtesy of A.M.B. and John Hornbuckle, respectively. (b) Circos plot of the reference genome assembly of A. flavipes. The outermost segment represents chromosome sequences, with the numbers on the external surface indicating genome size in Mb. The line plots, from outside to inside, represent the distribution of CDS density (from 0 to 0.15), GC content (from 0.30 to 0.65) and TE ratio (from 0.2 to 1.0). Frequencies were calculated in 500‐kb sliding windows. (c) Circos plot showing shared synteny of A. flavipes (chr1‐chr7) and Sarcophilus harrisii (NC_045426.1‐NC_045432.1). Aligned using lastz. The synteny blocks are linked using lines coloured in accordance with the A. flavipes chromosomes. Aligned blocks with length shorter than 20 kb are not shown. Chr7 in A. flavipes corresponds to the X chromosome of Sarcophilus harrisii

MATERIALS AND METHODS

Sample collection

“AdamAnt,” an adult male yellow‐footed antechinus (Antechinus flavipes), was sampled in Samford Valley (greater Brisbane), QLD, Australia (27°22′S, 152°52′E) in September 2018. A second male individual (“Ant2”), sampled at the same location in 2017, was not used for genome assembly but employed in demographic analysis. A total of 106 A. flavipes samples were collected for transcriptome sequencing (RNA‐seq). Ear tissue samples (i.e., ear clippings) for short‐insert whole‐genome sequencing were obtained from the black‐tailed dusky antechinus Antechinus arktos (Lamington National Park, QLD, Australia, 28°22′S, 153°15′E; specimen AA100 sampled by A.M.B. in 2015) and the silver‐headed antechinus Antechinus argentus (from Blackdown Tableland National Park, QLD, Australia; 23°82′S, 149°07′E; specimen BD‐17‐5A sampled by A.M.B. in 2017). Frozen liver tissue from the black‐tailed dasyure (Murexia melanurus; Wau, Morobe, Papua New Guinea, 7°20′S, 146°42′E; voucher specimen ABTC46020 sampled by the South Australian Museum in 1985) was also sequenced. All DNA and RNA were extracted by BGI Hong Kong.

Reference genome sequencing and assembly

We used a combination of long‐fragment sequencing, short‐insert library sequencing for error correction and gap filling, and chromatin conformation capture (Hi‐C) to generate chromosome‐level semelparous mammal reference genomes. High‐molecular‐weight (HMW) DNA extracted from the testis of the A. flavipes individual “AdamAnt” was used to generate long‐read (PacBio) sequencing data by Annoroad Gene Technology. Briefly, 8 µg of DNA was sheared using a g‐TUBE (Covaris) and concentrated with AMPure PB magnetic beads (PacBio). A total of 323.85 Gb (~100×) A. flavipes PacBio reads were assembled using canu version 1.7 (Koren et al., 2017) with the error correction module. The corrected subreads were used for initial draft assembly using wtdbg2 version 1.2.8 (Ruan & Li, 2020). Paired‐end (2 × 100 bp) BGI‐SEQ500 data were generated from cerebrum, liver, heart and lung tissue from the same individual by BGI‐Qingdao. To reduce base errors, the assembly was “polished” using pilon version 1.23 (Walker et al., 2014) and 151.43 Gb (50×) 100‐bp paired‐end BGISEQ‐500 reads (mapped to the initial PacBio assembly using minimap2 version 2.10: Li, 2018; and samtools version 1.9: Li et al., 2009).

Genome sizes were estimated by k‐mer frequency analysis (Liu et al., 2013). Briefly, 100‐bp paired‐end whole genome sequence (WGS) reads were used as input into gce (Genomic Charactor Estimator) version 1.0.0 (Marcais & Kingsford, 2011) to obtain the k‐mer frequency and the genome size was estimated using the equation “Genome size = k‐mer number/k‐mer depth,” where the “k‐mer number” is the total number of k‐mers and “k‐mer depth” denotes the peak frequency that occurred more than any other frequencies. Genome length was estimated on the basis of the total scaffold length of the assembly. Using the frequency distribution of 17‐mers of short paired‐end reads (Figure S1), the A. flavipes genome was estimated to be 3.2 Gb.

Assembly quality was assessed using busco (Benchmarking Universal Single‐Copy Orthologs) version 5.0.0_cv1 (Seppey et al., 2019), employing the gene predictor augustus version 3.2.1 (Stanke & Waack, 2003) and the 9226‐gene busco mammalian lineage data set (mammalia_odb10). Although, gene centric, the busco Score is a good predictor of genome completeness (Seppey et al., 2019).

Assignment of scaffolds to chromosomes

A. flavipes Hi‐C libraries from liver tissue were prepared as described previously (Lieberman‐Aiden et al., 2009). Hi‐C libraries were sequenced on the BGISEQ‐500 platform and quality controlled using the hic‐pro (version 2.8.0_devel) pipeline (Servant et al., 2015), resulting in 51.0 Gb uniquely aligned read pairs. Reads validated by hic‐pro were next used to scaffold contigs into seven chromosome clusters (2n‐14 karyotype) using the 3d‐dna version 1.12 pipeline (Dudchenko et al., 2017). The assembly was further improved by interactive correction using juicebox version 1.11.08 (Dudchenko et al., 2018; Durand et al., 2016). A chromatin contact map can be found in Figure S2.

Draft dasyurid genome assemblies

Specimens of the Australian dasyurids A. arktos and A. argentus and the New Guinean dasyurid M. melanurus were sequenced by short‐insert WGS (see Table S1). Reference‐based (i.e., “pseudo”) nuclear genomes were obtained by mapping to the A. flavipes reference genome (see Feigin et al., 2018 and Supporting Methods).

RNA‐seq and transcriptome data processing

Eight A. flavipes specimens in 2017 were sequenced by BGI Hong Kong. The RNA of the remainder (seven males and six females collected in 2018) was extracted by BGI Hong Kong and sequenced by BGI Qingdao. Only those samples collected during the 2018 breeding season were used for differential gene expression analysis (the samples collected in 2017 were used for gene annotation only). Raw data were filtered using flexbar version 3.4.0 (Dodt et al., 2012; Roehr et al., 2017) with default settings (removes reads with any uncalled bases). BGISEQ‐500 adapters were obtained from a comparative study of Illumina and BGI sequencing platforms (Mak et al., 2017). Any residual rRNA reads (the majority removed by poly(A) selection prior to sequencing library generation) were removed from the reads using sortmerna version 2.1b (Kopylova et al., 2012) against the SILVA version 119 ribosomal database (Quast et al., 2013). Tissue transcriptomes were de novo assembled using trinity version 2.8.4 (Grabherr et al., 2011; Haas et al., 2013; Henschel et al., 2012) with default parameters. The A. flavipes transcriptome assemblies were assessed using busco (Seppey et al., 2019). Details of transcriptome sequencing data used in study are listed in Tables S2–S4 (also see Supporting Methods).

Genome annotation

Repeat element annotation

We identified repetitive elements by integrating homology and de novo predictions. Homology‐based transposable elements (TEs) annotations were obtained by interrogating a genome assembly with known repeats in the Repbase database version 16.02 (Bao et al., 2015) using repeatmasker version 4.0.5 (DNA‐level) (Tarailo‐Graovac & Chen, 2009) and repeatproteinmask (protein‐level; implemented in repeatmasker). De novo TE predictions were obtained using repeatmodeler version 1.1.0.4 (Smit & Hubley, 2010) and ltrharvest version 1.5.8 (Ellinghaus et al., 2008) to generate database for a repeatmasker run. tandem repeat finder (version 4.07) (Benson, 1999) was used to find tandem repeats (TRs) in the genome. A nonredundant repeat annotation set was obtained by combining the above data. Genome sequences of other marsupials were also interrogated using the repeat element annotation pipeline. The following were downloaded from DNA Zoo (Virginia opossum [Didelphis virginiana; dv‐2k; Dudchenko et al., 2018]), the NCBI Reference Sequence Database (Release 86) (Pruitt et al., 2007) (gray short‐tail opossum [Monodelphis domestica; MonDom5; GCF_000002295.2], koala [Phascolarctos cinereus; phaCin_unsw_v4.1; GCA_002099425.1], Tasmanian devil [Sarcophilus harrisii; Devil_ref v7.0, also known as sarHar1; GCA_000189315.1], Tammar wallaby [Notamacropus eugenii; Meug_1.1; GCA_000004035.1] and common wombat [Vombatus ursinus; GCF_900497805.2]), and the GigaScience Database (brown antechinus [Antechinus stuartii; Brandies et al., 2020a; Brandies et al., 2020b]).

Protein‐coding gene annotation

Protein‐coding genes of A. flavipes and Didelphis virginiana were annotated using homology‐based prediction, de novo prediction and RNA‐seq‐assisted prediction methods. Sequences of homologous proteins from five mammals (human [Homo sapiens], M. domestica, P. cinereus, S. harrisii and V. ursinus]) were downloaded from NCBI. These protein sequences were aligned to the repeat‐masked genome using blat version 0.36 (Kent, 2002). genewise version 2.4.1 (Birney et al., 2004) was employed to generate gene structures based on the alignments of proteins to the genome assembly. De novo gene prediction was performed using augustus version 3.2.3 (Stanke et al., 2006), genscan version 1.0 (Burge & Karlin, 1997) and glimmerhmm version 3.0.1 (Majoros et al., 2004) with a human training set. Transcriptome data were mapped to the assembled genome using hisat2 version 2.1.0 (Kim et al., 2019) and samtools version 1.9 (Li et al., 2009), and coding regions were predicted using transdecoder version 5.5.0 (Grabherr et al., 2011; Haas et al., 2013). A final nonredundant reference gene set was generated by merging the three annotated gene sets using evidencemodeler version 1.1.1 (EVM) (Haas et al., 2008) and excluding EVM gene models with only ab initio support. The gene models were translated into amino acid sequences and used in local blastp (Camacho et al., 2009) searches against the public databases Kyoto Encyclopedia of Genes and Genomes (KEGG; version 89.1) (Kanehisa & Goto, 2000), Clusters of Orthologous Groups (COG) (Tatusov et al., 2000), NCBI nonredundant protein sequences (NR; version 20170924) (O'Leary et al., 2016), Swiss‐Prot (release‐2018_07) (UniProt Consortium, 2012), TrEMBL (TRanslation of EMBL [nucleotide sequences that are not in Swiss‐Prot]; release‐2018_07) (O'Donovan et al., 2002), and InterPro (version 69.0) (A. L. Mitchell et al., 2019). A total of 18,068 (93.5%) of A. flavipes genes could be functionally annotated. Where specific genes are named in this paper, human nomenclature assignments are used unless otherwise noted.

Phylogeny and divergence time estimation

Marsupials (metatherians) and eutherians diverged ~160 million years ago (Ma), long before the radiation of extant eutherian clades (~100 Ma) (Luo et al., 2011; Phillips et al., 2009). We employed the platypus (Ornithorhynchus anatinus) as an outgroup as this species is frequently used as an outgroup for both eutherian and marsupial mammals and it has a high‐continuity, well‐annotated genome (contig N50 15.1 Mb; scaffold N50 83.3 Mb). In contrast, the genomes of the most closely related eutherian mammals, the edentates (sloths, armadillos and anteaters) currently have relatively poor‐quality genome assemblies. These include, for example, the armadillo Dasypus novemcinctus (contig N50 0.03 Mb, scaffold N50 1.7 Mb) and the sloth Choloepus hoffmanni (contig N50 0.06 Mb; scaffold N50 0.4 Mb).

We identified 6090 high‐confidence 1:1 orthologues by interrogating the predicted proteins from the gene models of 11 species (ten marsupials and the platypus) using sonicparanoid version 1.3.0 (Cosentino & Iwasaki, 2019). Because the A. arktos, A. argentus and M. murexia pseudogenomes were derived from A. flavipes, their genes were extracted using the A. flavipes gene annotation file. The corresponding coding sequences (CDS) for each species were aligned using prank version 100802 (Loytynoja & Goldman, 2005) and filtered by gblocks version 0.91b (Talavera & Castresana, 2007) to identify conserved blocks (removing gaps, ambiguous sites and excluding alignments less than 300 bp in size), leaving 6090 genes. Maximum‐likelihood (ML) phylogenetic trees were generated using raxml version 7.2.8 (Stamatakis, 2006) and iq‐tree version 2.1.3 (Minh et al., 2020) with three CDS data sets: the whole coding sequence (whole‐CDS), first codon positions and fourfold degenerate (4d) sites. Identical topologies and similar support values were obtained (1000 bootstrap iterations were performed). The divergence time between species was estimated using mcmctree (a Bayesian molecular clock model implemented in paml version 4.7; Yang, 2007) with the GTR+G nucleotide substitution model, and the whole‐CDS ML tree and concatenated whole‐CDS supergenes as inputs. We used 100,000 iterations after a burn‐in of 10,000 iterations. mcmctree calibration points (million years ago; Ma) were obtained from timetree (Kumar et al., 2017): O. anatinus–P. cinereus (~167–192 Ma), V. ursinus–D. virginiana (~72–86 Ma), V. ursinus–M. melanurus (~56–64 Ma), V. ursinus–P. cinereus (~31–39 Ma), M. domestica–D. virginiana (~24–39 Ma) and S. harrisii–M. melanurus (~4–22 Ma). For comparison, phylogenetic trees of marsupials in the marsupial orders Dasyuromorphia and Didelphimorphia were obtained by querying the PHYLACINE (The Phylogenetic Atlas of Mammal Macroecology) resource (Faurby et al., 2018).

Genomic diversity

Genome‐wide heterozygosity (π) was calculated from (2 × 100 bp) WGS read mapped to the S. harrisii genome. For details, see Supporting Methods.

Demographic history

We inferred demographic histories by applying the Pairwise sequentially Markovian coalescent (PSMC) method (Li & Durbin, 2011). Short‐insert paired‐end WGS reads of A. flavipes (two individuals), A. arktos, A. argentus, A. stuartii (two individuals; see Brandies et al., 2020a, 2020b) and M. melanurus were mapped to a repeat‐masked S. harrisii genome. Putative X‐linked scaffolds in the S. harrisii assembly were excluded (see Supporting Methods). Geographical maps were generated using the R package “ggmap” (Kahle & Wickham, 2013).

Transcriptome profiling of Antechinus flavipes breeding

To search for transcripts associated with obligate semelparity, we performed differential gene expression analysis of a single A. flavipes population (Samford Valley, Queensland, Australia) from August to October 2018. Every 2 days, we microscopically examined female urine for cornified, anuclear vaginal epithelial cells, a sign of peak oestrus (Selwood, 1980). Concurrent visual inspections of captured males for behavioural changes (e.g., agitation) revealed evidence of stress beginning in the last week of September, followed by progressive loss of body condition (Figure S3). The practical and ethical challenges of longitudinal sampling from a wild population limited the number of samples that could be obtained. In this study, we present differential expression data of eight male and six female antechinuses (liver, kidney, skeletal muscle, cerebrum and reproductive tissues [testis and ovary]). See Supporting Methods and Figures S4 and S5 for data processing, assessment and differential expression procedures, and Appendix S1 for normalized counts.

RESULTS AND DISCUSSION

Genome assembly and annotation

We used PacBio sequencing and Hi‐C (Burton et al., 2013) to create a 3.2‐Gb (2n = 14) chromosome‐level male Antechinus flavipes nuclear genome assembly with a contig N50 of 51.8 Mb and scaffold N50 of 636.7 Mb (Table 1 and Figure 1). The seven chromosomes ranged in size from 86.6 to 727.3 Mb (Table 2) and covered ~99.7% of the assembly. We also generated A. flavipes reference‐based assemblies of Antechinus arktos (female), Antechinus argentus (female) and Murexia melanurus (male) from short‐insert WGS reads. Most A. flavipes genome regions had a GC‐depth of about 50×, while a small peak at 25× depth revealed heterozygous sex chromosomes (Figure S6). Hi‐C was able to recover an ~86.5‐Mb X chromosome (chromosome 7), while the Y chromosome could not be obtained (see Supporting Methods section S2.0). This is unsurprising, as marsupial Y chromosomes are small (~10–12 Mb; akin to microchromosomes of birds), repeat‐rich and refractory to assembly (Deakin & O'Neill, 2020). The chromosomes of A. flavipes were highly homologous to the related Tasmanian devil (Sarcophilus harrisii) (Figure 1c), supporting the quality of the assembly.

TABLE 1

Antechinus flavipes genome assembly statistics

Statistic	Assembly method
	PacBio (100×) + BGISEQ‐500 (50×) polishing	PacBio + BGISEQ‐500 polishing + Hi‐C
Contig number	1107	1107
Contig length (bp)	3,192,388,871	3,192,388,871
Contig N50 (bp)	51,843,704	51,843,704
Contig max length (bp)	168,884,005	168,884,005
Scaffold number	1105	487
Scaffold N50 (bp)	51,843,704	636,717,922
GC content (%)	36.21	36.21

TABLE 2

Summary of assembled Antechinus flavipes pseudochromosomes

Chr ID	Sequence number	Sequence length (bp)
Chr 1	88	727,308,865
Chr 2	83	686,337,178
Chr 3	146	636,717,922
Chr 4	50	486,195,733
Chr 5	104	300,479,852
Chr 6	63	264,590,327
Chr 7	91	86,570,782

The GC content of A. flavipes (~36%) (Table 1) is similar to the brown antechinus (A. stuartii) and S. harrisii. As indicated by 17‐mer frequency analysis (Figure S1), the A. flavipes genome is repeat‐rich. Repetitive elements account for 51.8% (~1.7 Gb) of the assembly, with long interspersed elements (LINEs; ~45%), long terminal retrotransposons (LTRs; 15.4%) and short interspersed elements (SINEs; 6.6%) being the major classes of TEs (Tables S5 and S6). We annotated 24,708 A. flavipes protein‐coding genes (82.2% supported by transcriptome data from 12 tissues). The A. flavipes reference assembly obtained a busco (Benchmarking Universal Single‐Copy Orthologs) (Seppey et al., 2019) genome completeness score of 92.4%—comparable to A. stuartii (92.2% and 92.3% for a female and male assembly, respectively) (Brandies et al., 2020a), the koala (Phascolarctos cinereus; 93.9%) (Johnson et al., 2018) and S. harrisii (91.6%; 2019 assembly mSarHar1.1) (Table S7). The scores of the three dasyurid reference‐based assemblies, here used for phylogenetic and demography analyses, ranged from 74.2% to 90.0% recovered complete BUSCO genes (Table S8).

Phylogeny and demographic history

To construct a time‐calibrated species tree (Figure 2a), we identified 6090 high‐confidence single‐copy orthologues from 10 marsupial species and a monotreme, the platypus (Ornithorhynchus anatinus). Extant Australasian and American marsupials diverged from a South American ancestor ~80 Ma (Kumar et al., 2017; Nilsson et al., 2010), emphasizing the deep evolutionary history of marsupials. We estimated that the lineages leading to extant Dasyuridae (e.g., Antechinus) and Didelphidae diverged 71.7 Ma (70.1–72.7 Ma 95% credibility interval [CI]). The dasyurid subfamily Phascogalinae comprises three genera: the semelparous Australian Antechinus and Phascogale, which have male die‐off, and the iteroparous New Guinean sister genus Murexia. Small dasyurids from mainland Australia are likely to have dispersed into New Guinea 9–11 Ma (Mitchell et al., 2014). We found the divergence time of Antechinus and Murexia to be about 11.0 Ma (95% CI 8.3–13.3). Our analysis shows that the four antechinus species split from a common ancestor ~6.8 Ma (95% CI 6.1–7.3), while A. flavipes and A. argentus + A. stuartii share an ancestor ~5.9 Ma (95% CI 5.3–6.5). Our phylogeny and divergence times concord with trees generated from the PHYLACINE database (Faurby et al., 2018) (Figure S7) and with recent estimates from a small set of mitochondrial and nuclear genes (Mutton et al., 2019).

FIGURE 2

Evolutionary history of semelparous mammals. (a) Inferred phylogeny of ten marsupials and the platypus (outgroup) based on whole‐coding sequences of 6090 1:1 orthologues. Blue numbers at nodes represent the estimated divergence time from present (million years ago; Ma) between lineages. Semelparous lineages are indicated in red. (b) Demography of the dasyurid subfamily Phascogalinae. Historical population sizes (N e) were estimated using the pairwise sequentially Markovian coalescent (PSMC) method (Li & Durbin, 2011) and diploid genome sequences. The x‐axis shows the time before present in years on a log scale; the dashed grey line shows human arrival in Australia 59.3–65 ka (Clarkson et al., 2017); red dots indicate the extinction times of 15 species of Pleistocene Sahul (Saltre et al., 2016). The y‐axis shows N e, with bootstrap estimates indicated by lighter lines. Plots were scaled using a mutation rate (μ) of 1.37 × 10−9 substitutions per nucleotide per generation and species‐specific generation times (g) of 1 year. Two Antechinus flavipes (coloured in orange and blue) and two Antechinus stuartii (in brown and black) individuals were interrogated. (c) Map showing the location (coloured circles) and elevation (m amsl; metres above sea level) of sampled dasyurids (see panel b)

The PSMC method (Li & Durbin, 2011), which can infer changes in the effective population size (N e) over the Pleistocene, was used to evaluate four sequenced antechinus species and M. melanurus (Figure 2b,c). Until the end of the Pleistocene, ~11,700 years ago, mainland Australia, New Guinea and Tasmania formed a single continent, Sahul (Barrows et al., 2002). Glaciation of Sahul was restricted to the Snowy Mountains of southeastern mainland Australia and the Tasmanian highlands, and the continent experienced cycles of aridity and sea‐level fluctuations (Barrows et al., 2002). Evidently, during the Early to Mid‐Pleistocene until about 250,000 years ago, all species had larger population sizes, which is probably a reflection of larger expanses of connected forest habitats (Baker & Van Dyck, 2013). After 250,000 years ago, however, the species occurring in higher‐elevation, wetter habitat, A. arktos, A. argentus and M. melanurus, apparently declined in population size, whereas A. flavipes and A. stuartii populations, which are tolerant of drier conditions, underwent expansion (Figure 2a,b). This pattern persisted into the Late Pleistocene ~11,700–126,000 years ago, a period when temperature fluctuations were large, and many mesic‐adapted species declined, were restricted or became extinct; more arid‐adapted species radiated into burgeoning desert landscapes (Saltre et al., 2016). Our results are generally in accordance with a recent study suggesting that aridity and geographical isolation promoted dasyurid speciation (García‐Navas et al., 2020).

While the PSMC method cannot infer population dynamics after ~20,000 years ago, it is noteworthy that A. flavipes had the largest N e in the most recent PSMC estimate. This species is one of the most common, arid‐adapted and broadly distributed species in the genus Antechinus today (Baker & Van Dyck, 2013; Lada et al., 2008; Mutton et al., 2019). There is an apparent reduction in population size of A. stuartii ~20,000–60,000 years ago. Interestingly, although A. stuartii is generally considered a common antechinus species, evidence is emerging that it should be considered a species complex, representing multiple lineages, and thus possesses a smaller distributional range than presently recognized (Viacava et al., 2021). Moreover, we note that the A. stuartii individuals included here were sampled from Lane Cove National Park, a site with a low current population size (NSW Department of Environment & Conservation, 2004). In conclusion, our analysis suggests that climate change during the Early (~2.5–0.78 Ma) (Lisiecki & Raymo, 2005) and Mid‐Pleistocene (0.78–0.126 Ma) (Elderfield et al., 2012) has shaped the evolution of the dasyurid subfamily Phascogalinae by isolating populations in different habitats, resulting in some species being now greatly restricted in their distribution.

We found that all small dasyurid taxa have much higher (100 orders of magnitude) genome‐wide heterozygosity (π from 0.048 to 0.096) than the endangered S. harrisii (π = 0.00032 (Miller et al., 2011)) (Figure S8). Moreover, the heterozygosity of A. flavipes is twofold higher than that of A. arktos, A. argentus and M. melanurus. This may reflect its larger historical and current population size (Figure 2b,c). Although we are comparing individuals from single populations, combined with factors such as environmental change, loss of habitat and population fragmentation, the lower genetic diversity of A. arktos and A. argentus may reflect their present‐day low population and threatened species status (Mutton et al., 2019; Riordan et al., 2020). Future studies including assessment of geographical population structures and the inclusion of historical and current specimens (Diez‐Del‐Molino et al., 2018) are required to assess the threat level of extant small dasyurids. Our high‐quality A. flavipes assembly serves as an ideal reference genome for population genomic studies of A. flavipes and the 14 other antechinus species.

A gene expression profile of male die‐off in Antechinus

We profiled gene expression during antechinus reproduction by surveying multiple tissues from an A. flavipes population in Samford Valley in Queensland, Australia, during the 2018 breeding season (Movie S1). We observed reproductive behaviour and peak oestrus on September 15, 2018 (see Supporting Methods for details). Adult male A. flavipes disappear from this site after the first week of October. Antechinus male die‐off is characterized by a dramatic increase in free cortisol (Figure 3a), mediated by a surge in testosterone which decreases the production of corticosteroid‐binding globulin (Boonstra, 2005; Naylor et al., 2008). The potent effects of cortisol mean that the associated gene expression and its consequences (Le Phuc et al., 2005; McDowell et al., 2018) become a suitable backdrop against which comparisons can be made. While corticosteroids regulate a large number of genes via the glucocorticoid receptor (GR) (Le Phuc et al., 2005; McDowell et al., 2018), complicating gene expression analyses, chronic (i.e., breeding male antechinus) exposure to free corticosteroids has distinct outcomes (Manoli et al., 2007). The steroid hormone mediates a stress response, alters metabolism and mobilizes glucose as energy from a series of stores, initially enhancing short‐term survival.

FIGURE 3

Gene expression profiling of the Antechinus flavipes breeding cycle. (a) Schematic overview of the yellow‐footed antechinus (A. flavipes) breeding season in Samford Valley, Queensland, Australia. Free‐CORT illustrates an increase in the level of free (unbound) stress hormones (chiefly the corticosteroid cortisol) in males during the breeding season of antechinuses (based on three studies illustrated in Naylor et al., 2008). (b) Rapid organismal senescence of male A. flavipes. Left, male (~11 months of age) prior to the breeding season. Right, late breeding‐stage male (~11.5 months of age). Note the loss of fur and deterioration of body condition. (c) In the moribund state, when other energy stores have been depleted, skeletal muscle catabolism (induced by the corticosteroid cortisol) supports whole‐body metabolic demands. (d) Venn diagram for the overlap of significantly up‐regulated (left) or down‐regulated (right) genes in moribund male antechinuses. Called by limma version 3.22.1 (≥1.5‐fold change and p ≤ .01). (e) Heat map of selected energy metabolism genes differentially expressed in the skeletal muscle of moribund male antechinuses. Moribund samples are indicated by a red bar below the map. Scaled TMM‐normalized read counts (denoted as the row Z score) are plotted, with red indicating high expression and blue indicating low expression. BPGM, 2,3‐bisphosphoglycerate mutase; FBP2, fructose‐bisphosphatase 2; GPD1, glycerol‐3‐phosphate dehydrogenase 1; GPD2, glycerol‐3‐phosphate dehydrogenase 2; GPI, glucose‐6‐phosphate isomerase; LDHA, lactate dehydrogenase A; PFKM, phosphofructokinase, muscle; PGK2, phosphoglycerate kinase 2; PKM, pyruvate kinase M1/2; PRKAG3, protein kinase AMP‐activated noncatalytic subunit gamma 3; SLC2A5, solute carrier family 2 member 5; SLC37A4, solute carrier family 37 member 4; STK11, serine/threonine kinase 11. (f) Expression pattern of the AMPK regulator LKB1 (STK11) in male (left) and female (right) tissues throughout the breeding season. Counts per million TMM‐normalized counts, with standard error bars

Because of a limited number of samples in our current data set, we considered genes that were differentially expressed at the end of the breeding season (early October), when all males die (are moribund; right panel in Figure 3b). Each sex was compared separately. There were more differentially expressed genes (limma‐modified Student's t‐test p ≤ .01) in the male liver and skeletal muscle than in the kidney and cerebrum (Figure 3c and Appendix S2), and most of these were downregulated. Genes upregulated in moribund male livers showed enrichment of biological process gene sets attributed to the innate immune response (bootstrap p ≤ .0001) (Table S9). However, several cytokine genes (PDGFRA, PDGFRB, PRNP, CD248, CD79A, CXCR5 and CXCR6) were downregulated in the moribund liver (Appendix S2). Characteristic liver pathology of the moribund male antechinus includes infections by bacteria such as Listeria monocytogenes (Barker et al., 1978; Bradley et al., 1980), viruses (e.g., dasyurid herpesvirus 1; Amery‐Gale et al., 2014) and parasitic infestation (Poskitt et al., 1984). Cortisol‐mediated immunosuppression in the preceding two weeks is likely to lead to poor late‐stage resistance to invading pathogens. We identified eight genes that were downregulated and two genes that were upregulated by both liver and skeletal muscle of moribund males (Figure 3d). This included decreased expression of the glycophagy receptor STBD1 (Johansen & Lamark, 2020) and increased expression of ASNS, a gene induced by amino acid starvation stress (Balasubramanian et al., 2013). Downregulation of SLC2A4 (encodes the glucose transporter GLUT4; Amoasii et al., 2019) and NUR4A2 (NUR11; a regulator of GLUT4 expression; Amoasii et al., 2019) indicates reduced glucose uptake by the skeletal muscle of moribund males. KEGG pathway and biological process enrichment analysis revealed metabolic dysregulation of skeletal muscle in moribund males, with marked downregulation of genes associated with glycolysis, gluconeogenesis and insulin signalling—indicating insulin resistance (bootstrap p ≤ .0001) (Figure 3e; Tables S9 and S10). These genes included PRKAG3 (−4.7‐fold; limma‐modified Student's t‐test, p = 3.7 × 10−4), which encodes a subunit of a key energy‐sensing molecule, AMP‐activated protein kinase, and its regulator LKB1 (STK11; −2.0‐fold, p = 2.4 × 10−3) (Appendix S2). AMPK responds to energetic stress to regulate multiple metabolic pathways. Reduced expression of LKB1 was limited to the skeletal muscle of moribund A. flavipes males (Figure 3f). LKB1 has a tissue‐specific role in energy metabolism. Complete or partial loss of Stk11 expression in mouse skeletal muscle prevents AMPK activation and causes metabolic defects such as acceleration of ageing‐induced myopathy (Bujak et al., 2015; Koh et al., 2006; Sakamoto et al., 2005; Thomson et al., 2007).

Our work is the first gene expression profile of a semelparous mammal. In agreement with observational data from A. flavipes and other antechinuses (Bradley et al., 1980; Naylor et al., 2008; Woollard, 1971), our results indicate near‐complete depletion of fat stores and skeletal muscle atrophy—leading to a catastrophic energy crisis and an inability to maintain cognitive and physical function in moribund males. Further experimental work on the endocrine system of semelparous mammal species is needed to characterize the role of testosterone and the dysregulation of the hypothalamic–pituitary–adrenal axis and negative feedback by cortisol. We were not able to capture gene expression changes in brain regions other than the whole cerebrum of A. flavipes during the breeding season, limiting our insights into the effects on major organs. Nevertheless, the A. flavipes reference genome permits future studies on gene expression and regulation (e.g., glucocorticoid receptor ChIP‐seq [chromatin immunoprecipitation sequencing]).

CONCLUSION

Our study and accompanying data sets, including the first chromosome‐level genome of a semelparous mammal, provide a critical resource for future studies on semelparity and its association with the universal ageing process. Complementary genomic data, including the genomes of additional Australasian and South American semelparous and iteroparous marsupials, are forthcoming and promise to shed light on whether semelparity in lineages with a common ancestor ~80 Ma evolved by independent, lineage‐specific molecular changes or by shared molecular convergence. The genomic resources also provide tools for broader studies on the ecology, evolution and conservation of marsupials.

AUTHOR CONTRIBUTIONS

I.S. and G.F. initiated and coordinated the project. C.P., K.M. and A.M.B. collected Antechinus flavipes samples. L.K.C. dissected A. flavipes. A.M.B. provided ear clippings from Antechinus arktos and Antechinus argentus. P.B.T. and L.K.C. provided frozen tissue samples for whole‐genome, Hi‐C and transcriptome sequencing. C.S. led the sequencing and assembly efforts with G.F. R.T., K.H. and C.S. contributed to genome assembly and annotation. R.T., K.H., C.S. and I.S. performed comparative genomics analyses. The de novo and reference‐based assemblies of A. arktos, A. argentus and M. melanurus were generated by I.S. I.S. performed population history analysis. K.H. performed diversity analysis. I.S. assembled transcriptomes and performed gene expression analysis. R.T. performed molecular evolution analyses. L.K.C., D.O.F., N.O.L. and A.M.B. provided direction on marsupial ecology and physiology. R.T., L.K.C., D.O.F., A.M.B. and I.S. wrote the paper with input from all authors.

OPEN RESEARCH BADGES

This article has earned an Open Data Badge for making publicly available the digitally‐shareable data necessary to reproduce the reported results. The data is available at https://github.com/sciseim/semelparity-genomics.

Supplementary Material

Click here for additional data file.

Movie S1

Click here for additional data file.

Appendix S1

Click here for additional data file.

Appendix S2

Click here for additional data file.