An environmentally based growth model that uses finite difference calculus with maximum likelihood method: its application to the brackish water bivalve Corbicula japonica in Lake Abashiri, Japan.

Abstract--We present a growth analysis model that combines large amounts of environmental data with limited amounts of biological data and apply it to Corbicula japonica. The model uses the maximum-likelihood method with the Akaike information criterion, which provides an objective criterion for model selection. An adequate distribution for describing a single cohort is selected from available probability density functions, which are expressed by location and scale parameters. Daily relative increase tales of the location parameter are expressed by a multivariate logistic function with environmental factors for each day and categorical variables indicating animal ages as independent variables. Daily relative increase rates of the scale parameter are expressed by an equation describing the relationship with the daily relative increase rate of the location parameter. Corbicula japonica grows to a modal shell length of 0.7 mm during the first year in Lake Abashiri. Compared with the attainable maximum size of about 30 mm, the growth of juveniles is extremely slow because their growth is less susceptible to environmental factors until the second winter. The extremely slow growth in Lake Abashiri could be a geographical genetic variation within C. japonica.

**********

Extreme fluctuations, both short-term and seasonal, in food availability (e.g. phytoplankton density) make it difficult to determine relationships between the growth of filter-feeding bivalves and environmental factors (Bayne, 1993). However, it is becoming easier to acquire large amounts of environmental data through the use of data loggers, submersible fluorometers, of remote-sensing satellites, which enable environmental monitoring at daily or subdaily intervals. The development of these devices could solve difficulties in data collection. However, analytical methods that combine large amounts of environmental data with limited amounts of biological data (e.g. shell length) are not yet well developed. We present an environmentally based growth model that combines such unbalanced data sets. This model is useful in elucidating relationships between environmental factors and growth of filter feeders from field data.

Complex box models, such as ecophysiological models, can derive the relationships between environmental factors and the growth of filter-feeding bivalves (Campbell and Newell, 1998; Grant and Bacher, 1998; Scholten and Smaal, 1998). These models are useful for estimating impacts of cultivated species on an ecosystem or the carrying capacity of a species (or both) (Dame, 1993; Heral, 1993; Grant el al., 1993). They are suitable for animals that have been widely studied, such as Mytilus edulis, because they are derived by integrating a huge amount of ecophysiological knowledge acquired mainly from laboratory experiments. However, extrapolation of such knowledge to natural conditions is still controversial (Jorgensen, 1996; Bayne, 1998). Our model treats complicated eco-physiological processes as a black box; we constructed the model directly from fluctuations in environmental factors and growth rates. Our approach is reasonable for animals for which ecophysiological knowledge is limited, especially when the main purpose of investigation is to derive the relationships between environment and growth.

We applied the model to a single cohort of Corbicula japonica juveniles spawned in August 1997. We did not consider any bias caused by adjacent cohorts because C. japonica failed to spawn in 1995, 1996, and 1998 in Lake Abashiri owing to low water temperatures during the spawning season (Baba et al., 1999). Such investigations provide important basic information, such as the shape of the distribution of a single cohort, and the relationship between growth rate and expansion rate of size variation in a single cohort.

Corbicula spp. are harvested commercially in Japan. The annual catch ranged from 19,000 to 27,000 metric tons in 1996 to 2000 (Ministry of Agriculture, Forestry and Fisherles (1)), of which C. japonica was the main species. Corbicula japonica is distributed in brackish lakes and tidal flats of rivers from the south of Japan to the south of Sakhalin (Kafanov, 1991), is a dominant macrozoobenthos in these lakes, and has important roles in bioturbation and energy flow (Nakamura et al., 1988; Yamamuro and Koike, 1993). Juvenile C. japonica growth is fast in southern habitats. Their spats collected in Lake Shinji, which lies in the southern part of its range, grow to a mean shell length of around 6.7 mm in natural conditions by the first winter (Yamane et al. (2)). In northern habitats, growth is also believed to be fast; Utoh (1981) reported that mean shell length at the first annual mark was around 5.7 mm in Lake Abashiri. In Utoh's study differences between the shell lengths at the first annual marks and the shell lengths of individuals aged to be one year were also reported. The purposes of the present study are to elucidate juvenile growth and its relationship to environmental factors in Lake Abashiri.

Materials and methods

Overview of the model

Our model expresses relative growth rate for C. japonica by a sigmoid function with environmental factors and animal ages as independent variables. Modeling processes in general follow five steps: 1) Shell lengths of a single cohort are summarized by ah adequate probability density function, which is expressed by a location parameter and a scale parameter; 2) Daily relative increase rate of the location parameter (dRIRL) is approximated by a sigmoid function with environmental factors and animal ages as independent variables; 3) Daily relative increase rate of the scale parameter is approximated by a simple function with the dRIRL as an independent variable; 4) The model is optimized by a maximum likelihood method; and 5) The best model is selected by Akaike information criterion (AIC). The AIC is an information-theoretic criterion extended from Fisher's likelihood theory and is useful for simultaneous comparison of models (Akaike, 1973; Burnham and Anderson, 1998).

Study site and sampling method

