You have count data. If the counts are "close to zero", ANOVA is clearly not recommended, although it may be acceptable when the sample size is large (what is your sample size, actually?). If the counts are far from zero, ANOVA is ok, but I'd recomment using log-transformed counts (again, for a large sample size this won't matter much).

If you simply want a p-value for the null hypothesis "the distribution of the consumption data is the same for both fish species" you can use Fisher's exact test. In your 2x2 case you will even get an odds ratio (plus CI). I think this is the simplest and best analysis for your data.

If you really want to use a GLM, then you can try a negative-binomial model, that nicely accounts for overdispersed and/or zero-inflated data. But: when there are far too many zeros, specialized tero-inflated model (hurdle models or something) will be more appropriate.

I further want to second Christian: 1) make sure you have real replications and 2) "Searching for significance" is suspect and no good scientific practice.

maybe the error bars are calculated on the basis of the normal distribution of errors, and therefore are not good for your data. As in my knowledge, anova is quite robust against normality but it is not against heteroskedasticity: being your data overdispersed,  have you tried to use a negative binomial GLM with log-link (I'm sorry but I do not know quasi-poisson)? As an alternaty you can try to log transform your data before the anova.

GLM p-values are messy and the frame work can be not understood in the same way as a standard p-value in a general linear model framework.

I have some questions: How did you exactly define your model was it number prey eaten ~ prey species * predator species. Did you do pairwise deletion of non-significant terms, i.e. first the interaction if not significant, then each main explanatory? How often did you replicate (not pseudo-replicate), i.e. what the N per level? How strong was the overdispersion scaling factor? Was the glm significant before correcting for overdispersion? Did you derive the error bar by using the glm parameter estimate output and upper and lower boundary estimates calculated from the model output, or did you calculate normal SE in a normal distribution framework and plotted that?

Finally, I am generally sceptical when somebody tries searching for significance. Count data is simply not suited to be modelled by a normal distribution, except it is well above 0 so that it is clear the fact that a count is bound at this value is not affecting your distribution.

Looking at the data itself (e.g. by plotting a simple stripplot with a mean and or median), and especially by looking at effect strength (i.e. what is the actiual difference in mean counts between groups) you often already get an idea. When you look at such a plot would you trust there is a difference? As a first subjective guess... that is often quite appealing.

In addition to the feedback you've already received, another aspect you might want to consider is if the overdispersion in your count data might be due to zero-inflation.  Perumean-Chaney, Morgan, McDowall, & Aban in the Journal of Statistical Computation and Simulation showed that not accounting for excess zeroes increased Type-II errors. Determining whether a zero-inflated model is preferable can be assessed with Vuong's test which is available in R. 

You have count data. If the counts are "close to zero", ANOVA is clearly not recommended, although it may be acceptable when the sample size is large (what is your sample size, actually?). If the counts are far from zero, ANOVA is ok, but I'd recomment using log-transformed counts (again, for a large sample size this won't matter much).

If you simply want a p-value for the null hypothesis "the distribution of the consumption data is the same for both fish species" you can use Fisher's exact test. In your 2x2 case you will even get an odds ratio (plus CI). I think this is the simplest and best analysis for your data.

If you really want to use a GLM, then you can try a negative-binomial model, that nicely accounts for overdispersed and/or zero-inflated data. But: when there are far too many zeros, specialized tero-inflated model (hurdle models or something) will be more appropriate.

I further want to second Christian: 1) make sure you have real replications and 2) "Searching for significance" is suspect and no good scientific practice.

It looks like the best distribution to describe your system is either poisson or negative binomial (in case of overdispersion). A zero-inflated or hurdle model might improve your data fit but should only be used if you have biological reasons to suspect a zero-generating process.

