| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
|
|
|
|||||||
|
|||||||||
ANIMAL GENETICS |
,3
* Norwegian Institute of Fisheries and Aquaculture Research, N-9291, Tromsø, Norway; and
and
Department of Animal Science, Iowa State University, Ames 50011-3150
 |
Abstract |
|---|
 Top Abstract INTRODUCTIONMATERIALS AND METHODSRESULTS AND DISCUSSIONIMPLICATIONSLITERATURE CITED |
|---|
Key Words: Atlantic cod • challenge test • Cox regression • disease resistance • heritability • Vibrio anguillarum
|
INTRODUCTION |
|---|
|
TopAbstractINTRODUCTIONMATERIALS AND METHODSRESULTS AND DISCUSSIONIMPLICATIONSLITERATURE CITED |
|---|
). Genetic improvement of disease resistance through selective breeding requires additive genetic variation. Considerable additive genetic variation found in disease resistance in various aquatic species (Gjøen et al., 1997; Henryon et al., 2005; Gitterle et al., 2006) gives promise that sufficient additive genetic variation may also exist in vibriosis resistance in Atlantic cod.
Heritability estimates for disease resistance have predominantly been obtained using mixed linear model approaches and considering challenge test survival as a binary trait (Gjedrem et al., 1991; Gjedrem and Gjøen, 1995; Gjøen et al., 1997). More recent studies have evaluated disease resistance as a binary and longitudinal trait (Henryon et al., 2002; Henryon et al., 2005; Gitterle et al., 2006). Partly because of differing definitions of heritability on different scales, considerable variation in the estimates of heritability was observed. Different approaches predicted highly correlated breeding values, indicating that different approaches may rank individuals identically (Henryon et al., 2005; Gitterle et al., 2006). Although survival analysis is becoming a standard approach for genetic analysis of survival traits (Ducrocq et al., 2000), the robustness of a linear model for genetic analysis of vibriosis resistance would simplify estimation of the genetic correlations between linear and survival traits.
The objective of this study was to estimate genetic variation in vibriosis resistance in Atlantic cod using 3 statistical approaches, all defining challenge test survival as a longitudinal trait. The effect of the 3 approaches, a Cox regression model (CRM), a linear model (LMO), and linear model accounting for censoring (LMC) on the sires’ predicted breeding values was evaluated.
|
MATERIALS AND METHODS |
|---|
|
TopAbstractINTRODUCTIONMATERIALS AND METHODSRESULTS AND DISCUSSIONIMPLICATIONSLITERATURE CITED |
|---|
).
Differences in PanI allele (Pogson et al., 1995) frequencies exist between northeast Arctic cod (NEAC) and Norwegian coastal cod (CC) populations (Fevolden and Pogson, 1997). The PanI was originally identified as a gene encoding for a microvesicle protein called synaptophysin (SypI; Fevolden and Pogson, 1997) but was later suggested by Pogson (2001) to represent an isoform of SypI called pantophysin. Allele PanIB was found in high frequency in samples grouped a priori in NEAC, whereas PanIA was predominant in CC samples (Fevolden and Pogson, 1997). Therefore, the brood-fish were grouped according to PanI genotypes, and only homozygous PanIAA and PanIBB genotypes were accepted in CC and NEAC brood-stock, respectively. For this study, the CC brood-stock was further divided into 2 groups, southern CC (CCS) and northern CC (CCN), depending on the origin of the fish.
The challenge test data represented a subset of the fish produced within the breeding program during the spring 2004. Data consisted of 3,576 individuals from 50 full-sib families, of which 24 (n = 1,743) were CCS, 4 were CCN (n = 309), and 22 were NEAC (n = 1,524; Table 1). Families were offspring of 39 sires and 49 dams, of which 11 sires and 1 dam had 2 families in the data set. Thus, due to the mating structure, there were some problems of confounding between sire and dam effects. The number of individuals tested per family varied from 19 to 100 (average of 72).
|
). For the challenge test, fish were randomly assigned to 1 of 3 experimental tanks. All populations and families were well represented in each of the test tanks. The proportion of individuals from different populations occupying the test tanks varied from 29.1 to 36.1%. The proportion of individuals from each family in the test tanks varied between 15.4 and 59.6%. Cod juveniles were from 248 to 305 d old when evaluated for survival after challenge with V. anguillarum at the Fish Health Unit of Tromsø Aquaculture Research Station (Tromsø, Norway).
Challenge Protocol
The V. anguillarum strain 2001/09/700-4299, isolated from wild Atlantic cod in Norway and characterised as serotype O2b (Va-O2b-4299), was provided by the Norwegian Veterinary Institute (Oslo, Norway). The bacteria were grown on blood-agar (Oxoid base No.2, Oxoid Ltd., Basingstoke, UK), supplemented with 1.5% NaCl and 2% human red blood cells, at 12°C for 2 d. Pure colonies were transferred to marine broth (MB-2216, Difco, Becton Dickinson, Sparks, MD) and incubated with shaking for 24 h at 12°C. Thereafter, a small volume of the culture was inoculated in an appropriate volume of marine broth and incubated with shaking for 24 h at 12°C.
In the prechallenge study, a dosage of 2.3 x 104 cfu/mL resulted in delayed mortality and thus was determined as the preferred challenge test dosage. Consequently, a volume of 13 mL of bacterial culture was initially mixed with 1 L of water and thereafter added into the test tank containing 1,410 L of water. The actual cfu of the challenge dose was determined to be 1.8 x 104 cfu/mL by spreading dilutions of the bacterial culture on blood agar.
Approximately 1,200 experimental fish were randomly assigned to 1 of the 3 test tanks, with each tank having a volume of 1,800 L. The average fish BW of 96 g was estimated from a random sample, resulting in an approximate biomass of 60 kg/m3. Filtered, UV-treated, normal salinity seawater (34 parts per thousand) of 10°C was provided. Before the bath challenge, the water flow was stopped, and the volume was reduced to 1,410 L in each test tank. Thereafter, the suspension (details above) containing the V. anguillarum was added into the tank, and the fish were exposed to the bacteria for 60 min with oxygenation (requirement of 80 to 100% saturation). The water flow and volume in the tank were reattained postchallenge.
Fish were fed according to appetite (Dana Feed AS, Horsens, Denmark) through the whole experiment, and mortality was recorded daily. Samples from head kidney of dead or moribund fish were collected regularly during the challenge period for reisolation of V. anguillarum. The growth of the V. anguillarum bacteria was confirmed (MonoVa, Bionor AS, Skien, Norway) from all samples from fish that were collected from each test tank. The day of challenge was defined as d 1, and the experiment was terminated either 30 or 31 d postchallenge, depending on the test tank. Moribund and surviving fish were euthanized with an overdose of benzocaine (0.09 g/L, Benzoak vet., Europharma AS, Leknes, Norway).
Definition of Disease Resistance.
In all analyses, resistance to vibriosis was defined as the number of days from challenge until death, or until the end of the test (censored observations). Longer survival times were assumed to be associated with greater resistance to vibriosis.
Statistical Analysis
Preliminary Analysis.
The Kaplan-Meier estimate of the survival function (Kaplan and Meier, 1958) was plotted until the full termination of the experiment at d 31. Estimated curves of survivor function, 
(t), for categorical predictors (test tank, population, and family) were investigated, and the null hypothesis of whether 2 or more samples could have arisen from identical survivor functions was tested with nonparametric, univariate log-rank and Wilcoxon tests. Furthermore, visual examination of plots of –ln[(t)] vs. time (t) and ln[–ln[(t)]] vs. ln(t) was an indicator of whether the data might have arisen from an exponential or Weibull distribution, respectively. In addition, investigation of the latter gives an indication of the proportionality of hazards (parallel curves) across strata (Ducrocq et al., 2000). Preliminary analyses were conducted with the LIFETEST procedure (SAS Inst. Inc., Cary, NC).
Variance Component and Breeding Value Estimation.
To estimate variance components and the effects of independent variables on the survival time after challenge, 3 alternative approaches were used.
First, a univariate CRM, sire-dam frailty model was fitted (Cox, 1972):
 |
