What model(s) are appropriate for unbounded species abundance data with 0's and environmental response variables in R? I am conducting a research project on the vegetation composition of an abandoned mine site to determine the environmental factors that influence species abundance (cover) and composition. 
Study Design: I have 10 sites, within each site 10 randomly placed 1x1m quadrats were sampled for plant species percent cover (abundance) along with a soil sample that was tested for heavy metals concentration(Pb, As, Cu, Zn), pH, NPK, and OM. Factors such as slope, aspect were also recorded. 
Analysis: From species abundance data I created the response variables: total cover, exotic cover, native cover, and species richness. I would like to investigate if my environmental variables are influencing the response variables I described above. I have tried using lm in R but my response variables are not normally distributed and therefore I believe violate the assumptions of linear modeling. I have also tried glm but since my cover abundance data contains many 0's and is not bounded at 100%, as there can be more than 100% cover in a quadrat, I don't think any of the error distribution families of glm can accommodate this.
So, I would appreciate any input on selecting an appropriate model for my response variables. If possible I would like to see the size of the effect or r2 rather than just a P value. Thank you very much for your time. If you would be interested in the data I could provide it. 
 A: You can try to transform your responses.
If it is continuous and right skeewed, try a cube root transformation since it deals with zeros.
Otherwise, check if the species abundance in counts are hightly correlated with their cover %, if they are, you can try a count distribution family, such as poisson or negative binomial in case of overdispersion.
How is it possible to have a cover greater than 100%? (curiosity)
Another detail, are you inspecting for normality in the model residuals or at response variable itself?
You should inspect for normality at the model residuals, such as:
plot(model) #residuals versus fitted should show homocedasticity
hist(resid(model)) #should looks like a bell shape
shapiro.test(resid(model)) # not significative= normality

Since your quadrats are nested in site (non indenpendence-pseudoreplication), it is worth to consider a mixed modeling approach, where you should use sites as random factors.
Good luck!
A: First, ordinary least squares regression does not make assumptions about the distribution of the variables, only about the errors - estimated by the residuals.
Second, assuming that your residuals are not normal, then, rather than transforming your variable, I suggest you use a method that makes fewer assumptions. One such is quantile regression - I have done a lot of this, but in SAS, not R, but there is an R package quantreg that may help.  Quantile regression makes no assumptions about the residuals. 