To collect juveniles of C. japonica spawned in August 1997, sediments were sampled with a 0.05-m2 Smith-McIntyre grab once or twice a month from September 1997 to July 1999 at a depth of 3.5-4.0 m in Lake Abashiri (Fig. 1). The habitat of C. japonica is restricted to areas shallower than 6-m depth because the deeper area, the lower stratum of the lake, is covered by anoxic polyhaline water. We assumed that the selectivity of the sampling gear on C. japonica juveniles was negligible because the gear grabs the juveniles with the sediment. Because the magnitude of spawning in 1997 was relatively small (Baba et al., 1999), we selected a sampling site where we found abundant settled juveniles in our preliminary investigations. Samples could not be obtained during winter because of ice cover. Sediments were washed with tap water on 2- mm and 0.125-mm mesh sieves from September 1997 to October 1998, and on 4.75-mm and 0.125-mm mesh sieves from April to July 1999. To separate the juveniles from the retained sediments, we treated the sediments with zinc chloride solution as described by Sellmer (1956). Then we sorted the juveniles under a binocular microscope. Identification of the cohort spawned in 1997 was quite easy because C. japonica failed to spawn in 1995, 1996, and 1998 owing to low water temperatures during the spawning season (Baba et al., 1999). We considered all the individuals that passed through the larger-mesh sieves and that were retained on the smaller-mesh sieve as the 1997 cohort. Shell lengths were measured under a profile projector (V-12, Nikon Ltd., Chiyoda, Tokyo) at 50x magnification with a digital caliper (Digimatic caliper, Mitsutoyo Ltd., Kawasaki, Kanagawa), which has a 0.02-mm precision.

Environmental factors

Values for water temperature ([degrees]C), water fluorescenee (fluorescence equivalent to uranin density, [micro]g/L), salinity (psu, practical salinity unit), and turbidity (equivalent to kaolin density, ppm) were obtained for 0.1-m intervals from unpublished data at the Abashiri Local Office of the Hokkaido Development Bureau. (3) The variables were measured by a submersible fluorometer (Memory Chlorotec, ACL-1180-OK, Alec Electronics Ltd., Kobe, Hyogo) at four sites in Lake Abashiri at intervals of about one week (Fig. 1). The average values of each variable between the depths of 1 m and 6 m were used for later analyses. Values between the measured dates were interpolated linearly for subsequent analysis with the environmentally based growth model. The water fluorescence reflects the density of phytoplankton.

[FIGURE 1 OMITTED]

Model structure

Modeling the distribution of a single sample Normal distribution is usually used to describe a single cohort in fishes and aquatic invertebrates (e.g. Pauly, 1987; Founier and Sibert, 1990; Yamakawa and Matsumiya, 1997). However, an adequate function to describe a single cohort of each animal should be selected to avoid biases caused by any inadequacies of the function. Probability density functions of many distributions are applicable for that purpose, and the appropriate can be selected among easily calculable functions to ensure convergence of the model. Characteristics of many distributions are well described by Evans et al. (1993). We used two distributions: normal distribution and largest extreme value distribution. The normal distribution is symmetric. The largest extreme value distribution is asymmetric with a longer tail toward the larger side. These ate expressed by a location parameter and a scale parameter.

To use all the information inherent in data, parameters of the distribution functions are estimated from raw data (e.g. lengths), not from summarized data such as length frequency. This estimation method is described by Sakamoto et al. (1983). The most adequate distribution is selected by AIC. Log-likelihood functions of the distributions take the following forms:

Normal distribution

(1) [log.sub.e][L.sub.normal](a,b) = [n.summation over i=1][log.sub.e]{1/[square root of (2[pi][b.sup.2])] exp[-[([l.sub.1] - a).sup.2]]},

Largest extreme value distribution

(2) [log.sub.e][L.sub.largest](a,b) = [n.summation over i=1][log.sub.e]{(1/b)exp[-([l.sub.1] - a)/b] x exp{-exp[-([l.sub.a] - a)/b]}},

where n = number of data;

[l.sub.1] = length of ith individual;

a = location parameter; and

b = scale parameter.

The location parameter is a mean in the normal distribution. The location parameter is a mode in the largest extreme distributions. The scale parameter is a standard deviation in the normal distribution.

The AIC is calculated by

(3) AIC = -2 [log.sub.e] (maximum likelihood) + 2m,

where m = number of parameters to be estimated.

The model with the minimum AIC is the best model. A difference of more than 1 or 2 is regarded as significant in terms of AIC (Sakamoto et al., 1983).

Modeling the change in the location Values of the location and scale parameters usually increase with the growth of an animal. The relative increase rate in a certain time step is defined as

(4) [r.sub.i] = ([P.sub.i] - [P.sub.i=1])/[P.sub.i=1],

where [r.sub.i] = relative increase rate of a parameter in the ith time step; and

[P.sub.i] = parameter value after the ith time step.

Relationships between the parameter value and the relative increase rate of the parameter can be expressed by

(5) [P.sub.1] = [P.sub.0](1 + [r.sub.1])

