> library(rsm)
> 
> tilapia.rsm1 <- rsm(mean_WBC ~ SO(temp,protein))
> summary(tilapia.rsm1)

Call:
rsm(formula = mean_WBC ~ SO(temp, protein))

Residuals:
    Min      1Q  Median      3Q     Max 
-12.238  -9.668  -5.398  10.541  14.683 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)   
(Intercept)  -615.12410  248.45437  -2.476  0.04247 * 
temp           47.20166   12.95047   3.645  0.00823 **
protein         8.89176    6.58767   1.350  0.21911   
temp:protein   -0.05245    0.15830  -0.331  0.75010   
temp^2         -0.75344    0.21237  -3.548  0.00937 **
protein^2      -0.09036    0.06644  -1.360  0.21602   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 13.65 on 7 degrees of freedom
Multiple R-squared: 0.8396,     Adjusted R-squared: 0.7251 
F-statistic:  7.33 on 5 and 7 DF,  p-value: 0.01050 

Analysis of Variance Table

Response: mean_WBC
                   Df Sum Sq Mean Sq F value   Pr(>F)
FO(temp, protein)   2 4319.9 2159.94 11.5900 0.006009
TWI(temp, protein)  1   20.5   20.46  0.1098 0.750101
PQ(temp, protein)   2 2490.2 1245.08  6.6809 0.023822
Residuals           7 1304.5  186.36                 
Lack of fit         3  708.9  236.29  1.5868 0.325086
Pure error          4  595.7  148.92                 

Stationary point of response surface:
    temp  protein 
29.91379 40.52278 

Eigenanalysis:
$values
[1] -0.08931963 -0.75447368

$vectors
            [,1]        [,2]
[1,]  0.03945517 -0.99922134
[2,] -0.99922134 -0.03945517