Rainfall effects on Anomalocardia flexuosa densities on the northeastern Brazilian coast using distributed lag models

In this study, the effects of precipitation on Anomalocardia flexuosa densities were analyzed from the perspective of temporal delays between the variables. The collections occurred bi-monthly between April 2016 and February 2018 at Mangue Seco beach, Pernambuco, Brazil. Total densities and den- sities according to the category of size (small, medium and large), and precipitation were analyzed with autocorrelation and cross-correlation functions, with a retrospective analysis of up to 12 months. Distributed lag models were used among the variables. The maximum density was found in June 2017 (378 ind.m 2 ) for the medium-size category of individuals. Significant autocorrelations occurred for densities (total and average size) and precipitation. In cross-correlation functions, high precipitation for a given month was followed by high densities of total and average individuals for up to two months. Distributed lag models also showed significant values between densities (total and medium-sized) and precipitation, which explained more than 50% of the variability of these two groups. The effect of precipitation was responsible for the increase in the density of for up to two months, mainly by increasing individuals with medium shell lengths. (378 ind.m²) médios até meses. Modelos de defasagem nificativos entre as densidades (total e médio) explicando mais de 50% da variabilidade desses dois grupos. O efeito da precipitação foi responsável pelo aumento da densidade de A. flexuosa por até dois meses, principalmente pelo incremento de indivíduos com comprimentos de concha medianos. Palavras-chave: precipitação; exploração de bivalves; planície de maré.


INTRODUCTION
Fisheries have great global importance and approximately half of their production comes from small-scale fisheries, which employ approximately 90% of the workforce of this activity (FAO, 2018b). Proportionally, small-scale extraction is significant and even greater in the production of bivalve mollusks, and this resource is an important source of income for the subsistence of traditional communities residing in coastal areas distributed throughout the planet (Castilla and Defeo, 2001;Silva-Cavalcanti and Costa, 2011;Oliveira et al., 2013;López-Rocha et al., 2020). In 2018, catches of marine and estuarine mollusks were estimated at approximately 6 million tons, and bivalves make up a substantial part of the total catch (FAO, 2018a).
Among the bivalve mollusks, those of the family Veneridae are the most traditionally exploited by commercial fisheries (Costa et al., 2020;López-Rocha et al., 2020), and in this family, Anomalocardia flexuosa (Linnaeus, 1767) stands out. This species occurs in the Atlantic Ocean from the West Indies down as far as Uruguay (Rios, 2009), but its commercial exploitation is concentrated on the Brazilian coast (Silva-Cavalcanti and Costa, 2011), where it inhabits different environments, including intertidal flats (Narchi, 1974). These areas are greatly influenced by rivers and estuaries (Rodrigues et al., 2013;Silva-Cavalcanti et al., 2018).
Species living in intertidal habitats, including A. flexuosa, have temporal distributions correlated with the interactions of physical (e.g. temperature and salinity, tidal amplitude) and biological factors (e.g. predation, competition) (Gerwing et al., 2016). However, in coastal ecosystems, these variations in densities are complex and not easily explained by linear dynamics associated with only one of these factors (Dame et al., 2002).
Among the environmental variables that are able to modify both the physical and biological factors of bivalves, rainfall is the variable that has the greatest influence (Silva-Cavalcanti et al., 2018;Almeida et al., 2021). Coupled with the absence of estimates of recruitment rates that can predict future densities, it is suggested that environmental variables that affect recruitment can be used in predictive analyses (Almeida et al., 2021). These variables can control the success of the recruitment of bivalve species with a short life cycle, and can be recruited for fishing during the first year of life (Gaspar et al., 2004).
Rainfall has a strong influence on the densities of A. flexuosa, as it can affect the fluctuations of salinity and fishing, among other factors (Belém et al., 2013;Rodrigues et al., 2013;Oliveira et al., 2014;Silva-Cavalcanti et al., 2018). In addition, it can also influence the different size classes of the species, which is a factor that is directly associated with fishing production (Oliveira et al., 2014;Rodrigues et al., 2018). However, in these previous studies, the relationships between rainfall and bivalve density were reported in a secondary way, in atypical periods of rain or with the period considered as a fixed factor, without the use of models that incorporated temporal lags between the variables.
Distributed lag models (DLM) are applied effectively in various environmental studies that consider the effect of a covariable being delayed over time (Rushworth et al., 2013). Therefore, monitoring of precipitation is important when estimating and forecasting trends of fishery stocks for the strategic management of bivalve mollusks (Almeida et al., 2021). Understanding the changes in density patterns correlated with precipitation through estimates that consider time lag is fundamental for understanding the population dynamics of A. flexuosa. Environmental influences on biota are common in nature, including aquatic organisms, so time series analyses have become powerful monitoring tools that can explain the increase in bivalve densities at the present time, though such densities have been stimulated as a result of natural phenomena occurring in the past (Olden and Neff, 2001). Thus, this study analyzes the effects of precipitation on the densities of the bivalve A. flexuosa from the perspective of time delays between the variables in an intertidal flat on the northeastern coast of Brazil.