[P.sub.2] = [P.sub.1](1 + [r.sub.2]) = [P.sub.0](1 + [r.sub.1])(1 + [r.sub.2])

[P.sub.3] = [P.sub.2](1 + [r.sub.3]) = [P.sub.0](1 + [r.sub.1])(1 + [r.sub.2])(1 + [r.sub.3])

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

where [P.sub.0] = parameter value at the first sampling;

[P.sub.i] = parameter value after the ith time step; and

[r.sub.i] = relative increase rate of the parameter in the ith time step.

We used one day us the time step in this study. In our environmentally based growth model, we assumed that the daily relative increase rate of location parameter (dRIRL) depends on the age of the animal and on environmental factors for each day. Sigmoid functions that take values between 0 and a certain maximum are empirically appropriate for expressing the relationships between the dRIRL and independent variables, especially for measures such as shell length that do not show negative growth. Therefore, using categorical variables indicating animal ages and environmental factors for each day as independent variables, we express the dRIRL by the multivariate logistic function

(6) [s.sub.i] = [s.sub.max]/{1 + exp[-([[n.sub.A].summation over j=1][[alpha].sub.j][A.sub.j] + [[n.sub.E].summation over j=1][[beta].sub.k][E.sub.ki])]},

where [s.sub.i] = dRIRL on the ith day from the first sampling;

[s.sub.max] = potential maximum dRIRL of the animal;

[[alpha].sub.j], [[beta].sub.k] = coefficients of each independent variable;

[A.sub.j] = categorical variable (a dummy variable indicating animal ages) that takes the value 1 or 0;

[E.sub.ki] = the kth environmental factor on the ith day from the first sampling;

[n.sub.A] = number of age categories; and

[n.sub.E] = number of environmental factors.

The categorical variable takes the value of 1 when the animal is the category, otherwise it takes 0. The multivariate logistic function with [s.sub.max] = 1 is used for logistic regressions (Sokal and Rohlf, 1995). A method of giving a value to the categorical variable is described by Zar (1999).

Modeling the change in scale The daily relative increase rate of scale parameter (dRIRS) and dRIRL must be correlated because the dRIRS is larger when the dRIRL is larger. Therefore, we estimated the dRIRS from an equation expressing the relationship to the dRIRL. We tested two functions,

(7) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

and

(8) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.],

where [t.sub.i] = dRIRS on the ith day from the first sampling;

[[gamma].sub.1], [[gamma].sub.2] = coefficients of the equations; and

[s.sub.i] = dRIRL on the ith day from the first sampling.

Model estimation

Likelihood function The location and scale parameters at the first sampling ([a.sub.0] and [b.sub.0]), the coefficients of Equation 6 ([s.sub.max], [[alpha].sub.a] and [[beta].sub.k]), and the coefficients of Equations 7 and 8 ([[gamma].sub.1] and [[gamma].sub.2]) are estimated as values that maximize total log-likelihood. The total log-likelihood is evaluated by the adequate probability density function selected in the first step. The log-likelihood functions take the following forms:

Normal distribution

(9) [log.sub.e][L.sub.normal]([a.sub.0], [b.sub.0], [s.sub.max], [[alpha].sub.j], [[beta].sub.k], [[gamma].sub.1], [[gamma].sub.2] = [N.summation over q=1][[n.sub.q].summation over i=1][log.sub.e]{1/[square root of (2[pi][b.sup.2.sub.q])]exp[-[l.sub.qi] - [a.sub.q])/2[b.sup.2.sub.q]]};

Largest extreme value distribution

