Comparative Genomic Characterization of the Multimammate Mouse Mastomys coucha.

Mastomys are the most widespread African rodent and carriers of various diseases such as the plague or Lassa virus. In addition, mastomys have rapidly gained a large number of mammary glands. Here, we generated a genome, variome, and transcriptomes for Mastomys coucha. As mastomys diverged at similar times from mouse and rat, we demonstrate their utility as a comparative genomic tool for these commonly used animal models. Furthermore, we identified over 500 mastomys accelerated regions, often residing near important mammary developmental genes or within their exons leading to protein sequence changes. Functional characterization of a noncoding mastomys accelerated region, located in the HoxD locus, showed enhancer activity in mouse developing mammary glands. Combined, our results provide genomic resources for mastomys and highlight their potential both as a comparative genomic tool and for the identification of mammary gland number determining factors.

Mastomys are the most widespread African rodent (Colangelo et al. 2013) (fig. 1) and are members of the Praomyini tribe, the closest clade to Murini, which includes Mus musculus (Lecompte et al. 2008). Mastomys are of immense research interest as they are carriers of various diseases such as the plague and Lassa virus (Bonwitt et al. 2017), are a good animal model for gastric cancer (Nilsson et al. 1992), Papillomavirus (Helfrich et al. 2004) and Schistosoma (Lurie and De Meillon 1956), and are a model for chromosomal diversity and speciation. They are also notable in that females carry a well-developed prostate gland (Snell and Stewart 1965).

. 1

Mastomys coucha phylogeny, range, and morphology. (a) Rodent phylogeny showing the number of mammary glands for each species. (b) Mastomys range and sampling locations in Africa obtained from the International Union of Conservation red list of threatened species. The range of M. coucha is marked in yellow and Mastomys natalensis in red. (c) A picture of a M. coucha from our colony. (d) The ventral side of a M. coucha showing 16 mammary glands.

An intriguing trait found in the genus Mastomys is the presence of an unusually high number of mammary glands, anywhere between 16 and 24, which is a rapid gain of 6–14 glands in the 10.2 ± 0.6 My since their divergence with Mu. musculus and 11.3 ± 0.5 from Rattus (Lecompte et al. 2008) (fig. 1). This is not correlated with their body size, which is on average 90 mm, compared with mouse and rat at 79 and 220 mm, respectively (Hayssen et al. 1993). Mammary glands form during early embryogenesis along the parallel lines on the ventral ectodermal surface. This line, also called the mammary line, is a region where Wnt family member 6 (Wnt6), Wnt10A, and Wnt10B (Chu 2004; Eblaghie et al. 2004; Veltmaat et al. 2004) and other signaling molecules interact to define the ectodermal placodes. Previous work on the five pairs of glands in the lab mouse found that knockouts of various developmental genes (Eda, Edar, Hoxc6, Hoxc8, Fgf10, Gli3, Tbx3, etc.) can lead to changes in the number of mammary glands (Veltmaat 2017). Enhancers that regulate these genes have also been associated with mammary gland development. For example, mammary expression of HoxD9 was found to be regulated by eutherian-specific changes in a nearby enhancer (Schep et al. 2016). Combined, this work suggests that each pair of mammary glands is the result of a distinct combination of signaling pathways. Despite having a partial gene and pathway list that is important in mammary development, the molecular components that determine mammary gland number remain largely unknown.

We established a colony of Mastomys coucha in our lab (fig. 1). This colony was derived from a colony that was founded in the National Institute of Health and has been inbred since 1983 (Modlin et al. 1988). Analysis of individuals from the colony found that they have eight pairs of mammary glands: three pectoral, three abdominal, and two inguinal with occasional supernumerary mammary glands in the anterior region (fig. 1). The average litter size we obtained from this colony was 3.2.

