Genomic Insights of Bruneian Malays

The Malays and their many sub-ethnic groups collectively make up one of the largest population groups in Southeast Asia. However, their genomes, especially those from Brunei, remain very much underrepresented and understudied. Here, we analysed the publicly available WGS and genotyping data of two and 39 Bruneian Malay individuals, respectively. NGS reads from the two individuals were first mapped against the GRCh38 human reference genome and their variants called. Of the total ∼5.28 million short nucleotide variants and indels identified, ∼217K of them were found to be novel; with some predicted to be deleterious and associated with risk factors of common non-communicable diseases in Brunei. Unmapped reads were next mapped against the recently reported novel Chinese and Japanese genomic contigs and de novo assembled. ∼227 Kbp genomic sequences missing in GRCh38 and a partial open reading frame encoding a potential novel small zinc finger protein were successfully discovered. Interestingly, although the Malays in Brunei, Singapore and Malaysia share >83% common variants, principal component and admixture analysis comparing the genetic structure of the local Malays against other Asian population groups suggested that they are genetically closer to some Filipino ethnic groups than the Malays in Malaysia and Singapore. Taken together, our work provides the first comprehensive insight into the genomes of the Bruneian Malay population.


Introduction
The Malays are an Austronesian ethnic group consisting of various subethnic populations residing mainly in the Southeast Asian countries of Brunei, Malaysia, Singapore and Indonesia. Despite them being one of the largest ethnic groups in the region, Malay genomes were not included in either the earlier HapMap or the later 1000 Genomes Project (1KGP) [1,2].
In fact, to date, only a handful of genomic studies have been conducted in Singapore, Malaysia and Indonesia [3][4][5][6][7][8][9][10]. Although valuable information in areas such as population genetic structure, rare or novel genetic variations, pharmacogenomics and genetic disease risk factors have been gained, much remain to be uncovered. This is especially true for the different Malay sub-ethnic groups in Brunei where their genomics landscape remained virtually unexplored till now.
With a population of about 430 thousand, Brunei is located on the northern running on default parameters. BCFtools was then used to filter the raw variants as previously described.
To minimise the limitations of each of the three variant callers [22,23], a consensus set of variants was generated by intersecting the three filtered VCF files using the "isec" command of the BCFtools. In addition, INDELs longer than 40 bp were also removed to further improve the confidence of the consensus variants.
On the other hand, variants of the other 37 Bruneian Malays were extracted from a PLINK-formatted genotyping dataset and transformed into a VCFformatted file using the program PLINK (version 1.9) [24]. A custom script was used to transform the variant tables from AncestryDNA ® to their respective VCF-formatted files. GRCh37-based genomic coordinates of all variants were finally lifted as described previously.
All variants called from the 41 Malays were next merged using BCFtools before being annotated against genomic features in various databases (Table S2) using the "annotate" command of BCFtools and AnnoVar [25].
PolyPhen-2 or classified as 'high impact' by ENSEMBL, were investigated by submitting them to the Gene Ontology online web-server (http://geneontology.org/). On the other hand, mitochondrial nonsynonymous variants were uploaded to the MITOMAP web-server (https://www.mitomap.org/) [26] and their functional impacts predicted using the APOGEE [27] function available there.
In addition, basic variant statistics on allele frequency, depth distribution, Het/Homo and Ts/Tv ratios were calculated using the statistical function of either BCFtools or RTG tool [28]. Variant density distributed across each chromosome was calculated based on a 1 Mbp sliding window using the SNPdensity command of VCFtool [29] and plotted using the CMplot R package [30].

CNV Calling
cn.MOPs was employed to identify CNVs in the two WGS datasets [31].

Uncovering Novel Sequences in the Bruneian Malay Genome
NGS reads from the two Bruneian Malays which failed to be mapped to the human reference genome were subjected to further investigation to search for novel sequences that may be unique to the Bruneian Malay genome. It is probable that some of these unmapped reads were in fact of microbial origins since salivary DNAs were used in the WGS. Therefore, two rounds of microbial mapping were applied to remove these microbial reads ( Figure   S2A). Firstly, they were BWA-mapped, as previously described, against the The presence of potential novel sequences in Bruneian Malay genomes was finally investigated by subjecting the final set of unmapped reads to a workflow involving mapping against East Asian sequences and de-novo assembly ( Figure S2B). A novel sequence is defined here as one that is missing from the human reference genome.
The recently reported ~12 and ~6 Mbp novel HX1 Chinese [12] and JRG Japanese [15] sequences, respectively, were chosen as the "reference genomes" for the mapping as comparative genomic analysis here showed that they are genetically the closest population groups, among those publicly available, to the Bruneian Malays. Again, a similar BWA mapping strategy used earlier was applied. However, the output BAM files were more stringently filtered (>Q30 mapping quality, >30× coverage depth, and mapped region spanning over the length of a read (i.e. >100 bp)) to increase the confidence that the mapped reads were indeed similar to the Chinese/Japanese sequences. Variant calling and filtering was next carried out as before using BCFtools' "mpileup"/"call" and "filter" (QUAL and GQ ≥30) commands. Consensus sequences representative of the Bruneian Malay were then generated using GATK's "FastaAlternateReferenceMaker" which replaced the reference alleles with the filtered variants.

Comparative Genomic Analysis
The Malay variant data files of (1) 89 and 25 genotyped individuals from SGVP and Morseburg's study, respectively, and (2) 96 whole genome sequenced individuals from SSMP were transformed, when necessary, and merged as previously described into two separate VCF files according to the country of origin. These were then compared against variants of the 41 Bruneian Malays using the "isec" command of BCFtools to identify shared and unique variants among the Malays in the three countries. Potential functional impacts of those that were unique to the Bruneian Malays were further investigated as previously described.
To gain more insights into the population genetic structure of the Malays in  (Table S3). As a control, the AncestryDNA ® 's genotyping variants of a European relative of the two mixed-race Bruneians were included. Again, whenever necessary, SNP files were transformed into GRCh38-based VCF files as previously described.
The PCA was conducted according to a method adopted from a previous study investigating the European population genetic structure [36]. Prior to calculating the principal components of these SNPs, a filtering step using PLINK was performed to increase the examining resolution of the SNPs while avoiding the removal of excessive number of true positives, i.e., SNPs with low linkage disequilibrium (LD). The filtering criteria included (1) missingness rate of ≥10%, i.e., below 90% genotyping rate to ensure all the SNPs were comparable across all the different population groups; and (2) genotypic r 2 , which is a common measurement for LD, ≥80% within a 50-SNP sliding window increasing at 5 SNPs per step to ensure the selected SNPs were all independent from one another as well as non-arbitrary. The first ten principal components of each individual with the filtered SNPs were then measured using the PLINK's command "pca" and those which could explain over 80% of the variance were selected and plotted against each other using the R package, ggplot2 [37].
An "unsupervised" admixture analysis on the same filtered SNP datasets was next conducted using the default settings of ADMIXTURE [38]. The analysis assumed no prior knowledge on ancestral origin and historical context of the observed genotypic profiles. The number of major ancestral origins (K) in each of the Asian population groups were first estimated and cross-validated before the hypothetical ancestral proportions (Q) of each of the 1,499 individuals were calculated at K=4, 5 and 6. A bar plot representing their hypothetical ancestral proportions was subsequently plotted using ggplot2.

Genomic landscape of the Bruneian Malays
The high-quality WGS paired-end reads from each of the two Bruneians could be mapped successfully across >94% the human reference genome at an average depth of 37x and 47x, respectively (Table S4). The remaining ~6% falling mainly within known gaps or challenging chromosomal regions, such as centromere and telomeres ( Figure S3).  (Table S7).

Novel Sequences in the Bruneian Malay Genome
A total of ~93 million unmapped or ~3% of the raw reads from the two WGS datasets were obtained after initial mapping against the human reference and microbial genomes. Of these, ~16 million were mapped to ~146 Kbp sequences making up of ~9.2 and ~137 Kbp novel Chinese and Japanese contigs, respectively. 170 short variants, including both SNPs and INDELs, were called from the mapping; suggesting the potential commonality of these sequences among the three population groups.
Furthermore, although more than 19K contigs with a cumulative length of ~53 Mbp could be assembled from the "leftover" reads, only 58 of them totalling ~82 Kbp were found to share >95% BLAST search sequence identity with either some regions of the GRCh38 human genome or fosmidsand BAC-cloned human sequences, i.e., they are of human origin. The "discarded" contigs were likely to belong to some other yet-to-be verified human or unculturable microbial or misassembled sequences.
The mapping and de-novo assembly approaches had, therefore, yielded 227,763 bp novel Bruneian sequences which are missing in the GRCh38 human reference genome. When the gene-coding potential of these sequences were further investigated, an open reading frame out of a total 473 was found to share 92% sequence identity across a 106-residue segment with a macaque's hypothetical small zinc finger protein ( Figure   S4). The protein is predicted to be involved in the homeostasis of zinc ions.

Distinct admixed genetic structure of the Bruneian Malays.
To gain a better understanding on the population genetics of Bruneian Malays in the context of the wider Asian population groups, the first ten While K=5 has the lowest cross-validation error, K=6 appears to be in best agreement with findings from the PCA ( Figure 3B).

Genomic Landscape of Bruneian Malays
The ~5.28 million variants identified in our study have allowed, for the first time, the genomic landscape of Bruneian Malays to be explored in greater details and made comparative genomic studies against other Malay groups in the region possible. Although more than 216K of these variants have never been reported in dbSNP, the slightly lower Het/Homo ratio observed here, when compared to other population groups, seems to suggest lower heterozygosity in the genetics of the Dusun Malays than other Malay groups in the region. Higher homozygosity has in fact been observed in smaller population and less admixed groups [41], which is likely true for the Bruneian Dusun Malays. Crucially, many of the Bruneian variants were found to have minor allele frequency <0.001. This subset of rare variants is likely to be unique to the Bruneian Malays. Similar pattern has also been observed in studies assessing the allele frequency of variants in different population groups [42]. Although most of these rare variants are expected to have benign functional impacts, the existence of population-specific disease risk alleles should not be ruled out. In fact, ~200 of the novel variants were found to be located on exonic regions of each of the Bruneian genome and majority of them are predicted to be non-synonymous which may impact genes known to be associated with genetic diseases, such as cancer, cardiovascular diseases, and diabetes mellitus, and pharmacogenetic markers. However, more genotypic and phenotypic data from larger cohorts will be needed to provide the necessary resolution power required in studies such as GWAS to establish the likely association between these rare variants and diseases.
Although majority of the variants were distributed fairly evenly across the genome, several variant-dense loci which have not been reported in other population groups were identified in the two Bruneian Malay individuals. For instance, chr10 and 13 were found to harbour the highest variant density whereas findings from 1KGP reported chr16 and 1 as having the highest and lowest SNP density, respectively [39]. However, it is important to take into consideration potential biases of such a comparison between two local individuals against the majority European samples in 1KGP.
In addition, a number of genes are predicted to be disrupted by the CNVs identified in the two local Malays. Although it may be challenging to reconstruct the exact haplotype structure of these CNVs using the data from short read sequencing, such disruptions do have the potential to produce new protein isoforms which may or may not be functional. In fact, 13 haplotypes/alleles of the HLA-H pseudogene were reported elsewhere to reach 19.6% in East Asian populations and they have been shown to be functional [43]. Therefore, further investigation of these CNV-impacted genes using such technology as long read sequencing may be warranted.

Novel sequences in the Bruneian Malay genome
Although the novel sequences identified here have not been validated, for example, by Sanger sequencing or PCR, the facts that they do share high sequence identity to partial human reference genome and fosmid-or BAC-  [15] were also found in the de-novo assembled Korean genome by Seo et al. (2016) [16].
Since both Chinese and Japanese studies have reported novel sequences in the Mbp range, it is probable that more Malay-specific genomic sequences have yet to be uncovered. In addition, the finding of an ORF encoding a novel human homologue of a primate's small zinc finger protein have added more weight to the importance of finding such populationspecific novel sequences.

Subtle Genetic Differences among Southeast Asian Malays
The Malays are one of the largest Austronesian population groups spreading over Island Southeast Asia and as far away as South Africa.
Although the different sub-ethnic groups share considerable cultural and linguistic ties, they are not genetically homogenous [6,11,44]. Adding to this genetic diversity, population genetic structure analysis conducted here has, for the first time, unveiled the differences in genetic ancestries among shown the existence of subtle differences in the population genetic structure among the different Malay groups in Southeast Asia. Hence, a more refined stratification of these groups using variants from larger cohorts will be necessary should the benefits of medical genetics and genomics are to be fully realised in the region.