BACKGROUND
the term t cell repertoire describes a collection of lymphocytes characterized by t cell receptor  expression, which plays a critical role in antigen recognition.
since alterations of the t cell repertoire provide a significant indication of immune status in physiological and disease conditions, t cell repertoire analyses have been conducted for the identification of antigen-specific t cells involved in the development of disease and for the diagnosis of t lymphocyte abnormalities.
comparison of variable-region usage by fluorescence-activated cell sorter analysis using a large panel of antibodies specific for tcr variable regions , polymerase chain reaction  with multiple primers  <cit>  or pcr-based enzyme-linked immunosorbent assay  <cit>  have been widely used to detect changes in t cell repertoire.
length distribution analysis known as cdr <dig> spectratyping is based on the addition of non-template nucleotides in v--j region and has been used to evaluate t cell clonality and diversity  <cit> .
to identify the antigen specificity of t cells further, pcr cloning of tcr clonotypes and subsequent sequence determination of the antigen recognition region, cdr <dig>  have been required.
these conventional approaches are commonly used but are time-consuming and a laborious way to study tcr repertoires.
in recent years, advances in high-throughput sequencing technologies known as next-generation sequencing  have rapidly progressed and enabled large-scale analysis of sequence data  <cit> .
although several ngs-based tcr repertoire analysis systems have been developed by other researches, many amplification techniques are based on multiple pcr with different primers specific for each variable region.
thus, bias during pcr amplification is unavoidable since bias is most commonly due to differential hybridization kinetics among variable region-specific primers to different target genes.
correction and additional computational normalization methods are therefore required to minimize pcr bias when using multiple pcr assays  <cit> .
the use of a single set of primers is a better way to achieve unbiased and quantitative amplification of all tcr genes including unknown variants where the 5’ ends of sequences are highly diverse.
a single strand oligonucleotide anchor ligation to the 3’ end of cdna with t <dig> rna ligase  <cit> , homopolymeric tailing of cdna, 5’ rapid amplification of cdna ends   <cit>  and template switching pcr   <cit>  have been used to analyze tcr repertoires  <cit> .
ts-pcr is simple and convenient but produces high levels of background amplification because ts primers non-specifically anneal to random regions in rna or allow the repeated addition of ts primers  <cit> .
thus, the current study describes an adaptor-ligation mediated pcr  developed by the addition of an adaptor to the 5’ end of double stranded  cdna from tcr transcripts and subsequent pcr amplification with the adaptor primer and constant region-specific primer, as first reported by tsuruta et al.
<cit> .
the adaptor ligation to blunt-ended ds cdna is less influenced by the sequence of a particular cdna while the efficiency of 5’ adaptor ligation with t <dig> rna ligase is sequence dependent  <cit> .
in addition, the ligation of dsdna by t <dig> ligase is more efficient than ssdna ligation with t <dig> rna ligase in ligation anchored pcr .
various sequencing technologies such as roche  <dig> , illumina , ion-torrent , solid , helicos  and pacbio  have been developed.
among these ngs platforms, the  <dig> dna sequencing produces sequence reads ranging from  <dig> to  <dig> base pairs  or more in length and sufficient read outputs, yet less reads per run than the illumina.
long read sequencing allows determination of the full or near-complete length of tcr genes including v, d, j and c regions.
furthermore, recombinant tcr proteins can be easily produced by subsequent pcr cloning of the tcr genes.
therefore, we applied an adaptor-ligation mediated pcr method to ngs with  <dig> dna sequencing.
natural killer t  cells are a distinct t cell population with an important role in innate and adaptive immunity.
nkt cells regulate a broad range of immune responses such as autoimmune diseases, tumor surveillance, and host defense against pathogenic infections.
nkt cells express an invariant tcrα consisting of vα <dig> and jα <dig> that recognizes glycolipids presented by a non-classical major histocompatibility complex class i-related protein, cd1d  <cit> .
recently, mucosal-associated invariant t  cells, which preferentially exist in mucosal tissues, were shown to be a unique t cell population expressing a semi-invariant tcrα consisting of vα <dig>  and jα <dig>  mait cells recognize microbial vitamin b metabolites presented by a non-classical mhc class i molecule, mhc-related protein  <dig>   <cit> .
these t cell populations bearing invariant tcrα play a pivotal role in immune regulation but it remains to be determined whether all invariant tcrα are expressed by these unique t cell populations.
in this study, we conducted ngs sequencing of tcr transcripts from  <dig> healthy individuals using a newly developed ngs-based tcr repertoire analysis.
initially, based on sequence read count, we examined usages of variable and joining regions, and further analyzed clonality and diversity in tcrα and β genes.
unique sequence reads identified using an originally developed gene analysis program were compared at a clonal level among healthy individuals.
these results showed a similar usage of trv and trj and similar extent of diversity of t cells among individuals.
interestingly, tcrβ reads were less shared among individuals while tcrα reads frequently contained shared sequences that overlapped between two or more individuals.
shared tcrα reads contained a high proportion of invariant tcrα indicating the presence of inkt cells or mait cells.
in this report, we demonstrated that analysis of tcr genes shared among multiple individuals from ngs data provided significant information on invariant tcrs expressed by nkt cells and mait cells.
methods
isolation of peripheral blood mononuclear cells and rna extraction
after obtaining written informed consent, whole blood was collected from  <dig> japanese healthy individuals .
the study was approved by the ethics committees of the clinical research center for rheumatology and allergy, sagamihara national hospital, national hospital organization.
ten ml of whole blood was collected into heparinized tubes.
peripheral blood mononuclear cells  were isolated with ficoll-paque plus™  density gradient centrifugation and washed with phosphate buffered saline .
cell numbers were counted and 1 ×  <dig> cells were used for rna extraction.
total rna was isolated and purified with rneasy lipid tissue mini kit  according to the manufacturer’s instructions.
rna amounts and purity were measured with agilent  <dig> bioanalyzer .
unbiased amplification of tcr genes
one microgram of total rna was converted to complementary dna  with superscript iii reverse transcriptase .
bsl-18e primer containing polyt <dig> and a noti site was used for cdna synthesis.
after cdna synthesis, double strand -cdna was synthesized with e. coli dna polymerase i , e. coli dna ligase , and rnase h .
ds-cdnas were blunted with t <dig> dna polymerase .
p10ea/p20ea adaptor was ligated to the 5’ end of the ds-cdna and then cut with noti restriction enzyme.
after removal of adaptor and primer with minelute reaction cleanup kit , pcr was performed using either tcr α-chain constant region-specific  or tcr β-chain constant region-specific primers  and p20ea .
pcr conditions were as follows: 95 °c , 55 °c , and 72 °c  for 20 cycles.
the second pcr was performed with either ca <dig> or cb <dig> and p20ea primers using the same pcr conditions.table  <dig> primers used in this study
ccatctcatccctgcgtgtctccgac
ccatctcatccctgcgtgtctccgac
v: a/c/g, n: a/c/g/t, adaptor a and b sequences were typed in bold and bold italic, respectively.
a key sequence  was underlined.
the following mid tag sequences were used for identification of sample source.
mid1: acgagtgcgt, mid2: acgctcgaca, mid3: agacgcactc, mid4: agcactgtag, mid5: atcagacacg, mid6: atatcgcgag, mid7: cgtgtctcta, mid8: ctcgcgtgtc, mid10: tctctatgcg, mid11: tgatacgtct, mid15: tacgacgta, mid16: tcacgtacta, mid17: cgtctagtac, mid18: tctacgtagc, mid19: tgtactactc, mid20: acgactacag, mid21: cgtagactag, mid22: tacgagtatg, mid23: tactctcgtg, mid24: tagagacgag
amplicon sequencing by roche  <dig> sequencing system
amplicons for ngs were prepared by amplification of the second pcr products using p20ea primer and fusion tag primer .
the fusion tag primers consisting of an a-adaptor sequence ,  <dig> base sequence key , multiple identifier  tag sequence , and tcr constant region-specific sequence were designed according to the manufacturer’s instructions.
after pcr amplification, amplicons were separated and evaluated by agarose gel electrophoresis.
the resulting fragment  was removed from the gel and purified with qiaex ii gel extraction kit .
the amount of purified amplicon was quantified by quant-it™ picogreen® dsdna assay kit .
each amplicon obtained with a different fusion tag primer from  <dig> healthy individuals was mixed at equal molar concentrations.
emulsion pcr  was performed using the amplicon mixtures with gs junior titanium empcr lib-l kit  according to the manufacturer’s instructions.
assignment of trv and trj segments
all sequence reads were classified by mid tag sequences.
artificially added sequences  and sequences with low quality scores were removed from both terminals of sequence reads using software installed on the  <dig> sequencing system.
the remaining sequences were used for assignment of trav and traj for tcr α sequences, and trbv and trbj for tcr β sequences.
assignment of sequences was performed by determining sequences with the highest identity in a data set of reference sequences for  <dig> trav,  <dig> traj,  <dig> trbv and  <dig> trbj genes including pseudogenes and open reading frame  reference sequences available from the international immunogenetics information system®  database .
data processing, assignment, and data aggregation were automatically performed using repertoire analysis software originally developed by our group .
rg implemented a program for sequence homology searches using blatn, an automatic aggregation program, a graphics program for trv and trj usage, and cdr <dig> length distribution.
sequence identities at the nucleotide level between query and entry sequences were automatically calculated.
parameters that increased sensitivity and accuracy  score) were carefully optimized for respective repertoire analysis.
data analyses
nucleotide sequences of cdr <dig> regions ranged from conserved cysteine at position  <dig>  of imgt nomenclature to conserved phenylalanine at position  <dig>  and the following glycine  was translated to deduced amino acid sequences.
a unique sequence read  was defined as a sequence read having no identity in trv, trj and deduced amino acid sequence of cdr <dig> with the other sequence reads.
the copy number of identical usr were automatically counted by rg software in each sample and then ranked in order of the copy number.
percentage occurrence frequencies of sequence reads with trav, traj, trbv and trbj genes in total sequence reads were calculated.
retrieval of shared usrs among samples
to retrieve shared sequences among samples, a concatenate string of “trv gene name”_” deduced amino acid sequence of cdr <dig> region”_” trj gene name” of individual usr  was used as a tcr identifier.
the tcr identifier in a sample was retrieved in read data sets from all the other samples.
diversity indices and similarity index
to estimate tcr diversity in deep sequence data, several diversity indices, simpson’s index and shannon-weaver index were calculated using a function “diversity” of the vegan package in the r program.
these indices were calculated based on the number of species per sample and the number of individuals per sample as measures for biological diversity in ecology.
in deep sequence data, usr and copy number were used for species and individuals, respectively.
simpson’s index  was defined as: 1−λ=1−∑i=1snini−1nn− <dig>  where n is the total number of sequence reads, ni is the copy number of ith usr, and s is the species number of usr.
this value ranges from  <dig> to  <dig>  where the maximum number  <dig> means high levels of diversity and  <dig> indicates low diversity.
simpson’s reciprocal index  was also calculated as the inverse of λ. the shannon-weaver index  was used for the diversity index and was defined as: h'=−∑i−1sninlnnin where n is the total number of sequence reads, ni is the number of ith usr, and s is the species number of usr.
these diversity indices should be biased by differences in read numbers among samples.
therefore, the number of sequence reads was standardized for each sample down to the smallest number of sequence reads  <cit> .
to standardize sample size, repeated random resampling  <dig> times without replacement and calculation of the diversity index were performed using the r program.
the median of their indices was used to determine the diversity index for the sample.
to estimate the similarity of tcr repertoires between healthy individuals, the morisita-horn index was calculated using the function “vegdist” in the vegan package of the r program.
the morisita-horn index  was defined as: ch=2∑i=1sxiyi∑i=1sxi2x2+∑i=1syi2y2xy where xi is the number of ith usr in the total x reads of one sample, yi is the number of ith usr in the total y reads of another sample, and s is the number of usr.
to standardize the sample size, repeated random resampling  <dig> times without replacement and calculation of similarity index were performed using the r program  <cit> .
median values were used for similarity indexes between a pair of samples.
statistics
statistical significances were tested by the nonparametric mann–whitney u-test using graphpad prism software .
a value of p <  <dig>  was considered statistically significant.
RESULTS
repertoire analysis software
the cloud-based software platform developed in this study, rg, is a high-speed, accurate and convenient computing system for tcr repertoire analysis.
rg provides an integrated software package for 1) assignment of v, d and j segments, 2) calculation of sequence identity, 3) extraction of cdr <dig> sequences, 4) counting identical reads, 5) amino acid translation, 6) frame analysis , and 7) cdr <dig> length analysis.
after uploading sequencing data from the ngs sequencer, v, d and j segments can be identified based on their sequence similarity with optimized parameters.
read counts are automatically summarized and then processing data, summary tables, and graphs can be easily downloaded.
read number, error rate, and unproductive reads
we performed high-throughput sequencing of tra and trb genes in pbmcs from  <dig> healthy individuals.
a total of  <dig>  and  <dig>  sequence reads were assigned for tra and trb repertoire analyses, respectively using the rg program .
a total of  <dig>  and  <dig>  unique sequence reads  were identified in tra and trb, respectively.
the number of nucleotide sequences per read obtained by roche  <dig> sequencing were ~400 bp in length , indicating these sequences were long enough to identify tcr genes ranged from v region to j regions.
to evaluate the accuracy and quality of ngs sequencing, we calculated the frequency of mismatched nucleotides between query and reference sequences as the error rate.
error rates were  <dig>  ±  <dig>  % for trav,  <dig>  ±  <dig>  % for traj,  <dig>  ±  <dig>  % for trbv and  <dig>  ±  <dig>  % for trbj .
these error rates were slightly lower than in a previous study reporting a mean error rate of  <dig>  % for 454-sequences  <cit> .
the error rates were significantly higher in v regions than in j regions , indicating higher sequence reliability in regions close to sequencing primers.
occurrence frequency of unproductive reads carrying a stop codon or a shift of reading frame in cdr <dig> regions  was calculated .
there was no significant difference in the percentage of occurrence frequency of unproductive unique sequence reads between tra and trb .
regarding the error rates, frequencies of mismatched nucleotides were significantly higher in out-of-frame reads than in in-frame reads in all regions of trav, traj.
trbv and trbj .
the differences were quite significant in j regions .
this suggested that sequencing errors occurred in j regions adjacent to cdr <dig> cause frequently frame-shift of coding sequence.
expression of tcr genes including pseudogenes and orf
to determine usages of trv and trj genes in tcr sequence reads, copy number  of usrs bearing respective trv or trj were counted.
individual usrs were ranked by order of copy number and the percentage frequency of respective trv and trj were calculated .
regarding tra repertoires, eight pseudogenes  were not expressed in healthy individuals.
av8- <dig>  classified as an orf , were slightly expressed .
expression of av <dig> and av <dig>  was not observed in healthy individuals.
in addition, the functional genes, av <dig> and av9- <dig>  were poorly expressed in  <dig> individual  and  <dig> individuals , respectively.
of eight aj genes  classified as orf genes, the expression of aj <dig> and aj <dig> were observed in all  <dig> individuals.
of these, aj <dig> and aj <dig> were slightly expressed in  <dig> individuals  and  <dig> individuals , respectively.
aj <dig>  aj <dig>  aj <dig> and aj <dig> were not present in any individuals.
there was no expression of the three pseudogenes, aj <dig>  aj <dig> and aj <dig>  in any individuals.
the functional gene aj <dig> was detected in only  <dig> reads from  <dig> individuals.fig.
<dig> usage of trav and traj in  <dig> healthy individuals.
the numbers of tcr sequences bearing respective trav and traj were counted.
the percentage frequencies of  <dig> trav and  <dig> traj were calculated and are shown as a scatter plot.
each dot indicates the percentage frequency of trav or traj in each individual.
bars in red indicate the mean values of  <dig> individuals.
: pseudogene, : open reading frame
fig.
<dig> usage of trbv and trbj in  <dig> healthy individuals.
the percentage frequencies of  <dig> trbv and  <dig> trbj are shown as a scatter plot.
each dot indicates the percentage frequency of trbv or trbj in each individual.
bars in red indicate the mean values.
: pseudogene, : open reading frame
for trb genes, there was no expression of  <dig> pseudogenes  in healthy individuals.
of five orf genes, bv5- <dig> , bv6- <dig>  and bv <dig>  were poorly expressed.
the bv7- <dig> orf gene was not observed in any individuals while bv23- <dig> was expressed in all  <dig> individuals.
regarding bj genes, there was no expression of the bj2-2p pseudogene.
infrequent recombination of trav and traj
gene recombination of  <dig> trav with  <dig> traj  were capable of generating a total of  <dig> av-aj recombinations , of which  <dig> av-aj recombinations  were detected in  <dig> individuals.
this indicated that almost all av-aj recombinations were used in tcr transcripts without restriction.
notably, av1-1–av <dig> genes were recombined less preferentially with aj50–aj <dig> genes and similarly, recombination of av35–av <dig> genes with aj3–aj <dig> was rarely observed.
given the chromosomal location of these gene segments, these results showed that av-aj recombinations occurred infrequently between proximal av genes and distal aj genes and between distal av genes and proximal aj genes.
for trb,  <dig> gene recombinations were generated with  <dig> bv  and  <dig> bj genes , of which  <dig> bv-bj  were used in  <dig> individuals.
there was no restriction for combinations of trbv with trbj.fig.
<dig> heatmap representation of gene recombination of trav with traj in pooled read data from  <dig> healthy individuals.
the numbers of tcr sequence reads bearing respective gene recombination of trav and traj were counted.
the occurrence tendency of recombination is visualized by heat map presentation of the number of each recombination.
color in each pixel indicates the number of each recombination.
for trav,  <dig> pseudogenes ,  <dig> orf , and poorly expressed genes  were excluded.
for traj,  <dig> pseudogenes ,  <dig> orfs , and poorly expressed genes  were excluded.
two orfs  found to be expressed were included.
a heat map representation of the  <dig> recombination events  is shown
preferential usage of trv and trj repertoires in healthy individuals
to reveal trv and trj usage in whole tcr transcripts, the occurrence frequency of usrs bearing respective trv or trj were calculated .
preferential usage in several trav genes was observed in trav <dig> , trav9- <dig> , trav13- <dig>  and trav <dig> .
these preferential usages were similar to a previous result obtained using a hybridization-based quantitation assay  <cit> .
several trbv genes were abundantly used in the trbv repertoire.
the top three trbv29- <dig> , trbv20- <dig>  and trbv <dig>  accounted for approximately one third of the total sequence reads.
this was similar to results from our previous study using the microplate hybridization assay   <cit> .
gene usage varied considerably among the trbj genes.
trbj2- <dig> and trbj 2– <dig> were highly expressed while the expressions of trbj1- <dig>  trbj1- <dig>  trbj1- <dig>  trbj2- <dig> and trbj2- <dig> were low.
next, to examine if the preferential usage of trv and trj are due to peripheral selection, we compared usages of tcr between in-frame and out-of-frame reads .
as a whole, we observed little difference in the usages of trav, traj, trbv and trbj between in-frame and out-of-frame reads.
on the other hand, particularly trav26- <dig>  traj <dig>  traj <dig>  trbv23- <dig>  trbj1- <dig> and trbj2- <dig> were more frequently used in out-of-frame reads than in-frame reads.
three-dimensional  view of tcr repertoire usage
to visualize the usage of tcrs bearing a combination of trv with trj genes, we produced 3d pictures of the tcr repertoires .
the advantage of 3d images is that the dominance of certain combination of trv with trj genes as well as the extent of diversity of tcr can easily be observed.
for trb, there was little preferential usage of recombination between trbj and trbv genes.
the frequency of each recombination depended on the usage of trbv or trbj.
bv29-1/bj2- <dig>  bv29-1/bj2- <dig>  bv29-1/bj2- <dig> and bv20-1/bj2- <dig> were frequently used in all combinations while others were expressed at a low frequency.
in contrast, 3d imaging of the tra repertoire indicated a low level of expression with a wide distribution of trav and traj.
the occupancy was lower than 1 % for all combinations.
of note, tcr reads bearing av1- <dig> and aj <dig> were highly expressed in all healthy individuals .fig.
<dig> 3d image of tra repertoires.
the numbers of tcr sequence reads bearing a given gene recombination of trav with traj were counted.
the mean percentage frequencies of  <dig>  in  <dig> healthy individuals are shown as a 3d bar graph.
x-axis and y-axis indicate trav and traj, respectively.
recombination of trav1- <dig> with traj <dig>  was the most expressed .
: pseudogene, : open reading frame
fig.
<dig> 3d image of trb repertoires.
the numbers of tcr sequence reads bearing a given gene recombination of trbv with trbj were counted.
the mean percentage frequencies of  <dig>  in  <dig> healthy individuals are shown as a 3d bar graph.
x-axis and y-axis indicate trbv and trbj, respectively.
: pseudogene, : open reading frame
digital cdr <dig> length distribution
the analysis of cdr <dig> length distribution, termed cdr <dig> size spectratyping  <cit>  or immunoscope analysis  <cit> , has been effectively used to estimate the diversity of the tcr repertoire.
this technique is based on actual peak distributions of pcr amplicons including cdr <dig> sequences by gel electrophoresis.
in the current study, the length of the determined nucleotide sequences of tcrs ranging from conserved cys <dig>  to conserved phenylalanine at position  <dig>  were calculated automatically.
this provided a visually easy way to estimate the diversity and clonality of tcrs by using ngs data.
rg can produce figures depicting the digital cdr <dig> length distribution for each v region.
cdr <dig> length distribution of both tra and trb were similar to a normal distribution but were not completely symmetrical .
cdr <dig> length was shorter in tra than in trb  and tra had more positive skewness than trb , indicating the distribution in tra was concentrated to the left.
in addition, tra had more positive kurtosis than in trb, showing a high degree of peakedness in tra .fig.
<dig> digital cdr <dig> length distribution of tra and trb.
lengths of cdr <dig> were determined in  <dig>  tra and  <dig>  trb sequence reads obtained from the pooled data of  <dig> individuals.
length of nucleotide sequences from conserved cysteine at position  <dig>  of imgt nomenclature to conserved phenylalanine at position  <dig>  were automatically calculated using rg software.
distribution of cdr <dig> length in tra  and trb  is shown as a histogram
diversity of tra and trb repertoires
to demonstrate diversity of the tcr repertoire, we calculated the mean copy number of usr and diversity indices such as the simpson index and shannon-weaver index .
the mean copy number of usr was not significantly different between tra and trb .
in addition, there was no significant difference in simpson inverse index  and shannon-weaver index  between tra and trb .
these results indicated that immune diversity for tcrα and β in healthy individuals was not different.
next, we examined whether the diversity was different between in-frame  read and out-of-frame  reads .
the result indicated that shannon diversity indices were significantly higher in in-frame reads than in out-of-frame reads .
however, there was no difference in the inverse simpson index between in-frame and out-of-frame reads.
furthermore, we examined if t cell diversity was correlated with age .
there was a significant correlation of shannon index in tra with age with a spearman’s correlation of − <dig>   but no significant correlation of the other indices with age.
overall, these results showed tendencies for inverse relationship of tcr diversity with age.fig.
<dig> diversity of tra and trb repertoires in healthy individuals.
copy number  of unique sequence reads  was calculated.
mean copy number per unique sequence read in each individual are shown as open circle .
inverse simpson index  and shannon-weaver index  were calculated with r program according to formulas described in materials and methods.
each open circle indicates the index of an individual.
there were no significant differences in mean copy number, inverse simpson index and shannon-weaver index between tra and trb
similarity of tra and trb repertoire among healthy individuals
to reveal correlations of gene usage between individuals, percentage frequencies of each trv and trj were plotted between all pairs of individuals by scatterplot .
spearman’s correlation coefficient between each pair was calculated.
the concordance correlation coefficients were lower in trav than in trbv  and the values were lower in traj than in trbj  .
these results indicated that the expression levels of trv and trj between healthy individuals were more similar among individuals in trb compared with tra.
to evaluate the potential similarity of tcr repertoires at the clonal level between healthy individuals, we retrieved tcr sequence reads that were shared among individuals.
the number of shared tcr reads between all pairs of individuals was counted and their occurrence frequencies were calculated .
the mean frequency was significantly higher in tra compared with trb  , indicating the tra repertoire contains more common tcr reads among individuals than trb does.
a similarity index, morisita-horn index, was significantly larger for tra than trb .
these results clearly indicated that tra repertoires were more similar between healthy individuals compared with trb repertoires.fig.
<dig> similarity of tra and trb repertoires in healthy individuals.
the occurrence frequency of tcr sequence reads shared between all pairs of  <dig> individuals was calculated .
mean percentage frequency of shared reads were compared between tra and trb .
similarity index, morisita-horn index, was calculated with the r program according to a formula described in materials and methods.
there were significant differences in the frequency of shared reads and similarity index between tra and trb
shared tcr sequences among healthy individuals
a small number of tcr sequences were shared between different healthy individuals.
in contrast, most tcr sequences were specific for each healthy individual.
to identify shared tcr sequences in  <dig> healthy individuals, we retrieved tra and trb reads shared between two or more healthy individuals.
in  <dig> individuals,  <dig> shared tra and  <dig> shared trb sequences were obtained from  <dig>  and  <dig>  usrs, respectively .
shared tra were more frequent in pbls from healthy individuals than trb were.
shared trb sequences were obtained from two to four individuals while the shared tra sequences were observed in  <dig> individuals.
these results indicated that shared tra sequences were more commonly used by individuals but that the trb repertoire was more specific for each individual.
furthermore, the occurrence frequencies per individual of tcr sequences shared between a pair of individuals were significantly higher for tra  than trb .
to characterize the shared tra sequences, we compared the cdr <dig> length between shared and unshared tra sequences and observed that shared tra had a shorter cdr <dig> length than unshared tra  .table  <dig> numbers of tra and trb sequences shared among multiple healthy individuals
numbers of identical tcr sequences observed in multiple healthy individuals  were counted
fig.
<dig> shared tcr had shorter cdr <dig> lengths than private tcr.
cdr <dig> lengths were calculated with  <dig> usrs of shared tcr  and  <dig>  usrs of unshared tcr .
percentage frequencies of usrs in each cdr <dig> length are plotted as a bar graph.
median values of cdr <dig> length in shared and unshared tcr are  <dig> and  <dig>  respectively
tcrs shared among multiple individuals frequently contain invariant tcrα chains
shared tra were frequently observed in pbls from healthy individuals.
to determine the origin of the shared tra, we investigated cdr <dig> sequences of shared tra reported previously.
interestingly, shared tra sequences that were common to multiple individuals contained a high proportion of tcrα related with invariant tcrα indicative of inkt cells or mait cells .
it was reported that mait cells express trav1- <dig> and traj <dig> while inkt express trav <dig> and traj <dig>  many shared tras used trav1- <dig> and traj <dig> with different cdr <dig> sequences.
the total percentage frequencies of mait tras bearing trav1- <dig> and traj <dig> and inkt tras bearing trav <dig> and traj <dig> were  <dig>  ±  <dig>  % and  <dig>  ±  <dig>  % per individual, respectively.
of  <dig> shared tra sequences,  <dig>  mait and one  inkt sequences were observed in six or more individuals .
the rate increased with the number of overlapping individuals.
germline-like cdr <dig> sequences that had no amino acid sequences altered from germline sequences were observed in  <dig>  of  <dig> shared tras except for mait  and nkt .table  <dig> invariant tcrs observed in shared tra sequences
shared individualsa
anumber of individuals in which respective shared tras were detected; bnon-germline amino acid sequences generated by nucleotide addition/deletion are underlined; cmait-like: mucosal-associated invariant t cells-like
discussion
high-throughput sequencing technologies have taken a great leap forward with the development of a wide variety of ngs platforms.
ngs facilitates the acquisition of an enormous amount of sequence data but still requires pcr amplification or gene enrichment to sequence genes of specific interest instead of the entire genome or gene library.
for heterogeneous tcr or bcr genes generated by rearrangement of many gene segments, multiplex pcr with many sets of gene-specific primers have been widely used.
however, the use of multiple primers causes amplification bias between respective genes, hampering the accurate estimation of gene frequency.
here, we used an unbiased pcr technique, an adaptor-ligation mediated pcr, for ngs-based tcr repertoire analysis.
the method uses a single set of primers to avoid pcr bias by competition between primers.
therefore, it is better suited to estimate accurately the abundances of respective tcr genes from a wide variety of samples.
we comprehensively examined tra and trb repertoires at the clonal level from a large number of individuals  and evaluated a large number amount of sequence data .
thus, this study precisely revealed the normal range of gene usage as well the extent of diversity and similarity of tcr repertoires in healthy individuals.
compared with the illumina ngs platform  <cit> , sample sequence reads were less numerous but were longer and of higher quality.
using the illumina platform, a different sequence depth among cdr <dig> contig generated from many shotgun reads may make it difficult to determine the frequency of tcr clonotypes.
however, all tcr sequences were determined from a single read and had long sequences that covered the entire region of cdr <dig>  v and j .
direct analysis from read sequences without assembly is likely to reflect accurately the actual frequencies of tcr clonotypes.
error rates in tcr sequences were slightly less than a previous report showing a mean error rate for 454-sequences of  <dig>  %  <cit> , suggesting high levels of accuracy and quality irrespective of nested pcr.
homopolymeric stretches within coding regions occur typically with the 454-sequence methodology, causing a frame shift in coding sequence.
this leads to higher rate of production of out-of-frame reads in the 454-sequence compared with other sequence platforms.
bolotin et al.
has previously reported that higher percentage of mismatch-containing sequencing in the illumina than in the roche  <dig> and ion torrent datasets   <cit> .
the error rate obtained in this study seems to be relatively lower than that in the previous report, even though our data showed higher error rate in out-of-frame than in-frame reads.
this supports that the error rates obtained are lower than in the illumina although our results did not provide a direct evidence.
furthermore, the assignment and aggregation software, rg, can rapidly summarize usage as well as recombination usage of trv and trj.
this integrated analysis easily allows the detection of preferential usage of a given trv and/or trj and therefore it will be useful for studying immune responses by antigen-specific t cells.
unlike widely used multiplex pcr methods that typically require compensation for pcr bias  <cit> , the al-pcr method is supposed to accurately estimate tcr repertoires without the compensation.
high expression levels of trbv <dig> , trbv <dig>  and trbv7- <dig>  as well as low expression levels of trbv20- <dig> , trbv <dig>  and trbv29- <dig>  were reported in cd4+ and cd8+ cells by multiplex pcr  <cit> .
however, flow cytometry analysis showed that trbv <dig> and trbv <dig> were abundantly expressed in pbls  <cit> .
to examine difference of accuracy between al-pcr and multiplex pcr in detail, we compared usage of trbv obtained with either al-pcr or multiplex pcr with facs data reported previously by van den beemd et al.
<cit> .
the result indicated that al-pcr method was better correlated with facs method than multiplex pcr method, suggesting that the al-pcr method with a set of universe primers enables us to accurately analyze tcr repertoires.
in addition, our results of tcr repertoires are similar to a previous report  <cit> .
therefore, this method will provide direct, accurate and dependable results of tcr repertoires.
by comparison to large number of healthy individuals, it has been disclosed that disease patients with x-linked agammaglobulinemia  <cit>  or common variable immune deficiency   <cit>  had skewed and contracted tcr repertoires.
it is important to clarify tcr repertoires of healthy individuals in considering how much disease patients differ from normal.
regarding usages of trv and trj repertoires, we observed preferential usages of trv and trj in peripheral bloods from healthy individuals.
similar usages between in-frame  and out-of-frame  reads suggests that the preferential usages are unlikely due to peripheral selection.
given the preferential usage was observed in immature t cells  <cit> , this is likely to be influenced by genetic factors such as recombination process.
recombination usage exhibited infrequent recombinations of aj-proximal 3’ av segment to av-distal 3’ aj segment and aj-distal 5’ av segment to av-proximal 5’ aj segments.
in gene rearrangement of the tcrαδ locus, activation of the tcrα enhancer  and the t early activation  promoter initiate primary rearrangement of proximal trav and traj segments.
subsequent secondary rearrangement occurs using 5’ distal trav and distal 3’ traj genes , resulting in the restricted usage of tra repertoires   <cit> .
however, all trav genes can recombine with traj genes in secondary rearrangement by the model of locus contraction and dna looping formation  <cit> .
although there was inefficient recombination of distal-proximal and proximal-distal trav-traj genes, traj usage was not limited over all trav but rather was equally distributed.
this suggests that the frequency of recombination varies dependent on the location of trav and probably depends on the ability of loop formation between trav and traj loci.
potential tcr diversity generated by recombination and nucleotide addition/deletion has been estimated to be up to  <dig>  <cit> .
by ngs-based estimation, trb diversity was estimated to be 3–4 ×  <dig>  <cit>  or approximately 1 ×  <dig> in humans  <cit> .
furthermore, diversity of tra is 50 % of that of trb in humans  <cit> .
in mice, tra diversity was suggested to be  <dig>  ×  <dig>  <cit>  or  <dig>  ×  <dig>  <cit>  and is 10-fold lower than trb diversity.
this lower diversity of tra might be caused by a difference in recombination processes between tra and trb.
however, our results showed a similar extent of diversity between tra and trb as evaluated by simpson and shannon-weaver indices.
similarly, wang et al.
reported that tcr diversity was estimated to be equal between tra and trb   <cit> .
contrary to previous reports obtained using limited number of sequences, large-scale sequencing suggests that the repertoire size for tra generated by v-j recombination is comparable with that for trb by v-d-j recombination.
as for tcr diversity, productive tcr had more diverse than unproductive one.
only a portion of t cells produce both productive and unproductive tcrs.
this difference might be depend on the number of reads obtained from the library.
also, there was a correlation between the diversity and age.
this is consistence with the previous report that age-related decrease in tcr repertoire was found  <cit> .
diverse t cells are generated from thymus and the thymic involution occurs with age.
the decrease in tcr diversity in periphery is likely due to the age-dependent decrease in thymic t cell regeneration.
of note, we found that tra repertoires were considerably similar between individuals.
this was mainly due to the presence of shared tcr sequences between two or more individuals.
it has been reported that shared tcrβ amino acid sequences have fewer additions in their nucleotide sequences  <cit> .
random nucleotide addition and deletion mediated by terminal deoxynucleotidyl transferase occurs during tcr rearrangement, resulting in a remarkable increase in diversity of the cdr <dig> region.
however, the shared tcrs appeared to have germline-like cdr <dig> sequences that did not undergo such modifications .
furthermore, the shared tcrs contained many tcr clonotypes with a shorter cdr <dig> length.
these results suggest that the frequent occurrence of shared tras is likely to be caused by a difference in the inherent recombination mechanism from trb .
it is noteworthy that the shared tra were present in a large number of individuals.
we unexpectedly found that the shared tra contained a high rate of tcrα related with invariant tcrα derived from mait cells or inkt cells.
these functionally important t cells have homogenous tcrα and diverse tcrβ.
mait cells express a canonical tcrα including trav1- <dig> -traj <dig>  and are preferentially localized in the gut lamina propria  <cit>  and a tcrα bearing trav1-2-traj <dig> and trav1-2-traj <dig>  <cit> .
mait cells recognize vitamin b <dig> metabolites presented by mr <dig>  non-classical mhc class i molecule  <cit> .
furthermore, cd1d-restricted inkt cells express an invariant trav <dig> -traj <dig>  chain and semi-invariant trbv25- <dig>   <cit>  and recognize glycolipids such as α-galactosyl ceramide, self-glycolipid, or isoglobotrihexosyl ceramide  <cit> .
both cell types play an essential role in the regulation of immune responses against infections, tumors, autoimmune diseases, and tolerance induction  <cit> .
frequencies of mait and inkt cells obtained in this study were consistent with previous reports showing mait cells expanded up to 1–4 % of peripheral blood t cells  <cit>  and that inkt cells accounted for  <dig>  % of total pbmcs  <cit> .
interestingly, different types of shared sequences bearing trav1- <dig>  and several shared tra sequences other than the well-known mait and inkt sequences exist.
therefore, ngs-based repertoire analysis is useful for both estimating the frequency of mait or inkt cells as well as identifying potential new invariant tcrα chains.
further identification and verification is required to identify potential novel invariant tcrα.
CONCLUSIONS
we developed a new ngs-based tcr repertoire analysis method and thereby clearly revealed comparable diversity and different interindividual similarity between tra and trb.
shared tra sequences contained frequent functionally significant t cell subpopulations, mait and inkt cells.
the approach to seeking shared tcr by ngs would be useful for identification of potential new invariant tcrα chains.
this useful technology for tcr repertoire analysis will enable us to reveal antigen-specific t cells relevant to the pathogenesis of human disease and contribute to studies of innate and adaptive immunity.
additional file
additional file 1:
table s <dig>  age, gender and chronic illness of  <dig> healthy individuals.
table s <dig>  numbers of unique reads, reads and nucleotides in tra reads obtained from pbmcs of  <dig> healthy individuals.
table s <dig>  numbers of unique reads, reads and nucleotides in trb reads obtained from pbmcs of  <dig> healthy individuals.
table s <dig>  percentage of mismatched nucleotides in in-frame and out-of-frame tcr sequences.
table s <dig>  occurrence frequency of out-of-frame unique sequence reads in tra and trb.
table s <dig>  percentage frequency of shared tra reads between all pairs of individuals.
table s <dig>  percentage frequency of shared trb reads between all pairs of individuals.
figure s <dig>  correlation of gene usage of trav, traj, trbv and trbj between healthy individuals.
figure s <dig>  concordance correlation coefficient in trav, traj, trbv and trbj.
figure s <dig>  comparison of tcr usages between in-frame and out-of-frame reads sequences.
figure s <dig>  diversity indices of in-frame and out-of-frame tra and trb.
figure s <dig>  correlation of tcr diversity with age.
figure s <dig>  correlation of tcr usage from a published facs data with al-pcr and multiplex pcr.
abbreviations
cdrcomplementarity determining region
trat cell receptor alpha
trajt cell receptor alpha joining
travt cell receptor alpha variable
trbt cell receptor beta
trbjt cell receptor alpha joining
trbvt cell receptor beta variable

 the term t cell repertoire describes a collection of lymphocytes characterized by t cell receptor expression, which plays a critical role in antigen recognition. since alterations of the t cell repertoire provide a significant indication of immune status in physiological and disease conditions, t cell repertoire analyses have been conducted for the identification of antigen-specific t cells involved in the development of disease and for the diagnosis of t lymphocyte abnormalities. comparison of variable-region usage by fluorescence-activated cell sorter analysis using a large panel of antibodies specific for tcr variable regions, polymerase chain reaction with multiple primers or pcr-based enzyme-linked immunosorbent assay have been widely used to detect changes in t cell repertoire. length distribution analysis known as cdr spectratyping is based on the addition of non-template nucleotides in v--j region and has been used to evaluate t cell clonality and diversity. to identify the antigen specificity of t cells further, pcr cloning of tcr clonotypes and subsequent sequence determination of the antigen recognition region, cdr have been required. these conventional approaches are commonly used but are time-consuming and a laborious way to study tcr repertoires. in recent years, advances in high-throughput sequencing technologies known as next-generation sequencing have rapidly progressed and enabled large-scale analysis of sequence data. although several ngs-based tcr repertoire analysis systems have been developed by other researches, many amplification techniques are based on multiple pcr with different primers specific for each variable region. thus, bias during pcr amplification is unavoidable since bias is most commonly due to differential hybridization kinetics among variable region-specific primers to different target genes. correction and additional computational normalization methods are therefore required to minimize pcr bias when using multiple pcr assays. the use of a single set of primers is a better way to achieve unbiased and quantitative amplification of all tcr genes including unknown variants where the 5’ ends of sequences are highly diverse. a single strand oligonucleotide anchor ligation to the 3’ end of cdna with t rna ligase, homopolymeric tailing of cdna, 5’ rapid amplification of cdna ends and template switching pcr have been used to analyze tcr repertoires. ts-pcr is simple and convenient but produces high levels of background amplification because ts primers non-specifically anneal to random regions in rna or allow the repeated addition of ts primers. thus, the current study describes an adaptor-ligation mediated pcr developed by the addition of an adaptor to the 5’ end of double stranded cdna from tcr transcripts and subsequent pcr amplification with the adaptor primer and constant region-specific primer, as first reported by tsuruta et al.. the adaptor ligation to blunt-ended ds cdna is less influenced by the sequence of a particular cdna while the efficiency of 5’ adaptor ligation with t rna ligase is sequence dependent. in addition, the ligation of dsdna by t ligase is more efficient than ssdna ligation with t rna ligase in ligation anchored pcr. various sequencing technologies such as roche, illumina, ion-torrent, solid, helicos and pacbio have been developed. among these ngs platforms, the dna sequencing produces sequence reads ranging from to base pairs or more in length and sufficient read outputs, yet less reads per run than the illumina. long read sequencing allows determination of the full or near-complete length of tcr genes including v, d, j and c regions. furthermore, recombinant tcr proteins can be easily produced by subsequent pcr cloning of the tcr genes. therefore, we applied an adaptor-ligation mediated pcr method to ngs with dna sequencing. natural killer t cells are a distinct t cell population with an important role in innate and adaptive immunity. nkt cells regulate a broad range of immune responses such as autoimmune diseases, tumor surveillance, and host defense against pathogenic infections. nkt cells express an invariant tcrα consisting of vα and jα that recognizes glycolipids presented by a non-classical major histocompatibility complex class i-related protein, cd1d. recently, mucosal-associated invariant t cells, which preferentially exist in mucosal tissues, were shown to be a unique t cell population expressing a semi-invariant tcrα consisting of vα and jα mait cells recognize microbial vitamin b metabolites presented by a non-classical mhc class i molecule, mhc-related protein. these t cell populations bearing invariant tcrα play a pivotal role in immune regulation but it remains to be determined whether all invariant tcrα are expressed by these unique t cell populations. in this study, we conducted ngs sequencing of tcr transcripts from healthy individuals using a newly developed ngs-based tcr repertoire analysis. initially, based on sequence read count, we examined usages of variable and joining regions, and further analyzed clonality and diversity in tcrα and β genes. unique sequence reads identified using an originally developed gene analysis program were compared at a clonal level among healthy individuals. these results showed a similar usage of trv and trj and similar extent of diversity of t cells among individuals. interestingly, tcrβ reads were less shared among individuals while tcrα reads frequently contained shared sequences that overlapped between two or more individuals. shared tcrα reads contained a high proportion of invariant tcrα indicating the presence of inkt cells or mait cells. in this report, we demonstrated that analysis of tcr genes shared among multiple individuals from ngs data provided significant information on invariant tcrs expressed by nkt cells and mait cells. methods isolation of peripheral blood mononuclear cells and rna extraction after obtaining written informed consent, whole blood was collected from japanese healthy individuals. the study was approved by the ethics committees of the clinical research center for rheumatology and allergy, sagamihara national hospital, national hospital organization. ten ml of whole blood was collected into heparinized tubes. peripheral blood mononuclear cells were isolated with ficoll-paque plus™ density gradient centrifugation and washed with phosphate buffered saline. cell numbers were counted and 1 ×  cells were used for rna extraction. total rna was isolated and purified with rneasy lipid tissue mini kit according to the manufacturer’s instructions. rna amounts and purity were measured with agilent bioanalyzer. unbiased amplification of tcr genes one microgram of total rna was converted to complementary dna with superscript iii reverse transcriptase. bsl-18e primer containing polyt and a noti site was used for cdna synthesis. after cdna synthesis, double strand -cdna was synthesized with e. coli dna polymerase i, e. coli dna ligase, and rnase h. ds-cdnas were blunted with t dna polymerase. p10ea/p20ea adaptor was ligated to the 5’ end of the ds-cdna and then cut with noti restriction enzyme. after removal of adaptor and primer with minelute reaction cleanup kit, pcr was performed using either tcr α-chain constant region-specific or tcr β-chain constant region-specific primers and p20ea. pcr conditions were as follows: 95 °c, 55 °c, and 72 °c for 20 cycles. the second pcr was performed with either ca or cb and p20ea primers using the same pcr conditions.table primers used in this study ccatctcatccctgcgtgtctccgac ccatctcatccctgcgtgtctccgac v: a/c/g, n: a/c/g/t, adaptor a and b sequences were typed in bold and bold italic, respectively. a key sequence was underlined. the following mid tag sequences were used for identification of sample source. mid1: acgagtgcgt, mid2: acgctcgaca, mid3: agacgcactc, mid4: agcactgtag, mid5: atcagacacg, mid6: atatcgcgag, mid7: cgtgtctcta, mid8: ctcgcgtgtc, mid10: tctctatgcg, mid11: tgatacgtct, mid15: tacgacgta, mid16: tcacgtacta, mid17: cgtctagtac, mid18: tctacgtagc, mid19: tgtactactc, mid20: acgactacag, mid21: cgtagactag, mid22: tacgagtatg, mid23: tactctcgtg, mid24: tagagacgag amplicon sequencing by roche sequencing system amplicons for ngs were prepared by amplification of the second pcr products using p20ea primer and fusion tag primer. the fusion tag primers consisting of an a-adaptor sequence, base sequence key, multiple identifier tag sequence, and tcr constant region-specific sequence were designed according to the manufacturer’s instructions. after pcr amplification, amplicons were separated and evaluated by agarose gel electrophoresis. the resulting fragment was removed from the gel and purified with qiaex ii gel extraction kit. the amount of purified amplicon was quantified by quant-it™ picogreen® dsdna assay kit. each amplicon obtained with a different fusion tag primer from healthy individuals was mixed at equal molar concentrations. emulsion pcr was performed using the amplicon mixtures with gs junior titanium empcr lib-l kit according to the manufacturer’s instructions. assignment of trv and trj segments all sequence reads were classified by mid tag sequences. artificially added sequences and sequences with low quality scores were removed from both terminals of sequence reads using software installed on the sequencing system. the remaining sequences were used for assignment of trav and traj for tcr α sequences, and trbv and trbj for tcr β sequences. assignment of sequences was performed by determining sequences with the highest identity in a data set of reference sequences for trav, traj, trbv and trbj genes including pseudogenes and open reading frame reference sequences available from the international immunogenetics information system® database. data processing, assignment, and data aggregation were automatically performed using repertoire analysis software originally developed by our group. rg implemented a program for sequence homology searches using blatn, an automatic aggregation program, a graphics program for trv and trj usage, and cdr length distribution. sequence identities at the nucleotide level between query and entry sequences were automatically calculated. parameters that increased sensitivity and accuracy score) were carefully optimized for respective repertoire analysis. data analyses nucleotide sequences of cdr regions ranged from conserved cysteine at position of imgt nomenclature to conserved phenylalanine at position and the following glycine was translated to deduced amino acid sequences. a unique sequence read was defined as a sequence read having no identity in trv, trj and deduced amino acid sequence of cdr with the other sequence reads. the copy number of identical usr were automatically counted by rg software in each sample and then ranked in order of the copy number. percentage occurrence frequencies of sequence reads with trav, traj, trbv and trbj genes in total sequence reads were calculated. retrieval of shared usrs among samples to retrieve shared sequences among samples, a concatenate string of “trv gene name”_” deduced amino acid sequence of cdr region”_” trj gene name” of individual usr was used as a tcr identifier. the tcr identifier in a sample was retrieved in read data sets from all the other samples. diversity indices and similarity index to estimate tcr diversity in deep sequence data, several diversity indices, simpson’s index and shannon-weaver index were calculated using a function “diversity” of the vegan package in the r program. these indices were calculated based on the number of species per sample and the number of individuals per sample as measures for biological diversity in ecology. in deep sequence data, usr and copy number were used for species and individuals, respectively. simpson’s index was defined as: 1−λ=1−∑i=1snini−1nn− where n is the total number of sequence reads, ni is the copy number of ith usr, and s is the species number of usr. this value ranges from to where the maximum number means high levels of diversity and indicates low diversity. simpson’s reciprocal index was also calculated as the inverse of λ. the shannon-weaver index was used for the diversity index and was defined as: h'=−∑i−1sninlnnin where n is the total number of sequence reads, ni is the number of ith usr, and s is the species number of usr. these diversity indices should be biased by differences in read numbers among samples. therefore, the number of sequence reads was standardized for each sample down to the smallest number of sequence reads. to standardize sample size, repeated random resampling times without replacement and calculation of the diversity index were performed using the r program. the median of their indices was used to determine the diversity index for the sample. to estimate the similarity of tcr repertoires between healthy individuals, the morisita-horn index was calculated using the function “vegdist” in the vegan package of the r program. the morisita-horn index was defined as: ch=2∑i=1sxiyi∑i=1sxi2x2+∑i=1syi2y2xy where xi is the number of ith usr in the total x reads of one sample, yi is the number of ith usr in the total y reads of another sample, and s is the number of usr. to standardize the sample size, repeated random resampling times without replacement and calculation of similarity index were performed using the r program. median values were used for similarity indexes between a pair of samples. statistics statistical significances were tested by the nonparametric mann–whitney u-test using graphpad prism software. a value of p <  was considered statistically significant.  repertoire analysis software the cloud-based software platform developed in this study, rg, is a high-speed, accurate and convenient computing system for tcr repertoire analysis. rg provides an integrated software package for 1) assignment of v, d and j segments, 2) calculation of sequence identity, 3) extraction of cdr sequences, 4) counting identical reads, 5) amino acid translation, 6) frame analysis, and 7) cdr length analysis. after uploading sequencing data from the ngs sequencer, v, d and j segments can be identified based on their sequence similarity with optimized parameters. read counts are automatically summarized and then processing data, summary tables, and graphs can be easily downloaded. read number, error rate, and unproductive reads we performed high-throughput sequencing of tra and trb genes in pbmcs from healthy individuals. a total of and sequence reads were assigned for tra and trb repertoire analyses, respectively using the rg program. a total of and unique sequence reads were identified in tra and trb, respectively. the number of nucleotide sequences per read obtained by roche sequencing were ~400 bp in length, indicating these sequences were long enough to identify tcr genes ranged from v region to j regions. to evaluate the accuracy and quality of ngs sequencing, we calculated the frequency of mismatched nucleotides between query and reference sequences as the error rate. error rates were  ±   % for trav,  ±   % for traj,  ±   % for trbv and  ±   % for trbj. these error rates were slightly lower than in a previous study reporting a mean error rate of  % for 454-sequences. the error rates were significantly higher in v regions than in j regions, indicating higher sequence reliability in regions close to sequencing primers. occurrence frequency of unproductive reads carrying a stop codon or a shift of reading frame in cdr regions was calculated. there was no significant difference in the percentage of occurrence frequency of unproductive unique sequence reads between tra and trb. regarding the error rates, frequencies of mismatched nucleotides were significantly higher in out-of-frame reads than in in-frame reads in all regions of trav, traj. trbv and trbj. the differences were quite significant in j regions. this suggested that sequencing errors occurred in j regions adjacent to cdr cause frequently frame-shift of coding sequence. expression of tcr genes including pseudogenes and orf to determine usages of trv and trj genes in tcr sequence reads, copy number of usrs bearing respective trv or trj were counted. individual usrs were ranked by order of copy number and the percentage frequency of respective trv and trj were calculated. regarding tra repertoires, eight pseudogenes were not expressed in healthy individuals. av8- classified as an orf, were slightly expressed. expression of av and av was not observed in healthy individuals. in addition, the functional genes, av and av9- were poorly expressed in individual and individuals, respectively. of eight aj genes classified as orf genes, the expression of aj and aj were observed in all individuals. of these, aj and aj were slightly expressed in individuals and individuals, respectively. aj aj aj and aj were not present in any individuals. there was no expression of the three pseudogenes, aj aj and aj in any individuals. the functional gene aj was detected in only reads from individuals.fig. usage of trav and traj in healthy individuals. the numbers of tcr sequences bearing respective trav and traj were counted. the percentage frequencies of trav and traj were calculated and are shown as a scatter plot. each dot indicates the percentage frequency of trav or traj in each individual. bars in red indicate the mean values of individuals.: pseudogene,: open reading frame fig. usage of trbv and trbj in healthy individuals. the percentage frequencies of trbv and trbj are shown as a scatter plot. each dot indicates the percentage frequency of trbv or trbj in each individual. bars in red indicate the mean values.: pseudogene,: open reading frame for trb genes, there was no expression of pseudogenes in healthy individuals. of five orf genes, bv5-, bv6- and bv were poorly expressed. the bv7- orf gene was not observed in any individuals while bv23- was expressed in all individuals. regarding bj genes, there was no expression of the bj2-2p pseudogene. infrequent recombination of trav and traj gene recombination of trav with traj were capable of generating a total of av-aj recombinations, of which av-aj recombinations were detected in individuals. this indicated that almost all av-aj recombinations were used in tcr transcripts without restriction. notably, av1-1–av genes were recombined less preferentially with aj50–aj genes and similarly, recombination of av35–av genes with aj3–aj was rarely observed. given the chromosomal location of these gene segments, these results showed that av-aj recombinations occurred infrequently between proximal av genes and distal aj genes and between distal av genes and proximal aj genes. for trb, gene recombinations were generated with bv and bj genes, of which bv-bj were used in individuals. there was no restriction for combinations of trbv with trbj.fig. heatmap representation of gene recombination of trav with traj in pooled read data from healthy individuals. the numbers of tcr sequence reads bearing respective gene recombination of trav and traj were counted. the occurrence tendency of recombination is visualized by heat map presentation of the number of each recombination. color in each pixel indicates the number of each recombination. for trav, pseudogenes, orf, and poorly expressed genes were excluded. for traj, pseudogenes, orfs, and poorly expressed genes were excluded. two orfs found to be expressed were included. a heat map representation of the recombination events is shown preferential usage of trv and trj repertoires in healthy individuals to reveal trv and trj usage in whole tcr transcripts, the occurrence frequency of usrs bearing respective trv or trj were calculated. preferential usage in several trav genes was observed in trav, trav9-, trav13- and trav. these preferential usages were similar to a previous result obtained using a hybridization-based quantitation assay. several trbv genes were abundantly used in the trbv repertoire. the top three trbv29-, trbv20- and trbv accounted for approximately one third of the total sequence reads. this was similar to results from our previous study using the microplate hybridization assay. gene usage varied considerably among the trbj genes. trbj2- and trbj 2– were highly expressed while the expressions of trbj1- trbj1- trbj1- trbj2- and trbj2- were low. next, to examine if the preferential usage of trv and trj are due to peripheral selection, we compared usages of tcr between in-frame and out-of-frame reads. as a whole, we observed little difference in the usages of trav, traj, trbv and trbj between in-frame and out-of-frame reads. on the other hand, particularly trav26- traj traj trbv23- trbj1- and trbj2- were more frequently used in out-of-frame reads than in-frame reads. three-dimensional view of tcr repertoire usage to visualize the usage of tcrs bearing a combination of trv with trj genes, we produced 3d pictures of the tcr repertoires. the advantage of 3d images is that the dominance of certain combination of trv with trj genes as well as the extent of diversity of tcr can easily be observed. for trb, there was little preferential usage of recombination between trbj and trbv genes. the frequency of each recombination depended on the usage of trbv or trbj. bv29-1/bj2- bv29-1/bj2- bv29-1/bj2- and bv20-1/bj2- were frequently used in all combinations while others were expressed at a low frequency. in contrast, 3d imaging of the tra repertoire indicated a low level of expression with a wide distribution of trav and traj. the occupancy was lower than 1 % for all combinations. of note, tcr reads bearing av1- and aj were highly expressed in all healthy individuals.fig. 3d image of tra repertoires. the numbers of tcr sequence reads bearing a given gene recombination of trav with traj were counted. the mean percentage frequencies of in healthy individuals are shown as a 3d bar graph. x-axis and y-axis indicate trav and traj, respectively. recombination of trav1- with traj was the most expressed.: pseudogene,: open reading frame fig. 3d image of trb repertoires. the numbers of tcr sequence reads bearing a given gene recombination of trbv with trbj were counted. the mean percentage frequencies of in healthy individuals are shown as a 3d bar graph. x-axis and y-axis indicate trbv and trbj, respectively.: pseudogene,: open reading frame digital cdr length distribution the analysis of cdr length distribution, termed cdr size spectratyping or immunoscope analysis, has been effectively used to estimate the diversity of the tcr repertoire. this technique is based on actual peak distributions of pcr amplicons including cdr sequences by gel electrophoresis. in the current study, the length of the determined nucleotide sequences of tcrs ranging from conserved cys to conserved phenylalanine at position were calculated automatically. this provided a visually easy way to estimate the diversity and clonality of tcrs by using ngs data. rg can produce figures depicting the digital cdr length distribution for each v region. cdr length distribution of both tra and trb were similar to a normal distribution but were not completely symmetrical. cdr length was shorter in tra than in trb and tra had more positive skewness than trb, indicating the distribution in tra was concentrated to the left. in addition, tra had more positive kurtosis than in trb, showing a high degree of peakedness in tra.fig. digital cdr length distribution of tra and trb. lengths of cdr were determined in tra and trb sequence reads obtained from the pooled data of individuals. length of nucleotide sequences from conserved cysteine at position of imgt nomenclature to conserved phenylalanine at position were automatically calculated using rg software. distribution of cdr length in tra and trb is shown as a histogram diversity of tra and trb repertoires to demonstrate diversity of the tcr repertoire, we calculated the mean copy number of usr and diversity indices such as the simpson index and shannon-weaver index. the mean copy number of usr was not significantly different between tra and trb. in addition, there was no significant difference in simpson inverse index and shannon-weaver index between tra and trb. these results indicated that immune diversity for tcrα and β in healthy individuals was not different. next, we examined whether the diversity was different between in-frame read and out-of-frame reads. the result indicated that shannon diversity indices were significantly higher in in-frame reads than in out-of-frame reads. however, there was no difference in the inverse simpson index between in-frame and out-of-frame reads. furthermore, we examined if t cell diversity was correlated with age. there was a significant correlation of shannon index in tra with age with a spearman’s correlation of − but no significant correlation of the other indices with age. overall, these results showed tendencies for inverse relationship of tcr diversity with age.fig. diversity of tra and trb repertoires in healthy individuals. copy number of unique sequence reads was calculated. mean copy number per unique sequence read in each individual are shown as open circle. inverse simpson index and shannon-weaver index were calculated with r program according to formulas described in materials and methods. each open circle indicates the index of an individual. there were no significant differences in mean copy number, inverse simpson index and shannon-weaver index between tra and trb similarity of tra and trb repertoire among healthy individuals to reveal correlations of gene usage between individuals, percentage frequencies of each trv and trj were plotted between all pairs of individuals by scatterplot. spearman’s correlation coefficient between each pair was calculated. the concordance correlation coefficients were lower in trav than in trbv and the values were lower in traj than in trbj. these results indicated that the expression levels of trv and trj between healthy individuals were more similar among individuals in trb compared with tra. to evaluate the potential similarity of tcr repertoires at the clonal level between healthy individuals, we retrieved tcr sequence reads that were shared among individuals. the number of shared tcr reads between all pairs of individuals was counted and their occurrence frequencies were calculated. the mean frequency was significantly higher in tra compared with trb, indicating the tra repertoire contains more common tcr reads among individuals than trb does. a similarity index, morisita-horn index, was significantly larger for tra than trb. these results clearly indicated that tra repertoires were more similar between healthy individuals compared with trb repertoires.fig. similarity of tra and trb repertoires in healthy individuals. the occurrence frequency of tcr sequence reads shared between all pairs of individuals was calculated. mean percentage frequency of shared reads were compared between tra and trb. similarity index, morisita-horn index, was calculated with the r program according to a formula described in materials and methods. there were significant differences in the frequency of shared reads and similarity index between tra and trb shared tcr sequences among healthy individuals a small number of tcr sequences were shared between different healthy individuals. in contrast, most tcr sequences were specific for each healthy individual. to identify shared tcr sequences in healthy individuals, we retrieved tra and trb reads shared between two or more healthy individuals. in individuals, shared tra and shared trb sequences were obtained from and usrs, respectively. shared tra were more frequent in pbls from healthy individuals than trb were. shared trb sequences were obtained from two to four individuals while the shared tra sequences were observed in individuals. these results indicated that shared tra sequences were more commonly used by individuals but that the trb repertoire was more specific for each individual. furthermore, the occurrence frequencies per individual of tcr sequences shared between a pair of individuals were significantly higher for tra than trb. to characterize the shared tra sequences, we compared the cdr length between shared and unshared tra sequences and observed that shared tra had a shorter cdr length than unshared tra.table numbers of tra and trb sequences shared among multiple healthy individuals numbers of identical tcr sequences observed in multiple healthy individuals were counted fig. shared tcr had shorter cdr lengths than private tcr. cdr lengths were calculated with usrs of shared tcr and usrs of unshared tcr. percentage frequencies of usrs in each cdr length are plotted as a bar graph. median values of cdr length in shared and unshared tcr are and respectively tcrs shared among multiple individuals frequently contain invariant tcrα chains shared tra were frequently observed in pbls from healthy individuals. to determine the origin of the shared tra, we investigated cdr sequences of shared tra reported previously. interestingly, shared tra sequences that were common to multiple individuals contained a high proportion of tcrα related with invariant tcrα indicative of inkt cells or mait cells. it was reported that mait cells express trav1- and traj while inkt express trav and traj many shared tras used trav1- and traj with different cdr sequences. the total percentage frequencies of mait tras bearing trav1- and traj and inkt tras bearing trav and traj were  ±   % and  ±   % per individual, respectively. of shared tra sequences, mait and one inkt sequences were observed in six or more individuals. the rate increased with the number of overlapping individuals. germline-like cdr sequences that had no amino acid sequences altered from germline sequences were observed in of shared tras except for mait and nkt.table invariant tcrs observed in shared tra sequences shared individualsa anumber of individuals in which respective shared tras were detected; bnon-germline amino acid sequences generated by nucleotide addition/deletion are underlined; cmait-like: mucosal-associated invariant t cells-like discussion high-throughput sequencing technologies have taken a great leap forward with the development of a wide variety of ngs platforms. ngs facilitates the acquisition of an enormous amount of sequence data but still requires pcr amplification or gene enrichment to sequence genes of specific interest instead of the entire genome or gene library. for heterogeneous tcr or bcr genes generated by rearrangement of many gene segments, multiplex pcr with many sets of gene-specific primers have been widely used. however, the use of multiple primers causes amplification bias between respective genes, hampering the accurate estimation of gene frequency. here, we used an unbiased pcr technique, an adaptor-ligation mediated pcr, for ngs-based tcr repertoire analysis. the method uses a single set of primers to avoid pcr bias by competition between primers. therefore, it is better suited to estimate accurately the abundances of respective tcr genes from a wide variety of samples. we comprehensively examined tra and trb repertoires at the clonal level from a large number of individuals and evaluated a large number amount of sequence data. thus, this study precisely revealed the normal range of gene usage as well the extent of diversity and similarity of tcr repertoires in healthy individuals. compared with the illumina ngs platform, sample sequence reads were less numerous but were longer and of higher quality. using the illumina platform, a different sequence depth among cdr contig generated from many shotgun reads may make it difficult to determine the frequency of tcr clonotypes. however, all tcr sequences were determined from a single read and had long sequences that covered the entire region of cdr v and j. direct analysis from read sequences without assembly is likely to reflect accurately the actual frequencies of tcr clonotypes. error rates in tcr sequences were slightly less than a previous report showing a mean error rate for 454-sequences of  %, suggesting high levels of accuracy and quality irrespective of nested pcr. homopolymeric stretches within coding regions occur typically with the 454-sequence methodology, causing a frame shift in coding sequence. this leads to higher rate of production of out-of-frame reads in the 454-sequence compared with other sequence platforms. bolotin et al. has previously reported that higher percentage of mismatch-containing sequencing in the illumina than in the roche and ion torrent datasets. the error rate obtained in this study seems to be relatively lower than that in the previous report, even though our data showed higher error rate in out-of-frame than in-frame reads. this supports that the error rates obtained are lower than in the illumina although our results did not provide a direct evidence. furthermore, the assignment and aggregation software, rg, can rapidly summarize usage as well as recombination usage of trv and trj. this integrated analysis easily allows the detection of preferential usage of a given trv and/or trj and therefore it will be useful for studying immune responses by antigen-specific t cells. unlike widely used multiplex pcr methods that typically require compensation for pcr bias, the al-pcr method is supposed to accurately estimate tcr repertoires without the compensation. high expression levels of trbv, trbv and trbv7- as well as low expression levels of trbv20-, trbv and trbv29- were reported in cd4+ and cd8+ cells by multiplex pcr. however, flow cytometry analysis showed that trbv and trbv were abundantly expressed in pbls. to examine difference of accuracy between al-pcr and multiplex pcr in detail, we compared usage of trbv obtained with either al-pcr or multiplex pcr with facs data reported previously by van den beemd et al.. the result indicated that al-pcr method was better correlated with facs method than multiplex pcr method, suggesting that the al-pcr method with a set of universe primers enables us to accurately analyze tcr repertoires. in addition, our results of tcr repertoires are similar to a previous report. therefore, this method will provide direct, accurate and dependable results of tcr repertoires. by comparison to large number of healthy individuals, it has been disclosed that disease patients with x-linked agammaglobulinemia or common variable immune deficiency had skewed and contracted tcr repertoires. it is important to clarify tcr repertoires of healthy individuals in considering how much disease patients differ from normal. regarding usages of trv and trj repertoires, we observed preferential usages of trv and trj in peripheral bloods from healthy individuals. similar usages between in-frame and out-of-frame reads suggests that the preferential usages are unlikely due to peripheral selection. given the preferential usage was observed in immature t cells, this is likely to be influenced by genetic factors such as recombination process. recombination usage exhibited infrequent recombinations of aj-proximal 3’ av segment to av-distal 3’ aj segment and aj-distal 5’ av segment to av-proximal 5’ aj segments. in gene rearrangement of the tcrαδ locus, activation of the tcrα enhancer and the t early activation promoter initiate primary rearrangement of proximal trav and traj segments. subsequent secondary rearrangement occurs using 5’ distal trav and distal 3’ traj genes, resulting in the restricted usage of tra repertoires. however, all trav genes can recombine with traj genes in secondary rearrangement by the model of locus contraction and dna looping formation. although there was inefficient recombination of distal-proximal and proximal-distal trav-traj genes, traj usage was not limited over all trav but rather was equally distributed. this suggests that the frequency of recombination varies dependent on the location of trav and probably depends on the ability of loop formation between trav and traj loci. potential tcr diversity generated by recombination and nucleotide addition/deletion has been estimated to be up to. by ngs-based estimation, trb diversity was estimated to be 3–4 ×  or approximately 1 ×  in humans. furthermore, diversity of tra is 50 % of that of trb in humans. in mice, tra diversity was suggested to be  ×  or  ×  and is 10-fold lower than trb diversity. this lower diversity of tra might be caused by a difference in recombination processes between tra and trb. however, our results showed a similar extent of diversity between tra and trb as evaluated by simpson and shannon-weaver indices. similarly, wang et al. reported that tcr diversity was estimated to be equal between tra and trb. contrary to previous reports obtained using limited number of sequences, large-scale sequencing suggests that the repertoire size for tra generated by v-j recombination is comparable with that for trb by v-d-j recombination. as for tcr diversity, productive tcr had more diverse than unproductive one. only a portion of t cells produce both productive and unproductive tcrs. this difference might be depend on the number of reads obtained from the library. also, there was a correlation between the diversity and age. this is consistence with the previous report that age-related decrease in tcr repertoire was found. diverse t cells are generated from thymus and the thymic involution occurs with age. the decrease in tcr diversity in periphery is likely due to the age-dependent decrease in thymic t cell regeneration. of note, we found that tra repertoires were considerably similar between individuals. this was mainly due to the presence of shared tcr sequences between two or more individuals. it has been reported that shared tcrβ amino acid sequences have fewer additions in their nucleotide sequences. random nucleotide addition and deletion mediated by terminal deoxynucleotidyl transferase occurs during tcr rearrangement, resulting in a remarkable increase in diversity of the cdr region. however, the shared tcrs appeared to have germline-like cdr sequences that did not undergo such modifications. furthermore, the shared tcrs contained many tcr clonotypes with a shorter cdr length. these results suggest that the frequent occurrence of shared tras is likely to be caused by a difference in the inherent recombination mechanism from trb. it is noteworthy that the shared tra were present in a large number of individuals. we unexpectedly found that the shared tra contained a high rate of tcrα related with invariant tcrα derived from mait cells or inkt cells. these functionally important t cells have homogenous tcrα and diverse tcrβ. mait cells express a canonical tcrα including trav1- -traj and are preferentially localized in the gut lamina propria and a tcrα bearing trav1-2-traj and trav1-2-traj. mait cells recognize vitamin b metabolites presented by mr non-classical mhc class i molecule. furthermore, cd1d-restricted inkt cells express an invariant trav -traj chain and semi-invariant trbv25- and recognize glycolipids such as α-galactosyl ceramide, self-glycolipid, or isoglobotrihexosyl ceramide. both cell types play an essential role in the regulation of immune responses against infections, tumors, autoimmune diseases, and tolerance induction. frequencies of mait and inkt cells obtained in this study were consistent with previous reports showing mait cells expanded up to 1–4 % of peripheral blood t cells and that inkt cells accounted for  % of total pbmcs. interestingly, different types of shared sequences bearing trav1- and several shared tra sequences other than the well-known mait and inkt sequences exist. therefore, ngs-based repertoire analysis is useful for both estimating the frequency of mait or inkt cells as well as identifying potential new invariant tcrα chains. further identification and verification is required to identify potential novel invariant tcrα.  we developed a new ngs-based tcr repertoire analysis method and thereby clearly revealed comparable diversity and different interindividual similarity between tra and trb. shared tra sequences contained frequent functionally significant t cell subpopulations, mait and inkt cells. the approach to seeking shared tcr by ngs would be useful for identification of potential new invariant tcrα chains. this useful technology for tcr repertoire analysis will enable us to reveal antigen-specific t cells relevant to the pathogenesis of human disease and contribute to studies of innate and adaptive immunity. additional file additional file 1: table s age, gender and chronic illness of healthy individuals. table s numbers of unique reads, reads and nucleotides in tra reads obtained from pbmcs of healthy individuals. table s numbers of unique reads, reads and nucleotides in trb reads obtained from pbmcs of healthy individuals. table s percentage of mismatched nucleotides in in-frame and out-of-frame tcr sequences. table s occurrence frequency of out-of-frame unique sequence reads in tra and trb. table s percentage frequency of shared tra reads between all pairs of individuals. table s percentage frequency of shared trb reads between all pairs of individuals. figure s correlation of gene usage of trav, traj, trbv and trbj between healthy individuals. figure s concordance correlation coefficient in trav, traj, trbv and trbj. figure s comparison of tcr usages between in-frame and out-of-frame reads sequences. figure s diversity indices of in-frame and out-of-frame tra and trb. figure s correlation of tcr diversity with age. figure s correlation of tcr usage from a published facs data with al-pcr and multiplex pcr. abbreviations cdrcomplementarity determining region trat cell receptor alpha trajt cell receptor alpha joining travt cell receptor alpha variable trbt cell receptor beta trbjt cell receptor alpha joining trbvt cell receptor beta variable 