(10) [log.sub.e][L.sub.largest]([a.sub.0], [b.sub.0], [s.sub.max], [[alpha].sub.j], [[beta].sub.k], [[gamma].sub.1], [[gamma].sub.2] = [N.summation over q=1][[n.sub.q].summation over i=1][log.sub.e]{(1/[b.sub.q])exp[-[l.sub.qi] - [a.sub.q])/[b.sub.q]] x exp{-exp[-([l.sub.qi] - [a.sub.q]]}},

where [a.sub.0], [b.sub.0] = values of the location and scale parameters, respectively, at the first sampling;

[s.sub.max], [[alpha].sub.j], [[beta].sub.k] = coefficients of Equation 6;

[[gamma].sub.1], [[gamma].sub.2] = coefficients of Equations 7 and 8;

N = number of samplings;

[n.sub.q] = number of data at the qth sampling;

[a.sub.] = location parameter at the qth sampling estimated by Equation 5 ([r.sub.i]=[s.sub.i]);

[b.sub.q] = scale parameter at the qth sampling estimated by Equation 5 ([r.sub.i]=[t.sub.i]); and

[l.sub.qi] = length of the ith individual at the qth sampling.

AIC is used to select significant environmental factors, the age categorization, and the equation to express the relationship between dRIRL and dRIRS, Le. Equation 7 or 8.

Confidence intervals To evaluate uncertainties of coefficient values and model selection, we estimated the 95% confidence intervals of all coefficients--i.e, [a.sub.0], [b.sub.0], [s.sub.max], [[alpha].sub.j], [[beta].sub.k], [[gamma].sub.1], and [[gamma].sub.2]--based on profile likelihood. For example, the 95% confidence interval of [a.sub.0]--[a.sub.0.95]--was estimated as an interval that suffices in the following equation:

(11) 2{max [log.sub.e] L([a.sub.0], [b.sub.0], [s.sub.max], [[alpha].sub.j], [[beta].sub.k], [[gamma].sub.1], [[gamma].sub.2])

- max [log.sub.e] L([a.sub.0], [b.sub.0], [s.sub.max], [[alpha].sub.j], [[beta].sub.k], [[gamma].sub.1], [[gamma].sub.2]

|[a.sub.0] = [[alpha].sub.0,96])}[less than or equal to][[chi square].sub.1](0.05),

where [x.sub.1.sub.2](0.05) = value of a chi-squared distribution at an upper probability of 0.05 with one degree of freedom, i.e. 3.84.

The characteristics of the interval are explained by Burnham and Anderson (1998).

We used Microsoft Excel (Microsoft Corp., Redmond, WA) as the analysis platform, and Solver (Microsoft Corp., Redmond, WA) as the nonlinear optimization tool.

Model selection

We used three procedures for model selection to achieve the best model. First, we constructed an a priori set of base models based on biological variables; then we selected the best base model. Fixation of the base model drastically decreases possible candidate models to be tested. To test all possible combinations of independent variables and model forms is quite impractical, Second, we excluded insignificant factors from the best base model. Third, we checked the significance of environmental factors that were not included in the base models. If one was significant, we included it in the best base model. All of these procedures were performed by AIC. The construction of the a priori set of candidate models is partially subjective, but it is an important part of the model construction (Burnham and Anderson, 1998).

Seasonal growth in bivalves is influenced by water temperature and food supply (Bayne and Newell, 1983). The growth rate of Corbicula fluminea changes with age (McMahon, 1983). Therefore, we constructed base models combining water temperature, water fluorescence, and categorical variables indicating age for the independent variables of Equation 6. We tested two types of categorization of age. The first segregates ages based on real age, i.e. two categories: 0+ or 1+. The second segregates ages in relation to winter, i.e. three categories: before the first winter, from the first to the second winter, and after the second winter. For the real-age categorization, age was segregated based on 1 September, because the spawning season was in August 1997. For the winter-base age categorization, we segregated ages based on 1 January. No biases should have occurred because of the segregation date of the winter base categorization and because the growth of C. japonica is negligible during winter. Four base models were constructed combining the two types of age categorization and two types of equations expressing the relationship between the dRIRL and the dRIRS, i.e. Equations 7 or 8. We selected the best base model by AIC.

To check the significance of each environmental factor and age categorization, we removed the independent variables one by one from the best base model and re-optimized the model. When the model was significantly improved by the removal in terms of AIC, the effect of the variable was insignificant on the model; therefore we excluded it.

To check the significance of salinity and turbidity, which were not included in the base models, we included them one at a time into the best base model and re-optimized the model. When the model was improved by the inclusion, the effect of the variable was significant on the model; therefore we included it.

Results

Modeling the distribution of a single sample

The largest extreme value distribution was the best in terms of AIC except for data sampled on 13 May 1998 (results are not shown). The exception is due probably to the small sample size (n=38) on that date. The largest extreme value distribution was therefore used to evaluate likelihood in later analyses: we selected Equation 10 from Equations 9 and 10. The result of fitting the two distributions to the shell lengths sampled on 22 April 1999 is shown in Figure 2 as a representative example. The largest extreme value distribution is apparently better than the normal distribution for describing the single cohort of C. japonica spawned in 1997.

[FIGURE 2 OMITTED]

Model selection and application

Model 4 was the best in terms of AIC among four base models (Table 1, models 1-4); ages were categorized in relation to winter; and the relationship between dRIRL and dRIRS was expressed by Equation 8.

Four models were made by removing each independent variable from model 4 (Table l, models 4.1 to 4.4). The effect of one age categorization--segregation of ages between the first and second winters--was insignificant on the model, because the model was significantly improved by its removal in terms of AIC. The effects of the other independent variables were significant on the model, because the model was significantly worse by their removal in terms of AIC. The effects of salinity and turbidity were insignificant on the model, because adding each variable made the model significantly worse in terms of AIC (Table 1, models 4.5 and 4.6). Consequently, model 4.1 was the best model to describe the relationships among environmental factors, ages, and growth of C. japonica juveniles spawned in 1997.

The coefficient value for age categorization of before the second winter (18.3) is much smaller than that of after the second winter (-10.0) (Table 1). This difference suggests that the growth response of C. japonica juveniles is much less susceptible to environmental factors before the second winter than after.

Peaks of the dRIRL corresponded with peaks of water fluorescence, when the water temperature was warmer than about 10[degrees]C, especially before the second winter (Fig. 3, B and C). Therefore, food supply is the most influential factor when the water temperature is above about 10[degrees]C. The slow growth or no growth during winter is due to the low water temperatures. The dRIRL reached a plateau after 30 May 1999. This was due to two factors: water fluorescence was relatively intense after 30 May 1999 (Fig. 3B); and the growth response of C. japonica to the environmental factors was more susceptible after the second winter than before.

[FIGURE 3 OMITTED]

The confidence limits of all the coefficients seem to be reasonably estimated by the profile likelihood method (Table 2). These results also guarantee the convergence of the model because the model was frequently optimized to seek each confidence limit with different starting values. We repeated the optimization at least 20 times to seek each confidence limit. On other models, we also confirmed the convergences as well.

The largest extreme value distributions estimated by model 4.1 fitted the shell lengths of C. japonica juveniles very well (Fig. 4).

[FIGURE 4 OMITTED]

Discussion

Model formulation and application

Largest extreme value distribution is apparently better than normal distribution to describe the single cohort of C. japonica that spawned in 1997. This distribution has a mode and a longer tail toward the larger side. If the shell length distribution becomes asymmetric during growth, skewness of the distribution would increase according to growth. However, there is no correlation between the skewness and the means of the shell lengths. Therefore, we thought that the shell length distribution of the cohort was already asymmetric just after settlement. Such a distribution might be influenced by fluctuations in larval settlement during the spawning season; and larval settlement would be influenced by fluctuations in larval supply from the water column. During the spawning season of 1997, the average planktonic larval density gradually increased from 26 ind/[m.sup.3] on 25 July to a maximum of 603 ind/[m.sup.3] on 13 August. Then it sharply decreased to 3 ind/[m.sup.3] on 19 August (Baba et al., 1999). Such a pattern of larval-density fluctuation might have caused the asymmetric distribution of shell lengths of the settled juveniles. Another possible factor that influenced the shapes of the shell length distributions and the relationship between dRIRL and dRIRS is size-dependent mortality, e.g. predations and fisheries. Size-dependent mortality has been reported in several marine bivalves (e.g. Nakaoka, 1996). Potential predators of C. japonica are fishes, such as Japanese dace (Tribolodon hakonensis) (also known as big-scaled Pacific redfin, FAO), Pacific redfin (Tribolodon brandtii), common carp (Cyprinus carpio), and the So-iny mullet (Liza haematocheila) (Kawasaki (4)), In our study, the size-dependent mortality was negligible because the range of the shell lengths observed in this study was very narrow.

The shape of the distribution to describe a single cohort should be determined from the data. In contrast, single cohorts are usually separated from multicohort data by assuming a normal distribution of lengths in a single cohort (e.g. Fournier and Sibert, 1990). Therefore, it is possible that multicohort analysis done without selection of an adequate distribution to describe a single cohort causes substantial bias in estimations of various stock features of animal populations, such as age composition, growth, mortality, and recruitment. In our preliminary analyses, we also tested smallest extreme value distribution, inverse Gaussian distribution, and lognormal distribution. The inverse Gaussian distribution was the best for two samples; the lognormal distribution, was the best for two samples; the largest extreme value distribution was the best for ten samples. Therefore, it is reasonable to select the largest extreme value distribution. We selected a single distribution for our analyses, otherwise a discontinuous point would have appeared in the growth curve.

Relatively large confidence intervals were obtained in the coefficients of the linear component of Equation 6, i.e. [[alpha].sub.j], and [[beta].sub.k] (Table 2). The relatively large confidence intervals may indicate that the number of estimated coefficients is somewhat larger than the number of samplings. Therefore, to estimate these coefficients more precisely, we may need to investigate more cohorts spawned in other years in future investigations.

Growth of C. japonica

We identified extremely slow growth in C. japonica juveniles, which grew to a modal shell length of 0.7 mm during the first year in Lake Abashiri, which lies at 43.7[degrees]N. Spats of C. japonica collected from 1992 to 1997 in Lake Shinji, which lies at 35.5[degrees]N, grew to a mean shell length of 6.7 mm in natural conditions by the first winter (Yamane et al. (2)). Using environmental factors measured in Lake Shinji from 1990 to 1998 at monthly intervals (Seike (5)), we simulated the growth of C japonica with model 4.1. Corbicula japonica grew to a mean shell length of 1.4 mm (standard error, 0.37) by the first winter in the simulations. Therefore, the large difference in juvenile growth between the two habitats cannot be explained by environmental differences because the results of the simulation were apparently an underestimate. We think that the extremely slow growth of the juveniles (prolonged phase of meiobenthic development) in Lake Abashiri is probably a geographical variation, which is genetically determined, within C. japonica. However, there remains a possibility that the juvenile growth differences depend on other environmental factors not measured in this study. Therefore, the geographical variation should be validated by reciprocal transplantations or laboratory experiments (or both) in future investigations, Prolonged phases of meiobenthic development have been reported in some marine bivalves (Nakaoka, 1992; Halwey and Gage, 1995). However, a prolonged phase of meiobenthic development as a geographical variation is rarely reported.

In many species of bivalve, populations from higher latitudes have a slower initial growth rate; but longevity and ultimate size in these populations are frequently greater than at lower latitudes (Newell, 1964; Seed, 1980). The extremely slow growth of C. japonica juveniles in Lake Abashiri may be an extreme example of this phenomenon. In Lake Abashiri, C. japonica failed to spawn in ten out of 21 years for which data were available because of low water temperatures during the summer spawning season (Baba et al., 1999). This means that a long life span is essential to sustain populations of C. japonica in northern habitats. We think that a long life span is the ultimate factor for the extremely slow growth rate of C. japonica juveniles in Lake Abashiri.

The growth response of C. japonica juveniles is much less susceptible to environmental factors before the second winter than after and is the proximate factor for an extremely slow growth rate. Nuculoma tennis, a detritus feeder, develops its palp proboscides, its feeding apparatus, during the prolonged phase of meiobenthic development (Harvey and Gage, 1995). The change of growth susceptibility to environmental factors in young ages may suggest that some functional morphological changes occur in C. japonica, also a filter feeder. In our preliminary analyses, we could not find a better model when we used different values of [s.sub.max] in Equation 6 between ages instead of categorical variables indicating ages. Therefore, we conclude that the difference in growth rates between ages is not due to a difference in potential maximum growth rate, at least in the range of the shell length observed in our study. When our model is applied to a wider range of the shell lengths of other species, it is best to examine the age dependence of [s.sub.max]

Table 1

Values of location and scale parameters at the first sampling,
coefficients, log-likelihood, and AIC of models constructed based on
the largest extreme value distribution. The best AIC among four base
models (models 1-4) is enclosed by a single line. The best AIC of all
models is enclosed by a double line. dRIRL = daily relative increase
rate of location parameter, dRIRS = daily relative increase rate of
scale parameter, Temp. = water temperature, WF = water fluorescence,
Sal. = salinity, Turb. = turbidity, C1 = before the first winter,
C2 = from the first to the second winter, C3 = after the second
winter.

              Parameters
                at 1st             Max.
               sampling            dRIRL
Model
no.     [a.sub.0]   [b.sub.0]   [s.sub.max]

 1        0.299       0.040        0.012
 2        0.297       0.040        0.011
 3        0.299       0.042        0.011
 4        0.299       0.042        0.011
4.1       0.299       0.042        0.011
4.2       0.297       0.038        0.005
4.3       0.295       0.037        0.008
4.4       0.299       0.041        0.013
4.5       0.299       0.042        0.011
4.6       0.299       0.042        0.011

                    Age categorization

Model     [A.sub.1]       [A.sub.2]       [A.sub.3]
no.     [alpha.sub.1]   [alpha.sub.2]   [alpha.sub.3]

               0+              1+
 1          -62.6           -23.7
 2          -56.1           -22.1
               C1              C2             C3
 3          -16.8           -16.7           -9.1
 4          -17.5           -17.6           -9.6
4.1               -18.3 (1)                -10.0
4.2         127.9                 -26.8 (1)
4.3         -47.3           -16.3           -8.8
4.4          -4.9            -8.9           -4.9
4.5               -16.7 (1)                 -9.1
4.6               -18.5 (1)                -10.2

                    Environmental factors

          Temp.          WF          Sal.        Turb.
Model    [[beta]-     [[beta]-     [[beta]-     [[beta]-
no.      .sub.1]      .sub.2]      .sub.3]      .sub.4]

 1         0.16         2.61
 2         0.20         2.44
 3         0.61         0.41
 4         0.65         0.42
4.1        0.68         0.44
4.2        0.34         4.15
4.3                     1.47
4.4        0.40
4.5        0.62         0.42        -0.25
4.6        0.68         0.44                     0.007

            Expressing relationship
                 betweeen dRIRS
                    and dRIRL

Model   [[gamma]-    [[gamma]-     Eq.
no.      .sub.1]      .sub.2]      no    Log-L     AIC

 1        0.0000       1.686       7     850.4   -1682.9
 2        0.0001       0.887       8     852.3   -1686.5
 3       -0.0076       2.902       7     950.4   -1880.9
 4        0.0034       0.760       8     952.2   -1884.4
4.1       0.0034       0.760       8     952.2   -1886.3
4.2       0.0000       0.895       8     735.0   -1451.9
4.3       0.0033       0.766       8     848.9   -1679.9
4.4       0.0020       0.806       8     909.6   -1801.1
4.5       0.0033       0.762       8     952.4   -1884.8
4.6       0.0034       0.760       8     952.2   -1884.4

(1) One common coefficient was used for the two categorical variables.

Table 2

95% confidence limits of location and scale parameters at the first
sampling and coefficients of the best model constructed based on the
largest extreme value distribution (models 4.1 in Table 1) estimated by
profile likelihood method. dRIRL = daily relative increase rate of
location parameter, dRIRS = daily relative increase rate of scale
parameter, Temp. = water temperature, WF = water fluorescence, Sal. =
salinity, Turb. = turbidity.

                   Parameters
                     at 1st              Max.
                    sampling            dRIRL

             [a.sub.0]   [b.sub.0]   [s.sub.max]

Lower 95 %     0.294       0.039        0.010
Upper 95 %     0.304       0.045        0.013

                     Age categorization

             [A.sub.1]   [A.sub.2]   [A.sub.3]
              [alpha-     [alpha-     [alpha-
              .sub.1]     .sub.2]     .sub.3]

Lower 95 %   -26.6 (1)                 -14.6
Upper 95 %   -11.5 (1)                  -6.4

                      Environmental factors

              Temp.        WF        Sal.      Turb.
             [[beta]-   [[beta]-   [[beta]-   [[beta]-
             .sub.1]    .sub.2]    .sub.3]    .sub.4]

Lower 95 %     0.41       0.27
Upper 95 %     1.00       0.64

                  Expressing
                 relationship
                betweeen dRIRS
                   and dRIRL

             [[gamma]-   [[gamma]-
              .sub.1]     .sub.2]

Lower 95 %    0.0027       0.734
Upper 95 %    0.0039       0.793

(1) One common coefficient for the two categorial variables.

Acknowledgments

We express our thanks to T. Kato, Vice-Head of the River Improvement Section in the Abashiri Local Office of the Hokkaido Development Bureau, for providing environmental data on Lake Abashiri. We also thank the reviewers of Fishery Bulletin for providing helpful suggestions on our manuscript.

(1) Ministry of Agriculture, Forestry and Fisheries. 1996 2002. Statistics on fisheries and water culture production. Association of Agriculture and Forestry, 1-2-1 Kasumigaseki, Chiyoda, Tokyo 100-0013, Japan.

(2) Yamane, K., M. Nakamura, T. Kiyokawa, H. Fukui, and E. Shigemoto. 1999. Experiment on the artificial spat collection. Bull. Shimane Pref. Fish. Exp. Stn., p. 232 234. Unpubl. rep. Shimane Prefectural Fisheries Experimental Station, 25-1 Setogashima, Hamada, Shimane 697-0051, Japan. [In Japanese.]

(3) Abashiri Local Office of the Hokkaido Development Bureau, 2-6-1 Shinmachi, Abashiri, Hokkaido 093-0046, Japan.

(4) Kawasaki, K. 1997. Lagoon structure and fish production fu Ogawara-ko Lagoon. In Final reports on fisheries in Ogawara-ko Lagoon (Tohoka Construction Corporation ed.), p. 4-33. Unpubl. rep. Construction Office for Takasegawa General Development of Tohoku Regional Construction Bureau, 3 Ishido, Hachinohe, Aomori 039-1165, Japan.

(5) Seike, Y. 1990-98. Gobiusu: monthly report of water quality in Lake Shinji and Lake Nakaumi. Unpubl. rep. Faculty of Science and Engineering, Shimane University, 1060 Nishikawatsu, Matsue, Shimane 690-0823, Japan.

Literature cited

Akaike, H.

1973. Information theory and an extension of the maximum likelihood principle. In Second international symposium on information theory (B. N. Petrov and F. Csaki, eds.), p. 267-281. Akademiai Kiado, Budapest.

Baba, K., T. Kawajiri, and Y. Kuwahara.

1999. Effects of water temperature and salinity on spawning of the brackish water bivalve Corbicula japonica in Lake Abashiri, Hokkaido, Japan. Mar. Ecol. Prog. Ser. 180: 213-221.

Bayne, B. L.

1993. Feeding physiology of bivalves: time dependence and compensation for changes in food availability. In Bivalve filter feeders and marine ecosystem processes (R. F. Dame ed.), p. 1-24. Springer Verlag, New York, NY.

1998. The physiology of suspension feeding by bivalve molluscs: an introduction to the Plymouth "TROPHEE" workshop. J. Exp. Mar. Biol. Ecol. 219:1 19. Bayne, B. L., and R. C. Newell.

1983. Physiological energetics of marine molluscs. In The Mollusca, 4(1) (A. S. M. Saleuddin and K. M. Wilbur eds.), p. 407-515. Academic Press, New York, NY.

Burnham, K. P., and D. R. Anderson.

1998. Model selection and inference, 349 p. Springer, New York, NY.

Campbell, D. E., and C. R. Newell.

1998. MUSMOD, a production model for bottom culture of the blue mussel, Mytilus edulis L. J. Exp. Mar. Biol. Ecol. 219:171 203.

Dame, R. F.

1993. The role of bivalve filter feeder material fluxes in estuarine ecosystems. In Bivalve filter feeders in estuarine and coastal ecosystem processes: NATO ASI Series (R. F. Dame ed.), p. 245-269. Springer-Verlag, Berlin, Heidelberg.

Evans, M., N. Hastings, and B. Peacock.

1993. Statistical distribution, 170 p. John Wiley & Sons, New York, NY.

Fournier, D. A., and J. R. Sibert.

1990. MULTIFAN a likelihood-based method for estimating growth parameters and age composition from multiple length frequency data sets illustrated using data for southern bluefin tuna (Thunnus rnaccoyii). Can. J. Fish. Aquat. Sci. 47:301-317.

Grant, J., and C. Bacher.

1998. Comparative models of mussel bioenergetics and their validation at fie]d culture sites. J. Exp. Mar. Biol. Ecol. 219:21-44.

