Evolution of sequence traits of prion-like proteins linked to amyotrophic lateral sclerosis (ALS)

Data

Sets of orthologs for the following protein families were downloaded from OrthoDB database version 10.1 at the Eukaryota level (Kriventseva et al., 2019): TDP43 (Group 1425837at2759), hnRNPA2 (Group 1202220at2759) and EWS/TAF15/FUS (Group 1539664at2759). These data also contain paralogs. For some analyses of the hnRNPA2 data, only proteins from a reduced list of representative eukaryotes used in a previous investigation were included (An & Harrison, 2016; Su & Harrison, 2020).

Multiple sequence alignment and phylogenetic trees

Multiple sequence alignments were constructed using Clustal Omega (Sievers & Higgins, 2018), using default parameters. Phylogenetic trees were made using PhyML with Bayesian information criterion and aBayes branch support, and other parameters at defaults (Guindon et al., 2005). The aBayes branch support is a Bayesian-like transformation of the standard approximate likelihood ratio test implemented in PhyML, that is very fast and has high statistical power (Anisimova et al., 2011). Trees were saved in Newick format for input into the phylogenetic tree drawing tool Evolview version 2 (Subramanian et al., 2019). Annotations for prion-like regions, structured protein domains and compositionally biased regions (annotation described below) were formatted for input into Evolview using AWK scripts.

Annotation of protein domains