where 
i(t | sj, dk) is the hazard function of a fish i, conditional on random effects sj and dk; 0,m(t) is the unspecified baseline hazard function stratified by test tank (m = 1, 2, 3); b is the regression coefficient for tagging weight, (tagwti); sj is the random breeding value of the sire j; and dk is the random effect of the dam k. The effects of s and d were assumed to be normally distributed with zero means and variances of var(s) = I
s2 and var(d) = Id2, where I is an identity matrix.
Parametric models widely used in survival analysis, such as exponential and Weibull models, assume that the failure time distribution is known (Kalbfleisch and Prentice, 2002). In the current study, evaluation of plots of –ln[(t)] against t and ln[–ln[(t)]] against ln(t) revealed that neither exponential nor Weibull distributions were appropriate for the data. The CRM is nonparametric in the sense that it involves an unspecified function in the form of an arbitrary baseline hazard function (Kalbfleisch and Prentice, 2002).
As a result of the model optimization, a test tank-stratified population-wise model was chosen for estimation of variance components. The model was fitted using The Survival Kit v3.01 (Ducrocq and Sölkner, 1998b), which allows only 1 random parameter to be estimated at a time. Therefore, sire and dam variances were estimated sequentially, in a manner in which the estimate of the sire(dam) variance from the previous analyses was set as a known value when the dam(sire) variance was estimated. This chain of consecutive Cox regression runs was repeated until the variance component estimates remained unaltered. The variance components were estimated as the modes of the marginal posterior distributions (Ducrocq and Sölkner, 1998a). The Bayesian approach used by this software is described in detail in Ducrocq and Casella (1996). The estimates of breeding values were obtained as a by-product at convergence (Ducrocq and Sölkner, 1998a). The heritability for the number of days survived after challenge on the logarithmic-time scale was calculated according to Ducrocq and Casella (1996) as 
where 
2/6 is the variance for the residuals on an extreme value distribution, and the estimate for the proportion of variance due to the full-sib family was 
Variance components and breeding values for vibriosis resistance were also estimated assuming a univariate LMO:
 |