to compare gene usage, diversity and similarity of expressed tcr repertoires among individuals, we conducted next-generation sequencing of tra and trb genes in peripheral blood mononuclear cells from healthy human individuals. preferential usage of several v and j genes were observed while some recombinations of trav with traj appeared to be restricted. here, we developed a new tcr repertoire analysis method using dna sequencing technology in combination with an adaptor-ligation mediated polymerase chain reaction. a publicly available tra had a near-germline tcr with a shorter cdr notably, shared tra sequences, especially those shared among a large number of individuals’, often contained tcrα related with invariant tcrα derived from invariant natural killer t cells and mucosal-associated invariant t cells. the extent of tcr diversity was not significantly different between tra and trb, while tra repertoires were more similar between individuals than trb repertoires were. high-throughput sequencing of t cell receptor genes is a powerful tool for analyses of antigen specificity, clonality and diversity of t lymphocytes. these results suggest that retrieval of shared tcrs by ngs would be useful for the identification of potential new invariant tcrα chains. this method allows the amplification of all tcr genes without pcr bias. t cell receptorrepertoirenext generation sequencingimmune profilinginvariant tcrαissue-copyright-statement© the author 2016 from a total of sequence reads from individuals, unique sequence reads were identified. this ngs method will enable the comprehensive quantitative analysis of tcr repertoires at a clonal level.
 high-throughput sequencing of t cell receptor genes is a powerful tool for analyses of antigen specificity, clonality and diversity of t lymphocytes. here, we developed a new tcr repertoire analysis method using dna sequencing technology in combination with an adaptor-ligation mediated polymerase chain reaction. this method allows the amplification of all tcr genes without pcr bias. to compare gene usage, diversity and similarity of expressed tcr repertoires among individuals, we conducted next-generation sequencing of tra and trb genes in peripheral blood mononuclear cells from healthy human individuals.  from a total of sequence reads from individuals, unique sequence reads were identified. preferential usage of several v and j genes were observed while some recombinations of trav with traj appeared to be restricted. the extent of tcr diversity was not significantly different between tra and trb, while tra repertoires were more similar between individuals than trb repertoires were. the interindividual similarity of tra depended largely on the frequent presence of shared tcrs among two or more individuals. a publicly available tra had a near-germline tcr with a shorter cdr notably, shared tra sequences, especially those shared among a large number of individuals’, often contained tcrα related with invariant tcrα derived from invariant natural killer t cells and mucosal-associated invariant t cells.  these results suggest that retrieval of shared tcrs by ngs would be useful for the identification of potential new invariant tcrα chains. this ngs method will enable the comprehensive quantitative analysis of tcr repertoires at a clonal level. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords t cell receptorrepertoirenext generation sequencingimmune profilinginvariant tcrαissue-copyright-statement© the author 2016 