From InterPro version 66.0, we extracted Pfam (version 34.0) protein domain positions (El-Gebali et al., 2019; Mitchell et al., 2019). These were checked against PROSITE and SMART domain annotations (Letunic & Bork, 2018; Sigrist et al., 2010), as in a previous study (Su & Harrison, 2020), and they diverge by <1% in the total number of domains labelled. Since there were some OrthoDB sequences that do not have Pfam annotations, we ran the Hmmsearch program using domain Hidden Markov Models (HMMs) downloaded from the Pfam website (http://pfam.xfam.org/) to label further domains (El-Gebali et al., 2019; Mitchell et al., 2019), with e-value threshold 0.01. Also, some further RRM RNA-binding domain annotations were made by searching for RRM domain protein structures in the PDB (Protein Data Bank) (Bittrich et al., 2021), and extracting the relevant protein domain sequences from the ASTRALSCOP database version 2.06 using PDB identifiers and residue ranges (Fox, Brenner & Chandonia, 2014). These sequences were then compared to the proteomes using BLASTP 2.9.0 (Altschul et al., 1997) with e-value threshold 1e−04. All these annotations were then reduced for overlap by giving HMM-derived annotations precedence and otherwise sorting them in decreasing order of domain length and progressively flagging overlappers further down the lists for deletion.

Prion-like regions, intrinsic disorder, and predicted stress-granule recruitment

Prion-like regions were annotated using the PLAAC program with default parameters (Lancaster et al., 2014). Experiments in S. cerevisiae wherein prion-forming sequences are scrambled yet maintain their prion nature have shown that prion-forming domains can be compositionally defined (Ross, Baxa & Wickner, 2004). PLAAC uses a Hidden Markov Model trained on known prion-forming domains from budding yeast that were detected through several tests for prion-like behaviour and amyloid formation (Alberti et al., 2009). Prion-like regions have compositions judged similar to those that can form prions. Although prion-like composition is typically identified as based on glutamine and asparagine bias, often regions that are biased for subsidiary residues of prion-like regions obtain high prion-like composition scores, e.g., regions rich in tyrosine, glycine and serine. Prion-like regions were labelled if the PLAAC PRD prion-like composition score was >0.0 (as in a previous analysis (Su & Harrison, 2019)), and colour-coded according to PRD score in the annotated trees. Higher scores >15.0 are observed for known prion-forming domains (Su & Harrison, 2019). The program fLPS2 was used to annotate compositionally-biased (CB) regions (Harrison, 2006; Harrison, 2017; Harrison, 2021). Expected background frequencies = 0.05 were used in running fLPS2, with a threshold P-value of 10^–6. CB regions were grouped together according to whether their bias signatures had a common form {Xyzr….}, where X is the main biasing residue and yzr…, etc. are subsidiary biases that can be in any order (i.e., allowing for permutation of any subsidiary biases). The ten most common such CB regions are labelled on each phylogenetic tree. To analyze milder biases of experimental relevance, the ‘-z thorough’ setting was used in fLPS2, with a threshold P = 1 × 10^–4 (Harrison, 2021).

Total proportions of intrinsic disorder were calculated using IUPRED3 with default parameters (Dosztanyi, 2018; Ward et al., 2004), as in a previous work (Su & Harrison, 2020).

The server SGNN was used to estimate whether sequences can be recruited to stress granules (Iglesias et al., 2021). This tool was derived through analysis of data for the budding yeast S. cerevisiae, and thus may miss some of the compositional characteristics that drive stress granule recruitment in other organisms. Also, it may not work for regions of proteins that do not have prion-like characteristics, since it was trained to identify prion-like regions that are stress-granule recruited.

Deep conservation of protein domains and prion-like regions since last common ancestor of eukaryotes

RNA-binding protein domain content is deeply conserved since the last common ancestor of eukaryotes for each of the protein families examined. For the hnRNPA2 family, the same number of RRM RNA-binding domains is largely maintained since the last common ancestor (LCA) of eukaryotes (>88% of orthologs) (File S1; Fig. 2C). Across the TDP43 family, all plant species have an equivalent RRM domain where metazoans have the NTD (N-terminal domain) (File S1). The total number of domains in TDP43 (N-terminal domain and RRM RNA-binding domains) is also deeply conserved (over >88% of sequences) (Fig. 2A). This NTD N-terminal domain has been shown to potentially have either DNA- or RNA-binding ability (Chang et al., 2012; Qin et al., 2014). Over the EWSR1/TAF15/FUS tree, a single RRM domain is maintained across 95% of sequences (Fig. 2B).

The prion-like regions (as defined by PLAAC) are generally conserved since the LCA of metazoans for all three protein families studied (File S1), with the percentage of metazoan orthologs maintaining some prion-like composition 74% for TDP43, 79% for hnRNPA2 and 90% for EWSR1/TAF15/FUS. The most notable exceptions are for TDP43 in some specific metazoan clades (File S1).

Distinct protein forms in metazoans and plants

For the protein families that have substantial involvement in both plants and metazoans (i.e., TDP43 and hnRNPA2), we can examine the distinct sequence characteristics in each of these kingdoms. These are summarized in Fig. 3. Most notably, the orthologs in plants have less intrinsic disorder and fewer prion-like regions with weaker prion-like compositions as defined using the PLAAC program (details in Fig. 2 and File S1). Interestingly however, the plant proteins are predicted to be more recruited to stress granules by the SGNN tool (details in Table 1). Zoom-ins of the plant and vertebrate parts of the TDP43 tree are shown in Fig. 4, to highlight these trends.

For TDP43 specifically, plant sequences do not at all assign to the NTD domain, but instead tend to a RRM domain assignment at the equivalent sequence position (File S1). Although they generally lack a prion-like region (with just a handful of exceptions; File S1), many of them (28%, 43/151) have a {GRY}-/{GY}-/{GS}-biased region instead at the same relative position in their sequences, and all but one of these is predicted to be recruited to stress granules by the SGNN tool (Iglesias et al., 2021). A {GMNQS}-rich prion-like region is predominant across vertebrates, with the exception of some fish clades (File S1, region of labelled by star symbols). Other animals demonstrate diverse bias patterns that are conserved in a clade-specific manner. One prominent feature of the zoom-ins of both the plant and vertebrate tree areas is a short ~30-residue {G}- or {Gx}-rich low-complexity region within the annotated prion-like regions, indicating a possible functional significance for such a domain (Fig. 4).

Similar trends are observed in the hnRNPA2 tree (File S1). A {GNY}-, {GNSY}-, or {GY}-rich prion-like region predominates in metazoans. However, few biased regions are observed in plant orthologs where they lack a prion-like region, and they maintain low levels of annotated intrinsic disorder (File S1; Fig. 3).

Origins of FUS, TAF15 and EWSR1

FUS likely became widely conserved in an ancestor of the Bilateria, i.e., animals with embryonic bilateral symmetry, although there are six orthologs in cnidarians, so the exact evolutionary timepoint of its emergence is difficult to discern, since there may have been some initial complex patterns of gene loss at such early evolutionary stages. In a similar manner, TAF15 and EWSR1 seem to have become conserved in an early ancestor of vertebrates, but there are a small number of diverse assigned orthologs outside of the vertebrate clade (18 non-vertebrates out of 250 total for TAF15, and 6 out of 280 for EWSR1), so here again the exact evolutionary timepoint of the emergence of these proteins is also difficult to discern. Generally, the small number of non-metazoan homologs in the tree are much shorter sequences (particularly in plants), and correspondingly have less intrinsic disorder (Fig. 5).

Specific evolutionary trends for EWSR1, TAF15 and FUS

Apart from the deep conservation of RNA-binding domains and prion-like regions (Figs. 2 and 5), the individual sequence evolution of EWSR1, TAF15 and FUS has unfolded idiosyncratically.

In particular, each of these protein families have characteristic amounts of deeply conserved intrinsic disorder (Fig. 5C). The EWSR1 family has the highest proportions of annotated intrinsic disorder, typically >0.8, that have been conserved across clades since the last common ancestral sequence (Fig. 5C, File S1). Comparatively, TAF15 has intermediate levels of intrinsic disorder (mode 0.4–0.6), with the lowest for FUS proteins (mode 0.2–0.4) (Fig. 5C). These tendencies do not correlate with PLAAC prion-like composition scores, which are notably higher for FUS compared to the other two proteins (Fig. 5D). All three families maintain some degree of prion-like compositions (score >0) over the vast majority of the family members (97% of members for EWSR1, 93% for TAF15 and 85% for FUS).

A single Zn-finger domain is a deeply conserved component of both the FUS and TAF15 protein families (File S1; Table 2). However, EWSR1 has greater variation in the number of Zinc-finger domains compared to the other two protein families (Table 2), with one in five orthologs having two domains, and one in 10 having none.

The default fLPS program was used to detect longer compositionally biased regions (CBRs). We discovered that CBRs are maintained across deep and diverse clades when we examined the detail of the annotated trees (File S1). TAF15 conserves a {GR}- or {GY}-rich bias across mammals (bias P-values < 1 × 10^–10) with the C-terminal R bias maintained in other vertebrates (labelled as a {CR} bias encompassing also the Zn-finger domains, P-values < 1 × 10^–6). Where a predominant {GY} bias is labelled, it stretches across the whole sequence, whereas the {GR} bias is beyond the Zn-finger domains at the C-terminus. For FUS, the fLPS2 program assigns a biased region across most of the sequence that varies between {GQSY} or {GQRSY} in a clade-specific manner (bias P-value < 1 × 10^–24). In EWSR1, a very specific {GMPRSTY} bias predominates across all vertebrates (bias P-value < 1 × 10^–52), with the exception of fish, which conserve a similar, but also very specific {GMPQRSY} bias (P-value < 1 × 10^–46) (File S1).

We also examined for milder compositional biases of experimental relevance, i.e., involved in recruitment to stress granules, such as {R} and {Y} bias (fLPS2 P-value ≤ 1 × 10^–4) (Wang et al., 2018). Both {R}- and {Y}-biased regions occur with high frequency, but also other positively-charged regions, i.e., {K}-biased regions, are quite common, particularly in TAF15 orthologs (File S3). {R}-biased regions are characteristically shorter than {Y}-biased regions (File S4), indicating that arginine residues corresponding to tyrosine residues that might interact with them are less dispersed along the sequence. The most prevalent bias across all protein types is {G}.

Low-complexity regions

We analyzed the distribution of low-complexity regions (LCRs) annotated using the program fLPS 2.0 and parameters recommended for shorter LCR discovery (Harrison, 2021). LCRs associated with RNA-binding, such as the {GR}-/{GPR}-/{G}-rich regions seen in EWSR, TAF15 and FUS, are sometimes termed RGG or RG regions, and associated functionally with RNA-binding (Chong, Vernon & Forman-Kay, 2018). In addition to these regions, we discovered that other LCRs are also deeply conserved across kingdoms for each of three protein families examined (File S5). For example, a short {P}-rich tract (as exemplified in EWSR1, Fig. 1) is deeply conserved across EWSR1 orthologs (219/249 cases) and some other non-metazoan orthologs (20/249 cases) (14.3 ± 5.4 residues in length, median P-value = 8.6e−08). A short {A}-rich LCR is deeply conserved across many clades for the TDP43 family, specifically within the putative prion-like domain in vertebrates (11.7 ± 8.9 residues in length, median P-value = 1.1e−06, visible in Fig. 4), and corresponds to an alpha-helical region that likely functions in phase separation (Conicella et al., 2016). In general, the top ten LCRs are conserved across kingdoms with few exceptions (non-grey bars in the charts, File S5).