Grant, J., M. Dowd, K. Thompson, C. Emerson, and A. Hatcher.

1993. Perspectives on field studies and related biological models of bivalve growth and carrying capacity. In Bivalve filter feeders in estuarine and coastal ecosystem processes: NATO ASI Series (R. F. Dame ed.), p. 371-420. Springer-Verlag, Berlin, Heidelberg.

Harvey, R., and J. D. Gage.

1995. Reproduction and recruitment of Nuculoma Tenuis (Bivalvia: Nuculoida) from Loch Etive, Scotland. J. Moll. Stud. 61:409-419.

Heral, M.

1993. Why carrying capacity models are useful tools for management of bivalve molluscs culture. In Bivalve filter feeders in estuarine and coastal ecosystem processes: NATO ASI Series (R. F. Dame ed.), p. 455-477. Springel-Verlag, Berlin, Heidelberg.

Jorgensen, C. B.

1996. Bivalve filter feeding revisited. Mar. Ecol. Prog. Ser. 142:287-302.

Kaihaov, A. I.

1991. Bivalves on continental shelf and continental slope of Northern Pacific, p. 81-82. Science Acad. USSR, Vladivostok.

MacMahon, R. F.

1983. Ecology of an invasive pest bivalve, Corbicula. In The Mollusca, 6(12) (W. D. Russell-Hunter, ed.), p. 505-561. Academic Press, Orlando, FL.

