DOI: 10.5061/dryad.vdncjsxqf
title: The expansion and diversification of pentatricopeptide repeat RNA editing factors in plants

description: The RNA-binding pentatricopeptide repeat (PPR) family comprises hundreds
 to thousands of genes in most plants, but only a few dozen in algae,
 evidence of massive gene expansions during land plant evolution. The
 nature and timing of these expansions has not been well-defined due to the
 sparse sequence data available from early-diverging land plant lineages.
 We exploit the comprehensive OneKP dataset of over 1,000 transcriptomes
 from diverse plants and algae to establish a clear picture of the
 evolution of this massive gene family, focusing on the proteins typically
 associated with RNA editing, which show the most spectacular variation in
 numbers and domain composition across the plant kingdom. We characterise
 over 2,250,000 PPR motifs in over 400,000 proteins. In lycophytes, polypod
 ferns and hornworts, nearly 10% of expressed protein-coding genes encode
 putative PPR editing factors, whereas they are absent from algae and
 complex-thalloid liverworts. We show that rather than a single expansion,
 most land plant lineages with high numbers of editing factors have
 continued to generate novel sequence diversity. We identify sequence
 variation that implies functional differences between PPR proteins in seed
 plants versus non-seed plants and which we propose to be linked to
 seed-plant-specific editing cofactors. Finally, using the sequence
 variation across the dataset, we develop a structural model of the
 catalytic DYW domain associated with C-to-U editing and identify a clade
 of unique DYW variants that are strong candidates as U-to-C RNA editing
 factors, given their phylogenetic distribution and sequence
 characteristics.
