There is no exact science behind including or omitting explanatory variables, as adding additional variables changes the meaning of your model. By adding additional variables, what you basically say is: How does x affect y, keeping all other explanatory variables constant?
Apart from paying attention to the p-values, you should also have a look at the $ R^2 $, which measures which portion of the variance in the population is explained by your model.
Unless you are interested in coding the regression algorithm yourself, you might want to have a look at the excellent Weka Java Library or R, an open source statistics software.
Since you have multiple dependent and independent variables, a multivariate analysis would be one way to proceed.
A multivariate analysis will attempt to model the relationship between your dependent and independent variables, and as an outcome you will be able to test if those factors are significant in your model.
This is useful if you want to assess the significance of the factors within such a model, but if you are interested in knowing the significance of the relationship between the covariates and one response you can run a regression the way you describe.
Suppose though, that you want to construct a model for both responses simultaneously, and assess the significance of the factors in $that$ model.
Then you can use multivariate analysis of covariance (MANCOVA). MANCOVA will provide you with the contribution to the variance in the responses made by each factor, as well as their significance. Note that MANCOVA will produce both type I, II, and III sums of squares (SS). Which one is appropriate depends on the balance of your data. See here for more information on the types of SS.
http://mcfromnz.wordpress.com/2011/03/02/anova-type-iiiiii-ss-explained/
If you are using R, you can determine the statistical significance of your factors by performing multivariate regression and using this as input in the manova function. The code would go something like:
#fit a multivariate regression model and then test the type I SS using MANCOVA.
fit = lm(formula = cbind(Abundance, Richness) ~ Temp_1 + Rain_1 + Sunlight_1 + Temp_2 + Rain_2 + Sunlight_2 + Temp_3 + Rain_3 + Sunlight_3 + Temp_4 + Rain_4 + Sunlight_4, data = yourData)
summary(manova(fit), test="Hotelling-Lawley")
A more thorough overview of how to perform such an analysis is provided here:
http://www.uni-kiel.de/psychologie/rexrepos/posts/multRegression.html
Note that separate regressions return the same slopes as multivariate regression, and also not that different tests besides the "Hotelling-Lawley" are possible for the MANCOVA test of type I SS, and that you can also test type II SS.
As you pointed out, PCA is another multivariate data analysis method. It may be interesting for performing exploratory data analysis, though you will not obtain the same sort of significance testing discussed above. If you perform PCA on your data, a bi-plot may be a good way to investigate interesting relationships.
SAS provides some rather clear discussion interpreting the biplot:
http://support.sas.com/documentation/cdl/en/imlsug/62558/HTML/default/viewer.htm#ugmultpca_sect2.htm
Best Answer
Of course, it all depends on your research question, but here are a few ideas.
Basic set of analyses:
Selection of more sophisticated options: