key: cord-0863320-wsryhzb8
authors: Seibert, Brittany; Joaquín Cáceres, C.; Cardenas-Garcia, Stivalis; Carnaccini, Silvia; Geiger, Ginger; Rajao, Daniela S.; Ottesen, Elizabeth; Perez, Daniel R.
title: Mild and severe SARS-CoV-2 infection induces respiratory and intestinal microbiome changes in the K18-hACE2 transgenic mouse model
date: 2021-04-21
journal: bioRxiv
DOI: 10.1101/2021.04.20.440722
sha: 7ffdc9db9a92143f9f43afa18873dd32468c78db
doc_id: 863320
cord_uid: wsryhzb8

Transmission of the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), has resulted in millions of deaths and declining economies around the world. K18-hACE2 mice develop disease resembling severe SARS-CoV-2 infection in a virus dose-dependent manner. The relationship between SARS-CoV-2 and the intestinal or respiratory microbiome is not fully understood. In this context, we characterized the cecal and lung microbiome of SARS-CoV-2 challenged K18-hACE2 transgenic mice in the presence or absence of treatment with the Mpro inhibitor GC376. Cecum microbiome showed decreased Shannon and Inv Simpson diversity index correlating with SARS-CoV-2 infection dosage and a difference of Bray-Curtis dissimilarity distances among control and infected mice. Bacterial phyla such as Firmicutes, particularly Lachnospiraceae and Oscillospiraceae, were significantly less abundant while Verrucomicrobiota, particularly the family Akkermansiaceae, were increasingly more prevalent during peak infection in mice challenged with a high virus dose. In contrast to the cecal microbiome, the lung microbiome showed similar microbial diversity among the control, low and high challenge virus groups, independent of antiviral treatment. Bacterial phyla in the lungs such as Bacteroidota decreased while Firmicutes and Proteobacteria were significantly enriched in mice challenged with a high dose of SARS-CoV-2. In summary, we identified changes in the cecal and lung microbiome of K18-hACE2 mice with severe clinical signs of SARS-CoV-2 infection. IMPORTANCE The COVID-19 pandemic has resulted in millions of deaths. The host’s respiratory and intestinal microbiome can affect directly or indirectly the immune system during viral infections. We characterized the cecal and lung microbiome in a relevant mouse model challenged with a low and high dose of SARS-CoV-2 in the presence or absence of an antiviral Mpro inhibitor, GC376. Decreased microbial diversity and taxonomic abundances of the phyla Firmicutes, particularly Lachnospiraceae, correlating with infection dosage was observed in the cecum. In addition, microbes within the family Akkermansiaceae were increasingly more prevalent during peak infection, which is observed in other viral infections. The lung microbiome showed similar microbial diversity to the control, independent of antiviral treatment. Decreased Bacteroidota and increased Firmicutes and Proteobacteria were observed in the lungs in a virus dose-dependent manner. These studies add to a better understanding of the complexities associated with the intestinal microbiome during respiratory infections.

Throughout 2020, the World Health Organization reported ~8 million confirmed COVID-55 19 cases and ~1.8 million confirmed deaths leading to a continuous increase of cases during 56 the early months of 2021 (1). The SARS-CoV-2 virus replicates and migrates to multiple tissues 57 including the airways and alveolar epithelial cells in the lungs, triggering a strong immune 58 response that may lead to exacerbation of inflammatory responses, a major complication in 59 SARS-CoV-2 patients (2-9). While many infected patients can present as asymptomatic, others 60 show clinical manifestations such as fever, shortness of breath, cough, headache, and 61 occasional gastrointestinal symptoms (10-12); however, there are still several aspects of the 62 host immune response that need to be elucidated.

The respiratory and intestinal microbiome can have direct impacts on host cells or an 64 indirect impact on the immune system during viral infections (13, 14) . Our knowledge of the SARS-CoV-2 virus and either receiving antiviral therapy with the M pro inhibitor GC-376 or vehicle 120 for 7 days post virus challenge. We performed 16S sequencing at 2-, 5-, and 14-days post 121 challenge (dpc) with a prototypic SARS-CoV-2 strain. The results of the intestinal microbiome 122 show microbial differences in alpha and beta diversity measures that are SARS-CoV-2 virus 123 dose-dependent and with little effect of GC-376 treatment on lung bacterial communities.

