[R-sig-ME] Modelling intraspecific differences with random slopes

[David R.]

Hi, 
I'm using mixed models (lme4 package) to analyze variability in 13 SPECIES of birds observed during 15 YEARS across 5 sites. All the species were observed in all the sites in most YEARS. 
My initial model was: 
response ~ a + b + c + d + e + (1 | YEAR) + (1 | SITE) + (1 | SPECIES) 
that after some LRT was simplified to: 
response ~ a + b + c + d + e + (1 | SPECIES) 
I was not interested in these species in their own right and treated them as being representative members of a population of similar species. But now I was asked about the possible intraspecific differences in the effect of a, b, c, d and e on the response. 
My question is: Is it appropriate a model of random intercept and slopes as initial full model to estimate these differences? 
For example: 
response ~ a + b + c + d + e + (1 | YEAR) + (1 | SITE) + (1 | SPECIES) + (1 + a | SPECIES) + (1 + b | SPECIES) + (1 + c | SPECIES) + (1 + d | SPECIES) + (1 + e | SPECIES) 
Or perhaps is it more appropriate to work with species as fixed effects like: 
response ~ (a + b + c + d + e) * SPECIES + (1 | YEAR) + (1 | SITE) 

Thanks in advance 
David