To generate a genome for M. coucha, we extracted DNA from the liver of an adult male and generated both short and mate-pair (2-, 5–6-, and 8–10-kb) Illumina sequencing libraries. The M. coucha genome, which we named mc1, was mapped to 26 pseudochromosomes. Mc1 is 2.48 Gb in size with 78× coverage and a contig and scaffold N50 of 29 kb and 3.6 Mb, respectively. The heterozygosity rate is 0.86 per kb, consistent with a highly inbred genome. The mitochondrial genome was assembled separately providing a 16.2-kb genome containing a complete set of mitochondrial genes. The genome was also annotated with transcriptomes from five tissues: adult testis, skin, spleen, embryonic day (E) 17.5 brain, and E11.5 developing mammary glands. These transcriptomes were de novo assembled and the mc1 assembly was assessed for completeness with BUSCO (Simão et al. 2015), finding it to contain 85% of the core vertebrate genes. Using RepeatMasker (Smit et al. 2015), we found that repeats make up 36.4% of the M. coucha genome.

As our colony was inbred, we also sequenced at low coverage (10×) a M. coucha collected from the wild to determine nucleotide variation. We found that the wild samples have a heterozygosity rate of 5.78 per kb, which is over five times higher than the lab strain, adding further support that the mastomys we selected for genome sequencing is inbred. We also sequenced at low-coverage (11×) Mastomys natalensis to determine genus level conservation. We found an even higher rate of heterozygosity, 13.41 per kb, compared with our M. coucha collected in the wild, suggesting that this species has a higher population size, as observed in a previous study (Sands et al. 2015).

As Mastomys diverged 10.2 ± 0.6 My from mouse and 11.3 ± 0.5 from rats (Lecompte et al. 2008) (fig. 1), it occupies a strategic position in the rodent phylogenetic tree and hence can serve as a robust comparative genomic tool for these commonly used animal models (fig. 2). Whole-genome alignment found that mc1 contains 90.1% and 85.5% of the mouse and rat genomes, respectively. The average length of synteny between mouse and mastomys is 3.8 Mb, comparable to the 2.5 Mb between mouse and rat (Bourque et al. 2004). We calculated the percent of mouse variation in the whole-genome alignment with all three species (mouse-mastomys-rats) and found that 84% of the mouse-rat divergence could be assigned to a branch of either the rat or the mouse lineage, providing a resource for rodent-specific evolution. Combined, our results show that mastomys provide a useful comparative genomics tool as a sister group to mouse and rat.

. 2

Comparative genomics. (a) Synteny comparisons between mastomys-mouse-rat. A circos plot showing chromosomal segments (in orange) that are conserved between mastomys (blue), mouse (red), and rat (green). (b) Pie chart showing the different genomic locations of MARs. (c) GREAT analysis for MARs showing their GO biological process. (d) Diagram showing the location of the coding MAR178 within the Klf5 protein. Green circles mark nonsynonymous protein changes and black circles mark amino acid insertions in mastomys. The orange line represents the ubiquitin binding domain (WWP1) and the green rectangles represent the zinc finger binding domains. The blue line shows the location of MAR178, which overlaps exon 2 of the Klf5 gene.

Accelerated regions represent a good starting point to search for genetic differences leading to lineage-specific traits and have been frequently found to function as enhancers in other species (Capra et al. 2013; Booker et al. 2016). We thus set out to identify sequences that are specifically accelerated in mastomys, as they could pose as potential drivers of mammary gland number increase in mastomys. A seven-way mammalian multiple genome alignment was created using the mastomys, mouse, rat, beaver, human, goat, and cat genomes. We used phastCons (Siepel 2005) to find highly conserved regions in the mammalian alignments (without including mc1). We then used phlyoP (Pollard et al. 2010) to identify accelerated evolution in Mastomys relative to a neutral model of evolution. Using a FDR cutoff of 0.05, we found 515 Mastomys accelerated regions (MARs), the majority (85%) of which are noncoding (fig. 2). Using the Genomic Regions Enrichment of Annotations Tool (GREAT [McLean et al. 2010]), we found that MARs are adjacent to genes enriched for several gene ontology (GO) biological processes including “chromatin modification” and “regulation of epithelial cell proliferation” (fig. 2).