Clinical outcomes of K18 hACE2 transgenic mice challenged with two different doses of 126 SARS-CoV-2 virus and samples for microbiome analyses.

Taking advantage of a study evaluating antiviral activity of GC-376 against SARS-CoV-2 128 virus in the K18-hACE2 mouse model, we evaluated the microbiome composition at different 129 times after SARS-CoV-2 challenge. We and others have shown that mice challenged with 10^3 130 TCID50/mouse of the SARS-CoV-2 virus (Low/Vehicle) presented with brief reduced activity 131 and clinical signs leading to ~60% survival (43). In contrast, mice challenged with 10^5 132 TCID50/mouse of SARS-CoV-2 (High/Vehicle) presented initially with relatively normal activity 133 followed by rapid weight loss and substantial deterioration of clinical outcomes (43). By 6 dpc, 134 mice in the high virus dose group showed ~20% weight loss and all mice died or had to be 135 euthanized by 8 dpc (43). Peak virus titers for the low and high dose groups were observed at 2 136 and then 5 dpc in the nasal turbinate's and lungs (43) . Antiviral GC-376 treatment resulted in 137 milder inflammation and reduced lesions and viral loads compared to the vehicle group, 138 although it did not improve clinical outcomes (43). We analyzed the changes in the intestinal 139 and respiratory microbiome by collecting cecum and lungs from mice of the following groups: 140 PBS/Vehicle, Low/Vehicle, and High/Vehicle (Fig. 1A) . Because the respiratory tract is the 141 primary site of replication for SARS-CoV-2, we also collected lung samples from the antiviral 142 GC-376 treated group (Mock/GC-376, Low/GC-376, High/GC-376) to evaluate whether antiviral 143 intervention would affect the residential respiratory microbiome (Fig. 1A) .

We performed 16S sequencing at 2-, 5-, and 14-dpc except for lung samples in the 145 PBS/Vehicle group due to limited DNA concentrations. Of the total ceca and lung combined 146 5,098,781 raw reads obtained, while 3,103,597 reads remained after dada2 trimming, filtering, 147 merging, and chimera removal. Two lung samples, one from the Low/Vehicle at 2 dpc and one 148 from the High/Vehicle at 5 dpc were removed from the analysis because of low coverage 149 (<10,000 reads). One lung sample from the Mock/GC-376 at 14 dpc, considered an outlier 150 according to Grubbs test on taxonomic abundance (p = 6.022e-07), was also removed from the 151 analysis. Because the microbiome of the lungs can be easily contaminated, we compared the sequencing coverage of the blank extractions and negative PCR controls to samples from the 153 lung and ceca (Fig. 1B ). Blank extraction samples had an average of 5,824 reads/sample and 154 the negative PCR controls had an average of 199 reads/sample. Meanwhile, cecum samples 155 had an average of 46,827 with a minimum of 12,892 reads/sample and the lung samples had an 156 average of 42,217 with a minimum of 14,808 reads/sample (Fig. 1B) . Since the number of reads 157 for the ceca and lung are notably greater than the blanks, the difference in coverage suggests 158 that the majority of the reads in the samples are not from cross-contamination.

Microbial diversity in the cecum of SARS-CoV-2 challenge K18 hACE2 mice.

To evaluate spatial differences in microbial diversity and community structure in the 161 cecum, we evaluated alpha diversity using count data from rarified ASVs to calculate the 

In order to assess the relationship between microbial community structure and SARS-

CoV-2 challenge during the course of infection, we analyzed the number of shared ASVs. The 180 three groups shared 76 ASVs after the count data was rarified with a detection limit of 0.001 in 181 at least 90% of the samples (Fig 2D) . Using the same criteria, PBS/Vehicle had 31 unique ASVs 182 while Low/Vehicle and High/Vehicle had 6 and 12 unique ASVs, respectively ( Fig 2D) . The

Low/Vehicle group shared 21 ASVs with PBS/Vehicle and 11 ASVs with High/Vehicle, 184 suggesting a more similar microbial composition between Low/Vehicle and PBS/Vehicle (Fig.  2D ). Next, we quantified changes of the cecum microbiome composition among different SARS-

CoV-2 infected groups by comparing weighted dissimilarity distances (Bray-Curtis) within and 187 across groups (Fig. 2E) 

Since the diversity metrics suggested a difference among the low and high virus dose 203 infected mice, we investigated those differences further by analyzing the relative abundance of 

