text:
A high-quality human reference panel reveals the complexity and distribution of genomic structural variants

Structural variation (SV) represents a major source of differences between individual human genomes and has been linked to disease phenotypes. However, the majority of studies provide neither a global view of the full spectrum of these variants nor integrate them into reference panels of genetic variation. Here, we analyse whole genome sequencing data of 769 individuals from 250 Dutch families, and provide a haplotype-resolved map of 1.9 million genome variants across 9 different variant classes, including novel forms of complex indels, and retrotransposition-mediated insertions of mobile elements and processed RNAs. A large proportion are previously under reported variants sized between 21 and 100 bp. We detect 4 megabases of novel sequence, encoding 11 new transcripts. Finally, we show 191 known, trait-associated SNPs to be in strong linkage disequilibrium with SVs and demonstrate that our panel facilitates accurate imputation of SVs in unrelated individuals.

Deletions have been stratified into four bins based on AF quartiles (panels a-d). Each plot shows the length distribution in the given AF range. Furthermore, thin bars indicate the approaches that have contributed to the discovery of the deletions in the respective size and AF range. The same deletion can have been discovered by multiple approaches, so each five thin bars can add up to more than the enclosing grey bar.   i.e. of a base-line random caller, are shown (black). Figure 11 Supplementary Figure 11. Distribution of imputed genotypes per 'ground truth' genotype and per SV type when only genotypes for which the gold standard genotype contains the rare allele are considered. Note that either the reference or the alternative allele can be the rare(r) allele. RR, RA, and AA are short for "homozygous in the reference allele", "heterozygous", and "homozygous in the alternative allele". Imputed genotypes have been filtered at a GL threshold of c=0.95.  Table 3 Overview of classes of novel genomic segments identified in the GoNL dataset.

Indel and structural variant calling
Variant identification was performed using a combination of 12 different algorithms, which represent 5 major approaches for indel and structural variation (SV) discovery: (i) gapped reads, (ii) split-read, (iii) read pair, (iv) read depth and (v) de novo assembly (Figure 1)  was filtered to only include reads of which one end could not be mapped to the reference genome, significantly reducing the file size for each region and resulting in less transfer overhead. In the second step of the analysis, one region was selected and analyzed across all samples. For each analysis 8 compute cores were used, consuming 5000 core hours in total. The confident variants are those observed in one child with at least 2 supporting reads and which also appeared in at least one of the parents.

