Don't trust a test like shapiro-wilks to tell you if your data meets model assumptions.  It will steer you wrong. 

But also, normality of residuals is not the only assumption in analysis of variance.  Homoscedasticity and other assumptions should be checked.

For example, on your Trial 1, the residuals really weren't too far from normal, but the heteroscedasticity was horrible.

I recommend using diagnostic plots.  Of course, you have to understand what they mean.

source("http://rcompanion.org/r_script/normal.histogram.r")

normal.histogram(x)

plot(fitted(m1),

residuals(m1))

qqnorm(x)

qqline(x, col="red")

plot(m1)

Trying a few models, it looks like doing the three-way anova approach with the log-transformed dependent variable is acceptable.

But another issue I see is that your concentrations are almost entirely different for each treatment.  I'm not sure, but you may not want to include any terms with conc:treat then.  It doesn't really make much sense to include them.  I.e.

m4 = aov(ll ~ conc + treat + variable + conc:variable + treat:variable, data=md2)

To answer your new questions:

1) You shouldn't reduce your model to a one-way anova, unless you really have to for some reason. There's no sense in designing a factorial (sort of) experiment if you are going to reduce the analysis to a one-way anova.

2) You shouldn't reduce your model as in trial 2, for the same reason.

3) Yes, but this is not the best way to analyze your data, for the same reason.

4) Yes, as long as "length" for each of them is really the same dependent variable, and it would make sense to compare root length to shoot length.

Your two-way anova doesn't really make a lot of sense as written. Either

A) Use the three-way anova proposed in #4, or

B) Consider root length to be a totally different dependent variable than shoot length, and do the analyses for root, shoot and seedling totally separately. In this case, you cannot compare the length of roots to the length of shoots.

General advice:

Find some good sources on the assumptions of parametric tests like analysis of variance.  Check all model assumptions.  You need to plot residuals to see if your model adequately explains the data.

Think about what hypotheses you are trying to test. For example, if you want to know if shoot length differed from root length, you need to use a three-way anova suggested in #4.  Also, does it make any sense to include conc:treat?  What would a significant result for this interaction tell you?  (From your design, I think the answer is, nothing useful).

Plot your data by groups and interactions.  Include standard errors of the mean or confidence intervals.  This way, you can see what the effects are of the terms in you model.

Question 1) WOULD IT BE MEANINGFUL TO USE THE "conc:variable" TERM (SINCE THE CONCENTRATIONS DIFFERS IN EACH TREATMENT GROUP)

IN THE CODE SUGGESTED BY YOU: "m4 = aov(ll ~ conc + treat + variable + conc:variable + treat:variable, data=md2)"?

Question 2) SHOULD I USE THE SAME TRANSFORMED DATA (USED FOR TWO-WAY ANALYSIS) FOR ONE-WAY MODEL REQUIRED DURING MULTIPLE  COMPARISON USING "glht" EVEN THOUGH THE UNTRANSFORMED ONE-WAY

MODEL MEETS THE ASSUMPTION?

Question 3) Actually I was trying to see how certain chemical treatments affect the growth of root, shoot and seedling and whether there is any interaction in the outcome. Since the "conc:treat" term did not make sense (as you have pointed out), I modified the analysis by splitting (as the levels- concentrations- were not comparable between treatments) the data by each treatment. Thus the dataset included the respose variables (root, shoot and seedling) at different concentrations for ONE treatment. Thus the ANOVA reduced to two-way (The codes are given below) and required log transformation (except for t2).

NOW WHAT TO DO WITH THE "t2" DATA WHICH COULD NOT BE FITTED EVEN AFTER LOG TRANSFORMATION?

##Now each treatment with its control and conc. and 3 response variables (factor)

read.csv("test.csv")->md2

rbind(md2[md2$treat=="control",],md2[md2$treat=="t1",])->t1

rbind(md2[md2$treat=="control",],md2[md2$treat=="t2",])->t2

rbind(md2[md2$treat=="control",],md2[md2$treat=="t3",])->t3

rbind(md2[md2$treat=="control",],md2[md2$treat=="t4",])->t4

##NONE of the data (t1,t2,t3) met normality, though homogeneity was met

aov(length~conc*variable,data=t1)->m1

shapiro.test(resid(m1))##ASSUMTION VIOLATED

library(nortest)

ad.test(resid(m1))##ASSUMTION VIOLATED

lillie.test(resid(m1))##ASSUMTION VIOLATED

##log transformation produced a mormal fit (EXCEPT FOR "t2")

aov(log(length+1)~conc*variable,data=t1)->m1

shapiro.test(resid(m1))##ASSUMTION MET

Question 1) If the hypothesis is of interest to you.  It tests if the different "variable"s behave the same way across the concentrations. 

Question 2) Usually if the interactions are significant, you want to do multiple comparisons on the interaction terms, and not on the main effects.  In any case, you should not redo your two-way anova as a one-way anova.  You may want to use lsmeans instead of glht.  See http://rcompanion.org/rcompanion/d_08.html for an example.  TukeyHSD can also be used, (It may not be appropriate for unbalanced designs).

Question 3)  I think that approach makes sense.  Once you have separated the treatments into different analyses, they are really separate analyses, and there is no need to use the same transformation on all of them.  It might be weird to present one as log and one as not transformed, tho.

Besides estimating the interaction, I wanted to know the-

No Observed Effect concentration (NOEC)for each variable as well.

I think the marginal means produced by pooling either all the concentration per variable OR all the variables per concentration make no sense. I would also like to know the back transformation method for the output produced.

Thanking you in advance

read.csv("test.csv")->md2 ; md2$conct2

rbind(md2[md2$treat=="control",],md2[md2$treat=="t3",])->t3

rbind(md2[md2$treat=="control",],md2[md2$treat=="t4",])->t4

1) IS THIS THE RIGHT METHOD OF TUKEY'S HSD COMPARISON ?

##Tukey

m1 = lm(log(length+1) ~ conc+ variable + conc:variable, data=t1)

 library(lsmeans)

lsm = lsmeans(m1, pairwise ~ conc:variable, adjust="tukey")

lsm$contrasts

cld(lsm, alpha=.05, Letters=letters)

2)  ARE THE RESULTS GIVEN BELOW  (A,B or C) ARE EQUIVALENT FOR DUNNET’S MULTIPLE COMPARISON?

A)

library(multcomp)

m1