Nakamura, M., M. Yamamuro, M. Ishikawa, and H. Nishimura.

1988. Role of the bivalve Corbicula japonica in the nitrogen cycle in a mesohaline lagoon. Mar. Biol. 99: 369-374.

Nakaoka, M.

1992. Age determination and growth analysis based on external shell rings of the protobranch bivalve Yoldia notabilis Yokoyama in Otsuchi Bay, Northeastern Japan. Benthos Res. 43:53-66.

1996. Size-dependent survivorship of the bivalve Yoldia notabilis (Yokoyama, 1920): the effect of crab predation. J. Shellfish Res. 15:355-362.

Newell, G. E.

1964. Physiological aspects of the ecology of intertidal molluscs. In Physiology of Mollusca (K. M. Wilbur and C. M. Yonge eds.), p. 59-81. Academic Press, London.

Pauly, D.

1987. A review of the ELEFAN system for analysis of length frequency data in fish and aquatic invertebrates. In Length-based methods in fisheries research (D. Pauly and G. R. Morgan eds.), p. 7-34. ICLARM conference proceedings 13, 688 International Center for Living Aquatic Resource Management, Manila, Philippines and Kuwait Institute for Scientific Research, Safet, Kuwait.

Sakamoto, Y., M. Ishiguro, and G. Kitagawa.