GATK Unified Genotyper
We used the GATK UnifiedGenotyper (UG; https://software.broadinstitute.org/gatk/) v2.1.8 for indel discovery and calculation of genotype likelihoods, running it on 5 Mb regions with 100 bp overlap on each side. The overlap ensured that an indel at the fringe of the 5 Mb region would be correctly called. In order to separate biallelic sites from multi-allelic ones, we used the GATK UG multi-allelic model for indel calling.
Only bi-allelic ones were kept for further processing and downstream analyses. We only reported indels with a Phred quality score > Q10.
We adopted a filtering strategy by using GATK VQSR on our initial calls. We used the iii.
Haplotype score iv.
Read position rank sum test v.

Inbreeding coefficient
We genotyped all autosomal calls passing VQSR with GATK PhaseByTransmission using a mutation prior of 10 -4 to obtain sensitive trio-aware genotypes.

1--2--3--SV
We used version 0.9 of the 1-2-3-SV software package (http://tools.genomes.nl/123sv.html) 5 for population-based calling of deletions, inversions, intra-and inter-chromosomal structural variants. The analysis was performed across the complete set of 2,238 libraries. We considered only uniquely mapped reads (BWA X0 flag = 1) that do not have a secondary hit (X1 flag = 0). A read pair was considered to be discordant if the two reads in a pair map to different chromosomes, display incorrect mapping orientation or exhibited too short (below 0.1 percentile) or too long (above 99.9 percentile of fragment size distribution characteristic for a given library) distance of the mapping positions. We only retained clusters supported by 4 or more discordant and independent (not PCR duplicates) read pairs.

Breakdancer
We used BreakDancer (http://breakdancer.sourceforge.net) version 1.1 with parameters -r 5c 7 6 . The rest of the parameters were set to default). Calling was performed per family. Calls from different families that have overlapping confidence intervals of both breakpoints were merged in an iterative fashion, using a custom Perl script. Clusters supported by at least five read pairs were retained for further analysis.

DWAC--Seq
Version 0.7 of DWAC-Seq (http://tools.genomes.nl/dwac-seq.html) was used to call segmental copy-number changes. A BAM file representing a common GoNL reference set was constructed by taking 0.13% (1/769) reads from each of the individual samples and merging them into a single file. Each GoNL sample was then compared against this common reference using dynamic windows that overlap 250 unambiguously mapped sequencing reads (BWA flag X0=1). We filtered out CNV candidates that showed a copy-number change of less than 30%. CNV calls from different families with 80% reciprocal overlap were merged together in an iterative manner using a custom Perl script.

CNVnator
The CNVnator (http://sv.gersteinlab.org/cnvnator/) v0.2.2 read-depth algorithm was used for calling CNVs 7 . The read-depth signal consists of readcounts in non-overlapping genomic bins of equal size. Subsequently, the read-depth signal is then segmented into regions with different copy numbers. CNVs calls are predicted on the basis of statistical significance tests performed on the segmented regions.
To determine the optimal bin size for our dataset, we analyzed a quartet with MZ twins using CNVnator with nine different bin sizes (ranging from 100 -500 bp, with an increment of 50 bp). We found that a binsize of 500 bp resulted in the highest concordance between the MZ twin samples and therefore used this setting for the full data set.
We applied a GC correction and included CNVs in the final call on the basis of the following criteria: (i) p < 0.05 (corrected for multiple testing), (ii) the proportion of reads in the CNV that map to > 1 location does not exceed 50%, and (iii) the CNV does not overlap with an assembly gap of the reference genome or centromere.

FACADE
Short reads were mapped using mrsFAST to a repeat-masked reference (hg19) 8 . A correction for GC bias was applied to the raw read depth. An initial quality control step involved calculation of the correlation (R 2 ) between copy numbers and coverage at known loci with a fixed genotype. One sample with R 2 < 0.85 was withdrawn from further analysis.
With read depth information, we could predict the copy number of a gene or a segment of the genome. As a second quality control step, we performed copy number calls with 3 kb windows across the whole genome. Samples that had > 120,000 windows with a copy number of 2 were discarded, leading to the exclusion of 30 samples. In total, 224 families including all 19 quartets met the criteria and were used for variant calling.
To discover CNV events, we first estimated copy numbers across the whole genome from read-depth analysis using windows of 1kb of unique sequence. Then, copy number log 2 ratios between offspring and parents as well as all samples versus the population median were calculated. All pairwise log 2 ratios were segmented using the previously published method FACADE to produce the initial call set 8 . To account for variability in call boundaries with events encompassing a small number of windows of variable genomic size, we applied a Multiscale French Tophat wavelet to refine the breakpoints for events between 2 and 7 windows in size. Each call was classified into four classes: (i) de novo, (ii) assortment, (iii) inherited, and (iv) common to the family. To filter out low-confident calls, only calls spanning ≥ 2 windows and with a log ratio greater than 0.31 or less than -0.25 for duplications and deletions, respectively, were retained.

Mate--Clever
For deletion discovery, we ran the discovery part of MATE-CLEVER 9 , with minor modifications that account for different insert size distributions in different libraries.
MATE-CLEVER is an integrated approach. Its major purpose in the frame of the project is to discover deletions of size 30--100 bp (sometimes termed the "twilight zone of NGS indels").
It incorporates CLEVER 10 , as an internal segment size based approach and LASER 9 , as a split-read aligner.
MATE-CLEVER uses CLEVER in a first step, in order to discover deletions of size 30--100 bp at extremely high sensitivity, and uses LASER in a second step, in order to refine the breakpoint annotations made by CLEVER. It also uses several auxiliary tools, as described below. In the following, we provide the full description of details and commands, by which to reproduce MATE-CLEVER's callset.
To run the CLEVER-based deletion discovery pipeline, revision 3097f2 from the git repository at https://bitbucket.org/tobiasmarschall/clever-toolkit has been used. It includes CLEVER, LASER, and several auxiliary tools described below. To run (and parallelize) the below pipeline, we used the Python-based workflow engine Snakemake 11 .
CLEVER was run on all individuals separately. In detail, that means the following steps were performed: 1. We used tool \texttt{bwabam-to-alignment-priors} to extract prior alignment probabilities. The output was split by chromosome. Commandline: "bwabam-toalignment-priors -m mean-and-sd.txt ref.fasta input.bam | split-priors-by-chromosome -z -s input.bam priors". For each chromosome and each sample in the input, one gzipped file prefixed with "priors" is created. All alternative alignments provided by BWA through XA tags are used. Furthermore, this step generates robust estimates of insert size mean and standard deviation, which are stored in the file mean-and-sd.txt.
2. All files with alignment probabilities are sorted by genomic position by means of a standard unix sort (sort -g -k7).
3. The CLEVER core engine is run on each of these files. Command line: "zcat input.probabilities.gz | clever -c 150 -v > clever.output". The option "-c 150" instructs CLEVER to limit local coverage to 150 to avoid long runtimes for regions with excess coverage (due to repeats). LASER was run on all regions with putative deletions as identified by CLEVER in Step 1: 1. To generate high-quality background distributions for insert size, insertion size, and deletion size from uniquely mapping reads, LASER was used to realign reads from randomly chosen regions. To this end, 5000 regions of length 10000 were sampled uniformly at random. Reads mapped to these regions by BWA were extracted and remapped using LASER as described previously 9 . The following command line was used for the LASER core step: "join-split-reads -XIS -A14 --anchor_distance 2000 -- 5. The scores of all alignment pairs were recalibrated as follows. Phred-scaled insertion and deletion probabilities were set according to the empiric distribution obtained in Step 2.1. Deletions in the candidate set generated in Step 2.4 incur a constant (length-independent) alignment cost of 1. This upweights alignments that support deletions that have also been reported by CLEVER in Step 1 and are thus in line with read pair evidence. For each read pair the alignment pair with the highest probability is reported as primary alignment. All secondary alignments were discarded at this stage.
6. The number of LASER (primary) alignments supporting each deletion were extracted from the BAM files for further processing.
All deletions of CLEVER and LASER were merged. In this step, the whole GoNL population is considered simultaneously. We iterated over all putative deletions as reported by LASER sorted (decreasingly) by total support of primary alignments in the whole population.
1. For each individual, the set of CLEVER predictions was searched for deletions with a center distance of at most 100 and a length difference of at most 20 from the currently processed deletion. All CLEVER calls matching the current deletion are marked and ignored in subsequent iterations. This ensures that each CLEVER deletion is assigned to at most one LASER prediction, giving precedence to predictions common in the population.
2. If a deletion is supported by CLEVER with support of at least 5 and by at least one LASER (split) alignment, we report it to be present in the individual in question.
3. All deletions present in at least one individual are included in the final VCF file.

SOAPdenovo de novo assembly
We selected discordantly mapped and unmapped read pairs from each family and performed de novo assembly per family using SOAPdenovo2 13 , version r240 using kmer size 63 (http://soap.genomics.org.cn/soapdenovo.html). We then obtained family-specific contigs and assembled all family-specific contigs for the entire GoNL sample set together in a second round of de novo assembly. The resultant non-redundant contigs were compared against the GRCh37/hg19 genome reference using BLAT software. Alignments that consisted of two distantly mapped HSP segments of at least 150 base pairs that correspond to left and right part of a contig (with no more than 10 bp overlap at contig-level) were parsed by a custom script (available upon request) to define breakpoints of structural genome rearrangements.

Mobster
Non-reference mobile element insertions (MEIs) were identified using Mobster v0.1.6 with default parameters 14 . Analysis was run separately for each family. We only retained candidate MEIs supported by 5 or more reads.
The sensitivity for detection of MEIs by Mobster was separately tested using a set of 134 MEIs based on a combination of both single and double cluster predictions.

Identification of novel genomic segments
Two types of variants called in section 3.1.11 were analysed separately: i) full contigs that exceed 150bp in size and did not map anywhere with >99% homology when mapped by NCBI BLAST against the GRCh37/hg19 assembly and ii) parts of contigs larger than 150 bases that did not map to any genomic location. We realigned the individual-specific sets of discordant reads using these new segments as a reference sequence, in order to determine their presence/absence in the libraries of each individual.
We divided these new segments into two sets: i) potential contaminations (segments that appear in some of the libraries, but not in all libraries of the same individual) and ii) true new Table 3 GRCh38/hg38 genome reference or a decoy dataset hg38d1. We required 99% identity in the alignment between assembled segment and latest genome reference to discard a segment as unreported by GRCh38.

Construction of simple indel set
A data set with simple indels (1-20 base pairs) was constructed by merging four individual callsets obtained by running GATK HaplotypeCaller, Pindel, Mate-Clever and SOAPdenovo assembly (as described in previous sections).

Construction of complex indel set
Genomic regions showing a high density of polymorphisms (distance between adjacent polymorphisms below 30 basepairs) were tested for being complex events or alleles that potentially appeared as part of the single mutational event, but called as separate adjacent events. GATK HaplotypeCaller was rerun for these regions using the 'mergeVariantsViaLD' option and non-simple indels were retained. Simple events that originally represented these complex events in the simple indel dataset were then excluded from the simple indel dataset.

Variant validation
We Here, we describe an extension of these validations, which is based on the GoNL indel and SV call set released with this paper. For all validation rounds described below, we performed PCR amplification of the variant region or breakpoint junction regions (Supplementary Data   1). PCR primers were designed using primer3 software based on breakpoint coordinates.
Subsequently, we performed sequencing of PCR products using Sanger sequencing or MiSeq (2*250bp). Sanger sequencing reads were mapped to the human reference genome

Indels of 1--20bp
Besides the validation assays as described previously for the GoNL release 1 dataset 2 , we

21--100bp deletions
We have previously described validation for 96 deletions of 21-100bp 2 . We extended and reanalyzed these data and performed verification PCR assays and MiSeq sequencing for an additional 48 midsize deletions. Altogether, a total of 142 out of 144 21-100bp indels were confirmed.

Validation of 100bp+ deletions
We have previously described validation for 48 deletions larger than 100bp and observed a confirmatory rate of 46/48 2 . Here, we performed validation assays for an additional 51 100+ deletions by PCR and MiSeq on a single individual from the GoNL sample set. Two assays failed and two deletions could not be confirmed.

Validation of tandem duplications
A total of 48 tandem duplications were tested using PCR and MiSeq on a single individual.
We could confirm 41 out of 48 variants and failed to obtain a PCR product for 7 of them, resulting in a confirmation rate of 85.4%.

Validation of inversions
We identified a total of 76 inversions in the GoNL dataset. These were all subjected to experimental verification using PCR and MiSeq on diverse set of GoNL samples. Out of 76 tested inversion candidates, we confirmed 49, while the PCR assay failed for 27 variants, resulting in a confirmation rate of 64.5%.

Validation of mobile element insertions
We tested a 96 mobile element insertion (all Mobster calls) using PCR of one of the two breakpoint junctions. The assays were performed on a single individual of the GoNL dataset.
We successfully validated 92 mobile element insertions, while 4 assays failed, resulting in a confirmation rate of 95.8%.

Validation of interchromosomal rearrangements
All 42 interchromosomal rearrangements detected in the GoNL cohort were independently tested by PCR and MiSeq. Based on 42 breakpoint PCR assays, we could confirm 35 interchromosomal changes, resulting in a confirmation rate of 83.3%. Recent data showed that interchromosomal rearrangements may occur as a result of gene retrotransposition insertion polymorphisms (GRIPs) 16 . We confirmed at least 10 GRIP events in our set of interchromosomal rearrangements (Supplementary Data 2).

Validation of novel genomic segments
We tested 96 novel genomic segments, which were not present in the GRCh38 reference build. A total of 86 variants could be confirmed based on DNA from a single GoNL sample, while 10 assays did not result in a PCR product, resulting in a confirmation rate of 89.6%.

Validation of complex indels
We tested 10 complex indels and could confirm 8 of these, resulting in a confirmation rate of 80%.

Functional annotation of Variants
The consensus events were annotated using a custom annotation pipeline as well as

Effects of SVs on gene expression
Whole blood paired end RNA-seq data is currently being generated for the Genome of The

Variant Hotspots
To identify genomic regions where deletions frequently occur, 1 Mb genomic bins across the primary sequence were created and the number of unique deletion events per bin counted based on the consensus set. Regions with significantly more events were identified based on a z-score greater than 2 and compared to those previously reported 4 . Likewise to identify hotspots across variant types, the analysis was repeated to include deletion, duplication,

mobile element insertions and inversions. A cross variant hotspot was defined as a genomic
region with a z-score greater than 2 across more than 3 variant classes.
Candidate hotspots were subsequently filtered based on accessibility for sequencing as previously reported 1 .

Variant Novelty
To determine which SVs were known previously, we compared our call set to dbSNP, DGV, 1000 Genomes Data, CHM1 calls, and our previous public GoNL release 2,[21][22][23][24] . Details are given below. In addition to the datasets named above, the novelty of the MEI call set was also compared to three additional studies 15,25,26 .

Comparison to dbSNP
We compared simple and complex indel calls to dbSNP release 146 (released on 2015-11-04) by searching for exact matches of reference and alternative alleles between our and dbSNP polymorphisms.

Comparison to DGV
We compared our call set to release 2015-07-23 of the database of genomic variants (DGV, see footnote 1 ) 24 . For deletions, we considered all DGV entries where "varianttype" was equal to "CNV" and "variantsubtype" was "deletion" or "loss". For duplications, we considered DGV entries with type "CNV" and subtype "duplication", "gain", or "tandem duplication". For MEIs, we used DGV entries with type "CNV" and subtype "mobile element insertion". For inversions, entries with type "OTHER" and subtype "inversion" were used. For inversions, duplications, and deletions longer than 100bp, we counted an event as a match if its center was at most 2000bp away from the center of the respective DGV event, and the length difference was at most 500bp. For the (more frequent) deletions 21-100bp, we used more stringent criteria to avoid false positives and counted matches with a center distance of at most 100bp and a length difference of at most 20bp. For MEIs (where a reliable length estimate is often not available), we counted DGV events as a match if the insertion point was at a distance of at most 500bp.

Comparison to 1000 Genomes Data
We compared our call set to the integrated SV release of phase 3 of the 1000 Genomes Project, see footnote 2 . We compared all deletions >20bp to 1000G events where SVTYPE matched DEL, CNV, DEL_ALU, DEL_HERV, DEL_LINE1, or DEL_SVA. Duplications were compared to SVTYPES matching DUP or CNV. Inversions were compared to records with SVTYPE=INV and MEIs were compared to SVTYPES matching ALU, LINE1, or SVA. As described above for DGV, matching was done based on center distance and length difference, using the same cutoffs.

Comparison to previous GoNL release
Previously 2 , we have reported on a set of SNPs, simple indels, and deletions derived from the same NGS data. This call set is publicly available 3 as "release 5". We also compared the call set reported here to this previous one and only report variants as novel that are not included in the previous set. For indels we required an exact match of coordinates and called reference and alternative alleles. For matching deletions >20bp, we used the same criteria as used above for DGV and 1000 Genomes.

Comparison to CHM1
We obtained calls from CHM1 22 provided as three sets of events (insertions, deletions, and inversions) and compared them to our call set. For deletions and inversions, we used the same matching criteria as above. For insertions, the CHM1 data provides intervals (in a bed file). We counted a MEI as matching an insertion from CHM1 when the distance of our estimated insertion point was included in or at most 100bp away from one of the provided intervals. SNPs whose description specified a non-Caucasian cohort and SNPs where the risk allele was not specified. Removing duplicates yielded 7,457 SNPs, all of which were shown to be associated with disease in Caucasian individuals, which in the following we will refer to as GWAS SNPs. In the following, we only consider data of the parents to ensure that the sample is (approximately) independent. Both the SNPs and the deletions were phased using MVNCall 27 . SNP and deletion variants with too low minor allele frequency (< 4%) were discarded, following the workflow described in Maurano et al. 18 . This filtering step reduced the number of autosomal GWAS, non-GWAS SNPs and deletions to 3,198, 1,259,007 and 21,685, respectively.

Identification of significant (non)--GWAS SNP--deletion pairs
First, we tested for linkage disequilibrium (non-random association between alleles) in every GWAS SNP and deletion that are 1) on the same autosome, and 2) less than 1 million base pairs removed from each other (one CentiMorgan on average [Lodish et al., "Molecular Cell Biology"]), that is, associations are likely to be missed only in 1% of cases. The total number of GWAS SNP-deletion pairs that satisfied these constraints were 55,250. In order to assess association between each pair, we computed Pearson's coefficient of determination (R2) and applied Fisher's exact test (two-sided). We accounted for multiple testing by employing Benjamini-Hochberg's false discovery rate (FDR) control procedure with a q-value (expected  Pearson's coefficient of determination (R 2 ) and applying Fisher's (two-sided) exact test on the basis of the data of the parents. The SNP with the highest R 2 (of at least .8) and a p-value lower than 5% is regarded to be the most appropriate tag SNP for that particular deletion.

Tagging GoNL deletions
The results can be found in Supplementary Data 7.

Genotyping indels and SVs
SV genotyping as performed for all SVs in the consensus call set described in Section 3.3. In the following we describe the performed steps in detail for each SV type.