The effect of GC-376 was analyzed from homogenized lung tissue samples and we specifically 

Considering the diversity metrics suggested a limited difference in infected and control 284 mice, we further analyzed the relative abundance of the microbial communities at the phylum and family levels. The most abundant phyla within the lungs were Bacteroidota, Firmicutes,

Proteobacteria, Actinobacteriota and Verrumicrobiota. In contrast to the ceca, Bacteroidota were 287 suppressed in GC-376-treated mice exposed to low and high dose virus, with the Mock/GC-376 288 exhibiting significantly higher abundance of Bacteroidota compared to High/GC-376 (p=0.017) 289 ( Fig 6A) . The High/GC-376 group had significantly higher abundance of Firmicutes than 

We analyzed the cecum and lung microbiome changes that occur in K18-hACE2 mice 

The most distinct differences in taxonomic relative abundance within the ceca of infected 348 mice is the overall lower abundance of Firmicutes, particularly the families Lachnospiraceae and

Oscillospiraceae, increased abundance of Proteobacteria in the low virus dose group, and the 350 increase of Verrucomicrobiota, particularly the family Akkermansiaceae, in the high virus dose group at 5 dpc (Fig. 3) . In addition, a significant difference among the F/B ratio was observed 352 among the control and high dose group (Fig 3B) . Previously, F/B ratio has been associated with 353 maintaining homeostasis and changes could be indicative of dysbiosis (55). The results in this 354 report showed a decrease in the F/B ratio within the ceca, which was also observed in patients were not observed besides decreased number of ASVs within the lungs (Fig. 5A ). Similar to 411 results within this study, a previous report also showed no significant differences in Shannon 412 diversity within the lung microbiome of mice infected with influenza and the nasopharynx of 413 negative and positive SARS-CoV-2 PCR patients (18, 28). Regarding beta diversity analysis, 414 the weighted Bray-Curtis NMDS did not show separation of clusters by treatment, but rather a 415 single overlapping cluster with outliers that primarily belong to the high virus group (Fig. 5F ). In 

The most distinct differences in taxonomic relative abundance within the lung of infected 424 mice is the overall lower abundance of Bacteroidota, higher abundance of Firmicutes, and 425 higher abundance of Proteobacteria in the low and high virus dose groups compared to the 426 mock control (Fig 6) , which is consistent with previous reports (66, 67). 

Negative controls including an extraction blank and a PCR blank were included in each 512 sequencing run (2 runs total). Due to limited DNA concentrations, we were unable to sequence 513 5 of the 6 PBS/Vehicle lung samples.

Sequence processing and analysis: Primer removal and de-multiplexing was performed using 515 Illumina Basespace using default settings. Sequence analysis was performed in R (72) with 516 open-source software package 'dada2' (Version 1.16.0) (73). Each sequencing batch was processed separately until chimera removal. For each batch, the quality of the raw pair-end 518 reads was visualized and used to determine appropriate truncation of read 1 (R1) by 10 bp and 519 read 2 (R2) by 50 bp. After truncation, reads were discarded if they contained more than 2 520 maxEE "expected errors" or a quality score of less than or equal to 2. Following, each quality-521 filtered and trimmed read was processed independently by applying the trained dada2 522 algorithm. The reads were then merged with a minimum overlap of 20 bp. After merging, both 523 sequencing batches that were previously processed separately were combined, and chimeras 524 were removed using the consensus method with default settings. Taxonomy was assigned in 525 'dada2' using the native implementation of the naïve Bayesian classifier using Silva v.38 526 database. A count table and taxonomy file were created and used for downstream analysis.

Prior to diversity analysis, potential sequence contaminants were identified using 528 package 'decontam' (Version 1.8.0) (74) in RStudio (Version 1.2.5042) (75) . Briefly, potential 529 contaminants were identified by using the prevalence-based contaminant identification, which 530 relies on the principle that sequences from contaminating taxa have a higher prevalence in 531 negative control samples (extraction and PCR blanks) than true samples (74). A threshold of 0.1 532 was used to identify contaminants. In total, 14 potential contaminants were identified by 533 package 'decontam' (Table S1) ; however, all contaminants had biological relevance to the 534 sample types collected except for one, Gemmobacter, which was removed from the data set.

