BACKGROUND
aquificales are gram-negative, non-sporulating bacteria that are thermophilic to hyperthermophilic
<cit> , living in terrestrial and marine hot springs.
they are autotrophs that primarily fix carbon by the tricarboxylic acid  cycle
<cit> .
the hyperthermophile a. aeolicus, living under extreme temperatures of up to 95°c, has been proposed to have adopted 10% of its protein-coding genes by horizontal gene transfer
<cit>  from archaea.
accumulation of all the special properties of thermophiles  are rarely understood.
special protein-protective mechanisms have been analyzed
<cit> , but we are far away from a comprehensive understanding of the molecular biology of extremophilic bacteria.
beyond idiosyncratic features of aquificales genomes, our interest focussed on their transcriptomes.
experimentally, we performed a deep sequencing analysis on the model hyperthermophile a. aeolicus with the primary goal of identifying novel ncrnas candidates.
ncrnas are known to have various functions in all domains of life.
apart from their general importance as gene expression regulators
<cit> , ncrnas are involved in processing
<cit>  and translation
<cit>  of other genes, in defending genomes from viral invasion
<cit> , in shaping and maintenance of bacterial chromosome architecture
<cit> , and they can even be multifunctional
<cit> .
according to 16s rrna analysis, the aquificales constitute the most deeply rooted bacterial group
<cit> .
however, protein-based phylogenetic reconstructions are not in line with this model
<cit> .
we compared the genomes of the three aquificales families, i.e.
aquificaceae, hydrogenothermaceae and desulfurobacteriaceae.
we have extended the protein annotation of the mentioned aquificales and reconstructed the phylogenetic position of these species based on 16s rrnas as well as on a set of orthologous proteins.
moreover, we have identified ncrnas based on known homologs and present a complete set of novel ncrna candidates based on sequence analyses and deep sequencing data obtained for a. aeolicus.
for selected ncrna loci, we provide independent experimental evidence for their expression.
methods
genomes
we analyzed the genomes of the following species split into their respective families:
– aquificaceae: aquifex aeolicus vf <dig> , hydrogenivirga sp.
128-5-r1- <dig> , hydrogenobacter thermophilus tk- <dig> , thermocrinis ruber , thermocrinis albus dsm  <dig> , hydrogenobaculum sp.
y04aas <dig> ,
– hydrogenothermaceae: sulfurihydrogenibium sp.
yo3aop <dig> , sulfurihydrogenibium azorense az-fu <dig> , persephonella marina ex-h <dig> , and
– desulfurobacteriaceae: desulfobacterium thermolithotrophum dsm  <dig> , and thermovibrio ammonificans hb- <dig> .
accession numbers and sources of genomes are listed in the electronic supplemental material
http://www.rna.uni-jena.de/supplements/aquificales/index.html.
whole-genome alignments were constructed using pomago
<cit>  and tba
<cit>  with default parameters.
pomago alignments were computed separately for each species as reference.
the tba alignment was projected to each of the reference genomes.
coverage, alignment quality  and gap ratio are given in figure
<dig>
<dig> and in the supplemental material.
tmrnas are classified into two types .
the phylogenetic tree shown at the top of the table is based on the whole genome as well as 16s rrna analysis of the  <dig> aquificales species.
it reproduces the results presented in
<cit> .
for further information, see supplemental material.
aae – a. aeolicus, hvi – hydrogenivirga sp., hth – h. thermophilus, hba – hydrogenobaculum sp., tal – t. albus, tru – t. ruber, pma – p. marina, saz – s. azorense, ssp – sulfurihydrogenibium sp., dth – d. thermolithotrophum, tam – t. ammonificans, rs – riboswitch, wsop – weighted sum-of-pairs score
<cit> , * denotes the hydrogenivirga sp.
genome of unfinished assembly.
extension of protein annotation
we used bacprot  to complement the present annotation of protein-coding genes for each aquificales genome.
it uses a database of groups of orthologous protein-coding genes present in most bacteria
<cit> .
matches in the genome of interest are annotated, and species-specific features like codon usage, shine-dalgarno sequences, pribnow box motifs and rho-independent terminators are used to predict additional protein-coding genes.
to actually achieve a de novo annotation, we excluded all aquificales genes from the reference database.
alternative start codons like att and ctg were considered as well
<cit> .
re-annotated and previously annotated proteins  and statistics  for each species are provided in the supplemental material.
annotation of ncrnas by homology
we used gorap  to annotate ncrnas in the following manner: transfer-rnas  were detected by trnascan-se
<cit>  with the option -b for bacteria.
split trnas were searched using splits
<cit> .
by applying aragorn , we searched for trnas containing introns
<cit> .
searches for rnase p rna were conducted with bcheck
<cit> .
for the detection of putative crispr loci, crt
<cit>  and crisprfinder <cit>  were used.
we searched for cas protein genes by blast
<cit>  based on known cas genes
<cit> ).
to find further ncrnas, we used blast and infernal
<cit> .
seed sequences from the rfam  database
<cit>  and european ribosomal rna database <cit>  were used as query with an e-value ≤ <dig>  for blast and the rfam-provided family specific noise cutoffa for infernal.
ncrnas expected to escape from detection  were searched in a second step with rnabob <cit>  for short motif search in combination with rnasubopt, rnaduplex, rnacofold, rnaalifold and rnaup from the rna vienna package
<cit> .
for verification, we aligned candidates with clustalw
<cit>  or locarnate
<cit> .
stockholm alignments were adjusted by hand in the emacs ralee mode <cit> .
resulting stockholm alignments are supplied in the supplemental material in the general feature format  as well as in fasta  and stockholm  formats.
phylogenetic reconstruction
protein-based phylogeny was performed based on the official ncbi <cit>  annotations for  <dig> bacteria shown in the supplemental material.
in addition to eleven aquificales species, we included two archaea as outgroup and a wide phylogenetic range of  <dig> bacterial species representing all bacterial clades.
protein sequences were clustered using proteinortho <cit>  in the blastp+-mode, thus performing a pairwise all-against-all comparison of sequences from different species to derive orthologous relationships.
whenever an orthologous group did not have a member in a certain species, we applied tblastn to the respective genome to complement for potentially incomplete annotations.
the highest scoring alignment to an orf above a fairly high e-value ≤ 10- <dig> was added to the initial protein annotation.
finally, proteinortho was applied again using the expanded annotation.
for a high resolution phylogeny within the aquificales, we created a whole genome alignment using pomago.
the alignment was analyzed using raxml
<cit>  with a gamma model of rate heterogeneity with an estimate on the proportion of invariable sites and  <dig> rapid bootstraps.
in an additional phylogenetic analysis we used single-copy orthologous proteins present in at least 50% of all species in the set .
each protein group was aligned separately using dialign-tx <cit> .
both ends of the group’s alignments were cropped to remove leading and tailing gaps.
the remaining sequences were concatenated resulting in a  <dig>  aa long alignment and applied to raxml using the lg substitution model
<cit>  as well as the gamma model of rate heterogeneity with  <dig> rapid bootstraps.
the 16s rrna-based phylogeny was computed with mafft
<cit>  using the l-ins-i method with  <dig> iterations.
we used different approaches:  neighbor joining with the kimura correction model
<cit>  ,  bayesian inference with mrbayes
<cit>  with default parameters,  maximum likelihood with raxml
<cit>   with the base substitution models  gtrgamma  and  gtrcat for the bootstrapping phase.
for all previously mentioned methods the archaea methanobacterium sp.
al- <dig> and archaeoglobus fulgidus were used as outgroup.
as state of the art, we have estimated a tree with  sate
<cit>  .
related sequences were aligned with mafft and subsequently merged by muscle
<cit> .
the tree was computed using raxml.
drna-seq of a. aeolicus total cellular rna
transcriptome analysis of a. aeolicus was based on cdna libraries from a differential deep sequencing approach
<cit> .
a. aeolicus cells, provided by m. thomm and r. huber , were grown for  <dig> day  and harvested as described
<cit> .
for preparation of total cellular rna, we used the hot phenol method
<cit> : cell pellets were resuspended in extraction buffer  and incubated for  <dig> min at room temperature with  <dig>  volumes of lysozyme .
sds was added to a final concentration of 1% followed by vigorous vortexing.
after addition of  <dig> volume phenol  and vortexing, the mixture was incubated for  <dig> min at 65°c, then cooled on ice for  <dig> min, and centrifuged for  <dig> min at 4°c and  <dig> g. phenol extraction was repeated, followed by chloroform  extraction and ethanol precipitation.
finally, the dna was digested with  <dig> u turbo dnase  for  <dig> min at 37°c, followed by addition of another  <dig> u dnase and incubation for another  <dig> min at 37°c.
subsequently, the rna was subjected to phenol/chloroform extraction and ethanol precipitation.
after redissolving the rna in double-distilled water, its concentration was determined by uv spectroscopy.
before cdna library construction, the rna was split into two fractions; one fraction was treated with terminator 5’ p-dependent exonuclease  for depletion of transcripts carrying a 5’-monophosphate.
both fractions were treated with tobacco acid phosphatase  before 5’-linker ligation, poly tailing and conversion into cdna .
the cdna libraries were then sequenced on a roche flx sequencer and resulted in the -library with  <dig>  reads and the -library  containing the enriched primary transcripts.
detection of novel ncrnas
we used the igb
<cit>  to visualize the following features of a. aeolicus:  nucleotide sequence;  local gc-content ;  protein genes annotated by ncbi <cit>  and bacprot;  locally stable secondary structures: calculation was performed with rnalfold with options -d <dig> and -l <dig> for both strands with a maximum base-pair span of  <dig> nucleotides.
sequences with local structures of fewer than  <dig> nt were discarded.
for the prediction of thermodynamically stable rna structures, each sequence was shuffled  <dig> times while preserving the dinucleotide frequencies; to classify extraordinarily stable rna secondary structures, we chose to use a z-score cutoff of - <dig>  ;  conserved regions among the aquificales: with default parameters of tba and pomago we aligned  <dig> genomes; the tba alignment was projected to each of the reference genomes; coverage, wsop and gap ratio are given in figure
1;  novel ncrnas: novel ncrna candidates were predicted using rnaz.
we used rnazwindow.pl –min-seqs= <dig> and rnaz -n -b -p  <dig>  on the alignments of pomago and tba.
as rnazwindow.pl assumes lower case nucleotides to be masked, the alignments were converted to upper case letters beforehand;  drna-seq: cdna libraries were mapped with segemehl
<cit>  applying the parameters -m  <dig> -d  <dig> -e  <dig> -p  <dig> -x  <dig> -a  <dig> -e  <dig> .
northern blot experiments
total rna preparation
total rna was prepared from cell pellets using the hot phenol method as described
<cit> .
positive and negative controls
the positive and the negative controls for the northern blot experiments were synthesized by in vitro transcription using the "transcriptaid t <dig> high yield transcription kit" , according to the protocol supplied by the manufacturer.
pcr products generated with the "long pcr enzyme mix"  served as templates for in vitro transcription.
as positive controls for the antisense trna blots, chemically synthesized rna oligonucleotides from "integrated dna technologies"  were used .
rna oligonucleotides were 5’-phosphorylated before gel electrophoresis.
the in vitro transcribed full-length sense trnas  were used as negative controls for the northern blots of antisense trnas.
digoxigenin and lna probes
for the northern blot detection internally digoxigenin-labeled probes were transcribed using the dig rna labeling mix  as described
<cit> .
the antisense trna transcripts were detected with chemically synthesized 5’-digoxigenin-labeled dna/lna mixmer probes .
5’-phosphorylation of rna oligonucleotides
<dig> ng/ μl rna oligonucleotide,  <dig>  mm dtt,  <dig> mm atp and  <dig> u t <dig> polynucleotide kinase  were incubated in  <dig> × t <dig> pnk buffer  in a volume of  <dig> μl for  <dig> h at 37°c, followed by transfer to and storage at -20°c.
electrophoresis
rnas were separated on 8% or 10% denaturing  paa gel with  <dig> × tbe as electrophoresis buffer
<cit> .
blotting, crosslinking, hybridization and detection
rna blotting, hybridization  and immunological detection were performed as described
<cit> , except that rna blotting was carried out at  <dig>  ma/cm  <dig> overnight.
prehybridization and hybridization were performed at 68°c  using  <dig> ml hybridization solution.
<dig>  μl of in vitro transcribed, internally digoxigenin-labeled probe were added for overnight hybridization.
<dig> pmol of chemically synthesized, 5’-digoxigenin-labeled dna/lna mixmer probe were used for northern detection of antisense trnas.
blotted membranes were stored at room temperature.
in vitro transcripts, probes and primers
further details on in vitro transcripts, probes and primers are listed in the supplemental material.
RESULTS
genome analysis – general observations
the genomes of the aquificales range from  <dig>  mb  to  <dig>  mb , thus being at the lower limit of bacterial genomes ranging in size from  <dig>  to  <dig>  mb with a mean of ∼  <dig> mb
<cit> .
the current annotation file of hydrogenivirga sp.
contains  <dig>  mb, which is considerably larger than the genome size of the other aquificales, which might be an assembly artefact as discussed later.
aquificales are known to be at-rich with a gc-content of about 43%
<cit> .
in hydrogenobaculum sp., sulfurihydrogenibium sp.
and s. azorense even only one-third of the nucleotides are guanine or cytosine.
for t. ammonificans an atypically high gc-content of more than 50% was observed.
between  <dig> %  and  <dig> %  of the genomes were found to be unique to each member bacterium .
the comparatively low coverage of hydrogenivirga sp.
is due to the currently assembled genome being almost twice as long as those of other aquificales.
<dig> % to  <dig> % of the pomago alignment, resp.
<dig> % to  <dig> % of the tba alignment, consist of gaps.
according to the wsop each nucleotide from the alignment is conserved on average in slightly less than half of the other  <dig> species  indicating that the genomes diverged relatively fast.
genomic rearrangements among the aquificales, underlining the diversity, can be seen in an overview of the pomago alignment in the supplemental material.
extended annotation of proteins
we extended the original ncbi annotation of proteins of the aquificalesde novo using bacprot, revealing a number of additional proteins .
since a large fraction of proteins are hypothetical or of unknown function, we added for each species a second row which exclusively depicts those with an associated function.
the annotations of ncbi and bacprot were merged to generate an extended annotation of protein genes in the aquificales.
annotations obtained with ncbi  and those identified with bacprot  lead to an extended current annotation of aquificales .
in the second lines, hypothetical proteins were removed.
equal – proteins equally identified by bacprot and ncbi; start/end shifted – proteins identified by bacprot and ncbi vary in length ; ncbi/bacprot only – proteins identified only by ncbi/bacprot.
all gff files are available in the supplemental material.
species abbreviations as in figure
<dig>
we added between  <dig> % of h. thermophilus  and  <dig> % of a. aeolicus  protein-coding genes to the ncbi annotation.
for all proteins annotated by bacprot, we extracted the shine-dalgarno and pribnow box  motifs  in order to facilitate the assignment of novel aquificales-specific proteins.
the shine-dalgarno sequence is rather conserved .
in contrast, the pribnow box is recognizable but less conserved, indicating more sequence variations among promoters.
with the appropriate covariance models we searched for species-specific novel proteins and listed them as predicted proteins in the supplemental material.
<dig>
an overview of the codon usage of a. aeolicus is shown in table
<dig>  complete data on all codon usage tables and mono/dinucleotide distributions are provided in the supplemental material.
we observe a disproportionate usage of certain triplets: isoleucine is mostly  encoded by aua, tyrosine by uac  and histidine by cac .
the four arginine codons with a cytosine at the first position of the triplet are rarely used, compared to the two adenine-containing triplets .
a. aeolicus
codon usage is based on  <dig>  protein-coding genes comprising  <dig>  codons.
codon usage of other aquificales can be viewed in the supplemental material.
aa – amino acid; the fraction of a particular amino acid encoded by the respective codon is given .
homology search and annotation of known ncrnas
a search for ncrna candidates with rnaz <cit>  predicted a relatively constant fraction of the genome to code for ncrnas .
besides the well-known and described rrnas and trnas, only a handful of other wide-spread ncrnas were detected .
rrna operons
most of the aquificales genomes have two rrna operons .
h. thermophilus and t. albus appear to harbor only one operon.
the genomes of t. ammonificans and sulfurihydrogenibium sp.
contain three operons, whereas hydrogenivirga sp.
appears to have two 16s, two 23s and three 5s rrna genes.
trnas
with the exception of hydrogenivirga sp.
, trnascan identified between  <dig>  and  <dig> trnas  per aquificales species.
with splits and aragorn no split trnas were found.
all possible codons are utilized in the aquificales , but the number of trna genes is reduced to a minimum in contrast to reference bacteria such as e. coli which encodes multiple copies of many trna isoacceptors.
figure
<dig> shows nearly no trna with 5’-a in the anticodon and only half of the aquificales have some anticodons with 5’-c, where the non-aquificaceae apparently favored the reduction of such trna genes .
important trna modification enzymes  are encoded in aquificales and x-ray structures of tada and tils from a. aeolicus have been reported
<cit> .
tada converts a residues in the 5’-position of certain trna anticodons to inosine to expand wobble decoding, and tils converts the 5’-c residue in the cau anticodon of specific trna-ile molecules to lysidine  to decode 5’-aua  codons instead of 5’-aug  codons
<cit> .
without this posttranscriptional modification, decoding of isoleucine aua codons would be impossible
<cit> .
selenocysteine-specific trnas decoding 5’-uga are present in the aquificaceae  and in the desulfurobacteriaceae , but are absent from the hydrogenothermaceae .
the aquificaceae , in contrast to the other aquificales or mesophiles such as e. coli or b. subtilis, encode the lysine isoacceptor with the anticodon 5’-cuu to decode the aag codon.
1; right: anticodons specified by the following colors: red – trna with this anticodon encoded in all aquificales; blue – trna encoded in the aquificaceae only.
other colors represent the absence or presence of a trna with this anticodon, as defined in the phylogenetic tree on the left.
for example, trna-sec is present in all aquificaceae except for hydrogenobaculum sp., and is additionally found in the non-aquificaceae species d. thermolithotrophum and t. ammonificans.
rnase p
the catalytic rna subunit of the trna processing endoribonuclease rnase p was previously identified in p. marina and s. azorense <cit> .
additionally, rnase p rnas were easily identified here with bcheck in sulfurihydrogenibium sp., t. ammonificans and d. thermolithotrophum.
in the aquificaceae, rnase p rna candidates were neither detected with bcheck, rnabob nor by manual in silico search methods using cdna libraries of a. aeolicus.
this is consistent with the negative results of previous searches for rnase p rna in a. aeolicus <cit> .
all identified rnase p rnas lack the p <dig> element, which appears to be a general feature of type a rnase p rnas in the hydrogenothermaceae and desulfurobacteriaceae.
the sulfurihydrogenibium sp., t. ammonificans and d. thermolithotrophum rnas differ from their p. marina and s. azorense counterparts by a weaker l9-p <dig> tertiary contact .
other differences are:  very short p <dig> stems in t. ammonificans and d. thermolithotrophum,  particularly weak p <dig> stems in sulfurihydrogenibium sp.
and d. thermolithotrophum,  a destabilized l8-p <dig> interaction, a destabilized p <dig> helix, but a stabilized l14-p <dig> interaction in t. ammonificans.
for details, see rnase p rna 2d structures in the supplemental material.
6s rna
bacterial 6s rnas, about  <dig> nt in length, form a rod-shaped secondary structure with a central bulge region flanked by largely helical arms on both sides.
their structure is thought to mimic the structure of an open dna promoter
<cit> .
6s rnas bind to the housekeeping rna polymerase holoenzyme to block transcription at dna promoters, primarily upon entry of cells into stationary growth phase.
when nutrients are resupplied , rna polymerase massively synthesizes transcripts  on 6s rna as template, which lead to a structural rearrangement of 6s rna and release of rna polymerase.
thus, 6s rna is a fast riboregulator that makes rna polymerase instantly available for a new exponential growth when nutrients are resupplied
<cit> .
in a. aeolicus the 6s rna was clearly identified via an experimental rnomics approach
<cit> .
6s rna candidates in the other aquificales were predicted computationally using the rfam covariance model and, as expected, vary substantially in primary, but less in secondary structure.
for hydrogenivirgia we found two copies.
predicted 6s rnas for t. ammonificans and d. thermolithotrophum remain candidates since they differ substantially from those of other aquificales.
the rnaalifold consensus structure for the 6s rna candidates from all other aquificales analyzed here is shown in the supplement.
individual rnafold predictions  support the notion that they are bona fide 6s rnas.
in the case of a. aeolicus 6s rna, we proposed that formation of a "central bulge collapse" helix  is the major component of the prna-induced rearrangement of this 6s rna structure
<cit> .
if at all, or to which extent, the adjacent hairpin structure forms in the prna-rearranged structure remains to be investigated.
for the eight other 6s rna candidates , we predicted rod-shaped structures with a destabilized central region that is not necessarily purely single-stranded .
according to our proposals, prnas would start with a g residue in the aquificaceae, whereas those of the hydrogenothermaceae  would initiate with an a residue.
<cit>  and the prna transcription start site identified by deep sequencing ; nucleotides of the central bulge region are marked in blue; during prna transcription on 6s rna as template, the endogenous helix is disrupted, leading to the formation of new base-pairing interactions.
here, a 6s rna hybrid with a prna 13-mer  is shown on the right; the proposed rearranged structure of the central 6s rna region
<cit>  has not yet been proven experimentally.
proposed structures of the central bulge regions and their prna-induced rearrangements of the other eight aquificales 6s rna candidates: rearranged structures upon duplex formation with putative prna 13-mers; the prna initiation sites are proposed on the basis of resemblance to a. aeolicus 6s rna.
for more details, see supplemental material.
tmrna
in bacteria, stalling of translating ribosomes on truncated mrnas is rescued through action of the dual-function transfer-messenger rnas
<cit> .
the trna-like domain is present and highly conserved in all aquificales.
an architectural feature of tmrnas is their intricate structure consisting of four pseudoknots.
interestingly, we found two different types of tmrnas, introduced here as type a  and b .
this classification is based on the observation that the lower stem of pseudoknot  <dig>  involves 4- <dig> bp in type a tmrnas, but only 2– <dig> bp in type b variants .
pk <dig> is critical for tmrna function and binds near the ribosomal decoding site
<cit> .
mutational analysis of e. coli tmrna revealed that mutations disrupting the upper stem of pk <dig> are not tolerated, whereas the outer two base pairs of the lower stem  can be disrupted  without loss of function
<cit> .
on the other hand, the tmrna of another thermophile, thermotoga maritima, has a lower pk <dig> stem expanded to  <dig> bp
<cit> .
this raises the question if the aquificales type b tmrnas, for which only a 2-bp lower pk <dig> stem is predicted , are still able to form this pseudoknot, or if the weakness or absence of this stem is compensated for by e.g.
tmrna ligand interactions that are idiosyncratic to the aquificales encoding a type b tmrna.
the messenger rna-like regions , which are in close vicinity of pk <dig>  encode tag preptides of  <dig> amino acids, with subphyla-specific signatures .
for example, all aquificaceae and hydrogenothermaceae tmrnas code for a proline at the second position, which is alanine in the desulfurobacteriaceae.
the genome of hydrogenivirga sp.
appears to encode both types of tmrnas .
whether this reflects a genuine tmrna gene duplication rather than a genome contamination or assembly artefact remains to be clarified .
furthermore, hydrogenobaculum sp.
carries a 78-nt hairpin-like insertion in the pseudoknot  <dig>  region, which however is compatible with formation of pk <dig> .
such a long extension within tmrnas has been not reported yet.
crispr system
for each member of the aquificales we could identify at least one locus of clustered interspaced short palindromic repeat sequences , which are involved in an immunity against viruses and plasmids
<cit> .
although the aquificales have very compact genomes, the number of identified crispr clusters varied from one to thirteen , indicating the presence of thermostable viruses in extreme environments as reported for archaea
<cit> .
the number of crispr clusters does not seem to be clade-specific.
also, the number of repeats in a cluster varies strongly.
for example, in t. albus we found in total four crispr systems containing  <dig>   <dig>   <dig> and  <dig> repeats, whereas in a. aeolicus the five crispr loci only had four to five repeats.
for some, but not all of the crispr clusters, we could detect associated cas genes .
the exact numbers of detected crispr clusters and cas protein cassettes can be seen in figure
<dig>  in this table we included only crispr clusters that were found by both approaches .
it has to be kept in mind that the genome of hydrogenivirga sp.
is in an unfinished state, so it is possible that some crispr loci and especially associated cas genes escaped detection.
other ncrna
srp rna was found once per genome being highly conserved in sequence and structure .
additionally, we show some riboswitch candidates: tpp, moco, cobalamin and crcb .
the moco riboswitch found in t. ammonificans and the two crcb riboswitches identified in hydrogenobaculum sp.
conform well to the rfam conservation model .
riboswitches were only found sporadically among the aquificales.
novel ncrnas in a. aeolicus
besides the annotation of ncrnas with known functions, we additionally aimed to detect novel ncrnas, as they often regulate transcription or play an important role as posttranscriptional regulators.
here we combined in silico analysis of the a. aeolicus genome and drna-seq data from the same organism to identify novel ncrna candidates, some of which were subsequently analyzed by northern blot analysis.
in the in silico search, small ncrnas  were distinguished from proteins by the following analysis steps:  the gc-content of the a. aeolicus genome is 43%.
however, the ncrnas described above show an average gc-content of 66%.
we associated each nucleotide with a local gc-value.
the function of small ncrnas, e.g.
6s rna, is often determined by their stable secondary structure.
to each position in the genome, we assigned the minimum free energy of the most stable local secondary structure including this nucleotide, using rnalfold.
most ncrnas are conserved among closely related organisms.
we calculated genomewide multiple sequence alignments  with tba and pomago of all aquificales genomes, which can be viewed in the supplemental material.
based on the msas we performed a novel ncrna prediction with rnaz and displayed their probability.
all ncrna candidates with a minimum length of  <dig> nt and not overlapping protein-coding sequences, rrna operons or trnas, were summarized in a full candidate table, containing all properties mentioned above .
a subset of these genes can be seen in table
<dig>  we identified  <dig> putative loci for ncrnas, abbreviated n <dig> to n <dig>  all above annotated ncrnas, such as tmrna  or srp rna  were mutually confirmed by our drna-seq and in silico approaches.
interestingly, known ncrnas as well as novel ncrna candidates show a significant level of antisense transcripts .
for unknown ncrnas the sense direction is not assignable.
putative ncrnas, referring to one genomic location and having comparable numbers of cdna read counts on both strands, are described with the same id.
a. aeolicus
the genomic locations and gc-content are listed in columns 2- <dig>  cdna – the maximal number of observed read counts in the - and -library; annotation – overlap to predicted proteins by ncbi and bacprot; rnalfold – energy in kcal/mol of locally stable rna secondary structure predicted by rnalfold; cons_p and cons_t – number of species with homologous regions aligned by pomago and tba; rnaz – probabilities > <dig>   or tba ).
further observations, for example that aae- <dig> was described earlier in
<cit> , are noted in the last column.
a complete list of novel ncrna candidates, and trnas can be found in the supplemental material.
<dig>  for northern blots, see methods.
for comparison reasons, we also added trnas to our table of ncrnas, which show the feature of sense-antisense  expression.
to exclude the possibility of mapping or other artefacts, we confirmed the presence of antisense transcripts exemplarily by northern blots for tmrna and trna  <dig>  .
furthermore, northern blots were conducted for the loci encoding candidates n <dig> and n <dig>  for which the drna-seq data indicated sense and antisense transcription each differing between the - and -library .
for n <dig>  we found most transcripts on the plus strand in the -library , whereas less than half as many transcripts  were detected in the -library.
interestingly, an inverse relation was observed for the minus strand .
for n <dig>  northern blot detection revealed a signal somewhat shorter than the one expected from the cdna read boundaries, whereas no signal could be detected for antisense transcripts .
this finding suggests that the sense transcript is the major one.
in the case of n <dig>  both sense and antisense transcripts of comparable intensity were detected, the major signals of the northern blot representing rnas larger and smaller than anticipated from the read boundaries .
thus, the polarity of the putative ncrna gene remains unclear.
interestingly, very high transcription levels are found in overlapping 5’-upstream regions of two protein-coding genes located on opposing strands .
beside these so-called head-to-head  transcripts we furthermore observed tail-to-tail overlaps  that are represented by very high read coverage .
if these are real transcripts with a certain function or artefacts remains unclear.
CONCLUSIONS
with the advent of a growing number of aquificales genome sequences in public databases, we have re-analyzed this group of bacteria thriving in extreme environments.
the aquificales share the feature of a small, compact genome with a reduced number of protein and ncrna genes.
the genes for trnas are reduced to a minimum but retain the capacity to decode all types of codons, and rrna genes are confined to 2– <dig> copies each.
several classical ncrnas are present, such as srp rna, tmrna, 6s rna, rnase p rna  and riboswitch candidates in some aquificales.
furthermore, by combining in silico analysis with drna-seq data of a. aeolicus, we were able to predict nearly  <dig> novel ncrna candidates, some of which might be specific to the aquificales.
finally, crispr systems of bacterial immunity were identified.
re-annotation of protein genes using bacprot revealed novel proteins with unknown function, some of which might turn out to be specific to the aquificales as well.
on average,  <dig> additional proteins were found that were missing in the respective original annotation.
in our cdna libraries of a. aeolicus, we observed massive amounts of antisense reads with similar patterns  at putative ncrna loci and terminal regions of mrnas.
examples of transcripts antisense to tmrna and trna are illustrated in figure
<dig>
we compared  <dig> bacterial and  <dig> archaeal genomes , and the presence or absence of proteins was used to determine their position in the phylogenetic tree of bacteria.
both archaea form a clear outgroup.
thermodesulfatator indicus branches first in the group of bacteria, followed immediately by the aquificales, while other bacterial branches diverge later.
in an additional protein-based analysis, we took the sequences of single-copy orthologs that were present in at least 50% of all species  .
in contrast to the protein presence/absence tree, neither the aquificales nor t. indicus were placed at a basal position here.
however, the two groups are still in close vicinity to each other.
this analysis not necessarily excludes the possibility of the aquificales being a basal clade.
the selection of orthologs being present in at least 50% of the species leads to a lower coverage of orthologs present in archaea species and therefore may favor long branch attraction
<cit> .
the idea behind selecting frequently occurring single-copy orthologs was to produce phylogenetic trees being less influenced by horizontal gene transfer.
however, proteins shared by archaea and aquificales only are not part of the selected "50% group" of proteins and are therefore not considered in this analysis.
both protein-based phylogenetic trees disagree with a previous study
<cit>  where desulfobacterium autotrophicum hrm <dig>  a δ-proteobacterium, was added to the desulfurobacteriaceae family based on 16s rrna analysis.
we assume that this was an artefact of the high gc-content of rrnas due to the high environmental temperatures.
regarding their proteomes, aquificales and d. autotrophicum are not significantly related.
the results of the 16s rrna phylogenetic analysis did not show a clear picture.
depending on the method used for reconstruction, the aquificales were either placed near the root of the bacterial tree  or not  .
in accordance with the results of
<cit> , the aquificales were always placed close to the thermotogales and thermales-deinococcales, archaea were more closely related to the aquificales than to the thermotogales.
we identified two 6s rna and two tmrna candidate genes in hydrogenivirga sp., rather than a single one as in the other aquificales.
likewise, hydrogenivirga sp.
has a comparatively high amount of trna copies and crispr loci and its genome is estimated to be of roughly double the size of the other aquificales genomes.
combined, these observations support the notion that the hydrogenivirga sp.
genome assembly is erroneous or two genomes of related bacteria  have entered the sequencing project, being in agreement with
<cit> .
based on the tmrna tag peptides identified in the hydrogenivirga sp.
assembly, the second one  matches the sequence exclusively found among the hydrogenothermaceae, although hydrogenivirga sp.
belongs to the aquificaceae .
this suggests that the hydrogenivirga sp.
assembly is a blend of sequences from a member of the aquificaceae and a member of the hydrogenothermaceae.
endnote
anoise cutoff is the highest observed false positive bit score for a potential gene which does not belong to the seed model.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
bioinformatical analysis: ml, sw, kr, bmb, nw, and mm.
experimental validation: ain and rkh.
analyzed data: all.
wrote, read and approved the final manuscript: all.

 aquificales are gram-negative, non-sporulating bacteria that are thermophilic to hyperthermophilic, living in terrestrial and marine hot springs. they are autotrophs that primarily fix carbon by the tricarboxylic acid cycle. the hyperthermophile a. aeolicus, living under extreme temperatures of up to 95°c, has been proposed to have adopted 10% of its protein-coding genes by horizontal gene transfer from archaea. accumulation of all the special properties of thermophiles are rarely understood. special protein-protective mechanisms have been analyzed, but we are far away from a comprehensive understanding of the molecular biology of extremophilic bacteria. beyond idiosyncratic features of aquificales genomes, our interest focussed on their transcriptomes. experimentally, we performed a deep sequencing analysis on the model hyperthermophile a. aeolicus with the primary goal of identifying novel ncrnas candidates. ncrnas are known to have various functions in all domains of life. apart from their general importance as gene expression regulators, ncrnas are involved in processing and translation of other genes, in defending genomes from viral invasion, in shaping and maintenance of bacterial chromosome architecture, and they can even be multifunctional. according to 16s rrna analysis, the aquificales constitute the most deeply rooted bacterial group. however, protein-based phylogenetic reconstructions are not in line with this model. we compared the genomes of the three aquificales families, i.e. aquificaceae, hydrogenothermaceae and desulfurobacteriaceae. we have extended the protein annotation of the mentioned aquificales and reconstructed the phylogenetic position of these species based on 16s rrnas as well as on a set of orthologous proteins. moreover, we have identified ncrnas based on known homologs and present a complete set of novel ncrna candidates based on sequence analyses and deep sequencing data obtained for a. aeolicus. for selected ncrna loci, we provide independent experimental evidence for their expression. methods genomes we analyzed the genomes of the following species split into their respective families: – aquificaceae: aquifex aeolicus vf, hydrogenivirga sp. 128-5-r1-, hydrogenobacter thermophilus tk-, thermocrinis ruber, thermocrinis albus dsm, hydrogenobaculum sp. y04aas, – hydrogenothermaceae: sulfurihydrogenibium sp. yo3aop, sulfurihydrogenibium azorense az-fu, persephonella marina ex-h, and – desulfurobacteriaceae: desulfobacterium thermolithotrophum dsm, and thermovibrio ammonificans hb-. accession numbers and sources of genomes are listed in the electronic supplemental material http://www.rna.uni-jena.de/supplements/aquificales/index.html. whole-genome alignments were constructed using pomago and tba with default parameters. pomago alignments were computed separately for each species as reference. the tba alignment was projected to each of the reference genomes. coverage, alignment quality and gap ratio are given in figure and in the supplemental material. tmrnas are classified into two types. the phylogenetic tree shown at the top of the table is based on the whole genome as well as 16s rrna analysis of the aquificales species. it reproduces the results presented in. for further information, see supplemental material. aae – a. aeolicus, hvi – hydrogenivirga sp., hth – h. thermophilus, hba – hydrogenobaculum sp., tal – t. albus, tru – t. ruber, pma – p. marina, saz – s. azorense, ssp – sulfurihydrogenibium sp., dth – d. thermolithotrophum, tam – t. ammonificans, rs – riboswitch, wsop – weighted sum-of-pairs score, * denotes the hydrogenivirga sp. genome of unfinished assembly. extension of protein annotation we used bacprot to complement the present annotation of protein-coding genes for each aquificales genome. it uses a database of groups of orthologous protein-coding genes present in most bacteria. matches in the genome of interest are annotated, and species-specific features like codon usage, shine-dalgarno sequences, pribnow box motifs and rho-independent terminators are used to predict additional protein-coding genes. to actually achieve a de novo annotation, we excluded all aquificales genes from the reference database. alternative start codons like att and ctg were considered as well. re-annotated and previously annotated proteins and statistics for each species are provided in the supplemental material. annotation of ncrnas by homology we used gorap to annotate ncrnas in the following manner: transfer-rnas were detected by trnascan-se with the option -b for bacteria. split trnas were searched using splits. by applying aragorn, we searched for trnas containing introns. searches for rnase p rna were conducted with bcheck. for the detection of putative crispr loci, crt and crisprfinder were used. we searched for cas protein genes by blast based on known cas genes ). to find further ncrnas, we used blast and infernal. seed sequences from the rfam database and european ribosomal rna database were used as query with an e-value ≤ for blast and the rfam-provided family specific noise cutoffa for infernal. ncrnas expected to escape from detection were searched in a second step with rnabob for short motif search in combination with rnasubopt, rnaduplex, rnacofold, rnaalifold and rnaup from the rna vienna package. for verification, we aligned candidates with clustalw or locarnate. stockholm alignments were adjusted by hand in the emacs ralee mode. resulting stockholm alignments are supplied in the supplemental material in the general feature format as well as in fasta and stockholm formats. phylogenetic reconstruction protein-based phylogeny was performed based on the official ncbi annotations for bacteria shown in the supplemental material. in addition to eleven aquificales species, we included two archaea as outgroup and a wide phylogenetic range of bacterial species representing all bacterial clades. protein sequences were clustered using proteinortho in the blastp+-mode, thus performing a pairwise all-against-all comparison of sequences from different species to derive orthologous relationships. whenever an orthologous group did not have a member in a certain species, we applied tblastn to the respective genome to complement for potentially incomplete annotations. the highest scoring alignment to an orf above a fairly high e-value ≤ 10- was added to the initial protein annotation. finally, proteinortho was applied again using the expanded annotation. for a high resolution phylogeny within the aquificales, we created a whole genome alignment using pomago. the alignment was analyzed using raxml with a gamma model of rate heterogeneity with an estimate on the proportion of invariable sites and rapid bootstraps. in an additional phylogenetic analysis we used single-copy orthologous proteins present in at least 50% of all species in the set. each protein group was aligned separately using dialign-tx. both ends of the group’s alignments were cropped to remove leading and tailing gaps. the remaining sequences were concatenated resulting in a aa long alignment and applied to raxml using the lg substitution model as well as the gamma model of rate heterogeneity with rapid bootstraps. the 16s rrna-based phylogeny was computed with mafft using the l-ins-i method with iterations. we used different approaches: neighbor joining with the kimura correction model, bayesian inference with mrbayes with default parameters, maximum likelihood with raxml with the base substitution models gtrgamma and gtrcat for the bootstrapping phase. for all previously mentioned methods the archaea methanobacterium sp. al- and archaeoglobus fulgidus were used as outgroup. as state of the art, we have estimated a tree with sate. related sequences were aligned with mafft and subsequently merged by muscle. the tree was computed using raxml. drna-seq of a. aeolicus total cellular rna transcriptome analysis of a. aeolicus was based on cdna libraries from a differential deep sequencing approach. a. aeolicus cells, provided by m. thomm and r. huber, were grown for day and harvested as described. for preparation of total cellular rna, we used the hot phenol method: cell pellets were resuspended in extraction buffer and incubated for min at room temperature with volumes of lysozyme. sds was added to a final concentration of 1% followed by vigorous vortexing. after addition of volume phenol and vortexing, the mixture was incubated for min at 65°c, then cooled on ice for min, and centrifuged for min at 4°c and g. phenol extraction was repeated, followed by chloroform extraction and ethanol precipitation. finally, the dna was digested with u turbo dnase for min at 37°c, followed by addition of another u dnase and incubation for another min at 37°c. subsequently, the rna was subjected to phenol/chloroform extraction and ethanol precipitation. after redissolving the rna in double-distilled water, its concentration was determined by uv spectroscopy. before cdna library construction, the rna was split into two fractions; one fraction was treated with terminator 5’ p-dependent exonuclease for depletion of transcripts carrying a 5’-monophosphate. both fractions were treated with tobacco acid phosphatase before 5’-linker ligation, poly tailing and conversion into cdna. the cdna libraries were then sequenced on a roche flx sequencer and resulted in the -library with reads and the -library containing the enriched primary transcripts. detection of novel ncrnas we used the igb to visualize the following features of a. aeolicus: nucleotide sequence; local gc-content; protein genes annotated by ncbi and bacprot; locally stable secondary structures: calculation was performed with rnalfold with options -d and -l for both strands with a maximum base-pair span of nucleotides. sequences with local structures of fewer than nt were discarded. for the prediction of thermodynamically stable rna structures, each sequence was shuffled times while preserving the dinucleotide frequencies; to classify extraordinarily stable rna secondary structures, we chose to use a z-score cutoff of -; conserved regions among the aquificales: with default parameters of tba and pomago we aligned genomes; the tba alignment was projected to each of the reference genomes; coverage, wsop and gap ratio are given in figure 1; novel ncrnas: novel ncrna candidates were predicted using rnaz. we used rnazwindow.pl –min-seqs= and rnaz -n -b -p on the alignments of pomago and tba. as rnazwindow.pl assumes lower case nucleotides to be masked, the alignments were converted to upper case letters beforehand; drna-seq: cdna libraries were mapped with segemehl applying the parameters -m -d -e -p -x -a -e. northern blot experiments total rna preparation total rna was prepared from cell pellets using the hot phenol method as described. positive and negative controls the positive and the negative controls for the northern blot experiments were synthesized by in vitro transcription using the "transcriptaid t high yield transcription kit", according to the protocol supplied by the manufacturer. pcr products generated with the "long pcr enzyme mix" served as templates for in vitro transcription. as positive controls for the antisense trna blots, chemically synthesized rna oligonucleotides from "integrated dna technologies" were used. rna oligonucleotides were 5’-phosphorylated before gel electrophoresis. the in vitro transcribed full-length sense trnas were used as negative controls for the northern blots of antisense trnas. digoxigenin and lna probes for the northern blot detection internally digoxigenin-labeled probes were transcribed using the dig rna labeling mix as described. the antisense trna transcripts were detected with chemically synthesized 5’-digoxigenin-labeled dna/lna mixmer probes. 5’-phosphorylation of rna oligonucleotides ng/ μl rna oligonucleotide, mm dtt, mm atp and u t polynucleotide kinase were incubated in × t pnk buffer in a volume of μl for h at 37°c, followed by transfer to and storage at -20°c. electrophoresis rnas were separated on 8% or 10% denaturing paa gel with × tbe as electrophoresis buffer. blotting, crosslinking, hybridization and detection rna blotting, hybridization and immunological detection were performed as described, except that rna blotting was carried out at ma/cm overnight. prehybridization and hybridization were performed at 68°c using ml hybridization solution. μl of in vitro transcribed, internally digoxigenin-labeled probe were added for overnight hybridization. pmol of chemically synthesized, 5’-digoxigenin-labeled dna/lna mixmer probe were used for northern detection of antisense trnas. blotted membranes were stored at room temperature. in vitro transcripts, probes and primers further details on in vitro transcripts, probes and primers are listed in the supplemental material.  genome analysis – general observations the genomes of the aquificales range from mb to mb, thus being at the lower limit of bacterial genomes ranging in size from to mb with a mean of ∼ mb. the current annotation file of hydrogenivirga sp. contains mb, which is considerably larger than the genome size of the other aquificales, which might be an assembly artefact as discussed later. aquificales are known to be at-rich with a gc-content of about 43%. in hydrogenobaculum sp., sulfurihydrogenibium sp. and s. azorense even only one-third of the nucleotides are guanine or cytosine. for t. ammonificans an atypically high gc-content of more than 50% was observed. between % and % of the genomes were found to be unique to each member bacterium. the comparatively low coverage of hydrogenivirga sp. is due to the currently assembled genome being almost twice as long as those of other aquificales. % to % of the pomago alignment, resp. % to % of the tba alignment, consist of gaps. according to the wsop each nucleotide from the alignment is conserved on average in slightly less than half of the other species indicating that the genomes diverged relatively fast. genomic rearrangements among the aquificales, underlining the diversity, can be seen in an overview of the pomago alignment in the supplemental material. extended annotation of proteins we extended the original ncbi annotation of proteins of the aquificalesde novo using bacprot, revealing a number of additional proteins. since a large fraction of proteins are hypothetical or of unknown function, we added for each species a second row which exclusively depicts those with an associated function. the annotations of ncbi and bacprot were merged to generate an extended annotation of protein genes in the aquificales. annotations obtained with ncbi and those identified with bacprot lead to an extended current annotation of aquificales. in the second lines, hypothetical proteins were removed. equal – proteins equally identified by bacprot and ncbi; start/end shifted – proteins identified by bacprot and ncbi vary in length; ncbi/bacprot only – proteins identified only by ncbi/bacprot. all gff files are available in the supplemental material. species abbreviations as in figure we added between % of h. thermophilus and % of a. aeolicus protein-coding genes to the ncbi annotation. for all proteins annotated by bacprot, we extracted the shine-dalgarno and pribnow box motifs in order to facilitate the assignment of novel aquificales-specific proteins. the shine-dalgarno sequence is rather conserved. in contrast, the pribnow box is recognizable but less conserved, indicating more sequence variations among promoters. with the appropriate covariance models we searched for species-specific novel proteins and listed them as predicted proteins in the supplemental material. an overview of the codon usage of a. aeolicus is shown in table complete data on all codon usage tables and mono/dinucleotide distributions are provided in the supplemental material. we observe a disproportionate usage of certain triplets: isoleucine is mostly encoded by aua, tyrosine by uac and histidine by cac. the four arginine codons with a cytosine at the first position of the triplet are rarely used, compared to the two adenine-containing triplets. a. aeolicus codon usage is based on protein-coding genes comprising codons. codon usage of other aquificales can be viewed in the supplemental material. aa – amino acid; the fraction of a particular amino acid encoded by the respective codon is given. homology search and annotation of known ncrnas a search for ncrna candidates with rnaz predicted a relatively constant fraction of the genome to code for ncrnas. besides the well-known and described rrnas and trnas, only a handful of other wide-spread ncrnas were detected. rrna operons most of the aquificales genomes have two rrna operons. h. thermophilus and t. albus appear to harbor only one operon. the genomes of t. ammonificans and sulfurihydrogenibium sp. contain three operons, whereas hydrogenivirga sp. appears to have two 16s, two 23s and three 5s rrna genes. trnas with the exception of hydrogenivirga sp., trnascan identified between and trnas per aquificales species. with splits and aragorn no split trnas were found. all possible codons are utilized in the aquificales, but the number of trna genes is reduced to a minimum in contrast to reference bacteria such as e. coli which encodes multiple copies of many trna isoacceptors. figure shows nearly no trna with 5’-a in the anticodon and only half of the aquificales have some anticodons with 5’-c, where the non-aquificaceae apparently favored the reduction of such trna genes. important trna modification enzymes are encoded in aquificales and x-ray structures of tada and tils from a. aeolicus have been reported. tada converts a residues in the 5’-position of certain trna anticodons to inosine to expand wobble decoding, and tils converts the 5’-c residue in the cau anticodon of specific trna-ile molecules to lysidine to decode 5’-aua codons instead of 5’-aug codons. without this posttranscriptional modification, decoding of isoleucine aua codons would be impossible. selenocysteine-specific trnas decoding 5’-uga are present in the aquificaceae and in the desulfurobacteriaceae, but are absent from the hydrogenothermaceae. the aquificaceae, in contrast to the other aquificales or mesophiles such as e. coli or b. subtilis, encode the lysine isoacceptor with the anticodon 5’-cuu to decode the aag codon. 1; right: anticodons specified by the following colors: red – trna with this anticodon encoded in all aquificales; blue – trna encoded in the aquificaceae only. other colors represent the absence or presence of a trna with this anticodon, as defined in the phylogenetic tree on the left. for example, trna-sec is present in all aquificaceae except for hydrogenobaculum sp., and is additionally found in the non-aquificaceae species d. thermolithotrophum and t. ammonificans. rnase p the catalytic rna subunit of the trna processing endoribonuclease rnase p was previously identified in p. marina and s. azorense. additionally, rnase p rnas were easily identified here with bcheck in sulfurihydrogenibium sp., t. ammonificans and d. thermolithotrophum. in the aquificaceae, rnase p rna candidates were neither detected with bcheck, rnabob nor by manual in silico search methods using cdna libraries of a. aeolicus. this is consistent with the negative results of previous searches for rnase p rna in a. aeolicus. all identified rnase p rnas lack the p element, which appears to be a general feature of type a rnase p rnas in the hydrogenothermaceae and desulfurobacteriaceae. the sulfurihydrogenibium sp., t. ammonificans and d. thermolithotrophum rnas differ from their p. marina and s. azorense counterparts by a weaker l9-p tertiary contact. other differences are: very short p stems in t. ammonificans and d. thermolithotrophum, particularly weak p stems in sulfurihydrogenibium sp. and d. thermolithotrophum, a destabilized l8-p interaction, a destabilized p helix, but a stabilized l14-p interaction in t. ammonificans. for details, see rnase p rna 2d structures in the supplemental material. 6s rna bacterial 6s rnas, about nt in length, form a rod-shaped secondary structure with a central bulge region flanked by largely helical arms on both sides. their structure is thought to mimic the structure of an open dna promoter. 6s rnas bind to the housekeeping rna polymerase holoenzyme to block transcription at dna promoters, primarily upon entry of cells into stationary growth phase. when nutrients are resupplied, rna polymerase massively synthesizes transcripts on 6s rna as template, which lead to a structural rearrangement of 6s rna and release of rna polymerase. thus, 6s rna is a fast riboregulator that makes rna polymerase instantly available for a new exponential growth when nutrients are resupplied. in a. aeolicus the 6s rna was clearly identified via an experimental rnomics approach. 6s rna candidates in the other aquificales were predicted computationally using the rfam covariance model and, as expected, vary substantially in primary, but less in secondary structure. for hydrogenivirgia we found two copies. predicted 6s rnas for t. ammonificans and d. thermolithotrophum remain candidates since they differ substantially from those of other aquificales. the rnaalifold consensus structure for the 6s rna candidates from all other aquificales analyzed here is shown in the supplement. individual rnafold predictions support the notion that they are bona fide 6s rnas. in the case of a. aeolicus 6s rna, we proposed that formation of a "central bulge collapse" helix is the major component of the prna-induced rearrangement of this 6s rna structure. if at all, or to which extent, the adjacent hairpin structure forms in the prna-rearranged structure remains to be investigated. for the eight other 6s rna candidates, we predicted rod-shaped structures with a destabilized central region that is not necessarily purely single-stranded. according to our proposals, prnas would start with a g residue in the aquificaceae, whereas those of the hydrogenothermaceae would initiate with an a residue. and the prna transcription start site identified by deep sequencing; nucleotides of the central bulge region are marked in blue; during prna transcription on 6s rna as template, the endogenous helix is disrupted, leading to the formation of new base-pairing interactions. here, a 6s rna hybrid with a prna 13-mer is shown on the right; the proposed rearranged structure of the central 6s rna region has not yet been proven experimentally. proposed structures of the central bulge regions and their prna-induced rearrangements of the other eight aquificales 6s rna candidates: rearranged structures upon duplex formation with putative prna 13-mers; the prna initiation sites are proposed on the basis of resemblance to a. aeolicus 6s rna. for more details, see supplemental material. tmrna in bacteria, stalling of translating ribosomes on truncated mrnas is rescued through action of the dual-function transfer-messenger rnas. the trna-like domain is present and highly conserved in all aquificales. an architectural feature of tmrnas is their intricate structure consisting of four pseudoknots. interestingly, we found two different types of tmrnas, introduced here as type a and b. this classification is based on the observation that the lower stem of pseudoknot involves 4- bp in type a tmrnas, but only 2– bp in type b variants. pk is critical for tmrna function and binds near the ribosomal decoding site. mutational analysis of e. coli tmrna revealed that mutations disrupting the upper stem of pk are not tolerated, whereas the outer two base pairs of the lower stem can be disrupted without loss of function. on the other hand, the tmrna of another thermophile, thermotoga maritima, has a lower pk stem expanded to bp. this raises the question if the aquificales type b tmrnas, for which only a 2-bp lower pk stem is predicted, are still able to form this pseudoknot, or if the weakness or absence of this stem is compensated for by e.g. tmrna ligand interactions that are idiosyncratic to the aquificales encoding a type b tmrna. the messenger rna-like regions, which are in close vicinity of pk encode tag preptides of amino acids, with subphyla-specific signatures. for example, all aquificaceae and hydrogenothermaceae tmrnas code for a proline at the second position, which is alanine in the desulfurobacteriaceae. the genome of hydrogenivirga sp. appears to encode both types of tmrnas. whether this reflects a genuine tmrna gene duplication rather than a genome contamination or assembly artefact remains to be clarified. furthermore, hydrogenobaculum sp. carries a 78-nt hairpin-like insertion in the pseudoknot region, which however is compatible with formation of pk. such a long extension within tmrnas has been not reported yet. crispr system for each member of the aquificales we could identify at least one locus of clustered interspaced short palindromic repeat sequences, which are involved in an immunity against viruses and plasmids. although the aquificales have very compact genomes, the number of identified crispr clusters varied from one to thirteen, indicating the presence of thermostable viruses in extreme environments as reported for archaea. the number of crispr clusters does not seem to be clade-specific. also, the number of repeats in a cluster varies strongly. for example, in t. albus we found in total four crispr systems containing and repeats, whereas in a. aeolicus the five crispr loci only had four to five repeats. for some, but not all of the crispr clusters, we could detect associated cas genes. the exact numbers of detected crispr clusters and cas protein cassettes can be seen in figure in this table we included only crispr clusters that were found by both approaches. it has to be kept in mind that the genome of hydrogenivirga sp. is in an unfinished state, so it is possible that some crispr loci and especially associated cas genes escaped detection. other ncrna srp rna was found once per genome being highly conserved in sequence and structure. additionally, we show some riboswitch candidates: tpp, moco, cobalamin and crcb. the moco riboswitch found in t. ammonificans and the two crcb riboswitches identified in hydrogenobaculum sp. conform well to the rfam conservation model. riboswitches were only found sporadically among the aquificales. novel ncrnas in a. aeolicus besides the annotation of ncrnas with known functions, we additionally aimed to detect novel ncrnas, as they often regulate transcription or play an important role as posttranscriptional regulators. here we combined in silico analysis of the a. aeolicus genome and drna-seq data from the same organism to identify novel ncrna candidates, some of which were subsequently analyzed by northern blot analysis. in the in silico search, small ncrnas were distinguished from proteins by the following analysis steps: the gc-content of the a. aeolicus genome is 43%. however, the ncrnas described above show an average gc-content of 66%. we associated each nucleotide with a local gc-value. the function of small ncrnas, e.g. 6s rna, is often determined by their stable secondary structure. to each position in the genome, we assigned the minimum free energy of the most stable local secondary structure including this nucleotide, using rnalfold. most ncrnas are conserved among closely related organisms. we calculated genomewide multiple sequence alignments with tba and pomago of all aquificales genomes, which can be viewed in the supplemental material. based on the msas we performed a novel ncrna prediction with rnaz and displayed their probability. all ncrna candidates with a minimum length of nt and not overlapping protein-coding sequences, rrna operons or trnas, were summarized in a full candidate table, containing all properties mentioned above. a subset of these genes can be seen in table we identified putative loci for ncrnas, abbreviated n to n all above annotated ncrnas, such as tmrna or srp rna were mutually confirmed by our drna-seq and in silico approaches. interestingly, known ncrnas as well as novel ncrna candidates show a significant level of antisense transcripts. for unknown ncrnas the sense direction is not assignable. putative ncrnas, referring to one genomic location and having comparable numbers of cdna read counts on both strands, are described with the same id. a. aeolicus the genomic locations and gc-content are listed in columns 2- cdna – the maximal number of observed read counts in the - and -library; annotation – overlap to predicted proteins by ncbi and bacprot; rnalfold – energy in kcal/mol of locally stable rna secondary structure predicted by rnalfold; cons_p and cons_t – number of species with homologous regions aligned by pomago and tba; rnaz – probabilities > or tba ). further observations, for example that aae- was described earlier in, are noted in the last column. a complete list of novel ncrna candidates, and trnas can be found in the supplemental material. for northern blots, see methods. for comparison reasons, we also added trnas to our table of ncrnas, which show the feature of sense-antisense expression. to exclude the possibility of mapping or other artefacts, we confirmed the presence of antisense transcripts exemplarily by northern blots for tmrna and trna. furthermore, northern blots were conducted for the loci encoding candidates n and n for which the drna-seq data indicated sense and antisense transcription each differing between the - and -library. for n we found most transcripts on the plus strand in the -library, whereas less than half as many transcripts were detected in the -library. interestingly, an inverse relation was observed for the minus strand. for n northern blot detection revealed a signal somewhat shorter than the one expected from the cdna read boundaries, whereas no signal could be detected for antisense transcripts. this finding suggests that the sense transcript is the major one. in the case of n both sense and antisense transcripts of comparable intensity were detected, the major signals of the northern blot representing rnas larger and smaller than anticipated from the read boundaries. thus, the polarity of the putative ncrna gene remains unclear. interestingly, very high transcription levels are found in overlapping 5’-upstream regions of two protein-coding genes located on opposing strands. beside these so-called head-to-head transcripts we furthermore observed tail-to-tail overlaps that are represented by very high read coverage. if these are real transcripts with a certain function or artefacts remains unclear.  with the advent of a growing number of aquificales genome sequences in public databases, we have re-analyzed this group of bacteria thriving in extreme environments. the aquificales share the feature of a small, compact genome with a reduced number of protein and ncrna genes. the genes for trnas are reduced to a minimum but retain the capacity to decode all types of codons, and rrna genes are confined to 2– copies each. several classical ncrnas are present, such as srp rna, tmrna, 6s rna, rnase p rna and riboswitch candidates in some aquificales. furthermore, by combining in silico analysis with drna-seq data of a. aeolicus, we were able to predict nearly novel ncrna candidates, some of which might be specific to the aquificales. finally, crispr systems of bacterial immunity were identified. re-annotation of protein genes using bacprot revealed novel proteins with unknown function, some of which might turn out to be specific to the aquificales as well. on average, additional proteins were found that were missing in the respective original annotation. in our cdna libraries of a. aeolicus, we observed massive amounts of antisense reads with similar patterns at putative ncrna loci and terminal regions of mrnas. examples of transcripts antisense to tmrna and trna are illustrated in figure we compared bacterial and archaeal genomes, and the presence or absence of proteins was used to determine their position in the phylogenetic tree of bacteria. both archaea form a clear outgroup. thermodesulfatator indicus branches first in the group of bacteria, followed immediately by the aquificales, while other bacterial branches diverge later. in an additional protein-based analysis, we took the sequences of single-copy orthologs that were present in at least 50% of all species. in contrast to the protein presence/absence tree, neither the aquificales nor t. indicus were placed at a basal position here. however, the two groups are still in close vicinity to each other. this analysis not necessarily excludes the possibility of the aquificales being a basal clade. the selection of orthologs being present in at least 50% of the species leads to a lower coverage of orthologs present in archaea species and therefore may favor long branch attraction. the idea behind selecting frequently occurring single-copy orthologs was to produce phylogenetic trees being less influenced by horizontal gene transfer. however, proteins shared by archaea and aquificales only are not part of the selected "50% group" of proteins and are therefore not considered in this analysis. both protein-based phylogenetic trees disagree with a previous study where desulfobacterium autotrophicum hrm a δ-proteobacterium, was added to the desulfurobacteriaceae family based on 16s rrna analysis. we assume that this was an artefact of the high gc-content of rrnas due to the high environmental temperatures. regarding their proteomes, aquificales and d. autotrophicum are not significantly related. the results of the 16s rrna phylogenetic analysis did not show a clear picture. depending on the method used for reconstruction, the aquificales were either placed near the root of the bacterial tree or not. in accordance with the results of, the aquificales were always placed close to the thermotogales and thermales-deinococcales, archaea were more closely related to the aquificales than to the thermotogales. we identified two 6s rna and two tmrna candidate genes in hydrogenivirga sp., rather than a single one as in the other aquificales. likewise, hydrogenivirga sp. has a comparatively high amount of trna copies and crispr loci and its genome is estimated to be of roughly double the size of the other aquificales genomes. combined, these observations support the notion that the hydrogenivirga sp. genome assembly is erroneous or two genomes of related bacteria have entered the sequencing project, being in agreement with. based on the tmrna tag peptides identified in the hydrogenivirga sp. assembly, the second one matches the sequence exclusively found among the hydrogenothermaceae, although hydrogenivirga sp. belongs to the aquificaceae. this suggests that the hydrogenivirga sp. assembly is a blend of sequences from a member of the aquificaceae and a member of the hydrogenothermaceae. endnote anoise cutoff is the highest observed false positive bit score for a potential gene which does not belong to the seed model. competing interests the authors declare that they have no competing interests. authors’ contributions bioinformatical analysis: ml, sw, kr, bmb, nw, and mm. experimental validation: ain and rkh. analyzed data: all. wrote, read and approved the final manuscript: all. 