We found that over 15% (N = 80) of our MARs were located in exons, with three of them located in exons of genes associated with mammary gland development (Baz2a, Flt4, and Klf5; supplementary table S1, Supplementary Material online). For instance, the MAR found in the Kruppel like factor 5 (Klf5), a gene that is known to promote epithelial growth in breast tissues (Chen et al. 2002), leads to 20 amino acid nonsynonymous changes compared with its mouse homolog (fig. 2). Although these changes are not within the zinc finger binding domains, they are located in regions that determine protein–protein interactions.

As the majority of MARs are within noncoding sequences (85%), we set out to investigate the possibility of MARs functioning as gene regulatory elements for mammary gland development in Mastomys. To this end, we looked at possible co-occurrence of MARs and expressed genes within topologically associated domains (TADs) previously identified in mouse (Dixon et al. 2012). TADs are conserved genomic regions where DNA sequences interact more frequently with each other than outside and represent functional genomic units. We found 242 MARs residing in TADs containing mammary gland expressed genes which were derived from our RNA-seq data set from the E11.5 developing mammary glands (supplementary table S1, Supplementary Material online). Interestingly, two reside within a TAD with genes that have a known mammary gland phenotype when mutated in mice: sclerostin domain containing 1 (Sostdc1) and integrin subunit alpha 1 (Itga1) (Veltmaat 2017) (supplementary fig. S1, Supplementary Material online). Sostdc1 is of particular interest as it is known to modulate Bmp and Wnt signaling pathways and control the shape and size of mammary buds (Närhi et al. 2012).

To specifically test whether MARs could function as enhancers within the developing mammary gland, we selected MAR456, which is located nearby the HoxD cluster (fig. 3), for mouse transgenic enhancer assays. Interestingly, this MAR was also found to be a therian-specific accelerated region, sequences that were found to be accelerated in eutherians and are associated with the evolution of mammalian-specific traits (Holloway et al. 2016). We cloned both the M. coucha and mouse sequence into a mouse enhancer assay vector and tested it for E11.5 and E13.5 enhancer activity in transgenic mice. At E11.5, we observed an overall stronger enhancer expression for the mastomys sequence compared with mouse (fig. 3 and supplementary fig. S2, Supplementary Material online). However, we should note that this assay is not quantitative, as each embryo can have a different transgene copy number. We also noted that several of the embryos having the mastomys sequence showed enhancer activity in the limb and brachial arch, which was not observed in the embryos injected with the mouse sequence. At E13.5, a stage where the mammary glands are visible, we observed staining for both sequences in the developing mammary glands and did not see any inherent expression differences between mastomys and mouse (fig. 3).

. 3

MAR456 enhancer expression. (a) UCSC Genome Browser snapshot showing the mouse HoxD genomic locus along with the mastomys and therian accelerated regions tracks. (b) Representative transgenic embryo showing MAR456 enhancer activity for either the mouse or the mastomys sequence at E11.5 and E13.5. Zoom in of the mammary gland expression (red arrow) for the E13.5 embryos are shown in the right panel.

In summary, using various genomic tools, we characterized the multimammate M. coucha. Mastomys are close relatives of mouse and rat and as such pose as a great comparative genomic tool for these commonly used animal models. Our mc1 genome showed extensive alignments and synteny blocks (average length 3.8 Mb) with these two genomes. In addition, it allowed assignment for 84% of the mouse-rat divergence to a branch of either the rat or mouse lineage, highlighting the comparative genomic advantage of this genome as a sister group to mouse and rat. Using low-coverage sequencing of M. coucha and M. natalensis collected from the wild, we provide additional information on nucleotide variation and diversity in these species. The generation of additional mastomys genomes could enhance the use of this species as a strong comparative genomics tool for the rodent clade in general and for the identification of sequences controlling mammary gland number.

Materials and Methods

Genome Assembly

