Fontsere et al. BMC Genomics (2021) 22:735 https://doi.org/10.1186/s12864-021-08025-y RESEARCH Open Access The genetic impact of an Ebola outbreak on a wild gorilla population Claudia Fontsere1*†, Peter Frandsen2,3†, Jessica Hernandez-Rodriguez1,4, Jonas Niemann5, Camilla Hjorth Scharff-Olsen5, Dominique Vallet6, Pascaline Le Gouar6, Nelly Ménard6, Arcadi Navarro1,7,8,9, Hans R. Siegismund3, Christina Hvilsom2, M. Thomas P. Gilbert5,10, Martin Kuhlwilm1,11*†, David Hughes1,12,13† and Tomas Marques-Bonet1,7,14,15*† Abstract Background: Numerous Ebola virus outbreaks have occurred in Equatorial Africa over the past decades. Besides human fatalities, gorillas and chimpanzees have also succumbed to the fatal virus. The 2004 outbreak at the Odzala- Kokoua National Park (Republic of Congo) alone caused a severe decline in the resident western lowland gorilla (Gorilla gorilla gorilla) population, with a 95% mortality rate. Here, we explore the immediate genetic impact of the Ebola outbreak in the western lowland gorilla population. Results: Associations with survivorship were evaluated by utilizing DNA obtained from fecal samples from 16 gorilla individuals declared missing after the outbreak (non-survivors) and 15 individuals observed before and after the epidemic (survivors). We used a target enrichment approach to capture the sequences of 123 genes previously associated with immunology and Ebola virus resistance and additionally analyzed the gut microbiome which could influence the survival after an infection. Our results indicate no changes in the population genetic diversity before and after the Ebola outbreak, and no significant differences in microbial community composition between survivors and non-survivors. However, and despite the low power for an association analysis, we do detect six nominally significant missense mutations in four genes that might be candidate variants associated with an increased chance of survival. Conclusion: This study offers the first insight to the genetics of a wild great ape population before and after an Ebola outbreak using target capture experiments from fecal samples, and presents a list of candidate loci that may have facilitated their survival. Keywords: Ebola, gorilla, non-invasive samples, candidate genes Background The Ebola virus (EBOV), discovered in 1976, causes a severe disease and often fatal hemorrhagic fever for which numerous human outbreaks have been reported throughout Africa [1]. The most virulent outbreak re- * Correspondence: claudia.fontsere@upf.edu; martin.kuhlwilm@univie.ac.at; tomas.marques@upf.edu ported to date was in West Africa in December 2013 † Claudia Fontsere and Peter Frandsen share co-first authorship. and lasted until 2016 with more than 28,000 confirmed † Martin Kuhlwilm, David Hughes and Tomas Marques-Bonet share co-last or suspected human cases and more than 11,000 human authorship. 1 Institut de Biologia Evolutiva, (CSIC-Universitat Pompeu Fabra), PRBB, Doctor deaths [2]. Since then, other outbreaks of Ebola have Aiguader 88, 08003 Barcelona, Catalonia, Spain been observed. In June 2020, when the 2018 outbreak Full list of author information is available at the end of the article © The Author(s). 2021 Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data. Fontsere et al. BMC Genomics (2021) 22:735 Page 2 of 12 that occurred in Ituri, North Kivu and South Kivu prov- marking the relevance of social dynamics for transmis- inces of Democratic Republic of Congo was declared sion [13]. A previous study using 17 microsatellites over by the World Health Organization (WHO), 3470 found no loss of genetic diversity after one EBOV out- cases had been reported with 2287 deaths (fatality rate break in Lossi sanctuary and Odzala-Kokoua National of 66%) [3]. Park, which could be explained by post-epidemic immi- EBOV belongs to the single-stranded RNA virus family gration, sufficiently large remnant effective population Filoviridae [4] with five distinct strains in the Ebola size or a short period of time after the decline [18]. The genus: Zaire, Sudan, Bundibugyo, Taï Forest and Reston. present study represents a continuation of this afore- The first three are responsible for the majority of human mentioned research since many aspects of EBOV infec- infections [5, 6]. The virus is highly infectious and can tion in wild gorilla populations are not yet explored, enter the body through direct contact of broken skin or such as genotypes of surviving individuals that might mucous membranes with infected blood or body fluids, contribute to resistance or higher chance of survival to causing symptoms including fever, vomiting, diarrhea, EVD. internal and external bleeding. Ebola hemorrhagic fever Furthermore, microbial organisms inhabiting the gut or Ebola virus disease (EVD) is an acute and severe dis- also play a potentially crucial role in training and main- ease with a fatality rate in humans around 50% [5–7]. taining the immune system [19–21]. While some com- Infectious diseases such as Ebola are considered to be mensal microbes are associated with priming the a threat to the survival of African great apes [8], together immune response or activating antiviral responses, with other threats such as habitat loss, climate change others facilitate the development of the infection or sup- and poaching [9]. In some cases, the human outbreaks press the immune response [22]. For instance, a recent have been linked to contact with infected bushmeat study reported a link between the gastrointestinal micro- from gorillas or chimpanzees [10] and several surveys biome of healthy humans and a predisposition to severe have reported dramatic declines in populations of great COVID-19 [23], with an abundance of Klebsiella, apes in parallel with human EVD outbreaks with labora- Streptococcus, and Ruminococcus being correlated with tory confirmation of Ebola virus infection in some car- elevated levels of proinflammatory cytokine. The associ- casses [10–12]. Gorilla populations from the Republic of ation between infectious diseases and the gut micro- Congo suffered severe die-offs during a human EVD biome of nonhuman primates is less well understood. outbreak near the Lossi sanctuary in 2002–2003 [12] While the gut microbiome of SIVgor-infected (gorilla and Odzala-Kokoua National Park in 2004 [13] with re- Simian Immunodeficiency virus) wild gorillas seems to ported mortality rates as high as 95%. In Lossi sanctuary be more robust to dysbiosis than those of chimpanzees alone, it was estimated that the Ebola virus killed 5000 and humans [24], it is unclear how shifts in the gorilla wild gorillas [12]. The severe population decline, as a re- gut microbiome can impact the severity of viral infec- sult of illegal hunting, disease and habitat loss, has con- tions, and in particular the immune response to Ebola tributed to the 2007 shift of the conservation status of virus infection. western gorillas from “endangered” to “critically endan- Long-term monitoring of a western lowland gorilla gered” by the International Union for Conservation of (Gorilla gorilla gorilla) population in Odzala-Kokoua Nature (IUCN) [14]. Furthermore, the recent outbreak National Park (Republic of Congo) has been ongoing of human EVD in North Kivu (Democratic Republic since 2001 until 2017 encompassing the Ebola outbreak Congo) [3] was in close proximity to the remnant popu- of 2004 [16–18, 25] which resulted in a mortality rate of lations of eastern gorilla species, hence a human to ape 95% [13]. Population monitoring involved the recording transmission of the virus could potentially mark the end of individual histories of hundreds of identified individ- of existence for this critically endangered species. uals, the determination of sex, age and social status and Threats such as infectious diseases are relevant for the collection of fecal samples in different time periods. conservation efforts, and efficient strategies are needed During the Ebola outbreak, which lasted six months in to reduce the effects of EVD on wild great ape popula- the monitored population, about 340 individuals disap- tions [15]. Understanding any genetic impact that EBOV peared and entire groups were missing. Knowing the oc- outbreaks might have on wild populations is vital, as currence of Ebola epidemics in humans, the occurrences EVD contributes to the fragmentation of gorilla popula- of gorilla carcasses infected by Ebola [10], and the detec- tions due to a heterogeneous spatial influence of the tion of Ebola virus antibodies in great ape fecal samples outbreak [16]. Social dynamics in gorillas are rapidly af- in the areas of gorilla population decline [26], the scien- fected by Ebola through a decrease in social cohesion, al- tific community has attributed the abnormal rate of though recovery after the outbreak has been observed gorilla disappearances to massive death caused by EVD. [17, 18]. One study reported that solitary individuals Individuals without any sign of injury, handicap or sen- were less affected than individuals living in groups, escence that went missing after the Ebola period were Fontsere et al. BMC Genomics (2021) 22:735 Page 3 of 12 considered as dead because of EVD, and individuals that samples were from gorillas identified before the epi- were seen before and after the outbreak were considered demic and still observed after the epidemic, and are as surviving the epidemic period. identified henceforth as ‘survivors’ [17]. Although no im- Here, focusing on a well described western lowland munological test could be performed on the dried fecal gorilla population [18, 27], we have obtained genetic samples, the survivors were known to live in the area data from gorilla fecal samples pre- and post-outbreak where the epidemic emerged and thus they were sus- to explore potential associations to survivorship. Particu- pected to have been in contact with the virus. These 15 larly, we compare the genetic variation at the single nu- samples were collected between 2005 and 2014 (Supple- cleotide level in 123 autosomal target-captured genes mentary Table S1). We used target capture enrichment with putative roles in virus immune response and the to sequence the genomic regions of 123 genes, which gut microbiome composition in a reduced panel of sur- had previous evidence of putative roles in immune re- vivors and non-survivors of this Zaire EBOV outbreak. sponse to EBOV or other viruses (Supplementary Table With that we show that targeted capture on non- S2). In addition, 15 neutral regions previously studied in invasive fecal samples and next-generation sequencing other human and non-human primate studies were also can be used to study the impact of this severe disease in targeted [28, 29] (Supplementary Table S2). Target de- a natural population. sign and all analyses were performed using the human reference genome due to the higher quality of annota- Results tions compared to the gorilla reference genome. We se- A total of 31 non-invasive fecal samples from identified quenced an average of 73 million paired reads per western lowland gorillas were collected between 2001 individual, 12% of which were unique (Supplementary and 2014 in Odzala-Kokoua National Park, Congo [18] Table S3). On average, 0.38% of the data mapped to the (Fig. 1A). Sixteen of these were from individuals de- target space, representing an on-target effective coverage clared missing after the Zaire Ebola virus outbreak in of 53.89-fold (range: 2.52-fold to 230.70-fold; Fig. 1B, 2004. These individuals did not show any sign of injury, Supplementary Table S3), with 72% of the target space handicap or senescence that could lead to natural death. covered by at least 4 reads per individual. Samples G282, Also, they were observed shortly before the epidemic G1392 and G638 performed poorly, with <50% of the period and were thus unlikely to disperse. Therefore, target space covered at a minimum depth of 4 reads they were suspected to have died because of the infec- (Supplementary Fig. S1). Overall performance can be tion, here termed ‘non-survivors’. These samples were assessed by calculating how well the capture resequen- collected between 2001 and 2004. The remaining 15 cing experiment went relative to expectations had we Fig. 1 Sample description. A) Geographical map of the extant range of gorillas and the Odzala-Kokoua National Park (Republic of Congo) where fecal samples were collected between 2001 and 2014, overlapping the Ebola outbreak in 2004. B) Average coverage reached in the target space per sample in both studied groups. C) Percentage of human contamination in each fecal sample Fontsere et al. BMC Genomics (2021) 22:735 Page 4 of 12 performed random shotgun sequencing. In that regard, we observe an average enrichment of 125-fold (88–346- fold) (Supplementary Table S3). Individuals with ex- tremely low proportions of target space covered by at least 4 reads (<30%) (G638 and G282) and high hetero- zygosity and high levels of human contamination were removed from further analysis (individuals G374 and G1392, Fig. 1C and Supplementary Fig. S2, Table S3). With the final dataset of 13 survivors and 14 non- survivors, we validated that the genotype information obtained was in concordance with previously published gorilla whole-genomes, as determined by a principal component analysis (PCA) (Supplementary Fig. S3). Since the sampling of the gorilla population was done to avoid related individuals by sampling from different groups (a total of 45 groups and 31 solitary males were present before the outbreak), we identified a single case of close relatedness (1st degree) among them (non-sur- vivor G374 and survivor G739; kinship coefficient = 0.4) Fig. 2 Genetic distance dendrogram and heterozygosity among (Supplementary Fig. S4). In addition, we observed no non-survivor and survivor gorilla groups. A) Clustering dendrogram stratification correlating with survivor/non-survivor clas- of pairwise genetic distance derived from genotype likelihoods (N = sification. Individuals appear to be dispersed randomly 5477). B) Mean heterozygosity (bp−1) in non-survivors and survivors; not significantly different (Student t-test, p-value = 0.34) across a dendrogram derived from shared genotype like- lihood dosage states (Fig. 2A), and a univariate linear re- gression of survivorship on the top 5 PCs identified no work with low coverage data that relies on genotype significant structure associated to survivorship (P-value likelihoods. With this method, we were able to recover >0.1; Supplementary Fig. S5). Hence, we determined no the majority of the genes found above (30 out of 36). genome-wide group structure differences between survi- However, ANGSD returned more hits and thus more vors and non-survivors. The overall level of genetic di- genes (Supplementary Figs. S7, S8 and Table S5), ren- versity within the target space of the studied gorillas dering the above approach more conservative. was, on average, lower than that of western lowland go- Next, we explored the potential functional impact of the rillas obtained from whole-genome sequencing (Supple- variants in the nominally significant candidate loci that mentary Fig. S6) [30, 31], an expected outcome differentiate survivors from non-survivors. While we following the target capture procedure. Moreover, there found no significant associations with gene ontology cat- are no statistically significant differences in heterozygos- egories, the analysis of predictions of functional conse- ity between survivors and non-survivors (Student’s t- quences pinpointed six missense mutations which differed test, p-value = 0.34; Fig. 2B). in frequency between the two groups (Supplementary In order to determine genetic differences between the Table S6), one each in the ATM, IGKV4–1 and RNF167 groups, we calculated three summary statistics on a genes (all lower in survivors) and three in the ACTB dataset of 6852 high-quality variants: (1) the difference (Actin Beta) gene (one unique to survivors, two lower in in allele frequency (ΔFrequency), (2) the fixation index survivors, Supplementary Table S4). All three missense (FST), and (3) the significance level (α) of each variant variants in ACTB are predicted as deleterious by both the for its association with the binary trait survivor/non-sur- PolyPhen [33] and SIFT [34] algorithms. Furthermore, the vivor (Fig. 3). We found 118 SNPs within the target derived variant in survivors in the immunoglobuline- space that surpassed our α threshold in the association encoding IGKV4–1 might be deleterious (C-score > 20 test, and we also reported their ΔFrequency and FST [35]), hence potentially functionally relevant. Since we values. However, after controlling for type I error none used the human genome for target design, mapping and of these remained significant. Out of these, seven genes variant calling, we caution that differences in exon usage have multiple nominally significant SNPs (CD1B, or pseudogenization on the gorilla lineage might confound IGKV4–1, HLA-A, ACTB, LYN, CD68 and MX1), while these inferences of protein-coding changes. In order to 10 neutral regions (~10 kb each) have at least 1 nomin- confirm the expression of these genes, and specifically the ally significant variant (Supplementary Table S4). For exons of interest, we mapped transcriptome data from six comparative validation, we repeated the association ana- tissues in gorillas [36] to the same reference genome, and lysis using ANGSD [32], a software explicitly built to quantified expression levels. We confirmed the expression Fontsere et al. BMC Genomics (2021) 22:735 Page 5 of 12 Fig. 3 Association analysis to detect SNPs and candidate genes related to survivorship to EBOV outbreak. A) Significance level (threshold set at α = 0.05, −log10(P-value) = 1.30). B) Difference in allele frequency (threshold set at ±0.2). C) Fixation index (threshold set at FST = 0.15). Dashed lines delineate the thresholds used of these genes and found high transcript abundance (log2- the nature of the samples, we first analyzed the value of counts > = 9) for ACTB, RNF167 and ATM, while microbiota present in each fecal extraction using a IGKV4–1 was only detected at low levels in these tissues 16S rRNA library (Methods). We obtained a total of (log2-value of counts = 5.84). We found ~10,000 RNA se- 96,928,400 reads with an average sequencing depth of quencing reads overlapping the three loci of interest 1,101,459 (SD ± 418,989) reads per gorilla fecal ex- in ACTB, 437 in RNF167, and 96 in IGKV4–1, sup- traction (Supplementary Table S8), and determined porting the expression of these specific loci in gorilla the abundance of taxa (Supplementary Figs. S9 and tissues, while for ATM only 10 reads overlapped S10). Firmicutes (53.79%), Bacteroidetes (12.02%) and (Supplementary Table S6). Chloroflexi (11.11%) were the predominant phyla Among the 118 significant SNPs, we found no direct (Supplementary Fig. S9 and Table S9), that include overlap with loci associated to 2385 traits in genome- the following most abundant orders: Clostridiales wide association studies (GWAS) for humans [37]. How- (39.37%), Bacillales (11.24%), Bacteroidales (10.87%) ever, we found 49 associated loci within close proximity and Anaerolineales (11.11%) (Supplementary Fig. S10 (5000 bp) to GWAS loci (Supplementary Table S7), and Table S10), concordant with previous findings among which the locus 19:1104078 near GPX4 stands [38, 39]. We found no taxa significantly differing in out for its association with five blood cell traits (counts relative abundance between survivor and non-survivor and percentages of leukocyte cell types). Furthermore, 6: gorillas (Bonferroni-corrected p-values >0.05; Supple- 29916885 in the HLA is associated to hemoglobin levels. mentary Fig. S11 and Table S11), and sample groups Five loci in HLA-DRA are associated with different auto- were not separated in a clustering analysis (Fig. S12). immune diseases like systemic lupus erythematosus or Since these results on the gut microbiome diversity multiple sclerosis (Supplementary Table S7). Further- did not support differences between both gorilla groups, more, nine loci in this region (HLA-DRB5 and HLA-A) we decided to perform deep sequencing on the fecal li- are associated with schizophrenia or autism spectrum braries. We generated a total of 801,132,281 sequences disorder (Supplementary Table S7). from DNA libraries (4,025,054-25,593,317 reads per It has been reported that intestinal microbiota play sample; Supplementary Table S8) and used MALT a critical role in immune response to infectious dis- (MEGAN Alignment Tool) [40] for an alternative eases [19–21]. Thus, the microbiome might be rele- characterization of the microbial profile of the samples vant for the survival of wild gorilla populations (Supplementary Fig. S13). We find that the majority of experiencing an Ebola outbreak. Taking advantage of the identified taxa are associated with the gut Fontsere et al. BMC Genomics (2021) 22:735 Page 6 of 12 microbiome (Supplementary Table S12 and S13). By far The latter is particularly true (2a) relative to the tim- the most abundant taxon is the gut bacterium Escheri- ing of any exposure, but also in respect to the (2b) chia coli, which could be detected in all samples, and dynamic nature of the gut microbiome [39]. makes up more than 50% of all assigned sequences in Given the limited sample size, we developed an ap- nine of the samples. Also in high abundance are species proach using differences in allele frequency, the fixation of the Bacteroidales order, such as Bacteroides cellulosi- index and the effect size to determine variants most lyticus and Prevotella spp., as well as members of the strongly associated to survivability in the studied popula- Clostridiales, Lactobacillales, and Bacillales orders, cor- tion, generally replicable using an association analysis roborating previous reports on the composition of west- with ANGSD. While 44 of the 118 nominally significant ern lowland gorilla gut microbiomes [38, 41]. SNPs (Supplementary Table S4) do fall within 10 of the Furthermore, we detected pathogenic taxa in high abun- 15 neutral regions included in the study, some SNPs dance in some of the samples, such as Clostridium botu- might be functionally relevant for surviving the EBOV linum, Acinetobacter baumannii, and Klebsiella outbreak. The non-synonymous variants in ACTB, pneumoniae, which have been previously found in the RNF167 and IGKV4–1 genes are obvious candidate loci, gorilla gut [42]. However, the microbial profiles of survi- and particularly the three deleterious missense muta- vors and non-survivors do not differ significantly from tions in Actin Beta appear to be strong candidates for a each other in this analysis either (Bonferroni-corrected higher survival rate. The actin cytoskeleton is important p-values >0.05, two-sided t-test, Supplementary Table for virus assembly [45], and a disturbed assembly S11), and the two groups do not form separate clusters process could have influenced the viral load in individ- in a Principal Coordinate analysis or a Neighbor Joining uals with changes in this protein. As expected for a gene Tree (Supplementary Figs. S14 and S15). encoding a structural protein, ACTB is highly expressed in gorilla tissues. Furthermore, the variant in IGKV4–1 Discussion might improve the immune response to viral infection We investigated non-invasive fecal samples from a through antigen recognition [46]. The missense muta- long-term monitored population of western lowland tion in ATM, which belongs to the PI3-kinase family, gorillas in the Republic of Congo, including individ- could interfere with the cellular entry specifically of the uals that most likely succumbed to the Zaire Ebola Ebola virus [47], although we could not confirm expres- virus outbreak in 2004, as well as surviving individ- sion of this locus in vivo in the available tissues. uals [17, 18]. We used targeted capture of 123 auto- We find other potentially relevant non-coding vari- somal genes with putative roles in immune response ants, 49 of which are in close proximity to SNPs as- to EBOV or other viruses (Supplementary Table S2) sociated to GWAS traits in humans, suggesting from fecal samples. This yielded an enrichment of possible regulatory functions (Supplementary Table more than 100-fold across samples, and a medium to S7). Among those, the association of GPX4 with high coverage of the target space across most individ- leukocyte cell type count might reflect differences in uals (Supplementary Fig. S1 and Table S3). Although leukocyte composition after viral infection. Differences a large proportion of reads were duplicates, the over- in hematocrit or hemoglobin levels might have con- all performance was high and these results demon- tributed to the survival of wild gorillas considering strate the great potential of capture experiments for that hemorrhage and internal bleeding are symptoms obtaining genotypes from fecal samples of wild great of Ebola infection. Since eight loci are associated to ape populations [43, 44], for which high-coverage se- the HLA-DRB gene, a direct involvement in the adap- quencing would be prohibitively expensive. We deter- tive immune system might cause the signature ob- mined that the studied individuals were not closely served at this locus, particularly given that human related, hence most likely representing a random sam- survivors of EVD show a lower frequency of HLA- pling of the wild gorilla population before and after DR-positive T cells [48]. the outbreak. We also investigated the microbial com- munity composition of survivors and non-survivors, Conclusion finding no significant differences in taxa abundance, By using fecal samples and targeted capture enrichment, neither using 16S rRNA or deep sequencing data non-invasive assessment of numerous individuals from (Supplementary Figs. S11, S13 and Table S11). Hence, wild populations is possible. Here, we demonstrate that we find no evidence that the gut microbiome of indi- this approach can be used to analyze temporal genetic viduals has an influence on the survival rate of wild changes in wild great ape populations in response to en- gorillas exposed to Ebola. However, these observations vironmental factors. Additionally, we present candidate are limited by (1) the sample size, and (2) the broad loci that may have facilitated the survival of gorilla indi- range of collection dates (Supplementary Table 1). viduals or groups after an outbreak of the Zaire Ebola Fontsere et al. BMC Genomics (2021) 22:735 Page 7 of 12 virus. Understanding putative adaptive responses to this and reverse primer (IS2_P7_16S_806rB: GTGACTGGAG pathogen in wild populations can help to advance our TTCAGACGTGTGCTCTTCCGATCTGGACTACNVG knowledge on the natural dynamics of this severe dis- GGTWTCTAAT) include the complementary sequences ease. Such a strategy might be useful in a broader con- necessary for the final indexing step [52]. Protocols are text, since these and other primates are susceptible to provided in [53]. other infectious diseases such as Covid-19 [49]. Target Design, Capture and Sequencing Methods RNA baits covering the target space were designed and Samples, DNA extraction, library preparation synthesized by Agilent with a minimum of 3x bait cover- Non-invasive fecal samples from western lowland go- age. The target space included specific autosomal genes rillas were collected between 2001 and 2014 near the (123 genes) and 15 neutral regions (~10 kb each) (Sup- Lokoué clearing in Odzala-Kokoua National Park, Re- plementary Table S2). For target enrichment, the fDNA public of Congo [18], following the standard practices in libraries were pooled into one equimolar batch and sub- the field. The population and group structure of this jected to two consecutive rounds of DNA capture with western lowland gorilla population has been described the RNA baits in 8 hybridizations. Captured fDNA li- previously [27]. Among collected fecal samples, we se- braries were sequenced on the Illumina system in four lected 31 samples from previously identified individuals HiSeq 2500 2x125 lanes and one HiSeq 2500 rapid run each from unique social groups (Supplementary Table at 2x250 bp. The 16S libraries were sequenced on one S1). Sixteen of those individuals were declared missing HiSeq 2500 rapid run 2x250 bp lane. In addition, we after the epidemic in 2004 [13]. Given the high mortality generated paired-end sequences from the fDNA libraries rate (95%) of Ebola virus in gorillas, they were presumed on four HiSeq 4000 lanes (2x150bp) to study the whole dead during the time span of the epidemic and were microbiome composition (Table S8). identified here as non-survivors. Fifteen individuals were observed before and after the epidemic, and were de- Mapping and Variant Discovery scribed here as survivors (Fig. 1A and Supplementary Prior to mapping, paired-end reads belonging to the same Table S1) [17]. Fresh samples were collected by field in- library but sequenced in different lanes were merged into vestigators wearing masks and gloves, and were dried a single FASTQ file. PCR duplicates were directly re- with silica beads and then stored at room temperature moved from FASTQ files using FASTuniq (v1.1) [57]. until arrival in a laboratory where they were stored at Overlapping reads were merged (minimum overlap of 10 4 °C until extraction. Sampling dates are recorded in bp, minimum length of final read to 50 bp) using PEAR Supplementary Table S1. (v0.9.6) [58]. Reads were mapped using BWA mem DNA was extracted from 10 mg of dried sample using (v0.7.12) [59] to the human reference genome Hg19 the 2CTAB/PCI protocol [50] using negative controls that (GRCh37 from the UCSC database). Assembled reads were checked for DNA contamination before subsequent were mapped considering single-end specifications and experiment. Three different extractions were carried out, unassembled reads considering paired-end specifications. except for samples G778, G374, G344, G498 and G372, Any remaining PCR duplicates were removed using where only two extractions were performed (Supplemen- PicardTools MarkDuplicates (v1.95) (http://broadinstitute. tary Table S1). A DNA library [51–53] and a 16S rRNA li- github.io/picard/). Non-primary alignments and reads brary [54–56] were prepared for each extract. Isolated with quality below 30 were filtered from the dataset with DNA samples were quantified with Qubit with a mean es- samtools (v1.5) [60]. Finally, single-end and paired-end timated concentration of 13.3 ng/μl (range: 0.90–74.7). reads were merged into a single BAM file using Picard- Whenever possible, a total of 250 ng of DNA was used to Tools MergeSamFiles (http://broadinstitute.github.io/ construct DNA libraries, but never more than a total vol- picard/). The percentage of aligned reads for each DNA ume of 33 μl was taken from any single sample. DNA was extraction and sample was calculated by dividing the sheared with a Covaris S2 instrument and 88 fecal DNA number of uniquely and high-quality mapped reads (with- (fDNA) libraries were prepared following a custom dual- out duplicates) by the total number of sequenced reads. indexing protocol with 25 cycles of amplification [51, 52]. The percentage of on-target aligned reads was calculated Subsequent to DNA library preparation, 88 16S rRNA li- for each sample by dividing the number of on-target fil- braries were prepared using 1 μl of total DNA. The V3 tered reads by the number of sequenced reads. The aver- and V4 regions of 16S rRNA were target amplified using age target effective coverage was calculated dividing the modified 341F and 806Rb primers [54–56], incorporated number of aligned bases by the total length of the targeted into the dual-indexing protocol [52]. The forward primer genomic space. Finally, the enrichment factor (ER) of the (IS1_P5_16S_341f: ACACTCTTTCCCTACACGACGC capture performance was calculated using the ratio be- TCTTCCGATCTNNNNCCTACGGGNGGCWGCAG), tween the on-target reads by the total mapped reads over Fontsere et al. BMC Genomics (2021) 22:735 Page 8 of 12 the targeted size by genomic size (ER = (On-Target representative of each know gorilla subspecies [30, 31]. Reads/Mapped Reads)/(Target Size/Genome Size)). The We also obtained a PCA using the GATK genotype calls coverage for each target region was retrieved using SAM- after keeping only variants with minor allele frequency tools bedcov [60]. of 0.02 with plink --pca option (N = 6051) [65]. For variant calling, all BAM files belonging to the same sample were merged into a single BAM file using Picard- Genetic distance, relatedness and heterozygosity Tools MergeSamFiles (v1.95) (http://broadinstitute. We used the genotype likelihood information for the github.io/picard/). Variant discovery was performed studied individuals to obtain the genetic distance by run- using GATK ‘Unified Genotyper’ [61] for each sample ning ngsDist in ANGSD [66] with the following parame- independently with the following parameters -out_mode ters: --n_sites 5477 --probs TRUE --pairwise_del. Then, EMIT_ALL_SITES -stand_call_conf 5.0 -stand_emit_ we constructed an Euclidean distance matrix based on conf 5.0 -A BaseCounts -A GCContent -A RMSMap- the genotypes and performed a hierarchical clustering pingQuality -A BaseQualityRankSumTest. Afterwards, using the R package ape [67]. We also run PCAngsd we merged each sample gvcf to a single one using [64] considering only the study gorillas to discard any GATK `CombineVariants’ [61] with the following pa- possible intra-group structure. rameters -genotypeMergeOptions UNIQUIFY –exclude- The theta coefficients of kinship (probability of a pair NonVariant. We also included in the gvcf the genotype of randomly sampled homologous alleles are identical by information of available whole genome data of six descent) were calculated using the NgsrelateV2 [68, 69] Gorilla beringei beringei, eight Gorilla beringei graueri, on the genotype likelihood obtained with ANGSD [32]. one Gorilla gorilla dielhi, and twenty-three Gorilla Note that all possible genotype likelihoods, even outside gorilla gorilla samples [30, 31]. The VCF was filtered the target space (N = 226,094), were used since the with VCFtools [62] to keep only biallelic positions with coverage of the kinship and neutral markers was DP >3 and quality >30 and without indels. insufficient. Genotype likelihoods were directly obtained from To assess global levels of heterozygosity, the unfolded BAM files with ANGSD [32] including four Gorilla ber- SFS (Site Frequency Spectrum) was calculated for each ingei beringei, four Gorilla beringei graueri, one Gorilla sample separately, including thirteen gorilla whole- gorilla dielhi, and four Gorilla gorilla gorilla, with the genomes representative of all gorilla subspecies [30, 31], following parameters and only in the target space: using ANGSD [32] and realSFS [70] only in the target -uniqueOnly 1 -remove_bads 1 -only_proper_pairs 1 space with the following quality filter parameters: -trim 0 -C 50 -baq 1 -minInd 21 -skipTriallelic 1 -GL 2 -uniqueOnly 1 -remove_bads 1 -only_proper_pairs 1 -minMapQ 30 -doGlf 2 -doMajorMinor 1 -doMaf 2 -trim 0 -C 50 -baq 1 -minMapQ 20 -minQ 20 -setMax- -minMaf 0.05 -SNP_pval 1e-6. Depth 200 -doCounts 1 -GL 1 -doSaf 1. We used the human genome (Hg19) to determine the ancestral state. Quality control We evaluated the amount of human contamination in Association analysis each fecal library using the HuConTest script [63], as The genotype calls obtained with GATK were further fil- described previously [44]. The majority of samples have tered with Plink [65] to exclude variants considering less than 2% of human contamination, but samples their missing rate (−−geno 0.05), minor allele frequency G348 and G1392 have estimates of human contamin- (−−maf 0.01) and Hardy-Weinberg equilibrium (−−hwe ation of 6.7 and 25.1%, respectively (Supplementary 0.00001). The final dataset consists of 27 samples (13 Table S3). These two samples also show extreme values survivors and 14 non-survivors) and 6852 high-quality of heterozygosity (deviating >1 s.d. from mean heterozy- variants. Nominal significance was set at an alpha of gosity; Fig. S2). For the identification of individuals and 0.05 and a Bonferroni corrected alpha threshold of markers with elevated missing data rates we used the 7.3x10–6 (0.05/6852) was defined to account for family- proportion of the target space covered by at least 4 wise error. Associations were tested for by a chi-square reads. Individuals with less than 30% of covered target allelic test with one degree of freedom and p-values were space (4 reads) were not used for further analysis (G638 estimated by permutation in plink (plink –assoc – and G282). mperm 10,000) [65].. The p-values were plotted in a A principal component analysis (PCA) was performed Manhattan plot in R (v3.4.1). to validate that the genotype information obtained for The allele frequency for each group (survivors and the case study gorillas was in concordance with previ- non-survivors) was obtained using the -freq2 option in ously published data. We used PCAngsd [64] with the VCFtools [62]. Then, we calculated the allele frequency genotype likelihoods obtained with ANGSD (N = 6484), difference per SNP by subtracting the allele frequency in including 13 previously published whole-genomes non-survivors from the allele frequency in survivors Fontsere et al. BMC Genomics (2021) 22:735 Page 9 of 12 (ΔFrequency). We chose a threshold of ±0.2, and plotted Microbiome sequencing the allele frequencies using R. We also reported the al- 16S RNA sequencing reads were processed using QIIME lele frequency differences per SNP considering all sam- (v1) (Quantitative Insights Into Microbial Ecology) [77] to ples (pre-Ebola outbreak: survivors and non-survivors) analyze the 16S rRNA. First, paired-end raw reads were and the survivors (post-Ebola outbreak), which results in merged using fastq-join from ea-utils package [78]. Then, more modest allele frequency differences. The fixation with usearch software [79], merged FASTQ reads were fil- index (FST) between both groups was calculated using tered (−fastq_trunclen 253 and –fastq_maxee 0.5). Using VCFtools –weir-fst-pop option (Weir and Cockerham) QIIME environment, the metadata mapping file was con- [62] with a threshold at 0.15, and results were plotted in structed and validated (validate_mapping_file.py) and R. We retrieved markers with α ≤ 0.05 in the association QIIME labels were added (add_qiime_labels.py). We ap- test and a p-value <0.05 in the permutation test. The plied open-reference OTUs (Operational Taxonomic ANGSD software [32] was used to perform a replication Units) picking (pick_open_reference_otus.py). Summary of the association analysis using the following parame- statistics were computed using biom summarize-table. ters -minQ 20 -minMapQ 30 -doAsso 1 -GL 1 -out The resulting dataset was rarefied to an even depth of assocGQ_filter -doMajorMinor 1 -doMaf 1 -SNP_pval 10,000 sequences per extract (6 extracts were excluded in 1e-6 -minInd 22 -minMaf 0.02. The output of the asso- diversity analysis: G191_5782, G344_5827, G314_5824, ciation analysis are LRT values (Likelihood Ratio Test), G489_5834, G489_5835, G344_5828). Finally, we ran di- which are chi square distributed with one degree of free- versity analysis with a sequence depth of 10,000 (core_di- dom. Since we set a threshold of significance at 95% versity_analysis.py). Taxa abundance quantification and confidence, the minimum score to be significant is LRT significance or relative taxa abundance (T-test and p- = 3.84. In both association analyses, we linked the nom- values adjusted for multiple testing with Bonferroni- inally significant SNPs with their genes (Supplementary correction) were computed in R. Table S4 and S5). Genes with multiple nominally signifi- Deep sequencing of the DNA library (pre-capture) was cant SNPs were considered to be potentially more rele- performed as stated above. To remove sequencing vant. Subsequently, we compared the overlap of adapters and merge the read pairs, we used Adapter- discovered genes between the datasets (Unified Genoty- Removal v2.2.4 with default settings [80]. We then per and ANGSD) in a Venn diagram (Supplementary aligned the merged sequences to the gorilla reference Fig. S8) using the R package VennDiagram [71]. genome Kamilah_GGO_v0 using bwa mem [59] to re- move host DNA. Subsequently we filtered out potential human contaminant DNA by aligning the unmapped se- Prediction of functional consequences of the significant quences to the human reference genome hg19, resulting markers in 724,738,878 filtered sequences. MALT v0.4.1 We used VEP (v91) [72] for the functional annotation of the (MEGAN Alignment Tool) [40] was used to characterize associated SNPs. We retrieved the predicted consequence of the microbial profile, using all archaeal, viral, and bacter- each significant marker found in the potentially related genes ial reference sequences downloaded from NCBI on to Ebola immune response, as well as PolyPhen-2 [33], Sift 06.05.2019. These were indexed using malt-build to [73] and C-scores [35]. Associated loci were intersected with build a custom database. Malt-run was then used with hits in the GWAS catalogue [37] within 5000 bp. We also minimum percent identity (−−minPercentIdentity) set to performed an overrepresentation test using Panther [74] to 95, the minimum support (−−minSupport) parameter test whether any of the potentially related genes are overrep- set to 10, and the top percent value (−−topPercent) set resented in biological or functional categories compared to as 1, other parameters were set to default. The resulting the rest of targeted genes with no apparent association with rma6 files were visualized with MEGAN6 [81] and clus- Ebola. In addition, we mapped previously published RNA se- tered in a Principal Coordinate analysis (PCoA) and quencing data from six tissues (brain, cerebellum, heart, kid- Neighbor Joining Tree analysis according to microbial ney, liver, testis) in two gorilla individuals [36] to the composition on the species level (Supplementary Figs. annotated genes in the human reference genome S14 and S15). (using the Ensembl Release 75 gene models) using Abbreviations Tophat2 [75], and estimated the gene expression with EBOV: Ebola virus; EVD: Ebola virus Disease; GWAS: Genome-Wide Association htseq-count [76]. Gene expression is reported log2- Study; IUCN: International Union for Conservation of Nature; normalized, and we counted the number of reads OTUs: Operational Taxonomic Units; PCA: Principal Component Analysis; SFS: Site Frequency Spectrum; SNP: Single Nucleotide Polymorphism overlapping the candidate missense mutations to con- firm their transcriptional activity in gorillas. Values presented are the cumulative sums of RNA sequen- Supplementary Information The online version contains supplementary material available at https://doi. cing reads across individuals and tissues. org/10.1186/s12864-021-08025-y. Fontsere et al. BMC Genomics (2021) 22:735 Page 10 of 12 Author details 1 Additional file 1. Institut de Biologia Evolutiva, (CSIC-Universitat Pompeu Fabra), PRBB, Doctor Aiguader 88, 08003 Barcelona, Catalonia, Spain. 2Research and Conservation, Additional file 2. Copenhagen Zoo, 2000 Frederiksberg, Denmark. 3Department of Biology, Section for Computational and RNA Biology, University of Copenhagen, Ole Acknowledgments Maaløes Vej 5, 2200 Copenhagen, Denmark. 4Genomics of Health Research We are grateful to S. Gatti, F. Levréro, C. Genton, R. Cristescu for their Group, Hospital Universitari Son Espases (HUSE) and Institut d’Investigacions assistance with sample collection on site. We thank the teams of the Sanitaries de Balears (IDISBA), Palma, Spain. 5Center for Evolutionary ECOFAC program (EU) and African Parks Networks for logistic assistance and Hologenomics, The GLOBE Institute, University of Copenhagen, 1353 permission to work in Odzala-Kokoua National Park. The storage and extrac- Copenhagen, Denmark. 6UMR 6553, ECOBIO: Ecosystems, Biodiversity, tion of fecal samples were performed in the molecular ecology platform Evolution, CNRS/University of Rennes 1, Station Biologique de Paimpont, (UMR 6553 Ecobio, CNRS/UR1) dedicated to non-invasive samples. 35380 Paimpont, France. 7Institucio Catalana de Recerca i Estudis Avançats (ICREA, 08010 Barcelona, Catalonia, Spain. 8CRG, Centre for Genomic Regulation, Barcelona Institute of Science and Technology (BIST), 08036 Research involving animals Barcelona, Spain. 9BarcelonaBeta Brain Research Center, Pasqual Maragall Non-invasive fecal samples were collected from the field following the Foundation, 08005 Barcelona, Spain. 10University Museum, NTNU, Trondheim, standard practices. Fecal samples are exempt by the Convention on the Norway. 11Department of Evolutionary Anthropology, University of Vienna, Trade in Endangered Species of Wild Fauna and Flora (CITES). Djerassiplatz 1, 1030 Vienna, Austria. 12Population Health Sciences, Bristol Medical School, University of Bristol, Bristol, UK. 13MRC Integrative Authors’ contributions Epidemiology Unit at University of Bristol, Bristol, UK. 14CNAG-CRG, Centre for D.H., T.M.-B., P.L.G., N.M., D.V., A.N. and M.T.P.G. conceived and designed the Genomic Regulation (CRG), Barcelona Institute of Science and Technology study; D.H., J.H.R. and P.F. selected the genes for the target (BIST), Baldiri i Reixac 4, 08028 Barcelona, Spain. 15Institut Català de capture design. D.H., D.V. and C.S.O. performed experiments; C.F., M.K., D.H., Paleontologia Miquel Crusafont, Universitat Autònoma de Barcelona, Edifici J.H.R., C.H.S.O and J.N. analyzed data; P.L.G., N.M. and D.V. collected genetic ICTA-ICP, c/ Columnes s/n, 08193 Cerdanyola del Vallès, Barcelona, Spain. and demographic data; P.F., H.R.S., C.H., provided analytical support. C.F., M.K., D.H., and T.M.-B. wrote the manuscript with input and approval from the Received: 31 May 2021 Accepted: 13 September 2021 other authors. Funding C.F. is supported by “la Caixa” PhD fellowship, fellowship code LCF/BQ/DE15/ References 10360006. M.K. is supported by “la Caixa” Foundation (ID 100010434), 1. Centers for Disease Control and Prevention. History of Ebola Virus Disease. fellowship code LCF/BQ/PR19/11700002. J. N is supported by the European U.S: Department of Health & Human Services; 2018. https://www.cdc.gov/ Union’s Horizon 2020 research and innovation programme under grant vhf/ebola/history/summaries.html. Accessed 4 Feb 2020 agreement no. 676154 (ArchSci2020) and an EMBO short-term fellowship 2. World Health Organization. Ebola virus disease situation report. 2016. STF-8036. P.F. is supported by the Innovation Fund Denmark. H.R.S is sup- https://apps.who.int/iris/handle/10665/208883. Accessed 3 Jul 2018. ported by The Danish Council for Independent Research | Natural Sciences. 3. World Health Organization. External Situation Report 98 - Ebola Virus A.N. is supported by BFU2015–68649-P (MINECO/FEDER, UE). M.T.P.G. is sup- Disease. 2020. https://www.who.int/publications/i/item/10665-332654. ported by the Danish Basic Research Foundation award DNRF143. T.M.-B is 4. Kuhn JH, Becker S, Ebihara H, Geisbert TW, Johnson KM, Kawaoka Y, et al. supported by funding from the European Research Council (ERC) under the Proposal for a revised taxonomy of the family Filoviridae: Classification, European Union's Horizon 2020 research and innovation programme (grant names of taxa and viruses, and virus abbreviations. Arch Virol. 2010;155: agreement No. 864203), BFU2017–86471-P (MINECO/FEDER, UE), "Unidad de 2083–103. https://doi.org/10.1007/s00705-010-0814-x. Excelencia María de Maeztu", funded by the AEI (CEX2018-000792-M), How- 5. Centers for Disease Control and Prevention. What is Ebola Virus Disease? | ard Hughes International Early Career and Secretaria d’Universitats i Recerca Ebola (Ebola Virus Disease) | CDC. U.S: Department of Health & Human and CERCA Programme del Departament d’Economia i Coneixement de la Services; 2019. https://www.cdc.gov/vhf/ebola/about.html. Accessed 4 Feb Generalitat de Catalunya (GRC 2017 SGR 880). P.L.G., N.M. and D.V. are sup- 2020 ported by the French National agency for research via the ANR-11-JVS7–015 6. World Health Organization. Ebola virus disease. 2018. http://www.who.int/ IDiPop project. D.H. is supported by Wellcome Investigator Award (202802/Z/ en/news-room/fact-sheets/detail/ebola-virus-disease. Accessed 3 Jul 2018. 16/Z) and works in the Medical Research Council Integrative Epidemiology 7. Feldmann H, Geisbert TW. Ebola haemorrhagic fever. Lancet. 2011;377:849– Unit at the University of Bristol, which is supported by the Medical Research 62. https://doi.org/10.1016/S0140-6736(10)60667-8. Council (MC_UU_00011/1–7). This long-term research on gorillas was funded 8. Smith KF, Acevedo-Whitehouse K, Pedersen AB. The role of infectious by the ECOsystèmes FORestiers program (Ministère de l’Ecologie et du Dével- diseases in biological conservation. Anim Conserv. 2009;12:1–12. https://doi. oppement Durable France), the Espèces-Phares program (DG Environnement, org/10.1111/j.1469-1795.2008.00228.x. UE) and Lundbeck Foundation Visiting Professorship R317–2019-5 grant to 9. Walsh PD, Abernethy KA, Bermejo M, Beyers R, De Wachter P, Akou ME, T.M.-B. and M.T.P.G. This work was supported by: AEI-PGC2018–101927- et al. Catastrophic ape decline in western equatorial Africa. Nature. 2003; BI00(FEDER/UE). 422:611–4. 10. Rouquet P, Froment JM, Bermejo M, Kilbourn A, Karesh W, Reed P, et al. Availability of data and materials Wild animal mortality monitoring and human ebola outbreaks, Gabon and The sequencing dataset generated in this study is publicly available Republic of Congo, 2001-2003. Emerg Infect Dis. 2005;11:283–90. https://doi. at ENA (project ID: PRJEB43265). org/10.3201/eid1102.040533. 11. Leroy EM, Rouquet P, Formenty P, Souquière S, Kilbourne A, Froment J-M, et al. Multiple Ebola Virus Transmission Events and Rapid Decline of Central Declarations African Wildlife. Science (80-). 2004;303:387–90. https://doi.org/10.1126/ science.1092528. Ethics approval and consent to participate Not applicable. 12. Bermejo M, Rodríguez-Teijeiro JD, Illera G, Barroso A, Vilà C, Walsh PD. Ebola outbreak killed 5000 gorillas. Science (80-). 2006;314:1564. https://doi.org/1 0.1126/science.1133105. Consent for publication 13. Caillaud D, Levréro F, Cristescu R, Gatti S, Dewas M, Douadi M, et al. Gorilla Not applicable. susceptibility to Ebola virus: the cost of sociality. Curr Biol. 2006;16:489–91. 14. Maisels F, Bergl R, Williamson E. Gorilla gorilla (errata version published in Competing interests 2016). In: The IUCN Red List of Threatened Species; 2016. http://dx.doi.org/ The authors declare that they have no competing interests. 10.2305/IUCN.UK.2016-2.RLTS.T9404A17963949.en. Fontsere et al. BMC Genomics (2021) 22:735 Page 11 of 12 15. Leendertz SAJ, Wich SA, Ancrenaz M, Bergl RA, Gonder MK, Humle T, et al. 38. Gomez A, Petrzelkova K, Yeoman CJ, Vlckova K, Mrázek J, Koppova I, et al. Ebola in great apes – current knowledge, possibilities for vaccination, and Gut microbiome composition and metabolomic profiles of wild western implications for conservation and human health. Mammal Rev. 2017;47:98–111. lowland gorillas (Gorilla gorilla gorilla) reflect host ecology. Mol Ecol. 2015; 16. Genton C, Cristescu R, Gatti S, Levréro F, Bigot E, Motsch P, et al. Using 24:2551–65. demographic characteristics of populations to detect spatial fragmentation 39. Hicks AL, Lee KJ, Couto-Rodriguez M, Patel J, Sinha R, Guo C, et al. Gut following suspected ebola outbreaks in great apes. Am J Phys Anthropol. microbiomes of wild great apes fluctuate seasonally in response to diet. Nat 2017;164:3–10. https://doi.org/10.1002/ajpa.23275. Commun. 2018;9:1786. https://doi.org/10.1038/s41467-018-04204-w. 17. Genton C, Cristescu R, Gatti S, Levréro F, Bigot E, Caillaud D, et al. Recovery 40. Vågene ÅJ, Herbig A, Campana MG, Robles García NM, Warinner C, Sabin S, potential of a western Lowland Gorilla population following a major Ebola et al. Salmonella enterica genomes from victims of a major sixteenth- outbreak: Results from a ten year study. PLoS One. 2012;7:e37106. https:// century epidemic in Mexico. Nat Ecol Evol. 2018;2:520–8. doi.org/10.1371/journal.pone.0037106. 41. Gomez A, Rothman JM, Petrzelkova K, Yeoman CJ, Vlckova K, Umaña JD, 18. le Gouar PJ, Vallet D, David L, Bermejo M, Gatti S, Levréro F, et al. How et al. Temporal variation selects for diet–microbe co-metabolic traits in the Ebola impacts genetics of western lowland gorilla populations. PLoS One. gut of Gorilla spp. ISME J. 2016;10:514–26. 2009;4:e8375. https://doi.org/10.1371/journal.pone.0008375. 42. Bittar F, Keita MB, Lagier JC, Peeters M, Delaporte E, Raoult D. Gorilla gorilla 19. Honda K, Littman DR. The microbiome in infectious disease and gorilla gut: A potential reservoir of pathogenic bacteria as revealed using inflammation. Annu Rev Immunol. 2012;30:759–95. https://doi.org/10.1146/a culturomics and molecular tools. Sci Rep. 2014;4:1–5. nnurev-immunol-020711-074937. 43. Hernandez-Rodriguez J, Arandjelovic M, Lester J, de Filippo C, Weihmann A, 20. Fiebiger U, Bereswill S, Heimesaat MM. Dissecting the interplay between Meyer M, et al. The impact of endogenous content, replicates and pooling intestinal microbiota and host immunity in health and disease: Lessons on genome capture from faecal samples. Mol Ecol Resour. 2018;18:319–33. learned from germfree and gnotobiotic animal models. Eur J Microbiol https://doi.org/10.1111/1755-0998.12728. Immunol. 2016;6:253–71. https://doi.org/10.1556/1886.2016.00036. 44. Fontsere C, Alvarez-Estape M, Lester J, Arandjelovic M, Kuhlwilm M, Dieguez 21. Zheng D, Liwinski T, Elinav E. Interaction between microbiota and immunity P, et al. Maximizing the acquisition of unique reads in noninvasive capture in health and disease. Cell Res. 2020;30:492–506. https://doi.org/10.1038/ sequencing experiments. Mol Ecol Resour. 2021;21:745–61. https://doi.org/1 s41422-020-0332-7. 0.1111/1755-0998.13300. 22. Domínguez-Díaz C, García-Orozco A, Riera-Leal A, Padilla-Arellano JR, Fafutis- 45. Kerviel A, Thomas A, Chaloin L, Favard C, Muriaux D. Virus assembly and Morris M. Microbiota and its role on viral evasion: Is it with us or against us? Front plasma membrane domains: which came first? Virus Res. 2013;171:332–40. Cell Infect Microbiol. 2019;9(JUL):256. https://doi.org/10.3389/fcimb.2019.00256. https://doi.org/10.1016/J.VIRUSRES.2012.08.014. 23. Gou W, Fu Y, Yue L, Chen GD, Cai X, Shuai M, et al. Gut microbiota may 46. Lefranc MP. Immunoglobulin and T cell receptor genes: IMGT® and the birth underlie the predisposition of healthy individuals to COVID-19. medRxiv. and rise of immunoinformatics. Front Immunol. 2014;5(FEB):22. https://doi. 2020. https://doi.org/10.1101/2020.04.22.20076091. org/10.3389/fimmu.2014.00022. 24. Moeller AH, Peeters M, Ayouba A, Ngole EM, Esteban A, Hahn BH, et al. 47. Saeed MF, Kolokoltsov AA, Freiberg AN, Holbrook MR, Davey RA. Stability of the gorilla microbiome despite simian immunodeficiency virus Phosphoinositide-3 kinase-akt pathway controls cellular entry of ebola infection. Mol Ecol. 2015;24:690–7. https://doi.org/10.1111/mec.13057. virus. PLoS Pathog. 2008;4:e1000141. https://doi.org/10.1371/journal.ppa 25. Genton C, Pierre A, Cristescu R, Lévréro F, Gatti S, Pierre JS, et al. How Ebola t.1000141. impacts social dynamics in gorillas: A multistate modelling approach. J 48. Ruibal P, Oestereich L, Lüdtke A, Becker-Ziaja B, Wozniak DM, Kerber R, et al. Anim Ecol. 2015;84:166–76. https://doi.org/10.1111/1365-2656.12268. Unique human immune signature of Ebola virus disease in Guinea. Nat. 26. Reed PE, Mulangu S, Cameron KN, Ondzie AU, Joly D, Bermejo M, et al. A 2016 5337601. 2016;533:100. https://doi.org/10.1038/nature17949. New Approach for Monitoring Ebolavirus in Wild Great Apes. PLoS Negl 49. Melin AD, Janiak MC, Marrone F, Arora PS, Higham JP. Comparative ACE2 Trop Dis. 2014;8(9):e3143. https://doi.org/10.1371/journal.pntd.0003143. variation and primate COVID-19 risk. Commun Biol. 2020;3:641. 27. Gatti S, Levréro F, Ménard N, Gautier-Hion A. Population and group structure 50. Vallet D, Petit EJ, Gatti S, Levréro F, Ménard N. A new 2CTAB/PCI method of western Lowland Gorillas (Gorilla gorilla gorilla) at Lokoué. Republic of improves DNA amplification success from faeces of Mediterranean (Barbary Congo Am J Primatol. 2004;63:111–23. https://doi.org/10.1002/ajp.20045 . macaques) and tropical (lowland gorillas) primates. Conserv Genet. 2008;9: 28. Fischer A, Prüfer K, Good JM, Halbwax M, Wiebe V, André C, et al. Bonobos 677–80. https://doi.org/10.1007/s10592-007-9361-8. fall within the genomic variation of Chimpanzees. PLoS One. 2011;6:1–10. 51. Meyer M, Kircher M. Illumina Sequencing Library Preparation for Highly 29. Voight BF, Adams AM, Frisse LA, Qian Y, Hudson RR, Di Rienzo A. Multiplexed Target Capture and Sequencing. Cold Spring Harb Protoc. 2010; Interrogating multiple aspects of variation in a full resequencing data set to 2010:pdb.prot5448. https://doi.org/10.1101/pdb.prot5448. infer human population size changes. Proc Natl Acad Sci U S A. 2005;102: 52. Kircher M, Sawyer S, Meyer M. Double indexing overcomes inaccuracies in 18508–13. https://doi.org/10.1073/pnas.0507325102. multiplex sequencing on the Illumina platform. Nucleic Acids Res. 2012;40:1–8. 30. Prado-Martinez J, Sudmant PH, Kidd JM, Li H, Kelley JL, Lorente-Galdos B, 53. Rohland N, Reich D. Cost-effective, high-throughput DNA sequencing et al. Great ape genetic diversity and population history. Nature. 2013;499: libraries for multiplexed target capture. Genome Res. 2012;22:939–46. 471–5. https://doi.org/10.1038/nature12228. https://doi.org/10.1101/gr.128124.111. 31. Xue Y, Prado-Martinez J, Sudmant PH, Narasimhan V, Ayub Q, Szpak M, et al. 54. Hugerth LW, Wefer HA, Lundin S, Jakobsson HE, Lindberg M, Rodin S, et al. Mountain gorilla genomes reveal the impact of long-term population DegePrime, a program for degenerate primer design for broad-taxonomic- decline and inbreeding. Science (80-). 2015;348:242–5. range PCR in microbial ecology studies. Appl Environ Microbiol. 2014;80: 32. Korneliussen TS, Albrechtsen A, Nielsen R. ANGSD: Analysis of Next 5116–23. https://doi.org/10.1128/AEM.01403-14. Generation Sequencing Data. BMC Bioinformatics. 2014;15:256.https://doi. 55. Klindworth A, Pruesse E, Schweer T, Peplies J, Quast C, Horn M, et al. org/10.1186/s12859-014-0356-4 . Evaluation of general 16S ribosomal RNA gene PCR primers for classical and 33. Adzhubei IA, Schmidt S, Peshkin L, Ramensky VE, Gerasimova A, Bork P, next-generation sequencing-based diversity studies. Nucleic Acids Res. 2013; et al. A method and server for predicting damaging missense mutations. 41:e1. https://doi.org/10.1093/nar/gks808. Nat Methods. 2010;7:248–9. https://doi.org/10.1038/nmeth0410-248. 56. Caporaso JG, Lauber CL, Walters WA, Berg-Lyons D, Huntley J, Fierer N, et al. 34. Vaser R, Adusumalli S, Leng SN, Sikic M, Ng PC. SIFT missense predictions for Ultra-high-throughput microbial community analysis on the Illumina HiSeq and genomes. Nat Protoc. 2016;11:1–9. https://doi.org/10.1038/nprot.2015.123. MiSeq platforms. ISME J. 2012;6:1621–4. https://doi.org/10.1038/ismej.2012.8. 35. Kircher M, Witten DM, Jain P, O’roak BJ, Cooper GM, Shendure J. A general 57. Xu H, Luo X, Qian J, Pang X, Song J, Qian G, et al. FastUniq: A Fast De Novo framework for estimating the relative pathogenicity of human genetic Duplicates Removal Tool for Paired Short Reads. PLoS One. 2012;7:1–6. variants. Nat Genet. 2014;46:310–5. https://doi.org/10.1038/ng.2892. 58. Zhang J, Kobert K, Flouri T, Stamatakis A. PEAR: A fast and accurate Illumina 36. Brawand D, Soumillon M, Necsulea A, Julien P, Csárdi G, Harrigan P, et al. Paired-End reAd mergeR. Bioinformatics. 2014;30:614–20. The evolution of gene expression levels in mammalian organs. Nature. 59. Li H, Durbin R. Fast and accurate short read alignment with Burrows- 2011;478. https://doi.org/10.1038/nature10532. Wheeler transform. Bioinformatics. 2009;25:1754–60. 37. MacArthur J, Bowler E, Cerezo M, Gil L, Hall P, Hastings E, et al. The new 60. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The NHGRI-EBI Catalog of published genome-wide association studies (GWAS Sequence Alignment/Map format and SAMtools. Bioinformatics. 2009;25: Catalog). Nucleic Acids Res. 2017;45:D896–901. 2078–9. Fontsere et al. BMC Genomics (2021) 22:735 Page 12 of 12 61. Depristo MA, Banks E, Poplin R, Garimella KV, Maguire JR, Hartl C, et al. A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nat Genet. 2011;43:491–8. 62. Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, et al. The variant call format and VCFtools. Bioinformatics. 2011;27:2156–8. 63. Kuhlwilm M, Fontsere C, Han S, Alvarez-Estape M, Marques-Bonet T. HuConTest: Testing Human Contamination in Great Ape Samples. Genome Biol Evol. 2021;13. https://doi.org/10.1093/gbe/evab117. 64. Meisner J, Albrechtsen A. Inferring population structure and admixture proportions in low-depth NGS data. Genetics. 2018;210:719–31. 65. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MAR, Bender D, et al. PLINK: A tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. 2007;81:559–75. 66. Vieira FG, Lassalle F, Korneliussen TS, Fumagalli M. Improving the estimation of genetic distances from Next-Generation Sequencing data. Biol J Linn Soc. 2016;117:139–49. 67. Paradis E, Claude J, Strimmer K. APE: Analyses of phylogenetics and evolution in R language. Bioinformatics. 2004;20:289–90. 68. Korneliussen TS, Moltke I. NgsRelate: A software tool for estimating pairwise e from next-generation sequencing data. Bioinformatics. 2015;31:4009–11. 69. Hanghøj K, Moltke I, Andersen PA, Manica A, Korneliussen TS. Fast and accurate relatedness estimation from high-throughput sequencing data in the presence of inbreeding. Gigascience. 2019;8(5):giz034. https://doi.org/1 0.1093/gigascience/giz034 . 70. Nielsen R, Korneliussen T, Albrechtsen A, Li Y, Wang J. SNP calling, genotype calling, and sample allele frequency estimation from new-generation sequencing data. PLoS One. 2012;7(7):e37558. https://doi.org/10.1371/journa l.pone.0037558 . 71. Chen H, Boutros PC. VennDiagram: A package for the generation of highly- customizable Venn and Euler diagrams in R; 2011. https://doi.org/10.1186/14 71-2105-12-35. 72. McLaren W, Gil L, Hunt SE, Riat HS, Ritchie GRS, Thormann A, et al. The Ensembl Variant Effect Predictor. Genome Biol. 2016;17:17–122. https://doi. org/10.1186/s13059-016-0974-4. 73. Sim NL, Kumar P, Hu J, Henikoff S, Schneider G, Ng PC. SIFT web server: Predicting effects of amino acid substitutions on proteins. Nucleic Acids Res. 2012;40:452–7. 74. Thomas PD, Campbell MJ, Kejariwal A, Mi H, Karlak B, Daverman R, et al. PANTHER: A library of protein families and subfamilies indexed by function. Genome Res. 2003;13:2129–41. 75. Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: Accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013;14:R36. https://doi.org/10.11 86/gb-2013-14-4-r36. 76. Anders S, Pyl PT, Huber W. HTSeq--a Python framework to work with high- throughput sequencing data. Bioinformatics. 2015;31:166–9. https://doi. org/10.1093/bioinformatics/btu638. 77. Caporaso JG, Kuczynski J, Stombaugh J, Bittinger K, Bushman FD, Costello EK, et al. QIIME allows analysis of high- throughput community sequencing data. Nat Methods. 2010;7:335–6. https://doi.org/10.1038/ nmeth0510-335. 78. Aronesty E. ea-utils : “Command-line tools for processing biological sequencing data”; 2011. 79. Edgar RC. Search and clustering orders of magnitude faster than BLAST. Bioinformatics. 2017;26:2460–1. 80. Schubert M, Lindgreen S, Orlando L. AdapterRemoval v2: Rapid adapter trimming, identification, and read merging. BMC Res Notes. 2016;9:88. https://doi.org/10.1186/s13104-016-1900-2 . 81. Huson DH, Auch AF, Qi J, Schuster SC. MEGAN analysis of metagenomic data. Genome Res. 2007;17:377–86. Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.