A maximum likelihood model for genotyping insertions and deletions
The following model has been adopted from [Marschall & Schönhuth, to appear], where it was described for midsize deletions. Here, we will make use of the model for deletions of arbitrary length, and also for insertions that can be detected via paired-end read and splitread data signals. The model described in the following has been implemented as part of the CLEVER Toolkit (that also contains MATE-CLEVER).
Let ! , ! , ! indicate the genotypes of an indel: ! denotes absence of the indel, ! denotes a heterozygous indel and ! a homozygous indel. In the following, we fix the region of interest, that is, the region that harbors the putative indel. Let A be a read alignment. Let R be a single read and R be the set of all reads that have an alignment with the region under consideration. For R in R, let A(R) be the alignment of R with the region we would like to genotype.
Let S be a subset of reads from R. In slight abuse of notation, we formally consider S as the event that precisely the alignments of reads from S in the region of interest are correct, while all others do not align in the region. Hence, is the corresponding probability. In the following, we consider a maximum likelihood (ML) setting, which in particular reflects that our prior belief in a genotype to be correct is the same for all types: This assumption is necessary for an efficient computation scheme, because a more flexible Bayesian approach would only yield a scheme that is exponential in |R|, the number of reads that align to the region of interest. If |R| exceeds 20, this becomes infeasible.
We are interested in maximizing (j=0,1,2) By considering probabilities P(S), we take alignment uncertainty due to the usual reasons like multiply mapped reads and alignment artifacts, which occur particularly often in indelaffected regions, into account. Let K:=|R| be the number of reads that align to the region to be genotyped. By Bayes' formula, equation (1) further implies that Note that we also assume that reads were generated independently of each other, which is a reasonable assumption. From equation (2), we conclude that where the second row results from expanding the third row. The last term, finally, can be computed in time linear in the number of reads R, which had been our goal. We use this formula for genotyping insertions and deletions of arbitrary length, that is, including insertion and deletion SVs.
It remains to compute reasonable probabilities ! |{ } and ! |{ } for read alignments A. While ! |{ } = ! because the read that underlies A does not stem from the region, hence has no influence on the respective genotype, computation of terms ! |{ } require further reasoning. One has to distinguish between two types of evidence for an indel given that A is correct: • Split-read evidence: A aligns with the region to be genotyped such that one read end stretches across the breakpoint(s) of the putative indel Let A be an alignment with a split/gap that agrees with D. By the above considerations, the read of A stems from a chromosomal copy that is affected by D. By Bayes' formula, and the formulas from above, In general, the assumption of no alignment uncertainty does not apply. In fact, split-read While these values may seem large, they are well supported by considerations on the uncertainty of (split-)alignments 9,30 .

Internal Segment--Based Evidence
Internal-segment-based evidence is provided by evaluating the empirical distribution on fragment length, which depends on the library from which the read stems. Here, we consider this distribution approximately Gaussian, which is a reasonable assumption for the vast majority of GoNL individuals. In general, all of the following considerations can be easily generalized to arbitrary, also non-Gaussian empirical distributions. Let D be the deletion to be genotyped and let C(D) be its centerpoint. Let R be the read of A where ! < < ! , that is, the alignment interval of A, which is the gap between ! , the alignment of the left end, and ! , the alignment of the right end, contains the counterpoint of the deletion D. Note that choosing centerpoints as anchors for the following considerations has advantages, as was explained in 9,30 . Let µ and σ be mean and standard deviation of the internal segment size distribution of the library R is from, which assumed to be Gaussian. So, internal segment length, as a random variable X, is distributed as the Gaussian distribution There are two cases: first, alignments A whose reads stem from a chromosomal copy that is not affected by D, and second, alignments B whose reads stem from a chromosomal copy that is affected by D. We obtain where the second case reflects that the alignment interval contains the deletion of length L.
We compute as appropriate densities for the cases of no variant and a homozygous variant. Let a normalization factor.
In analogy to considerations for the split-read case, we arrive at as an appropriate probability distribution for reads whose alignments span the breakpoints of deletions by their internal segments.
The procedure described above is implemented as part of MATE-CLEVER 9 , which can also use prior information in form of the Mendelian laws, if the input consists of multiple, ancestryrelated genomes. In order to genotype, MATE-CLEVER collects all (split-read and regular) alignments in a region, and uses them for carrying out the computations described above.

Mid--size deletions (20--10000bp)
For the consensus deletion calls in this length range, the respective regions (deletions +/-750bp) were extracted from the BAM files for all samples, converted to FASTQ format, and aligned using LASER 9 , version v2.0-rc1-20-gfea4980, with the following command line: The resulting VCFs for each family are merged using vcf-merge, which is part of vcftools 31 , to create a VCF for the full panel.
Phasing requires well-calibrated genotype likelihoods. We therefore conducted a simulation study to recalibrate GLs. We created three kinds of simulated diploid genomes, one that does not contain any midsize deletions (i.e. homozygous for the reference allele at all midsize deletion loci), one that is heterozygous for all midsize deletions (each one independently put at random onto one of the two haplotypes), and one that is homozygous for all midsize deletion alleles. We additionally added SNPs and indels to these three simulated genomes by randomly drawing SNP and indel alleles from the GoNL panel according to their AFs. We repeated this procedure 10 times to create a total of 30 simulated diploid genomes. For all these genomes, we simulated Illumina reads using the simseq program 32 , matching the coverage found in the GoNL. The resulting reads were aligned to the reference genome using LASER with the same parameters as described above.
Subsequent genotyping using the CLEVER Toolkit was also performed as described above.
For each deletion, we created a confusion matrix, counting the number of times we observed genotyped !"#$ (predicted by our genotype caller) when the true genotype was !"#$ , for all nine possible combinations of genotypes. Based on these counts, we estimated the empirical conditional probabilities ( !"#$ | !"#$ ), based on the assumption that ( !"#$ | !"#$ ) is proportional to ( !"#$ | !"#$ ). In the absence of observations (i.e. when the number of observed genotypes !"#$ is zero), we assume uniform conditional probabilities (corresponding to a flat prior). Using these deletion specific tables of conditional probabilities, we recalibrated all genotype likelihoods by computing posterior likelihoods.

Long deletions (10kbp and longer)
In principle, the genotyping of long deletions proceeded like the genotyping of mid-size deletions described above, with the following slight adaptations: • The analysis was done at a later point in time using the (then current) version v2.0-rc3-53 g7b13d29 of the CLEVER Toolkit. No major changes happened compared to the versions used for mid-size deletions, just slight fixes/usability improvements.
• Due to the size of the events, we did not extract the full-length deletions from the original BAM files, but used a +/-1000bp window around start position, center position, and end position of each deletion, which is sufficient for genotyping.

Inversions
To genotype the inversions in our consensus set, we used DELLY 33 , version 0.5.9. Since precise breakpoint coordinates were not available for all consensus inversions, we did not run delly in "genotyping mode", but used the "discovery mode" as follows:

Translocations
Genotyping of translocation was done analogously to genotyping inversions (see previous section).

Duplications
Genotyping of duplications proceeded mostly like for inversions as described above, with the 2. We matched consensus duplications to DELLY prediction, again using a maximum breakpoint distance of 500 and resolved ambiguities as follows: we prefer calls with a FILTER value of "PASS" over those with a "LowQual", we prefer PRECISE over IMPRECISE, we prefer calls with a lower distance to the original consensus duplication, and finally we prefer calls with a lower length difference.

Mobile Element Insertions
Genotyping of mobile element insertions (MEIs) was done via a logistic regression model,