Study area and data collection
Mangue Seco beach is located in the northern region of the state of Pernambuco, on the northeastern coast of Brazil, and has a warm and humid climate, with annual rainfall exceeding 1,000 mm. It registers rainfall throughout the year, and the periods of the rainfall regime are well-defined; September to January comprises the dry season and February to August comprises the rainy season (Moura, 2009). This is a traditional extractive region covering an area of approximately 2.9 km 2 . It is a tidal flat and has great influence due to its contribution of fresh water and nutrients from estuaries located to the South, at the mouth of the Timbó River, and to the north, at the Santa Cruz channel (Lima et al., 2020) (Figure 1). The beach line is two kilometers long, going from north to south along the coast, with shallow waters and small waves (Oliveira et al., 2014). The sediment found is sandy (composed of more than 90% sand) with a predominance of coarse (69%) or fine (27%) sand . This area is the largest global producer of wild A. flexuosa and had a maximum production record of 4,716 tons in 2007 (IBAMA, 2009) until the Brazilian governmental statistical program was interrupted in 2011.
The specimens were obtained bimonthly at 34 stations between April 2016 and February 2018 ( Figure 1). Bivalve samples were taken with a cylindrical corer (19 cm diameter and 10 cm height) at each site. The total area sampled was estimated with the sum of the corer considering A = ((3.14) × r 2 ) × 34, thus making 0.9635 m 2 . The total value was used to calculate the densities that were standardized in 1 m 2 . The introduction of the corer up to 10 meters deep occurred because individuals do not exceed this depth in the region (Lima et al., 2021).
The mollusks were measured (anteroposterior axis) using calipers and the average densities (individuals per square meter) estimated for three arbitrary size classes and also all individual aggregation (total). The first size class was small individuals up to 10 mm, which are of little interest for fishing. The second class was composed of sizes larger than 10 mm and smaller than 20 mm (medium-sized), which encompasses individuals that can already be caught (Lima et al., 2020). The third class, which has sizes equal to or greater than 20 mm, is the recommended size class for catching (Silva-Cavalcanti and Costa, 2011;Oliveira et al., 2014).
In addition, monthly rainfall data from April 2015 to February 2018 were obtained from Pernambuco State Water and Climate Agency, which operates a meteorological station close to the sampling site (APAC, 2021). The rainfall time series started twelve months before the mollusk sampling period in order 3/9 Rainfall effects on Anomalocardia fexuosa in Brazilian coast to permit the evaluation of lags in the relationship between estimations of precipitation and density of A. flexuosa.