1983. Statistics based on information amount (Jouhouryou Toukei Gaku), 236 p. Kyoritu Shuppan, Tokyo. [In Japanese.]

Scholten, H., and A. C. Smaal.

1998. Responses of Mytilus edulis L. to varying food concentrations: testing EMMY, an ecophysiological model. J. Exp. Mar. Biol. Ecol. 219:217-39.

Seed, R.

1980. Shell growth and form in the Bivalvia. In Skeletal growth of aquatic organisms (D. G. Rhoads ed.), p. 23-67. Plenum Press, New York, NY.

Sellmer, G. P.

1956. A method for the separation of small bivalve molluscs from sediments. Ecology 37:206.

Sokal, R. R. and F. J. Rohlf.

1995. Biometry, 3rd ed., 887 p. W.H. Freeman and Company, New York, NY.

Utoh, H.

1981. Growth of the brackish water bivalve, Corbicula japonica Prime, in Lake Abashiri. Sci. Rep. Hokkaido Fish. Exp. Stn. 23:65-81. [In Japanese.]

Yamakawa, T,. and Y. Matsumiya.

1997. Simultaneous analysis of multiple length frequency data sets when the growth rates fluctuate between years. Fish. Sci. 63:708-714.

Yamamuro, M., and I. Koike.

1993. Nitrogen metabolism of the filter-feeding bivalve Corbicula japonica and its significance in primary production of a brackish lake in Japan. Limnol. Oceanogr. 38: 997-1007.

Zar, J. H.

1999. Biostatistical analysis, 663 p. Prentice Hall, Upper Saddle River, NJ.

Katsuhisa Baba

Hokkaido Hakodate Fisheries Experiment Station

1-2-66, Yunokawa, Hakodate Hokkaido 042-0932, Japan

E-mail address: [e-mail]

Toshifumi Kawajiri

Nishiabashiri Fisheries Cooperative Association

1-7-1, Oomagari, Abashiri Hokkaido 093-0045, Japan

Yasuhiro Kuwahara

Hokkaido Abashiri Fisheries Experiment Station

31, Masuura, Abashiri Hokkaido 099-3119, Japan.

Shigeru Nakao

Graduate School of Fisheries Sciences

Hokaido University

3-1-1, Minato, Hakodate Hokkaido 041-8611, Japan

Manuscript approved for publication 14 August 2003 by Scientific Editor.

Manuscript received 20 October 2003 at NMFS Scientific Publications Office.