Phasing SVs using MVNCall
We used the same haplotype scaffolds as previously described 2 to phase SVs onto the already phased sets of SNPs and indels. The scaffold contains sites present on Omni2.5M chips. Refer to the supplement of Francioli et al 2 , Section 12, for details on how it was created. For phasing, we used MVNcall 27 version 1.1. Apart from PED files (containing the family/trio structure) and scaffolds, we used the genotype likelihoods (GLs) reported by the CLEVER Toolkit (for deletions) 10 , by DELLY (for inversions, duplications, and translocations) 33 , by Mobster 14 (for mobile element insertions), and by GATK's HaplotypeCaller (for complex indels) as input. Before phasing, the GLs were regularized so as to avoid too low probabilities as follows. We added a value of 10 5 to each of the three probabilities for genotypes (RR, AR, and AA) and renormalized the sum to one. The so transformed GLs are passed to MVNcall 27 for phasing using the following command line: "mvncall --sample-file samples.list --ped-file samples.ped --k 100 --iterations 50 --add-gls --int 0 300000000 -scaffold-file scaffolds.haps --glfs input.vcf --o output.vcf". We ran MVNcall in this way for each chromosome and for each SV type separately.
3.6.2 Imputing SVs in independent Dutch individuals using IMPUTE2 In order to assess the quality of our reference panel with respect to its structural variants (see 37 for an encompassing study on the quality of the panel in terms of SNPs), we ran IMPUTE2 38 on two independent Dutch individuals and analyzed the resulting genotypes. We aimed at a scenario where one could compare the imputed genotypes with true genotypes.
This differs from the evaluation in Deelen et al 2014 37 who focused on classical imputation metrics 39 as a basis for evaluation, which requires sufficiently many individuals, but does not use true genotypes. Since the application of imputation metrics for SVs has not been discussed before in the literature, we opted for a scenario where one could juxtapose imputed genotypes with (highly likely) true genotypes.
For genotyping structural variants in the two individuals, we analyzed the available wholegenome sequencing data of the individuals, which were Illumina HiSeq paired-end reads at 17x and 19.5x coverage respectively, aligned against reference genome GRCh37 using Bowtie2 v2.2.4 40 , followed by marking duplicated reads and indel realignment. For genotyping SVs, we used the same tools with the same setting as for the GoNL panel, see Section 3.5; that is, we used MATE-CLEVER 9 for genotyping deletions, DELLY 33 for genotyping inversions, duplications and translocations, Mobster 14 for genotyping mobile element insertions, and GATK's HaplotypeCaller for genotyping complex indels. We ran all these tools in genotyping mode -and not in discovery mode -on the structural variants that belong to the panel. This yielded a set of read-based genotyped SVs that is in 1-to-1correspondence with the SVs one can impute, in particular also in terms of breakpoint • lower and upper are integers to indicate the genomic region to be studied (we use the full chromosome).
This resulted in imputation-based genotype likelihoods (GL's) for each SV. We then filtered these GL's for SV's that correspond to SV's where gold standard genotypes, as described above, were available.
Comparing genotypes / allele dosages. Let V be a placeholder for single SV's, and let V be all SV's considered (stratified by type or panel-based allele frequency, for example). For each V, we compare the imputed (SNP-and panel-based) genotype likelihood !,! , !,! , !,! which correspond to the imputation-based probabilities that the V is not present (0), heterozygous (1) or homozygous (2), with the directly read-based, carefully filtered "gold standard genotype" likelihoods !,! , !,! , !,! which reflect the ground truth probabilities that V is absent (0), heterozygous (1) or homozygous (2); note that one of the !,! is at least 0.999 as per the filtering criteria (0.85 for MEI's). Let V be a placeholder for single SV's. We then compute the probability ! that the imputed genotypes agree with the ground truth as


extfieldsofstudy: Biology Medicine

s2fieldsofstudy: Biology Chemistry

text:
Distinct Activation Mechanisms of NF-κB Regulator Inhibitor of NF-κB Kinase (IKK) by Isoforms of the Cell Death Regulator Cellular FLICE-like Inhibitory Protein (cFLIP)*

The viral FLICE-like inhibitory protein (FLIP) protein from Kaposi sarcoma-associated herpesvirus activates the NF-κB pathway by forming a stable complex with a central region (amino acids 150–272) of the inhibitor of NF-κB kinase (IKK) γ subunits, thereby activating IKK. Cellular FLIP (cFLIP) forms are also known to activate the NF-κB pathway via IKK activation. Here we demonstrate that cFLIPL, cFLIPS, and their proteolytic product p22-FLIP all require the C-terminal region of NEMO/IKKγ (amino acids 272–419) and its ubiquitin binding function for activation of the IKK kinase (or kinase complex), but none form a stable complex with IKKγ. Our results further reveal that cFLIPL requires the linear ubiquitin chain assembly complex and the kinase TAK1 for activation of the IKK kinase. Similarly, cFLIPS and p22-FLIP also require TAK1 but do not require LUBAC. In contrast, these isoforms are both components of complexes that incorporate Fas-associated death domain and RIP1, which appear essential for kinase activation. This conservation of IKK activation among the cFLIP family using different mechanisms suggests that the mechanism plays a critical role in their function.

The first members of the FLICE 2 -like inhibitory protein (FLIP) family were identified as viral genome products (vFLIPs) following a bioinformatics search conducted to recognize death effector domain-containing proteins (1). A consensus sequence from the viral FLIPs was then used to screen human expressed sequence tags, and several human homologs were identified, now known collectively as cellular FLIPs (cFLIPs) (for a review, see Ref. 2). Cellular FLIPs are encoded by a single gene; in humans there are three main splice variants, long isoform (cFLIP L ; 55 kDa), short isoform (cFLIP S ; 26 kDa), and cFLIP Raji (cFLIP R ; 24 kDa) first detected in the human Burkitt lymphoma B cell line Raji (3). Mice do not express cFLIP S and contain only two isoforms cFLIP L and cFLIP R (4).
All cFLIP isoforms share an identical N terminus of 202 amino acids, which includes two death effector domains, but have different C termini. The structure of cFLIP S is similar to the vFLIPs, but it contains a 19-amino acid C terminus that is responsible for its ubiquitination and proteasomal degradation (5). In cFLIP L , the tandem death effector domains are followed by a caspase-like domain. In addition, two N-terminal cleavage products have been detected, p43-FLIP and p22-FLIP, both thought to be involved in stimulating NF-B signaling (6,7). The best described function of cellular FLIPs is in regulation of death receptor (DR) signaling. Procaspase-8 activation, which is pivotal for DR-induced apoptosis, occurs through an "induced proximity" mechanism where DR-mediated homodimerization of procaspase-8 leads to its self-cleavage and activation. Like procaspase-8, cFLIPs are recruited to the DR complex by FADD and can thus directly inhibit procaspase-8 activation (8). However, cFLIP L in particular also contributes to DR-induced NF-B signaling (9 -11), particularly when expressed at a moderate level (12).
The canonical NF-B pathway is activated by inhibitor of NF-B kinase (IKK), which phosphorylates IB, leading to its degradation. This releases NF-B transcription factors, which then activate gene expression. vFLIP is a latently expressed gene from Kaposi sarcoma-associated herpesvirus (KSHV), which activates NF-B signaling (13). NF-B activation by vFLIP leads to the characteristic phenotype and cytokine secretion of Kaposi sarcoma cells (14) and survival of KSHV-transformed B cells (15)(16)(17). vFLIP activates IKK (18) by forming a stable complex with the regulatory subunit IKK␥ (19). We have solved the crystal structure of two vFLIP molecules bound to a parallel ␣-helical coiled coil of two IKK␥ molecules and proposed a mechanism whereby this binding leads to activation of the IKK␣/␤ kinase subunits (20).
In this study, we wished to investigate how cFLIPs activate IKK and in particular whether they form a stable complex with IKK in a similar manner to vFLIP. We therefore engineered cells to express the various cFLIP isoforms and measured their ability to activate IKK in the absence of any additional signal. In contrast to previous reports, we did not find any direct interaction with IKK␥. The ubiquitin binding function of IKK␥ was required by all the cFLIP variants. In the case of cFLIP L , the linear ubiquitination complex LUBAC was necessary, presumably to generate a ubiquitinated substrate to interact with IKK␥. In contrast, cFLIP S and p22-FLIP formed a complex including FADD and RIP1 that led to RIP1 ubiquitination and activation of IKK. These data are relevant to a variety of human cancers as we have demonstrated that human tumor cells express a number of cFLIP variants at levels similar to those in our engineered cells.
Lentivectors and Cell Transductions-cDNA of the noncleavable cFLIP L was a kind gift from Prof. Inna Lavrik. cDNAs of wild-type or mutant FLIP variants, Tax, and N-terminally HA-tagged FADD were cloned into pDual (see Fig. 1B). Vesicular stomatitis virus G glycoprotein pseudotyped lentiviral particles were produced by co-transfection of HEK293T cells with pDual, pCMV8.91 (HIV gag-pol), and pMDG (vesicular stomatitis virus G glycoprotein). Supernatants were passed through 0.45-m filters and concentrated 100-fold by ultracentrifugation at 100,000 ϫ g at 4°C for 2 h in a Sorvall centrifuge (Beck-man Coulter). Viral pellets were resuspended in cold whole medium, incubated on ice for 1-2 h, aliquoted, and stored at Ϫ80°C. Titers were measured on HEK293T cells by fluorescence-activated cell sorting (FACS) analysis and TaqMan quantitative PCR (Applied Biosystems, Warrington, UK). For experiments, cells were transduced with test cFLIP, vFLIP, and Tax lentivectors at identical multiplicities of infection within experiments and in the range of 1-20 or 50 -500 for 293/293T or 70Z/3, respectively.
Generation of Stable Knock-down Cell Lines-shRNA-expressing lentivectors targeting HOIP, TAK1, Atg3, caspase-8, RIP1, and FADD were produced using pGIPZ lentivector plasmids (University College London Open Biosystems shRNA library, London, UK) that contain a puromycin resistance gene and a GFP marker. For silencing of HOIL-1 and SHARPIN, we designed customized shRNA-encoding DNA duplexes using the Clontech online RNAi designer tool and cloned them into the pSIREN Hygro lentivector backbone (modified from pSI-REN RetroQ; Clontech). The transduced cells were selected with hygromycin B (200 g/ml). The sequences targeted for each gene have been listed in Table 1.
Immunoblotting, Precipitation, and Kinase Assays-Cell extracts were prepared with lysis buffer (1% Triton X-100 in 20 mM Tris-HCl, pH7.5, 150 mM NaCl, and 5% glycerol supplemented with phosphatase and protease inhibitors). Lysates (20 -40 g/lane) were separated on 8 -12% SDS-polyacrylamide gels and transferred to PVDF membranes (GE Healthcare), which were then blocked with 5% skimmed milk in TBS-T (0.1% Tween 20 in 50 mM Tris-HCl, pH 7.6, and 150 mM NaCl) for 1 h. Membranes were incubated overnight at 4°C with suitable dilutions of primary antibodies in TBS-T buffer containing 5% BSA and 0.1% sodium azide. Bound antibodies were detected with appropriate peroxidase-conjugated secondary antibodies and visualized by the ECL system (LumiGLO, Cell Signaling Technology). For immunoprecipitation, lysates were first precleared by incubating with Protein A-or G-conjugated beads (Sigma) for 30 min at 4°C. The precleared extracts (300 -500 g) were then incubated with 1 g of primary antibody or its isotype-matched IgG control and 20 l of Sepharose beads for 2-4 h at 4°C. The beads were washed twice with high salt wash buffer (25 mM Tris-HCl, pH7.6, 500 mM NaCl, 1 mM DTT, and 5% glycerol) and two additional times with the lysis buffer and then were resuspended in 20 l of 2ϫ Laemmli buffer. IKK kinase assays were performed on immune complexes as described previously (19). The reaction mixtures were separated by 10% SDS-PAGE, transferred to a PVDF membrane, and exposed to a phosphor screen for 15-30 min. The radiolabeled phosphoproteins were then visualized using a Typhoon phosphorimaging system (GE Healthcare) and analyzed by ImageQuant TL Plus 7.0 software (GE Healthcare).
In Vitro IB␣ Phosphorylation Assay-Resting HEK293T cells were lysed in HTX lysis buffer (10 mM Hepes, pH 7.4, 10 mM MgCl 2 , 1 mM MnCl 2 , 0.1 mM EGTA, and 0.5% Triton X-100) supplemented with complete protease inhibitor mixture (Roche Applied Science), and the S100 fractions were prepared by sequential centrifugation at 1,000 and 100,000 ϫ g for 5 and 45 min, respectively. The cleared lysates were then incubated with recombinant FLIP and an ATP-generating system (10ϫ stock: 10 mM ATP, 20 mM Hepes, pH 7.2, 10 mM MgCl 2 , 350 M creatinine phosphate, and 500 g/ml creatinine kinase) for 10 min at 37°C. The reactions were terminated by adding 5ϫ Laemmli buffer and boiling at 95°C. The recombinant FLIP-induced IKK activation rate was detected by immunoblotting for pIB␣ and pIKK␣/␤.
Purification and Expression of Recombinant vFLIP, p22-FLIP, and GB1-p22-FLIP-Recombinant vFLIP was overexpressed and purified as described in Bagnéris et al. (20). For the production of recombinant p22-FLIP (and GB1-p22-FLIP), cFLIP S was first PCR-amplified from a human thymus cDNA library (Clontech) using the primers 5Ј-TTGCTAGCATGTCTGCT-GAAGTC and -TTCTCGAGTCACATGGAACAATTTC containing NheI and XhoI restriction sites, respectively (underlined). Following digestion, the product was then ligated into the pETM442 vector (described in Bagnéris et al. (20)) to provide an N-terminal His 6 -NusA tobacco etch virus proteasecleavable solubility tag to aid expression and purification. The resulting pETM442-cFLIP S vector was then used as the template for production of pETM442-p22-FLIP, which was generated by the introduction of a stop codon at Asp-196 using the primers TCCAAAAGAGTCTCAAGTAGCCTTCAAATAA-CTTCAGGAT and ATCCTGAAGTTATTTGAAGGCTAC-TTGAGACTCTTTTGGG. GB1-p22-FLIP was produced by subcloning p22-FLIP into the EcoRI and XhoI restriction sites of a modified pET15b vector (a kind gift from Prof. Paul Driscoll, National Institute for Medical Research) to produce pET15-GB1-p22-FLIP (containing a His 6 -GB1 solubility tag). The original pET15b vector was modified by inclusion of the coding sequence for the 56-amino acid immunoglobulin domain of Escherichia coli protein G (GB1) 5Ј to the multiple cloning site, the addition of an EcoRI restriction site directly 3Ј to GB1, and removal of the multiple cloning site EcoRI site. This enabled subcloning of p22-FLIP 3Ј to GB1 into the new EcoRI and existing XhoI sites. pETM442-p22-FLIP and pET15-GB1-p22-FLIP were then transformed into BL21(DE3) Star TM cells (Invitrogen). The cultures were then inoculated into LB medium containing ampicillin (100 g/ml) at a ratio of 1:100 and grown to an A 600 of between 0.6 and 0.8. Cultures were subsequently induced with 1 mM isopropyl 1-thio-␤-D-galactopyranoside, and the temperature was reduced to 16°C overnight. Cells were harvested by centrifugation, and the pellets were either frozen at Ϫ80°C or resuspended in a buffer comprising 25 mM Tris-HCl, pH 8.5, and 200 mM NaCl (buffer A) supplemented with an EDTA-free protease inhibitor complex tablet (Roche Applied Science) and DNase I (10 g/ml). The purification protocols used for p22-FLIP and GB1-p22-FLIP were nearly identical. Resuspended pellets from 8-liter cultures of both proteins were sonicated on ice, and the lysates were clarified by centrifugation (46,000 ϫ g for 1 h). Supernatants were subsequently filtered through a 0.45-m filter prior to loading on HisTrap FF columns (GE Healthcare) pre-equilibrated with buffer A. These were then washed with buffer B (buffer A and 20 mM imidazole) and eluted with buffer C (buffer A and 500 mM imidazole). For p22-FLIP, the His 6 -NusA tag originating from the pETM442-p22-FLIP construct was removed by incubating His 6 -NusA-p22-FLIP with tobacco etch virus protease overnight while dialyzing in buffer A (to remove the imidazole), and the solution was reapplied to a 5-ml His-Trap FF column pre-equilibrated in buffer A. p22-FLIP was eluted from the column using buffer A and 40 mM imidazole. Fractions having greater than 95% purity as judged by SDS-PAGE were subsequently pooled and frozen at Ϫ80°C. Following elution from the 5-ml HisTrap column, GB1-p22-FLIP appeared to be over 95% pure. To improve stability, however, GB1-p22-FLIP was subsequently dialyzed into a buffer containing 5 mM L-arginine, 300 mM NaCl, 5 mM tris(2-carboxyethyl)phosphine, and 25 mM imidazole prior to concentration using a 6-ml, 10-kDa-molecular mass cutoff Vivaspin concentrator (Millipore) and storage at Ϫ80°C.
NF-B Luciferase Reporter Assays-For some reporter assays, cells were transfected with 300 ng of NF-B firefly luciferase plasmid (pGL immunoglobin light chain (IgK) or pGL H2- Dk), 100 ng of pRL.TK as an internal control, and 500 ng of pcDNA3 vectors encoding NF-B activators. 24 h later, luciferase activity was measured using a Dual Glo Luciferase kit (Promega). Firefly luciferase activity levels were normalized to that of Renilla. Alternatively, HEK293/293T and 70Z/3 cells were transduced with NF-B firefly luciferase lentivectors at multiplicities of infection of 10 and 500, respectively. NF-B-induced luminescence was detected using a Varioskan Flash multimode reader (Thermo Scientific).
Statistical Analysis-Statistical differences between two groups were analyzed by the two-tailed unpaired Student's t test using the GraphPad Prism v4.03 software package (GraphPad Software, La Jolla, CA). The calculated p values are given in the figures. All experiments were performed in triplicates on at least three independent occasions. The quantitative data are presented as mean Ϯ S.D.

Results
cFLIP Isoforms Require the Ubiquitin Binding Function of IKK␥-To determine which functions of IKK␥ were necessary for IKK activation by cFLIP isoforms, we used the IKK␥-null cell line 1.3E2 (21) harboring an NF-B-responsive luciferase gene. These cells were reconstituted with a truncated IKK␥ lacking the ubiquitin binding domain (⌬271), a point mutant that does not bind linear ubiquitin (F312A (22)), or a mutant defective in vFLIP binding (F238R/D242R (23)) (Fig. 1A). To do this, we used lentivectors expressing wild-type or mutant IKK␥ together with an mCherry fluorescent protein (Fig. 1B). Cells were transduced so that about 50% were mCherry-positive. We then isolated cell clones, expanded those that were mCherrypositive, and performed immunoblotting to establish that the transduced cells expressed IKK␥ (Fig. 1C). Cells were then infected with a second lentivector expressing cFLIP variants, vFLIP, or Tax from human T lymphotropic virus, type 1, another known IKK activator that binds directly to IKK␥ (24), together with enhanced GFP. After 48 h, enhanced GFP was monitored, and cells with transduction rates over 80% were used to measure luciferase activity. Fig. 1E shows that both vFLIP and Tax could activate IKK in the absence of its ubiquitin binding function, unlike the control Toll-like receptor agonist lipopolysaccharide (Fig. 1D). Tax has been reported to interact with two regions of IKK␥ at both the N and C termini (Ref. 23 and shown in Fig. 1A). However, our data demonstrate that only the N-terminal interaction is necessary for IKK activation. It was also reported previously that vFLIP can activate IKK reconstituted with IKK␥ truncated at amino acid 251 (25). In contrast, Fig. 1F shows that cFLIP L , cFLIP S , and their proteolytic product p22-FLIP all require the ubiquitin binding domain of IKK␥ and its ubiquitin binding function to activate IKK but do not require the residues that interact with vFLIP.
To confirm these data, we examined the phosphorylation of IB␣ in the same cell lines. Fig. 2 shows that vFLIP and Tax could stimulate phosphorylation of IB␣ when IKK␥ lacked the ubiquitin binding domain in contrast to all the cFLIP isoforms. We also observed that Tax could stimulate phosphorylation of IB␣ in cells lacking IKK␥ (Fig. 2), and therefore Tax gave a higher stimulation of phosphorylation of IB␣ when IKK␥ lacked the ubiquitin binding domain than did vFLIP. This was observed previously, although the mechanism remains unexplained (21).
cFLIP Isoforms Generate an Active IKK without Stable Interaction with IKK␥-To examine the mechanism of activation of IKK by cFLIP isoforms, we transfected HEK293T cells with expression plasmids for cFLIP S , cFLIP L , a non-cleavable version of cFLIP L (cFLIP L D196E/D376N), and p22-FLIP as well as vFLIP as a positive control for IKK direct interaction. We then immunoprecipitated IKK using an anti-IKK␥ antibody and measured its ability to phosphorylate recombinant IB␣. Fig.  3A shows that in each case the cFLIP isoforms generated an activated IKK. The level of activation was comparable or greater than that generated by vFLIP although less than the transient activation observed following TNF-␣ treatment (Fig.  3A). In contrast to a previous report, we observed a similar level of activation with non-cleavable cFLIP L as with wild type (6). As another report had used an NF-B luciferase reporter with an H2-Dk response element to measure the activity of cFLIPs (7), we also tested this NF-B element. Fig. 4 shows that the noncleavable cFLIP L was no different in activity from wild type with either reporter. Note that cFLIP S generally gave a relatively lower level of stimulation in HEK293 and 293T cells, which we cannot explain. We then examined whether the cFLIP isoforms were found associated with the activated kinase and found no evidence for the presence of cFLIP isoforms in the IKK immunoprecipitates, although vFLIP was clearly present (Fig. 3B). We also demonstrated that recombinant p22-FLIP does not form a stable complex with recombinant IKK␥ (data not shown) unlike vFLIP (20). These data are in contrast to the previous reports that showed p22 and p43 fragments stably associated with IKK (7,12). Clearly, differences in the experimental conditions can explain these results; however, our data demonstrate that the cFLIP isoforms can generate an active IKK without remaining physically associated with the IKK complex.
All cFLIP Forms Require TAK1, and cFLIP L Requires LUBAC to Activate IKK-Given that the cFLIP isoforms require the ubiquitin binding domain of IKK␥ (Figs. 1 and 2), we sought to establish which of the ubiquitination pathways have a key role and whether the cFLIP variants have different requirements for activation. Linear ubiquitin chain binding to IKK␥ is crucial for IKK activation by TNF-␣ (22), and these are generated by the trimeric LUBAC composed of HOIL-1, HOIP, and SHARPIN (26). We therefore generated HEK293 cells that were stably transduced with lentivectors containing short hairpin RNAs targeting HOIL-1, HOIP, and SHARPIN (Fig. 5A). As reported previously, LUBAC was required for full IKK activation by TNF-␣ (26) and was dispensable for activation by vFLIP (27) and Tax (Fig. 5, B and C). A clear blockade of cFLIP L activation of IKK was observed in the LUBAC knockdown cells (Fig. 5C). However, cFLIP S and p22-FLIP were unaffected by LUBAC knockdown (Fig. 5C). We also overexpressed in HEK293 cells the deubiquitinase CYLD, which can remove ubiquitin chains with linear as well as KLys-63 linkages (28) (Fig. 5D) and which inhibited full IKK activation by TNF-␣ (Fig. 5E). In this case, cFLIP L was once again inhibited, whereas cFLIP S and p22-FLIP were unaffected (Fig. 5F). The kinase TAK1 has been identified as a critical activator of IKK, phosphorylating the activation loops of the kinase subunits in response to a variety of signals (29,30), and the kinase MEKK3 has also been reported to fulfil this role in some circumstances (31). We therefore also produced cells with stably knocked down TAK1 or MEKK3 (Fig. 5G). TAK1 was clearly necessary for activation of IKK by all the cFLIP variants, although both kinases were dispensable for vFLIP and Tax FIGURE 5. Only cFLIP L is dependent on LUBAC to induce NF-B, whereas all cFLIP isoforms require TAK1. A, HEK293 cells were transduced with lentivectors encoding either scrambled control shRNA (scr) or shRNAs targeting HOIL-1, HOIP, and SHARPIN sequentially. Generation of LUBAC knockdown (KD) cells was confirmed by immunoblotting against each component of the complex and the housekeeping protein GAPDH. ns, nonspecific. B, TNF-␣-induced NF-B activation which was measured 6 h after stimulation at 10 ng/ml concentration is inhibited in LUBAC knockdown cells. C, scrambled and LUBAC knockdown HEK293 cells were co-transfected with an NF-B firefly luciferase reporter construct (300 ng/well) and a Renilla luciferase reporter vector (normalization control, 100 ng/well) with empty or transactivator-expressing pcDNA3 vectors (500 ng/well). The luciferase reporter assay was performed 24 h post-transfection. D, immunoblot showing the expression of ectopically expressed CYLD in HEK293 cells. Cells transduced with a GFP-expressing lentivector were used as controls. Effects of CYLD overexpression on NF-B activation levels induced by TNF-␣ (E) or cFLIPs, vFLIP, and Tax (F) were measured as described above. G, shRNA-mediated stable silencing of TAK1 and MEKK3 in HEK293T cells was validated by immunoblotting. H, these cells were further transduced with NF-B luciferase lentivector to develop NF-B sensor cell lines. Luciferase reporter assays were then performed 6 h following stimulation with TNF-␣ (10 ng/ml) or 48 h following transduction with control (enhanced GFP) or test lentivector (I). Results shown are mean Ϯ S.D. (error bars) of a representative experiment of three, performed in triplicates. The two-tailed Student's t test was used to determine the statistical differences between two groups. RLU, relative luciferase units; LLV, luciferase lentivector.
activation. The lack of a requirement for TAK1 in vFLIP activation has been reported (27), although its role in Tax activation of IKK has been controversial (32, 33). The requirement for TAK1 in the cFLIP isoform-mediated pathways suggests that the active kinase identified in Fig. 3 has been phosphorylated and activated by TAK1. The lack of requirement for ubiquitin signaling by cFLIP S and p22-FLIP (Fig. 5, A and D) required further investigation as the ubiquitin binding domain of IKK␥ FIGURE 6. p22-FLIP and cFLIP S require a FADD⅐RIP1 complex to induce IKK. A, stable gene knockdown of Atg3, caspase-8 (Casp-8), FADD, and RIP1 in HEK293T cells was validated by immunoblotting. From five independent shRNAs used to silence each protein, two that yielded the most efficient knockdowns were chosen for the experiments. Targeted sequences for each gene are listed as in Table 1. scr, scrambled. B, control and knockdown cells were transduced with NF-B luciferase lentivector to develop stable NF-B reporter cell lines. TNF-␣-, cFLIPs-, vFLIP-, and Tax-induced NF-⌲B activation levels were then evaluated as described in Fig. 3, B and D. Values shown are the mean Ϯ S.D. (error bars) of the relative luciferase activity from five independent experiments. C, a stretch of hydrophobic amino acid (AA) residues on death effector domain 2 of cFLIP that is critical for interaction with FADD is highly conserved in the KSHV (KS) vFLIP. Sequence alignment was performed by PRALINE software. D, mutating the critical amino acids of this hydrophobic patch (Phe-114, Leu-115 in cFLIPs and Phe-115, Leu-116 in vFLIP) to alanine residues abrogated NF-B activation by cFLIP variants but not vFLIP. NF-B luciferase reporter assays were carried out in HEK293T cells as described in Fig. 5C. E, co-immunoprecipitation assays showing that F114A/L115A mutants of cFLIP variants are unable to associate with the FADD⅐RIP1 complex. HEK293T cells (seeded in 6-well plates) were transfected with 1 g of HA-FADD vector together with 1 g of pcDNA3 constructs expressing cFLIP variants or vFLIP. 48 h later, cell lysates were extracted and immunoprecipitated (IP) with 0.5 g of anti-HA or an isotype-matched control antibody (C). Cell lysates and the immunoprecipitates were then analyzed by Western blotting with the indicated antibodies. * indicates the position of p43-FLIP bands. and its ubiquitin binding function were clearly required (Figs. 1  and 2). cFLIP S and p22-FLIP Activate IKK via a FADD⅐RIP1 Complex-To further probe the mechanism of action of cFLIP S and p22-FLIP, we generated HEK293T cells that were knocked down for other molecules reported to be present in complexes with cFLIP isoforms or vFLIP (Fig. 6A). Atg3 has been reported to interact directly with vFLIP and cFLIPs (34), FADD interacts with cFLIP isoforms (for a review, see Ref. 2), and RIP1 interacts with FADD (35). Fig. 6B shows that FADD and RIP1 were both required for cFLIP S and p22-FLIP activation of IKK, whereas Atg3 and caspase-8 were dispensable. A hydrophobic stretch of amino acids in cFLIP has been identified as critical for its interaction with FADD (4), and a similar sequence can be found in vFLIP (Fig. 6C). We mutated critical amino acids in this hydrophobic stretch, which resulted in blocking IKK activation by cFLIP S , p22-FLIP, and cFLIP L but not vFLIP (Fig. 6D). We then demonstrated that this mutation prevented interaction of the cFLIP isoforms with FADD, whereas FADD and RIP1 were present in these cells as a preassembled complex (Fig. 6E). These data provide an explanation for the lack of inhibition of cFLIP S and p22-FLIP by CYLD as ubiquitinated RIP1 bound to IKK␥ has been reported to be protected from deubiquitination (36). Although mutation of the putative FADD binding domain in cFLIP L prevented activation of IKK (Fig. 6D), cFLIP L did not require FADD or RIP1 for its action (Fig. 6B), suggesting that this region interacts with a different effector to trigger LUBAC activity. This effector may be TRAF2, which has been reported to interact with the N-terminal region of cFLIP L and to be required for NF-B activation (6). Our finding that RIP1 was not required for cFLIP L activation differs from the report of its role in T cell receptor-triggered NF-B activation where RIP1 was required (11).
To further examine the role of RIP1 in IKK activation by cFLIP S and p22-FLIP, we then examined ubiquitinated RIP1 binding to IKK␥. Fig. 7 shows that both cFLIP S and p22-FLIP induced ubiquitinated RIP1 binding to IKK␥, whereas cFLIP L and vFLIP did not.
Recombinant vFLIP Can Directly Activate IKK When Added to Cell Lysate-As vFLIP and cFLIP isoforms had very different requirements for ubiquitination and TAK1 in their activation of IKK, we then produced recombinant vFLIP and p22-FLIP proteins and added them to HEK293T cell lysate. Fig. 8A shows that recombinant vFLIP induced phosphorylation of the kinase subunits of IKK and that IKK became activated in the lysate as IB␣ was phosphorylated. This occurred within 10 min of incubation with the lysate; in contrast, recombinant p22-FLIP could not activate IKK under these conditions (Fig. 8A). Similar direct activation of the IKK kinase in cell lysates supplemented with recombinant Tax has been reported (37).
The results from this study allowed us to formulate a model of how vFLIP and cFLIP isoforms activate IKK (Fig.  8B). We propose that vFLIP binds directly to IKK␥ and causes activation of the kinase subunits by stimulating autophosphorylation independently of upstream ubiquitination or kinases. In contrast, cFLIP S and p22-FLIP are found in a complex with FADD and RIP1. As cFLIP S and p22-FLIP could not be detected bound to IKK (Fig. 3), we suggest that ubiquitinated RIP1 is generated by and then released from this complex. Ubiquitinated RIP1 then binds directly to the ubiquitin binding domain of IKK␥, recruits TAK1, and activates IKK. However, cFLIP L requires LUBAC to generate a different, as yet unidentified, ubiquitinated substrate, which then activates IKK.
In many of these experiments, we had constitutively expressed cFLIP isoforms using lentivectors. To test whether this mimicked physiological expression levels of cFLIP, we examined the level of endogenous cFLIP isoform expression in an extensive panel of human tumor cells. Fig. 9A shows that all the cells we tested expressed cFLIP L and at least one shorter cFLIP variant; we did not detect p22-FLIP in any cell. For these experiments, we changed the cFLIP antibody to the monoclonal G11, which is considerably more sensitive than the monoclonal NF6 for the detection of short human cFLIP variants. Fig. 9B shows that the level of cFLIP S and cFLIP L achieved by lentivector transduction of HEK293T cells was a modest augmentation of the basal level expressed by these cells and similar to that seen in human tumor cells. These findings are compatible with the observation that 293T cells display a significant basal level of IKK activation (Figs. 3 and 8).

Discussion
These data show that all cFLIP variants can activate IKK via a transient interaction of a ubiquitinated protein with the ubiquitin binding domain of IKK␥. The cFLIP L isoform requires the LUBAC ubiquitination complex but does not require caspase-8 or cFLIP L cleavage. In contrast, cFLIP S and p22-FLIP form a stable complex with FADD and RIP1, which are required for IKK activation. cFLIP L was previously reported to activate IKK following CD95 stimulation by the binding of a p43 fragment to IKK (6), and another study reported that the p22-FLIP cleavage product interacts directly with NEMO/IKK␥ to activate IKK (7). We have demonstrated that IKK activation following cFLIP L overexpression does not require cFLIP L cleavage. Furthermore, our evidence for different pathways of NF-B activation demonstrates that p22-FLIP is not simply the effector of cFLIP L , although cFLIP L cleavage to generate p22-FLIP may switch the mechanism of NF-B activation.
The cFLIP S ⅐FADD⅐RIP1 or p22-FLIP⅐FADD⅐RIP1 complexes are reminiscent of a cytoplasmic complex, termed the ripoptosome, induced by treatments of human tumor cells that deplete or inhibit cellular inhibitor of apoptosis proteins (38,39). This depletion leads to the formation of a RIP1⅐FADD⅐procaspase-8 complex. Cellular inhibitors of apoptosis inhibit formation of this complex by ubiquitinating RIP1. In the absence of FLIP FIGURE 8. Recombinant KSHV vFLIP activates IKK when added to cell lysate. A, IKK activation by recombinant (rec) vFLIP or recombinant p22-FLIP was analyzed using an in vitro IB␣ phosphorylation assay. Briefly, HEK293T cells were lysed in HTX lysis buffer, and the S100 fractions were extracted. The lysates were then incubated with an ATP-generating buffer and increasing concentrations of recombinant vFLIP or p22-FLIP at 37°C for 10 min. Activation of the IKK complex was detected by immunoblotting to detect phospho-IB␣ and phospho-IKK␣/␤. The blot was reprobed for total IKK␣/␤ to ensure equal loading of protein. B, a model of IKK activation by cFLIP isoforms and KSHV vFLIP. U, ubiquitin.
proteins, procaspase-8 molecules homodimerize in the ripoptosome, promoting cell death by apoptosis. Recruitment of cFLIP L to this complex leads to partial caspase-8 activation, resulting in cell survival (40). In contrast, cFLIP S has been reported to block caspase-8 activation and facilitate cell death through necroptosis (39). A complex containing cFLIP L , RIP1, FADD, caspase-8, and TRADD (tumor necrosis factor receptor 1-associated death domain) is induced by TNF receptor signaling (41); a complex containing RIP1, FADD, caspase-8, and IKK␥ has been detected after TRAIL triggering (42), and a cFLIP⅐RIP1⅐FADD⅐caspase-8 complex has been found after Fas signaling (43). Our data show that constitutive expression of cFLIP S or p22-FLIP leads to formation of a FLIP⅐FADD⅐RIP1 complex that leads to NF-B activation. The relationship between other complexes and this NF-B-activating complex remains to be determined, but we have clearly shown that caspase-8 is not required for NF-B activation by cFLIP S or p22-FLIP.
The expression of cFLIP proteins is tightly regulated in normal cells. At the transcriptional level, cFLIP isoforms are up-regulated by NF-B activation so they positively regulate their own expression (44). cFLIP S is particularly more prone to ubiquitination and displays a considerably shorter halflife compared with cFLIP L (5). This is largely due to the unique 19-amino acid C-terminal sequence of cFLIP S that is ubiquitinated at lysine residues 192 and 195; the stability of cFLIP R is similar to cFLIP S (3). p22-FLIP is more stable than cFLIP S as it lacks this C-terminal sequence (7). By using overexpression of cFLIP in our experiments, we have overridden some of these control mechanisms; however, there is evidence for overexpression of cFLIP isoforms in tumor cells (for a review, see Ref. 45), which is the situation that we have mimicked.
Up-regulation of cFLIP expression has been reported in many types of human solid tumors and lymphomas. We have confirmed these finding and shown that diverse human tumor cell types express a variety of cFLIP isoforms. There is also evidence that the level of cFLIP may determine clinical outcome (46 -49). cFLIP inhibition or down-regulation has therefore been proposed as a potential cancer therapy with the assumption that this will operate by sensitizing tumor cells to Fas-, TRAIL-, or TNF-induced cell death (for a review, see Ref. 45). However, tumor cells may also be dependent on cFLIP activation of NF-⌲B for their survival (49). Our dissection of the different mechanisms of cFLIP activation of IKK will allow the development of reagents to inhibit these pathways selectively and thereby assess their role in cFLIP-driven tumorigenesis.


text:
Bacillus pumilus SAFR-032 Genome Revisited: Sequence Update and Re-Annotation

Bacillus pumilus strain SAFR-032 is a non-pathogenic spore-forming bacterium exhibiting an anomalously high persistence in bactericidal environments. In its dormant state, it is capable of withstanding doses of ultraviolet (UV) radiation or hydrogen peroxide, which are lethal for the vast majority of microorganisms. This unusual resistance profile has made SAFR-032 a reference strain for studies of bacterial spore resistance. The complete genome sequence of B. pumilus SAFR-032 was published in 2007 early in the genomics era. Since then, the SAFR-032 strain has frequently been used as a source of genetic/genomic information that was regarded as representative of the entire B. pumilus species group. Recently, our ongoing studies of conservation of gene distribution patterns in the complete genomes of various B. pumilus strains revealed indications of misassembly in the B. pumilus SAFR-032 genome. Synteny-driven local genome resequencing confirmed that the original SAFR-032 sequence contained assembly errors associated with long sequence repeats. The genome sequence was corrected according to the new findings. In addition, a significantly improved annotation is now available. Gene orders were compared and portions of the genome arrangement were found to be similar in a wide spectrum of Bacillus strains.

Introduction
The endospore forming species, Bacillus pumilus, naturally occurs in the plant root systems of tobacco, pepper, cucumber and tomato, as well as in the leaves of soybean crops [1,2]. In this context, it can confer systemic protection against several plant diseases [3]. It was therefore unexpected when the strain B. pumilus SAFR-032 was isolated from a spacecraft assembly facility, and its spores were found to be extremely resistant to radiation, desiccation, and hydrogen peroxide treatment [4,5].
Subsequently, SAFR-032 spores were found to exhibit high survival rates under exposure to outer space conditions in experiments on board the International Space Station (ISS) [6,7]. The spaceflight data generated in these studies pointed to the risks of forward contamination of celestial bodies by microbial contaminants carried by landing probes. As a result, B. pumilus SAFR-032 has been classified as an extreme microorganism according to the planetary protection standards [7], and established as an important reference organism in studies on bacterial spore viability [8][9][10][11][12][13]. The unique resistance of B. pumilus SAFR-032 attracted the attention of microbiologists. This interest culminated in the publication of a complete annotated genome sequence of the SAFR-032 strain in 2007 [14]. This preceded by several years the genome of the B. pumilus type strain, ATCC 7061 T . Because of its early availability and importance as a model organism, the SAFR-032 genome is often used as a source of genetic information that is considered likely to be representative of the entire B. pumilus species group [15][16][17][18][19][20][21][22]. The SAFR-032 genome has also been used as a reference for assembly, correction and annotation of closely related genomes [23][24][25].
In recent years, complete genome sequences were determined for several other B. pumilus strains [19,23,26,27]. Their comparisons with the SAFR-032 genome revealed discrepancies in gene order, which were possibly due to errors in the original SAFR-032 genome assembly. New studies, described herein, have corrected what were in fact errors. In addition, the SAFR-032 genome annotation was updated to include recently identified genes and correct boundaries of known coding sequences.

Materials and Methods
Media, cultivation, and genomic DNA extraction B. pumilus SAFR-032 culture growth was initiated from the original stock archived at Jet Propulsion Laboratory (JPL), Pasadena, CA. The cells from the stock were streaked on Tryptic Soy Agar and incubated at 30°C until colonies appeared. A single colony was inoculated into 5 mL of Tryptic Soy Broth, and the culture was grown overnight at 30°C and 200 rpm. The genomic DNA was extracted from 1 mL aliquot using Maxwell 16 DNA/RNA Extraction System (Promega, Madison, WI), and eluted with water [28].

Primer design
Since the putative sequence misconnections were associated with the repeating genome segments (ribosomal RNA operons or transposase genes), the amplification primers were designed to bind to the sites immediately adjacent to the repeat units. Primer picking was performed using Primer3 v. 2.3.6 [29], primer binding specificity was confirmed using Blast 2.2.29 + [30]. For each questionable junction, 2 primer pairs were selected according to the following criteria: melting temperature 54-56°C, amplicon size not exceeding 10,000 bp, no significant binding outside of the intended priming sites, no significant secondary structure or self-binding. Amplification primers were also used for sequencing the terminal parts of the corresponding amplicons. The internal regions of the rDNA amplicons were partially assessed using universal primers targeting intergenic 16S-23S and 23S-5S spacers. Sequences of all primers used in this study are presented in Table 1.
Amplicon synthesis, purification and sequencing PCR amplification of the genomic DNA fragments covering questionable B. pumilus SAFR-032 genome regions was performed using Q5 High Fidelity PCR kit (New England Biolabs, Ipswich, MA). The reaction mixtures contained 0.16 μg of genomic DNA, 1 μM of each amplification primer and 40 μL of Q5 High Fidelity Master Mix in a final volume of 80 μL. After initial incubation at 94°C for 90 sec, the reactions were carried out for 37 cycles at 94°C for 30 sec, 52°C for 30 sec, and 72°C for 5 min. Final extension was performed at 72°C for 10 min. Then, the entire reaction mixtures were separated by electrophoresis on 1% agarose gel containing 0.2 μg/mL ethidium bromide. The bands of relevant size were excised from the gel, and the DNA was isolated from the gel slices on silica spin columns [31]. The purified rDNA amplicons were partially sequenced by the Sanger method at SeqWright, Inc (Houston, TX) using 5 primers for each amplicon: 16s_bpum_seq, 23s_bpum_seq1, 23s_bpum_seq2, and 2 amplification primers. Transposase gene amplicons were sequenced using amplification primers only.

Annotation
The re-assembled new version of the B. pumilus SAFR-032 genome was annotated using the NCBI Prokaryotic Genome Annotation Pipeline (PGAP) [36] and the RAST pipeline [37]. The results were compared with those of the existing annotated genome (NC_009848.1/ CP000813.1). The comparison was conducted by manually aligning the predicted ORFs/genes (including the flanking 20-30 bases upstream) using the BioEdit sequence alignment software [38]. . Since some of the genomes were deposited in reverse-complement orientation to a standard layout, they were inverted to bring their orientation to normal. The first position of each genome was set at the first position of dnaA gene.

Identifying and Correcting Sequencing Problems
Six de novo assembled complete genomes of B. pumilus strains SAFR-032 (CP000813.1), WP8 (CP010075.1), MTCC B6033 (CP007436.1), NJ-M2 (CP012329.1), NJ-V2 (CP012482.1), and TUAT1 (AP014928.1) were aligned using the Mauve aligner to investigate conservation of B. pumilus genome layout. The original assembly of the SAFR-032 genome deviates from the consensus gene order near the origin of bacterial chromosome replication (Fig 1A). A conserved genome fragment encompassing 19 protein-encoding genes, 1 stand-alone tRNA gene and 1 rRNA operon was found displaced. This fragment is 522 kbp downstream of its location in the other genomes as related to the position of chromosomal replication initiation protein dnaA.
Another SAFR-032 genome segment containing 13 protein-encoding genes and 1 tRNA gene, and flanked by transposase coding sequences was found to be inverted as compared to the genomes of other strains ( Fig 1B). Both problematic fragments are flanked with long sequence repeats, rRNA operons in one case and transposase genes in the other. This makes it look like genome assembly errors were caused by incorrect contig joining over repeat regions. Better agreement with other B. pumilus strains can be achieved by moving a stretch of genes BPUM_0503-BPUM_0520 and the ribosomal operon 5 (BPUM_r0013, BPUM_r0013, BPUM_r0015) immediately downstream of the ribosomal operon 1 (BPUM_r0001, BPUM_r0002, BPUM_r0003), and by inverting a genome segment harboring genes from BPUM_0842 to BPUM_054 (Fig 2A). Validity of this rearrangement has been confirmed by selective PCR amplification of the affected SAFR-032 genome regions followed by dye dideoxy terminator sequencing of the amplicons (Fig 2B).

Comparison of New PGAP Annotation with Original Annotation and Resolution of Differences
Although the SAFR-032 genome rearrangement did not affect the existing gene boundaries, the genome sequence was re-annotated to bring it to the most up-to-date state. The primary annotation was performed with NCBI PGAP [36] and the results were compared with the original (CP000813.1) annotation. A surprisingly large number of differences were found and significant effort was required to resolve them.
To begin with, the automated PGAP annotation omitted 170 genes as well as all non-coding RNAs that have been present in the original genome annotation. These genes were checked for their veracity and those that were found to be genuine ORFs based on their homology to other genomes in the NCBI databases, were manually incorporated into the new annotation (S1 Table). Known non-coding RNAs were also added back. The revised version of the annotated B. pumilus SAFR-032 genome sequence was deposited in Genbank under accession number CP000813.4.
In addition, a total of 385 ORFs were found to differ between the new PGAP annotation and the original (CP000813.1) annotation in terms of the predicted location of the start codon  Table). These questionable gene boundaries were examined using the RAST pipeline [37]. In 197 cases RAST supported the new PGAP ORF coordinates. In 118 cases RAST placed the start codon position as it was in CP000813.1, and thus the earlier annotation was retained. Finally, in 70 cases RAST did not support either annotation. Hands on visual examination made it clear that the problematic ORFs do not have recognizable ribosome-binding sites in front of the candidate start codons. In the absence of experimental data, this makes the prediction of ORF start positions in these 70 cases unreliable, and does not allow a clear choice between several available alternatives. Therefore in the final updated annotation (CP000813.4), for consistency reasons, the ORFs with ambiguous coordinates were left as annotated by PGAP.
In the PGAP annotation, 58 entries were marked as pseudogenes (S3 Table). The two annotations agreed on 35 pseudogenes. The new PGAP version had 23 extra pseudogenes, 22 of which were previously regarded as valid genes, and 1 was newly identified. After examination of sequence alignments with close homologs, it was concluded that nine of these are indeed genuine pseudogenes with either in-frame stop codons or deletion of a substantial part of the coding sequence. The remaining fourteen examples rather correspond to properly expressed genes as it follows from homology analysis. Thus, only 44 pseudogenes were retained in the final CP000813.4 annotation.

BPUM_2060: A Special Case
BPUM_2060 was annotated by PGAP as a pseudogene of a "bifunctional GTP cyclohydrolase II/3,4-dihydroxy-2-butanone 4-phosphate synthase" (ribA). This gene is a part of the four gene operon ribTHAE that codes for enzymes involved in the riboflavin biosynthetic pathway [39]. The pseudogene status was assigned to BPUM_2060 because it maps to only the N-terminal half of the reference protein (WP_006636466.1, product of ribA gene from Bacillus sonorensis), which makes BPUM_2060 look like a truncated version of the reference. However, it appears that in some Bacillus species, GTP cyclohydrolase II and 3,4-dihydroxy-2-butanone 4-phosphate synthase are two distinct polypeptides encoded by separate genes. In B. pumilus SAFR-032, BPUM_2060 gene encodes 3,4-dihydroxy-2-butanone-4-phosphate synthase while two other genes, BPUM_0700 and BPUM_3531, each encodes GTP cyclohydrolase II. Besides B. pumilus strains, the same situation exists in B. safensis, B. altitudinis, B. xiamenensis, B. stratosphericus, B. isronensis, B. simplex, and B. butanolivorans. It can also be seen outside of the Bacillus genus, for example, in Planococcus donghaensis, P. halocryophilus, P. antarcticus, Solibacillus silvestris, Lysinibacillus fusiformis, L. sphaericus, L. varians, Paenibacillus polymyxa, Kurthia massiliensis, K. huakuii, and Oceanobacillus manasiensis. Therefore, BPUM_2060 was re-annotated as a valid gene. This example demonstrates that an improperly selected reference may cause an annotation artifact. Clearly, in the present day scenario of rapidly evolving sequencing technologies, automated annotation of genomes is likely to be error-prone, and manual curation is absolutely necessary to maintain better standards of quality of a genome and its various components.

Final Updated Annotation
Overall, the updated B. pumilus SAFR-032 genome annotation contains 3710 valid proteinencoding genes, 44 likely pseudogenes, 21 rRNA genes grouped in 7 rRNA operons, 72 tRNA genes, one tmRNA, 4 other non-coding RNA genes, and 16 riboswitches. The updated annotation includes 38 genes and 2 riboswitches that were not present in the original CP000813.1 genome annotation (S4 Table). Three of the newly identified genes, BPUM_03795, BPUM_04085, and BPUM_04135, are transposases, which brings a total number of mobile genetic elements in B. pumilus SAFR-032 genome to 20. Another two genes, BPUM_04115 (coding for the "capsular biosynthesis protein CpsH") and BPUM_04155 (coding for a "spore coat protein") with probable roles in cellular structural integrity and sporulation, respectively, may be an addition to the candidate genes involved in the unusual resistances exhibited by B. pumilus SAFR-032 [10,11]. Among the rest, seven genes are predicted to be involved in transport and secretion, two genes encode regulatory proteins, three genes encode enzymes, and five genes are predicted to encode stable cellular RNAs (three tRNAs, tmRNA, and 6S RNA). The final revised version of the annotated B. pumilus SAFR-032 genome sequence was deposited in Genbank under accession number CP000813.4.

Discussion
High-throughput shotgun sequencing methods generate a steady stream of genomic sequences filling public databases. Currently, the majority of deposited prokaryotic genomes have draft status. Thus, they represent a collection of contigs or scaffolds rather than a single circularly closed chromosome. While being acceptable for mutation mapping, strain typing, and metabolic pathways identification, the draft genome assemblies are not sufficient for studying largescale genome architecture, which is needed for better understanding of genome function and evolution. This creates a demand for completed genome sequences, even though the assembly finishing is costly and time-consuming. The major obstacle on the way towards a complete genome is posed by long sequence repeats, which may cause assembly errors because of ambiguities in contig order. The present study demonstrates that such errors can be revealed using multiple alignments of closely related genomes.
All de novo assembled complete B. pumilus genomes shows absolute conservation of gene order in the region immediately downstream of the replication origin, between dnaA and metS (metG) genes. This observed genetic layout is well preserved far beyond the B. pumilus strains. With minor modifications it can be seen in B. subtilis, B. mojavensis, B. atrophaeus, B. amyloliquefaciens, B. licheniformis, B. megaterium, B. endophyticus, B. smithii, B. cytotoxicus, in the entire B. cereus group of species, and throughout Geobacillus and Anoxybacillus genera (S1 Fig). Thus, the deviation of B. pumilus SAFR-032 genome sequence from such a pronounced trend strongly indicated a high possibility of an assembly error. Direct sequencing of the problematic locus confirmed the presence of an error. The corrected SAFR-032 sequence exhibits perfect agreement with the consensus gene distribution pattern near the replication origin.
The second problematic segment in the SAFR-032 genome covering genes from gpmB (BPUM_0834) to cspB (BPUM_0862) is significantly less conserved. Among six de novo assembled B. pumilus genomes, only strains WP8, NJ-M2, NJ-V2, and TUAT1 contain ssuBACD operon in this location while in SAFR-032 and MTCC B6033 it is missing. The SAFR-032 gene BPUM_0842 encoding an YncM-like protein maps to the QR42_04470 gene in the aligned fragment of the WP8 genome but is absent in the corresponding genomic regions of other strains. Two identical transposase genes are present only in the SAFR-032 genome fragment and flank the gene stretch from BPUM_0842 to BPUM_0854, which was inverted in the earlier B. pumilus SAFR-032 genome version, CP000813.1, as related to other B. pumilus genomes. The relatively low local conservation of the gene order has made the indications of sequence misassembly look weaker than in the previous case. Yet, the reversal of a 11,000 bp-long genome fragment was regarded as a sufficiently substantial deviation from the consensus gene pattern to justify its re-evaluation.

Conclusions
Possible assembly errors in the previously published B. pumilus SAFR-032 genome were evaluated and corrected. To accommodate these changes, re-annotation was required. Despite the fact no coding sequences were changed this led to significant changes in the annotation that required considerable effort to evaluate. This points to continuing difficulties in automated annotation and highlights the importance of being cautious when relying on it. The immediate consequences of the SAFR-032 genome sequence update reflect the organism's long-standing role as a reference for closely related species. The SAFR-032 genome was previously used as an assembly template in reference-driven assembly of B. pumilus W3 genome, which in its own turn was used as a reference for B. pumilus ku-bf1 genome assembly. Therefore, it is not surprising that both W3 and ku-bf1 genomes exhibit the same deviant layout that was seen in the CP000813.1 version of the SAFR-032 genome (Fig 3 and S2 Fig). As in the case of B. pumilus SAFR-032, these two genomes may need re-evaluation. Furthermore, all intergenomic synteny studies that have involved B. pumilus SAFR-032 strain may need to be critically revised.


extfieldsofstudy: Chemistry Medicine

text:
Datasets from a vapor diffusion mineral precipitation protocol for Dictyostelium stalks

Datasets from a slow carbonate vapor diffusion and mineral precipitation protocol for Dictyostelium ECM and cellulose stalks show examples for composite materials obtained by an in vitro approach, which differs substantially from the in vivo approach reported in The Journal of Structural Biology, doi: 10.1016/j.jsb.2016.03.015 [1]. Methods for obtaining the datasets include bright field transmitted light microscopy, fluorescence microscopy, LC-PolScope birefringence microscopy, variable pressure scanning electron microscopy (VP-SEM/ESEM), and Raman imaging spectroscopy.

Data
Datasets show fluorescence light microscopy, LC-PolScope/Abrio imaging, scanning electron microscopy (VP-SEM), and confocal Raman imaging spectrometry data of composite materials obtained from Dictyostelium extracellular matrix (ECM). Different parts of native and genetically modified ECM were subjected to slow precipitation of calcium carbonate in aqueous Ca 2 þ solution with ammonium carbonate in the vapor phase. Datasets of specimens modified with different protein derivatives, either from egg shell (OC-17) or from mollusc shells (n16N, perlucin) investigated in wet or dry conditions and without any further modification are presented.

Cell lines
Synthetic sequences of n16N, OC-17 and perlucin were cloned into modified Gateway expression systems for Dictyostelium Ax3-Orf þ under the control of the [ecmB] promoter [1]. The heterologous proteins were expressed with a C-terminal GFP domain for detection, and an N-terminal Dictyostelium signal peptide to facilitate their export into the extracellular matrix (ECM) which is formed during the multicellular stages of the Dictyostelium life cycle.

In vitro mineralization experiments
Dictyostelium ECM mineral composites were produced under slow ammonium carbonate vapor diffusion adapted from in vitro precipitation assays [2][3][4][5]. The 2-3 days old organic stalks of differentiated Dictyostelium cell lines were collected using tweezers and spore heads were gently removed.
The μ-slide was covered with gently perforated aluminum foil. Slides were incubated in humid atmosphere 12 h to 2 days in a closed incubator at room temperature containing 0.5 g of (NH 4 ) 2 CO 3 until the pH increased to 8.7 70.2. For drying the specimens, the Ca 2 þ containing salt solution with suspended organic stalks was quickly transferred onto nitrocellulose filters (0.2 μm pore diameter, Millipore) on top of a filter holder mounted on a vacuum pump, rinsed twice with 100 μl deionized water and immediately air dried under vacuum.

Transmitted light microscopy and fluorescence microscopy
Specimens were investigated using an inverted light microscope Cell Observer Z1 equipped with a mercury arc lamp and a Colibri LED fluorescence excitation system and the objectives A-Plan 10 Â / 0.25Ph1 and LD Plan Neofluar 40 Â /0.6Korr Ph2 (Zeiss, Göttingen, Germany). Electronic images were recorded with an Axiocam MRm Camera and Axiovision 4.8 software. Transmitted light microscopy images were taken at automatic exposure times. Fluorescence microscopy images were taken at a fixed exposure time of 10 s using a band pass filter (excitation 470/40, emission 525/50) at 50% Colibri LED intensity (Fig. 1).

LC-PolScope quantitative birefringence microscopy
The birefringence of specimens was analyzed using a LC-PolScope image processing system operated by the Abrio software (CRI, LOT-Oriel, Darmstadt, Germany) [6,7]. This system was mounted on the inverted light microscope Cell Observer Z1 (Zeiss Observer Z1, Göttingen, Germany). All background images were captured under identical conditions as the area of interest, with the only difference being the absence of the specimen. The measured birefringent retardance is encoded per pixel and displayed in gray scale or false colors. The retardance values represent the relative phase shift between two orthogonally polarized light waves after traversing an optically anisotropic material. Orientation values indicate the azimuth of the slow optical axis (Figs. 2 and 3).

Environmental scanning electron microscopy (ESEM)
Scanning electron microscopy was performed using a Quanta FEI 400F (FEI, Eindhoven, Netherlands) variable pressure SEM (VP-SEM) in either the environmental scanning electron microscopy mode (ESEM) or the low-vacuum SEM mode. Specimens were imaged in native state without additional coating. For ESEM, the stage was cooled to 3°C before sample application and the chamber was purged 8 times between 800 Pa and 1200 Pa during vacuum pumping. Images were captured with a secondary electron detector at an operating voltage between 10 kV and 30 kV, a spot size of 3, and a pressure of 760 Pa. The low-vacuum mode was commonly used with a secondary electron detector at 10-20 kV and 100 Pa pressure for dry specimens either directly on nitrocellulose membranes or mounted on a silicon wafer (Figs. 4 and 5).

Raman imaging spectroscopy
Raman spectra were collected using a confocal Raman spectrometer (Aramis Labram, Horiba Jobin Yvon, Bensheim, Germany) equipped with an peltier cooled detector. Samples were excited with 785 nm laser and captured using 50 Â and 100 Â objectives, a lattice of 600 over a scan range of 100-2000 cm À 1 and an exposure time of 10 s and 5 integrations. For confocal Raman spectra, the   7. Raman imaging microscopy data (left column) and low-vacuum VP-SEM images (overview, and region 3 of the stalk at higher magnification) of AX3_OC17GFP. Raman spectra (right) of mineral precipitates attached to region 1 (z1, z2) and region 3 (z3, z4) with corresponding images taken with the Raman microscope (left) are shown. A focal series of Raman spectra was collected manually from the mineral precipitate (z1, region 1) and underneath the mineral in the matrix (z2, region 1). Fig. 8. Low-vacuum VP-SEM images (top) and Raman imaging spectroscopy data (bottom) obtained from of a mineral particle precipitated in region 1 of AX3_n16NGFP (top right, dashed lines). The corresponding light microscopy image captured with the Raman microscope is shown in the inlet (top right). Raman spectra were collected in the confocal mode with 5 steps starting in the nitrocellulose filter (z1) through the mineral precipitate (z2 to z4) and ending slightly above the particle (z5). Spectra are adjusted to the height of the 1086 cm À 1 band. pinhole was closed to 100 μm. For each mineral precipitation experiment, spectra from at least three independent experiments were investigated (Figs. 6-8).