description: <i>Sources of data</i> Nucleotide sequences from the 1000 plants initiative (oneKP) (Matasci et al., 2014) were downloaded fromhttp://web.corral.tacc.utexas.edu/OneKP/. The files used for this study are the SOAPdenovo-Trans (Xie et al., 2014)assemblies (file names of the form ‘xxxx-SOAPdenovo-Trans-assembly.fa’ where xxxx is the four-letter sample code). <i>Prediction of PPR motif arrays</i> Assemblies were translated in all six frames, retaining open reading frames of at least 31 amino acids (the length of an S motif, the shortest of the PPR motif variants). <i>Hmmsearch</i> from the HMMER 3.2.1 package (Eddy, 2018) was used to search for PPR motifs, initially using the motif profiles developed by(Cheng et al., 2016), but subsequently we re-defined the DYW motif and defined the new DYW:KP variant using oneKP sequences and used these profiles in place of the DYW profile of (Cheng et al., 2016). <i>Hmmsearch</i> settings were ‘--domtblout --noali -E 0.1’. Our own in-house code (<i>PPRfinder</i>) was used to select the best-scoring non-overlapping chain of motifs from the domain table output of <i>hmmsearch</i>, with a motif score threshold of 0 for all motif types, except SS motifs (score &gt;= 10) and DYW motifs (score &gt;= 30). All motifs on the same strand, whatever frame they were in, were considered in this process. Amino acid sequences from different frames were joined where necessary to connect PPR motif arrays, with ‘X’ residues inserted to maintain the approximate length and to indicate the junction. For all the studies reported here, the protein sequences were filtered such that only sequences with either at least two adjacent PPR motifs and a sum of motif scores &gt;= 40, or a DYW motif, were retained. These are stringent thresholds that we believe eliminate almost all non-PPR sequences. <i>Alignment of PPR motifs</i> For most large alignments based on single motifs, sequences were aligned with <i>hmmalign</i> from the HMMER 3.2.1 package (Eddy, 2018). <i>Hmmalign</i> tends to give terminal gaps and unaligned tails where the sequence is a relatively poor match to the profile HMM and there is no surrounding sequence context provided. As all of the aligned motifs were originally identified using <i>hmmsearch</i> with the same profile HMM, we know that in its native sequence context each sequence does align to the HMM. Thus we wrote a script to replace terminal gaps at each end of the motif by the corresponding unaligned tails, as long as the gap and the tail were of equal length and thus there is no possible ambiguity in the placement of the amino acid residues. This procedure avoided overly sparse positions at the motif termini and a bias towards residues that best matched the HMM. For generating consensus sequences and structural models, the alignments were filtered to remove positions at which at least 50% of the sequences contained gaps. The alignments were then filtered a second time to remove sequences containing at least 5% gaps. For DYW alignments, the alignments were also filtered to remove sequences lacking the cytidine deaminase signature (HxEx<sub>25</sub>CxxC), to avoid including degenerate, non-functional sequences.  <i>Construction of phylogenetic trees</i> Large trees (thousands of sequences) were constructed using the approximate maximum likelihood approach of FastTree 2 (v. 2.1.11 double precision) (Price et al., 2010) using the JTT model of amino acid evolution (Jones et al., 1992), pseudocounts (-pseudo option) and a discrete gamma model with 20 rate categories ( -gamma option). For local support values, FastTree uses the Shimodaira-Hasegawa test on alternate nearest-neighbour interchanges for each split in the tree. <i>Modelling the DYW domain</i> Contact-derived <i>ab initio </i>modelling was performed using the RaptorX contact web server, available athttp://raptorx.uchicago.edu/ContactMap/, using a trimmed alignment of 15,177 DYW sequences as input for the contact prediction tool, and a consensus of the DYW domain as the subject for modelling (Ma et al., 2015; Wang et al., 2016; Wang et al., 2017). Models were visualised using the PyMOL Molecular Graphics System, Version 2.0 (Schrödinger, LLC), with zinc and 3-deazacytidine modelled into the best prediction based on the crystal structure of <i>E. coli </i>cytidine deaminase complexed to 3-deazacytidine, PDB ID 1ALN (Xiang et al., 1996).
description: Description of the files stored in this package<br> Figure
 numbers refer to Gutmann et al., Molecular Plant, 2019
 <br> DYW.pdb   538.32 kB    PDB format model of the DYW
 domain as depicted in Fig 6A and Fig S7A<br> DYW_contactmap.txt  
 182.63 kB    Contact map from Raptor X, depicted in Fig S7B<br>
 LL.pdb   448.41 kB    PDB format model of the LL motif as depicted in Fig
 3C<br> LL_contactmap.txt   131.78 kB    Contact map from Raptor X,
 used to build the model depicted in Fig 3C<br> PPRfinder_output.tgz
   193.58 MB    compressed output from PPRfinder for all 1KP samples; see
 https://github.com/ian-small/OneKP for details<br> all_PPR.hmm  
 295.33 kB    all PPR motif profiles built by hmmbuild and used in
 hmmsearch; the profiles are identical to those developed by Cheng et al.
 2016 except for the addition of the "KPAKA" profile (DYW:KP
 domain); see https://github.com/ian-small/OneKP for details<br>
 all_dyw_motifs.ali.fa   2.58 MB    Fasta format alignment of DYW domain
 sequences built using hmmalign and used for constructing trees and contact
 maps<br> all_dyw_motifs.nw   924.10 kB Newick format tree of DYW
 domain sequences, as depicted in Fig 2 and Fig 5<br>
 all_dyw_motifs_inc_algae.ali.fa   2.58 MB    Fasta format alignment of DYW
 domain sequences, including those from algae, built using
 hmmalign<br> all_dyw_motifs_inc_algae.nw   1 MB    Newick format
 tree of DYW domain sequences, including those from algae, as depicted in
 Fig S9B<br> dyw.hmm   64.25 kB    DYW domain profile built by
 hmmbuild and used in hmmsearch (is also included in all_PPR.hmm); see
 https://github.com/ian-small/OneKP for details<br> dyw_kp.hmm  
 63.79 kB    DYW:KP domain profile built by hmmbuild and used in hmmsearch
 (is also included in all_PPR.hmm); see https://github.com/ian-small/OneKP
 for details (the profile is named "KPAKA")<br>
 dywkp.ali.fa   291.55 kB    Fasta format alignment of DYW:KP domain
 sequences, built using hmmalign<br> dywnotkp.ali.fa   2.68 MB  
  Fasta format alignment of DYW domain sequences, not including DYW:KP
 sequences, built using hmmalign<br> ll_triplet.ali.fa   197.50 kB  
  Fasta format alignment of P1-LL-S1 triplet sequences, built using
 hmmalign<br> morf.hmm   51.23 kB    MORF domain profile built by
 hmmbuild and used in hmmsearch; used to construct Fig 4B<br>
 one_sample_per_clade_cter.nw   23.75 kB    Newick format tree of
 P2-L2-S2-E1-E2-DYW sequences, as depicted in Fig S3<br>
 pls_triplet.ali.fa   8.79 MB    Fasta format alignment of P1-L1-S1 triplet
 sequences, built using hmmalign<br> ss_triplet.ali.fa   1.10 MB  
  Fasta format alignment of SS-SS-SS triplet sequences, built using
 hmmalign