Data analysis
Autocorrelations were calculated for the five-time series, total mollusk density, the densities of the three size groups, and rainfall. Autocorrelative models can be widely used in time series to provide information on the density dependence of the populations (Bulmer, 1975), in addition to verifying relevant characteristics and possible relationships with other series. In addition, cross-correlation was calculated for all pairwise combinations. Cross-correlation analysis is the most valuable and widely used statistical tool for quantifying temporal relationships between variables of aquatic time series (Olden and Neff, 2001). The time lag with the highest cross-correlation coefficient is then typically taken as the true time lag between the two-time series (Chatfield, 2003). Up to 12 months of lag were considered when calculating autocorrelations and crosscorrelations. Analyses were conducted using a listwise deletion approach (also known as a complete case analysis) to cope with missing values. In addition, a DLM (Majid et al., 2018) was used to assess precipitation-density relationships. The DLM is a dynamic model in which the effect of a response variable y (densities) on the explanatory variable x (rainfall) may or may not be instantaneous due to the fact that the effect of change in an independent variable is not always completely exhausted within one time period, but can be distributed over several, and perhaps many, future periods (Majid et al., 2018). The selection of DLMs was conducted using the Bayesian information criterion (BIC) (Schwarz, 1978) and residual diagnostics. All the analyses were performed using the R 4.0.0 Statistical Environment software (R core team, 2020). The package dLagM (Demirhan, 2020) was used for the DLM analyses.

RESULTS
Rainfall and the average density time series are shown in Figure  2. Rainfall was generally higher from April to July, especially in May and June. The maximum value was reported for May 2016 (548.8 mm). The densities of the small shellfish were lower than the densities of the other size groups, except in April 2016, in which it had the highest density, and in February 2018, in which it presented higher density than the medium-sized specimens. In most of the months evaluated (seven collections), the densities of large shellfish (>20 mm) were higher than the values of the other groups. However, the shellfish belonging to the mediumsized class, exhibited maximum density values that occurred in the months of June (298 ind.m 2 ) in 2016 and April (301 ind.m 2 ) and June (378 ind.m 2 ) in 2017 (Figure 2).
There are significant negative autocorrelations for aggregate (total) and medium mollusk density, and also for the rainfall time series with a lag of six months ( Figures 3A, C and E). Hence, high density in a given month is followed by low density six months later. Similarly, high precipitation in a given month is followed by low precipitation six months later. There are no significant autocorrelations for small mollusk and large mollusk time series (Figures 3B and D).
The analysis of the cross-correlation shows that the rainfall was significantly associated with the total density of individuals, and that high rainfall in a given month is followed by a high total density for the same month (lag0) and up to the two months following that; though, with low density between seven and nine  months later. In addition, the data suggest that the decline in the total density can be caused by heavy rainfall during five months in the future (lead5) ( Figure 4A). The high precipitation in a month x is followed by the low density of small individuals with lags in the next three to five months, and with high density after 11 months, but an increase in the density of small individuals in rainy periods of one to three months in the future was also observed ( Figure 4B). Among the medium-sized individuals, high rainfall in a month x indicated an increase in density in lag0 and continued growth for up to two months; however it showed low seven months later, and also low density with rains of five months in the future ( Figure 4C). For the large-sized individuals, high precipitation in a given month was followed by high density five months later, but there was low density in future periods of rainfall lasting one to two months ( Figure 4D).
In the CCFs, for the relationship between densities ( Figure  4E, F and G), there were stronger significant seasonal correlations between the groups of individuals with the closest size, especially between small and medium sizes. High density of small individuals in a given month was followed by high density of medium individuals in the following two months, but the low density of medium-sized individuals may be caused by the presence of small individuals in the following months (two and four months) ( Figure 4E). Similarly, high density of small individuals two months in the future indicates low density of large individuals ( Figure 4F). The high density of medium-sized individuals in a given month is followed by the high density of large individuals four months later, but the low density of large individuals can be caused by medium-sized individuals that enter this population within two and four months ( Figure 4G).  Four DLMs were selected based on the BIC criterion for A. flexuosa densities and rainfall ). Lags of up to two months were considered in the analysis. The models for the total number of individuals (F = 4.723, DF = 3 and 8;p = 0.035) and medium individuals (F = 5.607, DF = 3 and 8;p = 0.023) were significant, while for the small (F = 1.999, DF = 3 and 8;p = 0.193) and large individuals (F = 0.656, DF = 3 and 8;p = 0.602), they were not significant ( Table 1). As for total density, there was an increase in the number of individuals with the stimulation of the rainfall for the same month (lag0) and two months later (lag2), but there was a small decrease for lag1. As for the density of the group of medium-sized individuals, there were successive increases in the number of individuals with the stimulus of the rainfall up to two months later (lag2). The best fits can be noted between the observed and predicted values of the models for the total number of individuals and medium-sized individuals ( Figure 5A and C) when compared to those of small and large individuals ( Figure 5B and D). In the best-fitting models for total densities and medium-sized individuals, 50.4 and 55.7% of the variances were explained by precipitation, respectively, which was observed in both cases in which there are higher densities between the months of April and August, with the highest peak occurring in June, and the lowest densities occurring from October to February of the following year.

