Natural variation in the consequences of gene overexpression and its implications for evolutionary trajectories
Abstract
Copy number variation through gene or chromosome amplification provides a route for rapid phenotypic variation and supports the long-term evolution of gene functions. Although the evolutionary importance of copy-number variation is known, little is understood about how genetic background influences its tolerance. Here, we measured fitness costs of over 4000 overexpressed genes in 15 Saccharomyces cerevisiae strains representing different lineages, to explore natural variation in tolerating gene overexpression (OE). Strain-specific effects dominated the fitness costs of gene OE. We report global differences in the consequences of gene OE, independent of the amplified gene, as well as gene-specific effects that were dependent on the genetic background. Natural variation in the response to gene OE could be explained by several models, including strain-specific physiological differences, resource limitations, and regulatory sensitivities. This work provides new insight on how genetic background influences tolerance to gene amplification and the evolutionary trajectories accessible to different backgrounds.
Introduction
Genetic variation that underlies phenotypic differences provides the material on which evolutionary selection acts. This variation includes single-nucleotide polymorphisms (SNPs), insertions and small deletions, and other structural rearrangements. DNA copy number variants (CNVs) can also serve as a powerful source of variation. CNVs span small tandem duplication of one or few genes to large segmental duplication and even chromosomal aneuploidy that amplifies many genes together (Hastings et al., 2009; Levasseur and Pontarotti, 2011). Although most duplications are likely lost shortly after creation (Ohno, 1970; Lynch and Conery, 2000; Voordeckers and Verstrepen, 2015), the functional redundancy afforded by gene duplication can support the long-term evolution of new functions (neofunctionalization) or a division of functional labor among the duplicated genes (subfunctionalization) (Graur and Wh, 2000). Neutral or nearly neutral variants can accumulate over time in a population, and this standing variation can accelerate evolution when organisms encounter a new environment (Hermisson and Pennings, 2005; Przeworski et al., 2005; Barrett and Schluter, 2008; Zheng et al., 2020). But CNVs can also produce immediate changes in cellular fitness due to the effective increase in gene expression, at least for some genes (Kondrashov, 2012). For example, yeast cultures challenged by low glucose, sulfate, or nitrogen levels benefit from amplifying genes encoding transporters of glucose (HXT6/7), sulfate (SUL1), and amino acids (GAP1) (Brown et al., 1998; Gresham et al., 2008; Gresham et al., 2010; Sanchez et al., 2017). A genome-wide study in Escherichia coli showed that 115 amplified genes, including efflux pumps/transporters, regulatory genes, and prophage genes, increased tolerance to numerous antibiotics and toxins when overexpressed (Soo et al., 2011). Consistently, amplification of transporters and resistance genes is an early event in the bacterial evolution of antibiotic resistance (Sandegren and Andersson, 2009). Furthermore, whole-chromosome duplication in human fungal pathogens can produce immediate resistance to anti-fungal agents, due to overexpression (OE) of drug efflux pumps and their regulators (Selmecki et al., 2006; Sionov et al., 2010; Ni et al., 2013; Berman and Krysan, 2020), and gene amplifications often underlies chemoresistance in cancer cells (Yasui et al., 2004; Mishra and Whetstine, 2016). These examples illustrate the importance of gene duplication in rapid phenotypic change, especially in response to drugs and environmental stresses where tolerant individuals can rapidly emerge.
While the potential evolutionary benefits afforded by gene amplification are well known, there is also a significant fitness cost (Adler et al., 2014; Moriya, 2015). The costs and consequences of gene OE are perhaps best studied in Saccharomyces cerevisiae. Protein overproduction can cause a shortage of resources, including nucleotides required for additional DNA synthesis, nucleosides consumed by transcriptional burden, and amino acids and ATP required for translation (Wagner, 2007). Which resources become limiting can depend on the environment: for example, transcription was shown to be limiting in yeast grown during phosphate starvation, whereas translation is likely limiting when cells are grown in minimal media with low amino acid availability (Kafri et al., 2016; Metzl-Raz et al., 2017; Metzl-Raz et al., 2020). Other cellular processes can become taxed as well. Several studies used the yeast gene-deletion library to investigate genes and processes required to accommodate OE of specific proteins such as GFP (Farkas et al., 2018; Kintaka et al., 2020). Although results varied somewhat by which protein was overproduced, collectively these studies showed that gene OE can put a burden on mRNA and protein export systems, protein folding chaperones, and protein degradation machinery.
Cells are also impacted by OE of specific genes and functional classes. Several studies have quantified the fitness consequences of gene-OE libraries in laboratory strains of budding yeast S. cerevisiae to reveal fitness consequences across many OE genes (Sopko et al., 2006; Ho et al., 2009; Magtanong et al., 2011; Makanae et al., 2013). Deleterious OE genes are enriched for those encoding proteins in multi-subunit complexes such as the ribosome. One model to explain this enrichment is that perturbing stoichiometric balances will perturb complex assembly and function (Papp et al., 2003; Veitia et al., 2008; Birchler and Veitia, 2012; Moriya, 2015). Cells have mechanisms to control the dosage of some of these proteins, yet high-copy OE can still overwhelm these mechanisms in the cell (Li et al., 1995; Li et al., 1996; Fewell and Woolford, 1999; Hose et al., 2015; Ascencio et al., 2021). Protein OE can also force promiscuous interactions, perturbing protein interaction networks. Proteins with intrinsically disordered regions (IDRs) are particularly susceptible to promiscuous interactions with other proteins, and deleterious OE genes are enriched for those with IDRs (Gsponer et al., 2008; Vavouri et al., 2009; Ma et al., 2010; Chakrabortee et al., 2016). In fact, proteins encoded by gene duplicates fixed in several species are under-enriched for those with IDRs (Banerjee et al., 2017), suggesting the impact on long-term evolution. Finally, OE of transcription factors, kinases, and other regulators can trigger broad downstream effects that in turn amplify the expression of other downstream proteins, further taxing proteostasis but potentially also producing phenotypes that could be beneficial (Sharifpoor et al., 2012; Moriya, 2015; Youn et al., 2017). Ultimately, the fitness benefit of gene OE must outweigh the fitness costs in a given environment in order for the duplication to be beneficial to the cell.
Although the evolutionary importance of gene duplication has long been appreciated, little is known about natural variation in the tolerance of duplication of specific genes. Variation in the cost of a gene’s duplication could have a significant influence on evolutionary trajectories that are accessible to different individuals. Anecdotal evidence shows that different individuals can vary widely in their response to OE of specific genes. For example, prior results from our lab showed that S. cerevisiae strains from different genetic lineages have unique fitness responses to overexpressed genes when strains are grown in the presence of toxins (Sardi et al., 2016). A major unanswered question is the degree to which reported trends in the fitness consequences of gene OE vary across natural isolates beyond lab strains and how natural variation influences the response to gene OE.
To investigate these questions, we expressed the same high-copy gene OE library applied previously to laboratory S. cerevisiae strains, in 15 different yeast isolates together representing four lineages and several admixed strains, to explore the variation in tolerance to gene OE. Our results distinguish universal effects common to many studied strains versus strain-specific effects, including global responses independent of the OE gene as well as gene-specific sensitivities. We present evidence for several general models explaining strain-specific variation in the response to gene OE. These results raise important implications for the accessibility of evolutionary trajectories afforded by gene OE depending on genetic background.
Results
Overview
We chose 15 genetically diverse S. cerevisiae strains for analysis, including strains from four defined genetic lineages (including European Wine, North American Oak, Asian, and West African), one commonly used lab strain (BY4743), and three strains that represent recently admixed ‘mosaic’ strains. These isolates were collected from diverse environments including soil, vineyards, sake production, sewage, and clinical samples (Capriotti, 1955; Gerke et al., 2006; Kurtzman, 1986; McCullough et al., 1998; Sniegowski et al., 2002; Strope et al., 2015; Supplementary file 1). In addition to genetic diversity, these strains display extensive phenotype variation, for example, in nitrogen and carbon utilization (Warringer et al., 2011), nutrient requirements (Liti et al., 2009; Warringer et al., 2011), and stress tolerance (Kvitek et al., 2008; Liti et al., 2009; Will et al., 2010; Warringer et al., 2011; Strope et al., 2015; Zheng and Wang, 2015; Sardi et al., 2016; Sardi et al., 2018).
Each strain was transformed with the MoBY 2.0 library that includes ~4900 open reading frames (ORFs) with their native upstream and downstream sequences, cloned into a high-copy 2-µm replicating plasmid (Ho et al., 2009; Magtanong et al., 2011). We chose the high-copy expression system to expose gene-specific fitness differences that may be too subtle to score when genes are merely duplicated. Although the lab strain replicates the empty vector at ~11 copies per haploid genome, most other strains maintain ~2–5 copies per haploid genome (see Figure 2—figure supplement 2). All strains were readily transformed with the library and grew at expected growth rates in selective media, indicating that all strains could maintain the plasmid. An aliquot of each library transformed culture was collected before and after 10 generations of competitive growth, and relative plasmid abundance was scored by quantitative sequencing of plasmid barcodes to measure changes in plasmid abundance in the population, in biological triplicate (Figure 1A, see Materials and methods).
Barcode abundance was normalized to the total number of reads per sample, thus producing a fitness score relative to the total set of genes expressed in each strain and accounting for strain-specific differences in library expression. We measured the log2(fold change) in relative barcode abundance after 10 generations of competitive growth, which we refer to as the relative fitness score. Genes that are detrimental when overexpressed will drop in frequency in the population because of reduced cell growth or because cells suppress the abundance of toxic plasmids (Makanae et al., 2013), both of which we interpret as a relative fitness defect. In contrast, beneficial plasmids will rise in frequency in the population over time. We focus on genes with a significant fitness effect when OE at a false discovery rate (FDR)<5% (see Materials and methods, Supplementary file 4).
We first validated our results by comparing fitness effects measured in lab strain BY4743 to a previous study using a similar library of yeast genes expressed from their native promoters in a similar strain background (Makanae et al., 2013). There was highly significant overlap between the 851 genes, we identified that produce a defect upon OE and previous results of Makanae et al. (p=8 ×10−45, Hypergeometric test) despite differences in media and experimental conditions (Makanae et al., 2013). Deleterious OE genes identified in both studies were enriched for genes involved in translation, including ribosomal proteins and essential genes, and a largely overlapping group of genes repressed as part of the yeast Environmental Stress Response (Gasch et al., 2000) (p<1.4×10−10, Hypergeometric test). Thus, our approach is robust to replication and comparable to previous studies, validating our methods.
We next quantified the library effects across strains (Figure 1B). We identified 4064 genes whose OE produced a reproducible fitness effect in at least one strain (FDR<0.05), with a median of 1726 genes per strain. However, there was a wide range in the number of consequential genes (Figure 2A). Mosaic strain Y2209 was affected by 635 OE genes, whereas Y12 (isolated from African palm wine but genotypically similar to Asian strains) was affected by 3060 OE genes. (We note that the low number of genes identified in YPS606 may be influenced by reduced statistical power since only duplicates of that strain were analyzed, and the lower abundance of the 2-µm plasmid in the case of YJM1592 and YJM978.) Most significant OE genes were detrimental, although there were some differences across strains. For example, whereas roughly half of the significant OE genes in the lab strain BY4743 caused a defect, over 95% of significant OE genes in the Y12 strain were detrimental (Figure 2—figure supplement 1).
In addition to the variable number of OE genes with fitness consequences, strains also varied in the severity of the defects (Figure 2B). While the median fitness cost of deleterious OE genes was not correlated overall with the number of deleterious genes per strain, strains with the most deleterious genes (NCYC3290, YJM1389, and Y12) did show an expanded range of fitness costs, with more genes showing very strong deleterious effects compared to other strains (Figure 2B). Importantly, there was no correlation between the number of deleterious OE genes and Moby 2.0 copy number maintained in the strains (Figure 2—figure supplement 2A), as expected since our normalization procedure reflects gene fitness effects relative to the overall library in that strain. Most significant genes produced fitness effects in only a subset of strains (Figure 2C), even though results were highly reproducible within strain replicates, with over half the 4064 significant genes producing a defect in only four or fewer strains. Together, these results highlight that different strains respond differently to gene OE, on both broad and gene-specific scales.
Genes whose overexpression is deleterious across many strains are functionally related
Before investigating strain-specific effects, we first characterized genes producing a defect in many strains (Figure 3). There were 431 OE genes that produced a significant defect in at least 66% of strains (FDR<0.05), and we refer to these as ‘commonly deleterious’ OE genes. This set was heavily enriched for genes involved in translation including ribosomal proteins, ribosome biogenesis factors, and other genes repressed during stress in the Environmental Stress Response. The group was also enriched for genes encoding helicases and ATP binding proteins, mitosis regulators, proteins that localize to the nucleus, and essential proteins (p<1E−4, Hypergeometric test). All of these categories remained significant if genes involved in translation were removed from the analysis (see Materials and methods), demonstrating that the enrichments are not due to overlap with translation factors.
The Balance Hypothesis (Birchler and Veitia, 2012) posits that genes encoding proteins in multi-subunit complexes or with many protein-protein interactions cause stoichiometric imbalances and thus toxicity when overexpressed. Indeed, we found that commonly deleterious genes had more protein interactions (Oughtred et al., 2021) (p=4.0×10−69, Wilcoxon test) and included more proteins that form complexes as defined by Pu et al., 2009 (p=6.6×10−12, Hypergeometric test) compared to all other genes measured in at least one strain. This confirms the result of Makanae et al. who also found that dosage-sensitive genes in the lab strain are enriched for proteins with many interactions and proteins in complexes (Makanae et al., 2013). Translation factors accounted for 163/431 (~38%) of the commonly deleterious genes and are known to participate in many interactions, raising the possibility that this functional group is driving the result. Even after removing genes involved in translation from consideration (see Supplementary file 2), the remaining commonly deleterious genes were still enriched for complex members and proteins with many physical interactions, strongly suggesting these features as driving factors in toxicity. We also noted that common deleterious OE genes contain a higher proportion of disordered proteins (as estimated by median IUPRED scores for each protein; Mészáros et al., 2018) than other genes measured in the experiment (p<4.7×10−12, Wilcoxon test), which is consistent with other analyses of deleterious OE genes (Vavouri et al., 2009; Ma et al., 2010).
Another hypothesis for dosage sensitivity is that already abundant proteins emerging from genes and transcripts that are highly transcribed and/or translated may be more subject to aggregation if further overexpressed. While the total set of common genes is expressed at higher mRNA and protein abundance, this trend was driven by translation factors and was not significant when translation factors were removed from the analysis. Thus, while high abundance of translation factors could contribute to their dosage effects, high expression alone is not a predictor of OE toxicity for other genes. Together, our work suggests that the propensity to disrupt protein interaction networks is likely a driving factor in OE gene toxicity.
Strain-specific responses to library overexpression may reflect global differences in resource allocation
Although we identified a common set of deleterious OE genes, strains clearly varied in their response to the library, in multiple distinct ways. Even for commonly deleterious genes, isolates varied in the severity of their responsiveness (Figure 3). In general, strains with a larger number of deleterious OE genes displayed more severe median relative fitness costs in response to common-gene OE than strains with fewer deleterious OE genes—the exception was North American oak-soil strains that showed aberrantly high fitness costs of common-gene OE even though they were not the most sensitive in terms of number of deleterious-gene effects (Figure 3B, purple boxes). These differences cannot be explained by gross differences in gene-OE levels, since strains with comparable Moby 2.0 copy numbers showed vastly different responses. For example, BC187 and YPS128 carry comparable plasmid copy numbers (Figure 2—figure supplement 2A), yet YPS128 is much more sensitive to commonly deleterious genes than BC187. Instead, these data suggest that some strains are more sensitive to gene amplification, even for OE genes that are commonly deleterious across many strains.
Several possibilities could explain these results. One is that cells have different capacities for tolerating protein overproduction, regardless of the OE gene. To test this, we measured growth rates in response to OE of a non-native yeast protein, GFP, expressed from the highly active TEF1 promoter on a low-copy CEN vector. Strains with the most deleterious OE genes in fact did not show higher sensitivity to GFP overexpression; however, three of the four North American oak-soil strains did, as indicated by significantly slower doubling times (Figure 4A). The reduced growth was not due to excessive GFP production as the oak strains express GFP at levels comparable to other strains tested (Figure 4—figure supplement 1A). One possibility is that these strains have reduced capacity to tolerate high protein production due to general amino acid shortage. To test, we measured growth rates during GFP OE when strains were grown in synthetic media with and without amino acids; however, growth rates of the oak strains as a group were not different than other strains, all of which grew comparably slower in the absence of amino acids (Figure 4B). The sensitivity of all strains to amino acid shortage is consistent with previous reports in the lab strain (Farkas et al., 2018). However, the degree of sensitivity to amino acid shortage did not correlate with the overall number of deleterious genes per strain, indicating that this is unlikely a driving factor explaining strain-specific effects.
Another possibility is that strains vary in the burden of protein overproduction in the context of the 2-µm plasmid, which may create a different type of stress on some strains. There was no overall correlation between the number of deleterious OE genes and the abundance of the strain’s native 2-µm plasmid (p>0.14) (although some strains with very low native 2-µm abundance also had a high number of deleterious genes [Figure 2—figure supplement 2B]). There was also no correlation between deleterious gene number and Moby 2.0 copy number. Nonetheless, we wondered if some strains may be more sensitive to the burden of the Moby 2.0 plasmid.
To test this, we measured growth rates of strains with and without the Moby 2.0 empty vector. Although many strains grew slower when expressing the empty vector under selection, some strains were more significantly affected (Figure 4C). While none of the tests passed an FDR<0.05, the trend was consistent across replicates for most strains. Indeed, the number of deleterious OE genes was correlated with the percent decrease in growth rate when strains expressed the empty Moby 2.0 vector (r=0.7, p=0.005, Figure 4D). Interestingly, strains with the greatest sensitivity to the empty Moby 2.0 vector were not the same as those most sensitive to GFP overproduction, revealing separable effects. Thus, the added stress of DNA/Moby 2.0 overproduction may render some strains more sensitive to gene OE (see Discussion).
Strain-specific responses to specific genes implicate models for variable OE tolerance
In addition to gene-independent differences in tolerating overproduction, strains also varied in their response to specific OE genes. Most overexpressed genes produced a relative fitness effect in a small number of strains, suggesting the widespread influence of genetic background (Figure 2C). To further investigate, we identified genes whose OE produced a significant fitness effect in each strain and not more than two others, which we defined as ‘strain-specific’ gene lists. The lists ranged from 41 genes in the Y2209 strain to 1763 genes in the Y12 strain (Figure 5A). It is important to note that some of these effects in strains overly sensitive to the empty vector could still reflect generalized strain sensitivities. For example, 60% of the deleterious genes identified by our criteria as ‘strain-specific’ in Y12 were shared with another of the top four strains most sensitive to the empty vector (DVBPG1373, YJM1592, YPS163, and YJM1389, Figure 4D). Thus, some of the identified genes may be deleterious if OE in other strains growing in suboptimal or stressful conditions.
Next, we investigated functional or biophysical features enriched in each strain’s list, which might implicate strain-specific constraints in tolerating gene OE. We compared, separately, genes that were specifically beneficial or detrimental in a given strain to a list of genes in that strain that were well-measured but produced no effect on fitness (FDR>0.1, see Materials and methods). To interpret the results, we also measured strain-specific gene expression differences through triplicated RNA-seq transcriptomic experiments (Supplementary file 5) to explore connections between native gene expression and the fitness consequences of gene OE.
Among the many sets of functional and biophysical enrichments (see Supplementary file 6), several themes emerged that suggest models for strain-specific responses to gene OE. The first model may be particular to our experimental design, in which the S288c allele is expressed in the library. Beneficial OE genes in BC187 and Y7568 harbored an overabundance of nonsynonymous SNPs between the overexpressed S288c allele and the strains’ native allele (p<0.0009, Wilcoxon test). One possibility is that the S288c allele may complement deleterious SNPs accumulated in those strains. For two other strains, YJM1273 and Y7568, deleterious OE genes showed a higher proportion of amino acid differences between native and expressed alleles (p<2×10−4, Wilcoxon test). Here, allelic conflict could explain strain-specific sensitivity to the S288c allele, if the focal strain evolved its own polymorphisms. These strains did not show higher genetic distances overall from S288c, raising the possibility that the high rate of allelic differences may reflect genes under accelerated evolution.
A second model explaining strain-specific gene OE effects is one in which the unique physiology of a strain causes a unique response to gene OE. This model is supported by functional enrichments for strain-specific OE gene lists, especially when those functions relate to differentially expressed genes in the same strains (Supplementary file 7). There were several functional categories that were enriched for multiple strain-specific lists. For example, OE genes beneficial to DBVPG1373 or Y12 were enriched for genes involved in the mitotic cell cycle (p<0.0004, Hypergeometric test). Both strains showed differential gene expression of cell-cycle genes: G2/M genes were expressed significantly lower in both strains and, in the case of DBVPG1373, S-phase genes were expressed significantly higher, raising the possibility of underlying differences in the strains’ cell-cycle regulation or timing. Another recurring example is reflected in nuclear-encoded mitochondrial genes. Beneficial OE genes affecting BC187 and Y389 were enriched for mitochondrial functions including respiration and cytochrome complex assembly, respectively; both strains showed altered expression of genes related to mitochondria. In contrast, genes whose OE was deleterious to Y12 were enriched for those encoding proteins in the mitochondrial matrix, and this functional category was enriched among genes expressed higher in Y12. Notably, genes related to mitochondrial function were also deleterious in a number of other strains sensitive to the empty vector (Figure 4D), raising the possibility that these genes may be generally deleterious in suboptimal or stressful conditions (or conversely that strains with inherent differences in mitochondrial function, suggested by the differential expression in the native strains, are sensitive to the vector). Although the exact connections in these cases will require experimentation to elucidate, that related functional categories are associated with a strain’s unique gene-OE susceptibilities and their unique expression differences points to physiological differences that could influence strain responses (see Discussion).
A third model is a specific example of physiological differences—strain-specific resource limitation. OE genes deleterious to vineyard strain DBVPG1373 encode proteins with a higher proportion of tryptophan compared to inconsequential genes (p=8.1×10−5, Wilcoxon test). That the encoded proteins are related by their composition hinted that limited tryptophan availability in this strain could sensitize cells to proteins high in tryptophan content. Interestingly, transcriptomic profiling revealed that genes repressed in this strain are enriched for genes encoding aromatic amino acid biosynthesis proteins (p=4.1×10−6, Hypergeometric test), and other amino acids. Tryptophan is a precursor for de novo synthesis of NAD+, and other genes in this pathway (BNA1, 4, 6, and 7) as well as genes in the nicotinic acid/nicotinamide salvage pathway (NTP1, NRK1, and TRP2) were also repressed in DBVPG1373. Together, these data raised the possibility that DBVPG1373 is sensitive to conditions that deplete tryptophan from the cell. One prediction is that DBVPG1373, but not other strains, should be especially sensitive to OE of tryptophan-containing proteins when grown in tryptophan-limiting media. To test this, we compared growth rates of DBVPG1373 and control strains BC187 and YPS128 overexpressing tryptophan-enriched genes COS1 (3.9% tryptophan) or MAL32 (3.4% tryptophan) in synthetic media with and without tryptophan. DBVPG1373 expressing Moby 2.0 vectors grew especially poorly in synthetic media compared to the other strains, perhaps obscuring the deleterious effects of COS1 and MAL32 OE (Figure 5B). Nonetheless, when overexpressing the genes, DBVPG1373 was reproducibly more sensitive in the absence of tryptophan (p≤0.05, one-tailed paired t-test), whereas the other strains were not. It is possible that this strain is more sensitive to any OE genes under these conditions, if the cellular system is taxed in the absence of tryptophan. Nonetheless, these results support the model that strains can vary in resource limitation (see Discussion).
A final model is that OE of regulators can perturb networks of downstream proteins, and strains sensitive to those networks will be overly sensitive to OE of the regulators. One possible example is seen in mosaic strain Y7568, whose deleterious OE gene list is enriched for DNA binding proteins (p=1.8×10−4, Hypergeometric test), including site-specific transcription factors (Flo8, War1, Pho2, Pdr1, Pdr8, and Hcm1). Collectively, the combined set of these factors’ targets is heavily enriched for genes encoding plasma membrane-localized proteins including drug and other transporters (p=5.7×10−5). Remarkably, the list of OE genes that are deleterious in Y7568 is also enriched for genes encoding plasma membrane proteins (p=1.8×10−4). Although not enriched above chance, 9% of the deleterious OE genes in Y7568 are direct targets of one of the six factors above (Monteiro et al., 2020). These connections give support to the model that OE of regulators can be deleterious via OE of downstream toxic genes, and that the effects can be strain-specific.
Genes highly beneficial to some strains may relate to 2-µm replication
We identified a unique cluster of 21 genes whose OE was strongly beneficial in over half (60%) of the strains, but notably not the lab strain. Although not enriched for any specific functions, the group included multiple genes involved in ribosome biogenesis/function (RRP6, RRP7, LOC1, RPL35B, and DOM34). Interestingly, over half the genes are located next to a centromere, and on closer inspection the plasmids would have cloned the centromere in the genes’ upstream regions. This was interesting because 2-µm segregation is closely coupled with chromosome segregation (Liu et al., 2014; Mehta et al., 2002). Past work showed that cloning centromeric sequences onto a 2-µm replicating plasmid reduces copy number to that of chromosome levels (Apostol and Greer, 1988). This raised the possibility that CEN sequences rather than the cloned genes could influence fitness effects in the cell, at least for a subset of these genes.
We selected two plasmids from the beneficial gene cluster that encode ERG26 (involved in ergosterol biosynthesis) or LOC1 (involved in mRNA localization and also ribosome biogenesis) and the adjacent CEN encompassed in their upstream regions. Oak-soil strain YPS128 carrying the ERG26 or LOC1 plasmids grew faster than the empty vector control, confirming the fitness benefit to this strain (Figure 6A). If the reduced copy number and growth benefit afforded to YPS128 is due only to the cloned CEN, then deleting the entire ORF or the start codon should retain the benefits provided by the CEN that remains on the plasmid. We generated derivatives of each plasmid in which the ORF (but not upstream sequence) was deleted or the start codon replaced with a stop codon (M*). Although the plasmid in which ERG26 was deleted showed some benefits, the other mutants did not, even though the ERG26 M* variant retained lower copy number (Figure 6). Thus, although half the plasmids in this cluster had cloned the CEN, the gene product may still be important. While future research will be required to disentangle why these plasmids provide a benefit, these results are yet another example in which strains respond differently to the same experimental environment.
Discussion
Our work shows that genetic background has a profound influence on how cells respond to gene OE. Out of the ~4000 genes whose OE impacted fitness in at least one isolate, only ~12% influenced fitness in 10 or more of the 15 strains. A hallmark of the 431 commonly deleterious OE genes is their potential to perturb protein-interaction networks, either because the proteins naturally display many connections or because their biophysical properties may force promiscuous interactions when overexpressed (Gsponer et al., 2008; Vavouri et al., 2009; Ma et al., 2010; Chakrabortee et al., 2016). But even for commonly deleterious genes, strains varied in the cost incurred by their OE. We suggest two main classes for strain-specific effects. One is general responses that may be independent of the overexpressed gene. For example, some strains became sensitized to gene OE in the context of the high-copy plasmid: those with greater sensitivity to the empty vector (independent of the vector copy number) showed proportionately more deleterious OE genes and with greater fitness costs. Whether this limitation is due to the burden of extra DNA or something related to the stress of 2-µm replication is not clear; nonetheless, the result unmasks strain-specific vulnerabilities that have a broad impact. We note that many of the genes scored as deleterious in these sensitive strains may cause fitness defects in other strains grown in suboptimal or stressful conditions. The second class of strain-specific effects pertains to gene-specific responses. Our results suggest several explanatory models, including strain-specific physiological differences, strain-specific resource limitation, and unique sensitivities to network perturbation by regulator amplification.
The implications of our study are several-fold. The first is that strains may have differential access to evolutionary trajectories if the cost of gene duplication varies across individuals. CNVs can produce immediate phenotypic gains for genes not subject to dosage control (Zhang et al., 2009; Kondrashov, 2012; Hose et al., 2015), but they can also produce standing genetic variation on which selection can later act. This standing variation is important for the long-term functional evolution (Ohno, 1970; Graur and Wh, 2000) and it can also accelerate evolution when selective pressures change (Zheng et al., 2020). If the cost of gene duplication, or simply increasing expression from a single-copy gene, is higher in some backgrounds, then those strains may be less likely to evolve through CNV mechanisms. An extreme example is whole-chromosomal aneuploidy, which is a potent mode of rapid evolution that is prevalent in some genetic backgrounds yet poorly tolerated in others (Torres et al., 2007; Gallone et al., 2018; Hose et al., 2020; Scopel et al., 2021). Differences in aneuploidy tolerance could be heavily influenced by different gene-specific sensitivities across strains. Consistent with the notion that these differences can affect evolutionary trajectories, several studies have found that different fungal strains evolve through different mechanisms when exposed to the same laboratory selections, where some genetic backgrounds leverage aneuploidy, polyploidy, and CNV while others do not (Filteau et al., 2015; Gerstein and Berman, 2020; Tung et al., 2021).
Another implication of our work is the interplay between gene OE, genetic background, and environment. Past work has shown that the cost of gene OE in a single strain can vary with nutrient limitation (Wagner, 2005; Kafri et al., 2016; Frumkin et al., 2017; Farkas et al., 2018; Kintaka et al., 2020). We propose that variation in environmental responses will further reveal variation in gene OE differences across strains, as hinted at by our studies. For example, strain Y7568 showed little sensitivity to GFP OE—unless amino acids were removed from the media in which case it grew among the worst across strains (Figure 4). The sensitivity of vineyard strain DBVPG1373 to tryptophan-containing proteins was exacerbated by tryptophan depletion, an environment that produced little added effect on other strains. Even the sensitivity to the Moby 2.0 empty vector may have unmasked strain sensitivities that are not evident if those strains are grown in other environments. Understanding gene-by-environment interactions is among the greatest challenges in genetics. Understanding how this interplay influences evolutionary potential is even more complicated but beginning to emerge through experimental studies (Filteau et al., 2015; Tung et al., 2021).
The results of our work also have broad application, from microbial engineering to human health. Many industrial processes use gene OE to improve microbial traits (Keasling, 1999; Xie and Fussenegger, 2018). Understanding (and ultimately predicting) how the response to engineering strategies will vary across host strains could accelerate engineering efforts (Steensels et al., 2014; Sardi and Gasch, 2017; Sardi and Gasch, 2018). Interpreting functional variants, from SNPs to CNVs, is also a major goal in human genetics and precision medicine. While already a colossal goal, incorporating genetic background interactions in such predictions will be fundamental. Elucidating mechanistic underpinnings in model organisms will continue to pave the way toward deeper understanding.
Materials and methods
Strains and growth conditions
Request a detailed protocolStrains used in this study are listed in Supplementary file 1. Unless otherwise indicated, strains were grown in rich YPD medium (10 g/L yeast extract, 20 g/L peptone, 20 g/L dextrose) in shake flasks at 30°C. Each strain was transformed with a pool of the molecular barcoded yeast ORF library (MoBY 2.0) containing 4871 pooled high copy number barcoded plasmids (Ho et al., 2009; Magtanong et al., 2011). At least 25,000 transformants were scraped from agar plates for fivefold replication of the library, and frozen stocks were made. All OE experiments were done in liquid YPD medium with G418 (200 mg/L) added for plasmid selection. Experiments interrogating single genes were performed via culture growths of yeast strains transformed with plasmids of interest, grown for 10 generations in YPD medium supplemented with G418 in shake flasks or test tubes at 30°C with shaking.
Moby 2.0 competitive growth
Request a detailed protocolThe competition experiments were performed as previously described (Ho et al., 2009; Magtanong et al., 2011; Piotrowski Jeff and Simpkins, 2015). Briefly, frozen glycerol stocks of library transformed cells were thawed into 100 ml of liquid YPD with G418 (200 mg/L) at a starting OD600 of 0.05. The remaining cells from the frozen stocks were pelleted by centrifugation and represented the starting pool (generation 0) for each strain. After precisely five generations, each pooled culture was diluted to an OD600 of 0.05 in fresh YPD containing G418, to maintain cells in log phase. At 10 generations, cells were harvested and cell pellets were stored at −80°C.
Library construction, sequencing, and analysis
Request a detailed protocolPlasmids were recovered from each pool using QIAprep spin miniprep kits (Qiagen, Hilden, Germany) after pretreatment with 1 μl R-Zymolyase (Zymo Research, Irvine, CA) and ~100 μl of glass beads, with vortexing for 5 min. Plasmid barcodes were amplified using primers containing Illumina multiplex adaptors as described in Magtanong et al., 2011; Piotrowski Jeff and Simpkins, 2015. Barcodes from three biological replicates pooled and split across three lanes on an Illumina HiSeq Rapid Run with single end 100 bp reads. Sequencing generated a median of 7,570,975 reads per barcode. Read data are available in the Short Read Archive under accession number GSE171586.
Moby normalization and analyses
Request a detailed protocolWe experimented with several normalization strategies, including TMM in the edgeR package (Robinson et al., 2010; McCarthy et al., 2012) and simple library size normalization, in which barcode reads were divided by the total barcode read count in the sample, multiplied by 1 million to rescale for edgeR analysis. The latter provided the most robust procedure with the fewest assumptions. To recapture genes that were clearly present in the starting pool but completely absent after 10 generation growth, we performed a data imputation: genes with at least 20 normalized read counts (>5th percentile of normalized reads) in all three replicates of the starting pool but missing reads from the end-point analysis received a pseudocount of 1 added to the barcode reads at 10 generations. Measured and imputed data were analyzed using edgeR version 3.22.1, using a linear model with generation (0 or 10) as a factor. Genes whose barcodes were significantly different after 10 generations of growth in each strain at an FDR<0.05 were taken as significant (Benjamini and Hochberg, 1995). Fitness scores were calculated by taking the ratio of normalized reads at generation 10 divided by reads at generation 0 (Supplementary file 4). Hierarchical clustering was performed on the log2(fold change) in normalized fitness scores using Cluster 3.0 (Eisen et al., 1998) and visualized using Java TreeView (Saldanha, 2004). We considered if differences in statistical power could explain differences in the number of significant genes. Figure 1B shows that the biological replicates were highly reproducible. The mean correlation among replicates per strain was generally high (0.74–0.89), aside of three triplicated strains (Y2209, YJM1592, and YJM978). The correlation was lower for these strains (0.55–0.65) even though their replicates agreed well (Figure 1B); the apparently lower correlation is almost certainly driven by noise in the nearly negligible fitness changes (i.e., log2 values close to 0). Strains with nearly identical correlation across replicates showed very different numbers of significant genes (e.g., Y12 and YJM1273). Thus, differences in statistical power cannot explain the differences in significant genes across strains.
Functional and biophysical enrichments were assessed using Wilcoxon rank-sum tests for continuous data (e.g., gene length, # of SNPs, and % amino acid content) and Hypergeometric tests for categorical terms, taking as the background data set the total number of measured genes (except for strain-specific gene lists, in which the background data set was a list of insignificant genes in that strain with FDR>0.1 and measured in at least two of the biological replicates). Because gene lists are heavily overlapping, standard FDR calculations over-correct p-values. We therefore took a stringent p-value of 5 ×10−4 as significant, but also cite FDR significance in data files. Genes involved in translation that were removed from several analyses are listed in Supplementary file 2. Functional and biophysical enrichments are available in Supplementary file 6. The background gene lists used for enrichments are available in Supplementary file 8.
Determining copy number using quantitative PCR
Request a detailed protocolWe measured Moby 2.0 plasmid copy numbers in strains grown 10 generations in log-phase as described above. Plasmid DNA was extracted from frozen cell pellets using phenol/chloroform and ethanol precipitation, which recovers both plasmid and genomic DNA. Quantitative PCR experiments were conducted using a Roche LightCycler 480 II and Roche LightCycler 480 SYBR Green I Master SYBR-Green (Bio‐Rad, Hercules, CA). Primers were designed to detect the KAN-MX resistance gene located on plasmids and genomic TUB1 (control) (Supplementary file 3). CT values for each sample were measured in technical triplicate with all experiments done in greater than three biological replicates. The CT values for KAN were internally normalized to TUB1 expressed from the genome and under an extreme constraint on copy number. KAN/TUB1 ratios measured for each isolate carrying the 2-µm plasmid were adjusted to BY4743 KAN/TUB1 ratios measured as an internal control in each experiment. Data were scaled to BY4743 values, which were adjusted relative to a KAN-MX marked CEN copy number measured in BY4743 (in the same way outlined above for Moby 2.0 plasmids).
2-µm copy number analysis
Request a detailed protocolWe determined the native copy number of the 2-µm gene, REP1, using publicly available DNA sequencing data for each strain (Bergström et al., 2014; Hose et al., 2015; Strope et al., 2015). (Note: BY4741 sequence was used instead of BY4743.) We mapped the sequencing data for each strain to a S. cerevisiae genome using BWA-MEM (version 0.7.12-r1039; Li and Durbin, 2010). Summed read counts for each gene were calculated by HT-Seq (version 0.6.0; Anders et al., 2015). Read counts were normalized using RPKM.
Cloning
Request a detailed protocolTo express GFP, 343 bp upstream of TEF1 (TEFPROM) was PCR amplified and sewn to a PCR product capturing the GFP ORF and ADH1 terminator (ADH1TERM) taken from the Yeast GFP Clone Collection (Thermo Fischer Scientific). PCR product was transformed into yeast with linearized Moby 1.0 empty vector (Ho et al., 2009) and homologous recombinants were selected and verified by sequencing.
Moby 2.0 plasmids expressing ERG26 and LOC1 were isolated from E. coli using a Qiagen Spin Miniprep Kit. ERG26 AND LOC1 deletions were generated by site-directed mutagenesis. The first methionine codon of each ORF was mutated to TAG using quick-change cloning (Wang and Malcolm, 1999). All constructs were verified with Sanger sequencing.
Western blot analysis
Request a detailed protocolYeast strains were grown in synthetic complete media to log phase (OD600 ~0.4). CEN-GFP was monitored by Western blot analysis, loading OD-normalized cells in sample buffer and using rabbit anti-GFP (Abcam) and mouse anti-PGK1 (Abcam) as a loading control, and imaging on the Licor Odyssey Infrared Imager.
Transcriptome profiling (RNA-Seq) and analysis
Request a detailed protocolYeast strains described in Supplementary file 1 were grown in biological triplicate in rich YPD medium with G418 at 30°C with shaking, for three generations to an OD600 ~0.5. Cultures were pelleted by centrifugation and flash frozen with liquid nitrogen and maintained at −80°C until RNA extraction. Total RNA was extracted by hot phenol lysis (Gasch, 2002), digested with Turbo DNase (Invitrogen) for 30 min at 37°C, and precipitated with 5 M lithium acetate for 30 min at −20°C. rRNA depletion was performed using the Ribo-Zero (Yeast) rRNA Removal Kit (Illumina, San Diego, CA) and libraries were generated according to the TruSeq Stranded Total RNA sample preparation guide (revision E). cDNA synthesis was performed using fragment prime finish mix (Illumina, San Diego, CA) and purified using Agencourt AMPure XP beads (Beckman Coulter, Indianapolis, IN). Illumina adaptors were ligated to DNA using PCR (10 cycles). The samples were pooled, resplit, and run across three lanes on an Illumina HiSeq 2500 sequencer, generating single-end 100 bp reads, with ~7,494,848 reads per sample. Data are available in GEO accession number GSE171585 and supplementary file 5.
Reads were processed using Trimmomatic version 0.3 (Bolger et al., 2014), and mapped to the S288c reference genome (version R64-1-1) with BWA-MEM (version 0.7.12-r1039; Li and Durbin, 2010). Read counts for each gene were calculated by HT-Seq (version 0.6.0; Anders et al., 2015). Differentially expressed genes were identified by edgeR (Robinson et al., 2010) using a linear model with strain background as a factor and paired replicates, identifying genes differentially expressed in each strain relative to the average of all strains using an FDR cutoff of 0.05 (Benjamini and Hochberg, 1995). Hierarchical clustering was performed by Cluster 3.0 (Eisen et al., 1998) and visualized using Java TreeView (Saldanha, 2004). There was a total of 4802 genes that were significant in at least one strain.
Data availability
Barcode sequencing data are available in the Short Read Archive under accession number GSE171586. RNA-Seq data are available in GEO accession number GSE171585.
-
NCBI Gene Expression OmnibusID GSE171586. Natural Variation in the Fitness Consequences of Gene Amplification in Wild Saccharomyces cerevisiae Isolates [Bar-seq].
-
NCBI Gene Expression OmnibusID GSE171585. Natural Variation in the Fitness Consequences of Gene Amplification in Wild Saccharomyces cerevisiae Isolates [RNA-seq].
References
-
High fitness costs and instability of gene duplications reduce rates of evolution of new genes by duplication-divergence mechanismsMolecular Biology and Evolution 31:1526–1535.https://doi.org/10.1093/molbev/msu111
-
Adaptation from standing genetic variationTrends in Ecology & Evolution 23:38–44.https://doi.org/10.1016/j.tree.2007.09.008
-
Controlling the false discovery rate: a practical and powerful approach to multiple testingJournal of the Royal Statistical Society 57:289–300.https://doi.org/10.1016/s0166-4328(01)00297-2
-
A high-definition view of functional genetic variation from natural yeast genomesMolecular Biology and Evolution 31:872–888.https://doi.org/10.1093/molbev/msu037
-
Drug resistance and tolerance in fungiNature Reviews Microbiology 18:319–331.https://doi.org/10.1038/s41579-019-0322-2
-
Trimmomatic: a flexible trimmer for illumina sequence dataBioinformatics 30:2114–2120.https://doi.org/10.1093/bioinformatics/btu170
-
Multiple duplications of yeast hexose transport genes in response to selection in a glucose-limited environmentMolecular Biology and Evolution 15:931–942.https://doi.org/10.1093/oxfordjournals.molbev.a026009
-
Yeasts in some netherlands soilsAntonie Van Leeuwenhoek 21:145–156.https://doi.org/10.1007/BF02543809
-
Gene architectures that minimize cost of gene expressionMolecular Cell 65:142–153.https://doi.org/10.1016/j.molcel.2016.11.007
-
Origins, evolution, domestication and diversity of Saccharomyces beer yeastsCurrent Opinion in Biotechnology 49:148–155.https://doi.org/10.1016/j.copbio.2017.08.005
-
Genomic expression programs in the response of yeast cells to environmental changesMolecular Biology of the Cell 11:4241–4257.https://doi.org/10.1091/mbc.11.12.4241
-
Yeast genomic expression studies using DNA microarraysMethods in Enzymology 350:393–414.https://doi.org/10.1016/s0076-6879(02)50976-9
-
Mechanisms of change in gene copy numberNature Reviews Genetics 10:551–564.https://doi.org/10.1038/nrg2593
-
Gene-expression tools for the metabolic engineering of BacteriaTrends in Biotechnology 17:452–460.https://doi.org/10.1016/S0167-7799(99)01376-1
-
Gene duplication as a mechanism of genomic adaptation to a changing environmentProceedings of the Royal Society B: Biological Sciences 279:5048–5057.https://doi.org/10.1098/rspb.2012.1108
-
The ARS culture collection: present status and new directionsEnzyme and Microbial Technology 8:328–333.https://doi.org/10.1016/0141-0229(86)90130-4
-
Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variationNucleic Acids Research 40:4288–4297.https://doi.org/10.1093/nar/gks042
-
Epidemiological investigation of vaginal Saccharomyces cerevisiae isolates by a genotypic methodJournal of Clinical Microbiology 36:557–562.https://doi.org/10.1128/JCM.36.2.557-562.1998
-
Gene transcription as a limiting factor in protein production and cell growthG3: Genes, Genomes, Genetics 10:3229–3242.https://doi.org/10.1534/g3.120.401303
-
Different facets of copy number changes: permanent, transient, and adaptiveMolecular and Cellular Biology 36:1050–1063.https://doi.org/10.1128/MCB.00652-15
-
YEASTRACT+: a portal for cross-species comparative genomics of transcription regulation in yeastsNucleic Acids Research 48:D642–D649.https://doi.org/10.1093/nar/gkz859
-
Quantitative nature of overexpression experimentsMolecular Biology of the Cell 26:3932–3939.https://doi.org/10.1091/mbc.E15-07-0512
-
BookChemical Genomic Profiling via Barcode Sequencing to Predict Compound Mode of ActionIn: Hempel Jonathan E, Williams C. H, editors. Chemical Biology: Methods and Protocols. New York: Springer. pp. 299–318.
-
Up-to-date catalogues of yeast protein complexesNucleic Acids Research 37:825–831.https://doi.org/10.1093/nar/gkn1005
-
Java treeview--extensible visualization of microarray dataBioinformatics 20:3246–3248.https://doi.org/10.1093/bioinformatics/bth349
-
Bacterial gene amplification: implications for the evolution of antibiotic resistanceNature Reviews Microbiology 7:578–588.https://doi.org/10.1038/nrmicro2174
-
Leveraging Genetic-Background effects in Saccharomyces cerevisiae to improve lignocellulosic hydrolysate toleranceApplied and Environmental Microbiology 82:5838–5849.https://doi.org/10.1128/AEM.01603-16
-
Improving industrial yeast strains: exploiting natural and artificial diversityFEMS Microbiology Reviews 38:947–995.https://doi.org/10.1111/1574-6976.12073
-
Energy constraints on the evolution of gene expressionMolecular Biology and Evolution 22:1365–1374.https://doi.org/10.1093/molbev/msi126
-
Energy costs constrain the evolution of gene expressionJournal of Experimental Zoology Part B: Molecular and Developmental Evolution 308B:322–324.https://doi.org/10.1002/jez.b.21152
-
Designing cell function: assembly of synthetic gene circuits for cell biology applicationsNature Reviews Molecular Cell Biology 19:507–525.https://doi.org/10.1038/s41580-018-0024-z
-
Copy number variation in human health, disease, and evolutionAnnual Review of Genomics and Human Genetics 10:451–481.https://doi.org/10.1146/annurev.genom.9.081307.164217
Article and author information
Author details
Funding
National Cancer Institute (R01CA229532)
- James Hose
- Adam Jochem
- Audrey P Gasch
U.S. Department of Energy (DE-SC0018409)
- Michael Place
- Audrey P Gasch
National Institutes of Health (GT32GM007133)
- DeElegant Robinson
National Human Genome Research Institute (5T32HG002760)
- DeElegant Robinson
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
The authors thank Peipei Wang and Shinhan Shiu for their input on data analyses, Auguste Dutcher for calculating IUPRED scores, and members of the Gasch Lab for useful feedback.
Copyright
© 2021, Robinson et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics
-
- 1,986
- views
-
- 191
- downloads
-
- 23
- citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading
-
- Chromosomes and Gene Expression
- Genetics and Genomics
The enhancer-promoter looping model, in which enhancers activate their target genes via physical contact, has long dominated the field of gene regulation. However, the ubiquity of this model has been questioned due to evidence of alternative mechanisms and the lack of its systematic validation, primarily owing to the absence of suitable experimental techniques. In this study, we present a new MNase-based proximity ligation method called MChIP-C, allowing for the measurement of protein-mediated chromatin interactions at single-nucleosome resolution on a genome-wide scale. By applying MChIP-C to study H3K4me3 promoter-centered interactions in K562 cells, we found that it had greatly improved resolution and sensitivity compared to restriction endonuclease-based C-methods. This allowed us to identify EP300 histone acetyltransferase and the SWI/SNF remodeling complex as potential candidates for establishing and/or maintaining enhancer-promoter interactions. Finally, leveraging data from published CRISPRi screens, we found that most functionally verified enhancers do physically interact with their cognate promoters, supporting the enhancer-promoter looping model.
-
- Genetics and Genomics
- Neuroscience
Continued methodological advances have enabled numerous statistical approaches for the analysis of summary statistics from genome-wide association studies. Genetic correlation analysis within specific regions enables a new strategy for identifying pleiotropy. Genomic regions with significant ‘local’ genetic correlations can be investigated further using state-of-the-art methodologies for statistical fine-mapping and variant colocalisation. We explored the utility of a genome-wide local genetic correlation analysis approach for identifying genetic overlaps between the candidate neuropsychiatric disorders, Alzheimer’s disease (AD), amyotrophic lateral sclerosis (ALS), frontotemporal dementia, Parkinson’s disease, and schizophrenia. The correlation analysis identified several associations between traits, the majority of which were loci in the human leukocyte antigen region. Colocalisation analysis suggested that disease-implicated variants in these loci often differ between traits and, in one locus, indicated a shared causal variant between ALS and AD. Our study identified candidate loci that might play a role in multiple neuropsychiatric diseases and suggested the role of distinct mechanisms across diseases despite shared loci. The fine-mapping and colocalisation analysis protocol designed for this study has been implemented in a flexible analysis pipeline that produces HTML reports and is available at: https://github.com/ThomasPSpargo/COLOC-reporter.