DNA was extracted from the liver of a male M. coucha using the Qiagen AllPrep DNA/RNA Mini Kit. Modified versions of the 4-μg protocol of the Nextera Mate-Pair Sample Preparation kit (Illumina) were used to generate libraries with insert sizes of 2, 5–6, and 8–10 kb. Multiple reactions were pooled before size selection, with five reactions of 4.5 μg grouped for the 5–6- and 8–10-kb size range and five reactions of 2 μg to enrich the 2-kb size range. The smaller-insert libraries, 220  and 350 bp, were generated with the TruSeq DNA PCR-Free Library Preparation kit (Illumina), following the manufacturer’s instructions. DNA was sonicated with the Diagenode Bioruptor Plus, end repaired, A-tailed, and ligated with barcoded adaptors. All libraries were sequenced on the Illumina HiSeq 2500 platform with 2× 125-bp high output mode. The 220- and 350-bp paired reads were trimmed on either side using Trimgalore (Krueger 2015). The mate-pair read pairs were trimmed with Skewer(Jiang et al. 2014). All libraries were error corrected with lighter. Trimmed and corrected short-fragment libraries were used to determine the optimal kmer size of 59 using KmerGenie (Chikhi and Medvedev 2014). We used Meraculous (Chapman et al. 2011) to assemble the genome using the short-fragment libraries for contig creation and the mate-pairs for scaffolding. Ragout (Kolmogorov et al. 2016) was used for pseudochromosome mapping. The mitochondrial genome was assembled separately using the short insert library with NOVOPlasty (Dierckxsens et al. 2016). The assembly was submitted as whole-genome sequencing data under BioProject PRJNA406979. Heterozygosity was estimated with BWA (Li and Durbin 2010) and samtools (Li et al. 2009).

Genome Annotation

The M. coucha genome was annotated using the MAKER2 pipeline (Holt and Yandell 2011). Repeat masking was done using RepeatMasker (Smit et al. 2015) and the RepBase annotation set (Bao et al. 2015) to softmask 36% of the genome. Preliminary ab initio gene models were created using SNAP (Korf 2004), AUGUSTUS (Stanke and Morgenstern 2005), and GeneMark-ES (Besemer and Borodovsky 2005). EST data from transcriptomes were included from mapping embryonic brain and mammary gland, as well as adult testis, skin, and spleen. Reads from these libraries were mapped to the draft genome using HISAT2 (Kim et al. 2015) and gene models were created with StringTie (Pertea et al. 2015). The transcriptomes were also de novo assembled with Trinity (Grabherr et al. 2011) and all models with TMP > 1 were retained. De novo assembled contigs from M. natalensis blood previously published (Andersen et al. 2015) were also used. Mouse evidence from Ensemble 89 cDNA was mapped with BlastN and proteins were mapped using TBlastN (Camacho et al. 2009). Noncoding annotation of tRNAs was performed with tRNAscan (Lowe and Eddy 1997). All de novo and functional evidence was merged in MAKER2 to create gene models. Gene models with annotation edit distance <0.9 were annotated with blast matches to UniProt/Sprot (UniProt Consortium 2017). Ensemble transcript IDs were mapped with reciprocal best blast hits to the MAKER2 transcriptomes. Assessment of completeness was performed against the BUSCO3 eukaryotic gene set (Simão et al. 2015).

RNA Extraction, Sequencing, and Assembly