Following, reads that did not identify as Bacteria, contained uncharacterized Phylum, identified 536 as chloroplast and/or mitochondria were removed using the 'phyloseq' package (Version 1.32) 537 (76). Subsequently, two samples with less than 10,000 reads/sample were removed (Lungs:

Low/Vehicle at 2 dpc and High/Vehicle at 5 dpc) and one sample (Lungs: Mock/GC-376 at 14 539 dpc), considered an outlier according to Grubbs test on taxonomic abundance using the 'outlier' 540 package (p = 6.022e-07) (77), was removed.

Alpha diversity metrics including observed number of amplicon sequence variants 542 (ASVs), Shannon diversity and inverse Simpson (Inv Simpson) indexes were calculated using 543 'phyloseq'. Briefly, samples were rarified to 12,000 using command rrarefy using the 'vegan' 544 package (Version 2.57) (78). Following, the rarified counts were imported into 'phyloseq' and 545 diversity indexes were calculated using command estimaterichness. Results were graphed 546 using 'ggplot2' (Version 3.3.2) (79) and 'ggpubr' package (Version 0.4) (80). Statistical pair-wise 547 comparison employing the Wilcox rank test was performed across groups and dpc. A Venn 548 diagram of unique and shared ASVs was created using the package 'microbiome' (Version 549 1.10) (81). Rarified count data was converted to relative abundances, and then ASVs that were common among groups were combined. ASVs with a limited detection of 0.001 in at least 90% 551 of the samples were included. The Venn diagram was graphed using package 'eulerr' (Version 552 6.1.0) (82). Regarding beta diversity, weighted Bray-Curtis dissimilarity matrix was calculated 553 with a minimum of 20 and maximum of 100 random starts using the rarified count data in 554 'vegan'. A Non-metric Multi-dimensional Scaling (NMDS) plot was used to graph the dissimilarity 555 matrix using 'ggplot2'. Ellipses were constructed using command stat_ellipse in 'ggplot2' with a 556 multivariate t-distribution. All distances displayed in boxplots for comparison of within and 557 across group Bray-Curtis dissimilarities were extracted from the same distance matrix as the 558 one used for the NMDS and graphed using 'ggplot2'. Hierarchal cluster analysis of the Bray-

Curtis distances was created using command hclust with agglomeration method "average" 

Relative abundances at the phylum and family level were generated using 'phyloseq'.

First, taxa were agglomerated at the phylum or family level and then transformed into relative 569 abundance. Taxa that had less than 1% (ceca) or 2% (lung) abundance across all samples 570 (separated by ceca and lungs) were grouped together. The box plots and bar plots of the 571 relative abundances were generated using 'ggplot2'. The three samples that were previously 572 removed were not included in the analysis. Statistical pair-wise comparison among groups was 573 performed using Wilcoxon signed-rank test. A p-value below 0.05 was considered significant. A 574 rough estimation of Firmicutes/Bacteroidota ratio was calculated by dividing the relative 575 abundance of the reads assigned to Firmicutes by the relative abundance of the reads assigned 576 to Bacteroidota. Scripts used for analysis can be found on githubt at: previously. Colored bars below the circles represent the different dpc (pink= 2 dpc, light blue = 5 867 dpc, green = 14 dpc). All statistical tests were performed using Kruskal-Wallis or Wilcox-rank test for 868 pair-wise comparisons using * = p<0.05 ; ** = p<0.005; *** = p<0.0005; **** = p<0.00005. 

Corresponding phyla are noted by the colored bar to the left of the graphs following legend in A. (D)

Relative abundances (%) of individuals were calculated by agglomerating at the family level and 879 then transformed into relative abundances. Taxa that had less than 1% abundance was grouped 880 together. Groups are indicated by the bars at the bottom of the graph and dpc at the top of the 881 graph. All statistical tests were performed using Wilcox-rank test using * = p<0.05 ; ** = p<0.005; *** 882 = p<0.0005; **** = p<0.00005. 

showing the relationship of different groups and dpc using Bray-Curtis dissimilarity distance.

Hierarchical cluster analysis was performed using hclust with agglomeration method average.

Shaded colors and circles correspond to the different groups described previously. Colored bars 906 below the circles represent the different dpc (pink= 2 dpc, light blue = 5 dpc, green = 14 dpc). All 907 statistical tests were performed using Kruskal-Wallis or Wilcox-rank test for pair-wise comparisons 908 using * = p<0.05 ; ** = p<0.005; *** = p<0.0005; **** = p<0.00005. 