by combining in silico analysis with drna-seq data of a. aeolicus we predict nearly novel ncrna candidates. we present several classical ncrnas and riboswitch candidates. in this work, we compare the available genomes, extend their protein annotation, identify regulatory sequences, annotate non-coding rnas of known function, predict novel ncrna candidates, show idiosyncrasies of the genetic decoding machinery, present two different types of transfer-messenger rnas and variations of the crispr systems. furthermore, we performed a phylogenetic analysis of the aquificales based on entire genome sequences, and extended this by a classification among all bacteria using 16s rrna sequences and a set of orthologous proteins. we have here re-analyzed the aquificales, a group of bacteria thriving in extreme environments, sharing the feature of a small, compact genome with a reduced number of protein and ncrna genes. combining several in silico features with an in vivo drna-seq transcriptome analysis of aquifex aeolicus, we predict roughly novel ncrna candidates in this bacterium.
 the aquificales are a diverse group of thermophilic bacteria that thrive in terrestrial and marine hydrothermal environments. they can be divided into the families aquificaceae, desulfurobacteriaceae and hydrogenothermaceae. although eleven fully sequenced and assembled genomes are available, only little is known about this taxonomic order in terms of rna metabolism.  in this work, we compare the available genomes, extend their protein annotation, identify regulatory sequences, annotate non-coding rnas of known function, predict novel ncrna candidates, show idiosyncrasies of the genetic decoding machinery, present two different types of transfer-messenger rnas and variations of the crispr systems. furthermore, we performed a phylogenetic analysis of the aquificales based on entire genome sequences, and extended this by a classification among all bacteria using 16s rrna sequences and a set of orthologous proteins. combining several in silico features with an in vivo drna-seq transcriptome analysis of aquifex aeolicus, we predict roughly novel ncrna candidates in this bacterium.  we have here re-analyzed the aquificales, a group of bacteria thriving in extreme environments, sharing the feature of a small, compact genome with a reduced number of protein and ncrna genes. we present several classical ncrnas and riboswitch candidates. by combining in silico analysis with drna-seq data of a. aeolicus we predict nearly novel ncrna candidates. aquificalesthermophilesncrnaaquificaceaedesulfurobacteriaceaehydrogenothermaceae 