RNA from the embryonic brain tissues and mammary glands was collected from samples stored in RNAlater (Thermo Fisher). Adult tissues were collected from mastomys and processed immediately. All RNA was extracted using the Qiagen RNeasy micro kit. RNA libraries for the mouse and Mastomys embryonic mammary glands were made using the NuGen V2 RNA library kits and the other Mastomys libraries were created with the NuGen RNA-Seq Strand-Specific kits. For prostate, liver, and adult mammary tissues, RNA libraries were constructed with NEBNext Ultra RNA Library Prep Kit. All libraries were sequenced on a HiSeq 4000 2× 150 bp. De novo transcriptomes were created using each tissue following the Oyster River Protocol best practices (MacManes 2018). Briefly, reads were trimmed using Trimmomatic (Bolger et al. 2014) to remove adapters and no quality trimming was performed. After read correction with Rcorrector (Song and Florea 2015), multiple assemblies were performed using Trinity (Grabherr et al. 2011), Shannon (Kannan et al. 2016), and Spades (Chikhi and Medvedev 2014), and the assemblies were merged using OrthoFuser (MacManes 2018). Contigs were filtered by mapping reads back to the transcriptomes with Kallisto (Pimentel et al. 2016) and all contigs with TMP > 1 were annotated with dammit (Scott 2016).

Comparative Genomics

The mastomys genome was aligned to the Mus. musculus genome (mm10) and Rattus norvegicus (rn6) using Mercator (Dewey 2007). Pairwise genome alignments were made using LAST (Kiełbasa et al. 2011) between M. coucha and human (hg38), cat (felCat8), goat (Bickhart et al. 2017), mouse (mm10), rat (rn6), and beaver (Lok et al. 2017). These genomes were chosen due to their high quality and ability to represent specific clades. The pairwise MAF alignments were chained and netted into a mammalian multiple species alignment using ROAST (Blanchette et al. 2004). A neutral model of evolution was constructed using 4-fold degenerate sites with the Ensembl 89 annotation of mm10. The mammalian alignment was used to find regions of conservation with phastCons (Siepel 2005) ignoring mastomys. phastCons tree input models of neutral evolution and conservation were created using 4D sites mapped onto the alignment from Ensemble gene models and created with phyloFit (Siepel and Haussler 2004). Regions of mammalian conservation were examined in mastomys for accelerated evolution using phyloP and those with significant acceleration (FDR < 0.05) were used (Pollard et al. 2010). Annotation of MARs was performed relative to the mouse using regions with greater than 40% identity to the mouse mm10 assembly using liftOver. HOMER v4.9 (Heinz et al. 2010) was used to obtain genomic annotation classes with RefSeq and GO term enrichment was performed with GREAT 3.0 (McLean et al. 2010).

Population Diversity and History

In addition to the lab strain of M. coucha, low-coverage sequencing was performed on a wild specimen of M. coucha and of M. natalensis. Samples collected in South Africa were provided by the Durban Natural Science Museum and DNA from liver tissue was extracted using DNeasy Blood and tissue kit (Qiagen). Libraries were created by Novogene Co., Ltd using the NEB Next Ultra DNA Library Prep Kit and sequenced to 10× coverage on the Illumina HiSeq X platform with HiSeq X Ten Reagent Kit v2.5 for paired-end reads with 150 bp each (2× 150 bp). After mapping with BWA, single nucleotide polymorphisms (SNPs) were called using the GATK3.6 pipeline (Van der Auwera et al. 2013) following best practices for nonmodel organisms. Briefly, an initial SNP and indel set was called and hard filters were applied to obtain a preliminary variant set. This set was used to recalibrate read quality scores and variants were called again. The second variant set was again filtered with hard cutoffs to obtain a high-quality variant set of SNPs and indels for both M. coucha and M. natalensis.

Mouse Enhancer Assays

MAR456 was amplified using polymerase chain reaction (PCR) on mastomys or mouse genomic DNA using the following primers: mastomys forward, agatccaaatcacatgcagc; mastomys reverse, gctgcatgtgatttggatct; mouse forward, gaagtaggctggcacaagtaga; mouse reverse, acttgttgaaattttgctagttcttg. PCR products were cloned into the Hsp68-LacZ vector and sequence verified. All transgenic mice were generated by Cyagen Biosciences using standard procedures (Nagy et al. 2002) and harvested and stained for LacZ expression at either E11.5 or E13.5 as previously described (Pennacchio et al. 2006). Pictures were obtained using an M165FC stereo microscope.

Supplementary Material

Supplementary data are available at Molecular Biology and Evolution online.

Click here for additional data file.