(B) Firmicutes/Bacteroidota ratio was calculated and graphed to analyze differences among groups.

(C) Relative abundances (%) of the most abundant families were compared via box plots.

Corresponding phyla are noted by the colored bar to the left of the graphs following legend in A. (D)

Relative abundances (%) of individuals were calculated by agglomerating at the family level and 919 then transformed into relative abundances. Taxa that had less than 2% abundance was grouped 920 together. Groups are indicated by the bars at the bottom of the graph and dpc at the top of the 921 graph. All statistical tests were performed using Wilcox-rank test for pair-wise comparisons using * = 922 p<0.05 ; ** = p<0.005; *** = p<0.0005; **** = p<0.00005. 

Since the sample size was limited (n=2 or 3) pair-wise comparisons were not performed. High/GC-376) separated by dpc from rarified ASV count table. All statistical tests were performed 940 using Kruskal-Wallis or Wilcox-rank test for pair-wise comparisons using * = p<0.05 ; ** = p<0.005; 941 *** = p<0.0005; **** = p<0.00005. 

Coronavirus disease (COVID-19) pandemic

Characteristics of SARS-CoV-2 and COVID-19

Clinical 593 features of severe patients infected with 2019 novel coronavirus: a systematic review 594 and meta-analysis

COVID-19: consider cytokine storm syndromes and 597 immunosuppression

Alterations of the Gut Microbiota in Patients with COVID-19 or H1N1 Influenza

Epidemiological, clinical and virological 607 characteristics of 74 cases of coronavirus-infected disease 2019 (COVID-19) with 608 gastrointestinal symptoms

Gastrointestinal and Other Clinical Characteristics of COVID-19

Direct 613 evidence of active SARS-CoV-2 replication in the intestine

The presence of SARS-CoV-2 RNA in the feces of COVID-19 617 patients

Targeting the gut-lung microbiota axis by means of a high-619 fibre diet and probiotics may have anti-inflammatory effects in COVID-19 infection

Epidemiological and clinical characteristics of 99 cases of 623 2019 novel coronavirus pneumonia in Wuhan, China: a descriptive study

Clinical Characteristics of 138 Hospitalized Patients With 627 2019 Novel Coronavirus-Infected Pneumonia in

Impact of the Respiratory Microbiome on Host 629 Responses to Respiratory Viral Infection. Vaccines (Basel)

Role of the intestinal microbiota in the 631 immunomodulation of influenza virus infection

Current challenges and best-practice protocols for 633 microbiome analysis

Dynamic Changes in the Microbiome and Mucosal Immune Microenvironment of the 636

Respiratory Bacteria Stabilize and Promote Airborne Transmission 639 of Influenza A Virus. mSystems 5

Influenza A virus infection impacts systemic microbiota dynamics and causes 642 quantitative enteric dysbiosis

Influenza Virus Affects Intestinal Microbiota and Secondary Salmonella Infection in 645 the Gut through Type I Interferons

Respiratory Disease following Viral Lung Infection Alters the Murine Gut Microbiota

Gut Dysbiosis during Influenza Contributes 653 to Pulmonary Pneumococcal Superinfection through Altered Short-Chain Fatty Acid 654 Production

Akkermansia 656 muciniphila Improves Host Defense Against Influenza Virus Infection

TLR5-660 mediated sensing of gut microbiota is necessary for antibody responses to seasonal 661 influenza vaccination

Alterations in Gut Microbiota of Patients With 665

COVID-19 During Time of Hospitalization

Gnotobiotic Rats Reveal 669

That Gut Microbiota Regulates Colonic mRNA of Ace2, the Receptor for SARS

Depicting SARS-CoV-2 faecal viral activity in association with gut 673 microbiota composition in patients with COVID-19

Metagenomic Next-Generation Sequencing of 676

Nasopharyngeal Specimens Collected from Confirmed and Suspect COVID-19 Patients

Nasopharyngeal Microbiota Profiling of SARS-CoV-2 Infected Patients

Gastrointestinal Manifestations of

SARS-CoV-2 Infection and Virus Load in Fecal Samples From a Hong Kong Cohort: 685 Systematic Review and Meta-analysis