in which censored observations were treated as real observations; i.e., censoring day was assigned as the day of death. The test tank was included as a fixed effect in the model. Other model factors were as described for CRM. The effects of s, d, and e were assumed random with zero means and var(s) = Is2, var(d) = Id2, and var(e) = Ie2. The covariances between s, d and e were assumed to be zero. Variance components were estimated with the REML-method using the average information algorithm (Johnson and Thompson, 1995) by the DMU package (Madsen and Jensen, 2002). The heritability was calculated as 
and the proportion of variance due to the full-sib family as 
The SE of the estimates were derived from the average information matrix.
In the last approach, disease resistance was analyzed using a LMC in Bayesian analysis. The statistical model was identical to the model of LMO. Data were augmented from a truncated normal distribution for fish still alive at the end of the experiment, corresponding to censoring of those observations. The blocked Gibbs sampler (García-Cortés and Sorensen, 1996) was run as a single chain of 1,100,000 cycles, where the first 100,000 cycles were treated as conservative burn-in. After the burn-in period, every 10th sample was stored; i.e., a total of 100,000 samples were kept to describe the posterior distributions. The blocked Gibbs sampler was developed using C++. The estimates of h2 and f2 were calculated as for LMO.
Comparison of the Statistical Approaches.
Predicted breeding values for sire from the different statistical approaches were evaluated using Spearman rank correlation coefficients. For LMC, the means of the marginal posterior distributions for the sire breeding values were used. In addition, the predicted breeding values from CRM were reversed so that larger breeding values indicated a greater resistance and a similar ranking was expressed by a positive correlation.
|
RESULTS AND DISCUSSION |
|---|
|
TopAbstractINTRODUCTIONMATERIALS AND METHODSRESULTS AND DISCUSSIONIMPLICATIONSLITERATURE CITED |
|---|
. Only a very low mortality was observed during the first 7 d postchallenge (1.1%) followed by a drastic increase in mortality at d 8. During the following 5-d period (between 8 and 12 d) 54.4% of the fish died. This corresponded to an approximate daily mortality of 11.0%. The following 7-d period resulted in an additional 25.4% mortality. Thereby, the cumulative mortality at d 19 postchallenge was 80.8%. Period from d 20 until the termination of the experiment at d 31 was responsible only for 5.6% of the total mortality.
|
). Similar pattern of mortality was observed in the first challenge test conducted with Atlantic cod within the breeding program (Kettunen and Fjalestad, 2006). None of the independent variables could explain the bimodality of the distribution. The second peak might be explained by mortality caused by more virulent bacteria or secondary infections. Nevertheless, the time lapse between the 2 peaks (Figure 1) suggests that fish withstanding the original challenge are reinfected by bacteria released from sick and especially moribund and dead fish. The authors speculate that the second peak most likely describes mortality after such a rechallenge. Great variation was observed in survival process of different families. At the end of the test, the lowest and greatest family-wise survival was 1.4% (n = 74) and 37.9% (n = 95), respectively. The homogeneity of the survival curves across families was tested, and log-rank and Wilcoxon test statistics indicated differences (P < 0.001) between families in the pattern of survival after challenge.
The CCS population showed superior vibriosis resistance over NEAC and CCN in this experiment. The survival percentage for CCS was 17.0, whereas the percentages for NEAC and CCN were 10.3 and 11.0%, respectively. The log-rank and Wilcoxon test statistics indicated that the survivor curves were different (P < 0.001) among CCS, NEAC, and CCN. Inspection of the graphs of the survivor curves revealed that this was mainly a consequence of the dissimilarity of CCS (POP1) from CCN and NEAC. Additional Kaplan-Meier analysis showed that CCN and NEAC were likely to arise from the same function (P = 0.38), and thus these populations were combined for the genetic analyses (POP2). High water temperatures increase the risk for vibriosis (Hellberg, 2005). Consequently, vibriosis outbreaks are more abundant in the southern parts of Norway, where POP1 has its origin. Therefore, it is possible that CCS fish are frequently exposed to V. anguillarum and as a result of natural selection have developed an increased resistance to vibriosis. In contrast, vibriosis is less common in the northernmost coastal waters of Norway. Especially, NEAC mainly live in Barents Sea and schools to the coastal waters only to spawn and are less likely to experience a vibriosis pressure similar to that of CCS.
The survival percentages at the end of the experiment for the 3 test tanks were 17.3, 15.7, and 18.0 for POP1 and 9.2, 9.8, and 12.1 for POP2. Survivor curves across test tanks were different for POP2 according to the log-rank and Wilcoxon test-statistics (P = 0.034 and P = 0.002), whereas only Wilcoxon test showed statistical significance (P = 0.025) for POP1. Nevertheless, plots of ln[–ln[(t)]] against ln(t) for test tank-stratified Kaplan-Meier analyses resulted in nonparallel lines suggesting that the proportional hazards assumption is violated. Consequently, conclusions based on log-rank test might be inappropriate, and fitting separate baseline functions for test tanks in CRM is necessary.
Fish were tagged and weighed approximately 3 mo before the challenge. It was likely that variation in the family means of BW recorded at tagging (minimum 15.3 g, maximum 32.6 g) would exist at the time of the experiment. The preliminary analysis indicated that tagging weight had an effect on the number of days at test and was therefore included in the statistical model.
Variance Component Estimation
Cox Regression Model.
A larger estimate of sire variance was obtained for POP1 than for POP2, whereas the contrary was true for the variance of the full-sib family (Table 2). The SD of the marginal posterior distributions express the variation connected to the estimates (Table 2). The resulting heritability estimates for days survived postchallenge on the logarithmic scale were 0.10 and 0.08 for POP1 and POP2, respectively (Table 2). Only a negligible proportion of the total variation was explained by the full-sib family (f2 from 0.004 to 0.008). Visual diagnostics of baseline hazard curves, as well as weight-dependent predicted hazard curves, were conducted (Figure 2). Baseline hazard curves across test tanks were not proportional. In both populations, the greatest hazard at the early stage of the test was associated to same test tank, whereas the ranking of the test tank hazards across the populations were different at the end of the test. Moreover, tagging weight had no effect on the hazard in POP1 (data not shown), but it had a notable effect on the probability of survival during the challenge in POP2 (Figure 2). These observations, together with unequal variance components, suggest that the choice of the model for Cox regression analysis was reasonable.
|
|
, 2005) reported genetic parameters from Weibull frailty model for challenge test survival of viral hemorrhagic septicaemia (VHS), enteric redmouth disease (ERM), and rainbow trout fry syndrome (RTFS) in rainbow trout (Onchorhynchus mykiss). Their heritability estimates for disease resistance on the logarithmic scale varied from 0.07 to 0.21 depending on the disease. Gitterle et al. (2006) estimated moderate heritability for resistance to White Spot Syndrome Virus (WSSV) in Pacific white shrimp (Penaeus vannamei) with Weibull frailty model (0.21), whereas only very low estimate (0.01) was obtained from CRM. Thus, our estimates are within the range of values with those in previous literature.
Linear Models.
With LMO, larger estimates for sire and residual variances were obtained for POP1 than for POP2 (Table 2). Again, the contrary pattern was true for the variance of full-fib family. Estimates of sire and full-sib family variance were connected with large asymptotic SD (Table 2). As with CRM, the estimates of heritability for the 2 populations were within the range of each other given their respective SE: 0.16 ± 0.04 for POP1 and 0.13 ± 0.04 for POP2. The estimates of f2 were only 0.006 and 0.018 for POP1 and POP2, respectively. The differences between the effects of the test tanks on the number of days survived postchallenge were slightly larger in POP2 than in POP1. In addition, the ranking of the test tank solutions were unequal in the 2 populations (Table 2). No appropriate comparisons of the test tank effect on the disease resistance can be made between CRM and LMO. This is due to the fact that in LMO, the effect of each test tank is assumed to stay constant over the test period. In contrary, the test tank-stratified CRM allows the relativity of the baseline hazards to change over time. In agreement with the results from CRM, in POP2 1 unit change in tagging weight had a 7-fold effect on the number of days survived postchallenge compared with POP1 (Table 2).
The variance components and genetic parameters as well as the test tank and tagging weight effect estimates from the LMC are expressed as modes of the marginal posterior distributions (Figure 3). In general, LMC estimated larger sire and residual variances and considerably smaller full-sib family variance than LMO (Table 2). As with the 2 previous approaches, the estimates of sire variance were larger, whereas full-sib family estimates were smaller for POP1 than for POP2. In addition, a much larger estimate of residual variance was obtained for POP1 than for POP2. Hence, contradictory to the results from CRM and LMO, slightly greater survival heritability was estimated for POP2 (0.17) than for POP1 (0.14). In addition, the estimate of f2 was low and nearly equal in POP1 and POP2 with LMC (Table 2). As for CRM and LMO, the estimates were connected with considerable amount of uncertainty (Table 2). The predictors for the test tank and tagging weight effect from LMC were in agreement with those from LMO (Table 2).
|
; Gjedrem and Gjøen, 1995; Gjøen et al., 1997). Specifically after challenge with V. anguillarum, Gjøen et al. (1997) reported relatively high heritability estimates in Atlantic salmon (Salmo salar): 0.38 and 0.69 on the observed and on the underlying continuous liability scale, respectively. In more recent studies, Gitterle et al. (2006) estimated low heritabilities with linear (from 0.03 to 0.07) and threshold model (0.04) for resistance to WSSV in Pacific white shrimp. Similarly, Henryon et al. (2005) used a threshold model and obtained high heritability estimates for VHS, ERM, and RTFS (between 0.42 and 0.57) in rainbow trout. Nevertheless, these estimates were combined with a substantial uncertainty (Henryon et al. 2005). It has to be noted that in all these studies resistance was assessed as a binary trait (i.e., survived or died), whereas in the current study, resistance was defined as the number of days survived post-challenge.
Correlations Between Breeding Values Estimated with Different Approaches
The ranking of the sires according to the predicted breeding values with the different statistical approaches was very similar, with the correlation ranging between 0.97 and 0.99. Spearman rank correlation coefficients were calculated between POP2 sire breeding values and were close to unity. For POP1, the rank correlation coefficient between LMO and CRM breeding values was 0.99, between LMO and LMC was 0.98, and between CRM and LMC was 0.97. Henryon et al. (2005) reported rank correlations (0.84 to 0.89) between the predicted breeding values for sires and dams when disease resistance of ERM, RTFS, and VHS in rainbow trout was assessed as longitudinal and binary trait. Similarly, rank correlations from 0.82 to 1.0 (as absolute values) between estimated breeding values for full-sibs were obtained with 5 different statistical models in the analyses of WSSV in Pacific white shrimp (Gitterle et al., 2006).
It was also of interest to evaluate how well the family phenotypic survival percentages correlate with sire predicted breeding values. Rank correlation coefficients varied between 0.75 and 0.80 for POP1, and from 0.86 to 0.87 for POP2, indicating reranking of sires, and the necessity for vibriosis challenge test survival genetic evaluation for Atlantic cod.
In this study, univariate genetic analyses of vibriosis resistance with 3 statistical models revealed additive genetic variation sufficient for selective breeding. Disease resistance was assessed as a longitudinal trait, and only observations of most resistant fish were censored. As a result, LMO ranking of the sires was similar to those from CRM and LMC. In the case of more complicated censoring, LMO will not necessarily perform equally well. In a linear model without censoring it is not possible to discriminate mortality because of the challenge from mortality caused by other reasons, e.g., accidents in production or other environmental factors. We can discard the data from fish that have died for other reasons, but this would undermine the use of the information before the incident. The alternative is to use all data and assign the same number of days for the fish dying the same day, independent of the cause of death. Nevertheless, biased estimates of breeding values would be obtained. In contrast, a model accounting for censoring would distinguish different causes of mortality, utilize all available information, and result in more accurate ranking of individuals.
In this study, partitioning of the total variance into variance components is somewhat complicated. Due to confounding between some sires and dams, it is likely that a part of the dams’ additive genetic variation appears in the sire component. Similarly, the effects of the dam and the full-sib family coincide in almost all cases. As a result the heritabilities reported here are likely to be overestimates. This may also be the case for the estimates of f2. The mating design and statistical models used mean that sire variance is likely to express from 0.25 to 0.5 of the additive genetic variance. Nevertheless, the calculation of heritability was performed identically for all models, and the multiplication factor does not affect the results of the model comparison.
Breeding program of Atlantic cod will start out with selecting for 2 traits only: growth (as BW after 2 summers at sea) and vibriosis resistance. As the breeding program proceeds, new traits (e.g., avoidance of early maturity, meat quality traits) will be aggregated in the overall breeding goal. When censoring structure of the challenge test data is simple, as it was in this study, a multivariate linear model seems adequate for estimation of genetic correlations. However, the robustness of LMO does not necessarily apply to challenge test data in general. For example, problems during the challenge test may cause unequal challenge for different contemporary groups. In such circumstances it would be important to augment censored observations, as was done in LMC. Alternatively, a 2-step approach (Ducrocq et al., 2001; Tarres et al., 2006) or a recently developed Bayesian approach (Damgaard and Korsgaard, 2006) could be used for the simultaneous analysis of linear and survival traits. Nevertheless, in the starting phase of the breeding program, a multivariate linear model for estimation of breeding values for growth and disease resistance is likely to be adequate.
|
IMPLICATIONS |
|---|
|
TopAbstractINTRODUCTIONMATERIALS AND METHODSRESULTS AND DISCUSSIONIMPLICATIONSLITERATURE CITED |
|---|
|
Footnotes |
|---|
3 MTT, Biotechnology and Food Research, Biometrical Genetics, FIN-31600, Jokioinen, Finland.
2 Corresponding author: anne.kettunen{at}fiskeriforskning.no
Received for publication February 28, 2006. Accepted for publication September 8, 2006.
|
LITERATURE CITED |
|---|
|
TopAbstractINTRODUCTIONMATERIALS AND METHODSRESULTS AND DISCUSSIONIMPLICATIONSLITERATURE CITED |
|---|
This article has been cited by other articles:
|
|
 |
|
  C. Y. Feng, S. C. Johnson, T. S. Hori, M. Rise, J. R. Hall, A. K. Gamperl, S. Hubert, J. Kimball, S. Bowman, and M. L. Rise Identification and analysis of differentially expressed genes in immune tissues of Atlantic cod stimulated with formalin-killed, atypical Aeromonas salmonicida Physiol Genomics, May 13, 2009; 37(3): 149 - 163. [Abstract] [Full Text] [PDF]
|
|
|
|
 |
|
 T. Serenius, K. J. Stalder, and R. L. Fernando Genetic associations of sow longevity with age at first farrowing, number of piglets weaned, and wean to insemination interval in the Finnish Landrace swine population J Anim Sci, December 1, 2008; 86(12): 3324 - 3329. [Abstract] [Full Text] [PDF]
|
|
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