DISCUSSION
The seasonality found in the autocorrelation of the rainfall was already expected because of the dry and rainy periods on the northeastern coast of Brazil. The total and medium-sized density also showed a seasonal pattern. The occurrence of autocorrelation in the densities of A. flexuosa is a reflection of the continuous reproductive cycle of the species, which leads to the constant entry of recruits and juveniles, though this varies throughout the year, depending on the climate and their habitat (Barreira and Araújo, 2005;Lavander et al., 2011;Corte et al., 2014). The non-occurrence of autocorrelation in the other size groups (small and large) of the bivalve may be a random effect of population parameters caused by natural factors (natural mortality -M) and anthropic factor (fishing mortality -F). These effects possibly affected the densities of medium-sized individuals less intensively, since they presented the maximum densities in the evaluated period and consequently this is reflected in the aggregate density.
Regarding the mortality of individuals from small and large groups, heavy rainfall and/or low salinities are variables that increase the mortality of A. flexuosa, since both factors limit the population growth of this species (Boehs et al., 2008;Belém et al., 2013). Our results indicated that precipitation at the site was not atypical and, therefore, mortality was not strongly affected by this variable. It is suspected that the decrease in small and large groups may have occurred due to other factors, such as morphological aspects and predation suffered by the species. A. flexuosa is a mollusk with a resistant shell, and this decreases its mortality rate when compared to other bivalves, such as Diplodonta punctata of the Ungulinidae family (Mattos and Cardoso, 2012), but these shells become thicker as the individuals increase in size, conferring greater protection, so it is assumed that M should occur gradually in the population of small, medium, and large individuals, respectively. Mortality in most aquatic organisms occurs mainly in smaller individuals, and also includes Veneridae bivalves (López-Rocha et al., 2018). Notwithstanding, predation also affects the population attributes of marine bivalves (Seitz et al., 2001), mainly in the group of A. flexuosa (Rodrigues et al., 2013).
Fishing activities also affect the densities of A. flexuosa, and are reflected in the high levels of exploitation in the coastal region, since the mollusk is the main fishing resource of the state in terms of production (Oliveira and Andrade, 2018); and there are already reports of growing overfishing by fishers in the region. The variation in density values by size class at the study site may be linked to the demand for the mollusk and also the fishing equipment used.
The greatest demand for the mollusk in Mangue Seco occurs at the end of the rainy period, at which time the value of this delicacy increases around 30% due to tourist activity (Silva-Cavalcanti and Costa, 2009;Oliveira et al., 2011). However, because it is a highly perishable product and rudimentarily caught through artisanal fishing (Bispo et al., 2004), there is no adequate storage and preservation logistics for the mollusk (e.g., cold stores). Consequently, bivalve exploitation increases to the extent that the demand also grows, directly influencing the decrease of the populations of A. flexuosa in this period of the year. Of course, large individuals are the most affected by the harvesting, but small individuals can also be affected by the use of fishing equipment. The use of dredges is quite widespread on site (Lima et al., 2020), and this may contribute to the mortality of individuals who escape, either by exposure to epibentonic predators or due to physical impact (Pezzuto et al., 2010;Eigaard et al., 2015), which may be more severe in small individuals who have more fragile shells. Therefore, both natural mortality and that caused by fishing affected small and large individuals to a greater degree, possibly covering up a seasonal pattern, as verified for the medium-sized group, which may have caused non-significant autocorrelation, which was presented in our results.
The results of the CCFs conducted in this study indicated that the precipitation showed a significant correlation with the mean densities of A. flexuosa, mainly with total densities and the medium-sized group, which were also significant for the DLMs, showing that rainfall plays a key role, but the role is less pronounced in the small-and large-sized groups, which indicates that these may be more influenced by other factors such as natural mortality and fishing.
Several studies on A. flexuosa carried out on the northeastern coast of Brazil have reported that the increase in rainfall caused a decrease in the density of the total number of individuals and, on Mangue Seco beach, the rainy period would have affected the large individuals to a lesser extent since they were more frequent (Oliveira et al., 2014;Rodrigues et al., 2018). The results of the models in this work do not corroborate those mentioned above, since the density of the total number of individuals increased along with rainfall, and this was promoted mainly by mediumsized mollusks, a class in which there would be an increase in densities that could occur up to two months later.
Additional fresh river water nutrient inflows can benefit bivalve populations (Cravo et al., 2006), and this a factor that coincides with the present study, since it is located between the mouth of two rivers, one of which is the Timbó River and the other the Santa Cruz channel, which is the largest estuarine complex in the state of Pernambuco (Moura, 2009). On the other hand, in the dry season, bivalve fishing focuses on catching large individuals, which when present at high densities could have less space avaliable. This factor corroborates the success of bivalve recruitment in local population stock (Boehs et al., 2008). The high density of the group of large-sized individuals that occurred in the dry season was also similar to the result found by Silva-Cavalcanti et al. (2018) for the species A. flexuosa in another extractive area in the state of Pernambuco.
However, it was observed that there are annual variations and positive and negative correlations for bivalve densities with precipitation, but these variations may also be linked to nutrient supply from rivers (Almeida et al., 2021). A. flexuosa densities can also vary significantly in the same area and according to period of study (Boehs et al., 2008;Oliveira et al., 2014), which could explain the different results obtained in the present study. Therefore, it is necessary to consider other variables such as chlorophyll-α concentrations. Nonetheless, it is worth mentioning that the present study is the first to use a statistical predictive model to investigate the relationship between density and rainfall for bivalves on the coast of Brazil, and evaluates precipitation stimuli with lags or leads over time. Therefore, the approach used here can be extrapolated to other datasets of precipitation-density relationships for A. flexuosa and other bivalves in intertidal flats.

CONCLUSIONS
The influence of rainfall on the densities of A. flexuosa, correlated with delays of this explanatory variable, is a step towards a better understanding of the populations of this bivalve on the Brazilian coast. The analyses proposed in this work considered the period as a non-fixed factor, and it was observed through the predictive models that, with the increase in rainfall, there would be an increase in total densities up to two months later, mainly with an increase in average individuals. The results of the models can serve as a reference for the adoption of management measures in the area that has the highest production of A. flexuosa in the world, considering the period and size of the individuals. Thus, it is possible to establish the best form of extraction using precipitation data, since this explained more than 50% of the variability of densities. In addition, the results of the models can serve as an initial parameter for other datasets that consider the period as a non-fixed factor in the precipitation-density relationships for A. flexuosa and other bivalves in intertidal environments.