Oral Bacteriotherapy in Patients With COVID-19: A Retrospective 689

Assessing gut microbial diversity from feces and rectal 692 mucosa

Comparison of the fecal, cecal, 697 and mucus microbiome in male and female mice after TNBS-induced colitis

Host genetic variation and its microbiome interactions within 701 the Human Microbiome Project

Temporal 703 and technical variability of human gut metagenomes

Effects of Psychological, Environmental and Physical Stressors on the 706

Unravelling the effects of the environment and host 708 genotype on the gut microbiome

Evaluation 710 of K18-hACE2 Mice as a Model of SARS-CoV-2 Infection

K18-hACE2 Mice for 714 Studies of COVID-19 Treatments and Pathogenesis Including Anosmia

Comparison of transgenic and 718 adenovirus hACE2 mouse models for SARS-CoV-2 infection

Publisher Correction: SARS-CoV-2 infection of human ACE2-transgenic mice 723 causes severe lung inflammation and impaired function

Efficacy of GC-376 against SARS-CoV-2 virus infection in the K18

Human angiotensin-converting enzyme 2 transgenic mice 738 infected with SARS-CoV-2 develop severe and fatal respiratory disease

COVID-19 treatments 741 and pathogenesis including anosmia in K18-hACE2 mice

K18-hACE2 mice develop 746 respiratory disease resembling severe COVID-19

Integrated gut 749 virome and bacteriome dynamics in COVID-19 patients

K18-hACE2 mice develop respiratory disease 753 resembling severe COVID-19

Considering the 755 Effects of Microbiome and Diet on SARS-CoV-2 Infection: Nanotechnology Roles

Dual role of commensal bacteria in viral 758 infections

Intestinal microbiota promote enteric virus replication and systemic 761 pathogenesis

Bacterial lipopolysaccharide binding 763 enhances virion stability and promotes environmental fitness of an enteric virus

Enteric viruses exploit the microbiota to promote infection

Viruses and the Microbiota

The Influence of Probiotics on the

the Treatment of Obesity and Inflammatory Bowel

The 772 Controversial Role of Human Gut Lachnospiraceae. Microorganisms 8

Influenza Infection on Mouse Gut Microbiome: An Exploratory Study of 775 the Role of Age-Related Microbiome Changes on Influenza Responses

Impaired type I interferon activity and inflammatory responses in severe COVID-782 19 patients

Respiratory Viral Infection-Induced 784

Microbiome Alterations and Secondary Bacterial Pneumonia

Akkermansia muciniphila in 786 the Human Gastrointestinal Tract: When, Where, and How? Microorganisms 6

Modulation 788 of Mucosal Immune Response, Tolerance, and Proliferation in Mice

Front Microbiol 2:166. 790 62

Allergic inflammation alters the lung microbiome and hinders 792 synergistic co-infection with H1N1 influenza virus and Streptococcus pneumoniae in 793 C57BL/6 mice

Modulation of Pulmonary 796

Microbiota by Antibiotic or Probiotic Aerosol Therapy: A Strategy to Promote 797 Immunosurveillance against Lung Metastases

Testing the Anna Karenina Principle in Human Microbiome

Stress and stability: applying the Anna 801 Karenina principle to animal microbiomes

Lung microbiome and coronavirus disease 2019 (COVID-803 19): Possible link and implications

Lung microbiota promotes tolerance to allergens in 806 neonates via PD-L1

Modulation of potential 808 respiratory pathogens by pH1N1 viral infection

Differences in airway microbiome and metabolome of single lung 812 transplant recipients

A simple method for estimating fifty percent endpoints

Ultra-high-817 throughput microbial community analysis on the Illumina HiSeq and MiSeq platforms

R: A language and environment for statistical computing, R Foundation 820 for Statistical Computing

DADA2: High-resolution sample inference from Illumina amplicon data

Simple statistical 825 identification and removal of contaminant sequences in marker-gene and metagenomics 826 data

RStudio: Integrated Development for R, PBC

phyloseq: an R package for reproducible interactive 830 analysis and graphics of microbiome census data

Sample Criteria for testing outlying observations

Eduard Szoecs, Wagner H. 2020. vegan: Community Ecology

ggplot2: Elegant Graphics for Data

ggpubr: 'ggplot2' Based Publication Ready Plots

_eulerr: Area-Proportional Euler and Venn Diagrams with Ellipses