While the vuong test is used to compare zero-inflated and non-inflated models there is recent proof that it is not adequate for the task (see http://dx.doi.org/10.1016/j.econlet.2014.12.029 for an explanation and possible alternatives).

I personally find selection of nested models via AICc to be easier to interpret when dealing with GLMs than p-values. The package MuMIn (http://cran.r-project.org/web/packages/MuMIn/) provides several functions to enable AICc based comparisons. 

Bear in mind that model-comparisons are always relative to models present in your total set, therefore it's always a good idea to check the absolute model fit via diagnostic plots (i.e. qqplots).

I recommend robust ANOVA in R (packages 'robust' or 'WRS' http://dornsife.usc.edu/labs/rwilcox/software/) 

I too have oft encountered the notion that ANOVA is robust to deviations from the normality assumption. I have yet to see a quantitative description of "robust." I am sure it depends a great deal on the data, and sample size.

     I suggest that you answer your own question using R. Use your data to take a guess as to the underlying distribution. Use the random number generator in R to create a population of integer values that come from the hypothesized distribution. Run the analyses -- many times. You can model power and see if the observed Type 1 error is close to the expected value. Since this is simulated data, you have set up the problem such that you know the correct answer. Hopefully, you can get the analysis to give the correct answer and thereby be more confident when applying the analysis to data where you don't know the answer.

   If you are up for a larger challenge you could simulate the entire thing. This would be creating a population of values from the underlying distribution, and then withdrawing a sample from this population.

as mentionned several times above I do not know quasi-Poisson error, what do you mean about it?

if you still decide to use GLM to analyze count data, maybe you should use Box-Cox transformation: it is a power transformation which computes the best lambda to apply to your data. Besides it could also resolve heteroscedasticity if encountered.

Céline,

to understand what is meant by "quasi-Poisson models" you need to understand what a likelihood function is. In statistics you have some given data and you have some general model of how this data may be "explained" by certain parameters in the model. Such parameters can be a center, a rate, a proportion or something alike. The model relates the values of the parameter(s) to the observed data. Now the idea is, that such parameter values are good that make the actually observed data most likely. So we need some part in the model that assigns a probability value to the observations. How this can be obtained is the subject of probability theory, and it is actually solved only for very simple scenarios that are based on some simple assumptions. For instance, if "events" happening in time with a constant expected frequency, but with no information about the size of any specific interval between two events, then the expectation about the numbers of observed events within a given interval can be expressed by a Poisson distribution with the parameter lambda (that has the meaning of a rate).

If I have some count data, for which I can resonably assume that the expected frequencies of events are constant and all the same, then I can use the Poisson distribution with any arbitraily selected lambda and see how likely I expect my observed data for this value of lambda. This is called the likelihood (of the given data, depending on the model [Poisson] and selected value for the parameter [lambda]).

One can now search the values for lambda for which the obtained likelihood is maximal. And based on (typically only approximated) properties of the likelhood function one can derive standard errors or confidence intervals for such "maximum lihelihood estimates" and perform hypothesis tests (which are then likelihood ratio tests).

The Poisson model is based on such simple and restrictive assumptions that there is a strict and neccesary relationship between the expectation and the variance of a variable according to a Poisson distribution (they are equal).

The Poisson model can be modified to allow for any other linear relationship between the expectation and the variance. This is achieved by adding a another parameter for the dispersion of the data to the Poisson model (actually just to allow the dispersion parameter in the general linear model form to take values different from 1). The resulting likelihood function is then not based on a specific probability model that was derived by theoretical consideration from simple assumptions. Therefore it is called a "quasi-likelihood" and the model is called a "quasi-Poisson model".

The original quasi-Poisson model estimates the dispersion parameter as some linear function of the relationship between the expectation ("mean") and the variance. The negative biniomial model is more flexible, since it estimates it from a quadratic relationship (what can include a linear relatinship). But one is not restricted to these. One could define other mean-variance relationships and use them in a more general setting of dealing with GLMs.

Regaring the Box-Cox transformation:

Advantage: it often works and gives more reliable error-estimates.

Disadvantage 1: transfroming data to make it suit a particular analysis renders the p-values and confidence intervals meaningless if they are obtained from the same data (split your dataset or use new data to test hypotheses).

Disadvantage 2: the meaning of the coefficients for transformed data is usually completely unclear. This strongly limits any scientific interpretation. Exceptionas are the log-transformation (changing a multiplicative model into an additive one) and the log-odds transformation (maybe there are more but I am only aware of these two leading to scientifically interpretable estimates).

If I well understood when you talk about "mean-variance relationship" you are talking about the link function, isn't it? If yes, why the link function regarding Poisson distrubution is log by default in many softwares? but maybe I totally misunderstand what you mean.

I understand concerning Box-Cox transformation, it was just in case still going with ANOVA in this kind of data.

No, the link function is the function that related the model predictions to a "linear predictor" (the eta-function). This is not related to the mean-variance relationship of the Poisson model. The link function is used to make sure that any possible predictor value(s) will give predictions that are within the support of the distribution. For the Poisson model the typical link function is "log". So the values predicted by the linear predictor can be positive or negative; since these predictions are on the log-scale (as specified by the link function), the prediction on the original scale are always neccesarily positive (there can not be negative counts).

The mean-variance relationship means:

E(X) = VAR(X) = lambda.

This is the case when X ~ Pois(lambda).

For instance, there is no such relationship for X ~ Norm(m,s²). The normal distribution has two parameters, one for the expectation (m) and the other for the variance (s²) - and they are independent (E(X)=m and VAR(x)=s²). Knowing m tells us nothing about s² (and vice versa). But when X ~ Pois(lambda) and this X has a given expectation, then we know that the variance must have the very same value.

ok now it is perfectly clear concerning especially the Relationship between mean and variance. a last question, you wrote in your in last mail:

"The original quasi-Poisson model estimates the dispersion parameter as some linear function of the relationship between the expectation ("mean") and the variance. The negative biniomial model is more flexible, since it estimates it from a quadratic relationship (what can include a linear relatinship). But one is not restricted to these. One could define other mean-variance relationships and use them in a more general setting of dealing with GLMs".

which kind of other mean-variance relationships do you mean? in what situations?

especially for counts data, probability of distribution of Poisson is involved. with continuous data, how do you manage deviation for normality and/or heteroscedasticity (without using data transformation of course)?