Heteroscedasticity, lme4 and nlme

I am working on my dissertation and using mixed models for quite a bit of
my analysis. Here are some details about my experiment and data:

I grew seedlings of 21 species in two sites, 10 blocks per site, and 2
plots per block either with or without competition ("Weeds").  I am also
including phylogenetic family (3 families) and invasiveness (3
categories).  I have site and family as fixed effects because I understand
there aren't enough categories to reliably estimate them as random
effects.  I have several potential continuous covariates as well.  My model
includes all possible two-way interactions, fixed and random, and I am
reducing the model using the likelihood ratio test.  My goal in reducing
the model is to improve clarity and avoid over-fitting.  I have both
continuous (transformed) and binomial response variables.

Here's a full (unreduced) model without covariates:
RGR.simple <- lmer(RGR ~ Inv + Weeds + Site + Fam + Inv:Weeds + Inv:Site +
Weeds:Site + Weeds:Fam + Site:Fam + (1|Blk) + (1|Plot) + (1|Sp) + (0 +
Site|Sp) + (0+Weeds|Sp), fdatRGR, family="gaussian")

My data are extremely unbalanced. I have heteroscedasticity in the
residuals of my full and reduced models, as seen both by plotting and using
leveneTest (package car).

I have several questions, which center largely around the issue of
heteroscedasticity:

I have both crossed random effects and serious heteroscedasticity.  I'd
like to implement a non-homogeneous variance structure, but as far as I can
tell, that can only be done in lme (from nlme), not in lmer (from lme4).
Is that right?  I have heard that there are some work-arounds to get
crossed random effects in lme.  

Would that be the best solution in this
situation?  If I specified an alternate variance structure by species,
should residuals by species be homogenous (ie non-significant Levene's
Test) or not?

I've also gotten a bit confused about random interactions.  lme4 allows
much more complex random effects - would that in some way help account for
heteroscedasticity?  If so, it hasn't worked for my data, but I want to
understand what the random effects are doing.  In other words, what is the
difference between specifying a variance structure (e.g. varIDent(form =
~1|Species) and including vector-valued random effects (e.g.
(Competition|Species))?  My model statistics (ie LRT and AIC) tell me that
many of the models are better with these complex random effects, but that
also restricts my movement from lme4 to nlme.

There is a weights function in lmer, but I know it is not the same as the
weights function in lme.  If I use the weights term in lmer to specify
weighting by the scaled reciprocal of species variance, that would give
species with smaller variances a greater "say" in the model - is that
right?  Weighting this way does fix my heteroscedasticity problems, but I'm
not sure if it is a valid approach.